Coverage for dragonfly_sim/core/RANS_model.py: 89%

1""" 

2Compressible RANS (and laminar Navier-Stokes) at p = 0, built on the Euler 

3model by composition. 

4 

5CompressibleRANSModel owns a CompressibleEulerModel (`self.flow`) whose state 

6carries the turbulence model's transported scalars after the mean flow 

7(rho*nu_tilde for Spalart-Allmaras), and plugs the viscous physics into it. 

8Everything else -- spaces, boundary conditions, the Newton solver and 

9positivity limiter, pseudo-transient continuation, time stepping, adjoint 

10and shape-derivative assembly -- is the Euler model's own code. Attributes 

11and methods the RANS model does not define itself are forwarded to 

12`self.flow` (reads and writes), so DG_windtunnel_model and DG_postprocessor 

13use it exactly like the Euler model. 

14 

15Discretization (cell-centred finite volume, as DG0) 

16--------------------------------------------------- 

17Inviscid: NS_utils.extended_hll_flux (HLL on the mean flow, the same fan 

18upwinding the scalars), with MUSCL-reconstructed interior-facet traces. 

19Without MUSCL the first-order HLL flux puts acoustic-speed dissipation on 

20the tangential momentum at wall-parallel facets, a numerical viscosity that 

21exceeds the turbulent one through the boundary layer (SA NACA0012 at M 0.7: 

22friction drag ~10x too high). 

23 

24Viscous: through each facet, F_v(U_f, grad U_f) . n with the reconstructed 

25facet gradient of core/reconstruction.py; no penalty term. Turbulence 

26sources use the Green-Gauss gradient G. G is condensed out of Newton 

27exactly, and is a derived field for adjoints and shape derivatives, as are 

28the cell centres and the wall distance. 

29 

30Nondimensionalization: rho_inf = p_inf = 1, chord 1, U_inf = M sqrt(gamma), 

31mu_inf = rho_inf U_inf / Re (Re held fixed under set_mach). 

32 

33Boundary conditions (kind: extra-variable values) 

34------------------------------------------------- 

35farfield the characteristic state; scalars from the free stream where it 

36 flows in, from the interior where it flows out 

37wall adiabatic no-slip (define_wall_bc): u = 0, interior density and 

38 internal energy, the turbulence model's wall values; zero heat flux 

39slip (symmetry) the Euler slip state, interior scalars, no viscous flux 

40 

41Solver defaults 

42--------------- 

43Pseudo-transient continuation from the free stream (use_ptc, CFL0 5), one 

44Newton step per pseudo step; FGMRES on the exact condensed Jacobian, 

45preconditioned with ASM(1)/ILU(2) of the compact gradient-frozen part 

46(jacobian_mode "compact_pc", utils/gradient_condensation.py). Set 

47direct_linear_solver for MUMPS on the exact Jacobian instead. 

48""" 

49import ufl 

50import dolfinx 

51 

52from dragonfly_sim.core.Euler_model import CompressibleEulerModel 

53from dragonfly_sim.core.reconstruction import GreenGaussReconstruction 

54from dragonfly_sim.core.turbulence_models import SpalartAllmarasNeg 

55from dragonfly_sim.core.wall_distance import WallDistance 

56from dragonfly_sim.utils.Euler_utils import hll_flux, pressure 

57from dragonfly_sim.utils.NS_utils import (mean_flow, extra_components, concat_state, viscous_flux, 

58 viscous_stress, primitive_gradients, sutherland_viscosity) 

59 

60 

61class CompressibleRANSModel: 

62 

63 # Attributes that live on the RANS object itself; every other attribute 

64 # read or written goes to self.flow. 

65 _OWN = frozenset(( 

66 "flow", "turbulence_model", "Re", "Pr", "Pr_t", "viscosity_law", "sutherland_S_ratio", 

67 "muscl", "limiter", "limiter_eps", "jacobian_mode", "reconstruction", "wall_distance", 

68 "wall_tags", "mu_inf_const")) 

69 

70 def __init__(self, mesh, Re, turbulence_model=None, gamma=1.4, Pr=0.72, Pr_t=0.9, 

71 viscosity_law="sutherland", T_inf_dim=300.0, sutherland_S_dim=110.4, 

72 muscl=True, limiter=None, limiter_eps=1e-2, jacobian_mode="compact_pc"): 

73 tm = SpalartAllmarasNeg() if turbulence_model is None else turbulence_model 

74 if viscosity_law not in ("sutherland", "constant"): 

75 raise ValueError("viscosity_law must be 'sutherland' or 'constant'") 

76 self.turbulence_model = tm 

77 self.Re = float(Re) 

78 self.Pr = float(Pr) 

79 self.Pr_t = float(Pr_t) 

80 self.viscosity_law = viscosity_law 

81 self.sutherland_S_ratio = float(sutherland_S_dim) / float(T_inf_dim) 

82 self.muscl = muscl 

83 self.limiter = limiter 

84 self.limiter_eps = limiter_eps 

85 self.jacobian_mode = jacobian_mode 

86 self.reconstruction = None 

87 self.wall_distance = None 

88 self.wall_tags = [] 

89 self.mu_inf_const = None 

90 

91 self.flow = flow = CompressibleEulerModel(mesh, poly_order=0, gamma=gamma, 

92 flux_function=hll_flux, extra_variables=tm.names) 

93 # solver defaults (see the module docstring) 

94 flow.use_ptc = True 

95 flow.ptc_settings = {"cfl0": 5.0} 

96 flow.fom_newton_inner_max_it = 1 

97 flow.krylov_type = "fgmres" 

98 flow.asm_overlap = 1 

99 flow.ilu_levels = 2 

100 # physics plug-ins 

101 flow.freestream_scalars = self._freestream_extras 

102 flow.boundary_scalars = self._boundary_scalars 

103 flow.residual_terms.append(self._interior_terms) 

104 flow.boundary_terms.append(self._viscous_boundary_terms) 

105 flow.wall_bc = self.define_noslip_wall_bc 

106 flow.wall_traction = self._wall_traction 

107 flow.time_scale_terms.append(self._viscous_time_scale) 

108 flow.extra_scales = self._extra_scales 

109 

110 # ------------------------------------------------------------------ 

111 # Forwarding to the Euler host 

112 # ------------------------------------------------------------------ 

113 def __getattr__(self, name): 

114 # only called for attributes not found on the RANS object 

115 if name.startswith("__") or name in CompressibleRANSModel._OWN: 

116 raise AttributeError(name) 

117 return getattr(self.__dict__["flow"], name) 

118 

119 def __setattr__(self, name, value): 

120 if name in CompressibleRANSModel._OWN: 

121 object.__setattr__(self, name, value) 

122 else: 

123 setattr(self.flow, name, value) 

124 

125 # ------------------------------------------------------------------ 

126 # Set-up (overrides of the host's steps) 

127 # ------------------------------------------------------------------ 

128 def define_elements_functionspaces(self): 

129 self.flow.define_elements_functionspaces() 

130 if self.reconstruction is None: 

131 self.reconstruction = GreenGaussReconstruction( 

132 self.flow, muscl=self.muscl, limiter=self.limiter, limiter_eps=self.limiter_eps, 

133 limiter_scales=self._limiter_scales) 

134 self.reconstruction.install() 

135 

136 def define_farfield_conditions(self, bc_dict): 

137 self.flow.define_farfield_conditions(bc_dict) 

138 ff = self.flow.farfield_conditions 

139 self.mu_inf_const = dolfinx.fem.Constant( 

140 self.flow.mesh_obj.mesh, ff['rho'] * ff['M'] * ff['c'] / self.Re) 

141 

142 def set_mach(self, M): 

143 # Re is held fixed, so mu_inf follows U_inf 

144 self.flow.set_mach(M) 

145 ff = self.flow.farfield_conditions 

146 self.mu_inf_const.value = ff['rho'] * ff['M'] * ff['c'] / self.Re 

147 

148 def set_up_solver(self): 

149 flow = self.flow 

150 was_set_up = flow.solver_is_set_up 

151 flow.set_up_solver() 

152 if was_set_up: 

153 return 

154 mode = "exact" if flow.direct_linear_solver else self.jacobian_mode 

155 if mode != "exact": 

156 flow.problem.configure_snes(flow.solver.snes, mode) 

157 flow.solver.krylov_solver.setFromOptions() 

158 

159 # ------------------------------------------------------------------ 

160 # Material laws and free stream 

161 # ------------------------------------------------------------------ 

162 def laminar_viscosity(self, U): 

163 if self.viscosity_law == "constant": 

164 return self.mu_inf_const 

165 ff = self.flow.farfield_conditions 

166 return sutherland_viscosity(U, self.flow.gamma, self.flow.n_mean, self.mu_inf_const, 

167 ff['p'] / ff['rho'], self.sutherland_S_ratio) 

168 

169 def eddy_viscosity(self, U): 

170 return self.turbulence_model.eddy_viscosity(U, self.laminar_viscosity(U), self.flow.n_mean) 

171 

172 def viscous_flux(self, U, G, heat_flux=True): 

173 flow = self.flow 

174 mu = self.laminar_viscosity(U) 

175 tm = self.turbulence_model 

176 return viscous_flux(U, G, flow.gamma, flow.n_mean, mu, 

177 tm.eddy_viscosity(U, mu, flow.n_mean), self.Pr, self.Pr_t, 

178 extra_fluxes=tm.diffusive_fluxes(U, mu, flow.n_mean), 

179 heat_flux=heat_flux) 

180 

181 def freestream(self): 

182 """The free-stream dict TurbulenceModel methods take.""" 

183 ff = self.flow.farfield_conditions 

184 return {"rho": ff['rho'], "p": ff['p'], "c": ff['c'], "M": ff['M'], 

185 "speed": ff['M'] * ff['c'], "mu": self.mu_inf_const} 

186 

187 def _freestream_extras(self): 

188 return self.turbulence_model.freestream_values(self.freestream()) 

189 

190 def _extra_scales(self): 

191 return list(self.turbulence_model.limiter_scales(self.freestream())) 

192 

193 def _limiter_scales(self): 

194 return self.flow.state_scale_expressions() 

195 

196 def _boundary_scalars(self, kind, U, inflow): 

197 nm = self.flow.n_mean 

198 if kind == "slip": 

199 return extra_components(U, nm) 

200 if kind == "farfield": 

201 return [ufl.conditional(inflow, s_inf, s) 

202 for s_inf, s in zip(self._freestream_extras(), extra_components(U, nm))] 

203 if kind == "wall": 

204 return self.turbulence_model.wall_values(U, nm) 

205 raise ValueError("unknown boundary kind {!r}".format(kind)) 

206 

207 # ------------------------------------------------------------------ 

208 # Residual terms (plugged into the host) 

209 # ------------------------------------------------------------------ 

210 def _interior_terms(self): 

211 flow, rec = self.flow, self.reconstruction 

212 U, v, n = flow.u_vec, flow.v_vec, flow.mesh_obj.n 

213 nm = flow.n_mean 

214 jump_v = ufl.outer(v('+') - v('-'), n('+')) 

215 F = - ufl.inner(self.viscous_flux(0.5 * (U('+') + U('-')), rec.face_gradient_interior()), 

216 jump_v) * flow.forms.dS 

217 if flow.n_extra: 

218 # the WallDistance field, or a UFL expression standing in for it 

219 # (manufactured-solution tests) 

220 d = self.wall_distance 

221 d = 1.0 if d is None else getattr(d, "function", d) 

222 mu = self.laminar_viscosity(U) 

223 S = self.turbulence_model.sources(U, rec.G, mu, d, nm, flow.gamma) 

224 F -= sum(S[k] * v[nm + k] for k in range(flow.n_extra)) * flow.forms.dx 

225 return F 

226 

227 def _viscous_boundary_terms(self, U_b, meshtag, kind): 

228 if kind == "slip": 

229 return None # symmetry: zero shear stress and heat flux 

230 flow = self.flow 

231 Fv = self.viscous_flux(U_b, self.reconstruction.face_gradient_boundary(U_b), 

232 heat_flux=(kind != "wall")) 

233 return - ufl.inner(ufl.dot(Fv, flow.mesh_obj.n), flow.v_vec) * flow.forms.ds(meshtag) 

234 

235 def _viscous_time_scale(self, U): 

236 rho = ufl.max_value(U[0], 1e-8) 

237 nu_eff = (self.flow.gamma / self.Pr) * (self.laminar_viscosity(U) 

238 + ufl.max_value(self.eddy_viscosity(U), 0.0)) / rho 

239 h = self.flow.mesh_obj.h_min 

240 return 2.0 * nu_eff / h**2 

241 

242 # ------------------------------------------------------------------ 

243 # No-slip wall 

244 # ------------------------------------------------------------------ 

245 def noslip_wall_state(self, U): 

246 d = self.flow.dimensions 

247 rho = U[0] 

248 rho_e = U[1 + d] - 0.5 * sum(U[1 + i]**2 for i in range(d)) / rho 

249 return concat_state(ufl.as_vector([rho] + [0.0] * d + [rho_e]), 

250 self._boundary_scalars("wall", U, None)) 

251 

252 def define_noslip_wall_bc(self, meshtag): 

253 """ 

254 Adiabatic no-slip wall. The inviscid flux is the exact flux of the 

255 wall state, (0, p n, 0, 0): no mass, energy or scalar flux through 

256 the wall, the interior pressure on it. Conduction is dropped from the 

257 viscous flux, which is the q.n = 0 condition. 

258 """ 

259 flow = self.flow 

260 flow.add_boundary_condition(meshtag, self.noslip_wall_state(flow.u_vec), "wall", 

261 inviscid_flux="exact") 

262 self.wall_tags.append(meshtag) 

263 if self.turbulence_model.needs_wall_distance: 

264 if self.wall_distance is None: 

265 self.wall_distance = WallDistance(flow.mesh_obj, self.wall_tags) 

266 flow.derived_fields.append(self.wall_distance) 

267 else: 

268 self.wall_distance.wall_tags = list(self.wall_tags) 

269 self.wall_distance.update_geometry() 

270 

271 def define_wall_bc(self, meshtag): 

272 self.define_noslip_wall_bc(meshtag) 

273 

274 def _wall_stress_arguments(self, U): 

275 """(state, gradient) the wall shear is evaluated with: the wall state 

276 and the reconstructed facet gradient, exactly what the residual's wall 

277 flux uses.""" 

278 U_w = self.noslip_wall_state(U) 

279 return U_w, self.reconstruction.face_gradient_boundary(U_w, U) 

280 

281 def _wall_traction(self, U, n): 

282 U_s, G_s = self._wall_stress_arguments(U) 

283 grad_u, _, _, _ = primitive_gradients(U_s, G_s, self.flow.gamma, self.flow.n_mean) 

284 tau = viscous_stress(grad_u, self.laminar_viscosity(U_s) + self.eddy_viscosity(U_s)) 

285 return - ufl.dot(tau, n) 

286 

287 def skin_friction_vector(self, U, n, q_inf): 

288 """tau_w / q_inf: the wall shear (tangential traction) coefficient 

289 vector. Call update_derived_fields() before assembling it.""" 

290 U_s, G_s = self._wall_stress_arguments(U) 

291 grad_u, _, _, _ = primitive_gradients(U_s, G_s, self.flow.gamma, self.flow.n_mean) 

292 t = ufl.dot(viscous_stress(grad_u, self.laminar_viscosity(U_s)), n) 

293 return (t - ufl.dot(t, n) * n) / q_inf 

294 

295 def eddy_viscosity_ratio(self): 

296 """mu_t / mu of the current state, as a UFL expression.""" 

297 U = self.flow.u_vec 

298 return self.eddy_viscosity(U) / self.laminar_viscosity(U)