Last Updated: June 23, 2019
DES (Detached Eddy Simulation) turbulence model is pioneered in 1997 by Spalart et al. based on the Spalart-Allmaras RANS model [1]. It is referred to as DES97 model in the CFD (Computational Fluid Dynamics) community and several improved versions (DDES and IDDES) have been developed so far that have continually widened the areas of application. They have been implemented in both commercial and open-source CFD softwares and used on a regular basis in both industry and academia.
In this blog post, I will try to describe in detail DES97 model that serves as a base for the subsequent models:
- DDES (Delayed DES)
- IDDES (Improved DDES)
Motivation and Aim of DES
The motivation and aim of DES models is to decrease the computational cost of massively separated turbulent flows compared to LES (Large Eddy Simulation). The reduction in the cost is achieved by modeling the boundary layers using unsteady RANS (Reynolds-Averaged Navier-Stokes) model [2].
Detached-Eddy Simulation (DES) is a hybrid technique first proposed by Spalart et al. (1997) for prediction of turbulent flows at high Reynolds numbers (see also Spalart 2000). Development of the technique was motivated by estimates which indicate that the computational costs of applying Large-Eddy Simulation (LES) to complete configurations such as an airplane, submarine, or road vehicle are prohibitive. The high cost of LES when applied to complete configurations at high Reynolds numbers arises because of the resolution required in the boundary layers, an issue that remains even with fully successful wall-layer modeling.
In Detached-Eddy Simulation (DES), the aim is to combine the most favorable aspects of the two techniques, i.e., application of RANS models for predicting the attached boundary layers and LES for resolution of time-dependent, three-dimensional large eddies. The cost scaling of the method is then favorable since LES is not applied to resolution of the relatively smaller-structures that populate the boundary layer.
The DES97 model accomplishes the switching between RANS and LES modes by changing the length scale as will be described later in this post.
The model is available in OpenFOAM. I will detail its implementations in various OpenFOAM versions.
OpenFOAM versions
* OpenFOAM v4
* OpenFOAM v1606+
* OpenFOAM v1612+
* OpenFOAM v1706
* OpenFOAM v1712
* OpenFOAM v1806
* OpenFOAM v1812
Governing Equation in v4
Source Code: src/TurbulenceModels/turbulenceModels/
LES/SpalartAllmarasDES/SpalartAllmarasDES.C
429 430 431 432 433 434 435 436 437 438 439 440 441 442 443 |
tmp<fvScalarMatrix> nuTildaEqn ( fvm::ddt(alpha, rho, nuTilda_) + fvm::div(alphaRhoPhi, nuTilda_) - fvm::laplacian(alpha*rho*DnuTildaEff(), nuTilda_) - Cb2_/sigmaNut_*alpha*rho*magSqr(fvc::grad(nuTilda_)) == Cb1_*alpha*rho*Stilda*nuTilda_ - fvm::Sp ( Cw1_*alpha*rho*fw(Stilda, dTilda)*nuTilda_/sqr(dTilda), nuTilda_ ) + fvOptions(alpha, rho, nuTilda_) ); |
Mathematical Expression
\begin{align}
&\;\frac{\partial (\rho \tilde{\nu})}{\partial t} + \frac{\partial (\rho u_j \tilde{\nu})}{\partial x_j} \\
-&\;\frac{1}{\sigma} \left[ \frac{\partial}{\partial x_j} \left( \rho \left( \tilde{\nu} + \nu \right) \frac{\partial \tilde{\nu}}{\partial x_j} \right) + C_{b2} \rho \frac{\partial \tilde{\nu}}{\partial x_i} \frac{\partial \tilde{\nu}}{\partial x_i} \right] \tag{1} \label{eq:sa4} \\
=&\;C_{b1} \rho \tilde{S} \tilde{\nu}\;-\;C_{w1} \rho f_{w} \left( \frac{\tilde{\nu}}{\tilde{d}} \right)^2
\end{align}
In the case of incompressible fluid flow, the density \(\rho = 1\) in Eq. \eqref{eq:sa4}. The third term in the l.h.s. of Eq. \eqref{eq:sa4} is the diffusion term. In the r.h.s, the terms starting from the left represent the production and destruction terms for the eddy viscosity \(\tilde{\nu}\). The destruction term is changed from the Spalart-Allmaras RANS model by replacing the distance to the wall \(d\) with a new length scale \(\tilde{d}\) [2]
\begin{equation}
\tilde{d} = {\rm min} \left( l_{RANS}, l_{LES} \right) = {\rm min} \left( d, C_{DES}\Delta \right) \tag{2} \label{eq:dTilda}
\end{equation}
where \(\Delta\) is the chosen measure of grid spacing and \(C_{DES}\) is a calibration constant and its default value is 0.65.
When balanced with the production term, this term adjusts the eddy viscosity to scale with the local deformation rate \(S\) and \(d\): \(\tilde{\nu} \propto Sd^2\). Subgrid-scale (SGS) eddy viscosities scale with \(S\) and the grid spacing \(\Delta\), i.e., \(\nu_{SGS} \propto S\Delta^2\). A subgrid-scale model within the S-A formulation can then be obtained by replacing \(d\) with a length scale \(\Delta\) directly proportional to the grid spacing.
Governing Equation in v1606+ – v1812
Source Code: src/TurbulenceModels/turbulenceModels/
DES/SpalartAllmarasDES/SpalartAllmarasDES.C
540 541 542 543 544 545 546 547 548 549 550 551 552 553 554 555 |
tmp<fvScalarMatrix> nuTildaEqn ( fvm::ddt(alpha, rho, nuTilda_) + fvm::div(alphaRhoPhi, nuTilda_) - fvm::laplacian(alpha*rho*DnuTildaEff(), nuTilda_) - Cb2_/sigmaNut_*alpha*rho*magSqr(fvc::grad(nuTilda_)) == Cb1_*alpha*rho*Stilda*nuTilda_*(scalar(1) - ft2) - fvm::Sp ( (Cw1_*fw(Stilda, dTilda) - Cb1_/sqr(kappa_)*ft2) *alpha*rho*nuTilda_/sqr(dTilda), nuTilda_ ) + fvOptions(alpha, rho, nuTilda_) ); |
Mathematical Expression
\begin{align}
&\;\frac{\partial (\rho \tilde{\nu})}{\partial t} + \frac{\partial (\rho u_j \tilde{\nu})}{\partial x_j} \\
-&\;\frac{1}{\sigma} \left[ \frac{\partial}{\partial x_j} \left( \rho \left( \tilde{\nu} + \nu \right) \frac{\partial \tilde{\nu}}{\partial x_j} \right) + C_{b2} \rho \frac{\partial \tilde{\nu}}{\partial x_i} \frac{\partial \tilde{\nu}}{\partial x_i} \right] \tag{3} \label{eq:sa1606} \\
=&\;C_{b1} \rho \tilde{S} \tilde{\nu} \left( 1 – \color{#dd4b39}{f_{t2}} \right)\;-\;\rho \left( C_{w1} f_{w} – \frac{C_{b1}}{\kappa^2} \color{#dd4b39}{f_{t2}} \right) \left( \frac{\tilde{\nu}}{\tilde{d}} \right)^2
\end{align}
The \(f_{t2}\) laminar suppression term is included in com versions (v1606+ – v1812) as against v4. The subscript \(t\) stands for “trip”. The production term is multiplied by \(1-f_{t2}\) so that \(\tilde{\nu} = 0\) is a stable solution when solving a transitional flow where both laminar and turbulent regions exist in one computational domain. The destruction term is also altered in accordance with the change in the production term in order to balance the energy budget near the wall [3]. For fully turbulent flows at high Reynolds numbers, the \(f_{t2}\) term takes small values and makes essentially no difference between the above two equations.
In addition to this difference, the low-Reynolds number correction option [4] is available only in com versions (v1606+ – v1812):
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 |
FoamFile { version 2.0; format ascii; class dictionary; location "constant"; object turbulenceProperties; } // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // simulationType LES; LES { turbulence on; LESModel SpalartAllmarasDES; // Optional entries SpalartAllmarasDESCoeffs { // Apply low-Reynolds number correction; default = yes lowReCorrection yes; } } |
152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 |
template<class BasicTurbulenceModel> tmp<volScalarField> SpalartAllmarasDES<BasicTurbulenceModel>::psi ( const volScalarField& chi, const volScalarField& fv1 ) const { tmp<volScalarField> tpsi ( new volScalarField ( IOobject ( type() + ":psi", this->time().timeName(), this->mesh(), IOobject::NO_READ, IOobject::NO_WRITE ), this->mesh(), dimensionedScalar("one", dimless, 1) ) ); if (lowReCorrection_) { volScalarField& psi = tpsi.ref(); const volScalarField fv2(this->fv2(chi, fv1)); const volScalarField ft2(this->ft2(chi)); psi = sqrt ( min ( scalar(100), (1 - Cb1_/(Cw1_*sqr(kappa_)*fwStar_)*(ft2 + (1 - ft2)*fv2)) /max(SMALL, (fv1*max(scalar(1e-10), 1 - ft2))) ) ); } return tpsi; } |
If the low-Re number correction option is deactivated (lowReCorrection false;) and the value of \(C_{t3}\) is set to 0 in v1606+ – v1812, the same calculation can be done as v4.
Low-Reynolds Number Term Correction \(\Psi\) Function
The low-Reynolds number terms in the S-A RANS model were designed to account for wall proximity, “measured” by the ratio of the eddy and molecular viscosities \(\nu_t/\nu\) (an expedient to avoid using the friction velocity). Most RANS models have similar terms, also based on various turbulent Reynolds numbers. However, in the LES mode of DES, the subgrid eddy viscosity (normalized by the molecular viscosity) decreases with grid refinement and decrease of the flow Reynolds number. At some point, standard DES will mis-interpret the situation as “wall proximity” and lower the eddy viscosity excessively, relative to the ambient velocity and length scales, through the \(f_v\) and \(f_t\) functions of the S-A model. Although this deficiency shows up only at relatively low Reynolds numbers and/or with very fine grids, it is still desirable to eliminate it, and to do so without zonal instructions from the user.
\begin{equation}
\tilde{d} = {\rm min} \left( d, \Psi C_{DES}\Delta \right) \tag{4} \label{eq:dTildaPsi}
\end{equation}
\begin{equation}
\Psi^2 = {\rm min} \left[ 10^2, \frac{1-\frac{C_{b1}}{C_{w1}\kappa^2 f_w^*} \left( f_{t2} + \left( 1 – f_{t2} \right) f_{v2} \right)}{f_{v1} {\rm max}\left[ 10^{-10}, 1-f_{t2} \right]} \right] \tag{5} \label{eq:Psi}
\end{equation}
Requirement for Computational Mesh
The DES97 model is based on the assumption that the wall-parallel grid spacing near the wall \(\Delta_{||}\) is in excess of the boundary layer thickness \(\delta\) by a good margin so that the entire boundary layer is handled by RANS model. As we can see from Eq. \eqref{eq:dTilda}, the new length scale of DES97 model \(\tilde{d}\) depends only on the geometric features of computational mesh and not on the computed flow fields at all.
If we choose the following measure of grid spacing that is maxDeltaxyz option in OpenFOAM
\begin{equation}
\Delta \equiv {\rm max} \left( \Delta x, \Delta y, \Delta z \right), \tag{6} \label{eq:delta}
\end{equation}
the wall-parallel grid spacing \(\Delta_{||}\) should meet the following inequality to satisfy the above mesh requirement
\begin{equation}
\Delta_{||} \approx \Delta \geq \frac{d}{C_{DES}} \geq \frac{\delta}{C_{DES}}, \tag{7} \label{eq:deltaInequality}
\end{equation}
where we assume that \(y\) is the wall normal direction and \(\Delta_{||} \approx \Delta_x \approx \Delta_z\).
There are several problems that DES97 model suffers from.
- Grey area
- Modeled-Stress Depletion (MSD)
Delayed DES and Zonal DES are the solutions to MSD.
- Grid-induced Separation (GIS)
MSD is the root cause of GIS.
Laminar Suppression Term ft2
In the excellent Turbulence Modeling Resource page [5] of NASA Langley Research Center, the various deformations of the Spalart-Allmaras turbulence (RANS) model are clearly documented.
According to the naming method, the “Standard” Spalart-Allmaras One-Equation Model (SA) is implemented in OpenFOAM v1606+ – v1812 and the Spalart-Allmaras One-Equation Model without \(f_{t2}\) Term (SA-noft2) is implemented in v4.
Spalart-Allmaras One-Equation Model without \(f_{t2}\) Term (SA-noft2)
Many implementations of Spalart-Allmaras ignore the \(f_{t2}\) term, which was a numerical fix in the original model in order to make zero a stable solution to the equation with a small basin of attraction (thus slightly delaying transition so that the trip term could be activated appropriately). It is argued that if the trip is not used, then \(f_{t2}\) is not necessary. The equations are the same as for the “standard” version (SA), except that the term \(f_{t2}\) does not appear at all (i.e., ct3=0 instead of 1.2).
Common Definitions between the Two Versions
38 39 40 41 42 |
template<class BasicTurbulenceModel> tmp<volScalarField> SpalartAllmarasDES<BasicTurbulenceModel>::chi() const { return nuTilda_/this->nu(); } |
\begin{equation}
\chi = \frac{\tilde{\nu}}{\nu} \tag{8} \label{eq:chi}
\end{equation}
45 46 47 48 49 50 51 52 53 |
template<class BasicTurbulenceModel> tmp<volScalarField> SpalartAllmarasDES<BasicTurbulenceModel>::fv1 ( const volScalarField& chi ) const { const volScalarField chi3("chi3", pow3(chi)); return chi3/(chi3 + pow3(Cv1_)); } |
\begin{equation}
f_{v1} = \frac{\chi^3}{\chi^3 + C_{v1}^3} \tag{9} \label{eq:fv1}
\end{equation}
56 57 58 59 60 61 62 63 64 |
template<class BasicTurbulenceModel> tmp<volScalarField> SpalartAllmarasDES<BasicTurbulenceModel>::fv2 ( const volScalarField& chi, const volScalarField& fv1 ) const { return 1.0 - chi/(1.0 + chi*fv1); } |
\begin{equation}
f_{v2} = 1\;-\;\frac{\chi}{1+\chi f_{v1}} \tag{10} \label{eq:fv2}
\end{equation}
77 78 79 80 81 82 83 84 |
template<class BasicTurbulenceModel> tmp<volScalarField> SpalartAllmarasDES<BasicTurbulenceModel>::Omega ( const volTensorField& gradU ) const { return sqrt(2.0)*mag(skew(gradU)); } |
\begin{align}
W_{ij} &= \frac{1}{2}\left( \nabla \boldsymbol{u} – (\nabla \boldsymbol{u})^T \right) = \frac{1}{2}\left( \frac{\partial u_j}{\partial x_i} – \frac{\partial u_i}{\partial x_j} \right) \tag{11a} \\
\Omega &= \sqrt{2 W_{ij} \; W_{ij}} \tag{11b} \label{eq:Omega}
\end{align}
87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 |
template<class BasicTurbulenceModel> tmp<volScalarField> SpalartAllmarasDES<BasicTurbulenceModel>::Stilda ( const volScalarField& chi, const volScalarField& fv1, const volScalarField& Omega, const volScalarField& dTilda ) const { return ( max ( Omega + fv2(chi, fv1)*nuTilda_/sqr(kappa_*dTilda), Cs_*Omega ) ); } |
\begin{equation}
\tilde{S} = {\rm max} \left[ \Omega + \frac{\tilde{\nu}}{\kappa^2 \tilde{d}^2} f_{v2}, C_{s} \Omega \right] \tag{12} \label{eq:Stilda}
\end{equation}
108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 |
template<class BasicTurbulenceModel> tmp<volScalarField> SpalartAllmarasDES<BasicTurbulenceModel>::r ( const volScalarField& nur, const volScalarField& Omega, const volScalarField& dTilda ) const { tmp<volScalarField> tr ( min ( nur /( max ( Omega, dimensionedScalar("SMALL", Omega.dimensions(), SMALL) ) *sqr(kappa_*dTilda) ), scalar(10) ) ); tr.ref().boundaryFieldRef() == 0.0; return tr; } |
\begin{equation}
r = {\rm min} \left[ \frac{\tilde{\nu}}{\Omega \kappa^2 \tilde{d}^2}, 10 \right] \tag{13} \label{eq:r}
\end{equation}
138 139 140 141 142 143 144 145 146 147 148 149 |
template<class BasicTurbulenceModel> tmp<volScalarField> SpalartAllmarasDES<BasicTurbulenceModel>::fw ( const volScalarField& Omega, const volScalarField& dTilda ) const { const volScalarField r(this->r(nuTilda_, Omega, dTilda)); const volScalarField g(r + Cw2_*(pow6(r) - r)); return g*pow((1 + pow6(Cw3_))/(pow6(g) + pow6(Cw3_)), 1.0/6.0); } |
\begin{align}
g &= r + C_{w2} \left( r^6 – r \right) \tag{14a} \\
f_{w} &= g \left[ \frac{1 + C_{w3}^6}{g^6 + C_{w3}^6} \right] ^{\frac{1}{6}} \tag{14b} \label{eq:fw}
\end{align}
Different Definitions between the Two Versions
152 153 154 155 156 157 158 159 160 161 162 163 |
template<class BasicTurbulenceModel> tmp<volScalarField> SpalartAllmarasDES<BasicTurbulenceModel>::dTilda ( const volScalarField& chi, const volScalarField& fv1, const volTensorField& gradU ) const { tmp<volScalarField> tdTilda(CDES_*this->delta()); min(tdTilda.ref().ref(), tdTilda(), y_); return tdTilda; } |
199 200 201 202 203 204 205 206 207 208 209 210 |
template<class BasicTurbulenceModel> tmp<volScalarField> SpalartAllmarasDES<BasicTurbulenceModel>::dTilda ( const volScalarField& chi, const volScalarField& fv1, const volTensorField& gradU ) const { tmp<volScalarField> tdTilda(psi(chi, fv1)*CDES_*this->delta()); min(tdTilda.ref().dimensionedInternalField(), tdTilda(), y_); return tdTilda; } |
67 68 69 70 71 72 73 74 |
template<class BasicTurbulenceModel> tmp<volScalarField> SpalartAllmarasDES<BasicTurbulenceModel>::ft2 ( const volScalarField& chi ) const { return Ct3_*exp(-Ct4_*sqr(chi)); } |
\begin{equation}
f_{t2} = C_{t3} e^{-C_{t4} \chi^2} \tag{15} \label{eq:ft2}
\end{equation}
Other Topics
- Function Objects – DESModelRegions
- Convection Scheme – DEShybrid
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] K. D. Squires, DETACHED-EDDY SIMULATION: CURRENT STATUS AND PERSPECTIVES.
[3] P. R. Spalart and S. R. Allmaras, A one-equation turbulence model for aerodynamic flows. La Rech. Aerospatiale 1, 5-21, 1994.
[4] 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.
[5] Turbulence Modeling Resource. Available at: https://turbmodels.larc.nasa.gov/index.html [Accessed: 10 October 2016].
[6] P. R. Spalart and C. L. Rumsey, Effective inflow conditions for turbulence models in aerodynamic calculations. AIAA Journal 45(10), 2544-2553, 2007.
[7] P. R. Spalart et al., THE USES OF DES: NATURAL, EXTENDED, AND IMPROPER. DESider, July 14-15 2005, Stockholm.