Pressure Hessian and viscous contributions to velocity gradient statistics based on Gaussian random fields
UUnder consideration for publication in J. Fluid Mech. Pressure Hessian and viscous contributionsto velocity gradient statistics based onGaussian random fields
Michael Wilczek † , and Charles Meneveau Department of Mechanical Engineering & Institute for Data Intensive Engineering andScience, The Johns Hopkins University, 3400 North Charles Street, Baltimore MD 21218(Received 3 September 2014)
Understanding the non-local pressure contributions and viscous effects on the small-scalestatistics remains one of the central challenges in the study of homogeneous isotropicturbulence. Here we address this issue by studying the impact of the pressure Hessian aswell as viscous diffusion on the statistics of the velocity gradient tensor in the frameworkof an exact statistical evolution equation. This evolution equation shares similarities withearlier phenomenological models for the Lagrangian velocity gradient tensor evolution,yet constitutes the starting point for a systematic study of the unclosed pressure Hessianand viscous diffusion terms. Based on the assumption of incompressible Gaussian velocityfields, closed expressions are obtained as the results of an evaluation of the characteristicfunctionals. The benefits and shortcomings of this Gaussian closure are discussed, anda generalization is proposed based on results from direct numerical simulations. Thisenhanced Gaussian closure yields, for example, insights on how the pressure Hessianprevents the finite-time singularity induced by the local self-amplification and how itsinteraction with viscous effects leads to the characteristic strain skewness phenomenon. † Email address for correspondence: [email protected] a r X i v : . [ phy s i c s . f l u - dyn ] S e p M. Wilczek and C. Meneveau
1. Introduction
The understanding of turbulence dynamics and statistics by theoretical means is ham-pered by two challenges: nonlinearity and non-locality. These challenges become espe-cially evident when studying the smallest scales of turbulence in terms of the velocitygradient tensor A = ∇ u (where u is the velocity vector). The velocity gradient ten-sor gives a comprehensive characterization of the small scales. Its symmetric part, therate-of-strain tensor, characterizes the local rates of deformation of fluid elements and,for example, determines the rate of kinetic energy dissipation. In addition to nonlinearadvection, viscous diffusion and the interaction with the vorticity field, the rate-of-straindynamics is also subject to non-local pressure effects through the pressure Hessian. Theisotropic part of the pressure Hessian preserves solenoidality of the velocity field, whereasthe deviatoric part communicates information between distant points in the velocity gra-dient field. The antisymmetric part of the velocity gradient tensor, the rate-of-rotationtensor, represents the vorticity and yields further insights into the small-scale coherentstructures of the field. Also the vorticity is subject to nonlinear advection through thevelocity field and viscous diffusion, but is locally stretched by the rate-of-strain tensor.With its wealth of information, the velocity gradient tensor has been subject to re-search for many decades. Whereas early works by Betchov (1956) revealed a number ofimportant kinematic constraints on its statistical properties, the dynamical properties ofthe velocity gradient tensor moved into focus by the works of Vieillefosse (1982, 1984)and Cantwell (1992). In these works, the local self-amplification of the velocity gradi-ent tensor along Lagrangian fluid trajectories has been studied in detail neglecting thenon-local pressure contributions and dissipative effects. This so-called restricted Eulerapproximation led to valuable insights into several features of small-scale turbulence. Forexample, these studies elucidated the preferential alignment of the vorticity vector withthe principal axes of the rate-of-strain tensor, which is also observed in direct numericalsimulations (DNS) studies (Ashurst et al. et al. (1998) that the inclusion of a linear diffusion term preventsthe singularity for initial conditions with moderate velocity gradient tensor values (ascompared to the damping rate), while for larger values the quadratic nonlinearity stilloverpowers the linear damping.These prior works revealed that the velocity gradient dynamics can be convenientlystudied in a Lagrangian frame. This idea has been picked up and considerably extendedby the works of Chertkov et al. (1999) and Naso & Pumir (2005), who argued that atetrad, whose corners are defined by four Lagrangian fluid particles, can be interpreted asa scale-dependent perceived velocity gradient tensor, giving insights into the statisticalproperties of not only the dissipative but also the inertial range of scales. The governingequations of motions have been closed by phenomenological arguments, where it hasbeen assumed that the main effect of the pressure Hessian is to deplete the nonlinearself-amplification of the gradients.Another related phenomenological approach, focusing on the smallest scales of fluidmotion, has been proposed by Chevillard & Meneveau (2006) who developed a closurebased on the deformation of an infinitesimal fluid element. By taking into account onlyits recent history, closed expressions for the pressure Hessian and the viscous term havebeen obtained. While it has been found that this model yields a good description ofmany observed features, it has also been shown in a critical comparison with DNS databy Chevillard et al. (2008) that the model misses some features of the pressure Hessian.Further limitations concern the behaviour for large Reynolds numbers. For a more ex-tensive literature review including other phenomenological models, see Meneveau (2011).The quality of the phenomenological models summarized here has reached a levelwhere satisfactory agreement of numerical and experimental data is achieved for mod-erate Reynolds numbers. The complexity of both the closure problem itself and thealready proposed models, however, renders progress based on further phenomenologi-cal refinements a challenge and motivates further theoretical investigations. Moreover, atheoretical justification of the proposed models remains elusive and motivates a deeperinvestigation of the pressure Hessian and viscous contributions.The current work aims at contributing to such investigations. To this end, we studythe statistical properties of the velocity gradient tensor in terms of an exact, yet unclosedstatistical evolution equation. This evolution equation describes the statistical propertiesof a class of fluid particles that share the same value of the velocity gradient tensor. Theadvantage of this approach is the rigorous formulation of the closure problem in terms ofrandom fields, which serves as a starting point to establish well-controlled closures basedon as few assumptions as possible. To obtain explicit expressions for the pressure Hessianand viscous terms, we then make the assumption that the velocity field can be representedby an incompressible Gaussian random field. This is known to be inaccurate due to, e.g.,small-scale intermittency, but this approach allows for a fully analytical treatment. Inparticular, insights into the formal structure of the pressure Hessian and the viscousterm are obtained without further ad-hoc assumptions. The merits and shortcomings ofthis Gaussian closure are discussed, and a generalization based on DNS observations isproposed. In this enhanced Gaussian closure, dimensionless parameters obtained fromthe analytical solutions are replaced by empirical fits that make use of DNS data. Thedynamical features of this new closure illustrate how the non-local pressure Hessiancontributions help to prevent a finite-time singularity and how the interaction with thediffusive ingredients of the dynamics allow for strain skewness and enstrophy production.Finally, predictions from this enhanced Gaussian closure are compared with DNS data.
2. Velocity gradient dynamics and statistical description
In this section we briefly review the basic equations of motion and then proceed to adetailed description of the statistical methods used in this paper.2.1.
Velocity gradient tensor dynamics and non-local pressure contributions
The dynamics of the velocity gradient A ij ( x , t ) = ∂u i ∂x j ( x , t ) is obtained by taking thespatial gradient of the incompressible three-dimensional Navier-Stokes equation, yielding ∂∂t A + u · ∇ A = − A − H + ν ∆ A + F . (2.1)Here u ( x , t ) is the velocity field, H ij ( x , t ) = ∂ p∂x i ∂x j ( x , t ) is the Hessian of the kine-matic pressure field, ν is the kinematic viscosity, and F ( x , t ) constitutes the gradientof a solenoidal large-scale forcing optionally included in the Navier-Stokes equation.This equation states that the velocity gradient tensor field is advected with the veloc-ity field while being subject to self-amplification or self-attenuation, as well as pressureand viscous diffusion effects and optionally an external forcing. The influence of theself-amplification term is well understood thanks to its locality, and has been analysedextensively in the context of the restricted Euler model, for example, by Cantwell (1992).The viscous term already incorporates some non-local information, because information M. Wilczek and C. Meneveau from neighbouring points in the field is needed to evaluate the Laplacian. The pressureHessian term is related to the velocity gradient tensor field by a Poisson equation,∆ p = Tr ( H ) = − Tr (cid:0) A (cid:1) , (2.2)which is obtained from (2.1) by taking the trace and using the fact that Tr ( A ) = 0 andTr ( F ) = 0. This shows that the isotropic part of the pressure Hessian is also local, suchthat the pressure Hessian can be decomposed into H = −
13 Tr (cid:0) A (cid:1) I + (cid:101) H (2.3)where I is the identity matrix and (cid:101) H is a traceless symmetric tensor containing all thenon-local contributions of the pressure field. This can be made more explicit by the factthat the non-local pressure Hessian contribution is obtained from the velocity gradientfield by a principal value integral (Ohkitani & Kishiba 1995) which (for unbounded flow)reads (cid:101) H ij ( x , t ) = − π (cid:90) P . V . d x (cid:48) (cid:20) δ ij | x − x (cid:48) | − x − x (cid:48) ) i ( x − x (cid:48) ) j | x − x (cid:48) | (cid:21) Tr (cid:0) A ( x (cid:48) , t ) (cid:1) . (2.4)This relation stresses that the deviatoric part of the pressure Hessian contains highlynon-local information from remote points in the fluid. Interestingly, its dependence onthe velocity gradient tensor occurs only through a scalar invariant, Q = − Tr (cid:0) A (cid:1) .For the discussion later on it is also useful to decompose the velocity gradient dy-namics into its symmetric and antisymmetric contributions, S = (cid:0) A + A T (cid:1) and W = (cid:0) A − A T (cid:1) , respectively. In terms of these tensors (for the moment considering a flowwithout body force) the velocity gradient dynamics (2.1) takes the form ∂∂t W + u · ∇ W = − SW − WS + ν ∆ W (2.5) ∂∂t S + u · ∇ S = − (cid:20) S −
13 Tr (cid:0) S (cid:1) I (cid:21) − (cid:20) W −
13 Tr (cid:0) W (cid:1) I (cid:21) − (cid:101) H + ν ∆ S . (2.6)Alternatively, the antisymmetric part is readily expressed in terms of the vorticity vector, ω i = − ε ijk W jk , where ε ijk is the Levi-Civita tensor. Taking S and ω as primary variables,the above equations can be written as ∂∂t ω + u · ∇ ω = S ω + ν ∆ ω (2.7) ∂∂t S + u · ∇ S = − (cid:20) S −
13 Tr (cid:0) S (cid:1) I (cid:21) − (cid:20) ωω T − ω I (cid:21) − (cid:101) H + ν ∆ S . (2.8)The velocity gradient dynamics becomes particularly clear from these equations whenassuming S as initially diagonal without loss of generality. According to (2.7), vorticityis stretched or attenuated depending on the sign of the specific eigenvalue of the rate-of-strain tensor, in addition to advection and viscous diffusion. As can be seen from(2.8), the rate-of-strain is advected and diffuses. More importantly, however, it can benoted that the first term on the right-hand side is diagonal in the eigenframe of the rate-of-strain tensor and causes a self-amplification or self-attenuation of velocity gradientssimilar to that observed in Burgers dynamics. The second term also contributes to thisamplification, but furthermore induces a rotation of the eigenframe depending on thevorticity vector. Understanding how the non-local pressure Hessian acts in this equationis one of the central topics in the present paper.In the following we will explicitly consider a stochastic large-scale forcing. Under theassumption that the small-scale statistics of turbulence is independent of the large-scaleforcing for sufficiently high Reynolds numbers, this particular choice does not appearto be a severe restriction, but it will turn out useful to make a connection to existingphenomenological models. Furthermore, it allows for a fully analytical treatment.2.2. Statistical evolution equation
We now turn to a statistical description of the problem, focusing on f ( A ; t ), the prob-ability density function (PDF) of the velocity gradient tensor at a single point in thefluid. This function, which depends upon the nine velocity gradient tensor elements andtime, contains rich information on the small-scale properties of turbulence, including thesingle-point statistics of vorticity, rate-of-strain and their mutual orientation.To obtain an evolution equation for the probability function, we use the statisticalframework of the Lundgren-Monin-Novikov hierarchy, which allows to derive exact, yetunclosed evolution equations for PDFs. A basic account on the methodology can be foundin Lundgren (1967); Dopazo (1994); Pope (2000); Wilczek et al. (2011) and Friedrich et al. (2012). Appendix A gives a detailed derivation of the relations used in this section,but we here rather focus on a discussion of the physical implications of non-localityon the small-scale statistics of turbulence. The closure problem arising in this statisticalframework can be cast in two different ways. One way is introducing conditional averagesof, e.g., the pressure Hessian, with respect to the velocity gradient, the other way is toexpress the unclosed terms in terms of multipoint statistics. The two formulations givecomplementary perspectives on the problem, and in fact both perspectives turn out tobe useful for the current work.The PDF equation for the velocity gradient tensor in homogeneous turbulence reads(see also Girimaji & Pope (1990)) ∂∂t f ( A ; t ) = − ∂∂ A ij (cid:18)(cid:20) − (cid:18) A ik A kj −
13 Tr (cid:0) A (cid:1) δ ij (cid:19) − (cid:10) (cid:101) H ij (cid:12)(cid:12) A (cid:11) + (cid:10) ν ∆ A ij (cid:12)(cid:12) A (cid:11)(cid:21) f ( A ; t ) (cid:19) + 12 Q ijkl ( ) ∂∂ A ik ∂∂ A jl f ( A ; t ) , (2.9)where Q ijkl ( ) denotes the two-point covariance tensor of the gradients of the large-scaleforcing evaluated at the origin. Without a stochastic forcing (2.9) is a Liouville equationdescribing the conservation of probability. Due to the fact that the self-amplification andthe isotropic part of the pressure Hessian are local, these terms appear closed. The onlyunclosed terms in this equation are the conditional averages of the non-local pressureHessian and the viscous diffusion. The stochastic forcing, which only enters the equationthrough its covariance tensor, turns the PDF equation into a Fokker-Planck equation.The Fokker-Planck equation can be associated to a Langevin equation which establishesthe connection to existing phenomenological models. This equation takes the formd A = (cid:20) − (cid:18) A −
13 Tr (cid:0) A (cid:1) I (cid:19) − (cid:10) (cid:101) H (cid:12)(cid:12) A (cid:11) + (cid:10) ν ∆ A (cid:12)(cid:12) A (cid:11)(cid:21) d t + d F . (2.10)Here, d F is a stochastic forcing which is determined by the large-scale forcing applied tothe velocity field (see appendix A.2 for more details). Again, the closure problem arisesin terms of the conditional averages of the pressure Hessian and the viscous diffusionof the velocity gradient. This exact, yet unclosed equation has clear similarities withthe dynamical phenomenological stochastic models for the Lagrangian velocity gradientevolution reported earlier (Girimaji & Pope 1990; Jeong & Girimaji 2003; Chevillard &Meneveau 2006; Meneveau 2011), which focus on modelling realizations of the process.The difference, however, is that this equation describes the statistical evolution of a class M. Wilczek and C. Meneveau of fluid particles which all share the identical value of the velocity gradient, whereas theLagrangian velocity gradient models express the pressure Hessian and the diffusive termin terms of the velocity gradient tensor along that trajectory.The fact that the closure problem in the current formulation arises in terms of con-ditional averages allows us to study the unclosed terms based on general random fields.While this certainly does not solve the central problem, it will yield some interestinginsights. Moreover, the assumption of Gaussian random fields allows to calculate theunclosed terms analytically. Both of these points will be discussed below in detail.2.3.
Conditional viscous diffusion and conditional pressure Hessian
To make contact with a formulation in terms of general random fields, we first have toestablish the relation of the unclosed terms to multipoint statistics. The unclosed terms in(2.9) and (2.10) can be expressed in terms of two-point statistics, which underscores theirnon-local nature. For the following we will consider homogeneous isotropic turbulence.The technical background on obtaining the relations discussed in this paragraph are givenin appendix A.3.The relation obtained for the viscous term reads (cid:10) ν ∆ x A ( x , t ) (cid:12)(cid:12) A (cid:11) = lim r → ν ∆ r (cid:10) A ( x , t ) (cid:12)(cid:12) A (cid:11) , (2.11)where we have introduced subscripts to discriminate the two points in space and thedistance vector r = x − x . The viscous term contains non-local information in terms ofthe average velocity gradient at point x conditional on the value of the velocity gradientat point x . Still, it can be considered as somewhat local because only the immediateneighbourhood of x has to be known to evaluate the derivative.For the conditional pressure Hessian we obtain (cid:10) (cid:101) H ij ( x , t ) (cid:12)(cid:12) A (cid:11) = 12 π (cid:90) P . V . d r (cid:20) δ ij r − r i r j r (cid:21) (cid:10) Q ( x , t ) (cid:12)(cid:12) A (cid:11) . (2.12)This result is interesting in two respects: First, the structure of the integral kernel as-sures that the conditional Hessian will be traceless and symmetric in i and j , no matterwhat is assumed for (cid:10) Q ( x , t ) (cid:12)(cid:12) A (cid:11) ; second, this unclosed tensorial conditional averagedepends only on a conditional scalar expression, which is a considerable simplification,especially for modelling. Still the conditional second invariant represents a rather com-plicated function, which may depend on all invariants that can be constructed from A and r .By expressing the conditional average in terms of two-point statistics, we have “shifted”the closure problem from an unclosed expression for the joint single-point statistics in-volving the pressure Hessian and the velocity gradient to the joint two-point statisticsinvolving the velocity gradient only. This allows for an evaluation in terms of general ran-dom fields, and we will pursue the special case of Gaussian random fields in the followingsections.We remark that up to now all of the discussed relations represent exact results. Oncea more general theory of random fields is available, they can be used to obtain moreelaborate closures.
3. A closure based on Gaussian random fields
Before presenting the central results of this paper, some words on why to choose aclosure based on Gaussian random fields are in order, especially given the fact thatit is generally known that the multipoint structure of the velocity field in turbulentflows exhibits important non-Gaussian features such as intermittency and skewness ofvelocity increments and gradients (see, e.g., Frisch (1995)). The justification is simple: Asa general theory for non-Gaussian random fields is currently lacking, Gaussian randomfields are the only available choice for an analytical treatment of the current closureproblem without involving further phenomenological assumptions.Furthermore, Gaussian fields serve as an important reference point for comparisonwith statistics from real turbulent flow, as for example discussed by Shtilman et al. (1993); Tsinober (1998, 2009). It is also worth pointing out that for pressure statistics,the assumption of Gaussianity of the velocity field has, in fact, been shown to leadto qualitatively correct results by Holzer & Siggia (1993). One might speculate thata possible reason for this is that some essential features of the non-locality might berobust to the details of the particular choice of random fields. The motivation for thecurrent work is the perspective that new insights on the structure of the unclosed terms,especially on the complex pressure Hessian, can be achieved this way. In fact, our resultswill show that nontrivial, but imperfect results can be obtained under the assumption ofGaussian velocity fields.3.1.
Gaussian characteristic functional for incompressible velocity fields
The following calculations rely exclusively on the spatial properties of Gaussian randomfields, such that we can suppress the time variable in our notation. The assumptionof Gaussian velocity fields can most comprehensively be captured on the level of thecharacteristic functional of the velocity field, which is defined in terms of φ u [ λ ( x )] = (cid:28) exp (cid:20) i (cid:90) d x λ i ( x ) u i ( x ) (cid:21)(cid:29) . (3.1)Here λ ( x ) denotes the Fourier field conjugate to the velocity field. The spatial coordi-nate x may be regarded as a continuous index of the functional Fourier transform. ForGaussian random fields with zero mean, the explicit expression for the characteristicfunctional reads φ u [ λ ( x )] = exp (cid:20) − (cid:90) d x (cid:90) d x (cid:48) λ i ( x ) R uij ( x , x (cid:48) ) λ j ( x (cid:48) ) (cid:21) (3.2)which depends on the velocity covariance tensor R uij ( x , x (cid:48) ) = (cid:104) u i ( x ) u j ( x (cid:48) ) (cid:105) . For homo-geneous isotropic turbulence in incompressible flows this tensor takes the form R uij ( x , x (cid:48) ) = R uij ( r ) = (cid:104) u (cid:105) (cid:20) f u ( r ) δ ij + 12 rf (cid:48) u ( r ) [ δ ij − ˆ r i ˆ r j ] (cid:21) . (3.3)Here, f u denotes the longitudinal velocity autocorrelation function, and ˆ r i = r i /r denotesa component of the direction of the difference vector r = x − x (cid:48) . We have used the factthat for homogenous turbulence the covariance tensor is a function of the distance vectoronly, R uij ( x , x (cid:48) ) = R uij ( r ). Both types of notation will be used in the following, dependingon the context. The statistics of the Gaussian velocity field is completely specified by thelongitudinal velocity autocorrelation function f u ( r ).Next, we consider the characteristic functional of the velocity gradient tensor, whichin analogy to (3.1) is defined as φ A [ Λ ( x )] = (cid:28) exp (cid:20) i (cid:90) d x Λ ij ( x ) A ij ( x ) (cid:21)(cid:29) . (3.4)Owing to A ij = ∂u i /∂x j , a simple relation between (3.1) and (3.4) can be established M. Wilczek and C. Meneveau by partial integration, φ A [ Λ ( x )] = φ u (cid:2) λ ( x ) = −∇ · Λ T ( x ) (cid:3) , (3.5)where the argument of the characteristic functional for the velocity field in componentnotation reads λ i ( x ) = − ∂ Λ ij ∂x j ( x ). This shows that the characteristic functional for thevelocity gradient is readily expressed in terms of the characteristic functional of thevelocity field. Using this result together with the Gaussian characteristic functional (3.2)leads, after partial integration, to φ A [ Λ ( x )] = exp (cid:20) − (cid:90) d x (cid:90) d x (cid:48) Λ ik ( x ) R ijkl ( x , x (cid:48) )Λ jl ( x (cid:48) ) (cid:21) . (3.6)Here we have introduced the velocity gradient covariance tensor, which is kinematicallyrelated to the velocity covariance tensor: R ijkl ( x , x (cid:48) ) = (cid:104) A ik ( x ) A jl ( x (cid:48) ) (cid:105) = (cid:28) ∂u i ∂x k ( x ) ∂u j ∂x (cid:48) l ( x (cid:48) ) (cid:29) = ∂ R uij ∂x k ∂x (cid:48) l ( x , x (cid:48) ) . (3.7)The general structure of this tensor is discussed in appendix B.1. Two conclusions can bedrawn from this result. First, it shows that if the velocity field is (multipoint) Gaussian,so is the velocity gradient tensor field (as cautioned already before, this is known to beunphysical due to intermittency, etc.). Second, because the velocity gradient covariancetensor is kinematically prescribed by the velocity covariance tensor, the full statisticaldescription of the velocity gradient tensor field is fixed by one scalar function, the longi-tudinal velocity autocorrelation function. It is crucial to note that the Reynolds-numberdependence as well as any length-scale information of the Gaussian field enter through thecorrelation function. Its precise shape, for example, especially near the origin, dependson the Reynolds number.For the calculations of the unclosed terms we will in particular make use of single-and two-point statistics. The characteristic functional contains the statistics for arbitrarynumbers of points, such that it can be conveniently projected to the single-point statisticsby evaluating (3.6) at Λ ( x ) = Λ δ ( x − x ): φ A ( Λ ) = φ A [ Λ δ ( x − x )] = exp (cid:20) − Λ T1 R ( ) Λ (cid:21) . (3.8)Here, we have introduced the short-hand notation Λ T1 R ( ) Λ = Λ ,ik R ijkl ( x , x )Λ ,jl .The characteristic function for two points is obtained analogously: φ A ( Λ , Λ ) = φ A [ Λ δ ( x − x ) + Λ δ ( x − x )]= exp (cid:20) − (cid:2) Λ T1 R ( ) Λ + Λ T1 R ( r ) Λ + Λ T2 R ( r ) Λ + Λ T2 R ( ) Λ (cid:3)(cid:21) . (3.9)It may be noted that both characteristic functions resemble a standard Gaussian charac-teristic function generalized to tensorial random variables. The velocity gradient tensorPDFs, whenever needed, can be obtained by Fourier transform. Incompressibility is ex-plicitly taken into account by the structure of the covariance tensor (see appendix B.1).3.2. Gaussian closure
Conditional Laplacian
We are now equipped with the technical prerequisites to evaluate the unclosed termsbased on the assumption of Gaussian random fields. To calculate the conditional viscousterm, we first obtain an expression for (cid:10) A ( x ) (cid:12)(cid:12) A (cid:11) and then evaluate the Laplacian byrelation (2.11). The details of the calculation are elaborated in appendix B.2, which yieldsthe simple result (cid:10) ν ∆ x A ( x ) (cid:12)(cid:12) A (cid:11) = δ A (3.10)with δ = ν f (4) u (0) f (cid:48)(cid:48) u (0) = − ν (cid:82) d k k E ( k ) (cid:82) d k k E ( k ) , (3.11)where E ( k ) is the energy spectrum function. That means, for incompressible Gaussianvelocity fields, the conditional Laplacian is a linear function of the velocity gradient witha prefactor depending on the kinematic viscosity as well as on derivatives of the longitu-dinal autocorrelation function of the velocity, or alternatively on integrals involving theenergy spectrum function. It is interesting to note that the Gaussian closure is consistentwith the linear diffusion model by Martin et al. (1998), thus giving theoretical underpin-ning for this phenomenological assumption: For the case of Gaussian velocity fields, theassumption of a linear dependence of the conditional Laplacian is exact. In addition tothe linear diffusion model, however, the Gaussian closure also fixes the coefficient by ex-pressing it in terms of the longitudinal velocity autocorrelation function or, equivalently,the energy spectrum function.The conditional Laplacian depends implicitly on the Reynolds number through theautocorrelation function. As a simple estimate of the Reynolds-number dependence ofthis term one may note that the coefficient δ defines the inverse of a time scale. If weassume, in accordance with Kolmogorov’s phenomenology, that this time scale dependson the mean rate of energy dissipation ε and viscosity ν , it is necessarily identified withthe Kolmogorov time scale τ η = ( ν/ε ) / . The same can be concluded from (3.11) sinceboth integrals in the numerator and denominator are dominated by viscous scales. Anorder-of-magnitude estimate based on a simple model spectrum described in appendixC yields δ τ η ≈ − .
65 for a range of Reynolds numbers. This value will turn out to beexcessive in magnitude because of an imperfectly modeled viscous range. We will discussa more accurate estimate based on the velocity derivative skewness in the context of theenhanced Gaussian closure later below.3.2.2.
Conditional pressure Hessian
To evaluate the non-local contribution to the pressure Hessian under the assumptionof incompressible Gaussian velocity fields, we first obtain the Gaussian expression for (cid:10) Q ( x ) (cid:12)(cid:12) A (cid:11) and insert this into (2.12). The details of this calculation are included inappendix B.3, so that we can here focus on a discussion of the main result. It takes theform (cid:10) (cid:101) H ( x ) (cid:12)(cid:12) A (cid:11) = α (cid:18) S −
13 Tr (cid:0) S (cid:1) I (cid:19) + β (cid:18) W −
13 Tr (cid:0) W (cid:1) I (cid:19) + γ ( S W − W S ) (3.12)0 M. Wilczek and C. Meneveau with α = −
27 (3.13 a ) β = −
25 (3.13 b ) γ = 625 + 1675 f (cid:48)(cid:48) u (0) (cid:90) d r f (cid:48) u f (cid:48)(cid:48)(cid:48) u r . (3.13 c )First, it is interesting to note that the non-local contribution to the conditional pres-sure Hessian consists of all possible combinations of the rate-of-strain tensor S and therate-of-rotation tensor W that are dimensionally consistent (quadratic), traceless andsymmetric. Second, the fact that the local nonlinear term in the stochastic evolutionequation (2.10) in terms of the rate-of-strain and the rate-of-rotation tensors takes theform − (cid:18) A −
13 Tr (cid:0) A (cid:1) I (cid:19) = − (cid:18) S −
13 Tr (cid:0) S (cid:1) I (cid:19) − (cid:18) W −
13 Tr (cid:0) W (cid:1) I (cid:19) − S W − W S (3.14)shows that the effect of the Gaussian non-local pressure Hessian can, in part, be un-derstood in terms of a “reduction of nonlinearity” as phenomenologically introduced byChertkov et al. (1999) and Naso & Pumir (2005) for the perceived velocity gradient ten-sor dynamics based on Lagrangian tetrads. This result from the Gaussian approximationgives further theoretical support for such an Ansatz. The Gaussian approximation, how-ever, yields differing coefficients for the rate-of-strain term and the rate-of-rotation term,thus suggesting a refinement of the reduction of nonlinearity.A remarkable property of the restricted Euler approximation is that the velocity gra-dient invariants dynamics can be reduced to only two invariants, R and Q , out of fivepossible invariants (Vieillefosse 1982; Cantwell 1992). This is a direct consequence of thefact that only the isotropic, local contribution of the pressure Hessian (which involvesthe trace of the squared velocity gradient) is taken into account in this approximation.In comparison to that, the Gaussian approximation also incorporates deviatoric contri-butions, which cannot be expressed in terms of the identity tensor or powers of A only.This results in a dynamical system that cannot be reduced to R and Q only. In that sensethe structure of the non-local pressure Hessian in the Gaussian approximation suggeststhe possibility of more complex dynamics compared to the restricted Euler model.Regarding the Reynolds-number dependence, it is interesting to note that the coef-ficients α and β do not depend on the autocorrelation function and thus are indepen-dent of the Reynolds number. The situation is more involved for the coefficient γ dueto its integral dependence on the velocity autocorrelation function. A closer look atthis integral relation, however, reveals that it depends on the non-dimensional function f (cid:48) u ( r ) f (cid:48)(cid:48)(cid:48) u ( r ) /f (cid:48)(cid:48) u (0) . Composed of derivatives of the velocity autocorrelation function, itis plausible to assume that this function is largely independent of the energy-containingrange of turbulent motion. If we additionally assume that the inertial and dissipativeranges of the energy spectrum exhibit a universal functional form at sufficiently highReynolds numbers, this implies a universal shape of this non-dimensional function whenproperly rescaled by the Kolmogorov length scale η . Under these assumptions the coef-ficient γ is also universal, i.e. Reynolds-number independent. The actual numeric valueof γ , however, must be obtained from a model autocorrelation function. We have usedthe model spectrum described in appendix C to estimate the parameter γ for a Reynoldsnumber of R λ ≈
432 and obtained γ ≈ .
08. Also the universality argument was verifiedfor a range of Reynolds numbers. Reynolds-number-dependent effects like intermittency1and bottleneck corrections as discussed, e.g., by Meyers & Meneveau (2008), however,will break the universality of the autocorrelation function and the coefficient γ .
4. Dynamics of the Gaussian closure
General considerations
With the results of the last section we have obtained a closure to the statistical evolutionequation (2.10), such that we arrive at a stochastic differential equation (SDE),d A = (cid:20) − (cid:18) A −
13 Tr (cid:0) A (cid:1) I (cid:19) − α (cid:18) S −
13 Tr (cid:0) S (cid:1) I (cid:19) − β (cid:18) W −
13 Tr (cid:0) W (cid:1) I (cid:19) − γ ( SW − WS ) + δ A (cid:21) d t + d F , (4.1)where the coefficients α to δ are fixed by the Gaussian approximation. We have droppedthe subscripts as all quantities are considered at a single point. This Langevin equationdescribes the evolution of a class of fluid particles which share the same value of thevelocity gradient.For the interpretation of the individual terms it turns out to be useful to decomposethe evolution equation for A into the dynamics of the rate-of-strain tensor S and rate-of-rotation tensor W :d S = (cid:20) − (1 + α ) (cid:18) S −
13 Tr (cid:0) S (cid:1) I (cid:19) − (1 + β ) (cid:18) W −
13 Tr (cid:0) W (cid:1) I (cid:19) − γ ( SW − WS ) + δ S (cid:21) d t + d F S (4.2)and d W = [ − SW − WS + δ W ] d t + d F W . (4.3)Here, d F S and d F W denote the contributions of the stochastic forcing to the rate-of-strain and rate-of-rotation dynamics. The rate-of-rotation dynamics written in terms ofthe vorticity takes the form d ω = [ S ω + δ ω ] d t + d F ω . (4.4)To elucidate the dynamics of the closure, we take the restricted Euler model as areference, whose dynamical properties have been discussed in detail by Cantwell (1992)and Nomura & Post (1998). For the vorticity equation, the well-known effect of vortexstretching is supplemented by a linear damping term and the random-force term.The same applies to the rate-of-strain equation. Additionally, the terms related to α and β modify the coefficients of the restricted Euler model, i.e. their relative strengthand correspondingly their relative time scale is changed. As we find typical values in therange − < α < − < β <
0, the effect of the restricted Euler terms is weakened.It is plausible to assume that this affects the occurrence of a finite-time singularity.It is furthermore instructive to consider the effect of the terms individually. The α -termin the rate-of-strain equation can be diagonalized along with the rate-of-strain tensor andthus affects only the eigenvalues of S and not the orientation of the eigenframe. It maybe understood as an algebraic growth or decay of the eigenvalues, which additionallymaintains the zero-trace condition of the rate-of-strain tensor. Considered on its own,it leads to a divergence in finite time, similar to that observed in the inviscid Burgersequation and the restricted Euler model.2 M. Wilczek and C. Meneveau
The β -term is independent of the rate-of-strain tensor, and thus depends linearly ontime for a fixed vorticity. Considered in a frame in which the e -axis coincides with thedirection of the vorticity, this term is diagonal and induces a linear decrease of the firstdiagonal component of S at a rate of − ( β + 1) ω / β +1) ω /
12, assuming β > −
1. As an accompanying effect,the eigendirection corresponding to the most negative eigenvalue tends to align with thedirection of the vorticity. It turns out that this term is crucial for preventing a blow-upof the SDE system. As a dynamically evolving vorticity tends to be amplified along thedirection of the most positive eigenvalue of S , the β -term counteracts an unboundedgrowth because it decreases the most positive eigenvalue with a rate ∼ ω , even allowingit to become negative eventually. This goes along with a tilting of the eigenframe of S .The γ -term is a new term not present in the restricted Euler model. Considered fora fixed vorticity, it leads to a rotation of the eigenframe of S without changing itseigenvalues. This can be seen by the fact that a rotation of S with a rotation matrix M takes the form S (cid:48) = M S M T . For an infinitesimal time interval the rate-of-rotationtensor induces a rotation of the form M = I + W d t , such that the infinitesimal evolutionof S takes the form d S = − [ SW − WS ] d t (4.5)which is precisely the term associated with γ . Combined with a vorticity vector thatundergoes vortex stretching with the tendency to align with the direction of the mostpositive eigenvector, the γ -term induces an accelerated rotation of the eigenframe of S about this axis. The full system (4.1) shows a complex dynamics, whose statistics willbe discussed below.Finally, we study the magnitudes of the different terms and their Reynolds-numberdependence. To this end it is useful to non-dimensionalize (4.1). The velocity gradienttensor being a small-scale quantity motivates a non-dimensionalization with the Kol-mogorov time scale τ η = (cid:0) (cid:104) Tr (cid:0) S (cid:1) (cid:105) (cid:1) − / . With A (cid:63) = A τ η and t (cid:63) = t/τ η , (4.1) takesthe formd A (cid:63) = (cid:20) − (cid:18) A (cid:63) −
13 Tr (cid:0) A (cid:63) (cid:1) I (cid:19) − α (cid:18) S (cid:63) −
13 Tr (cid:0) S (cid:63) (cid:1) I (cid:19) − β (cid:18) W (cid:63) −
13 Tr (cid:0) W (cid:63) (cid:1) I (cid:19) − γ ( S (cid:63) W (cid:63) − W (cid:63) S (cid:63) ) + δ (cid:63) A (cid:63) (cid:21) d t (cid:63) + d F (cid:63) , (4.6)where δ (cid:63) = δτ η and d F (cid:63) = d F τ η . The non-dimensional coefficients α , β and γ remainunchanged by this procedure. As discussed above, α and β are Reynolds-number inde-pendent and the same is expected for γ under the simple universality argument.As will be derived from the Kolmogorov equation in section 4.3, we have δ (cid:63) = 7 S / (6 √ S is the derivative skewness coefficient (typically S ≈ − . -4.0-2.00.02.04.0 -4.0 -2.0 0.0 2.0 4.0 Q ⋆ R ⋆ (cid:16) S − Tr(S )I (cid:17) (cid:16) W − Tr(W )I (cid:17) -4.0 -2.0 0.0 2.0 4.0 R ⋆ (cid:16) S − Tr(S )I (cid:17) (cid:16) W − Tr(W )I (cid:17) Figure 1. (cid:0) S − Tr (cid:0) S (cid:1) I (cid:1) and (cid:0) W − Tr (cid:0) W (cid:1) I (cid:1) projected to the RQ -plane (see text fordetails on the projection). The strain term is predominantly active in the strain-dominatedregion Q (cid:63) < Q (cid:63) >
0. All conditional averages have been non-dimensionalized by (cid:10) Tr (cid:0) S (cid:1) (cid:11) . Here, and forthe following plots, vectors have been scaled with a factor of 0 . Shortcomings of the Gaussian closure
The evolution of the full system (4.1), in which all of the terms are simultaneously activeand completed by the (stochastic) forcing term, is governed by a combination of alleffects discussed above, which can lead to complex temporal behaviour. It is plausiblethat whether or not the system diverges depends crucially on the particular choice ofthe parameters α to δ . If we assume that the singularity of the restricted Euler modelshould be regularized by the non-local pressure Hessian contributions, α , β and γ are thecrucial parameters. As we have seen, α and β are universal, whereas γ depends on thetwo-point correlation of the velocity field.To further elucidate the question of stability of the closure, let us consider first the pro-jections of the different terms to the RQ -plane. Here and in the following, we use velocitygradient tensor fields from DNS to evaluate the terms of the closure, contract them with − A and − A , respectively, and non-dimensionalize them with proper powers of (cid:10) Tr (cid:0) S (cid:1)(cid:11) .Finally, they are binned with respect to Q (cid:63) = Q/ (cid:10) Tr (cid:0) S (cid:1)(cid:11) and R (cid:63) = R/ (cid:10) Tr (cid:0) S (cid:1)(cid:11) / .With this procedure, the contribution of the different terms to the reduced RQ -dynamicsis revealed, which helps to compare our closure, for example, to the restricted Euler model.It has to be stressed, though, that the dynamics of the closure involves all five possibleinvariants, i.e. information is lost due to projection. In this context, we would like to referthe reader to the recent work by L¨uthi et al. (2009) for a study of the velocity gradientdynamics beyond the RQ -plane.We analyse data from the JHU turbulence database (Li et al. R λ ≈
433 is publicly available. For thefollowing estimates a single 1024 snapshot of the database was used. However, we alsoperformed further checks on a number of snapshots to ensure that the reported resultsare robust. Velocity gradients and the Laplacian of the velocity gradients have been cal-4 M. Wilczek and C. Meneveau -4.0-2.00.02.04.0 Q ⋆ Gaussian closureDNS n o n l o c a l p r e ss u r e H e ss i a n s u m Gaussian closureDNS n o n l o c a l p r e ss u r e H e ss i a n s u m -4.0-2.00.02.04.0 -4.0 -2.0 0.0 2.0 4.0 Q ⋆ R ⋆ Gaussian closureDNS n o n l o c a l p r e ss u r e H e ss i a n s u m -4.0 -2.0 0.0 2.0 4.0 R ⋆ Gaussian closureDNS n o n l o c a l p r e ss u r e H e ss i a n s u m Figure 2.
Top panels: Non-local pressure Hessian contributions from DNS and the Gaussian clo-sure, projected to the RQ -plane. Bottom panels: Sum of the restricted Euler, non-local pressureHessian and viscous contributions from DNS and the Gaussian closure. culated spectrally from the downloaded velocity fields. Due to limited resolution of thesmallest scales, the velocity field was low-pass filtered for the evaluation of the velocitygradients and the Laplacian of the velocity gradients. The cut-off was determined suchthat the unphysical pile-up of kinetic energy at the highest wavenumbers (i.e. beyond ap-proximately 0 . k max , where k max denotes the highest dynamically active wavenumber)is excluded. As a consequence we expect inaccuracies especially for the viscous contribu-tions. The pressure Hessian has been evaluated spectrally from the downloaded pressurefields.To discriminate the effects of the terms related to α and β , figure 1 shows the pro-jections of (cid:0) S − Tr (cid:0) S (cid:1) I (cid:1) and (cid:0) W − Tr (cid:0) W (cid:1) I (cid:1) to the RQ -plane, respectively. The α -term turns out to be predominantly active in the strain-dominated half-plane ( Q (cid:63) < β -term. Interestingly, the term associated with γ doesnot contribute to this projection, because SW − WS vanishes when contracted with ar-bitrary powers of A .5The top panel of figure 2 compares the non-local pressure Hessian contribution fromDNS with the Gaussian closure. While the overall topology is quite similar, differencesin magnitude can be observed. As in the DNS, the closure attenuates the singularityinduced along the right part of the Vieillefosse line, however less strongly. This becomeseven more clear when comparing the sum of the contributions from the restricted Eu-ler term, the pressure Hessian and the viscous contributions from DNS to the Gaussianclosure, as shown in the bottom panel of figure 2. For the Gaussian closure, the totaleffect is apparently not strong enough to counteract the singularity induced along theright Vieillefosse tail. Because projection to the RQ -plane does not allow to comprehen-sively judge the performance of this closure, we also ran numerical tests solving (4.6)in the unforced, inviscid case with the parameter values given by the Gaussian closure.This consistently led to blow-up, i.e. the non-local pressure Hessian contributions in theGaussian closure are not able to prevent singularities by themselves. Also including thelinear damping term for the parameter range discussed below does not solve this problem.Tests show that the system still diverges for certain initial conditions, which is consistentwith the a priori findings reported in figure 2. This is because the linear damping termbecomes subdominant for sufficiently large velocity gradients compared to the quadraticterms of the restricted Euler and non-local pressure contributions (Meneveau 2011).We conclude that, while the structure of the Gaussian closure yields promising insights,it fails to predict coefficients that lead to a non-divergent evolution in time. A furthernumerical analysis of how the system diverges shows that especially the reduction ofnonlinearity associated with the α -term is not sufficient to prevent the singularity, whichcan be regarded as the main shortcoming of the Gaussian closure. The spatial structureof Gaussian random fields is insufficient to counteract the local self-amplification effectsby the restricted Euler part of the dynamics.4.3. Estimation of the coefficients from DNS data
To overcome the shortcomings of the Gaussian approximation, we accept the tenso-rial structure of the non-local pressure Hessian contributions as well as the conditionalLaplacian as obtained analytically from the Gaussian closure, but obtain estimates forthe parameter values from DNS data. We will refer to this as the enhanced Gaussianclosure.To estimate the parameters from the DNS data, where the conditional averages of ten-sor elements can be directly measured, we recognize that in general we have an overdeter-mined system if we wish to obtain scalar parameters. The most practical direct manneris to consider tensor contractions with appropriate powers of the velocity gradient tensorsuch that non-vanishing moments can be constructed.For the linear Laplacian term one may consider the expectation value of the contractionof the velocity gradient with its Laplacian, (cid:90) d A (cid:10) ν Tr (cid:0) A T ∆ A (cid:1) (cid:12)(cid:12) A (cid:11) f ( A ) = δ (cid:90) d A Tr (cid:16) A T A (cid:17) f ( A ) (4.7)which yields the relation δ = (cid:10) ν Tr (cid:0) A T ∆ A (cid:1) (cid:11)(cid:10) Tr ( A T A ) (cid:11) . (4.8)It can be shown by a straightforward calculation from the velocity gradient covariancestructure (B 1) that this estimate of δ is equivalent to the result (3.11) from the Gaussianapproximation. Evaluation of the averages in the numerator and denominator from theDNS data leads to the result δτ η ≈ − . M. Wilczek and C. Meneveau
To provide further theoretical insight, δ can be related to the velocity derivative skew-ness. The relevant derivation starts from the Kolmogorov equation (Monin et al. (cid:68) ( u ( x + r e ) − u ( x )) (cid:69) − ν dd r (cid:68) ( u ( x + r e ) − u ( x )) (cid:69) = − εr . (4.9)An expansion of this expression about r = 0 yields (cid:42)(cid:18) ∂u ∂x (cid:19) (cid:43) r − ν (cid:42)(cid:18) ∂u ∂x (cid:19) (cid:43) r + 2 ν (cid:42)(cid:18) ∂ u ∂x (cid:19) (cid:43) r + h . o . t . = − εr . (4.10)Comparing the terms of order one yields the usual relation ε = 15 ν (cid:68) ( ∂u /∂x ) (cid:69) . Theterms of third order give (cid:42)(cid:18) ∂u ∂x (cid:19) (cid:43) = − ν (cid:42)(cid:18) ∂ u ∂x (cid:19) (cid:43) , (4.11)which can be used to relate the fourth-order derivative in (3.11) to the velocity gradientskewness S = (cid:10) ( ∂u /∂x ) (cid:11) / (cid:10) ( ∂u /∂x ) (cid:11) / . With (cid:68) ( ∂u /∂x ) (cid:69) = − ( (cid:104) u (cid:105) / f (cid:48)(cid:48) u (0)and (cid:68)(cid:0) ∂ u /∂x (cid:1) (cid:69) = ( (cid:104) u (cid:105) / f (4) u (0) the coefficient δ for isotropic turbulence obeyingthe Kolmogorov equation takes the form δ = 76 (cid:42)(cid:18) ∂u ∂x (cid:19) (cid:43)(cid:30)(cid:42)(cid:18) ∂u ∂x (cid:19) (cid:43) = 76 √ S τ η . (4.12)As a result the conditional Laplacian depends on the Reynolds number mainly through τ η , but also due to a weak dependence of S . The derivative skewness is known to havevalues near S ≈ − .
5. Recent DNS results by Ishihara et al. (2007) at R λ = 471, forexample, provide a value of S between − . − .
5, which leads to a value of δτ η between − .
15 and − .
18, which is in good agreement with our DNS findings. For furtheranalysis, the value of δτ η ≈ − .
15 will be assumed.To estimate the coefficients of the non-local pressure Hessian contribution (3.12) fromDNS, we obtain expectation values analogous to (4.7). The lowest-order contraction thatleads to a non-vanishing of the γ -term is the contraction with the quadratic expression SW . For consistency, we also use contractions with quadratic expressions of the form S and W to estimate α and β , which leads to a linear set of equations. This procedureis also interesting from a physical point of view, because it connects the parameterestimates for the non-local pressure Hessian with the flatness factors of the strain fieldand the vorticity field. The set of equations is readily solved by α = (cid:10) Tr (cid:0) W (cid:103) W (cid:1)(cid:11)(cid:10) Tr (cid:0) S (cid:101) H (cid:1)(cid:11) − (cid:10) Tr (cid:0) S (cid:103) W (cid:1)(cid:11)(cid:10) Tr (cid:0) W (cid:101) H (cid:1)(cid:11)(cid:10) Tr (cid:0) S (cid:102) S (cid:1)(cid:11)(cid:10) Tr (cid:0) W (cid:103) W (cid:1)(cid:11) − (cid:10) Tr (cid:0) S (cid:103) W (cid:1)(cid:11)(cid:10) Tr (cid:0) W (cid:102) S (cid:1)(cid:11) (4.13) β = (cid:10) Tr (cid:0) S (cid:102) S (cid:1)(cid:11)(cid:10) Tr (cid:0) W (cid:101) H (cid:1)(cid:11) − (cid:10) Tr (cid:0) W (cid:102) S (cid:1)(cid:11)(cid:10) Tr (cid:0) S (cid:101) H (cid:1)(cid:11)(cid:10) Tr (cid:0) S (cid:102) S (cid:1)(cid:11)(cid:10) Tr (cid:0) W (cid:103) W (cid:1)(cid:11) − (cid:10) Tr (cid:0) S (cid:103) W (cid:1)(cid:11)(cid:10) Tr (cid:0) W (cid:102) S (cid:1)(cid:11) (4.14) γ = (cid:10) Tr (cid:0) SW (cid:101) H (cid:1)(cid:11)(cid:10) Tr (cid:0) SW [ SW − WS ] (cid:1)(cid:11) . (4.15)We have used the short-hand notation (cid:103) W = W − Tr (cid:0) W (cid:1) I and (cid:102) S = S − Tr (cid:0) S (cid:1) I for denoting the traceless tensors. This means that the parameters α , β and γ can be7fixed by evaluating averages involving the rate-of-strain and rate-of-rotation tensors aswell as the pressure Hessian.From the DNS data set we obtain α ≈ − . β ≈ − .
65, and γ ≈ .
14, i.e. all valuesare larger in magnitude than predicted by the Gaussian approximation. In particular, thishas the consequence that the nonlinearity associated with the α -term is weaker comparedto the restricted Euler model and the Gaussian closure. As we will see below, this willallow for a non-divergent evolution of the system.Our tests with different snapshots of the DNS data indicate that the reported values areaccurate within a few per cent. They may, however, show a Reynolds-number dependence,whose study goes beyond the scope of the current work. We close this section withreiterating that the only difference between the Gaussian closure and its enhanced versionis in the choice of the constant coefficients of the non-local pressure Hessian.4.4. Comparison of enhanced Gaussian closure and DNS results
Figure 3 shows a comparison of the different terms of the velocity gradient dynamicsestimated from the JHU DNS database along with the terms of the enhanced Gaussianclosure closure projected to the RQ -plane. As the restricted Euler contribution appearsclosed in this description, the closure coincides with the DNS data.The non-local pressure Hessian contribution of the enhanced Gaussian closure sharesa number of qualitative similarities with the term from DNS data. Most importantly,and in contrast to the original Gaussian closure, it fully counteracts the singularity alongthe right branch of the Vieillefosse line. It also reproduces the tendency of the pressureHessian to push towards negative R (cid:63) for positive Q (cid:63) . The main difference is observed forpositive R (cid:63) and moderate values of Q (cid:63) , where the closure does not reproduce the verysmall amplitudes of the non-local pressure Hessian from DNS.The diffusive term of the closure compares quite well to the diffusive term from DNSin its overall damping influence, although minor differences in direction and amplitudecan be observed. We note that the enhanced Gaussian closure coincides with the originalone in this term.The total vector field of the enhanced closure compares qualitatively well with the DNSresults; it indicates a clockwise cyclic motion in the RQ -plane. However, the enhancedclosure seems to push towards the right Vieillefosse tail too early, such that the PDF of R (cid:63) and Q (cid:63) will show discrepancies in this region.4.5. Enhanced Gaussian closure in SDE
Next we investigate the statistics of the model by using it in the context of a SDE. This isdenoted as the ‘enhanced Gaussian closure SDE’ model and is based on (4.6). Althoughrooted in a Gaussian assumption for the unclosed term, this SDE will naturally pro-duce non-Gaussian statistics due to the joint nonlinear effects of local self-amplificationand non-local pressure Hessian contributions. The SDE of the enhanced Gaussian clo-sure is implemented with a fourth-order Runge-Kutta method for the deterministic part,combined with an Euler method for the noise term. The details of the noise term imple-mentation are given in appendix D.For the numerical solution we choose the values obtained from the DNS data, α = − . β = − . γ = 0 .
14 and δ (cid:63) = − .
15. To determine the amplitude of the forceterm, A F (cf. appendix D), we note that the model should fulfil (cid:10) Tr (cid:0) S (cid:63) (cid:1)(cid:11) = 1 /
2, which isan immediate consequence of the definition of the mean rate of kinetic energy dissipationand the non-dimensionalization leading to (4.6). We ran a number of numerical tests tochoose the forcing amplitude such that this constraint is fulfilled in good approximation;figure 4 shows (cid:10) Tr (cid:0) S (cid:63) (cid:1)(cid:11) as a function of the forcing amplitude A F . The resulting plot8 M. Wilczek and C. Meneveau -4.0-2.00.02.04.0 Q ⋆ DNS enhanced Gaussian closure r e s t r i c t e d E u l e r n o n l o c a l p r e ss u r e H e ss i a nd i ff u s i o n s u m DNS enhanced Gaussian closure r e s t r i c t e d E u l e r n o n l o c a l p r e ss u r e H e ss i a nd i ff u s i o n s u m -4.0-2.00.02.04.0 Q ⋆ DNS enhanced Gaussian closure r e s t r i c t e d E u l e r n o n l o c a l p r e ss u r e H e ss i a nd i ff u s i o n s u m DNS enhanced Gaussian closure r e s t r i c t e d E u l e r n o n l o c a l p r e ss u r e H e ss i a nd i ff u s i o n s u m -4.0-2.00.02.04.0 Q ⋆ DNS enhanced Gaussian closure r e s t r i c t e d E u l e r n o n l o c a l p r e ss u r e H e ss i a nd i ff u s i o n s u m DNS enhanced Gaussian closure r e s t r i c t e d E u l e r n o n l o c a l p r e ss u r e H e ss i a nd i ff u s i o n s u m -4.0-2.00.02.04.0 -4.0 -2.0 0.0 2.0 4.0 Q ⋆ R ⋆ DNS enhanced Gaussian closure r e s t r i c t e d E u l e r n o n l o c a l p r e ss u r e H e ss i a nd i ff u s i o n s u m -4.0 -2.0 0.0 2.0 4.0 R ⋆ DNS enhanced Gaussian closure r e s t r i c t e d E u l e r n o n l o c a l p r e ss u r e H e ss i a nd i ff u s i o n s u m Figure 3.
Projection of the terms governing the RQ -evolution from DNS data and theenhanced Gaussian closure. − − − h T r ( S ⋆ ) i A F Figure 4. (cid:10) Tr (cid:0) S (cid:63) (cid:1)(cid:11) as a function of the forcing amplitude A F for the enhanced Gaussianclosure SDE. The dashed blue line indicates the constraint (cid:10) Tr (cid:0) S (cid:63) (cid:1)(cid:11) = 1 /
2, which is used todetermine the amplitude.
Figure 5.
Joint PDFs of R (cid:63) and Q (cid:63) . Left: DNS data, right: SDE model. While the main featuresare captured accurately by the model, it overestimates the rotational and underestimates thestraining regions in the RQ -plane. evidences a strong dependence of the strain-rate variance resulting from the model asfunction of the forcing strength. Still, clearly it can be seen that choosing A F = 0 . (cid:10) Tr (cid:0) S (cid:63) (cid:1)(cid:11) ≈ /
2. Therefore, in the numerical solutions presented below, we set A F = 0 .
13. For the numerical integration, the time step is set to 10 − . For the statisticalevaluation, we draw 10 initial conditions from a Gaussian ensemble (that also fulfills (cid:10) Tr (cid:0) S (cid:63) (cid:1)(cid:11) ≈ /
2) and let them evolve for 5 × time steps. After this initial transient,we evolve the SDE for another 5 × time steps during which statistics are gathered.We have checked that all presented results are statistically well converged. However, wealso noted that higher- (e.g. fourth-)order moments did not fully converge, possibly dueto some rare trajectories visiting far-out regions of sample space.Figure 5 shows the PDFs of R and Q both from DNS and the enhanced Gaussian0 M. Wilczek and C. Meneveau P D F cosine of anglemost negative eigenvalueintermediate eigenvaluemost positive eigenvalue 0.01.02.03.04.05.00.0 0.2 0.4 0.6 0.8 1.0 P D F cosine of anglemost negative eigenvalueintermediate eigenvaluemost positive eigenvalue Figure 6.
PDFs of angle cosines (cid:98) ω · s i , where s i denotes the eigenvectors of the rate-of-straintensor associated to the three eigenvalues. Left: DNS data, right: SDE model. The geometrictrends of the DNS results are captured accurately by the model, although differences in theamplitudes of the PDFs can be seen. closure SDE. The PDFs share qualitative similarities like the characteristic tear-dropshape related to the non-vanishing enstrophy production and strain skewness. The modelPDF, however, overestimates the rotational regions in the RQ -plane, while the strainingregions, especially around the left branch of the Vieillefosse line, are underestimated. Thisresults in the fact that the model does not respect the Betchov relations (Betchov 1956) (cid:10) Tr (cid:0) S (cid:1)(cid:11) = (cid:10) ω (cid:11) / − (cid:10) Tr (cid:0) S (cid:1)(cid:11) = 3 (cid:104) ω i S ij ω j (cid:105) /
4. This discrepancy is plausiblefrom the model pressure Hessian, which tends to suppress the statistical evolution alongthe right Vieillefosse tail, as well as the linearity of the diffusive term, which does notcapture the details observed in the DNS. We also would like to note that Betchov’srelations represent constraints of averages over fields, which are inherently difficult toincorporate into single-point models.Focusing on geometric features of the enhanced Gaussian closure SDE, the alignmentPDFs of the vorticity vector with the eigenvectors of the rate-of-strain tensor are pre-sented in figure 6. Results are shown for both the DNS data as well as the enhancedGaussian closure SDE. As can be inferred from this figure, the geometrical trends ob-served in the DNS data are captured quite accurately by the model. While the qualitativebehaviour of the alignment with the most negative and intermediate eigendirections arecaptured satisfactorily, differences in amplitude of the PDF are observed.The topology of the velocity gradient tensor can be conveniently described using theparameter s (cid:63) = −√ (cid:0) S (cid:1) / Tr (cid:0) S (cid:1) / , introduced by Lund & Rogers (1994). The lim-iting case s (cid:63) = − s (cid:63) = 1corresponds to axisymmetric expansion. Figure 7 shows its PDF for the DNS data andthe enhanced Gaussian closure SDE results. For the DNS the PDF is a strictly increasingfunction showing that the preferred state of strain is axisymmetric expansion (Lund &Rogers 1994). This feature is qualitatively reproduced by the model, it however overes-timates states of axisymmetric contraction and underestimates states of axisymmetricexpansion somewhat. For a purely Gaussian field (without strain skewness, etc.), thecorresponding PDF is flat.These comparisons show that the enhanced Gaussian closure is capable of reproducingsome important qualitative features of homogeneous isotropic turbulence. For betterquantitative agreement, however, more accurate modelling of the unclosed terms will benecessary. Specifically, the coefficients α to δ need to be generalized to depend on the1 P D F s ⋆ DNSenhanced Gaussian closurereversible enhanced Gaussian closureGaussian
Figure 7.
PDFs of s (cid:63) = −√ (cid:0) S (cid:1) / Tr (cid:0) S (cid:1) / for DNS, the enhanced Gaussian closure aswell as its reversible version. The flat PDF of a Gaussian random field is shown for reference. P D F cosine of anglemost negative eigenvalueintermediate eigenvaluemost positive eigenvalue Figure 8.
Left: RQ -PDF for the time-reversible enhanced Gaussian closure ODE. The RQ -PDFis symmetric with respect to R → − R which implies vanishing strain skewness and enstrophyproduction. Right: alignment PDFs. Compared with the irreversible model (see figure 6), thealignment properties change significantly. invariants of the velocity gradient tensor. This, however, falls beyond the scope of thepresent work. 4.6. Time-reversal symmetry
One notable feature of the evolution equation (4.1) is the time-reversal symmetry inthe undamped ( δ = 0) and unforced (d F = 0) case. We have seen that the non-localpressure Hessian contributions are able to prevent a singularity for certain parameterchoices, thus leading to stationary statistics. Considering the undamped and unforcedcase, this has the interesting consequence that all odd-order moments of the velocitygradient tensor vanish, as these quantities change sign under time reversal. This impliesthat the restricted Euler part of the dynamics combined with the stabilizing non-localpressure Hessian contributions considered on its own yields vanishing strain skewness,2 M. Wilczek and C. Meneveau (cid:10) Tr (cid:0) S (cid:1) (cid:11) = 0, as well as vanishing enstrophy production, (cid:10) ω i S ij ω j (cid:11) = 0. We note inpassing that, although the restricted Euler model displays the time-reversal symmetry,too, the argument does not hold in this case: The restricted Euler model does not producestationary statistics.These considerations can be confirmed with numerical results from the enhanced Gaus-sian closure ordinary differential equation (ODE; since it has no forcing, it reduces to anODE). The PDF of Q and R for this case is displayed in figure 8. As expected, the PDF isfully symmetric with respect to the transformation R → − R . Along with the RQ -PDF,also the PDF of vorticity alignment with the principal strain axes is shown. Interestingly,the time-reversible ODE has dramatically different alignment properties. For example,the PDFs of vorticity alignment with the most negative and most positive strain-rateeigendirection now collapse. This can be understood from the fact that the rate-of-straintensor eigenvalues change sign under time reversal, i.e. the most positive and the mostnegative one change their role. As the statistics has to be invariant with respect to thattransformation, the alignment PDFs of the vorticity with the most positive and mostnegative eigenvalue then have to collapse.Also the PDF of s (cid:63) has been evaluated for the reversible case (see figure 7). It isfully symmetric with respect to the transformation s (cid:63) → − s (cid:63) , which also can be directlyinferred by the behaviour of this quantity under time reversal.Including linear damping and stochastic forcing, which both break this symmetry, thenproduces skewed statistics, non-vanishing enstrophy production and the familiar strainand alignment properties. These effects occur due to the combination of time-reversalsymmetry preserving and breaking terms. Consequently, some of the essential statisticalproperties of small-scale turbulence cannot exclusively be associated to the closure ofthe pressure Hessian alone, but depend on its interplay with the dissipative and energy-injecting terms.
5. Conclusions
We have evaluated the effects of non-local pressure Hessian contributions and viscousdiffusion in the framework of a statistical evolution equation for the velocity gradienttensor under the assumption of Gaussian incompressible velocity fields.In this scenario, the viscous term is obtained as a linear damping term, where thecoefficient is Reynolds-number dependent through its dependence on the velocity auto-correlation function.The non-local contributions to the pressure Hessian are found to be a combinationof quadratic, traceless and symmetric expressions of the rate-of-strain and the rate-of-rotation tensors. Two of these terms modify the original restricted Euler model withcoefficients that are independent of the Reynolds number. In addition, the Gaussianclosure yields a term which induces a rotation of the eigenframe of the rate-of-straintensor, whose coefficient depends on the details of the velocity two-point correlations.The simplicity of the Gaussian closure allowed to discuss the different dynamical ef-fects like stretching and tilting of the vorticity vector and rate-of-strain eigenvalues andeigenvectors, revealing how the various non-local contributions of the pressure Hessianhelp to attenuate the occurrence of the singularity.The closed dynamical system has then been investigated numerically, showing thatthe coefficients obtained in the Gaussian approximation are insufficient to prevent thesingularity caused by the restricted Euler part. Maintaining the overall structure of theGaussian approximation, physically more realistic values for coefficients were obtainedusing fitting to DNS data, in a mean-field approach. The pressure Hessian and diffu-3sive term in the enhanced Gaussian closure were shown to yield improved qualitativeagreement with the DNS results.The same applies to the statistical properties of the enhanced Gaussian closure SDE,which qualitatively captures main features such as strain skewness, enstrophy produc-tion and alignment of the vorticity with the strain eigenvectors. While the enhancedclosure leads to a non-divergent time evolution, some quantitative differences to DNSwere observed. A study of the closure in the time-reversible case of undamped, unforceddynamics showed some interesting results that pointed to the subtle interactions betweennon-local pressure contributions and the dissipative term that must be taking place inreal turbulence.With respect to possible generalizations it is interesting to note that (3.12) alreadyrepresents the most general structure of the non-local pressure Hessian, which is sym-metric, traceless and dimensionally consistent. In general, however, the coefficients maydepend on non-dimensional combinations of the five invariants of the velocity gradienttensor. In particular, the discrepancies of the closures compared with the DNS resultsmake coefficients that depend on, for example, strain skewness and enstrophy productionplausible, which will be the topic of future research.
Acknowledgments
We thank Gregory L. Eyink, Laurent Chevillard, Anton Daitche and Perry Johnsonfor insightful discussions. Cristian C. Lalescu is gratefully acknowledged for help with theJHU database. MW acknowledges support from the DFG under project WI 3544/2-1.CM is thankful for funding from the National Science Foundation (grants CBET 1033942and CDI, CMMI-0941530).
Appendix A. PDF equation for the velocity gradient
A.1.
Derivation of the deterministic part of the PDF equation
This appendix is devoted to derive the evolution equation for the velocity gradient PDFand discuss the different formulations of the closure problem. As usual, the PDF isconveniently introduced as an ensemble average over a fine-grained PDF (Lundgren 1967;Pope 2000) applied to the velocity gradient: f ( A ; x , t ) = (cid:10) δ ( A ( x , t ) − A ) (cid:11) . (A 1)Here, A denotes the sample-space variable corresponding to the velocity gradient tensorfield at position x and time t , and f is a probability density with respect to A and afunction with respect to x and t . We first consider the unforced case and postpone therather technical discussion of how to include a stochastic large-scale forcing. To obtainan evolution equation for the probability density, we take the partial derivative of (A 1)with respect to time and use the chain rule and a simple change of variables: ∂∂t f ( A ; x , t ) = (cid:28) ∂∂t δ ( A ( x , t ) − A ) (cid:29) (A 2 a )= − ∂∂ A ij (cid:28)(cid:20) ∂∂t A ij ( x , t ) (cid:21) δ ( A ( x , t ) − A ) (cid:29) . (A 2 b )Here and throughout the paper summation over double indices is implied. For the timederivative of the velocity gradient field we now can substitute the dynamical equation4 M. Wilczek and C. Meneveau (2.1) and obtain ∂∂t f ( A ; x , t ) = − ∂∂ A ij (cid:28)(cid:20) ∂∂t A ij ( x , t ) (cid:21) δ ( A ( x , t ) − A ) (cid:29) (A 3 a )= − ∂∂ A ij (cid:28)(cid:20) − u k ∂∂x k A ij − (cid:18) A ik A kj − A lk A kl δ ij (cid:19) − (cid:101) H ij + ν ∆ A ij (cid:21) δ ( A − A ) (cid:29) . (A 3 b )For brevity we have omitted the space-time dependence of all fields in the second line.This expression already shows in what way the closure problem enters this framework,namely in terms of joint averages of the various fields with the fine-grained distribution.The advective term can be treated further yielding − ∂∂ A ij (cid:28)(cid:20) − u k ∂∂x k A ij (cid:21) δ ( A − A ) (cid:29) = − ∂∂x k (cid:10) u k δ ( A − A ) (cid:11) . (A 4)Here we have used again the chain rule and incompressibility of the velocity field. Forhomogenous flow the average of the velocity field and the delta distribution is independentof the spatial variable and hence the derivative vanishes. That means for homogeneousturbulence the advective term is gone.Next we consider the term which is local in A . Using the sifting property of the deltadistribution, A ( x , t ) δ ( A ( x , t ) − A ) = A δ ( A ( x , t ) − A ), we obtain (cid:28)(cid:20) A ik ( x , t ) A kj ( x , t ) − A lk ( x , t ) A kl ( x , t ) δ ij (cid:21) δ ( A ( x , t ) − A ) (cid:29) = (cid:20) A ik A kj − A lk A kl δ ij (cid:21) f ( A ; x , t ) . (A 5)Again we have made use of the fact that sample-space variables can be pulled out of theaverage. The main point here is that the self-amplification term along with the isotropicpart of the pressure Hessian are closed and do not need to be modeled.For the pressure Hessian and the viscous term two options are available, namely ex-pressing them with the help of conditional averages or in terms of multipoint statistics.We start with the former option which leads to (cid:10) (cid:101) H ij ( x , t ) δ ( A ( x , t ) − A ) (cid:11) = (cid:10) (cid:101) H ij ( x , t ) (cid:12)(cid:12) A (cid:11) f ( A ; x , t ) (A 6) (cid:10) [ ν ∆ A ( x , t )] δ ( A ( x , t ) − A ) (cid:11) = (cid:10) ν ∆ A ( x , t ) (cid:12)(cid:12) A (cid:11) f ( A ; x , t ) . (A 7)By introduction of conditional averages the PDF can be isolated at the price of intro-ducing unknown functions. To close the expressions tensorial functions depending on atensorial argument have to be modeled. The relation to the two-point statistics will bediscussed later below. A.2. Inclusion of a stochastic force
To conclude the derivation of the PDF equation we consider the forcing term. If a de-terministic forcing term is considered, it also can be treated with the help of conditionalaverages. We here, however, consider a stochastic large-scale forcing, which eventuallyallows us to make direct contact with earlier phenomenological models. If the Reynoldsnumber is sufficiently high, the velocity gradient statistics should be independent of theparticular choice of large-scale forcing, making the choice of this forcing not a particularlystrong restriction.The inclusion of a stochastic forcing term turns the Navier-Stokes equation to a stochas-tic partial differential equation. The forcing term is specified as a homogeneous isotropic5Gaussian white-in-time random force F u ( x , t ) added as an additional acceleration termwith the properties (cid:104) F u ( x , t ) (cid:105) = and (cid:104) F ui ( x , t ) F uj ( x (cid:48) , t (cid:48) ) (cid:105) = Q uij ( r ) δ ( t − t (cid:48) ) (A 8)where r = x − x (cid:48) . Isotropy and solenoidality require the structure of its covariance to be Q uij ( r ) = σ F (cid:20) f F ( r ) δ ij + 12 rf (cid:48) F ( r ) [ δ ij − ˆ r i ˆ r j ] (cid:21) (A 9)where σ F denotes the forcing variance and f F denotes the longitudinal forcing autocor-relation function.We are interested in how the stochastic forcing for the velocity field translates to thevelocity gradient field. For the velocity gradient dynamics the gradient of this randomforce has to be added, as indicated in (2.1), which then has the properties (cid:104) F ij ( x , t ) (cid:105) = 0 and (cid:104) F ik ( x , t ) F jl ( x (cid:48) , t (cid:48) ) (cid:105) = Q ijkl ( r ) δ ( t − t (cid:48) ) . (A 10)The relations specifying the structure of the forcing resemble the results on the generalstructure of the velocity gradient tensor covariance discussed below in appendix B.1, sothat many of the technical results can be used here. The relation between the covariancesof the random forcing for the velocity field and the velocity gradient field analogouslyreads Q ijkl ( r ) = − ∂∂r k ∂∂r l Q uij ( r ) . (A 11)For the following we will especially need the case of r = , for which Q ijkl ( ) = − σ F f (cid:48)(cid:48) F (0) (cid:20) δ ij δ kl − δ ik δ jl − δ il δ jk (cid:21) . (A 12)To learn how the stochastic forcing carries through to our statistical description, weadapt the presentation by Haken (2004) for the derivation of the Fokker-Planck equationstarting from the Langevin equation and generalize it to partial differential equations.Specifically, we consider how the PDF evolves over a short time interval ∆ t . To this endwe introduce the short-hand notation∆ f = f ( A ; x , t + ∆ t ) − f ( A ; x , t ) and (A 13)∆ A = A ( x , t + ∆ t ) − A ( x , t ) (A 14)and consider furthermore∆ A = (cid:2) − u · ∇ x A − A − H + ν ∆ x A (cid:3) ∆ t + ∆ F + O (cid:0) ∆ t (cid:1) . (A 15)Note that, for the present discussion, ∆ refers to an increment rather than the Lapla-cian, which here is denoted as ∆ x to avoid ambiguities. The notation ∆ F for the forceincrement indicates that we are dealing with a stochastic force which is not differen-tiable in time. For the case that ∆ t is infinitesimally small, this resembles the notationof stochastic calculus, but for the moment considering a small, but finite ∆ t is sufficient.Expanding the evolution of the PDF up to second order we obtain∆ f = − ∂∂ A ij (cid:10) δ ( A − A ) ∆ A ij (cid:11) + 12 ∂∂ A ik ∂∂ A jl (cid:10) δ ( A − A ) ∆ A ik ∆ A jl (cid:11) + O (cid:0) ∆ A (cid:1) . (A 16)Next, (A 15) is inserted into this expression, and the individual terms are evaluatedwith respect to their dependence on the time increment ∆ t . Only terms linear in ∆ t areof interest, because we finally want to evaluate lim ∆ t → ∆ f / ∆ t .6 M. Wilczek and C. Meneveau
The first term in (A 16) corresponds to the deterministic contributions already dis-cussed, plus a term involving the joint average of the force increment with the deltadistribution. Because the force increment contains only contributions arising after thetime t , the fine-grained PDF and the force increment are statistically independent, thus (cid:10) δ ( A − A ) ∆ F ij (cid:11) = (cid:10) δ ( A − A ) (cid:11)(cid:10) ∆ F ij (cid:11) = 0 . (A 17)The last equality is due to the fact that we are considering stochastic forces with zeromean. Consequently the forcing term gives no first-order contribution. This is also thereason for why we could first consider the unforced case without loss of generality.Owing to the quadratic dependence of the second term on ∆ A , it contains contributionsproportional to ∆ t , ∆ t ∆ F and ∆ F . The terms proportional to ∆ t vanish in thelimit eventually taken, and so do the terms proportional to ∆ t ∆ F because of the aboveargument of statistical independence and vanishing mean forcing. The only remainingterm can be treated according to (cid:10) δ ( A − A ) ∆ F ik ∆ F jl (cid:11) = (cid:10) δ ( A − A ) (cid:11)(cid:10) ∆ F ik ∆ F jl (cid:11) (A 18 a )= (cid:10) δ ( A − A ) (cid:11) Q ijkl ( )∆ t + O (cid:0) ∆ t (cid:1) . (A 18 b )For the first equality we have again used the argument of statistical independence. Tosee that (cid:10) ∆ F ik ∆ F jl (cid:11) is linear in ∆ t , one needs to recall that the force is specified asdelta-correlated in time which cancels one of the integrations necessary to evaluate thefinite increment.By evaluating the limit lim ∆ t → ∆ f / ∆ t and combining the results of this and thepreceding section, we arrive at the PDF equation for the velocity gradient tensor inhomogeneous turbulence: ∂∂t f ( A ; t ) = − ∂∂ A ij (cid:18)(cid:20) − (cid:18) A ik A kj −
13 Tr (cid:0) A (cid:1) δ ij (cid:19) − (cid:10) (cid:101) H ij (cid:12)(cid:12) A (cid:11) + (cid:10) ν ∆ x A ij (cid:12)(cid:12) A (cid:11)(cid:21) f ( A ; t ) (cid:19) + 12 Q ijkl ( ) ∂∂ A ik ∂∂ A jl f ( A ; t ) . (A 19)A.3. The relation to multipoint statistics
Instead of treating the unclosed terms arising in the derivation of the PDF equation withthe help of conditional averages as done in (A 6) and (A 7), they can also be expressedin terms of two-point statistics. In the following subscripts on the position vectors andsample-space variables will be used to discriminate the two spatial points, and f and f will be used for the one- and two-point PDFs, respectively. We can evaluate the viscousterm according to (Lundgren 1967) (cid:10) [ ν ∆ x A ( x , t )] δ ( A ( x , t ) − A ) (cid:11) = lim | x − x |→ (cid:10) [ ν ∆ x A ( x , t )] δ ( A ( x , t ) − A ) (cid:11) (A 20 a )= lim | x − x |→ ν ∆ x (cid:10) A ( x , t ) (cid:12)(cid:12) A (cid:11) f ( A ; x , t ) . (A 20 b )Together with (A 7) this leads to the result (cid:10) ν ∆ x A ( x , t ) (cid:12)(cid:12) A (cid:11) = lim | x − x |→ ν ∆ x (cid:10) A ( x , t ) (cid:12)(cid:12) A (cid:11) . (A 21)For homogeneous turbulence, the conditional average is a function of the distance vector r = x − x only, such that we can also write (cid:10) ν ∆ x A ( x , t ) (cid:12)(cid:12) A (cid:11) = lim r → ν ∆ r (cid:10) A ( x , t ) (cid:12)(cid:12) A (cid:11) . (A 22)7For the non-local contributions to the pressure Hessian we have to make use of the Poissonrelation (2.4) and also of the identity (cid:90) d A δ ( A ( x , t ) − A ) = 1 . (A 23)First we consider (cid:10) Tr (cid:0) A ( x , t ) (cid:1) δ ( A ( x , t ) − A ) (cid:11) = (cid:90) d A (cid:10) Tr (cid:0) A ( x , t ) (cid:1) δ ( A ( x , t ) − A ) δ ( A ( x , t ) − A ) (cid:11) (A 24 a )= (cid:90) d A Tr (cid:0) A (cid:1) f ( A , A ; x , x , t ) (A 24 b )= (cid:10) Tr (cid:0) A ( x , t ) (cid:1) (cid:12)(cid:12) A (cid:11) f ( A ; x , t ) . (A 24 c )With this relation the non-local pressure Hessian can be expressed as (cid:10) (cid:101) H ij ( x , t ) δ ( A ( x , t ) − A ) (cid:11) = − π (cid:90) P . V . d x (cid:20) δ ij | x − x | − x − x ) i ( x − x ) j | x − x | (cid:21) (cid:10) Tr (cid:0) A ( x , t ) (cid:1) δ ( A ( x , t ) − A ) (cid:11) (A 25 a )= − π (cid:90) P . V . d x (cid:20) δ ij | x − x | − x − x ) i ( x − x ) j | x − x | (cid:21) (cid:10) Tr (cid:0) A ( x , t ) (cid:1) (cid:12)(cid:12) A (cid:11) f ( A ; x , t ) . (A 25 b )Together with (A 6) this leads to the result (cid:10) (cid:101) H ij ( x , t ) (cid:12)(cid:12) A (cid:11) = 12 π (cid:90) P . V . d x (cid:20) δ ij | x − x | − x − x ) i ( x − x ) j | x − x | (cid:21) (cid:10) Q ( x , t ) (cid:12)(cid:12) A (cid:11) , (A 26)and for homogeneous turbulence (cid:10) (cid:101) H ij ( x , t ) (cid:12)(cid:12) A (cid:11) = 12 π (cid:90) P . V . d r (cid:20) δ ij r − r i r j r (cid:21) (cid:10) Q ( x , t ) (cid:12)(cid:12) A (cid:11) . (A 27) Appendix B. Gaussian approximation
B.1.
Velocity gradient covariance tensor
The general structure of the velocity gradient covariance tensor is obtained by evaluatingthe kinematic relation (3.7). We also make use of homogeneity, which implies that thecovariance tensor depends on r = x − x (cid:48) only, and isotropy. As a result we obtain R ijkl ( r ) = a δ ij δ kl + a (cid:2) δ ik δ jl + δ il δ jk (cid:3) + a δ ij ˆ r k ˆ r l + a (cid:2) δ ik ˆ r j ˆ r l + δ il ˆ r j ˆ r k + δ jk ˆ r i ˆ r l + δ jl ˆ r i ˆ r k + δ kl ˆ r i ˆ r j (cid:3) + a ˆ r i ˆ r j ˆ r k ˆ r l (B 1)8 M. Wilczek and C. Meneveau where the coefficients depend on the longitudinal velocity autocorrelation function andare given by a = (cid:104) u (cid:105) (cid:20) − f (cid:48) u r − f (cid:48)(cid:48) u (cid:21) (B 2 a ) a = (cid:104) u (cid:105) (cid:20) f (cid:48) u r (cid:21) (B 2 b ) a = (cid:104) u (cid:105) (cid:20) f (cid:48) u r − f (cid:48)(cid:48) u − rf (cid:48)(cid:48)(cid:48) u (cid:21) (B 2 c ) a = (cid:104) u (cid:105) (cid:20) − f (cid:48) u r + f (cid:48)(cid:48) u (cid:21) (B 2 d ) a = (cid:104) u (cid:105) (cid:20) f (cid:48) u r − f (cid:48)(cid:48) u + rf (cid:48)(cid:48)(cid:48) u (cid:21) . (B 2 e )Some interesting observations can be made here. Incompressibility of the velocity fieldimplies A ii = 0. For the covariance tensor this implies R ijil = 0 and R ijkj = 0 which isreadily checked with the above results. We also need to explicitly evaluate the tensor for r = , in which case only the coefficients a and a should remain. Indeed, by makinguse of L’Hospital’s rule we obtain a = a = a = 0 and a (0) = lim r → (cid:104) u (cid:105) (cid:20) − f (cid:48) u r − f (cid:48)(cid:48) u (cid:21) = − (cid:104) u (cid:105) f (cid:48)(cid:48) u (0) = 215 εν (B 3 a ) a (0) = lim r → (cid:104) u (cid:105) (cid:20) f (cid:48) u r (cid:21) = (cid:104) u (cid:105) f (cid:48)(cid:48) u (0) = − εν . (B 3 b )For these calculations we have made use of the properties f (cid:48) u (0) = 0 and the relation(see, e.g., Pope (2000)) (cid:104) u (cid:105) f (cid:48)(cid:48) u (0) = − εν = − (cid:10) Tr (cid:0) S (cid:1)(cid:11) . (B 4)This leads to the result R ijkl ( ) = (cid:104) u (cid:105) f (cid:48)(cid:48) u (0) (cid:2) − δ ij δ kl + δ ik δ jl + δ il δ jk (cid:3) . (B 5)For the Fourier transform of the single-point characteristic function to the single-pointPDF of the velocity gradient tensor also the inverse of this expression is needed. Weconstruct the inverse by considering R ijkl ( ) R − jnlp ( ) A np = A ik . (B 6)In this context it is important to note that the velocity gradient covariance tensor issingular due to solenoidality of the velocity field. Consistently taking into account A ii =0, however, still allows to introduce the above definition of an inverse. It is readily checkedthat for traceless matrices the result reads R − jnlp ( ) = 25 1 (cid:104) u (cid:105) f (cid:48)(cid:48) u (0) (cid:2) − δ jn δ lp − δ jp δ ln (cid:3) . (B 7)B.2. First conditional moment and conditional Laplacian
To explicitly obtain the conditional Laplacian (2.11) in the Gaussian approximation, wefirst have to evaluate the first conditional moment in the Gaussian approximation. To9this end we consider the definition of the conditional moment (cid:10) A ( x ) (cid:12)(cid:12) A (cid:11) g ( A ) = (cid:90) d A A g ( A , A ) . (B 8)Here g and g denote the single- and two-point Gaussian distributions, respectively. Forthe following calculation we especially need the characteristic function of the two-pointGaussian distribution (3.9), which is related to the two-point PDF by inverse Fouriertransform: g ( A , A ) = (2 π ) − (cid:90) d Λ d Λ φ A ( Λ , Λ ) exp [ − i (Λ ,kl A ,kl + Λ ,kl A ,kl )] . (B 9)Inserting the last expression into (B 8) and noticing that A ,ij exp [ − i (Λ ,kl A ,kl + Λ ,kl A ,kl )] = i ∂∂ Λ ,ij exp [ − i (Λ ,kl A ,kl + Λ ,kl A ,kl )](B 10)we can make use of partial integration to obtain (cid:10) A ij ( x ) (cid:12)(cid:12) A (cid:11) g ( A ) = − i(2 π ) − (cid:90) d A d Λ d Λ (cid:20) ∂∂ Λ ,ij φ A ( Λ , Λ ) (cid:21) exp [ − i (Λ ,kl A ,kl + Λ ,kl A ,kl )] . (B 11)The derivative of the two-point characteristic function is readily obtained and reads ∂∂ Λ ,ij φ A ( Λ , Λ ) = − [ R ikjl ( r )Λ ,kl + R ikjl ( )Λ ,kl ] φ A ( Λ , Λ ) . (B 12)For this result we have also made use of the fact that R ijkl for homogeneous isotropicflows remains unchanged under a simultaneous change of indices i ↔ j and k ↔ l . Inanalogy to (B 10) we have the relationΛ ,ij exp [ − i (Λ ,kl A ,kl + Λ ,kl A ,kl )] = i ∂∂ A ,ij exp [ − i (Λ ,kl A ,kl + Λ ,kl A ,kl )](B 13)which together with (B 12) and (B 11) can be used to obtain the relation (cid:10) A ij ( x ) (cid:12)(cid:12) A (cid:11) g ( A ) = − R ikjl ( r ) ∂∂ A ,kl g ( A ) , (B 14)where we have carried out the integration with respect to A and identified the one-point PDF as the inverse Fourier transform of the one-point characteristic function. Thederivative of the Gaussian single-point PDF is readily obtained and yields ∂∂ A ,kl g ( A ) = − R − kmln ( ) A ,mn g ( A ) , (B 15)such that we obtain the final result (cid:10) A ij ( x ) (cid:12)(cid:12) A (cid:11) = R ikjl ( r ) R − kmln ( ) A ,mn . (B 16)By (2.11) we now have to evaluate the Laplacian of this expression (times the kinematicviscosity) which comes down to calculating lim r → ∆ r R ikjl ( r ). The calculation involvesagain a careful application of L’Hospital’s rule and yields the resultlim r → ∆ r R ikjl ( r ) = 718 (cid:104) u (cid:105) f (4) u (0) (cid:2) − δ ik δ jl + δ ij δ kl + δ il δ jk (cid:3) . (B 17)0 M. Wilczek and C. Meneveau
Together with the inverse (B 7) we obtain for the conditional Laplacian (cid:10) ν ∆ x A ij ( x ) (cid:12)(cid:12) A (cid:11) = ν (cid:104) lim r → ∆ r R ikjl ( r ) (cid:105) R − kmln ( ) A ,mn (B 18 a )= ν f (4) u (0) f (cid:48)(cid:48) u (0) A ,ij . (B 18 b )The prefactor of (B 18 b ) can also be conveniently expressed in terms of the energy spec-trum function. To establish the relation, we have to express the longitudinal velocityautocorrelation function in terms of the energy spectrum function, which can be achievedby considering R uij ( r ) = (cid:90) d k φ ij ( k ) exp[i k · r ] (B 19)where φ ij ( k ) = E ( k )4 πk (cid:18) δ ij − k i k j k (cid:19) (B 20)is the energy spectrum tensor and E ( k ) the energy spectrum function. This relation canbe used to evaluate lim r → ∂ ∂r l ∂r l R uii ( r ) = (cid:104) u (cid:105) f (cid:48)(cid:48) u (0) (B 21 a )= − (cid:90) d k k E ( k ) (B 21 b )leading to f (cid:48)(cid:48) u (0) = −
215 3 (cid:104) u (cid:105) (cid:90) d k k E ( k ) . (B 22)Similar calculations can be performed to derive the relation f (4) u (0) = 235 3 (cid:104) u (cid:105) (cid:90) d k k E ( k ) . (B 23)As a result, (B 18 b ) takes the simple form (cid:10) ν ∆ x A ( x ) (cid:12)(cid:12) A (cid:11) = − ν (cid:82) d k k E ( k ) (cid:82) d k k E ( k ) A . (B 24)B.3. Second conditional moment and conditional pressure Hessian
To evaluate the non-local contributions to the pressure Hessian in the Gaussian approx-imation following (2.12), we first need to evaluate the conditional second invariant (cid:10) Q ( x ) (cid:12)(cid:12) A (cid:11) = − (cid:10) A mn ( x ) A nm ( x ) (cid:12)(cid:12) A (cid:11) . (B 25)The evaluation of the second conditional moment (cid:10) A mn ( x ) A nm ( x ) (cid:12)(cid:12) A (cid:11) g ( A ) = (cid:90) d A A ,mn A ,nm g ( A , A ) (B 26)is analogous to the calculation of the first conditional moment and resembles many ofthe steps of the prior section. The result for the conditional second invariant eventually1reads (cid:10) Q ( x ) (cid:12)(cid:12) A (cid:11) = −
45 1 (cid:104) u (cid:105) f (cid:48)(cid:48) u (0) R mink ( r ) R nimk ( r ) −
15 1 (cid:104) u (cid:105) f (cid:48)(cid:48) u (0) R mink ( r ) R nkmi ( r ) − (cid:18) (cid:104) u (cid:105) f (cid:48)(cid:48) u (0) (cid:19) R mink ( r ) R njml ( r ) [4 A ,ik + A ,ki ] [4 A ,jl + A ,lj ] . (B 27)The result is a combination of tensor contractions of the velocity gradient covariancetensor and the velocity gradient tensor. It depends on the longitudinal velocity auto-correlation function through the velocity gradient covariance tensor, which can be seenwhen evaluating this expression for the general structure (B 1) of the velocity gradientcovariance tensor in homogeneous isotropic turbulence: (cid:10) Q ( x ) (cid:12)(cid:12) A (cid:11) = − f (cid:48)(cid:48) u (0) (cid:18) b (cid:104) u (cid:105) f (cid:48)(cid:48) u (0) + b Tr (cid:0) A (cid:1) + b Tr (cid:16) A A T1 (cid:17) + b (cid:16) ˆ r T A ˆ r (cid:17) + b (cid:16) ˆ r T A A T1 ˆ r (cid:17) + b (cid:16) ˆ r T A T1 A ˆ r (cid:17) + b (cid:16) ˆ r T A ˆ r (cid:17) (cid:19) (B 28)with b = − (cid:18) f (cid:48) u r (cid:19) + 2 f (cid:48) u r f (cid:48)(cid:48) u + 2 f (cid:48) u f (cid:48)(cid:48)(cid:48) u − f (cid:48)(cid:48) u − rf (cid:48)(cid:48) u f (cid:48)(cid:48)(cid:48) u (B 29 a ) b = 122 (cid:18) f (cid:48) u r (cid:19) + 86 f (cid:48) u r f (cid:48)(cid:48) u + 17 f (cid:48)(cid:48) u (B 29 b ) b = − (cid:18) f (cid:48) u r (cid:19) + 14 f (cid:48) u r f (cid:48)(cid:48) u + 8 f (cid:48)(cid:48) u (B 29 c ) b = − (cid:18) f (cid:48) u r (cid:19) + 308 f (cid:48) u r f (cid:48)(cid:48) u + 136 f (cid:48) u f (cid:48)(cid:48)(cid:48) u − f (cid:48)(cid:48) u − rf (cid:48)(cid:48) u f (cid:48)(cid:48)(cid:48) u (B 29 d ) b = 22 (cid:18) f (cid:48) u r (cid:19) − f (cid:48) u r f (cid:48)(cid:48) u + 32 f (cid:48) u f (cid:48)(cid:48)(cid:48) u − f (cid:48)(cid:48) u − rf (cid:48)(cid:48) u f (cid:48)(cid:48)(cid:48) u (B 29 e ) b = 22 (cid:18) f (cid:48) u r (cid:19) + 106 f (cid:48) u r f (cid:48)(cid:48) u + 32 f (cid:48) u f (cid:48)(cid:48)(cid:48) u − f (cid:48)(cid:48) u − rf (cid:48)(cid:48) u f (cid:48)(cid:48)(cid:48) u (B 29 f ) b = 50 (cid:18) f (cid:48) u r (cid:19) − f (cid:48) u r f (cid:48)(cid:48) u − f (cid:48) u f (cid:48)(cid:48)(cid:48) u + 500 f (cid:48)(cid:48) u + 50 rf (cid:48)(cid:48) u f (cid:48)(cid:48)(cid:48) u . (B 29 g )As expected these terms turn out to be invariants composed of A and ˆ r . It is in-teresting to check the limiting behaviour of this expression. In the limit of r → ∞ , theconditional second invariant vanishes due to the decay of correlations. In the limit r → b as well as b through b vanish, as a careful evaluation with L’Hospital’s rule shows. For b , however, we obtain b = 225 f (cid:48)(cid:48) u (0) in this limit, which yields the expected limitingbehaviour lim r → (cid:10) Q ( x ) (cid:12)(cid:12) A (cid:11) = Q . (B 30)Inserting expression (B 28) into (2.12), one realizes that the term independent of A doesnot give any contribution. A lengthy calculation reveals that the terms quadratic in A M. Wilczek and C. Meneveau can be grouped taking the form (cid:10) (cid:101) H ( x ) (cid:12)(cid:12) A (cid:11) = α (cid:18) S −
13 Tr (cid:0) S (cid:1) I (cid:19) + β (cid:18) W −
13 Tr (cid:0) W (cid:1) I (cid:19) + γ ( S W − W S ) (B 31)with α = − f (cid:48)(cid:48) u (0) (cid:90) d r (cid:18) f (cid:48) u r − f (cid:48) u f (cid:48)(cid:48) u r − f (cid:48) u f (cid:48)(cid:48)(cid:48) u r − f (cid:48)(cid:48) u r + f (cid:48)(cid:48) u f (cid:48)(cid:48)(cid:48) u (cid:19) (B 32 a ) β = − f (cid:48)(cid:48) u (0) (cid:90) d r (cid:18) f (cid:48) u r − f (cid:48) u f (cid:48)(cid:48) u r − f (cid:48) u f (cid:48)(cid:48)(cid:48) u r − f (cid:48)(cid:48) u r − f (cid:48)(cid:48) u f (cid:48)(cid:48)(cid:48) u (cid:19) (B 32 b ) γ = 475 f (cid:48)(cid:48) u (0) (cid:90) d r (cid:18) f (cid:48) u f (cid:48)(cid:48) u r − f (cid:48)(cid:48) u r − f (cid:48)(cid:48) u f (cid:48)(cid:48)(cid:48) u (cid:19) . (B 32 c )These terms can be significantly simplified by partial integration and identifying productrules, which then leads to α = −
27 (B 33 a ) β = −
25 (B 33 b ) γ = 625 + 1675 f (cid:48)(cid:48) u (0) (cid:90) d r f (cid:48) u f (cid:48)(cid:48)(cid:48) u r . (B 33 c )That means, the coefficients α and β become independent of the longitudinal velocityautocorrelation function, whereas this dependence remains for the term γ through anonlinear integral dependence. Appendix C. Parameter estimation from a model spectrum
The parameters γ and δ from the Gaussian approximation depend implicitly on theReynolds number through the longitudinal velocity autocorrelation function or, equiv-alently, through the energy spectrum function. To obtain approximate values for theparameters we use a simplified model spectrum similar to those proposed by Meyers& Meneveau (2008) or Pope (2000) (but unlike that proposed by Meyers & Meneveau(2008), without intermittency and bottleneck corrections). We use the basic form E ( k ) = C K ε / k − / F L ( kL ) F η ( kη ) with (C 1 a ) F L ( kL ) = (cid:34) c kL (cid:2) ( c kL ) / + 1 (cid:3) / (cid:35) / and (C 1 b ) F η ( kη ) = exp ( − c kη ) . (C 1 c )Here, C K is the Kolmogorov constant, ε is the dissipation rate, L denotes the integrallength scale, η denotes the Kolmogorov length scale and c and c are non-dimensionalparameters which are determined such that the length scales in the above expression areconsistent with the relations L = (2 E kin / / /ε and η = ( ν /ε ) / . The kinetic energy3and the rate of energy dissipation are given by the integrals E kin = (cid:90) d k E ( k ) (C 2) ε = 2 ν (cid:90) d k k E ( k ) . (C 3)For a given dissipation rate, the kinematic viscosity then is determined by (C 3), and theTaylor Reynolds number can be estimated as R λ = (cid:112) / E kin / √ εν . For a given energyspectrum function, the longitudinal velocity autocorrelation is given by f u ( r ) = 6 (cid:104) u (cid:105) (cid:90) d k E ( k ) (cid:20) sin( kr )( kr ) − cos( kr )( kr ) (cid:21) . (C 4)Also derivatives of the velocity autocorrelation function can be obtained straightfor-wardly, such that we can evaluate (3.13 c ) to obtain a numerical estimate for γ ; δ isreadily obtained from a numerical evaluation of (3.11) in terms of the energy spectrumfunction.For the present example we choose C K = 1 . ε = 1 . / s as well as c L = 1 . c η = 0 .
002 m, which results in a Reynolds number of R λ ≈ γ ≈ .
08 and δ τ η ≈ − . Appendix D. Implementation of noise term in the SDE
For the discussion of the implementation of the noise term we mainly follow Chevillard et al. (2008), to which we refer the reader for more background information. Note that ourpresentation differs with respect to some coefficients (which are a matter of convention).The noise term in (4.6) can be written as (Chevillard et al. F (cid:63)ij = A F D ijkl d W (cid:63)kl , (D 1)where A F is the forcing amplitude and d W (cid:63)kl in this section denotes an isotropic tensorialWiener process specified by (cid:104) d W (cid:63) (cid:105) = 0 and (cid:104) d W (cid:63)ik d W (cid:63)jl (cid:105) = δ ij δ kl d t (cid:63) . Here D is a tensorspecified such that the noise term complies with the tensorial structure of the forcingapplied to the velocity gradient tensor evolution equation discussed in appendix A.2. Werecall that this forcing is specified by (cid:104) d F (cid:63) (cid:105) = 0 and (cid:104) d F (cid:63)ik d F (cid:63)jl (cid:105) = Q (cid:63)ijkl ( )d t (cid:63) , wherethe Q (cid:63) = Q τ η is the non-dimensionalized version of the force covariance tensor (A 12).Now identifying the forcing amplitude as A F = (cid:0) − σ F f (cid:48)(cid:48) F (0) τ η (cid:1) / leads to the condition D ikmn D jlmn = 2 δ ij δ kl − δ ik δ jl − δ il δ jk . (D 2)The solution to this equation has been worked out by Chevillard et al. (2008) and reads D ijkl = 13 3 + √ √
10 + √ δ ij δ kl − √
10 + √ δ ik δ jl + 1 √
10 + √ δ il δ jk , (D 3)which concludes the specification of the stochastic forcing term. REFERENCESAshurst, W. T., Kerstein, A. R., Kerr, R. M. & Gibson, C. H.
Phys. Fluids (8), 2343–2353. M. Wilczek and C. Meneveau
Betchov, R.
J. Fluid Mech. (05), 497–504. Cantwell, B. J.
Phys. Fluids A (4), 782–793. Chertkov, M., Pumir, A. & Shraiman, B. I.
Phys. Fluids (8), 2394–2410. Chevillard, L. & Meneveau, C.
Phys. Rev. Lett. , 174501. Chevillard, L., Meneveau, C., Biferale, L. & Toschi, F.
Phys. Fluids (10), 101504. Dopazo, C.
Turbulent Reacting Flows II , chap. Recent developments in PDF methods.Springer.
Friedrich, R., Daitche, A., Kamps, O., Luelff, J., Vosskuhle, M. & Wilczek, M.
C. R. Phys. (9-10), 929 – 953. Frisch, U.
Turbulence: the legacy of A.N. Kolmogorov . Cambridge University Press.
Girimaji, S. S. & Pope, S. B.
Phys. Fluids (2), 242–256. Haken, H.
Synergetics: Introduction and Advanced Topics . Springer.
Holzer, M. & Siggia, E.
Phys. Fluids A (10), 2525–2532. Ishihara, T., Kaneda, Y., Yokokawa, M., Itakura, K. & Uno, A.
J. Fluid Mech. , 335–366.
Jeong, E. & Girimaji, S. S.
Theoret. Comput. Fluid Dynamics , 421–432. Li, Y., Perlman, E., Wan, M., Yang, Y., Meneveau, C., Burns, R., Chen, S., Szalay,A. & Eyink, G.
J. Turbul. p. N31.
Lund, T. S. & Rogers, M. M.
Phys. Fluids (5), 1838–1847. Lundgren, T. S.
Phys.Fluids (5), 969–975. L¨uthi, B., Holzner, M. & Tsinober, A.
J. Fluid Mech. , 497–507.
Martin, J., Dopazo, C. & Valino, L.
Phys. Fluids (8), 2012–2025. Meneveau, C.
Annu. Rev. Fluid Mech. (1), 219–245. Meyers, J. & Meneveau, C.
Physics of Fluids (6), 065109. Monin, A. S., Yaglom, A. M. & Lumley, J. L.
Statistical fluid mechanics: mechanicsof turbulence , Statistical Fluid Mechanics , vol. 2. Dover Publications.
Naso, A. & Pumir, A.
Phys. Rev. E , 056318. Nomura, K. K. & Post, G. K.
J. Fluid Mech. , 65–97.
Ohkitani, K. & Kishiba, S.
Phys. Fluids (2), 411–421. Pope, S. B.
Turbulent Flows . Cambridge University Press.
Shtilman, L., Spector, M. & Tsinober, A.
J. Fluid Mech. , 65–77.
Tsinober, A.
Eur. J. Mech. B/Fluids , 421–449. Tsinober, A.
An Informal Conceptual Introduction to Turbulence . Springer. Vieillefosse, P.
J. Phys. France (6), 837–842. Vieillefosse, P.
PhysicaA (1), 150 – 162.
Wilczek, M., Daitche, A. & Friedrich, R.
J. Fluid Mech.676