25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134 | def viscosity_term(eq: Equation, trial: Argument | Indexed | Function) -> Form:
r"""Viscosity term $-nabla * (mu nabla u)$ in the momentum equation.
Using the symmetric interior penalty method (Epshteyn & Rivière, 2007), the weak
form becomes
$$
{:( -int_Omega nabla * (mu grad u) phi dx , = , int_Omega mu (grad phi) * (grad u) dx ),
( , - , int_(cc"I" uu cc"I"_v) "jump"(phi bb n) * "avg"(mu grad u) dS
- int_(cc"I" uu cc"I"_v) "jump"(u bb n) * "avg"(mu grad phi) dS ),
( , + , int_(cc"I" uu cc"I"_v) sigma "avg"(mu) "jump"(u bb n) * "jump"(phi bb n) dS )
:}
$$
where σ is a penalty parameter.
Epshteyn, Y., & Rivière, B. (2007).
Estimation of penalty parameters for symmetric interior penalty Galerkin methods.
Journal of Computational and Applied Mathematics, 206(2), 843-872.
"""
dim = eq.mesh.geometric_dimension
identity = Identity(dim)
compressible_stress = eq.approximation.compressible
mu = eq.approximation.mu
stress = eq.stress
F = inner(nabla_grad(eq.test), stress) * eq.dx
sigma = interior_penalty_factor(eq)
sigma *= FacetArea(eq.mesh) / avg(CellVolume(eq.mesh))
if not is_continuous(eq.trial_space):
trial_tensor_jump = tensor_jump(eq.n, trial) + tensor_jump(trial, eq.n)
if compressible_stress:
trial_tensor_jump -= 2 / 3 * identity * jump(trial, eq.n)
F += (
sigma
* inner(tensor_jump(eq.n, eq.test), avg(mu) * trial_tensor_jump)
* eq.dS
)
F -= inner(avg(mu * nabla_grad(eq.test)), trial_tensor_jump) * eq.dS
F -= inner(tensor_jump(eq.n, eq.test), avg(stress)) * eq.dS
# NOTE: Unspecified boundaries result in free stress (i.e. free in all directions).
# NOTE: "un" can be combined with "stress" provided the stress component is
# tangential (e.g. no normal flow with wind)
for bc_id, bc in eq.bcs.items():
if "u" in bc and any(bc_type in bc for bc_type in ["stress", "un"]):
raise ValueError(
'"stress" or "un" cannot be specified if "u" is already given.'
)
if "normal_stress" in bc and any(bc_type in bc for bc_type in ["u", "un"]):
raise ValueError(
'"u" or "un" cannot be specified if "normal_stress" is already given.'
)
if "u" in bc:
trial_tensor_jump = 2 * sym(outer(eq.n, trial - bc["u"]))
if compressible_stress:
trial_tensor_jump -= (
2 / 3 * identity * (dot(eq.n, trial) - dot(eq.n, bc["u"]))
)
# Terms below are similar to the above terms for the DG dS integrals.
F += (
2
* sigma
* inner(outer(eq.n, eq.test), mu * trial_tensor_jump)
* eq.ds(bc_id)
)
F -= inner(mu * nabla_grad(eq.test), trial_tensor_jump) * eq.ds(bc_id)
F -= inner(outer(eq.n, eq.test), stress) * eq.ds(bc_id)
if "un" in bc:
trial_tensor_jump = 2 * outer(eq.n, eq.n) * (dot(eq.n, trial) - bc["un"])
if compressible_stress:
trial_tensor_jump -= 2 / 3 * identity * (dot(eq.n, trial) - bc["un"])
# Terms below are similar to the above terms for the DG dS integrals.
F += (
2
* sigma
* inner(outer(eq.n, eq.test), mu * trial_tensor_jump)
* eq.ds(bc_id)
)
F -= inner(mu * nabla_grad(eq.test), trial_tensor_jump) * eq.ds(bc_id)
# We only keep the normal part of stress; the tangential part is assumed to
# be zero stress (i.e. free slip) or prescribed via "stress".
F -= dot(eq.n, eq.test) * dot(eq.n, dot(stress, eq.n)) * eq.ds(bc_id)
if hasattr(eq.approximation, 'bulk_modulus'):
trial_tensor_jump = identity * (dot(eq.n, trial) - bc["un"])
trial_tensor_jump += transpose(trial_tensor_jump)
bulk = eq.approximation.bulk_modulus * eq.approximation.bulk_shear_ratio
# Terms below are similar to the above terms for the DG dS integrals.
F += (
2
* sigma
* inner(outer(eq.n, eq.test), bulk * trial_tensor_jump)
* eq.ds(bc_id)
)
F -= inner(bulk * nabla_grad(eq.test), trial_tensor_jump) * eq.ds(bc_id)
if "stress" in bc: # a momentum flux, a.k.a. "force"
# Here we need only the third term because we assume jump_u = 0
# (u_ext = trial) and stress = n . (mu . stress_tensor).
F -= dot(eq.test, bc["stress"]) * eq.ds(bc_id)
if "normal_stress" in bc:
F += dot(eq.test, bc["normal_stress"] * eq.n) * eq.ds(bc_id)
return F
|