Coverage for dragonfly_sim/core/turbulence_models.py: 97%
95 statements
« prev ^ index » next coverage.py v7.16.2, created at 2026-10-17 03:37 +0000
1"""
2Turbulence closures for CompressibleRANSModel.
4A model contributes n_vars conservative transported scalars rho*phi_k. They
5are appended after the mean flow in the state vector and solved fully
6coupled with it. The RANS model asks the closure for the following, all UFL
7expressions of the state U (mean flow in U[:n_mean], scalars after it):
9 eddy_viscosity(U, mu, n_mean) mu_t
10 diffusive_fluxes(U, mu, n_mean) -> f(grad_phi, grad_rho), one
11 viscous-flux row per scalar
12 sources(U, grad_U, mu, d, n_mean, gamma) one source per scalar (right-hand side)
13 freestream_values(fs) conservative rho*phi_inf per scalar
14 wall_values(U, n_mean) conservative values on a no-slip wall
15 limiter_scales(fs) a magnitude per scalar, for the MUSCL
16 limiter's smoothing parameter
18`fs` is a dict with the free stream: rho, p, c, M, speed (M c) and mu (a
19dolfinx Constant: mu_inf follows set_mach). Convection needs nothing from
20the closure: every scalar is carried with the mean flow, rho*phi*u
21(NS_utils.extended_flux), and upwinded with the mean flow's wave speeds. A
22closure whose sources need the distance to the wall sets
23needs_wall_distance.
25Laminar is the zero-equation case (pure Navier-Stokes). Adding a closure,
26k-omega SST say, is a new subclass with n_vars = 2 and names
27("rho_k", "rho_omega"). Nothing in CompressibleRANSModel changes: the state
28size, boundary values, limiter scales, the Jacobian and the adjoint all
29follow from this interface.
30"""
31import ufl
33from dragonfly_sim.utils.NS_utils import primitive_gradients, vorticity_magnitude
36class TurbulenceModel:
37 """Base class and interface (see the module docstring)."""
38 n_vars = 0
39 names = ()
40 needs_wall_distance = False
42 def eddy_viscosity(self, U, mu, n_mean):
43 return 0.0
45 def diffusive_fluxes(self, U, mu, n_mean):
46 return lambda grad_phi, grad_rho: []
48 def sources(self, U, grad_U, mu, d, n_mean, gamma):
49 return []
51 def freestream_values(self, fs):
52 return []
54 def wall_values(self, U, n_mean):
55 return []
57 def limiter_scales(self, fs):
58 """Default: 1000x the free-stream value, i.e. the scalar's size in a
59 boundary layer whose eddy viscosity is ~1e3 times the laminar one."""
60 return [1e3 * v for v in self.freestream_values(fs)]
63class Laminar(TurbulenceModel):
64 """No turbulence model: the laminar compressible Navier-Stokes equations."""
67class SpalartAllmarasNeg(TurbulenceModel):
68 """
69 The negative Spalart-Allmaras model, SA-neg (Allmaras, Johnson & Spalart,
70 ICCFD7-1902, 2012), in the conservative compressible form
72 div(rho nu~ u) = rho (P - D)
73 + (1/sigma) div[(mu + rho nu~ f_n) grad nu~]
74 + (c_b2/sigma) rho |grad nu~|^2
75 - (1/sigma) (nu + nu~ f_n) grad rho . grad nu~
77 For nu~ >= 0 this is standard SA (without the f_t2 term unless ft2=True),
78 with the S~ modification that keeps S~ positive. For nu~ < 0 it switches
79 to the negative-branch production, destruction and diffusion, which make
80 a negative nu~ decay instead of growing. Under-resolved boundary-layer
81 edges routinely undershoot below zero, which is why SA-neg is the usual
82 choice. mu_t is zero there.
84 The two branches meet continuously at nu~ = 0 (production and
85 destruction both vanish), so the conditionals add kinks, not jumps.
87 Free stream: nu~_inf = chi_inf nu_inf, chi_inf = 3 (NASA TMR). Wall:
88 nu~ = 0.
89 """
90 n_vars = 1
91 names = ("rho_nu_tilde",)
92 needs_wall_distance = True
94 cb1 = 0.1355
95 sigma = 2.0 / 3.0
96 cb2 = 0.622
97 kappa = 0.41
98 cw2 = 0.3
99 cw3 = 2.0
100 cv1 = 7.1
101 ct3 = 1.2
102 ct4 = 0.5
103 cn1 = 16.0
104 cv2 = 0.7
105 cv3 = 0.9
106 r_lim = 10.0
108 def __init__(self, chi_inf=3.0, ft2=False):
109 self.chi_inf = float(chi_inf)
110 self.ft2 = bool(ft2)
111 self.cw1 = self.cb1 / self.kappa**2 + (1.0 + self.cb2) / self.sigma
113 # --- pointwise closure functions (also what the unit tests check) ----
114 def fv1(self, chi):
115 chi3 = chi**3
116 return chi3 / (chi3 + self.cv1**3)
118 def fv2(self, chi):
119 return 1.0 - chi / (1.0 + chi * self.fv1(chi))
121 def ft2_fn(self, chi):
122 return self.ct3 * ufl.exp(-self.ct4 * chi**2) if self.ft2 else 0.0
124 def fn(self, chi):
125 chi3 = chi**3
126 return ufl.conditional(ufl.ge(chi, 0.0), 1.0,
127 (self.cn1 + chi3) / (self.cn1 - chi3))
129 def s_tilde(self, Omega, nu_t, chi, d):
130 S_bar = nu_t * self.fv2(chi) / (self.kappa**2 * d**2)
131 modified = Omega + Omega * (self.cv2**2 * Omega + self.cv3 * S_bar) / (
132 (self.cv3 - 2.0 * self.cv2) * Omega - S_bar)
133 return ufl.conditional(ufl.ge(S_bar, -self.cv2 * Omega), Omega + S_bar, modified)
135 def fw(self, r):
136 g = r + self.cw2 * (r**6 - r)
137 return g * ((1.0 + self.cw3**6) / (g**6 + self.cw3**6))**(1.0 / 6.0)
139 def production_destruction(self, nu_t, nu, Omega, d):
140 """
141 (P, D) per unit density, with both SA-neg branches. S~ is floored
142 before it divides r, which only matters on the unselected branch.
143 """
144 chi = nu_t / nu
145 St = self.s_tilde(Omega, nu_t, chi, d)
146 r = ufl.min_value(nu_t / ufl.max_value(St * self.kappa**2 * d**2, 1e-30), self.r_lim)
147 ft2 = self.ft2_fn(chi)
148 P_pos = self.cb1 * (1.0 - ft2) * St * nu_t
149 D_pos = (self.cw1 * self.fw(r) - self.cb1 * ft2 / self.kappa**2) * (nu_t / d)**2
150 P_neg = self.cb1 * (1.0 - self.ct3) * Omega * nu_t
151 D_neg = -self.cw1 * (nu_t / d)**2
152 positive = ufl.ge(nu_t, 0.0)
153 return (ufl.conditional(positive, P_pos, P_neg),
154 ufl.conditional(positive, D_pos, D_neg))
156 # --- the TurbulenceModel interface -----------------------------------
157 def _nu_tilde(self, U, n_mean):
158 return U[n_mean] / U[0]
160 def eddy_viscosity(self, U, mu, n_mean):
161 rho = U[0]
162 nu_t = self._nu_tilde(U, n_mean)
163 chi = nu_t * rho / mu
164 return ufl.conditional(ufl.ge(nu_t, 0.0), rho * nu_t * self.fv1(chi), 0.0)
166 def diffusive_fluxes(self, U, mu, n_mean):
167 rho = U[0]
168 nu_t = self._nu_tilde(U, n_mean)
169 coeff = (mu + rho * nu_t * self.fn(nu_t * rho / mu)) / self.sigma
170 return lambda grad_phi, grad_rho: [coeff * grad_phi[0]]
172 def sources(self, U, grad_U, mu, d, n_mean, gamma):
173 rho = U[0]
174 nu = mu / rho
175 nu_t = self._nu_tilde(U, n_mean)
176 grad_u, _, grad_phi, grad_rho = primitive_gradients(U, grad_U, gamma, n_mean)
177 grad_nt = grad_phi[0]
178 Omega = vorticity_magnitude(grad_u)
179 # d is zero only on the wall itself; quadrature points never sit there
180 d_safe = ufl.max_value(d, 1e-14)
181 P, D = self.production_destruction(nu_t, nu, Omega, d_safe)
182 fn = self.fn(nu_t / nu)
183 return [rho * (P - D)
184 + self.cb2 / self.sigma * rho * ufl.dot(grad_nt, grad_nt)
185 - (nu + nu_t * fn) / self.sigma * ufl.dot(grad_rho, grad_nt)]
187 def freestream_values(self, fs):
188 return [self.chi_inf * fs['mu']] # rho_inf * (chi_inf * mu_inf / rho_inf)
190 def wall_values(self, U, n_mean):
191 return [0.0]