Matrix Exponential-Based Closures for the Turbulent Subgrid-Scale Stress Tensor

0.33 MB·English

by  Yi Li

Matrix Exponential-Based Closures for the Turbulent Subgrid-Scale Stress Tensor
Yi Li1,2, Laurent Chevillard1,3, Gregory Eyink4 and Charles Meneveau1
1Department of Mechanical Engineering and Center of Environmental and Applied Fluid Mechanics,
The Johns Hopkins University, Baltimore, MD, 21218
2Department of Applied Mathematics, The University of Sheffield H23B,
Hicks Building, Hounsfield Road Sheffield, S3 7RH, UK
3Laboratoire de Physique de l’E´cole Normale Sup´erieure de Lyon,
CNRS, Universit´e de Lyon, 46 all´ee d’Italie F-69007 Lyon, France
4Department of Applied Mathematics & Statistics,
and Center of Environmental and Applied Fluid Mechanics,
9 The Johns Hopkins University, Baltimore, MD, 21218
0
0 Twoapproachesforclosingtheturbulencesubgrid-scalestresstensorintermsofmatrixexponen-
2 tials are introduced and compared. The first approach is based on a formal solution of the stress

n transport equation in which the production terms can be integrated exactly in terms of matrix
a exponentials. This formal solution of the subgrid-scale stress transport equation is shown to be
J useful to explore special cases, such as the response to constant velocity gradient, but neglecting
1 pressure-strain correlations and diffusion effects. The second approach is based on an Eulerian-
2 Lagrangian change of variables, combined with the assumption of isotropy for the conditionally
averaged Lagrangian velocity gradient tensor and with the ‘Recent Fluid Deformation’ (RFD) ap-

] proximation. It is shown that both approaches lead to the same basic closure in which the stress
n
tensorisexpressedastheproductofthematrixexponentialoftheresolvedvelocitygradienttensor
y
multipliedbyitstranspose. Short-timeexpansionsofthematrixexponentialsareshowntoprovide
d
an eddy-viscosity term and particular quadratic terms, and thus allow a reinterpretation of tradi-
-
u tional eddy-viscosity and nonlinear stress closures. The basic feasibility of the matrix-exponential
l closure is illustrated by implementing it successfully in Large Eddy Simulation of forced isotropic
f
turbulence. Thematrix-exponentialclosure employsthedrastic approximation of entirely omitting
.
s the pressure-strain correlation and other ‘nonlinear scrambling’ terms. But unlike eddy-viscosity
c
closures, the matrix exponential approach provides a simple and local closure that can be derived
i
s directlyfromthestresstransportequationwiththeproductionterm,andusingphysicallymotivated
y
assumptions about Lagrangian decorrelation and upstream isotropy.
h
p PACSnumbers:
[

3
I. INTRODUCTION Fluid Deformation (RFD) closure – has been proposed
v
1 [10, 11, 12]. In this approach, a change of variables is
8 madeexpressingspatialgradientsintermsofLagrangian
One of the most basic challenges in turbulence mod-
7 gradients(e.g. howdoesavariableatthepresentlocation
eling is the need for closures for the fluxes associated
3 vary if we change the initial position of the fluid particle
4. with unresolved turbulent fluctuations. In the context at an earlier time). Then the assumption of isotropy is
of Large Eddy Simulation (LES), closures are required
0 introduced for the Lagrangiangradienttensors. This as-
for the subgrid-scale (SGS) stress tensor [1, 2]. Tradi-
7 sumption allows simpler isotropic forms to be used, and
0 tional closures involve mostly algebraic expressions re- is argued to be justified based on Lagrangian decorre-
: latingthestresstensortopowersofthevelocitygradient
v lation effects. Deviations from isotropy at the present
tensor. More elaborate approaches using separate trans-
i location for the Eulerian gradient tensors develop as a
X port equations have sometimes also been employed, al-
resultoffluidmaterialdeformationalongtheLagrangian
r though these tend to be significantly more costly in the trajectory. More traditionally, the Lagrangian time evo-
a context of LES. Closures expressing the stress in terms
lution of the stress tensor following fluid particles can
ofthematrixexponentialfunctiondonotappeartohave
bederivedbytakingappropriatemomentsoftheNavier-
received much attention in the literature. The objective
Stokesequations. Inthispaperweexaminebothofthese
of the present work is to identify and discuss two sep-
approachestoformulatemodelsfortheSGSstresstensor
arate paths that lead to such closures. Both paths are
in turbulence in the context of LES.
based on the Lagrangiandynamics of turbulence, i.e. on
Adescriptionofsmall-scalestructureofturbulencebe-
an understanding of the evolution of turbulence as one
gins with the Navier-Stokes equations of an incompress-
follows fluid-particle paths in time.
ible fluid of velocity u:
The use of Lagrangianconcepts in turbulent flows has
du ∂u
alonghistory[3,4]and,inrecentyears,hasseenrenewed = +(u ∇)u= ∇p+ν∇2u , (1)
dt ∂t · −
interestformodeling[5,6,7,8,9]. Amongothers,anew
model for the pressure-Hessian tensor based on the re- whered/dtstandsfortheLagrangianmaterialderivative,
cent Lagrangianevolutionof fluid elements – the Recent p the pressure divided by the density of the fluid and
2
ν the kinematic viscosity. Because of incompressibility, resulting model is shown to be expressible compactly in
the velocity gradient tensor A =∂u /∂x must remain terms of matrix exponentials as well. Differences and
ij i j
trace-free, i.e. A = 0, and the pressure field is the similaritiesbetweenthe RFD andtransportequationso-
ii
solution of the Poissonequation 2p= A A . lutionsarediscussed. In V,thematrix-exponentialsolu-
lk kl
∇ − §
In the framework of LES, the SGS stress tensor is de- tionsareexpandedforshorttimes. Theexpansionsallow
fined using the filtering approach [13, 15], to establish relationships to traditional eddy-viscosity
and nonlinear closure models in turbulence. In VI the
τij =uiuj uiuj . (2) matrix-exponentialclosureisimplementedinamo§stsim-
−
ple flow to illustrate its feasibility and cost.
An overbar denotes spatial filtering at a scale ∆ and is
formallygivenbya convolutionwithanonnegative,spa-
tially well-localizedfiltering function (r) ofcharacteris- II. SOLUTION TO STRESS TRANSPORT
G
tic size ∆, with unit integral (r)dr=1, namely EQUATION USING MATRIX EXPONENTIALS
G
R
u(x,t)= (r)u(x+r,t)dr . Equation4 is of the formof the “time-dependent Lya-
G
Z punov equation”, if the tensor Φ’s implicit dependencies
TheSGStensorτ entersinthedynamicsofthefiltered upon the velocity fluctuations and the stresstensor were
velocity u as it can be seen when applying the filtering disregarded (in reality, Φ depends upon small-scale ve-
locity fluctuations and thus the full equation is highly
procedure to the Navier-Stokes equations (Eq. (1)),
non-linear and non-local). The formal solution of the
Du ∂u Lyapunov equation in terms of matrix exponentials has
= +(u ∇)u= ∇p+ν∇2u ∇ τ , (3)
Dt ∂t · − − · been found useful in a number of other fields: principal
oscillation pattern analysis [16], mechanics of finite de-
where D/Dt stands for the Lagrangian material deriva- formations [17], and fluctuation-dissipation theorems for
tive with u as the advecting velocity, and p the filtered stochastic linear systems [18]. In the context of the SGS
pressure divided by the density of the fluid. Because stress transport equation the solution at time t (start-
of incompressibility, the filtered velocity gradient tensor ing from an initial condition at time t′) may be written
Aij = ∂ui/∂xj must remain trace-free, i.e. Aii = 0, and formally as follows
the filtered pressurefield is the solutionofthe respective
Poisson equation 2p= A A ∂2τ /∂x ∂x . t
∇ − lk kl− ij i j τ(t)=H(t,t′)τ(t′)H⊤(t,t′)+ H(t,s)Φ(s)H⊤(t,s)ds ,
Next, we also consider the transport equation for the
SGS stress tensor τ [13, 14, 15], which follows from Eq. Zt′ (5)
(3): where
Dτ = ∂τ +(u ∇)τ = τA⊤ Aτ +Φ , (4) DH(t,t′) = A(t)H(t,t′) and H(t′,t′)=I . (6)
Dt ∂t · − − Dt −
wherethetermΦ Φp+Φν ∇ J includesthepressure Toourknowledge,thisapproachtosolvethestressequa-
≡ − ·
gradient-velocity correlation tion in RANS closures was first suggested in the turbu-
lence literature by [19] (see equation (4.4) in this ref-
Φp,ij = ui∂jp ui∂jp+uj∂ip uj∂ip , erence). For the general case of time-varying veloc-
− − −
ity gradient, we note that the auxiliary matrix H(t,t′)
the viscous term(cid:2), (cid:3)
can be written as a time-ordered exponential (see Refs.
[18,20,21]forbackgroundonthisbasicmatrixfunction)
Φ =ν(u ∇2u u ∇2u +u ∇2u u ∇2u ) ,
ν,ij i j i j j i j i
− −
t
and the generalized central third-order moment H(t,t′)=Texp+ A(s)ds .
(cid:20)−Zt′ (cid:21)
=u u u u τ u τ u τ u u u .
ijk i j k j ik i jk k ij i j k
J − − − − Equation(5)illustratesclearlythedistinctrolesplayed
In II it is shown that a formal solution for the stress bytheproductiontermandthecontributiongivenbyΦ.
transp§ort equation may be obtained by integrating the Evaluation of Eq. 5 requires the knowledge of the time
production term exactly. This solution, suggested by historyofA(s)aswellasaccurateclosuresforΦ(s)along
′
[19]but—toourknowledge—littlepursued,willbeshown the fluid history t <s<t.
to involve matrix exponentials. The developments pre- As a next step, one may consider the special case
sented require some assumptions of Lagrangian isotropy in which the velocity gradient is considered to be con-
′
anddecorrelation,and some empiricalevidence support- stant between the initial time t and t, and set equal to
ing these assumptions is provided in III based on re- (e.g.) A(t) and simply denoted by A. For this approx-
sults from Direct Numerical Simulation§s (DNS). In IV imate situation, the solution of Eq.(6) may be written
the RFD closure for the SGS stress is developed. T§he as an ordinary matrix exponential H(t,t′) = e−(t−t′)A,
3
′
where the matrixexponentialis definedin the usualway and t t will be chosen to be of the order of such a
eB = +∞ Bn/n!. To simplify further, consider Eq. 5 decorr−elation time-scale. Clearly, one must also assume
n=0
forthecaseΦ=0,i.e. nowretainingonlytheproduction local isotropy to hold for the statistics, and this is jus-
P
term. This step eliminates the important isotropization tified from the usual arguments in turbulence when ∆
effects of pressure-strainand also the nonlinear diffusion is sufficiently small compared to the integral scale. In-
effects of the transport terms. While clearly missing im- cidentally, it is expected that a decorrelation between
portantphysics,it is still instructive to observethat this τ(t′) and A(t) may occur due to pressure effects, tur-
simplification allows the solution (Eq. 5) to be written bulent diffusion, etc. Some numerical evidence for such
as: decorrelation and isotropization is provided in the next
section.
τ(t)= e−(t−t′)A τ(t′) e−(t−t′)A⊤ . (7) ThetraceoftheconditionalSGStensor, τ (t′)A(t) ,
kk
h | i
must still be specified. The simplest option that is con-
At this stage it is conceptually advantageous to make
sistent with a local evaluation of velocity and length-
connectionwithRefs. [22,23,24],whereitisproposedto scalesistochooseafactorproportionalto∆2 S2,where
use conditional statistics to capture the relevant statis- ⊤ | |
S (A+A )/2 is the filtered strain rate tensor and
tics of the SGS stress. For example, in Ref. [23], it is S≡ (2S S )1/2. Finally, replacing into Eq. 8 with
shown that the least-square-error best estimate for the | | ≡′ ij ij
t t =τ , we get
SGSstressis ofthe formofamulti-pointconditionalav- − a
erage,namely τ u ,u ,...,u . Themulti-pointcon-
ditioning variahbliejs|{u11,u22,...,uNNi} are,in principle,con- τ(o) =cexp∆2|S|2e−τaAe−τaA⊤ , (10)
stitutedby theentire(N-point)resolvedvelocityfieldat
where the parameter c is unknown and may be ob-
scale ∆. To simplify the conditioning, one may limit the exp
tained by empirical knowledge, or by generalizing the
information to the past time-history of the local veloc-
dynamic model [25].
ity structure. In particular, a good choice that captures
For completeness and clarity, we remark that the ma-
much of the local dynamics in a Galilean invariant fash-
trix exponential solution may equivalently be obtained
ion is the Lagrangianpasthistory of the filtered velocity
gradient tensor A. The dependence on the Lagrangian by solving the linearized equation for a turbulent fluc-
tuation that only keeps the Rapid Distortion term from
timehistoryalongafluidparticleadvectedbythefiltered
the large-scale velocity field, and neglects all other ef-
resolvedvelocityfieldisthusassumedtobe describedby
A(s) with s t (here and below, the dependence of A fects. That is to say, we solve formally the equation
≤ D u′ = u′A using the matrix exponential function.
on spatial position is omitted for clarity). According to t i − k ik
The solution is then multiplied by its transpose to form
these ideas, we define a quasi-optimal SGS stress ten-
′ ′
sor τ(o)(t) as the conditional averageτ(o)(t)= τ(t)A . ui(t)uj(t)whichisthenaveragedoverthefluctuatingini-
Noticing that the matrix exponential prefactorshente|rinig tialconditionu′(t′)u′(t′)(conditioned ona constantA).
i j
′ ′ ′ ′
in Eq. 7 are some deterministic functions of the velocity The averaging of the term u (t)u (t) yields the initial
i j
gradient tensor itself, they thus can be taken out from (‘upstream’)stresstensorτ(t′),andwiththeconditional
this conditional average and we get the following stress averaging,anexpressionequivalentto Eq. 8 is obtained.
tensor: This is similar to the equivalence between solving the
equationforco-variancesorforthefluctuationsandthen
τ(o)(t)= e−(t−t′)A τ(t′)A e−(t−t′)A⊤ . (8) averaging, as noted in the context of stochastic linear
h | i
systems in [18].
With this expression, the closure problem has been Equation (10) represents a closure for the SGS stress
changed from requiring a model for the local stress ten- expressed in terms of matrix exponentials instead of the
sor at time t to requiring a model for the conditional more commonly used algebraic closures [26]. In section
′
averageofthe ‘upstreaminitial condition’ attime t <t. IV,aconnectionisnotedbetweentheexpressionEq. (10)
The initial condition needed is a symmetric tensor. In andaphysicalclosureforthe subgridstresstensorbased
the absence of additional information, the simplest as- ontherecentfluiddeformationclosureintheLagrangian
sumptionis topostulatethatthis conditionallyaveraged frame.
‘upstream’ stress tensor is isotropic, namely
1
τ (t′)A τ (t′)A δ . (9) III. EMPIRICAL EVIDENCE OF LAGRANGIAN
ij kk ij
h | i ≈ 3h | i DECORRELATION AND ISOTROPY FROM DNS
The magnitude of the tensor is proportional to the
trace of the SGS tensor and has units of squared ve- In order to verify whether the decorrelation and
locity. The assumption of isotropy may be justified if isotropizationofconditionalaveragesofSGSstressesoc-
τ(t′) and A(t) become more and more de-correlated as cur in turbulence, we analyze a DNS dataset of forced
′
the elapsed time t t grows, then no locally strong and isotropic turbulence. The simulation is conducted using
statistically preferr−ed direction should exist. This step a pseudo-spectral method in a [0,2π)3 box. 1283 grid
introduces a characteristic decorrelation time-scale τ , points are used. Fourier modes in shells with k <2 are
a
| |
4
forced by a term added to the Navier-Stokes equations
0.15
which provides constant energy injection rate ǫ = 0.1.
f
The viscosityofthe fluidisν =0.0032. Datais collected
afterthesimulationreachesstatisticalsteadystate. Note 0.1
that τ(t′) is the SGS stress at a previous time t′ and at
the spatial location occupied by the fluid particle which )
is at the position x at time t (i.e. X(t′;x,t), in the no- (t-t’0.05
tations of IV). According to the transport equation for I
τ(t) (Eq. 4§), the fluidparticle isadvectedbythe filtered
velocityfield. Thus,thepositionofthefluidparticleatt′ 0
isfoundbybackwardparticletrackingstartingfromend-
time t inthe filteredvelocityfield. Toperformbackward
-0.05
particle tracking, the filtered velocity and SGS stress 0 1 2 3 4_ 5 6 7 8
fields are calculated and stored at every ∆t = 0.009, (t-t’)A
rms
corresponding to 1/20 of the Kolmogorov time scale. A
Gaussian filter is used with filter scale ∆ = 15η, where FIG. 1: Decay of the anisotropy factor I (Eq. 11) as func-
η is the Kolmogorov length scale. In order to quantify tion of normalized time lag (t−t′)Arms measured in DNS of
′
isotropy as function of t t, the ratio of off-diagonal to isotropic turbulence.
−
on-diagonal tensor elements of the conditional averaged
′
SGS stress at decreasing previous time t is computed.
According to the derivation, the averaging must be 0.2
conditioned on a particular value of the resolved veloc-
ity gradient A(t). There are a large number of pos-
sibilities, since A(t) has 8 independent elements. As )]0.15
(t
representative of an important class of velocity gradi- |A12
ent structure, we choose to consider regions where the ’), 0.1
A(t) is such that it has a large shear in one direction, -t
(t
whereas all other velocity gradient tensor elements are τ12
-
weak. We choose a particular shear direction, “12”, and ρ[0.05
define E (t) to be a “high-12-shear” events that oc-
12
cur at time t. These events are defined here as those
points where 1) A12(t) > Arms, i.e. large and posi- 0
tive 12-shear, 2) A (t) < A for other off-diagonal 0 1 2 3 4_ 5 6 7 8
| ij | rms (t-t’)A
components (i,j) = (1,3) and (i,j) = (2,3), and 3) rms
A (t) <A /√2 for the diagonalelements i=j. The
ij rms
| | FIG. 2: Decay of the correlation coefficient ρ (see Eq. 12)
gradientrmsArms isdefinedasArms ≡hAijAiji1/2. This betweenτ12(t−t′)and−A12(t)asfunctionofnormalizedtime
definition allows a sufficiently large number of events to lag (t−t′)Arms measured in DNSof isotropic turbulence.
be counted and thus help in reaching statistical conver-
gence. With this definition of a conditioning event, we
′
calculate the isotropy factor I(t t) according to:
− from the stored fields using 6th order Lagrangian inter-
τ (t′)E (t) polation. The conditional averages are then found by
′ 12 12
I(t t) h | i . (11) averaging over all tracked particles. Statistical sampling
− ≡−(1/3) τ (t′)E (t)
h kk | 12 i is increased by averaging over the trajectories starting
fromseveraldifferentendtimes t andalsooverthe other
′
I(t t) monitors the isotropization of the SGS stress ′
asso−ciated with large “12 shear events”, i.e. a particular two 13 and 23 off-diagonal elements (in both τij(t) and
E (t)).
anisotropic condition in the large-scale velocity gradient ij
′
The resulting ratio I(t t) is plotted in Fig. 1 as a
tensor. Sincetheturbulenceisstatisticallyisotropic,sim- − ′
function of the normalized time lag (t t)A . It is
ilarresultsareexpectediftheothertwoshearcomponent − rms
evident that the conditional average of the SGS stresses
ofτ ,namelyτ andτ ,hadbeenchoseninsteadofτ ,
ij 13 23 12
under conditioning based on events E13(t) or E23(t), re- b′ecomesmoreisotropicasthetimelagincreasesandI(t−
t)crosseszeroatabout0.7eddyturn-overtimes,namely
spectively.
−1
Backward particle tracking starts from spatial loca- t∼0.7Arms. Thenthereisnegativeundershoottoabout
tionswheretheconditionsinE (t)areverified,attimet. negative half of the initial value, before it is relaxed to
′ 12 −1
At each time t <t, the particle locations are calculated around zero (the isotropic value) at about t 6A .
∼ rms
from the stored filtered velocity fields using a second- The undershoot below zero is an interesting trend and
orderAdam-Bashforthscheme. The filteredvelocityand understanding the physics of this behavior would be an
SGS stresses at the particle locations are interpolated interesting goal for future studies.
5
As additional evidence for the Lagrangiantime decor- a Taylor expansion of uδ and evaluation of the filtering
relation between stress and large-scale velocity gradient, operation at scale ∆ in Eq. 2 leads to
in Fig. 2 we show the correlation coefficient between
τ (t t′) and A (t), namely ∂uδ ∂uδ ∆
12 − − 12 τij C2∆2 i j, δ = . (13)
≈ ∂x ∂x β
′ k k
τ (t t)A (t)
12 12
ρ= h − i . (12)
− τ (t t′)2 A (t)2 One observes that similarity-type models such as the
h 12 − ih 12 i standardnonlinearmodel[2, 32]correspondtousing the
q gradientofthelarge-scalevelocityfield(δ =∆orβ =1).
The correlation is near 15% at zero time-lag (similar
Nevertheless, it is the case β >1 which is physically rel-
to the correlation coefficient between SGS stresses and
evant since the true SGS stress includes scales smaller
strain-rate tensor often quoted in a-priori studies), but
−1 than ∆. However, for β > 1, the expression 13 does
then decays to nearly zero at times around t ∼ 2Arms. not constitute a closure since then uδ contains sub-grid
Taken together,the DNS analyses thus provide evidence
motions that are not known at the LES filter scale ∆.
for the isotropy assumption on τ (t′)A(t) , as long as
t−t′ &τa with τa =A−rm1s. h ij | i froAms ihner[e5]oann)d, aChLeavgilrlaanrdgia&nMlaebneelvepaousit(i2o0n06X–iCsMem06-
Note that due to the cost of storing the entire sim- ployedtoencodethetime-historyinformation. Usingthe
ulation for backward particle tracking, only moderate two-timeformulationof[3],the labelpositionsX(t′;x,t)
Reynolds numbers were considered in the analysis. The satisfy dX/dt′ = u(X(t′),t′) with X(t) = x. Thus
forcing length scale has been estimated to be about 50 X(t′;x,t) represents the position X at a prior time t′ of
times the Kolmogorov length scale, η, and the viscous thefluidparticlewhichisatpositionxattimet.Making
effects begin to significantly damp the motions at scales the Eulerian-Lagrangianchange of variables also used in
of about 10η and smaller. Therefore, using ∆ = 15η, CM06 leads to the following expression:
there may be some effects from the forcing and viscous
scales on the results. However, the observed tendency
towards isotropization is expected to become more, not τ =C ∆2 ∂Xp∂Xq ∂uδi ∂uδj . (14)
less,prevalentathigherReynoldsnumbers. Wepointout ij 2 ∂x ∂x ∂X ∂X
k k p q
that opportunities for much more in-depth future anal-
yses of such issues are provided by the availability of a Alltermsinthisexpressionarestronglyfluctuatingvari-
turbulence databaseat higher Reynolds number [29] (al- ables. But,asinpriorsection,themostrelevantinforma-
though this database could not be used for the present tion is retained by the conditional averaged expression.
data analysisdue to the fact thatit does notyetcontain We proposethe same conditionalaveragingbasedonthe
sufficiently efficient means of filtering the data). time-historyofthevelocitygradienttensoralongthepast
fluid particle trajectory. Therefore, combining the con-
ditional averaging and the change of variables one may
IV. STRESS TENSOR MODEL BASED ON THE write
RECENT FLUID DEFORMATION CLOSURE
∂X ∂X ∂uδ ∂uδ
τ(o)(t)=C ∆2 p q i j A(s);t′ <s t ,
ThisalternativeapproachisbasedonrelatingtheSGS ij 2 *∂xk ∂xk ∂Xp∂Xq | ≤ +
stress tensor to small-scale velocity gradients. To begin, (15)
one may recall the multiscale expansion [27, 28, 30] in wherethedependenceofstressτ(o)(t)oncurrentposition
which, among others, the exact subgrid stress (Eq. 2) is x is understood and not indicated to simplify the nota-
writtenintermsofuδ,thevelocityfieldcoarse-grainedat tion. The Jacobian matrix G (t′,t) = ∂X (t′;x,t)/∂x
ij i j
ascaleδ,butstillwithδ <<∆,i.e. containingsignificant satisfies (see for instance [17])
contributions from sub-grid scales. One may then define
the approximated stress tensor τiδj = uδiuδj −uδi uδj and DtG(t′,t)=−G(t′,t)A(t) with G(t,t)=I , (16)
naturally τ =limδ→0τδ. where I is the identity matrix. Thus,
Consistent with the Kolmogorov phenomenology, as
argued formally in [31], and also as used in various a-
t
priori analyses of experimental data (see e.g. [28]) the G(t′,t)=Texp− A(s)ds
SGS stress is relatively local in scale, stating that the (cid:20)−Zt′ (cid:21)
leading terms entering in its development are given by
is expressed as an “anti-time-orderedexponential”, with
the coarse-grainedvelocity at the resolutionscale ∆ and
matrices ordered from left to right for increasing times
including also the next range of length-scales between
[20,21,33]. Theonlydifferencewiththematrixfunction
δ ≈∆/β and∆(e.g. β ∼2). Asaconsequence,onemay H(t′,t) of the preceding section (Eq. (6)) is the sense of
use the approximation τij ≈ τiδj=∆/β. Furthermore as- time-ordering.
suming thatuδ=∆/β is sufficientlysmoothoverdistances SincethedeformationgradienttensorG =∂X /∂x
pk p k
∆(orusingthe‘coherentsubregionapproximation’[31]), is a deterministic function of the past velocity gradient
6
history,thesetensorscanbetakenoutsidetheconditional wherethe parameterc =C c, ina similarfashionas
exp 2
·
averagesin Eq. 15. So far one can thus write intheprecedingsection,isunknownandmaybeobtained
by empirical knowledge, or by generalizing the dynamic
∂X ∂X
τ(o)(t)=C ∆2 p q Y model [25].
ij 2 ∂x ∂x ijpq
k k As a final step, the Recent Fluid Deformation (RFD)
approximationisused(CM06)inwhichthetime-varying
with Y = ∂uδi ∂uδj A(s);t′ <s t velocity gradient A(s) between t′ and t is approximated
ijpq
*∂Xp∂Xq | ≤ + with a constant value (e.g.) equal to its value at t
and denoted by A. The initial condition for the fluid
where Y is a 4th rank Lagrangian gradient tensor. At
deformation (when the deformation gradient tensor is
this stage, it is now possible again to invoke Lagrangian
assumed to be the identity), is prescribed at the time
isotropy, following the approach of CM06 and of the ′
t < t. The solution to Eq. 16 can then be written as
preceding section. It is assumed that the tensor Y is G(t′,t) = e−(t−t′)A. Note that in this approximation,
isotropic due to loss of information caused by turbulent H(t,t′)=G(t′,t), since the sense of ordering of the ma-
′
dispersion, past pressure effects, etc. if t t is long
− trix products is no longer significant.
enough. Under the Lagrangian-isotropyclosure assump- ThenextstepistoreplacethesolutionforG(t′,t)into
tion, one may write
Eq. 18. And, as was done in the preceding section, to
assume that a characteristic de-correlation time-scale τ
a
Y =A′δ δ +B′δ δ +C′δ δ . (17) has elapsed between the time where the initial isotropy
ijpq ij pq ip jq iq jp
assumption is justifiable and the current time when the
While individual realizations of a small-scale gradient
stressclosureisrequired. Thismeansreplacingtheinitial
tensor in turbulence are of course not isotropic, statis- ′
time t witht τ . Finally,the closureforthe deviatoric
a
tical moments such as the conditional average can be −
part of the stress reads
morejustifiably approximatedas isotropic. The isotropy
assumptionstatesthattherateofchangeofturbulentve-
locitiesδuδ(x,t)(atthepresentlocationandtime(x,t)), τ(od) =c ∆2 S2 e−τaAe−τaA⊤ d (19)
withrespecttochangesinpastlocationsofthefluidpar- ij exp | |
ticles at time t′, is insensitive to orientation of δX. This h i
where all quantities are evaluated at (x,t). It is imme-
appears to be a plausible postulate, if sufficient time has
diately apparent that this closure is equivalent to the
′
elapsed, i.e. if t t is sufficiently large for decorrela-
− formal solution developed in the previous chapter (see
tionto takeplace. In the precedingsectionweuseddata
Eq. (10)).
from DNS to test the accuracy of such a de-correlation
and‘Lagrangianisotropyassumption’inacloselyrelated
context(directlybasedonthestressesratherthansmall-
V. EXPANSIONS
scalevelocitygradientstatistics). Still,itisimportantto
recognizethat this step is introducedhere as an‘ad-hoc’
Inprecedingsectionsithasbeenshownthatamatrix-
closure assumption and no claim is made that this is a
exponential closure for the deviatoric part of the SGS
formal step with controlled errors.
stress tensor may be written as in Eq. 19. As a next
While the assumption of isotropy eliminates the de-
step, the behavior of this closure is explored when τ is
pendence of Y upon A, the latter still affects the Ja- a
ijpq small enough so that the norm of τ A is much smaller
cobian matrix G = ∂X /∂x that enters in the closure a
for the SGS streijss. Nexit, wje focus attention only on than unity. Then e−τaA I τaA+(1/2)(τaA)2+....
≈ −
Up to second order one then obtains
the trace-free part of the modeled SGS stress tensor, i.e.
we subtract the trace of the stress. And, noticing that
tohnely‘rtihghetu’nCkanuocwhyn-BGr′e+enCt′eennstoerr∂skiXnpt∂hkeXreqsiusltsiynmgm‘qeutarsici-, τod ≈cexp∆2|S|2 −2 τa S
optimal’ model for the deviatoric part of the stress (su-
d
perscript od) model τ(od). Dimensionally, the parameter + τ2 A A⊤+ 1 A2+(A⊤)2 +... , (20)
B′+C′ has units of iinjverse time-scale squared, and de- a (cid:20) 2(cid:16) (cid:17)(cid:21) !
pends upon the turbulence statistics down to scales δ.
Crow,inRef. [19],derivedessentiallythesameresult(see
For a fixed ratio ∆/δ, and with both scales in the iner-
hiseq.(5.3))butwithunspecifiedcoefficientsobtainedas
tial range, for simplicity we assume that the parameter
moments of his memory kernel and an additional term
′ ′
B + C follows, as in the prior section, ‘Smagorinsky
proportional to the material derivative D S. It is imme-
scaling’, i.e. B′+C′ c S2. One thus obtains t
≈ | | diately apparent that if the time-scale τa is chosen as
τ = S−1,thenthefirsttermisthestandardSmagorin-
a
τ(od)(t)=c ∆2 S2 ∂Xi∂Xj 1∂Xm∂Xmδ , sky m|o|del with cexp = c2s (where cs is the Smagorinsky
ij exp | | ∂x ∂x − 3 ∂x ∂x ij coefficient). Furthermore, the second term, the term in
(cid:18) k k k k (cid:19)
(18) the square parentheses, is of the form of the ‘nonlinear
7
model’ [2, 28, 32] with a prefactor cexp∆2. Two differ- 100
ences with the standard ‘nonlinear model’ are apparent,
however. The firstis thatif c c2, then as coefficient
exp ∼ s 10-1
ofthenonlineartermthisissignificantlysmallerthanthe
coefficient for this term normally mentioned in the liter-
ature (whichrangestypically between1/12to 1/3). The 10-2
)
k
second difference is the presence of the additional term (
A2+(A⊤)2 /2. To make connections with standard E10-3
(cid:16)non-linear mo(cid:17)dels used more often in RANS (e.g. [34]),
the velocity gradient is decomposed into symmetric and 10-4
antisymmetric parts, A = S+Ω. The result is (again
with τ = S−1) 10-5
a | | 100 101 102
k
d
τod 2c2∆2 SS+c2∆2 S2+ 1 Ω S S Ω .
≈− s | | s 2 −
(cid:20) (cid:21) FIG. 3: Radial kinetic energy spectra of forced isotropic tur-
(cid:0) (cid:1) (21)
bulencefromLESusingthematrix-exponentialclosureofEq.
Itterims in(Ate2re+st(inAg⊤t)o2)n/o2tecatnhcaetlstheexaecxtplyantshieonΩinΩclupdairntgththaet 1γ9=wi2th, acnexdpd=as(h0e.d1)2li.neS:olγid=lin0e.:5.γ D=ot1t,eddalsihn-ed:otutnedivelirnsae:l
is included in the standardnon-linear model A A⊤. For Kolmogorov spectrum E(k)=1.6ǫ2f/3k−5/3.
detailed a-prioristudies of the various decompositionsof
the velocity gradient and non-linear terms see [35].
0
VI. MATRIX EXPONENTIAL CLOSURE IN -0.2
LES OF ISOTROPIC TURBULENCE
S -0.4
The expansion introduced in the last section is for-
mallyvalidonlyforsmallvaluesofthenormofτ A. For
a
morerealisticlargervalues,theexpansionmaybeinaccu- -0.6
rateandmanyadditionalhigher-orderterms areneeded.
Theycanallbeexpressedinterms ofexpansionsintoin-
-0.8
tegrity bases [26], but it is in general difficult to obtain 0 0.5 1 1.5 2 2.5 3
τ
t/
the coefficients of the expansion. Instead, it is proposed
L
here to utilize the matrix exponential directly in simu-
lations. Since the exponential involves the full velocity
FIG.4: LongitudinalderivativeskewnesscoefficientSasfunc-
gradient tensor, it appears more natural to choose the
time-scale τ according to τ = γ(A A )−1/2 γ A−1 tion of simulation time. Lines are the same as in Figure 3.
instead of uasing S−1. Theaparametijeriγj is an≡emp|iri|cal τL≈6 is theintegral time scale.
| |
coefficient of order unity.
As a first test, LES of forced isotropic turbulence is
pseudo-spectral scheme, the modeled SGS stress is eval-
performed. This flow is the simplest possible test-case
uated in physicalspace and is made trace-free (this only
and it is used here simply to determine whether simula-
affects the effective pressure, not the dynamics) before
tions using the matrix-exponential based closure are nu-
computing its divergence in Fourier space.
merically stable yielding realistic energy spectra, and to
Thematrixexponentialsareevaluatedusingtruncated
ascertain the associated computational cost. The gener-
Taylorexpansionwithscalingandsquaring[36]. Specifi-
alization to dynamic versions and tests in more complex
cally, we need to evaluate exp(B), where B= γA/A.
flowswillbeleftforfutureinvestigations. Thesimulation
− | |
For a matrix C in general, the Kth order truncated
usesthesamepseudo-spectralmethodaswasusedinthe
Taylor expansion uses matrix polynomial T (C) =
DNS outlined in II, with same grid resolution, forcing K
scheme and time§step size. Dealiasing is performed by K Cn/n! to approximate exp(C), incurring an error
n=0
zero-padding according to the two-third rule. The vis- bounded by C K+1 / [1 C /(K +2)](K +1)! .
cosity of the fluid is ν = 0.000137. The subgrid-scale TPhe error dekcreakses wit{h −thek norkm of the matrix C}.
model implemented is given by Eq. 19 and c =(0.1)2 Therefore,to evaluate exp(B), we firstdefine C=B/2j,
exp
is chosen (dynamic versions [25] of this model to deter- where the value of the integer j is chosen to ensure
mined cexp can be developed in the future). To specify k C k6 1/2. exp(C) is then approximated by TK(C)
τ , the values γ=0.5, 1 and 2 are tested (a dynamic ap- and finally exp(B) is given by [T (C)]2j. The cost of
a K
proachofdeterminingγ couldalsobedeveloped). Inthe calculating T (C) is reduced by using Cayley-Hamilton
K
8
theoremto expressCn (n>2)interms ofI, C, C2, and focusofattentionintheliterature. InLES,however,due
the invariants of C. Choosing K = 7 , we obtain the topracticalconstraints‘algebraic’closuresaremostoften
following equation for T (C) with an error smaller than preferred. Thepresentapproachshowsthattheeffectsof
7
10−8: production in the context of such closures may be taken
into account directly based on an exact solution of the
T7(C)=C0I+C1C+C2C2 (22) stress transport equation. A central step in the present
approachistouseisotropyforthe‘upstream’initialcon-
where
dition. Evidence for such isotropization of initial con-
R Q R Q2R dition, given present large-scale velocity gradients, has
C0 =1 C + C C C C, beenobtainedusingaDNSdatabase. Implementationof
− 3! 5! − 7!
the closure in LES of forced isotropic turbulence yielded
Q R Q2 2Q R R2 Q3
C =1 C C + C + C C + C − C, good results. The computational cost is significant, but
1
− 3! − 4! 5! 6! 7!
itisnotprohibitive. Sinceourcodewiththismodeltook
1 Q R Q2 2Q R
C = C C + C + C C . abouttwiceaslongtorunaswithatraditionalalgebraic
2
2 − 4! − 5! 6! 7! closure, LES with this model at a resolution of N3 has
similarCPUcostasLESwithatraditionalmodelrunat
Here Q = Tr(C2)/2 and R = Tr(C3)/3 are the
C − C − a resolution of (21/4N)3 (1.2N)3.
two non-zero invariants of C (note that Tr(C) = 0). ∼
In terms of cost, the above algorithm uses about (1 + It is crucial to stress that the additional, more subtle
j)N3 +5N2 +2N +37 flops to calculate exp(B) when physics of the remaining terms in Eq.4 (pressure effects,
B is given, where N is the dimension of the matrix. In turbulent diffusion, dissipation, etc..) are, in general,
our tests j = 1+floor(log γ), so j = 1 when γ = 1 and unlikely to be well-represented by the simple assump-
2
the cost is estimated at about 140 flops for each stress tion of ‘upstream’ isotropy. In addition, non-equilibrium
evaluation. Thiscanbecomparedwiththe singlematrix conditions in which A varies quickly along the particle
multiplication needed for the nonlinear model, which is trajectory are not included in the closure as written in
about N3 30 flops. Overallwith this closure,our code Eq. 19, in which the velocity gradient is assumed to
tookabout∼twiceaslongtorunascomparedtousingthe have remained constant over a time-scale τa. To explore
mixed model. non-equilibriumeffects,thefulltime-orderedexponential
Simulations were initialized with random Fourier functionmustbeused,althoughthiswouldstillleaveout
modes and evolved until statistical steady state was ob- thenon-equilibriumeffectsofΦ. Tocomparethepresent
tained. No numerical instabilities were observed for the approach to other closures will require more in-depth
three parameter cases considered (c = 0.01, γ =0.5, testinginmoredemanding,complexflows(e.g. whereef-
exp
1 and 2). In Figure 3 the energy spectra obtained from fects of anisotropy, non-equilibrium, and pressure-strain
thethreesimulationsasaveragedinthetimeintervalbe- correlations are expected to be important).
tweenoneandthreelarge-eddyturnovertimesareshown. It is also instructive to consider the case of two-
Figure4showsthetime-evolutionofthederivativeskew- dimensional (2D) turbulence. Nothing in the closure
ness coefficient = (∂1u1)3 / (∂1u1)2 3/2. As can be strategies pursued here limits their application to space
S h i h i
seen,thecaseγ =1appearstoyieldphysicallymeaning- dimension three, at least nothing very obvious. How-
ful results, which can be compared with the well-known ever, the expansions (Eqs. 20,21) show that this is not
results of the Smagorinsky model and the mixed model likely to be a qualitatively goodclosurefor space dimen-
(see,forexample,[40]). But,thereiscleardependenceon sion two, since one there expects an effective “negative
the parameter γ. The skewness coefficient quickly drops eddy-viscosity” corresponding to inverse energy cascade
tovaluesnear 0.3forγ =1and 0.36forγ =2. These [39]. Itisthusworthreflectingonsomeofthereasonsfor
− −
arerealisticvaluesforfilteredturbulence[37]. Theskew- the inaccuracy of the closures in 2D, since this may help
ness values for γ = 0.5, on the other hand, appear to pinpoint potential shortcomings in 3D as well. First, it
be too close to zero, consistent with some pile-up of the is known that the 2D inverse cascade is less local than
spectrum at high wave-numbers. the 3D forward cascade, with most of the flux coming
from triadic interactions for a scale-ratio β = 4 8.
∼
[38, 39]. However, the starting point of the RFD clo-
VII. DISCUSSION sure, Eq. 13, is not accurate for β 1. To get a qual-
≫
itatively reasonable alternative at β substantially larger
A new closure based on matrix exponentials and as- than 1—which involves only first-order gradients—one
sumptions about short-time Lagrangian dynamical evo- must instead use something like the “Coherent Subre-
lution has been proposed. Matrix exponentials as for- gions Approximation” of [39]. On the other hand, the
mal solution of the stress transport equation provides startingpointofthe closureapproximationinSectionII,
interesting insights into the effects of the production the stress transport equation 4, is exact in 2D just as in
(gradient-stretching) term. Historically, in the context 3D.Thefailureoftheclosureprocedurein2Disnowdue,
ofRANS modeling using additionaltransportequations, presumably, to the effects of the Φ source-terms in the
the (closed)productiontermhasjustifiably notbeenthe transport equation. Indeed, those terms are expected to
9
contribute as an effective “negative viscosity”, primarily is expressed in terms of subgrid-scale velocity gradients.
due to the pressure-Hessian rotating small-scale strain Perhapsitcanbeexpectedthatcastingthisnewlighton
matricesrelativetothelarge-scalestrain[39]. Notethat, the closure problem improves our understanding of this
strictlyspeaking,thisisprobablyalsotruein3D,sothat long-standing problem.
thematrix-exponentialclosuresarelikelytobeoverlydis- Finally, we remark that many transport equations for
sipative in every dimension. The main effect of the gra-
turbulencemomentshaveabasicstructuresimilartoEq.
dient stretching terms—which is a tendency to forward 4, including two production terms involving the veloc-
cascade, or positive eddy-viscosity—is well captured by
ity gradient and its transpose. Examples include higher-
thematrix-exponentialclosureinanydimension,butthe order moments of velocity, the spectral tensor encoun-
additional,moresubtlephysicsoftheremainingtermsin
teredinRapidDistortionTheorycalculations,etc... The
Eq. 4 are most likely not well-represented by the simple formal solution in terms of matrix exponentials provides
assumption of isotropy.
new possibilities of calculation and insights into the un-
As has been cautioned severaltimes in this paper, the derlying physics.
simplified matrix-exponential closure as written in Eq.
10 employs the drastic approximation of entirely omit-
ting the pressure-strain correlation and other ‘nonlin-
ear scrambling’ terms. But unlike eddy-viscosity based
Acknowledgments
closure assumptions, this expression can be derived di-
rectly from a relevant fluid dynamical equation, namely
the stress transport equation (with only the production WethankS.B.PopeandR.Rubinsteinforusefulcom-
term), and using physically motivated and straightfor- ments. L.C. thanks the Keck Foundation for financial
ward assumptions about Lagrangian decorrelation and support. The other authors thank the National Science
upstream isotropy. A similar result is obtained using an Foundationforfinancialsupport,andC.M.alsoacknowl-
Eulerian-Lagrangianchange of variables when the stress edges partial support from the Office of Naval Research.
[1] M. Lesieur, M and O. M´etais. New trends in large-eddy eling thepressureHessian and viscousLaplacian in Tur-
simulations of turbulence.Ann. Rev. Fluid Mech. 28, 45 bulence: comparisons with DNS and implications on
(1996). velocity gradients dynamics, Phys. Fluids 20, 101504
[2] C. Meneveau and J. Katz, Scale-Invariance and Turbu- (2008).
lenceModelsforLarge-EddySimulation,Ann.Rev.Fluid [13] A.Leonard, Energy cascadein large eddysimulations of
Mech. 32, 1 (2000). turbulent fluid flows, Adv. Geophys. 18, 237 (1974).
[3] R.H. Kraichnan, Lagrangian-History Closure Approxi- [14] J.W. Deardorff, Three-dimensional numerical study of
mation for Turbulence, Phys. Fluids 8, 575 (1965). the height and mean structure of a heated planetary
[4] H.TennekesandJ.L.LumleyAfirstcourseinturbulence, boundary layer, Bound. Layer Meteor. 7, 81 (1974).
MIT press, Cambridge (1972). [15] M.Germano,Turbulence: thefilteringapproach,J.Fluid
[5] S.S. Girimaji, Modeling Turbulent Scalar Mixing as En- Mech. 238, 325 (1992).
hanced Diffusion, Combust. Sci. Tech. 97, 85 (1994). [16] C. Penland, Random forcing and forecasting using prin-
[6] M. Chertkov, A. Pumir and B.I. Shraiman, Lagrangian cipal oscillation pattern analysis, Mon. Wea. Rev. 117,
tetrad dynamics and the phenomenology of turbulence 2165 (1989).
Phys. Fluids 11, 2394 (1999). [17] C.TruesdellandW.Noll,TheNon-LinearFieldTheories
[7] Y. Li and C. Meneveau, Intermittency trends and La- of Mechanics, Springer-Verlag, Berlin, 1992.
grangian evolution of non-Gaussian statistics in turbu- [18] G.L. Eyink, Linear stochastic models of nonlinear dy-
lent flow and scalar transport, J.Fluid Mech. 558, 133 namical systems, Phys. Rev. E. 58, 6975 (1998).
(2006). [19] S.C.Crow,Viscoelastic propertiesoffine-grainedincom-
[8] E.JeongandS.S.Girimaji,Velocity-GradientDynamics pressible turbulence,J. Fluid Mech. 33, 1 (1968).
in Turbulence: Effect of Viscosity and Forcing, Theor. [20] J.D.DollardandC.N.Friedman,ProductIntegration,in
Comput. Fluid Dyn.16, 421 (2003). Encyclopedia of Mathematics, Addison-Wesley, London
[9] A. Pumir and B. Shraiman, Lagrangian Particle Ap- (1979).
proachtoLargeEddySimulationsofHydrodynamicTur- [21] W.J.Rugh,LinearSystemTheory,2ndEd.PrenticeHall
bulence,J. Stat.Phys. 113, 693 (2003). , New Jersey, 1996.
[10] L. Chevillard and C. Meneveau, Lagrangian dynamics [22] R.J. Adrian, Stochastic estimation of sub-grid scale mo-
and statistical geometric structure of turbulence, Phys. tions, Appl. Mech. Rev. 43, 214 (1990).
Rev. Lett. 97, 174501 (2006). [23] J.LangfordandR.Moser,OptimalLESformulationsfor
[11] L. Chevillard and C. Meneveau, Intermittency and uni- isotropic turbulence, J. Fluid Mech. 398, 321 (1999).
versality in a Lagrangian model of velocity gradients in [24] S.B. Pope, Turbulent Flows, Cambridge Univ.
three-dimensionalturbulence,C.R.M´ecanique 335,187 Press,Cambridge, 2000.
(2007). [25] M. Germano, U. Piomelli, P. Moin and W.H. Cabot, A
[12] L.Chevillard,C.Meneveau,L.Biferale, F.Toschi,Mod- dynamicsubgrid-scaleeddyviscositymodel,Phys.Fluids
10
A, 3, 1760 (1991). of Reynolds-Stress Closures in Turbulence, Ann. Rev.
[26] S.B.Pope,Amoregeneraleffective-viscosityhypothesis, Fluid Mech. 23, 107 (1991).
J. Fluid Mech. 72, 331 (1975). [35] K. Horiuti, Roles of non-aligned eigenvectors of strain-
[27] R.H. Kraichnan, On Kolmogorov’s inertial-range theo- rateandsubgrid-scalestresstensorsinturbulencegener-
ries, J.Fluid Mech. 62, 305 (1974). ation J. Fluid Mech. 491, 65 (2003).
[28] S. Liu, C. Meneveau and J. Katz, On the properties of [36] C. Moler and C.F. Van Loan, Nineteen dubious ways to
similaritysubgrid-scalemodelsasdeducedfrommeasure- compute the exponential of a matrix, twenty-five years
mentsin a turbulentjet, J. Fluid Mech. 275, 83 (1994). later, SIAM Rev. 45, 3 (2003).
[29] Y.Li,E.Perlman,M.Wan,Y.Yang,R.Burns,C.Men- [37] S. Cerutti and C. Meneveau, Statistics of filtered veloc-
eveau, R. Burns, S. Chen, A. Szalay and G. Eyink, A ity in grid and wake turbulence, Phys. Fluids 12, 1143
public turbulence database cluster and applications to (2000).
studyLagrangianevolutionofvelocityincrementsintur- [38] S.Chen,R.E.Ecke,G.L.Eyink,M.Rivera,M.Wan,and
bulence,J. Turbulence 9, 31 (2008). Z. Xiao, Physical Mechanism of the Two-Dimensional
[30] G.L. Eyink, Locality of turbulent cascades, Physica D Inverse Energy Cascade, Phys. Rev. Lett. 96, 084502
207, 91 (2005). (2006).
[31] G.L.Eyink,Multi-scalegradientexpansionoftheturbu- [39] G.L. Eyink, A turbulent constitutive law for the two-
lent stress tensor, J. Fluid Mech. 549, 159 (2006). dimensional inverseenergy cascade, J. Fluid Mech. 549,
[32] R.A. Clark, J.H. Ferziger and W.C. Reynolds, Evalua- 191 (2006).
tion of subgrid-scale models using an accurately simu- [40] H.-S.Kang,S.Chester,andC.Meneveau,Decayingtur-
lated turbulentflow, J. Fluid Mech. 91, 1 (1979). bulenceinanactive-grid-generatedflowandcomparisons
[33] G. Falkovich, K. Gaw¸edzki and M. Vergassola, Particles with large-eddy simulation, J. Fluid Mech., 480, 129
and fields in fluid turbulence Rev. Mod. Phys. 73, 913 (2003).
(2001).
[34] C.G. Speziale, Analytical Methods for the Development