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 | 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 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 = 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
|