Editing Large eddy simulation (section)

==== Functional (eddy–viscosity) models ====

Functional models are simpler than structural models, focusing only on dissipating energy at a rate that is physically correct.  These are based on an artificial eddy viscosity approach, where the effects of turbulence are lumped into a turbulent viscosity.  The approach treats dissipation of kinetic energy at sub-grid scales as analogous to molecular diffusion.  In this case, the deviatoric part of <math>\tau_{ij}</math> is modeled as:

:<math>
\tau_{ij}^r - \frac{1}{3} \tau_{kk} \delta_{ij} = -2 \nu_\mathrm{t} \bar{S}_{ij}
</math>

where <math>\nu_\mathrm{t}</math> is the turbulent eddy viscosity and <math>\bar{S}_{ij} = \frac{1}{2} \left( \frac{\partial \bar{u}_i }{\partial x_j} + \frac{\partial \bar{u}_j}{ \partial x_i} \right)</math> is the rate-of-strain tensor.

Based on dimensional analysis, the eddy viscosity must have units of  <math>\left[ \nu_\mathrm{t} \right] = \frac{\mathrm{m^2}}{\mathrm{s}}</math>.  Most eddy viscosity SGS models model the eddy viscosity as the product of a characteristic length scale and a characteristic velocity scale.

===== Smagorinsky–Lilly model =====

The first SGS model developed was the Smagorinsky–Lilly SGS model, which was developed by [[Joseph Smagorinsky|Smagorinsky]]<ref name="Smagorinsky_1963" /> and used in the first LES simulation by Deardorff.<ref name="Deardorff_1970" />  It models the eddy viscosity as:

:<math>\nu_\mathrm{t} 
= C \Delta^2\sqrt{2\bar{S}_{ij}\bar{S}_{ij}} 
= C \Delta^2 \left| \bar{S} \right|
</math>

where <math>\Delta</math> is the grid size and <math>C</math> is a constant.

This method assumes that the energy production and dissipation of the small scales are in equilibrium - that is, <math>\epsilon = \Pi</math>.

===== The Dynamic Model (Germano et al. and beyond)=====

Germano et al.<ref name="Germano_1991">{{cite journal
|title=A dynamic subgrid-scale eddy viscosity model
|author1=Germano, M.
|author2=Piomelli, U.
|author3=Moin, P.
|author4=Cabot, W.
|journal=[[Physics of Fluids|Physics of Fluids A]]
|volume=3
|pages=1760–1765
|year=1991
|doi=10.1063/1.857955|bibcode = 1991PhFlA...3.1760G
|issue=7
|s2cid=55719851
|authorlink3=Parviz Moin }}
</ref> identified a number of studies using the Smagorinsky model that each found different values for the Smagorinsky constant <math>C</math> for different flow configurations.  In an attempt to formulate a more universal approach to SGS models, Germano et al. proposed a dynamic Smagorinsky model, which utilized two filters: a grid LES filter, denoted <math>\overline{f}</math>, and a test LES filter, denoted <math>\hat{f}</math> for any turbulent field <math>f</math>. The test filter is larger in size than the grid filter and adds an additional smoothing of the turbulence field over the already smoothed fields represented by the LES. Applying the test filter to the LES equations (which are obtained by applying the "grid" filter to Navier-Stokes equations) results in a new set of equations that are identical in form but with the SGS stress <math>\tau_{ij} = \overline{u_{i} u_{j}} - \bar{u}_{i} \bar{u}_{j}</math> replaced by <math>T_{ij} = \widehat{\overline{u_{i} u_{j}}} - \hat{\bar{u}}_{i} \hat{\bar{u}}_{j}</math>. Germano ''et'' al. noted that even though neither <math>\tau_{ij}</math> nor <math>T_{ij}</math> can be computed exactly because of the presence of unresolved scales, there is an exact relation connecting these two tensors. This relation, known as the Germano identity is 
<math>
L_{ij} = T_{ij} - \hat{\tau}_{ij}.
</math>
Here <math> L_{ij} = \widehat{\bar{u}_{i} \bar{u}_{j}} - \widehat{\bar{u}_{i}} \widehat{\bar{u}_{j}}</math> can be explicitly evaluated as it involves only the filtered velocities and the operation of test filtering. The significance of the identity is that if one assumes that turbulence is self similar so that the SGS stress at the grid and test levels have the same form <math>\tau_{ij} - (\tau_{kk}/3)\delta_{ij} = - 2 C \Delta^{2} |\bar{S}_{ij}| \bar{S}_{ij}</math> and <math>T_{ij} - (T_{kk}/3)\delta_{ij} = - 2 C \hat{\Delta}^{2} |\hat{\bar{S}}_{ij}| \hat{\bar{S}}_{ij}</math>, then the Germano identity provides an equation from which the Smagorinsky coefficient <math>C</math> (which is no longer a 'constant') can potentially be determined. 
[Inherent in the procedure is the assumption that the coefficient <math>C</math> is invariant of scale (see review
<ref name="MeneveauKatz_2000">{{cite journal
|title=Scale-Invariance and Turbulence Models for Large-Eddy Simulation
|author1=Meneveau, C.
|author2=Katz, J.
|journal=Annu. Rev. Fluid Mech.
|year=2000
|volume=32
|issue=1
|pages=1–32
|doi=10.1146/annurev.fluid.32.1.1
|bibcode = 2000AnRFM..32....1M }}</ref>)]. 
In order to do this, two additional steps were introduced in the original formulation. First, one assumed that even though <math>C</math> was in principle variable, the variation was sufficiently slow that it can be moved out of the filtering operation <math> \widehat{C (.)} = C \widehat{(.)} </math>. Second, since <math>C</math> was a scalar, the Germano identity was contracted with a second rank tensor (the rate of strain tensor was chosen) to convert it to a scalar equation from which <math>C</math> could be determined.
Lilly
<ref name="Lilly_1992">{{cite journal
|title=A proposed modification of the Germano subgrid-scale closure method
|journal=Physics of Fluids A
|year=1992
|volume=4
|issue=3
|pages=633–636
|doi=10.1063/1.858280
|last1=Lilly
|first1=D. K.|bibcode = 1992PhFlA...4..633L }}</ref>
found a less arbitrary and therefore more satisfactory approach for obtaining C from the tensor identity. He noted that the Germano identity required the satisfaction of nine equations at each point in space (of which only five are independent) for a single quantity <math>C</math>. The problem of obtaining <math>C</math> was therefore over-determined. He proposed therefore that <math>C</math> be determined using a least square fit by minimizing the residuals. This results in 
:<math>
C  = \frac{ L_{ij} m_{ij} }{ m_{kl} m_{kl} }.
</math>
Here 
: <math>
m_{ij} = \alpha_{ij} - \widehat{\beta}_{ij} 
</math> 
and for brevity 
<math> \alpha_{ij} = - 2 \hat{\Delta}^{2} | \hat{\bar{S}} | \hat{\bar{S}}_{ij}  </math>, 
<math> \beta_{ij} = - 2 \Delta^2 | \bar{S} | \bar{S}_{ij}  </math> 
Initial attempts to implement the model in LES simulations proved unsuccessful. First, the computed coefficient 
was not at all "slowly varying" as assumed and varied as much as any other turbulent field. Secondly, 
the computed <math>C</math> could be positive as well as negative. The latter fact in itself should not be regarded as a 
shortcoming as a priori tests using filtered DNS fields have shown that the local subgrid dissipation rate 
<math> - \tau_{ij} \bar{S}_{ij}</math> in a turbulent field is almost as likely to be negative as it is positive even though the integral over the fluid domain is always positive representing a net dissipation of energy in the large scales. A slight preponderance of positive values as opposed to strict positivity of the eddy-viscosity results in the observed net dissipation. This so-called "backscatter" of energy from small to large scales indeed corresponds to negative C values in the Smagorinsky model. Nevertheless, the Germano-Lilly formulation was found not to result in stable calculations. An ad hoc measure was adopted by averaging the numerator and denominator over homogeneous directions (where such directions exist in the flow) 
: <math>
C = \frac{ 
\left\langle L_{ij} m_{ij} \right\rangle
}{ 
\left\langle m_{kl} m_{kl} \right\rangle
}.
</math>
When the averaging involved a large enough statistical sample that the computed <math>C</math> was positive (or at 
least only rarely negative) stable calculations were possible. Simply setting the negative values to zero (a procedure called "clipping") with or without the averaging also resulted in stable calculations. 
Meneveau proposed <ref name="Meneveauetal_1996">{{cite journal
|title=A Lagrangian dynamic subgrid-scale model of turbulence
|author1=Meneveau, C.
|author2=Lund, T. S.
|author3=Cabot, W. H.
|journal=J. Fluid Mech.
|year=1996
|volume=319
|issue=1
|pages=353–385
|doi=10.1017/S0022112096007379
|bibcode = 1996JFM...319..353M |hdl=2060/19950014634
|s2cid=122183534
|hdl-access=free
}}</ref>
an averaging over Lagrangian fluid trajectories with an exponentially decaying "memory". This can be applied to problems lacking homogeneous directions and can be stable if the effective time over which the averaging is done is long enough and yet not so long as to smooth out spatial inhomogenieties of interest.

Lilly's modification of the Germano method followed by a statistical averaging or synthetic removal of negative viscosity regions seems ad hoc, even if it could be made to "work".  An alternate formulation of the least square minimization procedure known as the "Dynamic Localization Model" (DLM) was suggested by 
Ghosal et al.<ref name="Ghosal_1995">{{cite journal
|title=A dynamic localization model for large-eddy simulation of turbulent flows
|journal=Journal of Fluid Mechanics 
|year=1995
|volume=286
|issue=
|pages=229–255
|doi=10.1017/S0022112095000711
|last1=Ghosal
|first1=S.
|last2=Lund
|first2=T.S.
|last3=Moin
|first3=P.
|last4=Akselvoll
|first4=K.|bibcode=1995JFM...286..229G 
|s2cid=124586994 
}}</ref>
In this approach one first defines a quantity 
:<math>
E_{ij} = L_{ij} - T_{ij} + \hat{\tau}_{ij}
</math>
with the tensors <math>\tau_{ij}</math> and <math>T_{ij}</math> replaced by the appropriate SGS model. This tensor then represents the amount by which the subgrid model fails to respect the Germano identity at each spatial location. In Lilly's approach,  <math>C</math> is then pulled out of the hat operator 
:<math>
 \widehat{C (.)} =  C \widehat{(.)}
</math>
making <math>E_{ij}</math> an algebraic function of <math>C</math> which is then determined by requiring that
<math>E_{ij} E_{ij}</math> considered as a function of C have the least possible value.
However, since the <math>C</math> thus obtained turns out to be just as variable as any other fluctuating quantity in turbulence, the original assumption of the constancy of <math>C</math> cannot be justified a posteriori. In the DLM approach one avoids this inconsistency by not invoking the step of removing 
C from the test filtering operation. Instead, one defines a global error over the entire flow domain by the quantity 
:<math>
 E [ C ] = \int E_{ij} E_{ij} dV
</math>
where the integral ranges over the whole fluid volume. This global error <math>E[C(x,y,z,t)]</math> is then a functional of the spatially varying function <math>C(x,y,z,t)</math> (here the time instant, <math>t</math>, is fixed and therefore appears just as a parameter) which is determined so as to minimize this functional. The solution to this variational problem is that <math>C</math> must satisfy a Fredholm integral equation of the second kind 
:<math> 
C (\boldsymbol{x}) = f ( \boldsymbol{x} ) + 
\int K(\boldsymbol{x}, \boldsymbol{y}) C ( \boldsymbol{y} ) d\boldsymbol{y}
</math> 
where the functions <math>K(\boldsymbol{x}, \boldsymbol{y})</math> and <math>f ( \boldsymbol{x} )</math> are defined in terms of the resolved fields <math>L_{ij},\alpha_{ij},\beta_{ij}</math> and are therefore known at each time step and the integral ranges over the whole fluid domain. The integral equation is solved numerically by an iteration procedure and convergence was found to be generally rapid if used with a pre-conditioning scheme. Even though this variational approach removes an inherent inconsistency in Lilly's approach, the <math>C(x,y,z,t)</math> obtained from the integral equation still displayed the instability associated with negative viscosities. This can be resolved by insisting that <math>E[C]</math> be minimized subject to the constraint <math>C(x,y,z,t) \geq 0</math>. This leads to an equation for <math>C</math> that is nonlinear 
:<math> 
C (\boldsymbol{x}) = \left[ f ( \boldsymbol{x} ) + 
\int K(\boldsymbol{x}, \boldsymbol{y}) C ( \boldsymbol{y} ) d\boldsymbol{y} \right]_{+}
</math> 
Here the suffix + indicates the "positive part of" that is, <math> x_{+} = (x + |x|)/2 </math>. Even though this superficially looks like "clipping" it is not an ad hoc scheme but a bonafide solution of the constrained variational problem. This DLM(+) model was found to be stable and yielded excellent results for forced and decaying isotropic turbulence, channel flows and a variety of other more complex geometries. If a flow happens to have homogeneous directions (let us say the directions x and z) then one can introduce the ansatz  
<math> C = C(y,t) </math>. The variational approach then immediately yields Lilly's result with averaging over homogeneous directions without any need for ad hoc modifications of a prior result.

One shortcoming of the DLM(+) model was that it did not describe backscatter which is known to be a real "thing" from analyzing DNS data. Two approaches were developed to address this. In one approach due to Carati et al. 
<ref name="Carati_1995">{{cite journal
|title=On the representation of backscatter in Dynamic Localization models
|journal=Physics of Fluids
|year=1995
|volume=7
|issue=3
|pages=606–616
|doi= 10.1063/1.868585 
|last1=Carati
|first1=D.
|last2=Ghosal
|first2=S.
|last3=Moin
|first3=P.
|bibcode=1995PhFl....7..606C
|doi-access=free
}}</ref>
a fluctuating force with amplitude determined by the fluctuation-dissipation theorem is added in 
analogy to Landau's theory of fluctuating hydrodynamics. In the second approach, one notes that 
any "backscattered" energy appears in the resolved scales only at the expense of energy in the subgrid 
scales. The DLM can be modified in a simple way to take into account this physical fact so as to allow 
for backscatter while being inherently stable. This k-equation version of the DLM, DLM(k) replaces 
<math> \Delta | \bar{S} | </math> in the Smagorinsky eddy viscosity model by <math> \sqrt{k} </math> as an appropriate velocity scale. The procedure for determining <math>C</math> remains identical to the "unconstrained" version except that the tensors <math> \alpha_{ij} =  - 2 \hat{\Delta} \sqrt{K} \hat{\bar{S}}_{ij} </math>,
<math> \beta_{ij} =  - 2 \hat{\Delta} \sqrt{k} \bar{S}_{ij} </math> where the sub-test scale kinetic 
energy K is related to the subgrid scale kinetic energy k by <math> K = k + L_{ii}/2 </math> 
(follows by taking the trace of the Germano identity). To determine k we now use a transport equation 
:<math> 
\frac{\partial k}{\partial t} + u_{j} \frac{\partial k}{\partial x_{j}} = 
- \tau_{ij} \bar{S}_{ij} - \frac{C_{*}}{\Delta} k^{3/2} + 
\frac{\partial }{\partial x_j} \left( D \Delta \sqrt{k} \frac{\partial k }{\partial x_j} \right) 
+ \nu \frac{\partial^{2} k }{\partial x_j \partial x_j}
</math>
where <math>\nu</math> is the kinematic viscosity and <math>C_{*},D</math> are positive coefficients 
representing kinetic energy dissipation and diffusion respectively. These can be determined following the dynamic 
procedure with constrained minimization as in DLM(+). This approach, though more expensive to implement than the DLM(+) was found to be stable and resulted in good agreement with experimental data for a variety of flows
tested. Furthermore, it is mathematically impossible for the DLM(k) to result in an unstable computation as the sum of the large scale and SGS energies is non-increasing by construction. Both of these approaches incorporating backscatter works well. They yield models that are slightly less dissipative with somewhat improved performance over the DLM(+). The DLM(k) model additionally yields the subgrid kinetic energy, which may be a physical quantity of interest. These improvements are achieved at a somewhat increased cost in model implementation.

The Dynamic Model originated at the 1990 [https://ctr.stanford.edu/ctr-summer-program Summer Program] of the [[Center for Turbulence Research]] (CTR) at [[Stanford University]]. A series of "CTR-Tea" seminars celebrated the [https://ctr.stanford.edu/event-type/ctr-tea 30th Anniversary] {{Webarchive|url=https://web.archive.org/web/20221030050059/https://ctr.stanford.edu/event-type/ctr-tea |date=2022-10-30 }} of this important milestone in turbulence modeling.