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
\begin{equation}
\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}
\end{equation}
and the subgrid scale kinetic energy \(k_{sgs}\) is
\begin{equation}
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}
\end{equation}
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}
\begin{equation}
\tau_{ij} \; – \frac{1}{3} \tau_{kk} \delta_{ij} \approx \; – 2 \nu_{sgs} dev(\overline{D})_{ij}. \tag{1c}
\end{equation}
In OpenFOAM, the subgrid scale viscosity is computed as
\begin{equation}
\nu_{sgs} = C_{k} \Delta \sqrt{k_{sgs}}, \tag{4} \label{eq:nusgs}
\end{equation}
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)
\begin{equation}
\overline{D} : \tau_{ij} + C_{\epsilon} \frac{k_{sgs}^{1.5}}{\Delta} = 0, \tag{5} \label{eq:equilibrium}
\end{equation}
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
\begin{equation}
\vert \overline{D} \vert = \sqrt{2 \overline{D} : \overline{D}}. \tag{9} \label{eq:abs}
\end{equation}
By substituting these relations into Eq. \eqref{eq:ksgs2}, we can get
\begin{equation}
k_{sgs} = \frac{c}{a} = \frac{C_{k} \Delta^2 \vert \overline{D} \vert^2}{C_{\epsilon}}. \tag{10} \label{eq:ksgsInc}
\end{equation}
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}
\begin{equation}
\nu_{sgs} = C_{k} \sqrt{\frac{C_{k}}{C_{\epsilon}}} \Delta^2 \vert \overline{D} \vert. \tag{11} \label{eq:nusgsInc}
\end{equation}
By comparing it with the following common expression in the literature
\begin{equation}
\nu_{sgs} = \left( C_{s} \Delta \right)^2 \vert \overline{D} \vert, \tag{12} \label{eq:nusgstext}
\end{equation}
we can get the following relation for the Smagorinsky constant \(C_{s}\):
\begin{equation}
C_{s}^2 = C_{k} \sqrt{\frac{C_{k}}{C_{\epsilon}}}. \tag{13} \label{eq:smgConstInc}
\end{equation}

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

References

[1] (Free Access) J. Smagorinsky, General circulation experiments with the primitive equations: I. The basic experiment*. Monthly weather review, 91(3), 99-164, 1963.