https://mooseframework.inl.gov
kEpsilonViscosityAux.C
Go to the documentation of this file.
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 
19 {
21  params.addClassDescription(
22  "Calculates the turbulent viscosity according to the k-epsilon model.");
23  params.addRequiredParam<MooseFunctorName>("u", "The velocity in the x direction.");
24  params.addParam<MooseFunctorName>("v", "The velocity in the y direction.");
25  params.addParam<MooseFunctorName>("w", "The velocity in the z direction.");
26  params.addRequiredParam<MooseFunctorName>("k", "Coupled turbulent kinetic energy.");
27  params.deprecateParam("k", NS::TKE, "01/01/2025");
28  params.addRequiredParam<MooseFunctorName>(NS::TKED,
29  "Coupled turbulent kinetic energy dissipation rate.");
30  params.addRequiredParam<MooseFunctorName>(NS::density, "Density");
31  params.addRequiredParam<MooseFunctorName>(NS::mu, "Dynamic viscosity.");
32  params.addParam<Real>("C_mu", "Coupled turbulent kinetic energy closure.");
33  params.addParam<Real>("mu_t_ratio_max", 1e5, "Maximum allowable mu_t_ratio : mu/mu_t.");
34  params.addParam<std::vector<BoundaryName>>(
35  "walls", {}, "Boundaries that correspond to solid walls.");
36  params.addParam<bool>("bulk_wall_treatment", false, "Activate bulk wall treatment.");
37  MooseEnum wall_treatment("eq_newton eq_incremental eq_linearized neq", "eq_newton");
38  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  MooseEnum scale_limiter("none standard", "standard");
43  params.addParam<MooseEnum>("scale_limiter",
44  scale_limiter,
45  "The method used to limit the k-epsilon time scale."
46  "'none', 'standard'");
47  params.addParam<bool>("newton_solve", false, "Whether a Newton nonlinear solve is being used");
48  params.addParamNamesToGroup("newton_solve", "Advanced");
49 
50  return params;
51 }
52 
54  : AuxKernel(params),
55  _dim(_subproblem.mesh().dimension()),
56  _u_var(getFunctor<Real>("u")),
57  _v_var(params.isParamValid("v") ? &(getFunctor<Real>("v")) : nullptr),
58  _w_var(params.isParamValid("w") ? &(getFunctor<Real>("w")) : nullptr),
59  _k(getFunctor<Real>(NS::TKE)),
60  _epsilon(getFunctor<Real>(NS::TKED)),
61  _rho(getFunctor<Real>(NS::density)),
62  _mu(getFunctor<Real>(NS::mu)),
63  _C_mu(getParam<Real>("C_mu")),
64  _mu_t_ratio_max(getParam<Real>("mu_t_ratio_max")),
65  _wall_boundary_names(getParam<std::vector<BoundaryName>>("walls")),
66  _bulk_wall_treatment(getParam<bool>("bulk_wall_treatment")),
67  _wall_treatment(getParam<MooseEnum>("wall_treatment").getEnum<NS::WallTreatmentEnum>()),
68  _scale_limiter(getParam<MooseEnum>("scale_limiter")),
69  _newton_solve(getParam<bool>("newton_solve"))
70 {
71 }
72 
73 void
75 {
76  if (!_wall_boundary_names.empty())
77  {
83  }
84 }
85 
86 Real
88 {
89  // Convenient Arguments
90  const Elem & elem = *_current_elem;
91  const auto elem_arg = makeElemArg(_current_elem);
92  const Moose::StateArg state = determineState();
93  const auto rho = _rho(makeElemArg(_current_elem), state);
94  const auto mu = _mu(makeElemArg(_current_elem), state);
95  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  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  if (wall_bounded && _bulk_wall_treatment)
105  {
106  // Computing wall value for turbulent dynamic viscosity
107  const auto & elem_distances = _dist[&elem];
108  const auto min_wall_distance_iterator =
109  (std::min_element(elem_distances.begin(), elem_distances.end()));
110  const auto min_wall_dist = *min_wall_distance_iterator;
111  const size_t minIndex = std::distance(elem_distances.begin(), min_wall_distance_iterator);
112  const auto loc_normal = _face_infos[&elem][minIndex]->normal();
113 
114  // Getting y_plus
115  RealVectorValue velocity(_u_var(elem_arg, state));
116  if (_v_var)
117  velocity(1) = (*_v_var)(elem_arg, state);
118  if (_w_var)
119  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  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 
132  {
133  // Full Newton-Raphson solve to find the wall quantities from the law of the wall
134  const auto u_tau = NS::findUStar<Real>(mu, rho, parallel_speed, min_wall_dist);
135  y_plus = min_wall_dist * u_tau * rho / mu;
136  mu_wall = rho * Utility::pow<2>(u_tau) * min_wall_dist / parallel_speed;
137  mut_log = mu_wall - mu;
138  }
140  {
141  // Incremental solve on y_plus to get the near-wall quantities
142  y_plus = NS::findyPlus<Real>(mu, rho, std::max(parallel_speed, 1e-10), min_wall_dist);
143  mu_wall = mu * (NS::von_karman_constant * y_plus /
144  std::log(std::max(NS::E_turb_constant * y_plus, 1 + 1e-4)));
145  mut_log = mu_wall - mu;
146  }
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  (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  const auto u_tau = (-b_c + std::sqrt(std::pow(b_c, 2) + 4.0 * a_c * c_c)) / (2.0 * a_c);
156  y_plus = min_wall_dist * u_tau * rho / mu;
157  mu_wall = rho * Utility::pow<2>(u_tau) * min_wall_dist / parallel_speed;
158  mut_log = mu_wall - mu;
159  }
161  {
162  // Assign non-equilibrium wall function value
163  y_plus = min_wall_dist * std::sqrt(std::sqrt(_C_mu) * _k(elem_arg, state)) * rho / mu;
164  mu_wall = mu * (NS::von_karman_constant * y_plus /
165  std::log(std::max(NS::E_turb_constant * y_plus, 1 + 1e-4)));
166  mut_log = mu_wall - mu;
167  }
168  else
169  mooseAssert(false, "For `kEpsilonViscosityAux` , wall treatment should not reach here");
170 
171  if (y_plus <= 5.0)
172  // sub-laminar layer
173  mu_t = 0.0;
174  else if (y_plus >= 30.0)
175  // log-layer
176  mu_t = std::max(mut_log, NS::mu_t_low_limit);
177  else
178  {
179  // buffer layer
180  const auto blending_function = (y_plus - 5.0) / 25.0;
181  // the blending depends on the mut_log at y+=30
182  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  1.0);
185  mu_t = std::max(blending_function * mut_log, NS::mu_t_low_limit);
186  }
187  }
188  else
189  {
190  Real time_scale;
191  if (_scale_limiter == "standard")
192  {
193  time_scale = std::max(_k(elem_arg, state) / _epsilon(elem_arg, state),
194  std::sqrt(nu / _epsilon(elem_arg, state)));
195  }
196  else
197  {
198  time_scale = _k(elem_arg, state) / _epsilon(elem_arg, state);
199  }
200  // For newton solvers, epsilon might not be bounded
201  if (_newton_solve)
202  time_scale = _k(elem_arg, state) / std::max(NS::epsilon_low_limit, _epsilon(elem_arg, state));
203 
204  const Real mu_t_nl = _rho(elem_arg, state) * _C_mu * _k(elem_arg, state) * time_scale;
205  mu_t = mu_t_nl;
206  if (_newton_solve)
207  mu_t = std::max(mu_t, NS::mu_t_low_limit);
208  }
209  // Turbulent viscosity limiter
210  return std::min(mu_t, _mu_t_ratio_max * mu);
211 }
std::unordered_set< const Elem * > _wall_bounded
virtual void initialSetup() override
static constexpr Real von_karman_constant
Definition: NS.h:191
static const std::string mu_t
Definition: NS.h:125
const Moose::Functor< Real > & _mu
Dynamic viscosity.
void addParam(const std::string &name, const std::initializer_list< typename T::value_type > &value, const std::string &doc_string)
const Moose::Functor< Real > & _rho
Density.
Computes the turbuent viscosity for the k-Epsilon model.
const bool _bulk_wall_treatment
If the user wants to enable bulk wall treatment.
Moose::StateArg determineState() const
MeshBase & mesh
static const std::string density
Definition: NS.h:33
static const std::string TKE
Definition: NS.h:176
WallTreatmentEnum
Wall treatment options.
Definition: NS.h:182
template Real findUStar< Real >(const Real &mu, const Real &rho, const Real &u, const Real dist)
void getWallBoundedElements(const std::vector< BoundaryName > &wall_boundary_name, const FEProblemBase &fe_problem, const SubProblem &subproblem, const std::set< SubdomainID > &block_ids, std::unordered_set< const Elem *> &wall_bounded)
Map marking wall bounded elements The map passed in wall_bounded_map gets cleared and re-populated...
virtual const std::set< SubdomainID > & blockIDs() const
const Moose::Functor< Real > & _epsilon
Turbulent kinetic energy dissipation rate.
const Moose::Functor< Real > & _u_var
x-velocity
const Moose::Functor< Real > & _k
Turbulent kinetic energy.
void addRequiredParam(const std::string &name, const std::string &doc_string)
FEProblemBase & _c_fe_problem
Moose::ElemArg makeElemArg(const Elem *elem, bool correct_skewnewss=false) const
template Real findyPlus< Real >(const Real &mu, const Real &rho, const Real &u, Real dist)
void deprecateParam(const std::string &old_name, const std::string &new_name, const std::string &removal_date)
NS::WallTreatmentEnum _wall_treatment
Method used for wall treatment.
static const std::string mu
Definition: NS.h:123
const bool _newton_solve
Whether we are using a newton solve.
std::map< const Elem *, std::vector< const FaceInfo * > > _face_infos
void getWallDistance(const std::vector< BoundaryName > &wall_boundary_name, const FEProblemBase &fe_problem, const SubProblem &subproblem, const std::set< SubdomainID > &block_ids, std::map< const Elem *, std::vector< Real >> &dist_map)
Map storing wall ditance for near-wall marked elements The map passed in dist_map gets cleared and re...
std::map< const Elem *, std::vector< Real > > _dist
virtual Real computeValue() override
kEpsilonViscosityAux(const InputParameters &parameters)
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
const MooseEnum _scale_limiter
Method used to limit the k-e time scale.
static constexpr Real E_turb_constant
Definition: NS.h:192
const std::vector< BoundaryName > & _wall_boundary_names
Wall boundaries.
static const std::string TKED
Definition: NS.h:177
SubProblem & _subproblem
const Elem *const & _current_elem
void addClassDescription(const std::string &doc_string)
void getElementFaceArgs(const std::vector< BoundaryName > &wall_boundary_name, const FEProblemBase &fe_problem, const SubProblem &subproblem, const std::set< SubdomainID > &block_ids, std::map< const Elem *, std::vector< const FaceInfo *>> &face_info_map)
Map storing face arguments to wall bounded faces The map passed in face_info_map gets cleared and re-...
static InputParameters validParams()
static const std::string velocity
Definition: NS.h:45
registerMooseObject("NavierStokesApp", kEpsilonViscosityAux)
template Real computeSpeed< Real >(const libMesh::VectorValue< Real > &velocity)
static constexpr Real mu_t_low_limit
Definition: NS.h:194
static InputParameters validParams()
const Moose::Functor< Real > * _w_var
z-velocity
MooseUnits pow(const MooseUnits &, int)
static constexpr Real epsilon_low_limit
Definition: NS.h:196
const Real _C_mu
C-mu closure coefficient.
const Moose::Functor< Real > * _v_var
y-velocity
void addParamNamesToGroup(const std::string &space_delim_names, const std::string group_name)