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