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 "INSFVTKEDSourceSink.h"
11 : #include "NonlinearSystemBase.h"
12 : #include "NavierStokesMethods.h"
13 : #include "libmesh/nonlinear_solver.h"
14 :
15 : registerMooseObject("NavierStokesApp", INSFVTKEDSourceSink);
16 :
17 : InputParameters
18 1005 : INSFVTKEDSourceSink::validParams()
19 : {
20 1005 : InputParameters params = FVElementalKernel::validParams();
21 1005 : params.addClassDescription("Elemental kernel to compute the production and destruction "
22 : " terms of turbulent kinetic energy dissipation (TKED).");
23 2010 : params.addRequiredParam<MooseFunctorName>("u", "The velocity in the x direction.");
24 2010 : params.addParam<MooseFunctorName>("v", "The velocity in the y direction.");
25 2010 : params.addParam<MooseFunctorName>("w", "The velocity in the z direction.");
26 2010 : params.addRequiredParam<MooseFunctorName>("k", "Coupled turbulent kinetic energy.");
27 2010 : params.deprecateParam("k", NS::TKE, "01/01/2025");
28 1005 : params.addRequiredParam<MooseFunctorName>(NS::density, "fluid density");
29 1005 : params.addRequiredParam<MooseFunctorName>(NS::mu, "Dynamic viscosity.");
30 1005 : params.addRequiredParam<MooseFunctorName>(NS::mu_t, "Turbulent viscosity.");
31 2010 : params.addParam<std::vector<BoundaryName>>(
32 : "walls", {}, "Boundaries that correspond to solid walls.");
33 2010 : params.addParam<bool>(
34 : "linearized_model",
35 2010 : true,
36 : "Boolean to determine if the problem should be used in a linear or nonlinear solve");
37 2010 : MooseEnum wall_treatment("eq_newton eq_incremental eq_linearized neq", "neq");
38 2010 : 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 2010 : params.addParam<Real>("C1_eps", 1.44, "First epsilon coefficient");
43 2010 : params.addParam<Real>("C2_eps", 1.92, "Second epsilon coefficient");
44 2010 : params.addParam<Real>("C_mu", 0.09, "Coupled turbulent kinetic energy closure.");
45 2010 : params.addParam<Real>("C_pl", 10.0, "Production limiter constant multiplier.");
46 1005 : params.set<unsigned short>("ghost_layers") = 2;
47 2010 : params.addParam<bool>("newton_solve", false, "Whether a Newton nonlinear solve is being used");
48 2010 : params.addParamNamesToGroup("newton_solve", "Advanced");
49 1005 : return params;
50 1005 : }
51 :
52 541 : INSFVTKEDSourceSink::INSFVTKEDSourceSink(const InputParameters & params)
53 : : FVElementalKernel(params),
54 541 : _dim(_subproblem.mesh().dimension()),
55 1082 : _u_var(getFunctor<ADReal>("u")),
56 1623 : _v_var(params.isParamValid("v") ? &(getFunctor<ADReal>("v")) : nullptr),
57 541 : _w_var(params.isParamValid("w") ? &(getFunctor<ADReal>("w")) : nullptr),
58 541 : _k(getFunctor<ADReal>(NS::TKE)),
59 541 : _rho(getFunctor<ADReal>(NS::density)),
60 541 : _mu(getFunctor<ADReal>(NS::mu)),
61 541 : _mu_t(getFunctor<ADReal>(NS::mu_t)),
62 1082 : _wall_boundary_names(getParam<std::vector<BoundaryName>>("walls")),
63 1082 : _linearized_model(getParam<bool>("linearized_model")),
64 1082 : _wall_treatment(getParam<MooseEnum>("wall_treatment").getEnum<NS::WallTreatmentEnum>()),
65 1082 : _C1_eps(getParam<Real>("C1_eps")),
66 1082 : _C2_eps(getParam<Real>("C2_eps")),
67 1082 : _C_mu(getParam<Real>("C_mu")),
68 1082 : _C_pl(getParam<Real>("C_pl")),
69 1623 : _newton_solve(getParam<bool>("newton_solve"))
70 : {
71 541 : if (_dim >= 2 && !_v_var)
72 0 : paramError("v", "In two or more dimensions, the v velocity must be supplied!");
73 :
74 541 : if (_dim >= 3 && !_w_var)
75 0 : paramError("w", "In three or more dimensions, the w velocity must be supplied!");
76 541 : }
77 :
78 : void
79 1997 : INSFVTKEDSourceSink::initialSetup()
80 : {
81 3994 : NS::getWallBoundedElements(
82 1997 : _wall_boundary_names, _fe_problem, _subproblem, blockIDs(), _wall_bounded);
83 1997 : NS::getWallDistance(_wall_boundary_names, _fe_problem, _subproblem, blockIDs(), _dist);
84 1997 : NS::getElementFaceArgs(_wall_boundary_names, _fe_problem, _subproblem, blockIDs(), _face_infos);
85 1997 : }
86 :
87 : ADReal
88 482584 : INSFVTKEDSourceSink::computeQpResidual()
89 : {
90 482584 : ADReal residual = 0.0;
91 482584 : ADReal production = 0.0;
92 482584 : ADReal destruction = 0.0;
93 482584 : const auto elem_arg = makeElemArg(_current_elem);
94 482584 : const auto state = determineState();
95 : const auto old_state =
96 482584 : _linearized_model ? Moose::StateArg(1, Moose::SolutionIterationType::Nonlinear) : state;
97 482584 : const auto mu = _mu(elem_arg, state);
98 482584 : const auto rho = _rho(elem_arg, state);
99 : const auto TKE_old =
100 567464 : _newton_solve ? std::max(_k(elem_arg, old_state), 1e-10) : _k(elem_arg, old_state);
101 : ADReal y_plus;
102 :
103 482584 : if (_wall_bounded.find(_current_elem) != _wall_bounded.end())
104 : {
105 : std::vector<ADReal> y_plus_vec;
106 :
107 112008 : Real tot_weight = 0.0;
108 :
109 224016 : ADRealVectorValue velocity(_u_var(elem_arg, state));
110 112008 : if (_v_var)
111 112008 : velocity(1) = (*_v_var)(elem_arg, state);
112 112008 : if (_w_var)
113 0 : velocity(2) = (*_w_var)(elem_arg, state);
114 :
115 112008 : const auto & face_info_vec = libmesh_map_find(_face_infos, _current_elem);
116 112008 : const auto & distance_vec = libmesh_map_find(_dist, _current_elem);
117 : mooseAssert(distance_vec.size(), "Should have found a distance vector");
118 : mooseAssert(distance_vec.size() == face_info_vec.size(),
119 : "Should be as many distance vectors as face info vectors");
120 :
121 231376 : for (unsigned int i = 0; i < distance_vec.size(); i++)
122 : {
123 119368 : const auto distance = distance_vec[i];
124 : mooseAssert(distance > 0, "Should be at a non-zero distance");
125 :
126 119368 : if (_wall_treatment == NS::WallTreatmentEnum::NEQ) // Non-equilibrium / Non-iterative
127 12960 : y_plus = distance * std::sqrt(std::sqrt(_C_mu) * TKE_old) * rho / mu;
128 : else
129 : {
130 : // Equilibrium / Iterative
131 : const auto parallel_speed = NS::computeSpeed<ADReal>(
132 345144 : velocity - velocity * face_info_vec[i]->normal() * face_info_vec[i]->normal());
133 :
134 115048 : y_plus = NS::findyPlus<ADReal>(mu, rho, std::max(parallel_speed, 1e-10), distance);
135 : }
136 :
137 119368 : y_plus_vec.push_back(y_plus);
138 :
139 119368 : tot_weight += 1.0;
140 : }
141 :
142 231376 : for (const auto i : index_range(y_plus_vec))
143 : {
144 119368 : const auto y_plus = y_plus_vec[i];
145 :
146 119368 : if (y_plus < 11.25)
147 : {
148 78956 : const auto fi = face_info_vec[i];
149 78956 : const bool defined_on_elem_side = _var.hasFaceSide(*fi, true);
150 78956 : const Elem * const loc_elem = defined_on_elem_side ? &fi->elem() : fi->neighborPtr();
151 78956 : const Moose::FaceArg facearg = {
152 78956 : fi, Moose::FV::LimiterType::CentralDifference, false, false, loc_elem, nullptr};
153 78956 : destruction += 2.0 * TKE_old * _mu(facearg, state) / rho /
154 157912 : Utility::pow<2>(distance_vec[i]) / tot_weight;
155 : }
156 : else
157 40412 : destruction += std::pow(_C_mu, 0.75) * std::pow(TKE_old, 1.5) /
158 80824 : (NS::von_karman_constant * distance_vec[i]) / tot_weight;
159 : }
160 :
161 224016 : residual = _var(makeElemArg(_current_elem), state) - destruction;
162 112008 : }
163 : else
164 : {
165 370576 : const auto & grad_u = _u_var.gradient(elem_arg, state);
166 370576 : const auto Sij_xx = 2.0 * grad_u(0);
167 370576 : ADReal Sij_xy = 0.0;
168 370576 : ADReal Sij_xz = 0.0;
169 370576 : ADReal Sij_yy = 0.0;
170 370576 : ADReal Sij_yz = 0.0;
171 370576 : ADReal Sij_zz = 0.0;
172 :
173 370576 : const auto grad_xx = grad_u(0);
174 370576 : ADReal grad_xy = 0.0;
175 370576 : ADReal grad_xz = 0.0;
176 370576 : ADReal grad_yx = 0.0;
177 370576 : ADReal grad_yy = 0.0;
178 370576 : ADReal grad_yz = 0.0;
179 370576 : ADReal grad_zx = 0.0;
180 370576 : ADReal grad_zy = 0.0;
181 370576 : ADReal grad_zz = 0.0;
182 :
183 370576 : auto trace = Sij_xx / 3.0;
184 :
185 370576 : if (_dim >= 2)
186 : {
187 370576 : const auto & grad_v = (*_v_var).gradient(elem_arg, state);
188 370576 : Sij_xy = grad_u(1) + grad_v(0);
189 741152 : Sij_yy = 2.0 * grad_v(1);
190 :
191 370576 : grad_xy = grad_u(1);
192 370576 : grad_yx = grad_v(0);
193 370576 : grad_yy = grad_v(1);
194 :
195 741152 : trace += Sij_yy / 3.0;
196 :
197 370576 : if (_dim >= 3)
198 : {
199 0 : const auto & grad_w = (*_w_var).gradient(elem_arg, state);
200 :
201 0 : Sij_xz = grad_u(2) + grad_w(0);
202 0 : Sij_yz = grad_v(2) + grad_w(1);
203 0 : Sij_zz = 2.0 * grad_w(2);
204 :
205 0 : grad_xz = grad_u(2);
206 0 : grad_yz = grad_v(2);
207 0 : grad_zx = grad_w(0);
208 0 : grad_zy = grad_w(1);
209 0 : grad_zz = grad_w(2);
210 :
211 0 : trace += Sij_zz / 3.0;
212 : }
213 : }
214 :
215 : const auto symmetric_strain_tensor_sq_norm =
216 370576 : (Sij_xx - trace) * grad_xx + Sij_xy * grad_xy + Sij_xz * grad_xz + Sij_xy * grad_yx +
217 370576 : (Sij_yy - trace) * grad_yy + Sij_yz * grad_yz + Sij_xz * grad_zx + Sij_yz * grad_zy +
218 370576 : (Sij_zz - trace) * grad_zz;
219 :
220 370576 : ADReal production_k = _mu_t(elem_arg, state) * symmetric_strain_tensor_sq_norm;
221 : // Compute production limiter (needed for flows with stagnation zones)
222 : const auto eps_old =
223 370576 : _newton_solve ? std::max(_var(elem_arg, old_state), 1e-10) : _var(elem_arg, old_state);
224 370576 : const ADReal production_limit = _C_pl * rho * eps_old;
225 : // Apply production limiter
226 370576 : production_k = std::min(production_k, production_limit);
227 :
228 370576 : const auto time_scale = raw_value(TKE_old) / raw_value(eps_old);
229 741152 : production = _C1_eps * production_k / time_scale;
230 741152 : destruction = _C2_eps * rho * _var(elem_arg, state) / time_scale;
231 :
232 370576 : residual = destruction - production;
233 : }
234 :
235 482584 : return residual;
236 : }
|
