diff --git a/modules/navier_stokes/doc/content/source/auxkernels/RANSYPlusAux.md b/modules/navier_stokes/doc/content/source/auxkernels/RANSYPlusAux.md index ea3d063dd113..89e47a631ee8 100644 --- a/modules/navier_stokes/doc/content/source/auxkernels/RANSYPlusAux.md +++ b/modules/navier_stokes/doc/content/source/auxkernels/RANSYPlusAux.md @@ -2,7 +2,7 @@ Computes the dimensionless wall distance $y^+$ for RANS models. -Four different formulations are supported as defined by the "wall_treatment" parameter. +Four different formulations are supported as defined by the [!param](/AuxKernels/RANSYPlusAux/wall_treatment) parameter. Details on each of the four formulations can be found in [INSFVTurbulentViscosityWallFunction](source/fvbcs/INSFVTurbulentViscosityWallFunction.md). diff --git a/modules/navier_stokes/doc/content/source/fvbcs/INSFVTurbulentViscosityWallFunction.md b/modules/navier_stokes/doc/content/source/fvbcs/INSFVTurbulentViscosityWallFunction.md index ec6f135dfa1c..0df2d34f9ece 100644 --- a/modules/navier_stokes/doc/content/source/fvbcs/INSFVTurbulentViscosityWallFunction.md +++ b/modules/navier_stokes/doc/content/source/fvbcs/INSFVTurbulentViscosityWallFunction.md @@ -17,14 +17,23 @@ $\mu_w = \mu + \mu_t $ , such that the wall shear stress $\tau_w$ is accurately without the need of fully resolving the gradients at the near wall region. \begin{equation} - \tau_w = /frac{ \mu_w $u_p$}{y_p} + \tau_w = /frac{ \mu_w u_p}{y_p} \end{equation} -To obtain the relationship between the wall shear stress and the dimensionless wall distance, -four different formulations are supportedas defined by the -[!param](/FVBCs/INSFVTurbulentViscosityWallFunction/wall_treatment) parameter. +where: -To define a proper grid spacing during the meshing process, we recommend using the Auxiliary Kernel +- $\mu_w = \mu + \mu_t$ is the total viscosity evaluated at the wall face +- $\mu_t$ is the turbulent viscosity, evaluated at the wall for the purpose of this boundary condition +- $\mu$ is the kinematic viscosity, evaluated at the wall for the purpose of this boundary condition +- $\tau_w$ is the wall-shear stress +- $u_p$ is the wall-parallel velocity at the centroid +- $y_p$ is the wall normal distance to the centroid + +To impose a correct boundary condition for $\mu_t$, as seen in the Equation above, we need to compute $\tau_w$ using analytical +relationships between the wall shear stress and the dimensionless wall distance $y^+$. For this purpose, four different +formulations are supported as defined by the [!param](/FVBCs/INSFVTurbulentViscosityWallFunction/wall_treatment) parameter. + +To set the grid spacing for the first cell near the wall in your mesh, we recommend using the Auxiliary Kernel [RANSYPlusAux.md] to estimate the dimensionless wall distance $y^+$. @@ -46,6 +55,7 @@ for the turbulent viscosity. where: - $\rho$ is the density +- $\mu$ is the kinematic viscosity - $u_{\tau} = \sqrt{\frac{\tau_w}{\rho}}$ is the friction velocity and $\tau_w$ is the wall friction - $y_p$ is the distance from the boundary to the center of the near-wall cell - $u_p$ is the parallel velocity to the boundary computed at the center of the near-wall cell diff --git a/modules/navier_stokes/doc/content/source/fvkernels/INSFVTKEDSourceSink.md b/modules/navier_stokes/doc/content/source/fvkernels/INSFVTKEDSourceSink.md index b05d43244729..ee2b89dd0534 100644 --- a/modules/navier_stokes/doc/content/source/fvkernels/INSFVTKEDSourceSink.md +++ b/modules/navier_stokes/doc/content/source/fvkernels/INSFVTKEDSourceSink.md @@ -8,7 +8,7 @@ A different treatment is used for the bulk and the near wall regions. ## Bulk formulation: -The turbulent production $G_\epsilon$ is modeled as follows: +The production of turbulent kinetic energy dissipation $G_\epsilon$ is modeled as follows: \begin{equation} G_{\epsilon} = C_{1,\epsilon} \frac{\epsilon}{k} G_k \,, @@ -29,6 +29,7 @@ where: - $C_{2,\epsilon} = 1.92$ is a closure parameter, - $\epsilon$ is the solution variable, i.e., the dissipation rate of the turbulent kinetic energy, +- $k$ is the turbulent kinetic energy, - $t_k = \frac{k}{\epsilon}$ is the turbulent time scale; if the [!param](/FVKernels/INSFVTKEDSourceSink/linearized_model) is `true`, this timescale is computed from the previous iteration; if [!param](/FVKernels/INSFVTKEDSourceSink/linearized_model) is `false`, in a nonlinear solve, this timescale is aded to the Jacobian. ## Wall formulation: diff --git a/modules/navier_stokes/doc/content/source/fvkernels/INSFVTKESourceSink.md b/modules/navier_stokes/doc/content/source/fvkernels/INSFVTKESourceSink.md index 44853034fdba..25e99d25deec 100644 --- a/modules/navier_stokes/doc/content/source/fvkernels/INSFVTKESourceSink.md +++ b/modules/navier_stokes/doc/content/source/fvkernels/INSFVTKESourceSink.md @@ -32,7 +32,7 @@ G_k = min \left( G_k , C_{PL} \rho \epsilon \right) \,, where: -- $C_{PL}$ it the limiter constant, and set to a recommended value of 10 . +- $C_{PL}$ it the limiter constant, and set by default to a recommended value of 10 \cite{durbin1996k}. ## Wall formulation: @@ -59,6 +59,8 @@ incremental fixed-point search algorithm. The cells belonging to the `sub-laminar` boundary layers are defined as those for which $y^+ < 11.25$. The ones belonging to the `logarithmic` boundary layer are those for which $y^+ \ge 11.25$. +The imposed threshold of $y^+ = 11.25$ is given by the value of $y^+$ at which the `sub-laminar` +and `logarithmic` boundary profiles intersect. In the `sub-laminar` region production of turbulent kinetic energy is negligible, therefore, if $y^+ \lt 11.25$: @@ -69,7 +71,7 @@ G_k = 0.0 \,, In the `logarithmic` boundary layers the production term is no longer negligible and is defined as: \begin{equation} -G_k = \tau_w ||\nabla \vec{u}|| = \left( \mu_t + \mu \right) ||\nabla \vec{u}|| \frac{ C_{\mu}^{0.25} \sqrt(k)}{\kappa y_p} \,, +G_k = \tau_w ||\nabla \vec{u}|| = \mu_w ||\nabla \vec{u}|| \frac{ C_{\mu}^{0.25} \sqrt(k)}{\kappa y_p} \,, \end{equation} where: @@ -84,7 +86,7 @@ The formulation assumes that the near wall value is already imposed in the $\mu_ When solving a linear problem, instead of the nonlinear formulation, the production term is formulated as: \begin{equation} -G_k = \left( \mu_t + \mu \right) ||\nabla \vec{u}|| \frac{ C_{\mu}^{0.25} k}{\sqrt{k_{old}} \kappa y_p} \,. +G_k = \mu_w ||\nabla \vec{u}|| \frac{ C_{\mu}^{0.25} k}{\sqrt{k_{old}} \kappa y_p} \,. \end{equation} where: diff --git a/modules/navier_stokes/include/auxkernels/kEpsilonViscosityAux.h b/modules/navier_stokes/include/auxkernels/kEpsilonViscosityAux.h index ad96f5bde298..1062bba8bedd 100644 --- a/modules/navier_stokes/include/auxkernels/kEpsilonViscosityAux.h +++ b/modules/navier_stokes/include/auxkernels/kEpsilonViscosityAux.h @@ -79,90 +79,3 @@ class kEpsilonViscosityAux : public AuxKernel std::map> _face_infos; ///@} }; - -// WORKING VERSION // - -//* This file is part of the MOOSE framework -//* https://www.mooseframework.org -//* -//* All rights reserved, see COPYRIGHT for full restrictions -//* https://github.com/idaholab/moose/blob/master/COPYRIGHT -//* -//* Licensed under LGPL 2.1, please see LICENSE for details -//* https://www.gnu.org/licenses/lgpl-2.1.html - -// #pragma once - -// #include "AuxKernel.h" -// #include "INSFVVelocityVariable.h" - -// /** -// * Computes the turbuent viscosity for the k-Epsilon model. -// * Implements two near-wall treatments: equilibrium and non-equilibrium wall functions. -// */ -// class kEpsilonViscosityAux : public AuxKernel -// { -// public: -// static InputParameters validParams(); - -// virtual void initialSetup() override; - -// kEpsilonViscosityAux(const InputParameters & parameters); - -// protected: -// virtual Real computeValue() override; - -// /// Local method to find friction velocity. -// /// This method may need to be reimplemented for each new turbulence model -// ADReal findUStarLocalMethod(const ADReal & u, const Real & dist); - -// /// The dimension of the domain -// const unsigned int _dim; - -// /// x-velocity -// const Moose::Functor & _u_var; -// /// y-velocity -// const Moose::Functor * _v_var; -// /// z-velocity -// const Moose::Functor * _w_var; - -// /// Turbulent kinetic energy -// const Moose::Functor & _k; -// /// Turbulent kinetic energy dissipation rate -// const Moose::Functor & _epsilon; - -// /// Density -// const Moose::Functor & _rho; -// /// Dynamic viscosity -// const Moose::Functor & _mu; - -// /// C-mu closure coefficient -// const Real _C_mu; - -// /// Wall boundaries -// const std::vector & _wall_boundary_names; - -// /// If the user wants the linearized computation of y_plus -// const bool _linearized_yplus; - -// /// If the user wants to enable bulk wall treatment -// const bool _bulk_wall_treatment; - -// /// IF the user requested non-equilibrium wall treatment -// const bool _non_equilibrium_treatment; - -// // -- Parameters of the wall function method - -// /// Maximum number of iterations to find the friction velocity -// static constexpr int _MAX_ITERS_U_TAU{50}; - -// /// Relative tolerance to find the friction velocity -// static constexpr Real _REL_TOLERANCE{1e-4}; - -// ///@{ -// /// Maps for wall bounded elements -// std::map _wall_bounded; -// std::map> _dist; -// std::map> _face_infos; -// ///@} -// }; diff --git a/modules/navier_stokes/include/base/NS.h b/modules/navier_stokes/include/base/NS.h index 25d3e4a05a5e..06ce16f900be 100644 --- a/modules/navier_stokes/include/base/NS.h +++ b/modules/navier_stokes/include/base/NS.h @@ -164,6 +164,8 @@ static const std::string TKED = "epsilon"; // Turbulence constants static constexpr Real von_karman_constant = 0.4187; static constexpr Real E_turb_constant = 9.793; +// Lower limit for mu_t +static constexpr Real mu_t_low_limit = 1.0e-12; } namespace NS_DEFAULT_VALUES diff --git a/modules/navier_stokes/include/fvbcs/INSFVTurbulentViscosityWallFunction.h b/modules/navier_stokes/include/fvbcs/INSFVTurbulentViscosityWallFunction.h index 257988b8b99b..bd92dc76efd5 100644 --- a/modules/navier_stokes/include/fvbcs/INSFVTurbulentViscosityWallFunction.h +++ b/modules/navier_stokes/include/fvbcs/INSFVTurbulentViscosityWallFunction.h @@ -51,57 +51,3 @@ class INSFVTurbulentViscosityWallFunction : public FVDirichletBCBase /// Method used for wall treatment const MooseEnum _wall_treatment; }; - -// //* This file is part of the MOOSE framework -// //* https://www.mooseframework.org -// //* -// //* All rights reserved, see COPYRIGHT for full restrictions -// //* https://github.com/idaholab/moose/blob/master/COPYRIGHT -// //* -// //* Licensed under LGPL 2.1, please see LICENSE for details -// //* https://www.gnu.org/licenses/lgpl-2.1.html - -// #pragma once - -// #include "FVDirichletBCBase.h" -// #include "FVFluxBC.h" - -// /** -// * Applies a wall function to the turbulent viscosity field -// */ -// class INSFVTurbulentViscosityWallFunction : public FVDirichletBCBase -// { -// public: -// INSFVTurbulentViscosityWallFunction(const InputParameters & parameters); - -// static InputParameters validParams(); - -// ADReal boundaryValue(const FaceInfo & fi) const override; - -// private: -// /// the dimension of the domain -// const unsigned int _dim; - -// /// x-velocity -// const Moose::Functor & _u_var; -// /// y-velocity -// const Moose::Functor * _v_var; -// /// z-velocity -// const Moose::Functor * _w_var; - -// /// Density -// const Moose::Functor & _rho; -// /// Dynamic viscosity -// const Moose::Functor & _mu; -// /// Turbulent dynamic viscosity -// const Moose::Functor & _mu_t; - -// /// Turbulent kinetic energy -// const Moose::Functor & _k; - -// /// C_mu turbulent coefficient -// const Real _C_mu; - -// /// Method used for wall treatment -// const MooseEnum _wall_treatment; -// }; diff --git a/modules/navier_stokes/include/fvkernels/INSFVTKEDSourceSink.h b/modules/navier_stokes/include/fvkernels/INSFVTKEDSourceSink.h index 16dbeefd1eb5..4358db073bb1 100644 --- a/modules/navier_stokes/include/fvkernels/INSFVTKEDSourceSink.h +++ b/modules/navier_stokes/include/fvkernels/INSFVTKEDSourceSink.h @@ -88,97 +88,3 @@ class INSFVTKEDSourceSink : public FVElementalKernel std::map> _face_infos; ///@} }; - -// //* This file is part of the MOOSE framework -// //* https://www.mooseframework.org -// //* -// //* All rights reserved, see COPYRIGHT for full restrictions -// //* https://github.com/idaholab/moose/blob/master/COPYRIGHT -// //* -// //* Licensed under LGPL 2.1, please see LICENSE for details -// //* https://www.gnu.org/licenses/lgpl-2.1.html - -// #pragma once - -// #include "FVElementalKernel.h" -// #include "MathFVUtils.h" -// #include "INSFVMomentumResidualObject.h" -// #include "INSFVVelocityVariable.h" - -// /** -// * Computes the source and sink terms for the turbulent kinetic energy dissipation rate. -// */ -// class INSFVTKEDSourceSink : public FVElementalKernel -// { -// public: -// static InputParameters validParams(); - -// virtual void initialSetup() override; - -// INSFVTKEDSourceSink(const InputParameters & parameters); - -// protected: -// ADReal computeQpResidual() override; - -// protected: -// /// The dimension of the simulation -// const unsigned int _dim; - -// /// x-velocity -// const Moose::Functor & _u_var; -// /// y-velocity -// const Moose::Functor * _v_var; -// /// z-velocity -// const Moose::Functor * _w_var; - -// /// Turbulent kinetic energy -// const Moose::Functor & _k; - -// /// Density -// const Moose::Functor & _rho; - -// /// Dynamic viscosity -// const Moose::Functor & _mu; - -// /// Turbulent dynamic viscosity -// const Moose::Functor & _mu_t; - -// /// Wall boundaries -// const std::vector & _wall_boundary_names; - -// /// Maximum mixing length allowed for the domain -// const Real _max_mixing_length; - -// /// If the user wants to use the linearized model -// const bool _linearized_model; - -// /// No equilibrium treatement -// const bool _non_equilibrium_treatment; - -// /// Value of the first epsilon closure coefficient -// const Real _C1_eps; - -// /// Value of the second epsilon closure coefficient -// const Real _C2_eps; - -// /// C_mu constant -// const Real _C_mu; - -// // Production Limiter Constant -// const Real _C_pl; - -// /// Stored strain rate -// std::map _symmetric_strain_tensor_norm_old; -// /// Map for the previous destruction field -// std::map _old_destruction; - -// /// Map for the previous nonlienar iterate -// std::map _pevious_nl_sol; - -// ///@{ -// /** Maps for wall treatment */ -// std::map _wall_bounded; -// std::map> _dist; -// std::map> _face_infos; -// ///@} -// }; diff --git a/modules/navier_stokes/include/fvkernels/INSFVTKESourceSink.h b/modules/navier_stokes/include/fvkernels/INSFVTKESourceSink.h index f61b494a3a54..c7dabd1d8714 100644 --- a/modules/navier_stokes/include/fvkernels/INSFVTKESourceSink.h +++ b/modules/navier_stokes/include/fvkernels/INSFVTKESourceSink.h @@ -72,81 +72,3 @@ class INSFVTKESourceSink : public FVElementalKernel std::map> _face_infos; ///@} }; - -// //* This file is part of the MOOSE framework -// //* https://www.mooseframework.org -// //* -// //* All rights reserved, see COPYRIGHT for full restrictions -// //* https://github.com/idaholab/moose/blob/master/COPYRIGHT -// //* -// //* Licensed under LGPL 2.1, please see LICENSE for details -// //* https://www.gnu.org/licenses/lgpl-2.1.html - -// #pragma once - -// #include "FVElementalKernel.h" -// #include "INSFVVelocityVariable.h" - -// /** -// * Computes source the sink terms for the turbulent kinetic energy. -// */ -// class INSFVTKESourceSink : public FVElementalKernel -// { -// public: -// static InputParameters validParams(); - -// virtual void initialSetup() override; - -// INSFVTKESourceSink(const InputParameters & parameters); - -// protected: -// ADReal computeQpResidual() override; - -// protected: -// /// The dimension of the domain -// const unsigned int _dim; - -// /// x-velocity -// const Moose::Functor & _u_var; -// /// y-velocity -// const Moose::Functor * _v_var; -// /// z-velocity -// const Moose::Functor * _w_var; - -// /// epsilon - dissipation rate of TKE -// const Moose::Functor & _epsilon; - -// /// Density -// const Moose::Functor & _rho; - -// /// Dynamic viscosity -// const Moose::Functor & _mu; - -// /// Turbulent dynamic viscosity -// const Moose::Functor & _mu_t; - -// /// Wall boundaries -// const std::vector & _wall_boundary_names; - -// /// Maximum mixing length allowed for the domain -// const Real _max_mixing_length; - -// /// Linearized model? -// const bool _linearized_model; - -// /// No equilibrium treatement -// const bool _non_equilibrium_treatment; - -// /// C_mu constant -// const Real _C_mu; - -// // Production Limiter Constant -// const Real _C_pl; - -// ///@{ -// /// Maps for wall treatement -// std::map _wall_bounded; -// std::map> _dist; -// std::map> _face_infos; -// ///@} -// }; diff --git a/modules/navier_stokes/src/auxkernels/RANSYPlusAux.C b/modules/navier_stokes/src/auxkernels/RANSYPlusAux.C index 8fa3d877becb..362a7ac7cd44 100644 --- a/modules/navier_stokes/src/auxkernels/RANSYPlusAux.C +++ b/modules/navier_stokes/src/auxkernels/RANSYPlusAux.C @@ -19,12 +19,12 @@ InputParameters RANSYPlusAux::validParams() { InputParameters params = AuxKernel::validParams(); - params.addClassDescription("Calculates y+ value"); + params.addClassDescription("Calculates non-dimensional wall distance (y+) value."); params.addRequiredCoupledVar("u", "The velocity in the x direction."); params.addCoupledVar("v", "The velocity in the y direction."); params.addCoupledVar("w", "The velocity in the z direction."); - params.addRequiredParam(NS::TKE, "Coupled turbulent kinetic energy."); - params.addRequiredParam(NS::density, "fluid density"); + params.addRequiredParam(NS::TKE, "Turbulent kinetic energy functor."); + params.addRequiredParam(NS::density, "Fluid density."); params.addRequiredParam(NS::mu, "Dynamic viscosity."); params.addParam>( "walls", {}, "Boundaries that correspond to solid walls."); @@ -33,7 +33,7 @@ RANSYPlusAux::validParams() wall_treatment, "The method used for computing the wall functions " "'eq_newton', 'eq_incremental', 'eq_linearized', 'neq'"); - params.addParam("C_mu", 0.09, "Coupled turbulent kinetic energy closure."); + params.addParam("C_mu", 0.09, "Coupled turbulent kinetic energy closure coefficient."); return params; } @@ -74,7 +74,6 @@ RANSYPlusAux::computeValue() const auto elem_arg = makeElemArg(_current_elem); const auto rho = _rho(elem_arg, state); const auto mu = _mu(elem_arg, state); - const auto TKE = _k(elem_arg, state); ADReal y_plus; if (_wall_bounded.find(_current_elem) != _wall_bounded.end()) @@ -99,19 +98,18 @@ RANSYPlusAux::computeValue() if (_wall_treatment == "neq") { // Non-equilibrium / Non-iterative - y_plus = std::pow(_C_mu, 0.25) * distance * std::sqrt(TKE) * rho / mu; + y_plus = std::pow(_C_mu, 0.25) * distance * std::sqrt(_k(elem_arg, state)) * rho / mu; } else { // Equilibrium / Iterative y_plus = NS::findyPlus(mu, rho, std::max(parallel_speed, 1e-10), distance); } + y_plus_vec.push_back(raw_value(y_plus)); } - // return *std::max_element(y_plus_vec.begin(), y_plus_vec.end()); - return raw_value(y_plus); + // Return average of y+ for cells with multiple wall faces + return 1.0 * std::accumulate(y_plus_vec.begin(),y_plus_vec.end(),0.0) / y_plus_vec.size(); } else - { - return 0.; - } + return 0.; } diff --git a/modules/navier_stokes/src/auxkernels/kEpsilonViscosityAux.C b/modules/navier_stokes/src/auxkernels/kEpsilonViscosityAux.C index 50826aadb8b0..b6dab8c0c2f5 100644 --- a/modules/navier_stokes/src/auxkernels/kEpsilonViscosityAux.C +++ b/modules/navier_stokes/src/auxkernels/kEpsilonViscosityAux.C @@ -36,7 +36,7 @@ kEpsilonViscosityAux::validParams() MooseEnum wall_treatment("eq_newton eq_incremental eq_linearized neq", "eq_newton"); params.addParam("wall_treatment", wall_treatment, - "The method used for computing the wall functions " + "The method used for computing the wall functions." "'eq_newton', 'eq_incremental', 'eq_linearized', 'neq'"); MooseEnum scale_limiter("none standard", "standard"); params.addParam("scale_limiter", @@ -76,15 +76,6 @@ kEpsilonViscosityAux::initialSetup() _wall_boundary_names, _c_fe_problem, _subproblem, blockIDs(), _face_infos); } } -/////////////////////////////////////////////// -// TO BE REMOVED DURING THE PR REVIEW PROCESS// -/////////////////////////////////////////////// -// Removed the method findUstar to make this file more compact and similar to -// the viscosity wall function bc. I also incorporated the initial guess for -// the Newton method to obtain Ustar in the findUstar method, previously present -// in this file, in the findUstar routine in NavierStokesMethods.C. -// If there was a specific reason I missed why it was there in the first place let's -// bring the findUstar method in this file back. Real kEpsilonViscosityAux::computeValue() @@ -126,7 +117,6 @@ kEpsilonViscosityAux::computeValue() // Switch for determining the near wall quantities // wall_treatment can be: "eq_newton eq_incremental eq_linearized neq" - ADReal u_tau; ADReal y_plus; ADReal mut_log; // turbulent log-layer viscosity ADReal mu_wall; // total wall viscosity to obtain the shear stress at the wall @@ -134,7 +124,7 @@ kEpsilonViscosityAux::computeValue() if (_wall_treatment == "eq_newton") { // Full Newton-Raphson solve to find the wall quantities from the law of the wall - u_tau = NS::findUStar(mu, rho, parallel_speed, min_wall_dist); + const auto u_tau = NS::findUStar(mu, rho, parallel_speed, min_wall_dist); y_plus = min_wall_dist * u_tau * rho / mu; mu_wall = rho * Utility::pow<2>(u_tau) * min_wall_dist / parallel_speed; mut_log = mu_wall - mu; @@ -155,7 +145,7 @@ kEpsilonViscosityAux::computeValue() (std::log(NS::E_turb_constant * std::max(min_wall_dist, 1.0) / mu) + 1.0); const ADReal c_c = parallel_speed; - u_tau = (-b_c + std::sqrt(std::pow(b_c, 2) + 4.0 * a_c * c_c)) / (2.0 * a_c); + const auto u_tau = (-b_c + std::sqrt(std::pow(b_c, 2) + 4.0 * a_c * c_c)) / (2.0 * a_c); y_plus = min_wall_dist * u_tau * rho / mu; mu_wall = rho * Utility::pow<2>(u_tau) * min_wall_dist / parallel_speed; mut_log = mu_wall - mu; @@ -169,13 +159,15 @@ kEpsilonViscosityAux::computeValue() std::log(std::max(NS::E_turb_constant * y_plus, 1 + 1e-4))); mut_log = mu_wall - mu; } + else + mooseAssert(false,"Should not reach here"); if (y_plus <= 5.0) // sub-laminar layer mu_t = 0.0; else if (y_plus >= 30.0) // log-layer - mu_t = std::max(mut_log.value(), 1e-12); + mu_t = std::max(mut_log.value(), NS::mu_t_low_limit); else { // buffer layer @@ -184,7 +176,7 @@ kEpsilonViscosityAux::computeValue() const auto mut_log = mu * (NS::von_karman_constant * 30.0 / std::log(std::max(NS::E_turb_constant * 30.0, 1 + 1e-4)) - 1.0); - mu_t = std::max(raw_value(blending_function * mut_log), 1e-12); + mu_t = std::max(raw_value(blending_function * mut_log), NS::mu_t_low_limit); } } else diff --git a/modules/navier_stokes/src/executioners/SIMPLENonlinearAssembly.C b/modules/navier_stokes/src/executioners/SIMPLENonlinearAssembly.C index d5db50f280d1..516bd162c7a0 100644 --- a/modules/navier_stokes/src/executioners/SIMPLENonlinearAssembly.C +++ b/modules/navier_stokes/src/executioners/SIMPLENonlinearAssembly.C @@ -479,6 +479,8 @@ SIMPLENonlinearAssembly::execute() size_t residual_index = 0; // Execute all objects tagged as nonlinear + // This will execute everything in the problem at nonlinear, including the aux kernels. + // This way we compute the aux kernels before the momentum equations are solved. _problem.execute(EXEC_NONLINEAR); // We clear the caches in the momentum and pressure variables diff --git a/modules/navier_stokes/src/fvbcs/INSFVTurbulentViscosityWallFunction.C b/modules/navier_stokes/src/fvbcs/INSFVTurbulentViscosityWallFunction.C index 97ebb55d16e2..36a48b4742fe 100644 --- a/modules/navier_stokes/src/fvbcs/INSFVTurbulentViscosityWallFunction.C +++ b/modules/navier_stokes/src/fvbcs/INSFVTurbulentViscosityWallFunction.C @@ -75,16 +75,14 @@ INSFVTurbulentViscosityWallFunction::boundaryValue(const FaceInfo & fi) const // Switch for determining the near wall quantities // wall_treatment can be: "eq_newton eq_incremental eq_linearized neq" - ADReal u_tau; ADReal y_plus; - const ADReal mut_visc = mu; // laminar sublayer viscosity ADReal mut_log; // turbulent log-layer viscosity ADReal mu_wall; // total wall viscosity to obtain the shear stress at the wall if (_wall_treatment == "eq_newton") { // Full Newton-Raphson solve to find the wall quantities from the law of the wall - u_tau = NS::findUStar(mu, rho, parallel_speed, wall_dist); + const auto u_tau = NS::findUStar(mu, rho, parallel_speed, wall_dist); y_plus = wall_dist * u_tau * rho / mu; mu_wall = rho * Utility::pow<2>(u_tau) * wall_dist / parallel_speed; mut_log = mu_wall - mu; @@ -105,7 +103,7 @@ INSFVTurbulentViscosityWallFunction::boundaryValue(const FaceInfo & fi) const (std::log(NS::E_turb_constant * std::max(wall_dist, 1.0) / mu) + 1.0); const ADReal c_c = parallel_speed; - u_tau = (-b_c + std::sqrt(std::pow(b_c, 2) + 4.0 * a_c * c_c)) / (2.0 * a_c); + const auto u_tau = (-b_c + std::sqrt(std::pow(b_c, 2) + 4.0 * a_c * c_c)) / (2.0 * a_c); y_plus = wall_dist * u_tau * rho / mu; mu_wall = rho * Utility::pow<2>(u_tau) * wall_dist / parallel_speed; mut_log = mu_wall - mu; @@ -119,12 +117,15 @@ INSFVTurbulentViscosityWallFunction::boundaryValue(const FaceInfo & fi) const std::log(std::max(NS::E_turb_constant * y_plus, 1 + 1e-4))); mut_log = mu_wall - mu; } + else + mooseAssert(false,"Should not reach here"); + if (y_plus <= 5.0) // sub-laminar layer return 0.0; else if (y_plus >= 30.0) // log-layer - return std::max(mut_log, 1e-12); + return std::max(mut_log, NS::mu_t_low_limit); else { // buffer layer @@ -133,6 +134,6 @@ INSFVTurbulentViscosityWallFunction::boundaryValue(const FaceInfo & fi) const const auto mut_log = mu * (NS::von_karman_constant * 30.0 / std::log(std::max(NS::E_turb_constant * 30.0, 1 + 1e-4)) - 1.0); - return blending_function * std::max(mut_log, 1e-12); + return blending_function * std::max(mut_log, NS::mu_t_low_limit); } } diff --git a/modules/navier_stokes/src/fvkernels/INSFVTKEDSourceSink.C b/modules/navier_stokes/src/fvkernels/INSFVTKEDSourceSink.C index e27e39a8a20f..3a53e5638eec 100644 --- a/modules/navier_stokes/src/fvkernels/INSFVTKEDSourceSink.C +++ b/modules/navier_stokes/src/fvkernels/INSFVTKEDSourceSink.C @@ -118,7 +118,7 @@ INSFVTKEDSourceSink::computeQpResidual() if (_wall_treatment == "neq") { // Non-equilibrium / Non-iterative - y_plus = std::pow(_C_mu, 0.25) * distance * std::sqrt(TKE) * rho / mu; + y_plus = distance * std::sqrt(std::sqrt(_C_mu) * TKE) * rho / mu; } else { @@ -145,14 +145,11 @@ INSFVTKEDSourceSink::computeQpResidual() const Elem * const loc_elem = defined_on_elem_side ? &fi->elem() : fi->neighborPtr(); const Moose::FaceArg facearg = { fi, Moose::FV::LimiterType::CentralDifference, false, false, loc_elem}; - const ADReal wall_mu = _mu(facearg, state); - destruction += 2.0 * TKE * wall_mu / rho / Utility::pow<2>(distance_vec[i]) / tot_weight; + destruction += 2.0 * TKE * _mu(facearg, state) / rho / Utility::pow<2>(distance_vec[i]) / tot_weight; } else - { destruction += std::pow(_C_mu, 0.75) * rho * std::pow(TKE, 1.5) / (NS::von_karman_constant * distance_vec[i]) / tot_weight; - } } residual = _var(makeElemArg(_current_elem), state) - destruction; diff --git a/modules/navier_stokes/src/fvkernels/INSFVTKESourceSink.C b/modules/navier_stokes/src/fvkernels/INSFVTKESourceSink.C index 15f7897f4064..22b0d35526a4 100644 --- a/modules/navier_stokes/src/fvkernels/INSFVTKESourceSink.C +++ b/modules/navier_stokes/src/fvkernels/INSFVTKESourceSink.C @@ -118,7 +118,7 @@ INSFVTKESourceSink::computeQpResidual() if (_wall_treatment == "neq") { // Non-equilibrium / Non-iterative - y_plus = std::pow(_C_mu, 0.25) * distance * std::sqrt(_var(elem_arg, old_state)) * rho / mu; + y_plus = distance * std::sqrt(std::sqrt(_C_mu) * _var(elem_arg, old_state)) * rho / mu; } else { @@ -171,9 +171,7 @@ INSFVTKESourceSink::computeQpResidual() const auto tau_w = (wall_mut + wall_mu) * velocity_grad_norm_vec[i]; if (y_plus < 11.25) - { destruction += destruction_visc; - } else { destruction += destruction_log;