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