**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:

\begin{equation}

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

\end{equation}

where \(\nu_{sgs}\) is the subgrid scale eddy viscosity, \(\overline{D}\) is the resolved-scale strain rate tensor 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:Dij}

\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 eddy viscosity \(\nu_{sgs}\) is computed using \(k_{sgs}\)

\begin{equation}

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

\end{equation}

where \(C_{k}\) is a model constant whose default value is \(0.094\).

1 2 3 4 5 6 7 8 9 |
template<class BasicTurbulenceModel> void kEqn<BasicTurbulenceModel>::correctNut() { this->nut_ = Ck_*sqrt(k_)*this->delta(); this->nut_.correctBoundaryConditions(); fv::options::New(this->mesh_).correct(this->nut_); BasicTurbulenceModel::correctNut(); } |

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.

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 |
template<class BasicTurbulenceModel> void kEqn<BasicTurbulenceModel>::correct() { if (!this->turbulence_) { return; } // Local references const alphaField& alpha = this->alpha_; const rhoField& rho = this->rho_; const surfaceScalarField& alphaRhoPhi = this->alphaRhoPhi_; const volVectorField& U = this->U_; volScalarField& nut = this->nut_; fv::options& fvOptions(fv::options::New(this->mesh_)); LESeddyViscosity<BasicTurbulenceModel>::correct(); volScalarField divU(fvc::div(fvc::absolute(this->phi(), U))); tmp<volTensorField> tgradU(fvc::grad(U)); volScalarField G(this->GName(), nut*(tgradU() && dev(twoSymm(tgradU())))); tgradU.clear(); tmp<fvScalarMatrix> kEqn ( fvm::ddt(alpha, rho, k_) + fvm::div(alphaRhoPhi, k_) - fvm::laplacian(alpha*rho*DkEff(), k_) == alpha*rho*G - fvm::SuSp((2.0/3.0)*alpha*rho*divU, k_) - fvm::Sp(this->Ce_*alpha*rho*sqrt(k_)/this->delta(), k_) + kSource() + fvOptions(alpha, rho, k_) ); kEqn.ref().relax(); fvOptions.constrain(kEqn.ref()); solve(kEqn); fvOptions.correct(k_); bound(k_, this->kMin_); correctNut(); } |

1 2 3 4 5 6 7 8 |
//- Return the effective diffusivity for k tmp<volScalarField> DkEff() const { return tmp<volScalarField> ( new volScalarField("DkEff", this->nut_ + this->nu()) ); } |

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}

Thanks, very presented

Hi Sharad,

Thank you for your compliment!

Why the source code contains two terms instead of one for the production. In other words why do we need to separate it into two? (Refer your equation 6.)

Hi Khedar,

Thanks for your question.

I think the reason why the production term is divided into two part is that the SGS kinetic energy \(k_{sgs}\) is only in the 1st term of the R.H.S. of Eq. (6). This term can take both positive and negative values so the SuSp operator is used in order to create diagonally dominant matrix.

Oh I understand now. Thanks for your response. I have one request. I found one of your presentations on cyclic, cyclicAMI and cylicACMI boundary conditions which was really helpful even though it was not totally in English. Do you have an English translation for that.