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 1109 : INSFVTKESourceSink::validParams()
19 : {
20 1109 : InputParameters params = FVElementalKernel::validParams();
21 1109 : params.addClassDescription("Elemental kernel to compute the production and destruction "
22 : " terms of turbulent kinetic energy (TKE).");
23 2218 : params.addRequiredParam<MooseFunctorName>("u", "The velocity in the x direction.");
24 2218 : params.addParam<MooseFunctorName>("v", "The velocity in the y direction.");
25 2218 : params.addParam<MooseFunctorName>("w", "The velocity in the z direction.");
26 1109 : params.addRequiredParam<MooseFunctorName>(NS::TKED,
27 : "Coupled turbulent kinetic energy dissipation rate.");
28 1109 : params.addRequiredParam<MooseFunctorName>(NS::density, "fluid density");
29 1109 : params.addRequiredParam<MooseFunctorName>(NS::mu, "Dynamic viscosity.");
30 1109 : params.addRequiredParam<MooseFunctorName>(NS::mu_t, "Turbulent viscosity.");
31 2218 : params.addParam<std::vector<BoundaryName>>(
32 : "walls", {}, "Boundaries that correspond to solid walls.");
33 2218 : params.addParam<bool>(
34 : "linearized_model",
35 2218 : true,
36 : "Boolean to determine if the problem should be use in a linear or nonlinear solve.");
37 2218 : MooseEnum wall_treatment("eq_newton eq_incremental eq_linearized neq", "neq");
38 2218 : 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 2218 : params.addParam<Real>("C_mu", 0.09, "Coupled turbulent kinetic energy closure.");
43 2218 : params.addParam<Real>("C_pl", 10.0, "Production Limiter Constant Multiplier.");
44 1109 : params.set<unsigned short>("ghost_layers") = 2;
45 2218 : params.addParam<bool>("newton_solve", false, "Whether a Newton nonlinear solve is being used");
46 2218 : params.addParamNamesToGroup("newton_solve", "Advanced");
47 :
48 1109 : return params;
49 1109 : }
50 :
51 257 : INSFVTKESourceSink::INSFVTKESourceSink(const InputParameters & params)
52 : : FVElementalKernel(params),
53 257 : _dim(_subproblem.mesh().dimension()),
54 514 : _u_var(getFunctor<ADReal>("u")),
55 771 : _v_var(params.isParamValid("v") ? &(getFunctor<ADReal>("v")) : nullptr),
56 257 : _w_var(params.isParamValid("w") ? &(getFunctor<ADReal>("w")) : nullptr),
57 257 : _epsilon(getFunctor<ADReal>(NS::TKED)),
58 257 : _rho(getFunctor<ADReal>(NS::density)),
59 257 : _mu(getFunctor<ADReal>(NS::mu)),
60 257 : _mu_t(getFunctor<ADReal>(NS::mu_t)),
61 514 : _wall_boundary_names(getParam<std::vector<BoundaryName>>("walls")),
62 514 : _linearized_model(getParam<bool>("linearized_model")),
63 514 : _wall_treatment(getParam<MooseEnum>("wall_treatment").getEnum<NS::WallTreatmentEnum>()),
64 514 : _C_mu(getParam<Real>("C_mu")),
65 514 : _C_pl(getParam<Real>("C_pl")),
66 771 : _newton_solve(getParam<bool>("newton_solve"))
67 : {
68 257 : if (_dim >= 2 && !_v_var)
69 0 : paramError("v", "In two or more dimensions, the v velocity must be supplied!");
70 :
71 257 : if (_dim >= 3 && !_w_var)
72 0 : paramError("w", "In three or more dimensions, the w velocity must be supplied!");
73 257 : }
74 :
75 : void
76 929 : INSFVTKESourceSink::initialSetup()
77 : {
78 1858 : NS::getWallBoundedElements(
79 929 : _wall_boundary_names, _fe_problem, _subproblem, blockIDs(), _wall_bounded);
80 929 : NS::getWallDistance(_wall_boundary_names, _fe_problem, _subproblem, blockIDs(), _dist);
81 929 : NS::getElementFaceArgs(_wall_boundary_names, _fe_problem, _subproblem, blockIDs(), _face_infos);
82 929 : }
83 :
84 : ADReal
85 351348 : INSFVTKESourceSink::computeQpResidual()
86 : {
87 : using std::max, std::sqrt, std::pow, std::min;
88 :
89 351348 : ADReal residual = 0.0;
90 351348 : ADReal production = 0.0;
91 351348 : ADReal destruction = 0.0;
92 :
93 351348 : const auto state = determineState();
94 351348 : const auto elem_arg = makeElemArg(_current_elem);
95 : const auto old_state =
96 351348 : _linearized_model ? Moose::StateArg(1, Moose::SolutionIterationType::Nonlinear) : state;
97 351348 : const auto rho = _rho(elem_arg, state);
98 351348 : const auto mu = _mu(elem_arg, state);
99 : // To prevent negative values & preserve sparsity pattern
100 1033064 : auto TKE = _newton_solve ? max(_var(elem_arg, old_state), ADReal(0) * _var(elem_arg, old_state))
101 351348 : : _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 351348 : if (_newton_solve)
105 82592 : TKE = max(TKE, 1e-10);
106 :
107 351348 : if (_wall_bounded.find(_current_elem) != _wall_bounded.end())
108 : {
109 : std::vector<ADReal> y_plus_vec, velocity_grad_norm_vec;
110 :
111 83194 : Real tot_weight = 0.0;
112 :
113 166388 : ADRealVectorValue velocity(_u_var(elem_arg, state));
114 83194 : if (_v_var)
115 83194 : velocity(1) = (*_v_var)(elem_arg, state);
116 83194 : if (_w_var)
117 0 : velocity(2) = (*_w_var)(elem_arg, state);
118 :
119 83194 : const auto & face_info_vec = libmesh_map_find(_face_infos, _current_elem);
120 83194 : const auto & distance_vec = libmesh_map_find(_dist, _current_elem);
121 :
122 172110 : for (unsigned int i = 0; i < distance_vec.size(); i++)
123 : {
124 : const auto parallel_speed = NS::computeSpeed<ADReal>(
125 266748 : velocity - velocity * face_info_vec[i]->normal() * face_info_vec[i]->normal());
126 88916 : const auto distance = distance_vec[i];
127 :
128 : ADReal y_plus;
129 88916 : if (_wall_treatment == NS::WallTreatmentEnum::NEQ) // Non-equilibrium / Non-iterative
130 8640 : y_plus = distance * sqrt(sqrt(_C_mu) * TKE) * rho / mu;
131 : else // Equilibrium / Iterative
132 86036 : y_plus = NS::findyPlus<ADReal>(mu, rho, max(parallel_speed, 1e-10), distance);
133 :
134 88916 : 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 88916 : velocity_grad_norm_vec.push_back(velocity_grad_norm);
155 :
156 88916 : tot_weight += 1.0;
157 : }
158 :
159 172110 : for (unsigned int i = 0; i < y_plus_vec.size(); i++)
160 : {
161 88916 : const auto y_plus = y_plus_vec[i];
162 :
163 88916 : const auto fi = face_info_vec[i];
164 88916 : const bool defined_on_elem_side = _var.hasFaceSide(*fi, true);
165 88916 : const Elem * const loc_elem = defined_on_elem_side ? &fi->elem() : fi->neighborPtr();
166 88916 : const Moose::FaceArg facearg = {
167 88916 : fi, Moose::FV::LimiterType::CentralDifference, false, false, loc_elem, nullptr};
168 88916 : const ADReal wall_mut = _mu_t(facearg, state);
169 88916 : const ADReal wall_mu = _mu(facearg, state);
170 :
171 177832 : const auto destruction_visc = 2.0 * wall_mu / Utility::pow<2>(distance_vec[i]) / tot_weight;
172 88916 : const auto destruction_log = pow(_C_mu, 0.75) * rho * pow(TKE, 0.5) /
173 177832 : (NS::von_karman_constant * distance_vec[i]) / tot_weight;
174 88916 : 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 88916 : if (y_plus < 11.25)
178 : {
179 61452 : destruction += destruction_visc;
180 61452 : if (_newton_solve)
181 85152 : destruction += 0 * destruction_log + 0 * tau_w;
182 : }
183 : else
184 : {
185 27464 : destruction += destruction_log;
186 27464 : if (_newton_solve)
187 0 : destruction += 0 * destruction_visc;
188 27464 : production += tau_w * pow(_C_mu, 0.25) / sqrt(TKE) /
189 54928 : (NS::von_karman_constant * distance_vec[i]) / tot_weight;
190 : }
191 : }
192 :
193 166388 : 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 83194 : if (_newton_solve)
196 51184 : residual += 0 * _epsilon(elem_arg, old_state);
197 83194 : }
198 : else
199 : {
200 268154 : const auto subdomain_id = _current_elem->subdomain_id();
201 268154 : const auto coord_sys = _subproblem.getCoordSystem(subdomain_id);
202 : const auto rz_radial_coord =
203 268154 : coord_sys == Moose::COORD_RZ ? _subproblem.getAxisymmetricRadialCoord() : 0;
204 : const auto symmetric_strain_tensor_sq_norm = NS::computeShearStrainRateNormSquared<ADReal>(
205 268154 : _u_var, _v_var, _w_var, elem_arg, state, coord_sys, rz_radial_coord);
206 :
207 536308 : production = _mu_t(elem_arg, state) * symmetric_strain_tensor_sq_norm;
208 :
209 268154 : const auto tke_old_raw = raw_value(TKE);
210 268154 : const auto epsilon_old = _epsilon(elem_arg, old_state);
211 :
212 268154 : if (MooseUtils::absoluteFuzzyEqual(tke_old_raw, 0))
213 0 : destruction = rho * epsilon_old;
214 : else
215 536308 : destruction = rho * _var(elem_arg, state) / tke_old_raw * raw_value(epsilon_old);
216 :
217 : // k-Production limiter (needed for flows with stagnation zones)
218 : const ADReal production_limit =
219 536308 : _C_pl * rho * (_newton_solve ? max(epsilon_old, ADReal(0)) : epsilon_old);
220 :
221 : // Apply production limiter
222 268154 : production = min(production, production_limit);
223 :
224 268154 : residual = destruction - production;
225 :
226 : // Additional 0-value terms to make sure new derivative entries are not added during the solve
227 268154 : if (_newton_solve)
228 114000 : residual += 0 * _epsilon(elem_arg, state);
229 : }
230 :
231 351348 : return residual;
232 : }
|
