Merge pull request #25752 from tanoret/k_epsilon_polished

K epsilon model

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

@@ -0,0 +1,23 @@
+# TurbulentConductivityAux
+
+!syntax description /AuxKernels/TurbulentConductivityAux
+
+## Overview
+
+This is the auxiliary kernel used to compute the thermal effective turbulent conductivity
+
+\begin{equation}
+k_t = \frac{c_p \mu_t}{Pr_t} \,,
+\end{equation}
+
+where:
+
+- $c_p$ is the specific heat at constant pressure,
+- $\mu_t$ is the dynamic turbulent viscosity,
+- $Pr_t$ is the turbulent Prandtl number, which usually ranges between 0.3 and 0.9.
+
+!syntax parameters /AuxKernels/TurbulentConductivityAux
+
+!syntax inputs /AuxKernels/TurbulentConductivityAux
+
+!syntax children /AuxKernels/TurbulentConductivityAux

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

@@ -0,0 +1,34 @@
+# kEpsilonViscosityAux
+
+This is the auxiliary kernel used to compute the dynamic turbulent viscosity
+
+\begin{equation}
+\mu_t = \rho C_{\mu} \frac{k^2}{\epsilon} \,,
+\end{equation}
+
+where:
+
+- $\rho$ is the density,
+- $C_{\mu} = 0.09$ is a closure parameter,
+- $k$ is the turbulent kinetic energy,
+- $\epsilon$ is the turbulent kinetic energy dissipation rate.
+
+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
+[!param](/AuxKernels/kEpsilonViscosityAux/walls) to the dynamic turbulent viscosity predicted
+from the law of the wall or non-equilibrium wall functions.
+See [INSFVTurbulentViscosityWallFunction](INSFVTurbulentViscosityWallFunction.md) for more
+details about the near-wall implementation.
+
+!alert note
+If the boundary conditions for the dynamic turbulent viscosity are already set via [INSFVTurbulentViscosityWallFunction](INSFVTurbulentViscosityWallFunction.md),
+there is no need to add bulk wall treatment, i.e., we should set
+[!param](/AuxKernels/kEpsilonViscosityAux/bulk_wall_treatment) to `false`.
+This type of bulk wall treatment is mainly designed for porous media formulations
+with large computational cells.
+
+!syntax parameters /AuxKernels/kEpsilonViscosityAux
+
+!syntax inputs /AuxKernels/kEpsilonViscosityAux
+
+!syntax children /AuxKernels/kEpsilonViscosityAux

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

@@ -0,0 +1,22 @@
+# INSFVInletIntensityTKEBC
+
+This object wraps [`FVFunctionDirichletBC`](FVFunctionDirichletBC.md),
+to impose a precomputed value for the turbulent kinetic energy.
+
+The value set for the turbulent kinetic energy is:
+
+\begin{equation}
+k = \frac{3}{2} (I |\vec{u} \cdot \vec{n}|)^2 \,,
+\end{equation}
+
+where:
+
+- $I$ is the turbulent intensity, which can be set by the user or computed via correlations, e.g., $I = 0.16 Re^{-\frac{1}{8}}$
+- $\vec{u}$ is the velocity vector at a bounary face,
+- $\vec{n}$ is the normal vector for a boundary face.
+
+!syntax parameters /FVBCs/INSFVInletIntensityTKEBC
+
+!syntax inputs /FVBCs/INSFVInletIntensityTKEBC
+
+!syntax children /FVBCs/INSFVInletIntensityTKEBC

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

@@ -0,0 +1,22 @@
+# INSFVMixingLengthTKEDBC
+
+This object wraps [`FVFunctionDirichletBC`](FVFunctionDirichletBC.md),
+to impose a precomputed value for the turbulent kinetic energy dissipation rate.
+
+The value set for the turbulent kinetic energy is:
+
+\begin{equation}
+\epsilon = C_{\mu}^{0.75} \frac{k^{\frac{3}{2}}}{0.07 L} \,,
+\end{equation}
+
+where:
+
+- $C_{\mu} = 0.09$ is a closure parameter,
+- $k$ is the turbulent kinetic energy,
+- $L$ is a characterisitc length, e.g., the diameter of a pipe, diameter of an inlet jet, or the side length of the lid-driven cavity.
+
+!syntax parameters /FVBCs/INSFVMixingLengthTKEDBC
+
+!syntax inputs /FVBCs/INSFVMixingLengthTKEDBC
+
+!syntax children /FVBCs/INSFVMixingLengthTKEDBC

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

@@ -0,0 +1,12 @@
+# INSFVTKEDWallFunctionBC
+
+This boundary condtion should only be used if no wall treatment is added.
+Implements wall boundary conditions for the turbulent kinetic energy dissipation rate.
+A separate treatment is used for the laminar and logarithmic layers.
+A linear blending functions is used to interpolate between both layers.
+
+!syntax parameters /FVBCs/INSFVTKEDWallFunctionBC
+
+!syntax inputs /FVBCs/INSFVTKEDWallFunctionBC
+
+!syntax children /FVBCs/INSFVTKEDWallFunctionBC

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

@@ -0,0 +1,86 @@
+# INSFVTurbulentTemperatureWallFunction
+
+The function sets up the equivalent heat flux of the near-wall
+boundary layer when the wall temperature is set by the `T_w` parameter.
+
+The boundary conditions are different depending on whether the centroid
+of the cell near the identified boundary lies in the wall function profile.
+Taking the non-dimensional wall distance as $y^+$, the three regions of the
+boundary layer are identified as follows:
+
+- Sub-laminar region: $y^+ \le 5$
+- Buffer region: $y^+ \in (5, 30)$
+- Logarithmic region: $y^+ \ge 30$
+
+For the procedure of determining the non-dimensional wall distance as $y^+$
+and the friction velocity $u_{\tau}$ plese see
+[INSFVTurbulentViscosityWallFunction](INSFVTurbulentViscosityWallFunction.md).
+
+For the sub-laminar and logarithmic layer, the thermal diffusivity is defined
+as follows:
+
+\begin{equation}
+    \alpha =
+    \begin{cases}
+        \alpha_l = \frac{k}{\rho c_p} & \text{if } y^+ \le 5 \\
+        \alpha_t = \frac{u_{\tau} y_p}{Pr_t w_s} & \text{if } y^+ \ge 30
+    \end{cases}
+\end{equation}
+
+where:
+
+- $k$ is the thermal conductivity
+- $\rho$ is the density
+- $c_p$ is the specific heat at constant pressure
+- $u_{\tau}$ is the friction velocity defined by law of the wall
+- $y_p$ is the distance from the boundary to the centroid of the near-wall cell
+- $Pr_t$ is the turbulent Prandtl number, which typically ranges between 0.3 and 0.9
+- $w_s$ is a near-wall scaling factor that is defined as follows:
+
+\begin{equation}
+  w_s = \frac{1}{\kappa} \operatorname{ln}(E y^+) + J_k \,,
+\end{equation}
+
+where:
+
+- $\kappa = 0.4187$ is the von Kármán constant
+- $E = 9.793$ is a closure parameter
+- $J_k$ is the Jayatilleke wall functions defined as follows:
+
+\begin{equation}
+  J_k = 9.24 \left[ \left( \frac{Pr}{Pr_t} \right)^{0.75} - 1 \right] \left( 1 + 0.28 e^{-0.007 \frac{Pr}{Pr_t}} \right) \,,
+\end{equation}
+
+where:
+
+- $Pr$ is the Prandtl number
+
+For the buffer layer, i.e., in $y^+ \in (5, 30)$, the thermal diffusivity
+is defined via a linear blending function as follows:
+
+\begin{equation}
+  \alpha_b = \alpha_t \frac{(y^+ - 5)}{25} + \alpha_l \left[ 1 - \frac{(y^+ - 5)}{25} \right]
+\end{equation}
+
+Finally, using the thermal diffusivity, the heat flux at the wall is defined
+as follows:
+
+\begin{equation}
+  q''_w = - \rho c_p \alpha \frac{T_w - T_p}{y_p} \,,
+\end{equation}
+
+where:
+
+- $T_p$ is the temperature at the centroid of the near-wall cell
+- $y_p$ is the distance from the wall to the centroid of the near-wall cell
+
+!alert note
+The thermal wall functions are only valid for regions in which the equilibrium
+momentum wall functions are valid, i.e., no flow detachment, recirculation, etc.
+For resolving non-equilibrium phenomena, we recommend refining the mesh.
+
+!syntax parameters /FVBCs/INSFVTurbulentTemperatureWallFunction
+
+!syntax inputs /FVBCs/INSFVTurbulentTemperatureWallFunction
+
+!syntax children /FVBCs/INSFVTurbulentTemperatureWallFunction

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

@@ -0,0 +1,147 @@
+# INSFVTurbulentViscosityWallFunction
+
+Implements wall function boundary condition for the turbulent
+dynamic viscosity $\mu_t$.
+
+The boundary conditions are different depending on where the centroid
+of the cell near the identified boundary lies in the wall function profile.
+Taking the non-dimensional wall distance as $y^+$, the three regions of the
+boundary layer are identified as follows:
+
+- Sub-laminar region: $y^+ \le 5$
+- 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.
+
+## Equilibrium wall functions using a Newton solve
+
+This treatment can be enabled by setting the parameter
+[!param](/FVBCs/INSFVTurbulentViscosityWallFunction/wall_treatment) to `eq_newton`.
+The treatment uses equilibrium wall functions where the following formulation is used
+for the turbulent viscosity.
+
+\begin{equation}
+    \mu_t =
+    \begin{cases}
+        0 & \text{if } y^+ \le 5 \\
+        \frac{\rho u_{\tau}^2 y_p}{u_p} & \text{if } y^+ \ge 30
+    \end{cases}
+\end{equation}
+
+where:
+
+- $\rho$ is the density
+- $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}
+\end{equation}
+
+Note that for $y^+ = 5$ and $y^+ = 30$ we recover the sub-laminar and logarithmic profiles, respectively.
+
+Here the standard or equilibrium law of the wall defines $y^+$ and $u_{\tau}$ as follows:
+
+\begin{equation}
+  \frac{u_p}{u_{\tau}} = \frac{1}{\kappa} \operatorname{ln}(E y^+) \,,
+\end{equation}
+
+\begin{equation}
+  y^+ = \frac{\rho u_{\tau} y_p}{\mu} \,,
+\end{equation}
+
+where:
+
+- $\mu$ is the molecular dynamic viscosity
+- $E = 9.793$ is a closure parameter
+- $\kappa = 0.4187$ is the von Kármán constant
+
+In this method, we iterate on the wall function and $y^+$ to find
+$u_{\tau}$ via a Newton solve. Once $u_{\tau}$ is defined, $y^+$ is
+computed followed by the determination of the boundary turbulent viscosity.
+
+!alert note
+`eq_newton` solve will converge the fastest for simple flow geometries but it
+may diverge for more complicated flows. Also, the code will run if the center
+of the near wall cells are in the buffer layer. However, using a mesh that
+contains nodes in the buffer layer is not recommended.
+
+
+## Equilibrium wall functions using a fixed-point solve
+
+This treatment is enabled by setting parameter
+[!param](/FVBCs/INSFVTurbulentViscosityWallFunction/wall_treatment) to `eq_incremental`.
+The method uses the same equilibrium wall treatement than the Newton solve.
+However, the main difference is that, instead of computing $u_{\tau}$ for the
+near wall cells, a fixed point iteration is performed in the wall functions
+to find $y^+$.
+
+!alert note
+`eq_incremental` has a larger convergence radius than the Newton solve and
+internal controls are added to avoid issues converging the wall function
+at the buffer layer. However, it will take more iterations than the Newton
+solve to converge. Using a mesh that contains nodes in the buffer layer is
+not recommended.
+
+
+## Equilibrium wall functions using linearized wall function
+
+This treatment is enabled by setting parameter
+[!param](/FVBCs/INSFVTurbulentViscosityWallFunction/wall_treatment) to `eq_linearized`.
+The treatement uses a linearized version of the wall function, in which
+a linear Taylor approximation is used for the natural logarithm.
+This approximation results in a quadratic equation that is solved directly for $u_{\tau}$.
+Then, $y^+$ is computed from $u_{\tau}$.
+
+!alert note
+`eq_linearized` will work fast as there is no nonlinear solve at
+the near-wall region. However, the method may introduce significant
+near-wall errors. The method is designed to be used in conjunction
+with porous media treatement and not necessairly for free flow.
+
+## Non-equilibrium wall functions
+
+This treatment is enabled by setting parameter
+[!param](/FVBCs/INSFVTurbulentViscosityWallFunction/wall_treatment) to `neq`.
+In this case, the non-dimensional wall distance is computed from the
+turbulent kinetic energy near the wall as follows:
+
+\begin{equation}
+  y^+ = \frac{\rho y_p C_{\mu}^{0.25} \sqrt{k_p}}{\mu} \,,
+\end{equation}
+
+where:
+
+- $C_{\mu} = 0.09$ is a fitting parameter
+- $k_p$ is the turbulent kinetic energy at the centroid of the near-wall cell
+
+Then, the turbulent viscosity is defined as follows:
+
+\begin{equation}
+    \mu_t =
+    \begin{cases}
+        0 & \text{if } y^+ \le 5 \\
+        \mu \left[ \frac{\kappa y^+}{\operatorname{ln}(E y^+)} - 1.0 \right] & \text{if } y^+ \ge 30
+    \end{cases}
+\end{equation}
+
+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}
+\end{equation}
+
+!alert note
+`neq` should mainly be used for dettached flow or other cases for which equilibrium wall
+functions are not valid. One should try to use equilibrium wall functions when possible.
+
+!syntax parameters /FVBCs/INSFVTurbulentViscosityWallFunction
+
+!syntax inputs /FVBCs/INSFVTurbulentViscosityWallFunction
+
+!syntax children /FVBCs/INSFVTurbulentViscosityWallFunction

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

@@ -9,6 +9,19 @@ Navier-Stokes momentum equation, e.g.
 
 where $\mu$ is the dynamic viscosity, and $\bm{v}$ is the velocity.
 
+The object also takes a parameter
+[!param](/FVKernels/INSFVMomentumDiffusion/complete_expansion) which is
+`false` by default. If [!param](/FVKernels/INSFVMomentumDiffusion/complete_expansion)
+is activated, the following complete formulation is used for the momentum viscous stress:
+
+\begin{equation}
+- \left[ \nabla \cdot \mu \left( \nabla \bm{v} +  (\nabla \bm{v})^T \right) \right]
+\end{equation}
+
+!alert note
+The term $\nabla \cdot \mu (\nabla \bm{v})^T = 0$ for incompressible flow if a constant
+dynamic viscosity is used.
+
 !syntax parameters /FVKernels/INSFVMomentumDiffusion
 
 !syntax inputs /FVKernels/INSFVMomentumDiffusion