# 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}.

## Author: fumiya

CFD engineer in Japan

## 3 thoughts on “Smagorinsky SGS model in OpenFOAM”

1. MUHAMMED Manaf says:

Dear sir

If we want to change Cs . Then how can i change Ck and Ce, which one is more likely to change if you want to change Cs

2. Hi,

Is trace part of the SGS stress is added separately to the pressure at each time step or directly modified pressure is only solved.

1. fumiya says:

Hi NAVEEN,

Thanks for your comment.
The modified pressure is directly solved in OpenFOAM.

This site uses Akismet to reduce spam. Learn how your comment data is processed.