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 2102 : INSFVTKESourceSink::validParams()
19 : {
20 2102 : InputParameters params = FVElementalKernel::validParams();
21 2102 : params.addClassDescription("Elemental kernel to compute the production and destruction "
22 : " terms of turbulent kinetic energy (TKE).");
23 4204 : params.addRequiredParam<MooseFunctorName>("u", "The velocity in the x direction.");
24 4204 : params.addParam<MooseFunctorName>("v", "The velocity in the y direction.");
25 4204 : params.addParam<MooseFunctorName>("w", "The velocity in the z direction.");
26 2102 : params.addRequiredParam<MooseFunctorName>(NS::TKED,
27 : "Coupled turbulent kinetic energy dissipation rate.");
28 2102 : params.addRequiredParam<MooseFunctorName>(NS::density, "fluid density");
29 2102 : params.addRequiredParam<MooseFunctorName>(NS::mu, "Dynamic viscosity.");
30 2102 : params.addRequiredParam<MooseFunctorName>(NS::mu_t, "Turbulent viscosity.");
31 4204 : params.addParam<std::vector<BoundaryName>>(
32 : "walls", {}, "Boundaries that correspond to solid walls.");
33 4204 : params.addParam<bool>(
34 : "linearized_model",
35 4204 : true,
36 : "Boolean to determine if the problem should be use in a linear or nonlinear solve.");
37 4204 : MooseEnum wall_treatment("eq_newton eq_incremental eq_linearized neq", "neq");
38 4204 : 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 4204 : params.addParam<Real>("C_mu", 0.09, "Coupled turbulent kinetic energy closure.");
43 4204 : params.addParam<Real>("C_pl", 10.0, "Production Limiter Constant Multiplier.");
44 2102 : params.set<unsigned short>("ghost_layers") = 2;
45 4204 : params.addParam<bool>("newton_solve", false, "Whether a Newton nonlinear solve is being used");
46 4204 : params.addParamNamesToGroup("newton_solve", "Advanced");
47 :
48 2102 : return params;
49 2102 : }
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 482560 : ADReal residual = 0.0;
88 482560 : ADReal production = 0.0;
89 482560 : ADReal destruction = 0.0;
90 :
91 482560 : const auto state = determineState();
92 482560 : const auto elem_arg = makeElemArg(_current_elem);
93 : const auto old_state =
94 482560 : _linearized_model ? Moose::StateArg(1, Moose::SolutionIterationType::Nonlinear) : state;
95 482560 : const auto rho = _rho(elem_arg, state);
96 482560 : const auto mu = _mu(elem_arg, state);
97 : // To prevent negative values & preserve sparsity pattern
98 482560 : auto TKE = _newton_solve
99 1304640 : ? std::max(_var(elem_arg, old_state), ADReal(0) * _var(elem_arg, old_state))
100 482560 : : _var(elem_arg, old_state);
101 : // Prevent computation of sqrt(0) with undefined automatic derivatives
102 : // This is not needed for segregated solves, as TKE has minimum bound in the solver
103 482560 : if (_newton_solve)
104 84880 : TKE = std::max(TKE, 1e-10);
105 :
106 482560 : if (_wall_bounded.find(_current_elem) != _wall_bounded.end())
107 : {
108 : std::vector<ADReal> y_plus_vec, velocity_grad_norm_vec;
109 :
110 112000 : Real tot_weight = 0.0;
111 :
112 224000 : ADRealVectorValue velocity(_u_var(elem_arg, state));
113 112000 : if (_v_var)
114 112000 : velocity(1) = (*_v_var)(elem_arg, state);
115 112000 : if (_w_var)
116 0 : velocity(2) = (*_w_var)(elem_arg, state);
117 :
118 112000 : const auto & face_info_vec = libmesh_map_find(_face_infos, _current_elem);
119 112000 : const auto & distance_vec = libmesh_map_find(_dist, _current_elem);
120 :
121 231360 : for (unsigned int i = 0; i < distance_vec.size(); i++)
122 : {
123 : const auto parallel_speed = NS::computeSpeed<ADReal>(
124 358080 : velocity - velocity * face_info_vec[i]->normal() * face_info_vec[i]->normal());
125 119360 : const auto distance = distance_vec[i];
126 :
127 : ADReal y_plus;
128 119360 : if (_wall_treatment == NS::WallTreatmentEnum::NEQ) // Non-equilibrium / Non-iterative
129 12960 : y_plus = distance * std::sqrt(std::sqrt(_C_mu) * TKE) * rho / mu;
130 : else // Equilibrium / Iterative
131 115040 : y_plus = NS::findyPlus<ADReal>(mu, rho, std::max(parallel_speed, 1e-10), distance);
132 :
133 119360 : y_plus_vec.push_back(y_plus);
134 :
135 : const ADReal velocity_grad_norm = parallel_speed / distance;
136 :
137 : /// Do not erase!!
138 : // More complete expansion for velocity gradient. Leave commented for now.
139 : // Will be useful later when doing two-phase or compressible flow
140 : // ADReal velocity_grad_norm_sq =
141 : // Utility::pow<2>(_u_var->gradient(elem_arg, state) *
142 : // _normal[_current_elem][i]);
143 : // if (_dim >= 2)
144 : // velocity_grad_norm_sq +=
145 : // Utility::pow<2>(_v_var->gradient(elem_arg, state) *
146 : // _normal[_current_elem][i]);
147 : // if (_dim >= 3)
148 : // velocity_grad_norm_sq +=
149 : // Utility::pow<2>(_w_var->gradient(elem_arg, state) *
150 : // _normal[_current_elem][i]);
151 : // ADReal velocity_grad_norm = std::sqrt(velocity_grad_norm_sq);
152 :
153 119360 : velocity_grad_norm_vec.push_back(velocity_grad_norm);
154 :
155 119360 : tot_weight += 1.0;
156 : }
157 :
158 231360 : for (unsigned int i = 0; i < y_plus_vec.size(); i++)
159 : {
160 119360 : const auto y_plus = y_plus_vec[i];
161 :
162 119360 : const auto fi = face_info_vec[i];
163 119360 : const bool defined_on_elem_side = _var.hasFaceSide(*fi, true);
164 119360 : const Elem * const loc_elem = defined_on_elem_side ? &fi->elem() : fi->neighborPtr();
165 119360 : const Moose::FaceArg facearg = {
166 119360 : fi, Moose::FV::LimiterType::CentralDifference, false, false, loc_elem, nullptr};
167 119360 : const ADReal wall_mut = _mu_t(facearg, state);
168 119360 : const ADReal wall_mu = _mu(facearg, state);
169 :
170 238720 : const auto destruction_visc = 2.0 * wall_mu / Utility::pow<2>(distance_vec[i]) / tot_weight;
171 119360 : const auto destruction_log = std::pow(_C_mu, 0.75) * rho * std::pow(TKE, 0.5) /
172 238720 : (NS::von_karman_constant * distance_vec[i]) / tot_weight;
173 119360 : const auto tau_w = (wall_mut + wall_mu) * velocity_grad_norm_vec[i];
174 :
175 : // Additional 0-value terms to make sure new derivative entries are not added during the solve
176 119360 : if (y_plus < 11.25)
177 : {
178 78956 : destruction += destruction_visc;
179 78956 : if (_newton_solve)
180 88080 : destruction += 0 * destruction_log + 0 * tau_w;
181 : }
182 : else
183 : {
184 40404 : destruction += destruction_log;
185 40404 : if (_newton_solve)
186 0 : destruction += 0 * destruction_visc;
187 40404 : production += tau_w * std::pow(_C_mu, 0.25) / std::sqrt(TKE) /
188 80808 : (NS::von_karman_constant * distance_vec[i]) / tot_weight;
189 : }
190 : }
191 :
192 224000 : residual = (destruction - production) * _var(elem_arg, state);
193 : // Additional 0-value term to make sure new derivative entries are not added during the solve
194 112000 : if (_newton_solve)
195 52760 : residual += 0 * _epsilon(elem_arg, old_state);
196 : }
197 : else
198 : {
199 370560 : const auto & grad_u = _u_var.gradient(elem_arg, state);
200 370560 : const auto Sij_xx = 2.0 * grad_u(0);
201 370560 : ADReal Sij_xy = 0.0;
202 370560 : ADReal Sij_xz = 0.0;
203 370560 : ADReal Sij_yy = 0.0;
204 370560 : ADReal Sij_yz = 0.0;
205 370560 : ADReal Sij_zz = 0.0;
206 :
207 370560 : const auto grad_xx = grad_u(0);
208 370560 : ADReal grad_xy = 0.0;
209 370560 : ADReal grad_xz = 0.0;
210 370560 : ADReal grad_yx = 0.0;
211 370560 : ADReal grad_yy = 0.0;
212 370560 : ADReal grad_yz = 0.0;
213 370560 : ADReal grad_zx = 0.0;
214 370560 : ADReal grad_zy = 0.0;
215 370560 : ADReal grad_zz = 0.0;
216 :
217 370560 : auto trace = Sij_xx / 3.0;
218 :
219 370560 : if (_dim >= 2)
220 : {
221 370560 : const auto & grad_v = (*_v_var).gradient(elem_arg, state);
222 370560 : Sij_xy = grad_u(1) + grad_v(0);
223 741120 : Sij_yy = 2.0 * grad_v(1);
224 :
225 370560 : grad_xy = grad_u(1);
226 370560 : grad_yx = grad_v(0);
227 370560 : grad_yy = grad_v(1);
228 :
229 741120 : trace += Sij_yy / 3.0;
230 :
231 370560 : if (_dim >= 3)
232 : {
233 0 : const auto & grad_w = (*_w_var).gradient(elem_arg, state);
234 :
235 0 : Sij_xz = grad_u(2) + grad_w(0);
236 0 : Sij_yz = grad_v(2) + grad_w(1);
237 0 : Sij_zz = 2.0 * grad_w(2);
238 :
239 0 : grad_xz = grad_u(2);
240 0 : grad_yz = grad_v(2);
241 0 : grad_zx = grad_w(0);
242 0 : grad_zy = grad_w(1);
243 0 : grad_zz = grad_w(2);
244 :
245 0 : trace += Sij_zz / 3.0;
246 : }
247 : }
248 :
249 : const auto symmetric_strain_tensor_sq_norm =
250 370560 : (Sij_xx - trace) * grad_xx + Sij_xy * grad_xy + Sij_xz * grad_xz + Sij_xy * grad_yx +
251 370560 : (Sij_yy - trace) * grad_yy + Sij_yz * grad_yz + Sij_xz * grad_zx + Sij_yz * grad_zy +
252 370560 : (Sij_zz - trace) * grad_zz;
253 :
254 741120 : production = _mu_t(elem_arg, state) * symmetric_strain_tensor_sq_norm;
255 :
256 370560 : const auto tke_old_raw = raw_value(TKE);
257 370560 : const auto epsilon_old = _epsilon(elem_arg, old_state);
258 :
259 370560 : if (MooseUtils::absoluteFuzzyEqual(tke_old_raw, 0))
260 0 : destruction = rho * epsilon_old;
261 : else
262 741120 : destruction = rho * _var(elem_arg, state) / tke_old_raw * raw_value(epsilon_old);
263 :
264 : // k-Production limiter (needed for flows with stagnation zones)
265 : const ADReal production_limit =
266 741120 : _C_pl * rho * (_newton_solve ? std::max(epsilon_old, ADReal(0)) : epsilon_old);
267 :
268 : // Apply production limiter
269 370560 : production = std::min(production, production_limit);
270 :
271 370560 : residual = destruction - production;
272 :
273 : // Additional 0-value terms to make sure new derivative entries are not added during the solve
274 370560 : if (_newton_solve)
275 117000 : residual += 0 * _epsilon(elem_arg, state);
276 : }
277 :
278 482560 : return residual;
279 : }
|
