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

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

## Abstract

Two approaches for closing the turbulence subgrid-scale stress tensor in terms of matrix exponentials are introduced and compared. The first approach is based on a formal solution of the stress transport equation in which the production terms can be integrated exactly in terms of matrix exponentials. This formal solution of the subgrid-scale stress transport equation is shown to be useful to explore special cases, such as the response to constant velocity gradient, but neglecting pressure-strain correlations and diffusion effects. The second approach is based on an Eulerian-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) approximation. It is shown that both approaches lead to the same basic closure in which the stress tensor is expressed as the product of the matrix exponential of the resolved velocity gradient tensor multiplied by its transpose. Short-time expansions of the matrix exponentials are shown to provide an eddy-viscosity term and particular quadratic terms, and thus allow a reinterpretation of traditional eddy-viscosity and nonlinear stress closures. The basic feasibility of the matrix-exponential closure is illustrated by implementing it successfully in Large Eddy Simulation of forced isotropic turbulence. The matrix-exponential closure employs the drastic approximation of entirely omitting the pressure-strain correlation and other ‘nonlinear scrambling’ terms. But unlike eddy-viscosity closures, the matrix exponential approach provides a simple and local closure that can be derived directly from the stress transport equation with the production term, and using physically motivated assumptions about Lagrangian decorrelation and upstream isotropy.

## I Introduction

One of the most basic challenges in turbulence modeling is the need for closures for the fluxes associated with unresolved turbulent fluctuations. In the context of Large Eddy Simulation (LES), closures are required for the subgrid-scale (SGS) stress tensor (1); (2). Traditional closures involve mostly algebraic expressions relating the stress tensor to powers of the velocity gradient tensor. More elaborate approaches using separate transport equations have sometimes also been employed, although these tend to be significantly more costly in the context of LES. Closures expressing the stress in terms of the matrix exponential function do not appear to have received much attention in the literature. The objective of the present work is to identify and discuss two separate paths that lead to such closures. Both paths are based on the Lagrangian dynamics of turbulence, i.e. on an understanding of the evolution of turbulence as one follows fluid-particle paths in time.

The use of Lagrangian concepts in turbulent flows has a long history (3); (4) and, in recent years, has seen renewed interest for modeling (5); (6); (8); (9); (7). Among others, a new model for the pressure-Hessian tensor based on the recent Lagrangian evolution of fluid elements – the Recent Fluid Deformation (RFD) closure – has been proposed (10); (11); (12). In this approach, a change of variables is made expressing spatial gradients in terms of Lagrangian gradients (e.g. how does a variable at the present location vary if we change the initial position of the fluid particle at an earlier time). Then the assumption of isotropy is introduced for the Lagrangian gradient tensors. This assumption allows simpler isotropic forms to be used, and is argued to be justified based on Lagrangian decorrelation effects. Deviations from isotropy at the present location for the Eulerian gradient tensors develop as a result of fluid material deformation along the Lagrangian trajectory. More traditionally, the Lagrangian time evolution of the stress tensor following fluid particles can be derived by taking appropriate moments of the Navier-Stokes equations. In this paper we examine both of these approaches to formulate models for the SGS stress tensor in turbulence in the context of LES.

A description of small-scale structure of turbulence begins with the Navier-Stokes equations of an incompressible fluid of velocity u:

 dudt=∂u∂t+(u⋅\boldmath∇)u=−\boldmath∇p+ν% \boldmath∇2u , (1)

where stands for the Lagrangian material derivative, the pressure divided by the density of the fluid and the kinematic viscosity. Because of incompressibility, the velocity gradient tensor must remain trace-free, i.e. , and the pressure field is the solution of the Poisson equation .

In the framework of LES, the SGS stress tensor is defined using the filtering approach (13); (15),

 τij=¯¯¯¯¯¯¯¯¯¯uiuj−¯¯¯ui¯¯¯uj . (2)

An overbar denotes spatial filtering at a scale and is formally given by a convolution with a nonnegative, spatially well-localized filtering function of characteristic size , with unit integral , namely

 ¯¯¯u(x,t)=∫G(r)u(x+{r},t)dr .

The SGS tensor enters in the dynamics of the filtered velocity as it can be seen when applying the filtering procedure to the Navier-Stokes equations (Eq. (1)),

 D¯¯¯uDt=∂¯¯¯u∂t+(¯¯¯u⋅\boldmath∇)¯¯¯u=−\boldmath∇¯¯¯p+ν\boldmath∇2¯¯¯u−\boldmath∇⋅\boldmathτ , (3)

where stands for the Lagrangian material derivative with as the advecting velocity, and the filtered pressure divided by the density of the fluid. Because of incompressibility, the filtered velocity gradient tensor must remain trace-free, i.e. , and the filtered pressure field is the solution of the respective Poisson equation .

Next, we also consider the transport equation for the SGS stress tensor (13); (14); (15), which follows from Eq. (3):

 D\boldmathτDt=∂\boldmathτ∂t+(¯¯¯u⋅\boldmath∇)% \boldmathτ=−\boldmathτ¯¯¯A⊤−¯¯¯A\boldmathτ+\boldmathΦ , (4)

where the term includes the pressure gradient-velocity correlation

 Φp,ij=−[¯¯¯¯¯¯¯¯¯¯¯¯ui∂jp−¯¯¯¯¯ui∂j¯¯¯p+¯¯¯¯¯¯¯¯¯¯¯¯uj∂ip−¯¯¯¯¯uj∂i¯¯¯p] ,

the viscous term,

 Φν,ij=ν(¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯ui\boldmath∇2uj−¯¯¯ui\boldmath∇2¯¯¯uj+¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯uj\boldmath∇2ui−¯¯¯uj\boldmath∇2¯¯¯ui) ,

and the generalized central third-order moment

 Jijk=¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯uiujuk−¯¯¯ujτik−¯¯¯uiτjk−¯¯¯ukτij−¯¯¯¯¯ui ¯¯¯¯¯uj¯¯¯uk .

In §II it is shown that a formal solution for the stress transport equation may be obtained by integrating the production term exactly. This solution, suggested by (19) but—to our knowledge—little pursued, will be shown to involve matrix exponentials. The developments presented require some assumptions of Lagrangian isotropy and decorrelation, and some empirical evidence supporting these assumptions is provided in §III based on results from Direct Numerical Simulations (DNS). In §IV the RFD closure for the SGS stress is developed. The resulting model is shown to be expressible compactly in terms of matrix exponentials as well. Differences and similarities between the RFD and transport equation solutions are discussed. In §V, the matrix-exponential solutions are expanded for short times. The expansions allow to establish relationships to traditional eddy-viscosity and nonlinear closure models in turbulence. In §VI the matrix-exponential closure is implemented in a most simple flow to illustrate its feasibility and cost.

## Ii Solution to stress transport equation using matrix exponentials

Equation 4 is of the form of the “time-dependent Lyapunov equation”, if the tensor ’s implicit dependencies upon the velocity fluctuations and the stress tensor were disregarded (in reality, depends upon small-scale velocity fluctuations and thus the full equation is highly non-linear and non-local). The formal solution of the Lyapunov equation in terms of matrix exponentials has been found useful in a number of other fields: principal oscillation pattern analysis (16), mechanics of finite deformations (17), and fluctuation-dissipation theorems for stochastic linear systems (18). In the context of the SGS stress transport equation the solution at time (starting from an initial condition at time ) may be written formally as follows

 \boldmathτ(t)=H(t,t′)\boldmathτ(t′)H⊤(t,t′)+∫tt′H(t,s)\boldmathΦ(s)H⊤(t,s)ds , (5)

where

 DH(t,t′)Dt=−¯¯¯A(t)H(t,t′) and H(t′,t′)=I . (6)

To our knowledge, this approach to solve the stress equation in RANS closures was first suggested in the turbulence literature by (19) (see equation (4.4) in this reference). For the general case of time-varying velocity gradient, we note that the auxiliary matrix can be written as a time-ordered exponential (see Refs. (20); (18); (21) for background on this basic matrix function)

 H(t,t′)=Texp+[−∫tt′¯¯¯A(s)ds] .

Equation (5) illustrates clearly the distinct roles played by the production term and the contribution given by . Evaluation of Eq. 5 requires the knowledge of the time history of as well as accurate closures for along the fluid history .

As a next step, one may consider the special case in which the velocity gradient is considered to be constant between the initial time and , and set equal to (e.g.) and simply denoted by . For this approximate situation, the solution of Eq.(6) may be written as an ordinary matrix exponential , where the matrix exponential is defined in the usual way . To simplify further, consider Eq. 5 for the case , i.e. now retaining only the production term. This step eliminates the important isotropization effects of pressure-strain and also the nonlinear diffusion effects of the transport terms. While clearly missing important physics, it is still instructive to observe that this simplification allows the solution (Eq. 5) to be written as:

 \boldmathτ(t)= e−(t−t′)¯¯¯¯A %\boldmath$τ$(t′) e−(t−t′)¯¯¯¯A⊤ . (7)

At this stage it is conceptually advantageous to make connection with Refs. (22); (23); (24), where it is proposed to use conditional statistics to capture the relevant statistics of the SGS stress. For example, in Ref. (23), it is shown that the least-square-error best estimate for the SGS stress is of the form of a multi-point conditional average, namely . The multi-point conditioning variables are, in principle, constituted by the entire (-point) resolved velocity field at scale . To simplify the conditioning, one may limit the information to the past time-history of the local velocity structure. In particular, a good choice that captures much of the local dynamics in a Galilean invariant fashion is the Lagrangian past history of the filtered velocity gradient tensor . The dependence on the Lagrangian time history along a fluid particle advected by the filtered resolved velocity field is thus assumed to be described by with (here and below, the dependence of on spatial position is omitted for clarity). According to these ideas, we define a quasi-optimal SGS stress tensor as the conditional average . Noticing that the matrix exponential prefactors entering in Eq. 7 are some deterministic functions of the velocity gradient tensor itself, they thus can be taken out from this conditional average and we get the following stress tensor:

 \boldmathτ(o)(t)= e−(t−t′)¯¯¯¯A ⟨\boldmathτ(t′)|¯¯¯¯¯A⟩ e−(t−t′)¯¯¯¯A⊤ . (8)

With this expression, the closure problem has been changed from requiring a model for the local stress tensor at time to requiring a model for the conditional average of the ‘upstream initial condition’ at time . The initial condition needed is a symmetric tensor. In the absence of additional information, the simplest assumption is to postulate that this conditionally averaged ‘upstream’ stress tensor is isotropic, namely

 ⟨τij(t′)|¯¯¯¯¯A⟩ ≈ 13⟨τkk(t′)|¯¯¯¯¯A⟩ δij. (9)

The magnitude of the tensor is proportional to the trace of the SGS tensor and has units of squared velocity. The assumption of isotropy may be justified if and become more and more de-correlated as the elapsed time grows, then no locally strong and statistically preferred direction should exist. This step introduces a characteristic decorrelation time-scale , and will be chosen to be of the order of such a decorrelation time-scale. Clearly, one must also assume local isotropy to hold for the statistics, and this is justified from the usual arguments in turbulence when is sufficiently small compared to the integral scale. Incidentally, it is expected that a decorrelation between and may occur due to pressure effects, turbulent diffusion, etc. Some numerical evidence for such decorrelation and isotropization is provided in the next section.

The trace of the conditional SGS tensor, , must still be specified. The simplest option that is consistent with a local evaluation of velocity and length-scales is to choose a factor proportional to , where is the filtered strain rate tensor and . Finally, replacing into Eq. 8 with , we get

 \boldmathτ(o)=cexpΔ2|¯¯¯¯S|2e−τa¯¯¯¯Ae−τa¯¯¯¯A⊤ , (10)

where the parameter is unknown and may be obtained by empirical knowledge, or by generalizing the dynamic model (25).

For completeness and clarity, we remark that the matrix exponential solution may equivalently be obtained by solving the linearized equation for a turbulent fluctuation that only keeps the Rapid Distortion term from the large-scale velocity field, and neglects all other effects. That is to say, we solve formally the equation using the matrix exponential function. The solution is then multiplied by its transpose to form which is then averaged over the fluctuating initial condition (conditioned on a constant ). The averaging of the term yields the initial (‘upstream’) stress tensor , and with the conditional averaging, an expression equivalent to Eq. 8 is obtained. This is similar to the equivalence between solving the equation for co-variances or for the fluctuations and then averaging, as noted in the context of stochastic linear systems in (18).

Equation (10) represents a closure for the SGS stress expressed in terms of matrix exponentials instead of the more commonly used algebraic closures (26). In section IV, a connection is noted between the expression Eq. (10) and a physical closure for the subgrid stress tensor based on the recent fluid deformation closure in the Lagrangian frame.

## Iii Empirical Evidence of Lagrangian decorrelation and Isotropy from DNS

In order to verify whether the decorrelation and isotropization of conditional averages of SGS stresses occur in turbulence, we analyze a DNS dataset of forced isotropic turbulence. The simulation is conducted using a pseudo-spectral method in a box. grid points are used. Fourier modes in shells with are forced by a term added to the Navier-Stokes equations which provides constant energy injection rate . The viscosity of the fluid is . Data is collected after the simulation reaches statistical steady state. Note that is the SGS stress at a previous time and at the spatial location occupied by the fluid particle which is at the position at time (i.e. , in the notations of §IV). According to the transport equation for (Eq. 4), the fluid particle is advected by the filtered velocity field. Thus, the position of the fluid particle at is found by backward particle tracking starting from end-time in the filtered velocity field. To perform backward particle tracking, the filtered velocity and SGS stress fields are calculated and stored at every , corresponding to of the Kolmogorov time scale. A Gaussian filter is used with filter scale , where is the Kolmogorov length scale. In order to quantify isotropy as function of , the ratio of off-diagonal to on-diagonal tensor elements of the conditional averaged SGS stress at decreasing previous time is computed.

According to the derivation, the averaging must be conditioned on a particular value of the resolved velocity gradient . There are a large number of possibilities, since has 8 independent elements. As representative of an important class of velocity gradient structure, we choose to consider regions where the is such that it has a large shear in one direction, whereas all other velocity gradient tensor elements are weak. We choose a particular shear direction, “12”, and define to be a “high-12-shear” events that occur at time . These events are defined here as those points where 1) , i.e. large and positive 12-shear, 2) for other off-diagonal components and , and 3) for the diagonal elements . The gradient rms is defined as . This definition allows a sufficiently large number of events to be counted and thus help in reaching statistical convergence. With this definition of a conditioning event, we calculate the isotropy factor according to:

 I(t−t′)≡−⟨τ12(t′)|E12(t)⟩(1/3)⟨τkk(t′)|E12(t)⟩. (11)

monitors the isotropization of the SGS stress associated with large “12 shear events”, i.e. a particular anisotropic condition in the large-scale velocity gradient tensor. Since the turbulence is statistically isotropic, similar results are expected if the other two shear component of , namely and , had been chosen instead of , under conditioning based on events or , respectively.

Backward particle tracking starts from spatial locations where the conditions in are verified, at time . At each time , the particle locations are calculated from the stored filtered velocity fields using a second-order Adam-Bashforth scheme. The filtered velocity and SGS stresses at the particle locations are interpolated from the stored fields using 6th order Lagrangian interpolation. The conditional averages are then found by averaging over all tracked particles. Statistical sampling is increased by averaging over the trajectories starting from several different end times and also over the other two and off-diagonal elements (in both and ).

The resulting ratio is plotted in Fig. 1 as a function of the normalized time lag . It is evident that the conditional average of the SGS stresses becomes more isotropic as the time lag increases and crosses zero at about 0.7 eddy turn-over times, namely . Then there is negative undershoot to about negative half of the initial value, before it is relaxed to around zero (the isotropic value) at about . The undershoot below zero is an interesting trend and understanding the physics of this behavior would be an interesting goal for future studies.

As additional evidence for the Lagrangian time decorrelation between stress and large-scale velocity gradient, in Fig. 2 we show the correlation coefficient between and , namely

 ρ=−⟨τ12(t−t′)¯¯¯¯A12(t)⟩√⟨τ12(t−t′)2⟩⟨¯¯¯¯A12(t)2⟩ . (12)

The correlation is near 15% at zero time-lag (similar to the correlation coefficient between SGS stresses and strain-rate tensor often quoted in a-priori studies), but then decays to nearly zero at times around . Taken together, the DNS analyses thus provide evidence for the isotropy assumption on , as long as with .

Note that due to the cost of storing the entire simulation for backward particle tracking, only moderate Reynolds numbers were considered in the analysis. The forcing length scale has been estimated to be about 50 times the Kolmogorov length scale, , and the viscous effects begin to significantly damp the motions at scales of about and smaller. Therefore, using , 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 less, prevalent at higher Reynolds numbers. We point out that opportunities for much more in-depth future analyses of such issues are provided by the availability of a turbulence database at higher Reynolds number (29) (although this database could not be used for the present data analysis due to the fact that it does not yet contain sufficiently efficient means of filtering the data).

## Iv Stress Tensor Model Based On The Recent Fluid Deformation Closure

This alternative approach is based on relating the SGS stress tensor to small-scale velocity gradients. To begin, one may recall the multiscale expansion (27); (28); (30) in which, among others, the exact subgrid stress (Eq. 2) is written in terms of , the velocity field coarse-grained at a scale , but still with , i.e. containing significant contributions from sub-grid scales. One may then define the approximated stress tensor and naturally .

Consistent with the Kolmogorov phenomenology, as argued formally in (31), and also as used in various a-priori analyses of experimental data (see e.g. (28)) the SGS stress is relatively local in scale, stating that the leading terms entering in its development are given by the coarse-grained velocity at the resolution scale and including also the next range of length-scales between and (e.g. ). As a consequence, one may use the approximation . Furthermore assuming that is sufficiently smooth over distances (or using the ‘coherent subregion approximation’ (31)), a Taylor expansion of and evaluation of the filtering operation at scale in Eq. 2 leads to

 τij≈C2Δ2 ∂uδi∂xk∂uδj∂xk,     δ=Δβ . (13)

One observes that similarity-type models such as the standard nonlinear model (32); (2) correspond to using the gradient of the large-scale velocity field ( or ). Nevertheless, it is the case which is physically relevant since the true SGS stress includes scales smaller than . However, for , the expression 13 does not constitute a closure since then contains sub-grid motions that are not known at the LES filter scale .

As in (5) and Chevillard & Meneveau (2006 – CM06 from here on), a Lagrangian label position is employed to encode the time-history information. Using the two-time formulation of (3), the label positions satisfy with . Thus represents the position at a prior time of the fluid particle which is at position at time Making the Eulerian-Lagrangian change of variables also used in CM06 leads to the following expression:

 τij=C2Δ2 ∂Xp∂xk∂Xq∂xk∂uδi∂Xp∂uδj∂Xq. (14)

All terms in this expression are strongly fluctuating variables. But, as in prior section, the most relevant information is retained by the conditional averaged expression. We propose the same conditional averaging based on the time-history of the velocity gradient tensor along the past fluid particle trajectory. Therefore, combining the conditional averaging and the change of variables one may write

 τ(o)ij(t)=C2Δ2⟨∂Xp∂xk∂Xq∂xk∂uδi∂Xp∂uδj∂Xq | ¯¯¯¯¯A(s);t′

where the dependence of stress on current position is understood and not indicated to simplify the notation. The Jacobian matrix satisfies (see for instance (17))

 DtG(t′,t)=−G(t′,t)¯¯¯¯¯A(t) % with G(t,t)=I , (16)

where I is the identity matrix. Thus,

 G(t′,t)=Texp−[−∫tt′¯¯¯A(s)ds]

is expressed as an “anti-time-ordered exponential”, with matrices ordered from left to right for increasing times (20); (21); (33). The only difference with the matrix function of the preceding section (Eq. (6)) is the sense of time-ordering.

Since the deformation gradient tensor is a deterministic function of the past velocity gradient history, these tensors can be taken outside the conditional averages in Eq. 15. So far one can thus write

 τ(o)ij(t)=C2Δ2 ∂Xp∂xk∂Xq∂xk Yijpq
 with~{} Yijpq=⟨∂uδi∂Xp∂uδj∂Xq | ¯¯¯¯¯A(s);t′

where is a 4th rank Lagrangian gradient tensor. At this stage, it is now possible again to invoke Lagrangian isotropy, following the approach of CM06 and of the preceding section. It is assumed that the tensor is isotropic due to loss of information caused by turbulent dispersion, past pressure effects, etc. if is long enough. Under the Lagrangian-isotropy closure assumption, one may write

 Yijpq=A′δijδpq+B′δipδjq+C′δiqδjp. (17)

While individual realizations of a small-scale gradient tensor in turbulence are of course not isotropic, statistical moments such as the conditional average can be more justifiably approximated as isotropic. The isotropy assumption states that the rate of change of turbulent velocities (at the present location and time ), with respect to changes in past locations of the fluid particles at time , is insensitive to orientation of . This appears to be a plausible postulate, if sufficient time has elapsed, i.e. if is sufficiently large for decorrelation to take place. In the preceding section we used data from DNS to test the accuracy of such a de-correlation and ‘Lagrangian isotropy assumption’ in a closely related context (directly based on the stresses rather than small-scale velocity gradient statistics). Still, it is important to recognize that this step is introduced here as an ‘ad-hoc’ closure assumption and no claim is made that this is a formal step with controlled errors.

While the assumption of isotropy eliminates the dependence of upon , the latter still affects the Jacobian matrix that enters in the closure for the SGS stress. Next, we focus attention only on the trace-free part of the modeled SGS stress tensor, i.e. we subtract the trace of the stress. And, noticing that the ‘right’ Cauchy-Green tensor is symmetric, only the unknown enters in the resulting ‘quasi-optimal’ model for the deviatoric part of the stress (superscript ) model . Dimensionally, the parameter has units of inverse time-scale squared, and depends upon the turbulence statistics down to scales . For a fixed ratio , and with both scales in the inertial range, for simplicity we assume that the parameter follows, as in the prior section, ‘Smagorinsky scaling’, i.e. . One thus obtains

 τ(od)ij(t)=cexp Δ2|¯¯¯¯S|2 (∂Xi∂xk∂Xj∂xk−13∂Xm∂xk∂Xm∂xkδij), (18)

where the parameter , in a similar fashion as in the preceding section, is unknown and may be obtained by empirical knowledge, or by generalizing the dynamic model (25).

As a final step, the Recent Fluid Deformation (RFD) approximation is used (CM06) in which the time-varying velocity gradient between and is approximated with a constant value (e.g.) equal to its value at and denoted by . The initial condition for the fluid deformation (when the deformation gradient tensor is assumed to be the identity), is prescribed at the time . The solution to Eq. 16 can then be written as . Note that in this approximation, since the sense of ordering of the matrix products is no longer significant.

The next step is to replace the solution for into Eq. 18. And, as was done in the preceding section, to assume that a characteristic de-correlation time-scale has elapsed between the time where the initial isotropy assumption is justifiable and the current time when the stress closure is required. This means replacing the initial time with . Finally, the closure for the deviatoric part of the stress reads

 \boldmathτ(od)ij=cexpΔ2|¯¯¯¯S|2[e−τa¯¯¯¯Ae−τa¯¯¯¯A⊤]d (19)

where all quantities are evaluated at . It is immediately apparent that this closure is equivalent to the formal solution developed in the previous chapter (see Eq. (10)).

## V Expansions

In preceding sections it has been shown that a matrix-exponential closure for the deviatoric part of the SGS stress tensor may be written as in Eq. 19. As a next step, the behavior of this closure is explored when is small enough so that the norm of is much smaller than unity. Then . Up to second order one then obtains

 \boldmathτod ≈cexpΔ2|¯¯¯S|2(−2 τa ¯¯¯S + τ2a [¯¯¯A ¯¯¯A⊤+12(¯¯¯A2+(¯¯¯A⊤)2)]d +...), (20)

Crow, in Ref. (19), derived essentially the same result (see his eq.(5.3)) but with unspecified coefficients obtained as moments of his memory kernel and an additional term proportional to the material derivative It is immediately apparent that if the time-scale is chosen as , then the first term is the standard Smagorinsky model with (where is the Smagorinsky coefficient). Furthermore, the second term, the term in the square parentheses, is of the form of the ‘nonlinear model’ (32); (28); (2) with a prefactor . Two differences with the standard ‘nonlinear model’ are apparent, however. The first is that if , then as coefficient of the nonlinear term this is significantly smaller than the coefficient for this term normally mentioned in the literature (which ranges typically between 1/12 to 1/3). The second difference is the presence of the additional term . To make connections with standard non-linear models used more often in RANS (e.g. (34)), the velocity gradient is decomposed into symmetric and antisymmetric parts, . The result is (again with )

 \boldmathτod≈−2c2sΔ2|¯¯¯S|¯¯¯S+c2sΔ2[¯¯¯S2+12(¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯\boldmathΩ ¯¯¯S−¯¯¯S ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯\boldmathΩ)]d . (21)

It is interesting to note that the expansion including the term cancels exactly the part that is included in the standard non-linear model . For detailed a-priori studies of the various decompositions of the velocity gradient and non-linear terms see (35).

## Vi Matrix exponential closure in LES of isotropic turbulence

The expansion introduced in the last section is formally valid only for small values of the norm of . For more realistic larger values, the expansion may be inaccurate and many additional higher-order terms are needed. They can all be expressed in terms of expansions into integrity bases (26), but it is in general difficult to obtain the coefficients of the expansion. Instead, it is proposed here to utilize the matrix exponential directly in simulations. Since the exponential involves the full velocity gradient tensor, it appears more natural to choose the time-scale according to instead of using . The parameter is an empirical coefficient of order unity.

As a first test, LES of forced isotropic turbulence is performed. This flow is the simplest possible test-case and it is used here simply to determine whether simulations using the matrix-exponential based closure are numerically stable yielding realistic energy spectra, and to ascertain the associated computational cost. The generalization to dynamic versions and tests in more complex flows will be left for future investigations. The simulation uses the same pseudo-spectral method as was used in the DNS outlined in §II, with same grid resolution, forcing scheme and time step size. Dealiasing is performed by zero-padding according to the two-third rule. The viscosity of the fluid is . The subgrid-scale model implemented is given by Eq. 19 and is chosen (dynamic versions (25) of this model to determined can be developed in the future). To specify , the values =0.5, 1 and 2 are tested (a dynamic approach of determining could also be developed). In the pseudo-spectral scheme, the modeled SGS stress is evaluated in physical space and is made trace-free (this only affects the effective pressure, not the dynamics) before computing its divergence in Fourier space.

The matrix exponentials are evaluated using truncated Taylor expansion with scaling and squaring (36). Specifically, we need to evaluate , where . For a matrix in general, the th order truncated Taylor expansion uses matrix polynomial to approximate , incurring an error bounded by . The error decreases with the norm of the matrix . Therefore, to evaluate , we first define , where the value of the integer is chosen to ensure . is then approximated by and finally is given by . The cost of calculating is reduced by using Cayley-Hamilton theorem to express () in terms of , , , and the invariants of . Choosing , we obtain the following equation for with an error smaller than :

 T7(C)=C0I+C1C+C2C2 (22)

where

 C0=1−RC3!+QCRC5!−Q2CRC7!, C1=1−QC3!−RC4!+Q2C5!+2QCRC6!+R2C−Q3C7!, C2=12−QC4!−RC5!+Q2C6!+2QCRC7! .

Here and are the two non-zero invariants of (note that ). In terms of cost, the above algorithm uses about flops to calculate when is given, where is the dimension of the matrix. In our tests , so when and the cost is estimated at about 140 flops for each stress evaluation. This can be compared with the single matrix multiplication needed for the nonlinear model, which is about flops. Overall with this closure, our code took about twice as long to run as compared to using the mixed model.

Simulations were initialized with random Fourier modes and evolved until statistical steady state was obtained. No numerical instabilities were observed for the three parameter cases considered (, 0.5, 1 and 2). In Figure 3 the energy spectra obtained from the three simulations as averaged in the time interval between one and three large-eddy turnover times are shown. Figure 4 shows the time-evolution of the derivative skewness coefficient . As can be seen, the case appears to yield physically meaningful results, which can be compared with the well-known results of the Smagorinsky model and the mixed model (see, for example, (40)). But, there is clear dependence on the parameter . The skewness coefficient quickly drops to values near for and for . These are realistic values for filtered turbulence (37). The skewness values for , on the other hand, appear to be too close to zero, consistent with some pile-up of the spectrum at high wave-numbers.

## Vii Discussion

A new closure based on matrix exponentials and assumptions about short-time Lagrangian dynamical evolution has been proposed. Matrix exponentials as formal solution of the stress transport equation provides interesting insights into the effects of the production (gradient-stretching) term. Historically, in the context of RANS modeling using additional transport equations, the (closed) production term has justifiably not been the focus of attention in the literature. In LES, however, due to practical constraints ‘algebraic’ closures are most often preferred. The present approach shows that the effects of production in the context of such closures may be taken into account directly based on an exact solution of the stress transport equation. A central step in the present approach is to use isotropy for the ‘upstream’ initial condition. Evidence for such isotropization of initial condition, given present large-scale velocity gradients, has been obtained using a DNS database. Implementation of the closure in LES of forced isotropic turbulence yielded good results. The computational cost is significant, but it is not prohibitive. Since our code with this model took about twice as long to run as with a traditional algebraic closure, LES with this model at a resolution of has similar CPU cost as LES with a traditional model run at a resolution of .

It is crucial to stress that the additional, more subtle physics of the remaining terms in Eq.4 (pressure effects, turbulent diffusion, dissipation, etc..) are, in general, unlikely to be well-represented by the simple assumption of ‘upstream’ isotropy. In addition, non-equilibrium conditions in which varies quickly along the particle trajectory are not included in the closure as written in Eq. 19, in which the velocity gradient is assumed to have remained constant over a time-scale . To explore non-equilibrium effects, the full time-ordered exponential function must be used, although this would still leave out the non-equilibrium effects of . To compare the present approach to other closures will require more in-depth testing in more demanding, complex flows (e.g. where effects of anisotropy, non-equilibrium, and pressure-strain correlations are expected to be important).

It is also instructive to consider the case of two-dimensional (2D) turbulence. Nothing in the closure strategies pursued here limits their application to space dimension three, at least nothing very obvious. However, the expansions (Eqs. V,21) show that this is not likely to be a qualitatively good closure for space dimension two, since one there expects an effective “negative eddy-viscosity” corresponding to inverse energy cascade (39). It is thus worth reflecting on some of the reasons for the inaccuracy of the closures in 2D, since this may help pinpoint potential shortcomings in 3D as well. First, it is known that the 2D inverse cascade is less local than the 3D forward cascade, with most of the flux coming from triadic interactions for a scale-ratio (39); (38). However, the starting point of the RFD closure, Eq. 13, is not accurate for To get a qualitatively reasonable alternative at substantially larger than 1—which involves only first-order gradients—one must instead use something like the “Coherent Subregions Approximation” of (39). On the other hand, the starting point of the closure approximation in Section II, the stress transport equation 4, is exact in 2D just as in 3D. The failure of the closure procedure in 2D is now due, presumably, to the effects of the source-terms in the transport equation. Indeed, those terms are expected to contribute as an effective “negative viscosity”, primarily due to the pressure-Hessian rotating small-scale strain matrices relative to the large-scale strain (39). Note that, strictly speaking, this is probably also true in 3D, so that the matrix-exponential closures are likely to be overly dissipative in every dimension. The main effect of the gradient stretching terms—which is a tendency to forward cascade, or positive eddy-viscosity—is well captured by the matrix-exponential closure in any dimension, but the additional, more subtle physics of the remaining terms in Eq. 4 are most likely not well-represented by the simple assumption of isotropy.

As has been cautioned several times in this paper, the simplified matrix-exponential closure as written in Eq. 10 employs the drastic approximation of entirely omitting the pressure-strain correlation and other ‘nonlinear scrambling’ terms. But unlike eddy-viscosity based closure assumptions, this expression can be derived directly from a relevant fluid dynamical equation, namely the stress transport equation (with only the production term), and using physically motivated and straightforward assumptions about Lagrangian decorrelation and upstream isotropy. A similar result is obtained using an Eulerian-Lagrangian change of variables when the stress is expressed in terms of subgrid-scale velocity gradients. Perhaps it can be expected that casting this new light on the closure problem improves our understanding of this long-standing problem.

Finally, we remark that many transport equations for turbulence moments have a basic structure similar to Eq. 4, including two production terms involving the velocity gradient and its transpose. Examples include higher-order moments of velocity, the spectral tensor encountered in Rapid Distortion Theory calculations, etc… The formal solution in terms of matrix exponentials provides new possibilities of calculation and insights into the underlying physics.

###### Acknowledgements.
We thank S.B. Pope and R. Rubinstein for useful comments. L.C. thanks the Keck Foundation for financial support. The other authors thank the National Science Foundation for financial support, and C.M. also acknowledges partial support from the Office of Naval Research.

### References

1. M. Lesieur, M and O. Métais. New trends in large-eddy simulations of turbulence. Ann. Rev. Fluid Mech. 28, 45 (1996).
2. C. Meneveau and J. Katz, Scale-Invariance and Turbulence Models for Large-Eddy Simulation, Ann. Rev. Fluid Mech. 32, 1 (2000).
3. R.H. Kraichnan, Lagrangian-History Closure Approximation for Turbulence, Phys. Fluids 8, 575 (1965).
4. H. Tennekes and J.L. Lumley A first course in turbulence, MIT press, Cambridge (1972).
5. S.S. Girimaji, Modeling Turbulent Scalar Mixing as Enhanced Diffusion, Combust. Sci. Tech. 97, 85 (1994).
6. M. Chertkov, A. Pumir and B.I. Shraiman, Lagrangian tetrad dynamics and the phenomenology of turbulence Phys. Fluids 11, 2394 (1999).
7. Y. Li and C. Meneveau, Intermittency trends and Lagrangian evolution of non-Gaussian statistics in turbulent flow and scalar transport, J.Fluid Mech. 558, 133 (2006).
8. E. Jeong and S. S. Girimaji, Velocity-Gradient Dynamics in Turbulence: Effect of Viscosity and Forcing, Theor. Comput. Fluid Dyn. 16, 421 (2003).
9. A. Pumir and B. Shraiman, Lagrangian Particle Approach to Large Eddy Simulations of Hydrodynamic Turbulence, J. Stat. Phys. 113, 693 (2003).
10. L. Chevillard and C. Meneveau, Lagrangian dynamics and statistical geometric structure of turbulence, Phys. Rev. Lett. 97, 174501 (2006).
11. L. Chevillard and C. Meneveau, Intermittency and universality in a Lagrangian model of velocity gradients in three-dimensional turbulence, C. R. Mécanique 335, 187 (2007).
12. L. Chevillard, C. Meneveau, L. Biferale, F. Toschi, Modeling the pressure Hessian and viscous Laplacian in Turbulence: comparisons with DNS and implications on velocity gradients dynamics, Phys. Fluids 20, 101504 (2008).
13. A. Leonard, Energy cascade in large eddy simulations of turbulent fluid flows, Adv. Geophys. 18, 237 (1974).
14. J.W. Deardorff, Three-dimensional numerical study of the height and mean structure of a heated planetary boundary layer, Bound. Layer Meteor. 7, 81 (1974).
15. M. Germano, Turbulence: the filtering approach, J. Fluid Mech. 238, 325 (1992).
16. C. Penland, Random forcing and forecasting using principal oscillation pattern analysis, Mon. Wea. Rev. 117, 2165 (1989).
17. C. Truesdell and W. Noll, The Non-Linear Field Theories of Mechanics, Springer-Verlag, Berlin, 1992.
18. G.L. Eyink, Linear stochastic models of nonlinear dynamical systems, Phys. Rev. E. 58, 6975 (1998).
19. S.C. Crow, Viscoelastic properties of fine-grained incompressible turbulence, J. Fluid Mech. 33, 1 (1968).
20. J.D. Dollard and C.N. Friedman, Product Integration, in Encyclopedia of Mathematics, Addison-Wesley, London (1979).
21. W.J. Rugh, Linear System Theory, 2nd Ed. Prentice Hall , New Jersey, 1996.
22. R.J. Adrian, Stochastic estimation of sub-grid scale motions, Appl. Mech. Rev. 43, 214 (1990).
23. J. Langford and R. Moser, Optimal LES formulations for isotropic turbulence, J. Fluid Mech. 398, 321 (1999).
24. S.B. Pope, Turbulent Flows, Cambridge Univ. Press,Cambridge, 2000.
25. M. Germano, U. Piomelli, P. Moin and W.H. Cabot, A dynamic subgrid-scale eddy viscosity model, Phys. Fluids A, 3, 1760 (1991).
26. S.B. Pope, A more general effective-viscosity hypothesis, J. Fluid Mech. 72, 331 (1975).
27. R.H. Kraichnan, On Kolmogorov’s inertial-range theories, J.Fluid Mech. 62, 305 (1974).
28. S. Liu, C. Meneveau and J. Katz, On the properties of similarity subgrid-scale models as deduced from measurements in a turbulent jet, J. Fluid Mech. 275, 83 (1994).
29. Y. Li, E. Perlman, M. Wan, Y. Yang, R. Burns, C. Meneveau, R. Burns, S. Chen, A. Szalay and G. Eyink, A public turbulence database cluster and applications to study Lagrangian evolution of velocity increments in turbulence, J. Turbulence 9, 31 (2008).
30. G.L. Eyink, Locality of turbulent cascades, Physica D 207, 91 (2005).
31. G.L. Eyink, Multi-scale gradient expansion of the turbulent stress tensor, J. Fluid Mech. 549, 159 (2006).
32. R.A. Clark, J.H. Ferziger and W.C. Reynolds, Evaluation of subgrid-scale models using an accurately simulated turbulent flow, J. Fluid Mech. 91, 1 (1979).
33. G. Falkovich, K. Gawȩdzki and M. Vergassola, Particles and fields in fluid turbulence Rev. Mod. Phys. 73, 913 (2001).
34. C.G. Speziale, Analytical Methods for the Development of Reynolds-Stress Closures in Turbulence, Ann. Rev. Fluid Mech. 23, 107 (1991).
35. K. Horiuti, Roles of non-aligned eigenvectors of strain-rate and subgrid-scale stress tensors in turbulence generation J. Fluid Mech. 491, 65 (2003).
36. C. Moler and C.F. Van Loan, Nineteen dubious ways to compute the exponential of a matrix, twenty-five years later, SIAM Rev. 45, 3 (2003).
37. S. Cerutti and C. Meneveau, Statistics of filtered velocity in grid and wake turbulence, Phys. Fluids 12, 1143 (2000).
38. S. Chen, R.E. Ecke, G.L. Eyink, M. Rivera, M. Wan, and Z. Xiao, Physical Mechanism of the Two-Dimensional Inverse Energy Cascade, Phys. Rev. Lett. 96, 084502 (2006).
39. G.L. Eyink, A turbulent constitutive law for the two-dimensional inverse energy cascade, J. Fluid Mech. 549, 191 (2006).
40. H.-S. Kang, S. Chester, and C. Meneveau, Decaying turbulence in an active-grid-generated flow and comparisons with large-eddy simulation, J. Fluid Mech., 480, 129 (2003).
You are adding the first comment!
How to quickly get a good reply:
• Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
• Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
• Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minumum 40 characters