## Turbulence Models in OpenFOAM – Hybrid Methods

Last Updated: May 4, 2019

Large Eddy Simulation (LES) can be used to accurately simulate unsteady flow behaviors but its numerical cost is higher compared to the conventional Reynolds Averaged Navier-Stokes (RANS) approach.

Several methods have been developed in order to save computational time and cost for unsteady flow computations compared to LES.

Keywords
Detached-Eddy Simulation (DES), Scale Adaptive Simulation (SAS), Grey area problem

##### Detached Eddy Simulation (DES)

Detached-Eddy Simulation (DES) is a hybrid Reynolds Averaged Navier-Stokes/Large-Eddy Simulation model. The DES model was first proposed in 1997 by Spalart et al. [1] based on the Spalart-Allmaras RANS model and it is commonly referred to as DES97.

Spalart-Allmaras Based DES Formulation (DES97)

Its formulation is briefly described in [2]:

The driving length scale of the RANS S-A model is the distance to the closest wall, $$d$$. This makes a modification to this model for DES mode quite straightforward (exactly for this reason it was used as a basis of DES in the first publication [1]). The modification consists in substituting for $$d$$, everywhere in the equations, the new DES length scale, $$\tilde{d}$$. This length is also based on the grid spacing $$\Delta$$ and is defined as:

\tilde{d} = {\rm min}\left( d, C_{DES}\Delta \right), \tag{1} \label{eq:dTilda}

where $$C_{DES}$$ is the only new adjustable model constant, and $$\Delta$$ is based on the largest dimension of the local grid cell

\Delta = {\rm max}\left( \delta_{x}, \delta_{y}, \delta_{z} \right). \tag{2} \label{eq:delta}

Here we assume for simplicity that the grid is structured and that the coordinates $$\left(x, y, z\right)$$ are aligned with the grid cell, but the generalizations are obvious.

For wall-bounded separated flows, the above formulation results in a bybrid model that functions as the standard RANS S-A model inside the whole attached boundary layer, and as its subgrid-scale version in the rest of the flow including the separated regions and near wake. Indeed, in the attached boundary layer, due to the significant grid anisotropy $$\left(\delta_{x} \approx \delta_{z} \gg \delta_{y}\right)$$ typical of this flow region, in accordance with \eqref{eq:dTilda}, $$\tilde{d} = d$$, and the model reduces to the standard S-A RANS model. Otherwise, once a field point is far enough from walls $$\left( d > C_{DES}\Delta \right)$$, the length scale of the model becomes grid-dependent, i.e., the model performs as a subgrid-scale version of the S-A model. Note that at “equilibrium” (meaning a balance of production and destruction terms) this model reduces to an algebraic miximg-length Smagorinski-like subgrid model.

The “DES limiter” defined by Eq. \eqref{eq:dTilda} that is used in DES97 switches between the RANS length scale and LES length scale so that the model behaves in RANS-like and LES-like manners as illustrated in Figure 1. The length scale \eqref{eq:dTilda} depends only on the grid used in the simulation and it is solution-independent.

Menter’s SST Based DES Formulation

• k-$$\omega$$ SST DES model
• k-$$\omega$$ SST DDES model
• k-$$\omega$$ SST IDDES model

There is a project called DESider (Detached Eddy Simulation for Industrial Aerodynamics) and many other DES models have been developed based on different RANS models, including RSM ones.

Other Models

• $$v^2-f$$ DES model

We can check the RANS and LES regions using the DESModelRegions function objects in OpenFOAM.

• k-$$\omega$$ SST SAS model
• FLOMANIA (Flow Physics Modelling An Integrated Approach) Project: 2002-2004
• DESider (Detached Eddy Simulation for Industrial Aerodynamics) Project: 2004-2007
• ATAAC (Advanced Turbulence Simulation for Aerodynamic Application Challenges) Project
##### References

[1] P. R. Spalart, W.-H. Jou, M. Strelets and S. R. Allmaras, Comments on the feasibility of LES for wings, and on a hybrid RANS/LES approach. 1st AFOSR Int. Conf. on DNS/LES, Aug. 4-8, 1997, Ruston, LA. In “Advances in DNS/LES”, C. Liu and Z. Liu Eds., Greyden Press, Columbus, OH.
[2] M. Strelets, Detached Eddy Simulation of Massively Separated Flows. AIAA, 2001-0879.
[3] P. R. Spalart, S. Deck, M. L. Shur, K. D. Squires, M. Kh. Strelets, and A. Travin, A new version of detached-eddy simulation, resistant to ambiguous grid densities. Theoretical and Computational Fluid Dynamics, 20(3), 181-195, 2006.
[4] M. L. Shur, P. R. Spalart, M. Kh. Strelets, and A. Travin, A hybrid RANS-LES approach with delayed-DES and wall-modelled LES capabilities. International Journal of Heat and Fluid Flow, 29, 1638–1648, 2008.
[5] Detached eddy simulation (DES). Available at: http://www.cfd-online.com/Wiki/Detached_eddy_simulation_(DES) [Accessed: 03 May 2019].
[6] OpenFOAM® v3.0+: New Solver and Physical Modelling Functionality. Available at: https://www.openfoam.com/releases/openfoam-v3.0+/solvers-and-physics.php#physics-kOmegaSSTDES [Accessed: 03 May 2019].
[7] C. M. Winkler, A. J. Dorgan and M. Mani, Scale Adaptive Simulations of Turbulent Flows on Unstructured Grids. AIAA, 2011-3559.

## Calculation of Reynolds stress field in OpenFOAM

Keywords
Reynolds stress, shear flow, fieldAverage, variance, covariance

In the simulations of wall-bounded turbulent flows, such as turbulent channel flow [1], the distributions of the Reynolds stress components in the wall-normal direction are matters of interest to researchers and engineers.

For a statistically steady turbulent flow, a prime-squared mean of the velocity field (UPrime2Mean in OpenFOAM) gives the Reynolds stress field

\begin{eqnarray}
\frac{1}{N} \sum_{k=1}^{N} \left( U_i^{(k)} – \overline{U_i} \right)^2
&=& \frac{1}{N} \sum_{k=1}^{N} {u_{i}^{(k)}}^2 \\
&=& \frac{1}{N} \sum_{k=1}^{N}
\left[
\begin{array}{ccc}
u_{1}^{(k)}u_{1}^{(k)} & u_{1}^{(k)}u_{2}^{(k)} & u_{1}^{(k)}u_{3}^{(k)} \\
u_{2}^{(k)}u_{1}^{(k)} & u_{2}^{(k)}u_{2}^{(k)} & u_{2}^{(k)}u_{3}^{(k)} \\
u_{3}^{(k)}u_{1}^{(k)} & u_{3}^{(k)}u_{2}^{(k)} & u_{3}^{(k)}u_{3}^{(k)}
\end{array}
\right] \\
&=&
\left[
\begin{array}{ccc}
\overline{u_{1}u_{1}} & \overline{u_{1}u_{2}} & \overline{u_{1}u_{3}} \\
\overline{u_{2}u_{1}} & \overline{u_{2}u_{2}} & \overline{u_{2}u_{3}} \tag{1} \\
\overline{u_{3}u_{1}} & \overline{u_{3}u_{2}} & \overline{u_{3}u_{3}}
\end{array}
\right], \\
\end{eqnarray}

where $$U_i^{(k)}$$ denotes the instantaneous velocity field $$U_i(\boldsymbol{x}, k\Delta t)$$ and it is decomposed into the time-averaged mean $$\overline{U_i}$$ and fluctuating components $${u_{i}}^{(k)}$$,

{U_{i}}^{(k)} = \overline{U_i} + {u_{i}}^{(k)}. \tag{2}

The diagonal components of the Reynolds stress tensor are the variances of the velocity components and the off-diagonal elements are the covariances of them.

 Function Object – fieldAverage

 Turbulent Channel Flow

Turbulent channel flow is one of the most fundamental wall-bounded shear flows and it has been widely used to study the structure of near-wall turbulence. Many DNS calculations have been carried out and produced a lot of informative data, which has contributed considerably to the improvements of other turbulence models, such as RANS and LES.

Moser et al. have released the statistical data obtained from their DNS of turbulent channel flow on their web site [2]. The data set contains, among other information, the components of the Reynolds stress tensor at three Reynolds numbers

Re_{\tau} = \frac{u_{\tau}h}{\nu} \approx 180,\; 395,\; 590, \tag{3}

where $$u_{\tau}$$ is the friction velocity, $$h$$ is the channel half-height and $$\nu$$ is the kinematic viscosity.

I will post my DNS results, which will help us to understand the flow behaviors in the near-wall region. The channel395 tutorial in OpenFOAM is the case of large eddy simulation (LES) instead of DNS at $$Re_{\tau}=395$$.

 References

[1] R. D. Moser, J. Kim and N. N. Mansour, Direct numerical simulation of turbulent channel flow up to $${\rm Re}_{\tau}$$=590. Physics of Fluids 11(4), 943-945, 1999.
[2] DNS Data for Turbulent Channel Flow. Available at: http://turbulence.ices.utexas.edu/data/MKM/ [Accessed: 9 April 2017].
[3] turbulence.ices.utexas.edu file server. Available at: http://turbulence.ices.utexas.edu/ [Accessed: 4 May 2017].

## Filter Width – maxDeltaxyz

Several options are available in OpenFOAM for calculation of the filter width used in the large eddy simulation (LES) and detached eddy simulation (DES). In this blog post, the maxDeltaxyz option is covered in some detail.

OpenFOAM Version: OpenFOAM-dev, OpenFOAM v1612+

 Implementation in OpenFOAM

The maxDeltaxyz option calculates the filter width of the $$i-$$th cell $$\Delta_i$$ by taking the maximum distance between the cell center $$P_i$$ and each face center $$F_j$$

\Delta_i = {\rm deltaCoeff} \times \max_{1 \le j \le n_i} \left\{ \overline{P_iF_j} \right\}, \tag{1} \label{eq:deltaxyz}

where $${\rm deltaCoeff}$$ is a constant of proportion (user input) and $$n_i$$ is the number of the faces of the $$i-$$th cell.

For a regular hexahedral cell shown in Fig. 1, the computed $$\Delta$$ using Eq. \eqref{eq:deltaxyz} equals one-half of the maximum cell width $$\Delta x$$, so the deltaCoeff coefficient should be set to 2 in the turbulenceProperties file as shown below.

The filter width is calculated in the following function.

 Implementation in OpenFOAM v1612+

In OpenFOAM v1612+, the face normal vectors are considered in order to take into account the mesh non-orthogonality.

This definition \eqref{eq:deltaxyz} is often used in the detached eddy simulation (DES) where some anisotropic grid cells such as “book”, “pencil” and “ribbon”-shaped cells exist in a boundary layer mesh. The position where switching between the RANS and LES modes occurs in the DES97 model [1] depends on how to calculate the filter width $$\Delta$$

\tilde{d} \equiv {\rm min}\left( d, C_{DES}\Delta \right), \tag{2} \label{eq:dTilda}

where $$d$$ is the distance to the closest wall and $$C_{DES}$$ is a calibration constant.
Several modified length scales have been developed to prevent a delay of the Kelvin-Helmholtz instability in free and separated shear layers [2].

 References

[1] P. R. Spalart, W.-H. Jou, M. Strelets and S. R. Allmaras, Comments on the feasibility of LES for wings, and on a hybrid RANS/LES approach. 1st AFOSR Int. Conf. on DNS/LES, Aug. 4-8, 1997, Ruston, LA. In “Advances in DNS/LES”, C. Liu and Z. Liu Eds., Greyden Press, Columbus, OH.
[2] P. R. Spalart, Detached-Eddy Simulation. Annu. Rev. Fluid Mech. 41, 181-202, 2009.

## WALE SGS model in OpenFOAM

Last Updated: May 6, 2019

The WALE (Wall-Adapting Local Eddy-viscosity) SGS model is one of the major SGS models. It is an algebraic eddy viscosity model (0-equation model) as with the Smagorinsky SGS model, but it has some excellent features that the Smagorinsky model does not have.

##### Implementation in OpenFOAM

In the WALE model, the subgrid scale eddy viscosity $$\nu_{sgs}$$ is evaluated as

\nu_{sgs} = C_{k} \Delta \sqrt{k_{sgs}}, \tag{1} \label{eq:nusgs}

where $$C_{k}$$ is a model constant and $$k_{sgs}$$ is the subgrid scale kinetic energy. You can find its definition in this post.

The traceless symmetric part of the square of the velocity gradient tensor $$S^d$$ is calculated in the following function.
\begin{eqnarray}
S_{ij}^d = \frac{1}{2} \left( \frac{\partial \overline{u}_k}{\partial x_i}\frac{\partial \overline{u}_j}{\partial x_k} + \frac{\partial \overline{u}_k}{\partial x_j}\frac{\partial \overline{u}_i}{\partial x_k} \right) – \frac{1}{3} \delta_{ij} \frac{\partial \overline{u}_k}{\partial x_l}\frac{\partial \overline{u}_l}{\partial x_k}, \tag{2} \label{eq:sd}
\end{eqnarray}
where $$\delta_{ij}$$ is the Kronecker delta.

These tensor operations used above are summarized in the following post:

The subgrid scale kinetic energy $$k_{sgs}$$ is

k_{sgs} = \left( \frac{C_w^2 \Delta}{C_k} \right)^2 \frac{\left( S_{ij}^d S_{ij}^d \right)^3}{\left( \left( \overline{S}_{ij} \overline{S}_{ij} \right)^{5/2} + \left( S_{ij}^d S_{ij}^d \right)^{5/4} \right)^2}, \tag{3} \label{eq:ksgs}

where

\overline{S}_{ij} = \frac{1}{2} \left( \frac{\partial \overline{u}_{i}}{\partial x_{j}} + \frac{\partial \overline{u}_{j}}{\partial x_{i}}\right), \tag{4} \label{eq:s}

is the resolved-scale strain rate tensor.

Finally, we can get the following expression by substituting Eq. \eqref{eq:ksgs} into Eq. \eqref{eq:nusgs}:
\begin{eqnarray}
\nu_{sgs} = \left( C_w \Delta \right)^2 \frac{\left( S_{ij}^d S_{ij}^d \right)^{3/2}}{\left( \overline{S}_{ij} \overline{S}_{ij} \right)^{5/2} + \left( S_{ij}^d S_{ij}^d \right)^{5/4}}, \tag{5} \label{eq:nusgs2}
\end{eqnarray}
It is the same as Eq. (13) in [1].

##### Features
• Algebraic eddy viscosity model (0-equation model)
• The rotation rate is taken into account in the calculation of $$\nu_{sgs}$$
• To be able to handle transition
• Damping is Not necessary for $$\nu_{sgs}$$ in the near-wall region
##### References

Final expression of the SGS eddy viscosity is relatively simple and the implementation into OpenFOAM is not so complicated but the process of deriving the expression described in [1] is not easy to understand. I want to develop my SGS model someday 🙂

## One equation eddy-viscosity SGS model in OpenFOAM

OpenFOAM Version: OpenFOAM-dev, OpenFOAM-1606+

 Implementation in OpenFOAM

Assumption 1: Eddy viscosity approximation
As in the case of the Smagorinsky SGS model, the one equation eddy viscosity SGS model uses the eddy viscosity approximation as the name suggests, so the subgrid scale stress tensor $$\tau_{ij}$$ is approximated as follows:

\tau_{ij} \approx \frac{2}{3} k_{sgs} \delta_{ij} – 2 \nu_{sgs} dev(\overline{D})_{ij}, \tag{1} \label{eq:tauij}

where $$\nu_{sgs}$$ is the subgrid scale eddy viscosity, $$\overline{D}$$ is the resolved-scale strain rate tensor defined as

\overline{D}_{ij} = \frac{1}{2} \left( \frac{\partial \overline{u}_{i}}{\partial x_{j}} + \frac{\partial \overline{u}_{j}}{\partial x_{i}}\right), \tag{2} \label{eq:Dij}

and the subgrid scale kinetic energy $$k_{sgs}$$ is

k_{sgs} = \frac{1}{2} \tau_{kk} = \frac{1}{2} \left( \overline{u_{k}u_{k}} – \overline{u}_{k}\overline{u}_{k} \right). \tag{3} \label{eq:ksgs}

The subgrid scale eddy viscosity $$\nu_{sgs}$$ is computed using $$k_{sgs}$$

\nu_{sgs} = C_{k} \sqrt{k_{sgs}} \Delta, \tag{4} \label{eq:nusgs}

where $$C_{k}$$ is a model constant whose default value is $$0.094$$.

The procedure so far is the same as the Smagorinsky SGS model but there is a difference in what follows. These models are different in terms of how to compute $$k_{sgs}$$. The Smagorinsky model assumes the local equilibrium to compute $$k_{sgs}$$ but the one equation eddy viscosity model solves a transport equation of $$k_{sgs}$$ [1].

The second category of SGS model is one-equation eddy viscosity models. The main reason to develop the one-equation SGS models is to overcome the deficiency of local balance assumption between the SGS energy production and dissipation adopted in algebraic eddy viscosity models. Such a phenomenon may occur in high Reynolds number flows and/or in the cases of coarse grid resolution. The first one-equation eddy viscosity SGS model was theoretically derived by Yoshizawa and Horiuti [2], in which the SGS kinetic energy is defined as $$k_{sgs} = \frac{1}{2}(\overline{u_k^2} – \overline{u}_k^2)$$, and the SGS eddy viscosity, $$\nu_{sgs}$$, is computed using $$k_{sgs}$$ as $$\nu_{sgs} = C_{k} k_{sgs}^{1/2} \Delta$$. A transportation equation is derived to account for the historic effect of $$k_{sgs}$$ due to production, dissipation and diffusion:

Transport Equation of $$k_{sgs}$$
\begin{align}
&\frac{\partial (\rho k_{sgs})}{\partial t} + \frac{\partial (\rho \overline{u}_j k_{sgs})}{\partial x_{j}} – \frac{\partial}{\partial x_{j}} \left[ \rho \left(\nu + \nu_{sgs}\right) \frac{\partial k_{sgs}}{\partial x_{j}} \right] \\
&= \; – \rho \tau_{ij} : \overline{D}_{ij} \; – \; C_{\epsilon} \frac{\rho k_{sgs}^{3/2}}{\Delta}, \tag{5} \label{eq:ksgsEqn}
\end{align}
where the operator : is a double inner product, $$\rho$$ is the density, $$\nu$$ is the kinetic viscosity and $$C_{\epsilon}$$ is another model constant. The terms in Eq. \eqref{eq:ksgsEqn} are, starting from the left, the time derivative term, convective term, diffusion term, production term and dissipation term. In the case of the Smagorinsky SGS model, only the production and dissipation terms are taken into account with the assumption of the local equilibrium.

The production term in Eq. \eqref{eq:ksgsEqn} can be rearranged to yield the expression used in the source code as follows:
\begin{align}
– \rho \tau_{ij} : \overline{D}_{ij} &= \left( – \frac{2}{3} \rho k_{sgs} \delta_{ij} + 2 \rho \nu_{sgs} dev(\overline{D})_{ij} \right) : \overline{D}_{ij} \\
&=\;- \frac{2}{3} \rho k_{sgs} \frac{\partial \overline{u}_{k}}{\partial x_{k}} \\
&\;\;\;+ \rho \nu_{sgs} \frac{\partial \overline{u}_i}{\partial x_{j}} \left( 2\overline{D}_{ij} – \frac{1}{3} tr(2\overline{D}) \delta_{ij} \right). \tag{6} \label{eq:production}
\end{align}
Please see the appendix below for more detailed process of this deformation.

 References

[1] S. Huang and Q. S. Li, A new dynamic one-equation subgrid-scale model for large eddy simulations. INTERNATIONAL JOURNAL FOR NUMERICAL METHODS IN ENGINEERING, 81, 835-865, 2010.
[2] A. Yoshizawa and K. Horiuti, A Statistically-Derived Subgrid-Scale Kinetic Energy Model for the Large-Eddy Simulation of Turbulent Flows. Journal of the Physical Society of Japan, 54(8), 2834-2839, 1985.

 Appendix: Deformation of the Production Term

\begin{align}
&\; 2 \rho \nu_{sgs} dev(\overline{D})_{ij} : \overline{D}_{ij} \\
=&\; \rho \nu_{sgs} \frac{1}{2} \left( \frac{\partial \overline{u}_i}{\partial x_j} + \frac{\partial \overline{u}_j}{\partial x_i} – \frac{2}{3} tr(\overline{D}) \delta_{ij} \right) \left( \frac{\partial \overline{u}_i}{\partial x_j} + \frac{\partial \overline{u}_j}{\partial x_i} \right) \\
=&\; \rho \nu_{sgs} \frac{1}{2} \left( 2 \frac{\partial \overline{u}_i}{\partial x_j}\frac{\partial \overline{u}_i}{\partial x_j} + 2 \frac{\partial \overline{u}_i}{\partial x_j}\frac{\partial \overline{u}_j}{\partial x_i} – \frac{4}{3} tr(\overline{D}) \delta_{ij} \frac{\partial \overline{u}_i}{\partial x_j} \right) \\
=&\; \rho \nu_{sgs} \frac{\partial \overline{u}_i}{\partial x_j} \left( \frac{\partial \overline{u}_i}{\partial x_j} + \frac{\partial \overline{u}_j}{\partial x_i} – \frac{2}{3} tr(\overline{D}) \delta_{ij} \right) \\
=&\; \rho \nu_{sgs} \frac{\partial \overline{u}_i}{\partial x_j} \left( 2 \overline{D}_{ij} \; – \frac{1}{3} tr(\overline{2D}) \delta_{ij} \right) \\
=&\; \rho G \tag{7} \label{eq:productionTerm}
\end{align}

## Smagorinsky SGS model in OpenFOAM

Last Updated: May 4, 2019

The Smagorinsky subgrid scale (SGS) model was developed by Joseph Smagorinsky in the meteorological community in the 1960s. It is based on the eddy viscosity assumption, which postulates a linear relationship between the SGS shear stress and the resolved rate of strain tensor. This model serves as a base for other SGS models.

OpenFOAM Version: OpenFOAM-dev, OpenFOAM-1606+

##### Implementation in OpenFOAM

The Smagorinsky SGS model is the oldest and best known subgrid scale model. The subgrid scale stress tensor $$\tau_{ij}$$ is
\begin{align}
\tau_{ij} &= \overline{u_{i}u_{j}} – \overline{u}_{i}\overline{u}_{j} \tag{1a}\\
&= \frac{1}{3} \tau_{kk} \delta_{ij} + \left( \tau_{ij} \; – \frac{1}{3} \tau_{kk} \delta_{ij} \right) \tag{1b} \label{eq:smg2}\\
&\approx \frac{1}{3} \tau_{kk} \delta_{ij} – 2 \nu_{sgs} dev(\overline{D})_{ij} \tag{1c} \label{eq:smg3}\\
&= \frac{2}{3} k_{sgs} \delta_{ij} – 2 \nu_{sgs} dev(\overline{D})_{ij}, \tag{1d} \label{eq:smg4}
\end{align}
where $$\nu_{sgs}$$ is the subgrid scale eddy viscosity and the resolved-scale strain rate tensor $$\overline{D}_{ij}$$ is defined as

\overline{D}_{ij} = \frac{1}{2} \left( \frac{\partial \overline{u}_{i}}{\partial x_{j}} + \frac{\partial \overline{u}_{j}}{\partial x_{i}}\right), \tag{2} \label{eq:D}

and the subgrid scale kinetic energy $$k_{sgs}$$ is

k_{sgs} = \frac{1}{2} \tau_{kk} = \frac{1}{2} \left( \overline{u_{k}u_{k}} – \overline{u}_{k}\overline{u}_{k} \right). \tag{3} \label{eq:ksgs}

The subgrid scale stress tensor $$\tau_{ij}$$ is split into an isotropic part $$\frac{1}{3} \tau_{kk} \delta_{ij}$$ and an anisotropic part $$\tau_{ij} \; – \frac{1}{3} \tau_{kk} \delta_{ij}$$ in Eq. \eqref{eq:smg2}. I think $$dev(\overline{D})$$ is used instead of $$\overline{D}$$ because the anisotropic part is a traceless tensor in Eq. \eqref{eq:smg3}.

Assumption 1: Eddy viscosity approximation
In analogy with the molecular viscous stress in laminar flows, the anisotropic part is approximated by relating it to the resolved rate of strain tensor $$\overline{D}_{ij}$$ as in Eq. \eqref{eq:smg3}

\tau_{ij} \; – \frac{1}{3} \tau_{kk} \delta_{ij} \approx \; – 2 \nu_{sgs} dev(\overline{D})_{ij}. \tag{1c}

In OpenFOAM, the subgrid scale viscosity is computed as

\nu_{sgs} = C_{k} \Delta \sqrt{k_{sgs}}, \tag{4} \label{eq:nusgs}

where $$C_{k}$$ is a model constant whose default value is $$0.094$$ and $$\Delta$$ is the grid size that defines the subgrid length scale. The method for calculating $$\Delta$$ is specified in the turbulenceProperties file and the available options are as follows:

• cubeRootVol
• maxDeltaxyz
• maxDeltaxyzCubeRoot
• smooth
• vanDriest
• Prandtl
• IDDESDelta

For the Smagorinsky SGS model, the vanDriest option is the first choice.

Our remaining task is to evaluate the distribution of subgrid scale kinetic energy $$k_{sgs}$$.

Assumption 2: Local equilibrium
The SGS kinetic energy $$k_{sgs}$$ is computed with the assumption of the balance between the subgrid scale energy production and dissipation (local equilibrium)

\overline{D} : \tau_{ij} + C_{\epsilon} \frac{k_{sgs}^{1.5}}{\Delta} = 0, \tag{5} \label{eq:equilibrium}

where the operator : is a double inner product of two second-rank tensors that can be evaluated as the sum of the 9 products of the tensor components. We can compute $$k_{sgs}$$ by solving Eq. \eqref{eq:equilibrium} as shown below:
\begin{align}
&\overline{D} : \left( \frac{2}{3} k_{sgs} I -2 \nu_{sgs} dev(\overline{D}) \right) + C_{\epsilon} \frac{k_{sgs}^{1.5}}{\Delta} = 0 \\
\Leftrightarrow &\quad \overline{D} : \left( \frac{2}{3} k_{sgs} I -2 C_{k} \Delta \sqrt{k_{sgs}} dev(\overline{D}) \right) + C_{\epsilon} \frac{k_{sgs}^{1.5}}{\Delta} = 0 \\
\Leftrightarrow &\quad \sqrt{k_{sgs}} \left( \frac{C_{\epsilon}}{\Delta} k_{sgs} + \frac{2}{3} tr(\overline{D}) \sqrt{k_{sgs}} -2 C_{k} \Delta \left( dev(\overline{D}) : \overline{D} \right) \right) = 0 \\
\Leftrightarrow &\quad a k_{sgs} + b \sqrt{k_{sgs}} – c = 0 \\
\Rightarrow &\quad k_{sgs} = \left( \frac{-b + \sqrt{b^2 + 4ac}}{2a} \right)^2, \tag{6} \label{eq:ksgs2}
\end{align}
where
\begin{eqnarray}
\left\{
\begin{array}{l}
a = \frac{C_{\epsilon}}{\Delta} \\
b = \frac{2}{3} tr(\overline{D}) \tag{7} \label{eq:coeffs}\\
c = 2 C_{k} \Delta \left( dev(\overline{D}) : \overline{D} \right).
\end{array}
\right.
\end{eqnarray}

In the case of incompressible flows, Eq. \eqref{eq:coeffs} reduces to
\begin{eqnarray}
\left\{
\begin{array}{l}
b = \frac{2}{3} tr(\overline{D}) = 0 \\
c = 2 C_{k} \Delta \left( dev(\overline{D}) : \overline{D} \right) = C_{k} \Delta \vert \overline{D} \vert^2, \tag{8} \label{eq:coeffsInc}
\end{array}
\right.
\end{eqnarray}
where

\vert \overline{D} \vert = \sqrt{2 \overline{D} : \overline{D}}. \tag{9} \label{eq:abs}

By substituting these relations into Eq. \eqref{eq:ksgs2}, we can get

k_{sgs} = \frac{c}{a} = \frac{C_{k} \Delta^2 \vert \overline{D} \vert^2}{C_{\epsilon}}. \tag{10} \label{eq:ksgsInc}

We can get the following expression for the SGS eddy viscosity in the case of incompressible flows by substituting the above equation into Eq. \eqref{eq:nusgs}

\nu_{sgs} = C_{k} \sqrt{\frac{C_{k}}{C_{\epsilon}}} \Delta^2 \vert \overline{D} \vert. \tag{11} \label{eq:nusgsInc}

By comparing it with the following common expression in the literature

\nu_{sgs} = \left( C_{s} \Delta \right)^2 \vert \overline{D} \vert, \tag{12} \label{eq:nusgstext}

we can get the following relation for the Smagorinsky constant $$C_{s}$$:

C_{s}^2 = C_{k} \sqrt{\frac{C_{k}}{C_{\epsilon}}}. \tag{13} \label{eq:smgConstInc}

##### Features
• Algebraic eddy viscosity model (0-equation model)
• The rotation rate is not taken into account in the calculation of $$\nu_{sgs}$$
• Based on the local equilibrium assumption
• Appropriate value of the Smagorinsky constant is problem dependent
• Not to be able to handle transition
• Not to be able to deal with energy backscatter
• Van-Driest damping is necessary for the SGS eddy viscosity in the near-wall region

The last three items arise from the fact that $$\nu_{sgs} \ge 0$$ in Eq. \eqref{eq:nusgs}.

## Large Eddy Simulation (LES)

Large Eddy Simulation (LES) is one of the most promising methods for computing industry-relevant turbulent flows. It is used to predict unsteady flow behaviors with lower computational cost as compared to Direct Numerical Simulation (DNS).

The essential idea of LES dates back to Joseph Smagorinsky (1963) in meteorology and to James W. Deardorff (1970) in engineering.

Keywords
Inertial subrange, Wall-modeled LES, Wall-resolved LES, Smagorinsky, Taylor microscale

##### Nomenclature – Abbreviations
• LES: Large Eddy Simulation
• ILES: Implicit Large Eddy Simulation
• MILES: Monotone Integrated Large Eddy Simulation
• VLES: Very Large Eddy Simulation
• WMLES: Wall-Modeled Large Eddy Simulation
• SGS: subgrid-scale
##### What is Large Eddy Simulation?

Professor Parviz Moin of Stanford University briefly explains what LES is in his interview. I’m surprised that he received his Masters and Ph.D. degrees in Mathematics and Mechanical Engineering from Stanford in the same year (1978)!

The verification of LES simulations is difficult because both the error induced by the SGS model and the numerical discretization error are dependent on the grid resolution. In the following presentation (37:40-42:06), he introduces one method using explicit filtering to distinguish the errors [6]:

If the grid-independent solution of the explicit filtered LES equations fails to accurately predict the filtered DNS flow field, its failure can be solely attributed to the capability of the subgrid stress model employed.

##### Categorization of LES

LES can be classified in terms of the following points as shown in Fig. 1:

• How the effect of the SGS stress is modeled
• How the spatial filtering is applied in the solution procedure.

In OpenFOAM, LES with the implicit filtering is implemented, in which only the filter width is specified and the filter shape is not. There exists this ambiguity in the definition of the filter. Lund [7] provides a clear description of the implicit filtering in LES:

The nearly universal approach is to simply write down the filtered Navier-Stokes equations together with an assumed model for the subgrid-scale stresses and then apply the desired spatial discretization to this “filtered” system. Although it is rarely mentioned, what one is doing by adopting this procedure is to imagine that the finite support of the computational mesh together with the low-pass characteristics of the discrete differentiating operators act as an effective filter. One then directly associates the computed velocity field with the filtered velocity. This procedure will be referred to as implicit filtering since an explicit filtering operation never appears in the solution procedure.

The explicit filtering methodology aims to damp the error-prone length scales smaller than the filter width that can be generated by the nonlinear convective term. In the case of the implicit filtering, we hope that this error has little effect on the resolved scales.

##### SGS models in OpenFOAM
• Smagorinsky

Smagorinsky SGS model

• kEqn

One equation eddy-viscosity model

• dynamicLagrangian

Dynamic SGS model with Lagrangian averaging

• dynamicKEqn

Dynamic one equation eddy-viscosity model

• WALE

Wall-adapting local eddy-viscosity (WALE) SGS model

• DeardorffDiffStress

Differential SGS Stress Equation Model

For the dynamic SGS models, the spatial averaging operations of the coefficients are often performed to stabilize the calculation. The homogeneousDynSmagorinsky model that had been implemented in older versions takes the average of the coefficient in the whole computational domain.

##### Other SGS models
• Vreman SGS model [3]

\nu_{sgs} = c \sqrt{\frac{B_{\beta}}{\alpha_{ij}\alpha_{ij}}}

##### Calculation of filter width in OpenFOAM

The method for calculating the filter width $$\Delta$$ is specified in the turbulenceProperties file. The available options in OpenFOAM are as follows:

• cubeRootVol
• maxDeltaxyz
• maxDeltaxyzCubeRoot
• smooth
• vanDriest
• Prandtl
• IDDESDelta
##### Grid resolution measures

The choice of the computational grid size has sensible impact on the quality of the LES simulations and there are several techniques to assess the grid resolution in the LES computations [4].

• Estimations based on prior RANS results
• Single-grid estimators
• Multi-grid estimators
##### Principal use
• Aeroacoustics
• Aerodynamics
• Combustion

## Direct Numerical Simulation (DNS)

Keywords
Inertial subrange, Kolmogorov microscale, Taylor microscale, Görtler vortex

 What is Direct Numerical Simulation?
 Channel Flow
 Taylor-Couette Flow
 Taylor-Green Vortex
 Backward-Facing Step
 DNS using iconCFD
 List of Laboratories (Japan)
 References

## Reynolds-Averaged Navier-Stokes (RANS)

Last Updated: May 11, 2019

Keywords
Boussinesq approximation, closure problem, Reynolds averaging

##### Reynolds Average

It is computationally expensive to resolve the wide range of time and length scales observed in turbulent flows. We now consider decomposing a flow property $$f$$, such as velocity and pressure, into a mean component $$\overline{f}$$ and a fluctuating component $$f’$$.

f(\boldsymbol{x}, t) = \overline{f}(\boldsymbol{x}, t) + f'(\boldsymbol{x}, t), \tag{1} \label{eq:decomposition}

where $$\boldsymbol{x}$$ is the position vector and $$t$$ is time.

The Reynolds-averaged Navier-Stokes (RANS) turbulence models aim to solve the mean flow $$\overline{f}$$ that changes more slowly in time and space than the original variable $$f$$. The governing equations of the mean component will be derived later.

There are many averaging operations defined in mathematics but the RANS models use the Reynolds average. It is briefly described in the newly published textbook by Kajishima and Taira [1].

For the discussion in this chapter, let us redefine the averaging operation such that it satisfies

\overline{f’} = 0,\;\;\overline{f’ \overline{f}} = 0,\;\;\overline{\overline{f}} = \overline{f}. \tag{7.2} \label{eq:reave}

These relations in Eq. \eqref{eq:reave} are referred to as the Reynolds-averaging laws. The ensemble average that satisfies these laws is called the Reynolds average. This conceptual averaging operation conveniently removes fluctuating components from the flow field variables without explicitly defining the spatial length scale used in the averaging operation.

The ensemble average that appears in the above definition is defined as (and usually denoted as $$\langle f \rangle$$)

\langle f \rangle(\boldsymbol{x}, t) \equiv \lim_{N \to \infty}\frac{1}{N}\sum_{i=1}^{N}f_{i}(\boldsymbol{x}, t), \tag{2} \label{eq:ensembleAve}

where $$f_i$$ are the samples of $$f$$ and $$N$$ is the number of samples. In other words, it is the average of the instantaneous values of the property at a given point in space $$\boldsymbol{x}$$ and time $$t$$ over a large number of repeated identical experiments. In general, this ensemble average varies with space and time (time-dependent).

For the stationary random processes, we can define the time average $$f_T$$:

f_{T}(\boldsymbol{x}) \equiv \frac{1}{T} \int_{0}^{T} f(\boldsymbol{x}, t) dt, \tag{3} \label{eq:timeAve}

where $$T$$ is the integration time. In the case of stationary random processes, the time averages equal the ensemble averages as stated in [3]:

if the signal is stationary, the time average defined by equation \eqref{eq:timeAve} is an unbiased estimator of the true average $$\langle f \rangle$$. Moreover, the estimator converges to $$\langle f \rangle$$ as the time becomes infinite; i.e., for stationary random processes

\langle f \rangle(\boldsymbol{x}) = \lim_{T \to \infty} \frac{1}{T} \int_{0}^{T} f(\boldsymbol{x}, t) dt, \tag{8.28} \label{eq:aveRelation}

Thus the time and ensemble averages are equivalent in the limit as $$T \to \infty$$, but only for a stationary random process.

Interested readers might want to search by the keyword ”ergodic hypothesis” on the relation between the ensemble and time averages.

To Be Updated

##### Closure Problem – Reynolds Stress

The linear eddy viscosity models (LEVM) assume the linear stress-strain relationship and employ the eddy-viscosity concept (Boussinesq approximation) introduced by Joseph Valentin Boussinesq

-\rho\overline{u_i u_j} = \mu_t \left(\frac{\partial \overline{U}_i}{\partial x_j} + \frac{\partial \overline{U}_j}{\partial x_i} \right) -\frac{2}{3}\delta_{ij}\rho k. \tag{4} \label{eq:BoussinesqApprox}

##### RANS Models in OpenFOAM

Linear Eddy Viscosity Model (LEVM)

Nonlinear Eddy Viscosity Model (NLEVM)

Reynolds Stress Model (RSM)

##### Differential Reynolds Stress model
• SSG/LRR-$$\omega$$
• JH-$$\omega^h$$
##### References

[1] T. Kajishima and K. Taira, Computational Fluid Dynamics: Incompressible Turbulent Flows. Springer, 2016.
[2] H. K. Versteeg and W. Malalasekera, An introduction to Computational Fluid Dynamics: The Finite Volume Method. Person Prentice Hall, 1995.
[3] W. K. George, Lectures in Turbulence for the 21st Century.

## Spalart-Allmaras DDES model in OpenFOAM

The formulation of Spalart-Allmaras DDES (Delayed Detached Eddy Simulation) model is described in [1].

Keywords:
delayed DES, shielding function, modeled-stress depletion

 Modifications Compared to Original DES97 Model

The DES length scale $$\tilde{d}$$ is modified compared to that of DES97 model to prevent a premature switch of DES to LES mode within boundary layers. This incursion can be a cause of modeled-stress depletion (MSD) issue that the DES97 formulation suffers from.

Two functions $$r_d$$ and $$f_d$$ are introduced in new definition \eqref{eq:dTilda} and it now depends not only on the grid but also the time-dependent eddy-viscosity field as Eq. \eqref{eq:rd} shows. The new length scale \eqref{eq:dTilda} detects the boundary layer and delays the switch to LES mode until outside of it. For this reason, the method was named delayed DES.

 Implementation in OpenFOAM

r_d \equiv \frac{\nu_t + \nu}{\sqrt{U_{i,j}U_{i,j}}\kappa^2 d^2} \tag{1} \label{eq:rd}

where $$\nu_{t}$$ is the kinematic eddy viscosity, $$\nu$$ the molecular viscosity, $$U_{i,j}$$ the velocity gradients, $$\kappa$$ the Kármán constant and $$d$$ the distance to the wall.

The quantity $$r_{d}$$ is used in the shielding function $$f_{d}$$:

f_d \equiv 1 – {\rm tanh} \left[ \left( 8r_d \right)^3 \right] \tag{2} \label{eq:fd}

which is designed to be 1 in the LES region, where $$r_{d}\ \ll 1$$, and 0 elsewhere.

\tilde{d} \equiv d \; – f_d {\rm max}(0, d \; – \; C_{DES}\Delta) \tag{3} \label{eq:dTilda}

 References

[1] P. R. Spalart, S. Deck, M. L. Shur, K. D. Squires, M. Kh. Strelets and A. Travin, A new version of detached-eddy simulation, resistant to ambiguous grid densities. Theoretical and Computational Fluid Dynamics 20(3), 181-195, 2006.