Merge pull request #27889 from freiler/k-epsilon-mods

K epsilon mods

modules/navier_stokes/doc/content/bib/navier_stokes.bib

@@ -467,3 +467,23 @@ @article{hibiki2000interface
   year={2000},
   publisher={Elsevier}
 }
+@article{durbin1996k,
+  title={On the k-3 stagnation point anomaly},
+  author={Durbin, Paul A},
+  journal={International journal of heat and fluid flow},
+  volume={17},
+  number={1},
+  pages={89--90},
+  year={1996},
+  publisher={Elsevier}
+}
+
+@article{menter1994two,
+  title={Two-equation eddy-viscosity turbulence models for engineering applications},
+  author={Menter, Florian R},
+  journal={AIAA journal},
+  volume={32},
+  number={8},
+  pages={1598--1605},
+  year={1994}
+}

modules/navier_stokes/doc/content/source/auxkernels/RANSYPlusAux.md

@@ -0,0 +1,13 @@
+# RANSYPlusAux
+
+Computes the dimensionless wall distance $y^+$. 
+
+Four different formulations are supported for computing $y^+$ 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).
+
+!syntax parameters /AuxKernels/RANSYPlusAux
+
+!syntax inputs /AuxKernels/RANSYPlusAux
+
+!syntax children /AuxKernels/RANSYPlusAux

modules/navier_stokes/doc/content/source/auxkernels/kEpsilonViscosityAux.md

@@ -3,7 +3,7 @@
 This is the auxiliary kernel used to compute the dynamic turbulent viscosity
 
 \begin{equation}
-\mu_t = \rho C_{\mu} \frac{k^2}{\epsilon} \,,
+\mu_t = \rho C_{\mu} k T_e \,,
 \end{equation}
 
 where:
@@ -12,6 +12,7 @@ where:
 - $C_{\mu} = 0.09$ is a closure parameter,
 - $k$ is the turbulent kinetic energy,
 - $\epsilon$ is the turbulent kinetic energy dissipation rate.
+- $T_e = max( \frac{k}{\epsilon} , \sqrt(\frac{\nu}{\epsilon}) )$.
 
 By setting parameter [!param](/AuxKernels/kEpsilonViscosityAux/bulk_wall_treatment) to `true`, the
 kernel allows us to set the value of the cells on the boundaries specified in

modules/navier_stokes/doc/content/source/fvbcs/INSFVTurbulentViscosityWallFunction.md

@@ -12,8 +12,29 @@ boundary layer are identified as follows:
 - Buffer region: $y^+ \in (5, 30)$
 - Logarithmic region: $y^+ \ge 30$
 
-Four different formulations are supported
-as defined by the [!param](/FVBCs/INSFVTurbulentViscosityWallFunction/wall_treatment) parameter.
+The wall function goal is to set the total viscosity at the wall $\mu_w$, decomposed as
+$\mu_w = \mu + \mu_t$, such that the wall shear stress $\tau_w$ is accurately captured 
+without the need of fully resolving the gradients at the near wall region. 
+
+\begin{equation}
+    \tau_w = \frac{ \mu_w u_p}{y_p} \,,
+\end{equation}
+
+where:
+
+- $\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 dynamic 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 [RANSYPlusAux.md] auxiliary kernel. 
+to estimate the dimensionless wall distance $y^+$.
 
 ## Equilibrium wall functions using a Newton solve
 
@@ -26,21 +47,22 @@ for the turbulent viscosity.
     \mu_t =
     \begin{cases}
         0 & \text{if } y^+ \le 5 \\
-        \frac{\rho u_{\tau}^2 y_p}{u_p} & \text{if } y^+ \ge 30
+        \frac{\rho u_{\tau}^2 y_p}{u_p} - \mu & \text{if } y^+ \ge 30 \,,
     \end{cases}
 \end{equation}
 
 where:
 
 - $\rho$ is the density
+- $\mu$ is the dynamic 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
 
 For the buffer layer, a linear blending method is used that defines the turbulent viscosity as follows:
 
 \begin{equation}
-    \mu_t = \frac{\rho u_{\tau}^2 y_p}{u_p} \frac{(y^+ - 5)}{25}
+    \mu_t = \mu_t(y^+=30) \frac{(y^+ - 5)}{25}
 \end{equation}
 
 Note that for $y^+ = 5$ and $y^+ = 30$ we recover the sub-laminar and logarithmic profiles, respectively.
@@ -133,7 +155,7 @@ Then, the turbulent viscosity is defined as follows:
 For the buffer layer, a linear blending method is used that defines the turbulent viscosity as follows:
 
 \begin{equation}
-    \mu_t = \mu \left[ \frac{\kappa y^+}{\operatorname{ln}(E y^+)} - 1.0 \right] \frac{(y^+ - 5)}{25}
+    \mu_t = \mu_t(y^+=30) \frac{(y^+ - 5)}{25}
 \end{equation}
 
 !alert note

modules/navier_stokes/doc/content/source/fvkernels/INSFVTKEDSourceSink.md

@@ -8,19 +8,16 @@ A different treatment is used for the bulk and the near wall regions.
 
 ## Bulk formulation:
 
-The turbulent production $G_k$ is modeled as follows:
+The production of turbulent kinetic energy dissipation $G_\epsilon$ is modeled as follows:
 
 \begin{equation}
-G_{\epsilon} = C_{1,\epsilon} \rho C_{\mu} S^2 k \,,
+G_{\epsilon} = C_{1,\epsilon} \frac{\epsilon}{k} G_k \,,
 \end{equation}
 
 where:
 
 - $C_{1,\epsilon} = 1.44$ is a closure parameter,
-- $\rho$ is the density
-- $C_{\mu} = 0.09$ is another closure parameter,
-- $S$ is the shear strain tensor internal norm, defined as $S = \sqrt{2\mathbf{S}:\mathbf{S}}$ with the shear strain tensor defined as $\mathbf{S} = \frac{1}{2} [\nabla \vec{u} + (\nabla \vec{u})^T]$,
-- $k$ is the turbulent kinetic energy, which generally comes from a coupled transport equation or can be set by the user for reproducing canonical cases.
+- $G_k$ is the limited turbulent kinetic energy production. For more details please refer to [INSFVTKESourceSink](INSFVTKESourceSink.md).
 
 The destruction of the turbulent kinetic energy dissipation rate is modeled as follows:
 
@@ -32,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:
@@ -41,12 +39,12 @@ treatment in which the equilibrium value for the $\epsilon = \epsilon_{eq}$ is s
 A separate formulation is used for the `sub-laminar` and `logarithmic` boundary layers.
 The determination of whether the near-wall cell lies in the laminar or logarithmic region
 is performed via the non-dimensional wall distance $y^+$.
-The non-dimensional wall distance can be as defined differently according to the
-[!param](/FVKernels/INSFVTKEDSourceSink/non_equilibrium_treatment) parameter.
+The non-dimensional wall distance can be defined differently according to the
+[!param](/FVKernels/INSFVTKEDSourceSink/wall_treatment) parameter. 
 
-If [!param](/FVKernels/INSFVTKEDSourceSink/non_equilibrium_treatment) is `false`, the
-standard wall function formulations is used in
-which $y^+$ is found by an incremental fixed point algorithm as follows:
+The four formulations are described in more detail in [INSFVTurbulentViscosityWallFunction.md]. 
+
+If an equilibrium [!param](/FVKernels/INSFVTKEDSourceSink/wall_treatment) is defined, i.e. `eq_newton`,`eq_incremental` or `eq_linearized`, the standard wall function formulations are used in which $y^+$ is found:
 
 \begin{equation}
 y^+ = \frac{\rho y_p u_{\tau}}{\mu} \,,
@@ -59,11 +57,11 @@ where:
 - $u_{\tau}$ is the friction velocity, defined as $u_{\tau} = \sqrt{\frac{\tau_w}{\rho}}$ with $\tau_w$ the shear stress at the wall for which the condition is applied,
 - $\mu$ is the dynamic molecular viscosity.
 
-If [!param](/FVKernels/INSFVTKEDSourceSink/non_equilibrium_treatment) is `true`,
-non-equilibrium wall function formulations is used in which the $y^+$ is defined as follows:
+If a non-equilibrium [!param](/FVKernels/INSFVTKEDSourceSink/wall_treatment) is defined, i.e. `neq`,
+the $y^+$ is defined non-iteratively as follows:
 
 \begin{equation}
-y^+ = \frac{C_{\mu}^{0.25} y_p \sqrt{k}}{\mu} \,,
+y^+ = \frac{y_p \sqrt{\sqrt{C_{\mu}}k}}{\mu} \,,
 \end{equation}
 
 !alert note
@@ -76,7 +74,7 @@ A different value is used for $\epsilon_{eq}$ in each of the two regions.
 For the `sub-laminar` boundary layer, the equilibrium value is determined as follows:
 
 \begin{equation}
-\epsilon_{eq} = 2 \frac{k \mu_t}{y_p^2}\,,
+\epsilon_{eq} = 2 \frac{k \mu}{y_p^2}\,,
 \end{equation}
 
 where:

modules/navier_stokes/doc/content/source/fvkernels/INSFVTKESourceSink.md

@@ -24,6 +24,16 @@ transport equation for $\epsilon$.
 However, for canonical or measured cases, e.g., isotropic decaying turbulence,
 the user can utilize predefined fields through functors in MOOSE.
 
+To avoid the overproduction of turbulent kinetic energy in stagnation zones \cite{durbin1996k}, a production limiter is imposed in relation to the dissipation using the formulation in \cite{menter1994two}:
+
+\begin{equation}
+G_k = min \left( G_k , C_{PL} \rho \epsilon \right) \,,
+\end{equation}
+
+where:
+
+- $C_{PL}$ it the limiter constant, and set by default to a recommended value of 10 \cite{durbin1996k}.
+
 ## Wall formulation:
 
 All cells in contact with a boundary identified in the [!param](/FVKernels/INSFVTKESourceSink/walls) list are applied a different
@@ -49,11 +59,19 @@ 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$:
+
+\begin{equation}
+G_k = 0.0 \,,
+\end{equation}
 
-For both the `sub-laminar` and `logarithmic` boundary layers the production term is defined as:
+In the `logarithmic` boundary layers the production term is no longer negligible and is defined as:
 
 \begin{equation}
-G_k = C_{\mu}^{0.25} \mu_t \frac{k^{\frac{3}{2}} ||\nabla \vec{u}||}{\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:
@@ -63,12 +81,12 @@ where:
 - $||\nabla \vec{u}||$ is the near wall velocity gradient norm, which is defined as $||\nabla \vec{u}|| = (\nabla \vec{u} \cdot \hat{n}) \cdot \hat{n}$,
 - $\kappa = 0.41$ is the von Kármán constant.
 
-The formulation assumes that the near wall value is already imposed in the $\mu_t$ functor.
+The formulation assumes that the near wall value is already imposed in the $\mu_t$ functor. 
 
 When solving a linear problem, instead of the nonlinear formulation, the production term is formulated as:
 
 \begin{equation}
-G_k = C_{\mu}^{0.25} \mu_t \frac{k \sqrt{k_{old}}||\nabla \vec{u}||}{\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:
@@ -82,7 +100,7 @@ For the destruction, formulation is different for the `sub-laminar` and `logarit
 For the `sub-laminar` layer, the destruction is defined as follows:
 
 \begin{equation}
-\epsilon = \frac{2 \mu_t k}{y_p ^2} \,.
+\epsilon = \frac{2 \mu k}{y_p ^2} \,.
 \end{equation}
 
 For the `logarithmic` layer, the destruction is defined as follows:

modules/navier_stokes/include/auxkernels/RANSYPlusAux.h

@@ -0,0 +1,65 @@
+//* 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"
+#include "NS.h"
+
+/**
+ * Computes wall y+ based on wall functions.
+ */
+class RANSYPlusAux : public AuxKernel
+{
+public:
+  static InputParameters validParams();
+
+  virtual void initialSetup() override;
+
+  RANSYPlusAux(const InputParameters & parameters);
+
+protected:
+  virtual Real computeValue() override;
+
+  /// the dimension of the simulation
+  const unsigned int _dim;
+
+  /// x-velocity
+  const Moose::Functor<ADReal> & _u_var;
+  /// y-velocity
+  const Moose::Functor<ADReal> * _v_var;
+  /// z-velocity
+  const Moose::Functor<ADReal> * _w_var;
+
+  /// Turbulent kinetic energy
+  const Moose::Functor<ADReal> * _k;
+
+  /// Density
+  const Moose::Functor<ADReal> & _rho;
+
+  /// Dynamic viscosity
+  const Moose::Functor<ADReal> & _mu;
+
+  /// Wall boundary names
+  const std::vector<BoundaryName> & _wall_boundary_names;
+
+  /// Method used for wall treatment
+  NS::WallTreatmentEnum _wall_treatment;
+
+  /// C_mu constant
+  const Real _C_mu;
+
+  ///@{
+  /// Maps for wall treatement
+  std::map<const Elem *, bool> _wall_bounded;
+  std::map<const Elem *, std::vector<Real>> _dist;
+  std::map<const Elem *, std::vector<const FaceInfo *>> _face_infos;
+  ///@}
+};

modules/navier_stokes/include/auxkernels/kEpsilonViscosityAux.h

@@ -1,5 +1,4 @@
 
-
 //* This file is part of the MOOSE framework
 //* https://www.mooseframework.org
 //*
@@ -13,6 +12,7 @@
 
 #include "AuxKernel.h"
 #include "INSFVVelocityVariable.h"
+#include "NS.h"
 
 /**
  * Computes the turbuent viscosity for the k-Epsilon model.
@@ -30,10 +30,6 @@ class kEpsilonViscosityAux : public AuxKernel
 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;
 
@@ -57,17 +53,20 @@ class kEpsilonViscosityAux : public AuxKernel
   /// C-mu closure coefficient
   const Real _C_mu;
 
+  // Maximum allowable mu_t_ratio : mu/mu_t
+  const Real _mu_t_ratio_max;
+
   /// Wall boundaries
   const std::vector<BoundaryName> & _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;
+  /// Method used for wall treatment
+  NS::WallTreatmentEnum _wall_treatment;
+
+  /// Method used to limit the k-e time scale
+  const MooseEnum _scale_limiter;
 
   // -- Parameters of the wall function method
 

modules/navier_stokes/include/base/NS.h

@@ -167,9 +167,22 @@ static const std::string mass_flux = "mass_flux";
 static const std::string TKE = "k";
 static const std::string TKED = "epsilon";
 
+/**
+ * Wall treatment options
+ */
+enum class WallTreatmentEnum
+{
+  EQ_NEWTON = 0,
+  EQ_INCREMENTAL = 1,
+  EQ_LINEARIZED = 2,
+  NEQ = 3
+};
+
 // 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-8;
 }
 
 namespace NS_DEFAULT_VALUES