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.
Source code in g-adopt/gadopt/momentum_equation.py
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 | def viscosity_term(
eq: Equation, trial: Argument | ufl.indexed.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 are equivalent to free stress (i.e. free in all
# directions).
# NOTE: "un" can be combined with "stress" provided the stress force 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 = outer(eq.n, trial - bc["u"])
if compressible_stress:
trial_tensor_jump -= (
2 / 3 * identity * (dot(eq.n, trial) - dot(eq.n, bc["u"]))
)
trial_tensor_jump += transpose(trial_tensor_jump)
# 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 = 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"])
trial_tensor_jump += transpose(trial_tensor_jump)
# 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 "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
|