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 925 : kEpsilonViscosityAux::validParams()
19 : {
20 925 : InputParameters params = AuxKernel::validParams();
21 925 : params.addClassDescription(
22 : "Calculates the turbulent viscosity according to the k-epsilon model.");
23 1850 : params.addRequiredParam<MooseFunctorName>("u", "The velocity in the x direction.");
24 1850 : params.addParam<MooseFunctorName>("v", "The velocity in the y direction.");
25 1850 : params.addParam<MooseFunctorName>("w", "The velocity in the z direction.");
26 925 : params.addRequiredParam<MooseFunctorName>(NS::TKE, "Coupled turbulent kinetic energy.");
27 925 : params.addRequiredParam<MooseFunctorName>(NS::TKED,
28 : "Coupled turbulent kinetic energy dissipation rate.");
29 925 : params.addRequiredParam<MooseFunctorName>(NS::density, "Density");
30 925 : params.addRequiredParam<MooseFunctorName>(NS::mu, "Dynamic viscosity.");
31 1850 : params.addParam<Real>("C_mu", "Coupled turbulent kinetic energy closure.");
32 1850 : params.addParam<Real>("mu_t_ratio_max", 1e5, "Maximum allowable mu_t_ratio : mu/mu_t.");
33 1850 : params.addParam<std::vector<BoundaryName>>(
34 : "walls", {}, "Boundaries that correspond to solid walls.");
35 1850 : params.addParam<bool>("bulk_wall_treatment", false, "Activate bulk wall treatment.");
36 1850 : MooseEnum wall_treatment("eq_newton eq_incremental eq_linearized neq", "eq_newton");
37 1850 : 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 1850 : MooseEnum scale_limiter("none standard", "standard");
42 1850 : params.addParam<MooseEnum>("scale_limiter",
43 : scale_limiter,
44 : "The method used to limit the k-epsilon time scale."
45 : "'none', 'standard'");
46 1850 : params.addParam<bool>("newton_solve", false, "Whether a Newton nonlinear solve is being used");
47 1850 : params.addParamNamesToGroup("newton_solve", "Advanced");
48 :
49 925 : return params;
50 925 : }
51 :
52 494 : kEpsilonViscosityAux::kEpsilonViscosityAux(const InputParameters & params)
53 : : AuxKernel(params),
54 494 : _dim(_subproblem.mesh().dimension()),
55 988 : _u_var(getFunctor<Real>("u")),
56 1482 : _v_var(params.isParamValid("v") ? &(getFunctor<Real>("v")) : nullptr),
57 494 : _w_var(params.isParamValid("w") ? &(getFunctor<Real>("w")) : nullptr),
58 494 : _k(getFunctor<Real>(NS::TKE)),
59 494 : _epsilon(getFunctor<Real>(NS::TKED)),
60 494 : _rho(getFunctor<Real>(NS::density)),
61 494 : _mu(getFunctor<Real>(NS::mu)),
62 988 : _C_mu(getParam<Real>("C_mu")),
63 988 : _mu_t_ratio_max(getParam<Real>("mu_t_ratio_max")),
64 988 : _wall_boundary_names(getParam<std::vector<BoundaryName>>("walls")),
65 988 : _bulk_wall_treatment(getParam<bool>("bulk_wall_treatment")),
66 988 : _wall_treatment(getParam<MooseEnum>("wall_treatment").getEnum<NS::WallTreatmentEnum>()),
67 988 : _scale_limiter(getParam<MooseEnum>("scale_limiter")),
68 1482 : _newton_solve(getParam<bool>("newton_solve"))
69 : {
70 494 : }
71 :
72 : void
73 494 : kEpsilonViscosityAux::initialSetup()
74 : {
75 494 : if (!_wall_boundary_names.empty())
76 : {
77 900 : NS::getWallBoundedElements(
78 450 : _wall_boundary_names, _c_fe_problem, _subproblem, blockIDs(), _wall_bounded);
79 450 : NS::getWallDistance(_wall_boundary_names, _c_fe_problem, _subproblem, blockIDs(), _dist);
80 900 : NS::getElementFaceArgs(
81 450 : _wall_boundary_names, _c_fe_problem, _subproblem, blockIDs(), _face_infos);
82 : }
83 494 : }
84 :
85 : Real
86 4442748 : kEpsilonViscosityAux::computeValue()
87 : {
88 : // Convenient Arguments
89 4442748 : const Elem & elem = *_current_elem;
90 4442748 : const auto elem_arg = makeElemArg(_current_elem);
91 4442748 : const Moose::StateArg state = determineState();
92 4442748 : const auto rho = _rho(makeElemArg(_current_elem), state);
93 4442748 : const auto mu = _mu(makeElemArg(_current_elem), state);
94 4442748 : 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 4442748 : 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 4442748 : if (wall_bounded && _bulk_wall_treatment)
104 : {
105 : // Computing wall value for turbulent dynamic viscosity
106 3960 : const auto & elem_distances = _dist[&elem];
107 : const auto min_wall_distance_iterator =
108 : (std::min_element(elem_distances.begin(), elem_distances.end()));
109 3960 : const auto min_wall_dist = *min_wall_distance_iterator;
110 : const size_t minIndex = std::distance(elem_distances.begin(), min_wall_distance_iterator);
111 3960 : const auto loc_normal = _face_infos[&elem][minIndex]->normal();
112 :
113 : // Getting y_plus
114 3960 : RealVectorValue velocity(_u_var(elem_arg, state));
115 3960 : if (_v_var)
116 3960 : velocity(1) = (*_v_var)(elem_arg, state);
117 3960 : 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 3960 : 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 3960 : 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 3960 : const auto u_tau = NS::findUStar<Real>(mu, rho, parallel_speed, min_wall_dist);
134 3960 : y_plus = min_wall_dist * u_tau * rho / mu;
135 3960 : mu_wall = rho * Utility::pow<2>(u_tau) * min_wall_dist / parallel_speed;
136 3960 : 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 3960 : if (y_plus <= 5.0)
171 : // sub-laminar layer
172 3960 : 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 3960 : }
187 : else
188 : {
189 : Real time_scale;
190 4438788 : if (_scale_limiter == "standard")
191 : {
192 4425828 : time_scale = std::max(_k(elem_arg, state) / _epsilon(elem_arg, state),
193 8851656 : std::sqrt(nu / _epsilon(elem_arg, state)));
194 : }
195 : else
196 : {
197 12960 : time_scale = _k(elem_arg, state) / _epsilon(elem_arg, state);
198 : }
199 : // For newton solvers, epsilon might not be bounded
200 4438788 : if (_newton_solve)
201 46309 : time_scale = _k(elem_arg, state) / std::max(NS::epsilon_low_limit, _epsilon(elem_arg, state));
202 :
203 4438788 : const Real mu_t_nl = _rho(elem_arg, state) * _C_mu * _k(elem_arg, state) * time_scale;
204 4438788 : mu_t = mu_t_nl;
205 4438788 : if (_newton_solve)
206 25940 : mu_t = std::max(mu_t, NS::mu_t_low_limit);
207 : }
208 : // Turbulent viscosity limiter
209 4442748 : return std::min(mu_t, _mu_t_ratio_max * mu);
210 : }
|
