Line data Source code
1 : //* This file is part of the MOOSE framework
2 : //* https://mooseframework.inl.gov
3 : //*
4 : //* All rights reserved, see COPYRIGHT for full restrictions
5 : //* https://github.com/idaholab/moose/blob/master/COPYRIGHT
6 : //*
7 : //* Licensed under LGPL 2.1, please see LICENSE for details
8 : //* https://www.gnu.org/licenses/lgpl-2.1.html
9 :
10 : #include "INSFVTKESourceSink.h"
11 : #include "NonlinearSystemBase.h"
12 : #include "NavierStokesMethods.h"
13 : #include "libmesh/nonlinear_solver.h"
14 :
15 : registerMooseObject("NavierStokesApp", INSFVTKESourceSink);
16 :
17 : InputParameters
18 2121 : INSFVTKESourceSink::validParams()
19 : {
20 2121 : InputParameters params = FVElementalKernel::validParams();
21 2121 : params.addClassDescription("Elemental kernel to compute the production and destruction "
22 : " terms of turbulent kinetic energy (TKE).");
23 4242 : params.addRequiredParam<MooseFunctorName>("u", "The velocity in the x direction.");
24 4242 : params.addParam<MooseFunctorName>("v", "The velocity in the y direction.");
25 4242 : params.addParam<MooseFunctorName>("w", "The velocity in the z direction.");
26 2121 : params.addRequiredParam<MooseFunctorName>(NS::TKED,
27 : "Coupled turbulent kinetic energy dissipation rate.");
28 2121 : params.addRequiredParam<MooseFunctorName>(NS::density, "fluid density");
29 2121 : params.addRequiredParam<MooseFunctorName>(NS::mu, "Dynamic viscosity.");
30 2121 : params.addRequiredParam<MooseFunctorName>(NS::mu_t, "Turbulent viscosity.");
31 4242 : params.addParam<std::vector<BoundaryName>>(
32 : "walls", {}, "Boundaries that correspond to solid walls.");
33 4242 : params.addParam<bool>(
34 : "linearized_model",
35 4242 : true,
36 : "Boolean to determine if the problem should be use in a linear or nonlinear solve.");
37 4242 : MooseEnum wall_treatment("eq_newton eq_incremental eq_linearized neq", "neq");
38 4242 : params.addParam<MooseEnum>("wall_treatment",
39 : wall_treatment,
40 : "The method used for computing the wall functions "
41 : "'eq_newton', 'eq_incremental', 'eq_linearized', 'neq'");
42 4242 : params.addParam<Real>("C_mu", 0.09, "Coupled turbulent kinetic energy closure.");
43 4242 : params.addParam<Real>("C_pl", 10.0, "Production Limiter Constant Multiplier.");
44 2121 : params.set<unsigned short>("ghost_layers") = 2;
45 4242 : params.addParam<bool>("newton_solve", false, "Whether a Newton nonlinear solve is being used");
46 4242 : params.addParamNamesToGroup("newton_solve", "Advanced");
47 :
48 2121 : return params;
49 2121 : }
50 :
51 541 : INSFVTKESourceSink::INSFVTKESourceSink(const InputParameters & params)
52 : : FVElementalKernel(params),
53 541 : _dim(_subproblem.mesh().dimension()),
54 1082 : _u_var(getFunctor<ADReal>("u")),
55 1623 : _v_var(params.isParamValid("v") ? &(getFunctor<ADReal>("v")) : nullptr),
56 541 : _w_var(params.isParamValid("w") ? &(getFunctor<ADReal>("w")) : nullptr),
57 541 : _epsilon(getFunctor<ADReal>(NS::TKED)),
58 541 : _rho(getFunctor<ADReal>(NS::density)),
59 541 : _mu(getFunctor<ADReal>(NS::mu)),
60 541 : _mu_t(getFunctor<ADReal>(NS::mu_t)),
61 1082 : _wall_boundary_names(getParam<std::vector<BoundaryName>>("walls")),
62 1082 : _linearized_model(getParam<bool>("linearized_model")),
63 1082 : _wall_treatment(getParam<MooseEnum>("wall_treatment").getEnum<NS::WallTreatmentEnum>()),
64 1082 : _C_mu(getParam<Real>("C_mu")),
65 1082 : _C_pl(getParam<Real>("C_pl")),
66 1623 : _newton_solve(getParam<bool>("newton_solve"))
67 : {
68 541 : if (_dim >= 2 && !_v_var)
69 0 : paramError("v", "In two or more dimensions, the v velocity must be supplied!");
70 :
71 541 : if (_dim >= 3 && !_w_var)
72 0 : paramError("w", "In three or more dimensions, the w velocity must be supplied!");
73 541 : }
74 :
75 : void
76 1997 : INSFVTKESourceSink::initialSetup()
77 : {
78 3994 : NS::getWallBoundedElements(
79 1997 : _wall_boundary_names, _fe_problem, _subproblem, blockIDs(), _wall_bounded);
80 1997 : NS::getWallDistance(_wall_boundary_names, _fe_problem, _subproblem, blockIDs(), _dist);
81 1997 : NS::getElementFaceArgs(_wall_boundary_names, _fe_problem, _subproblem, blockIDs(), _face_infos);
82 1997 : }
83 :
84 : ADReal
85 482560 : INSFVTKESourceSink::computeQpResidual()
86 : {
87 : using std::max, std::sqrt, std::pow, std::min;
88 :
89 482560 : ADReal residual = 0.0;
90 482560 : ADReal production = 0.0;
91 482560 : ADReal destruction = 0.0;
92 :
93 482560 : const auto state = determineState();
94 482560 : const auto elem_arg = makeElemArg(_current_elem);
95 : const auto old_state =
96 482560 : _linearized_model ? Moose::StateArg(1, Moose::SolutionIterationType::Nonlinear) : state;
97 482560 : const auto rho = _rho(elem_arg, state);
98 482560 : const auto mu = _mu(elem_arg, state);
99 : // To prevent negative values & preserve sparsity pattern
100 1304640 : auto TKE = _newton_solve ? max(_var(elem_arg, old_state), ADReal(0) * _var(elem_arg, old_state))
101 482560 : : _var(elem_arg, old_state);
102 : // Prevent computation of sqrt(0) with undefined automatic derivatives
103 : // This is not needed for segregated solves, as TKE has minimum bound in the solver
104 482560 : if (_newton_solve)
105 84880 : TKE = max(TKE, 1e-10);
106 :
107 482560 : if (_wall_bounded.find(_current_elem) != _wall_bounded.end())
108 : {
109 : std::vector<ADReal> y_plus_vec, velocity_grad_norm_vec;
110 :
111 112000 : Real tot_weight = 0.0;
112 :
113 224000 : ADRealVectorValue velocity(_u_var(elem_arg, state));
114 112000 : if (_v_var)
115 112000 : velocity(1) = (*_v_var)(elem_arg, state);
116 112000 : if (_w_var)
117 0 : velocity(2) = (*_w_var)(elem_arg, state);
118 :
119 112000 : const auto & face_info_vec = libmesh_map_find(_face_infos, _current_elem);
120 112000 : const auto & distance_vec = libmesh_map_find(_dist, _current_elem);
121 :
122 231360 : for (unsigned int i = 0; i < distance_vec.size(); i++)
123 : {
124 : const auto parallel_speed = NS::computeSpeed<ADReal>(
125 358080 : velocity - velocity * face_info_vec[i]->normal() * face_info_vec[i]->normal());
126 119360 : const auto distance = distance_vec[i];
127 :
128 : ADReal y_plus;
129 119360 : if (_wall_treatment == NS::WallTreatmentEnum::NEQ) // Non-equilibrium / Non-iterative
130 12960 : y_plus = distance * sqrt(sqrt(_C_mu) * TKE) * rho / mu;
131 : else // Equilibrium / Iterative
132 115040 : y_plus = NS::findyPlus<ADReal>(mu, rho, max(parallel_speed, 1e-10), distance);
133 :
134 119360 : y_plus_vec.push_back(y_plus);
135 :
136 : const ADReal velocity_grad_norm = parallel_speed / distance;
137 :
138 : /// Do not erase!!
139 : // More complete expansion for velocity gradient. Leave commented for now.
140 : // Will be useful later when doing two-phase or compressible flow
141 : // ADReal velocity_grad_norm_sq =
142 : // Utility::pow<2>(_u_var->gradient(elem_arg, state) *
143 : // _normal[_current_elem][i]);
144 : // if (_dim >= 2)
145 : // velocity_grad_norm_sq +=
146 : // Utility::pow<2>(_v_var->gradient(elem_arg, state) *
147 : // _normal[_current_elem][i]);
148 : // if (_dim >= 3)
149 : // velocity_grad_norm_sq +=
150 : // Utility::pow<2>(_w_var->gradient(elem_arg, state) *
151 : // _normal[_current_elem][i]);
152 : // ADReal velocity_grad_norm = sqrt(velocity_grad_norm_sq);
153 :
154 119360 : velocity_grad_norm_vec.push_back(velocity_grad_norm);
155 :
156 119360 : tot_weight += 1.0;
157 : }
158 :
159 231360 : for (unsigned int i = 0; i < y_plus_vec.size(); i++)
160 : {
161 119360 : const auto y_plus = y_plus_vec[i];
162 :
163 119360 : const auto fi = face_info_vec[i];
164 119360 : const bool defined_on_elem_side = _var.hasFaceSide(*fi, true);
165 119360 : const Elem * const loc_elem = defined_on_elem_side ? &fi->elem() : fi->neighborPtr();
166 119360 : const Moose::FaceArg facearg = {
167 119360 : fi, Moose::FV::LimiterType::CentralDifference, false, false, loc_elem, nullptr};
168 119360 : const ADReal wall_mut = _mu_t(facearg, state);
169 119360 : const ADReal wall_mu = _mu(facearg, state);
170 :
171 238720 : const auto destruction_visc = 2.0 * wall_mu / Utility::pow<2>(distance_vec[i]) / tot_weight;
172 119360 : const auto destruction_log = pow(_C_mu, 0.75) * rho * pow(TKE, 0.5) /
173 238720 : (NS::von_karman_constant * distance_vec[i]) / tot_weight;
174 119360 : const auto tau_w = (wall_mut + wall_mu) * velocity_grad_norm_vec[i];
175 :
176 : // Additional 0-value terms to make sure new derivative entries are not added during the solve
177 119360 : if (y_plus < 11.25)
178 : {
179 78956 : destruction += destruction_visc;
180 78956 : if (_newton_solve)
181 88080 : destruction += 0 * destruction_log + 0 * tau_w;
182 : }
183 : else
184 : {
185 40404 : destruction += destruction_log;
186 40404 : if (_newton_solve)
187 0 : destruction += 0 * destruction_visc;
188 40404 : production += tau_w * pow(_C_mu, 0.25) / sqrt(TKE) /
189 80808 : (NS::von_karman_constant * distance_vec[i]) / tot_weight;
190 : }
191 : }
192 :
193 224000 : residual = (destruction - production) * _var(elem_arg, state);
194 : // Additional 0-value term to make sure new derivative entries are not added during the solve
195 112000 : if (_newton_solve)
196 52760 : residual += 0 * _epsilon(elem_arg, old_state);
197 112000 : }
198 : else
199 : {
200 370560 : const auto & grad_u = _u_var.gradient(elem_arg, state);
201 370560 : const auto Sij_xx = 2.0 * grad_u(0);
202 370560 : ADReal Sij_xy = 0.0;
203 370560 : ADReal Sij_xz = 0.0;
204 370560 : ADReal Sij_yy = 0.0;
205 370560 : ADReal Sij_yz = 0.0;
206 370560 : ADReal Sij_zz = 0.0;
207 :
208 370560 : const auto grad_xx = grad_u(0);
209 370560 : ADReal grad_xy = 0.0;
210 370560 : ADReal grad_xz = 0.0;
211 370560 : ADReal grad_yx = 0.0;
212 370560 : ADReal grad_yy = 0.0;
213 370560 : ADReal grad_yz = 0.0;
214 370560 : ADReal grad_zx = 0.0;
215 370560 : ADReal grad_zy = 0.0;
216 370560 : ADReal grad_zz = 0.0;
217 :
218 370560 : auto trace = Sij_xx / 3.0;
219 :
220 370560 : if (_dim >= 2)
221 : {
222 370560 : const auto & grad_v = (*_v_var).gradient(elem_arg, state);
223 370560 : Sij_xy = grad_u(1) + grad_v(0);
224 741120 : Sij_yy = 2.0 * grad_v(1);
225 :
226 370560 : grad_xy = grad_u(1);
227 370560 : grad_yx = grad_v(0);
228 370560 : grad_yy = grad_v(1);
229 :
230 741120 : trace += Sij_yy / 3.0;
231 :
232 370560 : if (_dim >= 3)
233 : {
234 0 : const auto & grad_w = (*_w_var).gradient(elem_arg, state);
235 :
236 0 : Sij_xz = grad_u(2) + grad_w(0);
237 0 : Sij_yz = grad_v(2) + grad_w(1);
238 0 : Sij_zz = 2.0 * grad_w(2);
239 :
240 0 : grad_xz = grad_u(2);
241 0 : grad_yz = grad_v(2);
242 0 : grad_zx = grad_w(0);
243 0 : grad_zy = grad_w(1);
244 0 : grad_zz = grad_w(2);
245 :
246 0 : trace += Sij_zz / 3.0;
247 : }
248 : }
249 :
250 : const auto symmetric_strain_tensor_sq_norm =
251 370560 : (Sij_xx - trace) * grad_xx + Sij_xy * grad_xy + Sij_xz * grad_xz + Sij_xy * grad_yx +
252 370560 : (Sij_yy - trace) * grad_yy + Sij_yz * grad_yz + Sij_xz * grad_zx + Sij_yz * grad_zy +
253 370560 : (Sij_zz - trace) * grad_zz;
254 :
255 741120 : production = _mu_t(elem_arg, state) * symmetric_strain_tensor_sq_norm;
256 :
257 370560 : const auto tke_old_raw = raw_value(TKE);
258 370560 : const auto epsilon_old = _epsilon(elem_arg, old_state);
259 :
260 370560 : if (MooseUtils::absoluteFuzzyEqual(tke_old_raw, 0))
261 0 : destruction = rho * epsilon_old;
262 : else
263 741120 : destruction = rho * _var(elem_arg, state) / tke_old_raw * raw_value(epsilon_old);
264 :
265 : // k-Production limiter (needed for flows with stagnation zones)
266 : const ADReal production_limit =
267 741120 : _C_pl * rho * (_newton_solve ? max(epsilon_old, ADReal(0)) : epsilon_old);
268 :
269 : // Apply production limiter
270 370560 : production = min(production, production_limit);
271 :
272 370560 : residual = destruction - production;
273 :
274 : // Additional 0-value terms to make sure new derivative entries are not added during the solve
275 370560 : if (_newton_solve)
276 117000 : residual += 0 * _epsilon(elem_arg, state);
277 : }
278 :
279 482560 : return residual;
280 : }
|
