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 "kEpsilonViscosityAux.h"
11 : #include "NavierStokesMethods.h"
12 : #include "NonlinearSystemBase.h"
13 : #include "libmesh/nonlinear_solver.h"
14 :
15 : registerMooseObject("NavierStokesApp", kEpsilonViscosityAux);
16 :
17 : InputParameters
18 907 : kEpsilonViscosityAux::validParams()
19 : {
20 907 : InputParameters params = AuxKernel::validParams();
21 907 : params.addClassDescription(
22 : "Calculates the turbulent viscosity according to the k-epsilon model.");
23 1814 : params.addRequiredParam<MooseFunctorName>("u", "The velocity in the x direction.");
24 1814 : params.addParam<MooseFunctorName>("v", "The velocity in the y direction.");
25 1814 : params.addParam<MooseFunctorName>("w", "The velocity in the z direction.");
26 1814 : params.addRequiredParam<MooseFunctorName>("k", "Coupled turbulent kinetic energy.");
27 1814 : params.deprecateParam("k", NS::TKE, "01/01/2025");
28 907 : params.addRequiredParam<MooseFunctorName>(NS::TKED,
29 : "Coupled turbulent kinetic energy dissipation rate.");
30 907 : params.addRequiredParam<MooseFunctorName>(NS::density, "Density");
31 907 : params.addRequiredParam<MooseFunctorName>(NS::mu, "Dynamic viscosity.");
32 1814 : params.addParam<Real>("C_mu", "Coupled turbulent kinetic energy closure.");
33 1814 : params.addParam<Real>("mu_t_ratio_max", 1e5, "Maximum allowable mu_t_ratio : mu/mu_t.");
34 1814 : params.addParam<std::vector<BoundaryName>>(
35 : "walls", {}, "Boundaries that correspond to solid walls.");
36 1814 : params.addParam<bool>("bulk_wall_treatment", false, "Activate bulk wall treatment.");
37 1814 : MooseEnum wall_treatment("eq_newton eq_incremental eq_linearized neq", "eq_newton");
38 1814 : 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 1814 : MooseEnum scale_limiter("none standard", "standard");
43 1814 : params.addParam<MooseEnum>("scale_limiter",
44 : scale_limiter,
45 : "The method used to limit the k-epsilon time scale."
46 : "'none', 'standard'");
47 1814 : params.addParam<bool>("newton_solve", false, "Whether a Newton nonlinear solve is being used");
48 1814 : params.addParamNamesToGroup("newton_solve", "Advanced");
49 :
50 907 : return params;
51 907 : }
52 :
53 485 : kEpsilonViscosityAux::kEpsilonViscosityAux(const InputParameters & params)
54 : : AuxKernel(params),
55 485 : _dim(_subproblem.mesh().dimension()),
56 970 : _u_var(getFunctor<Real>("u")),
57 1455 : _v_var(params.isParamValid("v") ? &(getFunctor<Real>("v")) : nullptr),
58 485 : _w_var(params.isParamValid("w") ? &(getFunctor<Real>("w")) : nullptr),
59 485 : _k(getFunctor<Real>(NS::TKE)),
60 485 : _epsilon(getFunctor<Real>(NS::TKED)),
61 485 : _rho(getFunctor<Real>(NS::density)),
62 485 : _mu(getFunctor<Real>(NS::mu)),
63 970 : _C_mu(getParam<Real>("C_mu")),
64 970 : _mu_t_ratio_max(getParam<Real>("mu_t_ratio_max")),
65 970 : _wall_boundary_names(getParam<std::vector<BoundaryName>>("walls")),
66 970 : _bulk_wall_treatment(getParam<bool>("bulk_wall_treatment")),
67 970 : _wall_treatment(getParam<MooseEnum>("wall_treatment").getEnum<NS::WallTreatmentEnum>()),
68 970 : _scale_limiter(getParam<MooseEnum>("scale_limiter")),
69 1455 : _newton_solve(getParam<bool>("newton_solve"))
70 : {
71 485 : }
72 :
73 : void
74 485 : kEpsilonViscosityAux::initialSetup()
75 : {
76 485 : if (!_wall_boundary_names.empty())
77 : {
78 882 : NS::getWallBoundedElements(
79 441 : _wall_boundary_names, _c_fe_problem, _subproblem, blockIDs(), _wall_bounded);
80 441 : NS::getWallDistance(_wall_boundary_names, _c_fe_problem, _subproblem, blockIDs(), _dist);
81 882 : NS::getElementFaceArgs(
82 441 : _wall_boundary_names, _c_fe_problem, _subproblem, blockIDs(), _face_infos);
83 : }
84 485 : }
85 :
86 : Real
87 4058172 : kEpsilonViscosityAux::computeValue()
88 : {
89 : // Convenient Arguments
90 4058172 : const Elem & elem = *_current_elem;
91 4058172 : const auto elem_arg = makeElemArg(_current_elem);
92 4058172 : const Moose::StateArg state = determineState();
93 4058172 : const auto rho = _rho(makeElemArg(_current_elem), state);
94 4058172 : const auto mu = _mu(makeElemArg(_current_elem), state);
95 4058172 : const auto nu = mu / rho;
96 :
97 : // Determine if the element is wall bounded
98 : // and if bulk wall treatment needs to be activated
99 4058172 : const bool wall_bounded = _wall_bounded.find(_current_elem) != _wall_bounded.end();
100 : Real mu_t;
101 :
102 : // Computing wall value for near-wall elements if bulk wall treatment is activated
103 : // bulk_wall_treatment should only be used for very coarse mesh problems
104 4058172 : if (wall_bounded && _bulk_wall_treatment)
105 : {
106 : // Computing wall value for turbulent dynamic viscosity
107 3960 : const auto & elem_distances = _dist[&elem];
108 : const auto min_wall_distance_iterator =
109 : (std::min_element(elem_distances.begin(), elem_distances.end()));
110 3960 : const auto min_wall_dist = *min_wall_distance_iterator;
111 : const size_t minIndex = std::distance(elem_distances.begin(), min_wall_distance_iterator);
112 3960 : const auto loc_normal = _face_infos[&elem][minIndex]->normal();
113 :
114 : // Getting y_plus
115 3960 : RealVectorValue velocity(_u_var(elem_arg, state));
116 3960 : if (_v_var)
117 3960 : velocity(1) = (*_v_var)(elem_arg, state);
118 3960 : if (_w_var)
119 0 : velocity(2) = (*_w_var)(elem_arg, state);
120 :
121 : // Compute the velocity and direction of the velocity component that is parallel to the wall
122 : const auto parallel_speed =
123 3960 : NS::computeSpeed<Real>(velocity - velocity * loc_normal * loc_normal);
124 :
125 : // Switch for determining the near wall quantities
126 : // wall_treatment can be: "eq_newton eq_incremental eq_linearized neq"
127 : Real y_plus = 0;
128 : Real mut_log; // turbulent log-layer viscosity
129 : Real mu_wall; // total wall viscosity to obtain the shear stress at the wall
130 :
131 3960 : if (_wall_treatment == NS::WallTreatmentEnum::EQ_NEWTON)
132 : {
133 : // Full Newton-Raphson solve to find the wall quantities from the law of the wall
134 3960 : const auto u_tau = NS::findUStar<Real>(mu, rho, parallel_speed, min_wall_dist);
135 3960 : y_plus = min_wall_dist * u_tau * rho / mu;
136 3960 : mu_wall = rho * Utility::pow<2>(u_tau) * min_wall_dist / parallel_speed;
137 3960 : mut_log = mu_wall - mu;
138 : }
139 0 : else if (_wall_treatment == NS::WallTreatmentEnum::EQ_INCREMENTAL)
140 : {
141 : // Incremental solve on y_plus to get the near-wall quantities
142 0 : y_plus = NS::findyPlus<Real>(mu, rho, std::max(parallel_speed, 1e-10), min_wall_dist);
143 0 : mu_wall = mu * (NS::von_karman_constant * y_plus /
144 0 : std::log(std::max(NS::E_turb_constant * y_plus, 1 + 1e-4)));
145 0 : mut_log = mu_wall - mu;
146 : }
147 0 : else if (_wall_treatment == NS::WallTreatmentEnum::EQ_LINEARIZED)
148 : {
149 : // Linearized approximation to the wall function to find the near-wall quantities faster
150 : const Real a_c = 1 / NS::von_karman_constant;
151 : const Real b_c = 1 / NS::von_karman_constant *
152 0 : (std::log(NS::E_turb_constant * std::max(min_wall_dist, 1.0) / mu) + 1.0);
153 : const Real c_c = parallel_speed;
154 :
155 0 : const auto u_tau = (-b_c + std::sqrt(std::pow(b_c, 2) + 4.0 * a_c * c_c)) / (2.0 * a_c);
156 0 : y_plus = min_wall_dist * u_tau * rho / mu;
157 0 : mu_wall = rho * Utility::pow<2>(u_tau) * min_wall_dist / parallel_speed;
158 0 : mut_log = mu_wall - mu;
159 : }
160 0 : else if (_wall_treatment == NS::WallTreatmentEnum::NEQ)
161 : {
162 : // Assign non-equilibrium wall function value
163 0 : y_plus = min_wall_dist * std::sqrt(std::sqrt(_C_mu) * _k(elem_arg, state)) * rho / mu;
164 0 : mu_wall = mu * (NS::von_karman_constant * y_plus /
165 0 : std::log(std::max(NS::E_turb_constant * y_plus, 1 + 1e-4)));
166 0 : mut_log = mu_wall - mu;
167 : }
168 : else
169 : mooseAssert(false, "For `kEpsilonViscosityAux` , wall treatment should not reach here");
170 :
171 3960 : if (y_plus <= 5.0)
172 : // sub-laminar layer
173 3960 : mu_t = 0.0;
174 0 : else if (y_plus >= 30.0)
175 : // log-layer
176 0 : mu_t = std::max(mut_log, NS::mu_t_low_limit);
177 : else
178 : {
179 : // buffer layer
180 0 : const auto blending_function = (y_plus - 5.0) / 25.0;
181 : // the blending depends on the mut_log at y+=30
182 0 : const auto mut_log = mu * (NS::von_karman_constant * 30.0 /
183 : std::log(std::max(NS::E_turb_constant * 30.0, 1 + 1e-4)) -
184 0 : 1.0);
185 0 : mu_t = std::max(blending_function * mut_log, NS::mu_t_low_limit);
186 : }
187 3960 : }
188 : else
189 : {
190 : Real time_scale;
191 4054212 : if (_scale_limiter == "standard")
192 : {
193 4041252 : time_scale = std::max(_k(elem_arg, state) / _epsilon(elem_arg, state),
194 8082504 : std::sqrt(nu / _epsilon(elem_arg, state)));
195 : }
196 : else
197 : {
198 12960 : time_scale = _k(elem_arg, state) / _epsilon(elem_arg, state);
199 : }
200 : // For newton solvers, epsilon might not be bounded
201 4054212 : if (_newton_solve)
202 46309 : time_scale = _k(elem_arg, state) / std::max(NS::epsilon_low_limit, _epsilon(elem_arg, state));
203 :
204 4054212 : const Real mu_t_nl = _rho(elem_arg, state) * _C_mu * _k(elem_arg, state) * time_scale;
205 4054212 : mu_t = mu_t_nl;
206 4054212 : if (_newton_solve)
207 25940 : mu_t = std::max(mu_t, NS::mu_t_low_limit);
208 : }
209 : // Turbulent viscosity limiter
210 4058172 : return std::min(mu_t, _mu_t_ratio_max * mu);
211 : }
|
