A multifractal model for the velocity gradients dynamics in turbulent flows

TL;DR

This paper introduces a multifractal stochastic model for velocity gradient dynamics in turbulence that accurately captures small-scale statistics and intermittency effects across different Reynolds numbers.

Contribution

The authors develop a novel non-linear, non-Markovian model incorporating multifractal intermittency, improving the realism of velocity gradient statistics in turbulent flows.

Findings

Stationary regime with gradients variance proportional to Reynolds number

Gradients correlated over Kolmogorov time scale, norms over integral time scale

Reproduces characteristic teardrop shape of joint PDF of invariants Q and R

Abstract

We develop a stochastic model for the velocity gradients dynamics along a Lagrangian trajectory. Comparing with different attempts proposed in the literature, the present model, at the cost of introducing a free parameter known in turbulence phenomenology as the intermittency coefficient, gives a realistic picture of velocity gradient statistics at any Reynolds number. To achieve this level of accuracy, we use as a first modelling step a regularized self-stretching term in the framework of the Recent Fluid Deformation (RFD) approximation that was shown to give a realistic picture of small scales statistics of turbulence only up to moderate Reynolds numbers. As a second step, we constrain the dynamics, in the spirit of Girimaji & Pope (1990), in order to impose a peculiar statistical structure to the dissipation seen by the Lagrangian particle. This probabilistic closure uses as a building block a random field that fulfils the statistical description of the intermittency, i.e. multifractal, phenomenon. To do so, we define and generalize to a statistically stationary framework a proposition made by Schmitt (2003). These considerations lead us to propose a non-linear and non-Markovian closed dynamics for the elements of the velocity gradient tensor. We numerically integrate this dynamics and observe that a stationary regime is indeed reached, in which (i) the gradients variance is proportional to the Reynolds number, (ii) gradients are typically correlated over the (small) Kolmogorov time scale and gradients norms over the (large) integral time scale (iii) the joint probability distribution function of the two non vanishing invariants $Q$ and $R$ reproduces the characteristic teardrop shape, (iv) vorticity gets preferentially aligned with the intermediate eigendirection of the deformation tensor and (v) gradients are strongly non-Gaussian and intermittent.

Equations

$$A_{ij}=\frac{\partial u_i}{\partial x_j}.$$

$$\frac{d{A}_{ij}}{dt}=-A_{ik}{A}_{kj}-\frac{{\partial }^{2}p}{\partial x_i\partial x_j}+\nu \Delta A_{ij},$$

$$d{A}_{ij}=V_{ij}dt+D_{ijkl}{W}_{kl}(dt),$$

$$\text{mathsfbi}{V}^{\text{mbox}RE}=-\text{mathsfbi}{A}^2+\frac{\text{mbox}tr(\text{mathsfbi}{A}^2)}{3}\text{mathsfbi}I,$$

$$\text{mathsfbi}{V}^{\text{mbox}RFD}=-\text{mathsfbi}{A}^2+\frac{\text{mbox}tr(\text{mathsfbi}{A}^2)}{\text{mbox}tr(\text{mathsfbi}{C}_{{\tau }_{\eta }}^{-1})}\text{mathsfbi}{C}_{{\tau }_{\eta }}^{-1}-\frac{\text{mbox}tr(\text{mathsfbi}{C}_{{\tau }_{\eta }}^{-1})}{3T}\text{mathsfbi}A,$$

$$\text{mathsfbi}{C}_{{\tau }_{\eta }}=e^{{\tau }_{\eta }\text{mathsfbi}A}{e}^{{\tau }_{\eta }\text{mathsfbi}{A}^{\top }}.$$

$$\mathbb{E}(\phi )=\mathbb{E}\left[\text{mbox}tr(\text{mathsfbi}A\text{mathsfbi}{A}^{\top })\right]=\frac{1}{{\tau }_{\eta }^2},$$

$$d{A}_{ij}=\left[V_{ij}^{\text{mbox}RFD}+f(t)A_{ij}\right]dt+\frac{1}{2}\sqrt{\frac{{\mu }^l\phi }{{\tau }_{\eta }}}\left[W_{ij}-\frac{1}{3}\text{mbox}tr(\text{mathsfbi}W){\delta }_{ij}\right].$$

$${\phi }_{\tau }(t)=\frac{1}{\tau }{\int }_{s\in [t-\tau ,t]}\phi (s)ds.$$

$$\phi (\mathbf{x})=\frac{1}{{\tau }_{\eta }^2}{e}^{\sqrt{{\mu }^e}{X}_{\eta }(\mathbf{x})-\frac{{\mu }^e}{2}\mathbb{E}\left[X_{\eta }^2\right]},$$

$$X_{\eta }(\mathbf{x})=\frac{1}{4\pi }{\int }_{|\mathbf{x}-\mathbf{y}|\le L}\frac{1}{|\mathbf{x}-\mathbf{y}{|}_{\eta }^{3/2}}W(dy),$$

$$\mathbb{E}\left[X_{\eta }^2\right]\underset{\eta \to 0}{\sim }\log \left(\frac{L}{\eta }\right).$$

$$\underset{\eta \to 0}{\lim }\mathbb{E}\left[X_{\eta }(\mathbf{x})X_{\eta }(\mathbf{x}+\mathbf{h})\right]={\log }^{+}\left(\frac{L}{|\mathbf{h}|}\right)+g(|\mathbf{h}|),$$

$$\mathbb{E}\left[{\phi }^q\right]\underset{\eta \to 0}{\sim }\frac{e^{\frac{{\mu }^e}{2}q(q-1)g(0)}}{{\tau }_{\eta }^{2q}}{\left(\frac{L}{\eta }\right)}^{\frac{{\mu }^e}{2}q(q-1)}.$$

$${\phi }_{\ell }(\mathbf{x})=\frac{1}{\frac{4}{3}\pi {\ell }^3}{\int }_{|\mathbf{x}-\mathbf{y}|\le \ell }\phi (\mathbf{y})dy,$$

$$\underset{\eta \to 0}{\lim }\mathbb{E}\left[{\phi }_{\ell }^q\right]\underset{\ell \to 0}{\sim }{C}_q{\left(\frac{\ell }{L}\right)}^{\frac{{\mu }^e}{2}q(1-q)},$$

$$\phi (t)=\frac{1}{{\tau }_{\eta }^2}{e}^{\sqrt{{\mu }^l}{X}_{{\tau }_{\eta }}(t)-\frac{{\mu }^l}{2}\mathbb{E}\left[X_{{\tau }_{\eta }}^2\right]},$$

$$d{X}_{{\tau }_{\eta }}(t)=\left[-\frac{1}{T}{X}_{{\tau }_{\eta }}(t)+{\beta }_{{\tau }_{\eta }}(t)\right]dt+\frac{1}{\sqrt{{\tau }_{\eta }}}W(dt),$$

$${\beta }_{{\tau }_{\eta }}(t)=-\frac{1}{2}{\int }_{s=-\infty }^t\frac{1}{(t-s+{\tau }_{\eta }{)}^{3/2}}W(ds).$$

$$X_{{\tau }_{\eta }}(t)={\int }_{s=-\infty }^t{e}^{-\frac{t-s}{T}}{\beta }_{{\tau }_{\eta }}(s)ds+\frac{1}{\sqrt{{\tau }_{\eta }}}{\int }_{s=-\infty }^t{e}^{-\frac{t-s}{T}}W(ds),$$

$$\mathbb{E}\left[{\phi }^q\right]\underset{{\tau }_{\eta }\to 0}{\propto }{\left(\frac{T}{{\tau }_{\eta }}\right)}^{\frac{{\mu }^l}{2}q(q-1)},$$

$$\underset{{\tau }_{\eta }\to 0}{\lim }{\tau }_{\eta }^{2q}\mathbb{E}\left[{\phi }_{\tau }^q\right]\underset{\tau \to 0}{\propto }{\left(\frac{T}{\tau }\right)}^{\frac{{\mu }^l}{2}q(q-1)}{e}^{\frac{{\mu }^l}{2}q(q-1)\tilde{g}(0)}{\int }_{[0,1{]}^q}\prod _{i<j}^q\frac{1}{|s_i-s_j{|}^{{\mu }^l}}\prod _{i=1}^{q}d{s}_i,$$

$$d{A}_{ij}(t)=V_{ij}(t)dt+D_{ijkl}(t)W_{kl}(dt),$$

$$d\phi =d\text{mbox}tr(\text{mathsfbi}A\text{mathsfbi}{A}^{\top })=\left[2A_{ij}{V}_{ij}+D_{ijkl}{D}_{ijkl}\right]dt+2A_{ij}{D}_{ijkl}{W}_{kl}(dt).$$

$$D_{ijkl}=a{\delta }_{ij}{\delta }_{kl}+b{\delta }_{ik}{\delta }_{jl}+c{\delta }_{il}{\delta }_{kj},$$

$$3a+b+c=0.$$

$$d\phi =d\text{mbox}tr(\text{mathsfbi}A\text{mathsfbi}{A}^{\top })=\left[2A_{ij}{V}_{ij}+9(a^2+b^2+c^2)+6(ab+ac+bc)\right]dt+2(b{A}_{kl}+c{A}_{lk})W_{kl}(dt).$$

$$d{\phi }_{{\tau }_{\eta }}$$

$$+\sqrt{\frac{{\mu }^l}{{\tau }_{\eta }}}{\phi }_{{\tau }_{\eta }}(t)W(dt).$$

$$4\left[(b^2+c^2)\phi +2bc\text{mbox}tr(\text{mathsfbi}{A}^2)\right]=\frac{{\mu }^l}{{\tau }_{\eta }}{\phi }^2.$$

$$a=-\frac{1}{6}\sqrt{\frac{{\mu }^l\phi }{{\tau }_{\eta }}}\text{mbox},b=\frac{1}{2}\sqrt{\frac{{\mu }^l\phi }{{\tau }_{\eta }}}\text{mbox}andc=0.$$

Full text

††volume: ..

1 Introduction

The study of the statistical properties of the velocity gradients dynamics has shed a new light on hydrodynamics turbulence. In the past years (Tsinober, 2001; Wallace, 2009; Meneveau, 2011), quantitative progress has been made while focusing on the small scales of turbulence, i.e. below the Kolmogorov length scale, in particular on the statistical and geometrical properties of the velocity gradient tensor $\mathsfbi{A}$, defined as the $3\times 3$ matrix made up of the gradients of the three components of velocity $u_{i}$, for $1\leq i\leq 3$ along the three spatial directions $x_{j}$, with $1\leq j\leq 3$, namely

[TABLE]

Henceforth, as it is classically used in fluid mechanics, we will denote by $\mathsfbi{S}$ and $\mathsfbi{\Omega}$ the symmetric and antisymmetric decompostion of $\mathsfbi{A}$, respectively, such that $\mathsfbi{A}=\mathsfbi{S}+\mathsfbi{\Omega}$. It has been recognized that many of the geometrical properties of $\mathsfbi{A}$, such as the relative amplitudes and signs of the eigenvalues of $\mathsfbi{S}$ and the alignments of the vorticity vector $\boldsymbol{\omega}=\boldsymbol{\nabla}\wedge\boldsymbol{u}$ (related to $\mathsfbi{\Omega}$ according to $\mathsfbi{\Omega}\boldsymbol{h}=\frac{1}{2}\boldsymbol{\omega}\wedge\boldsymbol{h}$ for all vectors $\boldsymbol{h}\in\mathbb{R}^{3}$) with respect to the eigenframe of $\mathsfbi{S}$, can be understood in a kinematic sense while considering the Lagrangian evolution of $\mathsfbi{A}$. This dynamics is obtained by taking a spatial derivative of the Navier-Stokes equations and reads

[TABLE]

where $d/dt\equiv\partial/\partial t+\boldsymbol{u}\cdot\boldsymbol{\nabla}$ stands for the material derivative, $p$ is the pressure field determined by the incompressibility condition, hence solution of the respective Poisson equation $\Delta p=-\mbox{tr}(\mathsfbi{A}^{2})$, and $\nu$ the kinematic viscosity (Tsinober, 2001; Wallace, 2009; Meneveau, 2011).

The transport equation (2) of the velocity gradient tensor $\mathsfbi{A}$ states that its time variation along a Lagrangian trajectory is governed by the action of three terms: the self-stretching term $-\mathsfbi{A}^{2}$, the pressure Hessian and the viscous diffusion. Let us mention that (2) is not closed in terms of the temporal profile of the nine elements of $\mathsfbi{A}$ along the trajectory of the fluid particle under consideration since both the pressure and viscous terms require the knowledge of the whole spatial field. Indeed, the viscous term comprises the gradients of velocity gradients, which are not known along a single trajectory, and it is known (Ohkitani, 1993; Chevillard et al., 2008; Wilczek & Meneveau, 2014) that the pressure Hessian can be expressed as a convolution over space of $\mbox{tr}(\mathsfbi{A}^{2})$, as a consequence of the Poisson equation. It is nonetheless tempting to focus on a trajectory of a single particle, disregarding all the others, thus decreasing drastically the number of degrees of freedom. This implies studying an approximative dynamics based on simplistic closures of the pressure Hessian and viscous terms that are, hopefully, realistic enough to reproduce the observed statistics of $\mathsfbi{A}$ in experimental and numerical flows.

In this spirit, by neglecting the anisotropic part of the pressure Hessian entering (2) and taking $\nu=0$ one builds the simplest closure that preserves incompressibility. This is the so-called restricted Euler (RE) approximation (Vieillefosse, 1982; Cantwell, 1992, see also Meneveau 2011 for a review). It is known that this closure leads to a finite time singularity of the elements of $\mathsfbi{A}$ for any non-vanishing initial condition. This predicted singularity is unphysical, although this closure allows to understand, among others, one experimental observation: vorticity has the tendency to align with the eigendirection of the deformation $\mathsfbi{S}$ associated to the intermediate eigenvalue. The RE approximation is thus an appealing starting point to design more realistic closures of pressure and diffusivity able, at least, to regularize the finite time singularity implied by the self-stretching term, while keeping track of the underlying non-linear dynamics. Closures can then be compared against experimental measurements and numerical simulations. Usually, this closed dynamics is associated with a tensorial random forcing, delta-correlated in time, to maintain a statistically steady state (Meneveau, 2011). In the stationary regime, the statistics of $\mathsfbi{A}$ can then be compared to those obtained, for instance, in numerical simulations of homogeneous and isotropic turbulence, available at several databases in the world (see for instance Li et al., 2008). Recently, several such closures for the pressure Hessian and the viscous term have been proposed in the literature (Chevillard & Meneveau, 2006; Wilczek & Meneveau, 2014; Johnson & Meneveau, 2016). They all show realistic statistics of, among other properties, the joint probability distribution function (PDF) of the two invariants $Q$ and $R$, and the peculiar preferential alignment of vorticity with the eigenframe of deformation (we will define precisely these quantities later in the article). Moreover, with different levels of success, they even compared well against more precise estimations of the pressure field, such as conditional averages of the pressure Hessian given an instance of $\mathsfbi{A}$, and more precisely, of a joint instance of $R$ and $Q$ (Chevillard et al., 2008; Meneveau, 2011; Chevillard et al., 2011; Wilczek & Meneveau, 2014; Johnson & Meneveau, 2016). Whereas all these closed dynamics have in common to preserve the self stretching term $-\mathsfbi{A}^{2}$, they unfortunately, to our knowledge, give unrealistic statistics when the Reynolds number $\mathcal{R}_{e}$, introduced as a natural free parameter of the closures, becomes too large, preventing the model from reaching the asymptotic limit of infinite Reynolds numbers. They nonetheless provide realistic closures at moderate Reynolds numbers, and they are at least able to regularize the finite-time singularity implied by the self-stretching term $-\mathsfbi{A}^{2}$. To this regards, they appear as good starting points to build up more sophisticated models that can describe flows at arbitrarily large Reynolds Numbers. Let us now define more precisely the underlying dynamics that we will study. Consider the following stochastic differential equation for the velocity gradient tensor $\mathsfbi{A}$

[TABLE]

where $V_{ij}$ and $D_{ijkl}$ are called respectively the drift and diffusion terms, and $W_{kl}(dt)$ is a realization of a white noise stochastic process, defined such that all components are independent and, loosely speaking, delta-correlated in time. Let us assume at this stage that, following this dynamics (3), $\mathsfbi{A}$ reaches a statistically stationary regime, characterized for example by the covariance of its elements, supposed thus finite. We are then asking the tensor terms $V_{ij}$ and $D_{ijkl}$ from (3) to be, as a necessary condition, also of finite variance. In the framework of additive noise, the diffusion term is deterministic, and the covariance of its elements can be chosen such that it is consistent with a fourth-order isotropic tensor, ensuring furthermore that the trace-free condition imposed by incompressibility is respected (Chevillard et al., 2008; Johnson & Meneveau, 2016). As for the drift term $V_{ij}$, let us for example consider the one predicted by the restricted Euler (RE) approximation, that reads in matrix form

[TABLE]

where $\mathsfbi{I}$ stands for the identity matrix, and the one constructed under the Recent Fluid Deformation (RFD) approximation, namely

[TABLE]

where $\tau_{\eta}$ and $T$ are respectively the Kolmogorov and integral time scales, the two characteristic scales of turbulence at a given Reynolds number $\mathcal{R}_{e}$, and $\mathsfbi{C}_{\tau_{\eta}}$ the short-time Cauchy-Green tensor

[TABLE]

The RFD approximation and its consequences have been discussed in Chevillard & Meneveau (2006); Chevillard et al. (2008); Meneveau (2011), and we will not recall them here. We will keep in mind that whereas the RE approximation leads to finite time singularity and is independent of the Reynolds number, the RFD approximation regularizes this finite-time singularity through the joint action of the modelled pressure Hessian and viscous term using the short time Cauchy-Green tensor $\mathsfbi{C}_{\tau_{\eta}}$ (6). Furthermore, the Reynolds number becomes explicit through the ratio of the two time scales $T/\tau_{\eta}\propto\sqrt{\mathcal{R}_{e}}$. Similar drift terms can be derived from the approaches developed in Wilczek & Meneveau (2014) and Johnson & Meneveau (2016). As commented before, such closures give a realistic picture of the statistics of $\mathsfbi{A}$ at moderate Reynolds numbers, i.e. when the ratio $\tau_{\eta}/T$ does not become too small. As an example, the RFD approximation gives realistic results of homogeneous and isotropic turbulence (of Taylor-based Reynolds number $\mathcal{R}_{\lambda}=140$) when $\tau_{\eta}/T=0.1$, and deteriorate for $\tau_{\eta}/T<0.01$ (Chevillard et al., 2008). Subsequent numerical and theoretical developments in the framework of the RFD approximation indeed show that the model pointed in (5) fails at giving a realistic picture of turbulence at high Reynolds numbers (Afonso & Meneveau, 2010; Moriconi et al., 2014; Grigorio et al., 2017). We may wonder whether it is possible to constrain the drift $V_{ij}$ and diffusion $D_{ijkl}$ terms in a different way to ensure, for instance, a basic property of turbulence (Frisch, 1995) that is the finiteness of the average dissipation rate $2\nu\mathbb{E}\left[\mbox{tr}(\mathsfbi{S}^{2})\right]$ when $\mathcal{R}_{e}\to\infty$.

Actually, such an approach has already been explored in this context by Girimaji & Pope (1990). They end up with a very different dynamics for $\mathsfbi{A}$ (in the spirit of (3)), where a diffusion term $D_{ijkl}$ that depends explicitly on $\mathsfbi{A}$ appears, thus lying in the class of multiplicative noise models. More precisely, $V_{ij}$ and $D_{ijkl}$ are constrained to ensure that a particular contraction of $\mathsfbi{A}$, namely the pseudo-dissipation $\varphi=\mbox{tr}(\mathsfbi{A}\mathsfbi{A}^{\top})$ (in units of viscosity $\nu$), follows the dynamics of an exponentiated Ornstein-Uhlenbeck (OU) process. They thus impose that the dynamics of $\mathsfbi{A}$ (3), in particular, fulfils the dimensional Kolmogorov’s prediction

[TABLE]

which thus defines precisely what we mean by the Kolmogorov time scale $\tau_{\eta}$. Consequently, imposing the condition (7), they are able to reach any Reynolds numbers, as far as the variance of the elements of $\mathsfbi{A}$ are concerned. But assuming that $\varphi$ is an exponentiated OU process has further consequences. In particular, $\log(\varphi)$ is thus a Gaussian process, whose average is set by fulfilling the Kolmogorov prediction given in (7), and whose covariance decays exponentially along Lagrangian trajectories over a characteristic time scale. This was already checked in early direct numerical simulations of Lagrangian turbulence (Yeung & Pope, 1989), where it is observed that indeed, to a good approximation, $\log(\varphi)$ is a Gaussian random process whose covariance is consistent with an exponential decrease over the large time scale $T$ (or integral turnover time scale) of turbulence. This stochastic modelling of $\varphi$ compares thus well at this stage to numerical flows and is defined up to a free parameter, called $\hat{a}$ in Girimaji & Pope (1990), that has to be determined from empirical data. Unfortunately, higher order moments of $\varphi$ are poorly predicted, and do not compare appropriately with the $\mathcal{R}_{e}$ dependence of the fluctuations observed in experimental (Van Atta & Antonia, 1980; Sreenivasan & Antonia, 1997) and numerical (Ishihara et al., 2007) flows. We will develop these ideas in section 2.1. Let us just keep in mind that this modelling needs to be improved in order to be realistic in face of observed fluctuations. To this regard, it turns out that the $\mathcal{R}_{e}$-dependence of moments of velocity derivatives observed experimentally and numerically (Sreenivasan & Antonia, 1997; Ishihara et al., 2007) are well described by the multifractal formalism, as depicted, for instance, in Borgas (1993) and Frisch (1995). We need thus to develop a new stochastic model for pseudo-dissipation able to reproduce the observed fluctuations of the velocity gradients, and furthermore to be consistent with the multifractal picture.

Such a random field, able to reproduce the statistics of the pseudo-dissipation, is at the core of the work of Kolmogorov (1962) and Obukhov (1962), where the dissipation field, and alternatively, the pseudo-dissipation field, determine the statistical features of the velocity increments through the refined similarity hypothesis (Frisch, 1995). In an Eulerian description of the flow, the spatial field $\varphi$ is modelled in a homogeneous way by taking the exponential of a Gaussian field, whose covariance decreases logarithmically over the integral length scale $L$. This stochastic construction was first proposed by Mandelbrot (1972) following a construction due to Yaglom (1966), and formalised by Kahane (1985) in the context of Multiplicative Chaos Theory (Rhodes & Vargas, 2014). Going back to the velocity gradient dynamics where the Lagrangian point of view has to be adopted, we need a temporal form of such a multiplicative chaos. This was considered by Schmitt (2003), where, as a first step, it is shown how to adapt the approach of Kahane (1985) in a Lagrangian context that requires, in particular, a causal construction. Schmitt (2003) advances also some heuristics in order to build a dynamics, i.e. a stochastic differential equation, for which the statistically stationary solution is indeed such a causal multifractal process. We will develop these ideas in section 2.2.

The purpose of the present article is to include the multifractal picture given by the formerly described multiplicative chaos approach into the dynamics of the velocity gradient tensor (3). In the sequel, we will develop this idea and finally propose the following stochastic model for $\mathsfbi{A}$ along a Lagrangian trajectory

[TABLE]

Let us now comment on the various terms entering the proposed dynamics (8). First, the drift term proportional to $dt$ is made up of the dynamics given by the RFD closure (6), supplemented by an additional damping term proportional to $\mathsfbi{A}$ itself. The scalar function $f(t)$, given in (34), depends explicitly on $\mathsfbi{A}$ at present time $t$ and on an additional random variable (given in (32)) built on the past values of both $\mathsfbi{A}$ and the white noise tensor field $\mathsfbi{W}$. Its presence is crucial to ensure multifractal properties in the stationary regime, and makes the overall dynamics (8) non-Markovian. Secondly, the diffusion term is very similar to the one used by Girimaji & Pope (1990) that, as mentioned, have a multiplicative nature. It is given by the product of the trace-free part of the white tensor noise $\mathsfbi{W}(dt)$ and the pseudo-dissipation $\varphi=\mbox{tr}(\mathsfbi{A}\mathsfbi{A}^{\top})$ evaluated at present time $t$, and its intensity is moderated by a free parameter $\mu^{l}$, interpreted as the intermittency coefficient from turbulence phenomenology, for reasons that will become clearer.

In section 2, we review and recall the two existing stochastic models for the pseudo-dissipation seen by a fluid particle. This includes the proposition made by Pope (1990), that is to take the exponential of an Ornstein-Uhlenbeck process, and the one of Schmitt (2003) that we define rigorously and extend in order to deal with statistically stationary processes. In section 3, following the idea of Girimaji & Pope (1990), we show how to include the multifractal picture into the dynamics of the velocity gradient tensor. We end up with a proposition for a closed tensorial stochastic differential equation (8), which is non-linear and non-Markovian. Section 4 is devoted to the description of the numerical procedure used to simulate the trajectories of our process (8). We present then in section 5 the results of numerical simulations of the trajectories for various Reynolds numbers and discuss their comparison with known empirical facts of turbulence. We gather in section 6 our conclusions and perspectives.

2 Stochastic models of pseudo-dissipation

2.1 Pseudo-dissipation as an exponentiated Ornstein-Uhlenbeck process

As a building block of their stochastic model for the velocity gradients, Girimaji & Pope (1990) have considered for pseudo-dissipation $\varphi$ the exponential of an Ornstein-Uhlenbeck process. This model is justified by the observations made in early numerical simulations of Yeung & Pope (1989) that the shape of the one point PDF of $\log\varphi$ is close to a Gaussian and, furthermore, the autocorrelation of $\log\varphi(t)$ follows an exponential decrease over a typical time scale given by the integral time scale $T$ of turbulence (Pope, 1990). At this point, we may remark that the small scale quantity $\varphi$ is correlated over the large energy containing scale of turbulence, which is characteristic of the absence of scale decoupling of turbulence. This was already recognized in the Eulerian framework where it is observed that dissipation is correlated in space over the large integral length scale $L$ (Gagne & Hopfinger, 1979; Antonia et al., 1981). From the theoretical side, this long-range correlated nature of the small scales is at the heart of the models of intermittency of Kolmogorov (1962) and Obukhov (1962) (see also Novikov 1989, 1990 for discussions in the lagrangian framework), although, at that stage, nothing is said about how the correlation function decreases over $T$, and as we will see, multifractal phenomenology requires a logarithmic decrease (see subsection 2.2) which is not reproduced by the proposition of Pope (1990).

The statistical properties of this model are reviewed in appendix A. Its causal dynamics is given in (51) and its unique solution, given an initial condition, is given explicitly in (52). Let us call $\hat{a}$, as in Girimaji & Pope (1990), the free parameter of this model. It can be conveniently related to the variance of the logarithm of pseudo-dissipation as $\mathbb{E}\log^{2}\varphi-\mathbb{E}^{2}\log\varphi=2\hat{a}^{2}T$, or equivalently, to the moments $\mathbb{E}(\varphi^{q})=\frac{1}{\tau_{\eta}^{2q}}e^{\hat{a}^{2}Tq(q-1)}$ (53). Whatever the value of $\hat{a}$, this model fulfils the basics of Kolmogorov’s phenomenology, namely $\mathbb{E}(\varphi)=1/\tau_{\eta}^{2}$ (7). As it is explained and observed in Yeung & Pope (1989), $\hat{a}^{2}$ is expected to depend logarithmically on the Reynolds number (see (55)) in order to be consistent with both empirical observations and the phenomenology of Kolmogorov (1962) and Obukhov (1962) such as $\hat{a}^{2}\sim\log\mathcal{R}_{e}$ when $\mathcal{R}_{e}\to\infty$.

To build a more precise connection with Kolmogorov (1962) and Obukhov (1962) phenomenology, and more generally with the multifractal (i.e. intermittent) phenomenology (Borgas, 1993; Frisch, 1995), let us consider the average value of pseudo-dissipation over a time interval of extension $\tau$, i.e.

[TABLE]

This quantity is the Lagrangian analogue of the (pseudo-)dissipation averaged over a ball of size $\ell$ as considered by Kolmogorov (1962) and Obukhov (1962) in the refined similarity hypothesis relating velocity structure functions and dissipation field (Borgas, 1993; Frisch, 1995). It is shown in appendix A (see equation 54) how to compute the high-order moments of $\varphi_{\tau}$, and it is easy to be convinced that if $\hat{a}^{2}T$ is taken proportional to $\log\mathcal{R}_{e}$, as it is required by the intermittent phenomenology formerly described, then, for any $\tau>0$, we get diverging high-order moments of the coarse-grained pseudo-dissipation, i.e., for $q>1$, $\tau_{\eta}^{2q}\mathbb{E}[\varphi_{\tau}^{q}]\to\infty$ when $\mathcal{R}_{e}\to\infty$, which is not consistent with the refined similarity hypothesis (see the discussion on the Reynolds number dependence of the model in appendix A). Thus, taking as a stochastic model for pseudo-dissipation the exponential of an Ornstein-Uhlenbeck process prevents us from fulfilling the statistical properties required by the phenomenology of Kolmogorov (1962) and Obukhov (1962). We will see in the following section that, on the other hand, a pseudo-dissipation model given by the multiplicative chaos allows one to obtain (i) moments $\mathbb{E}[\varphi^{q}]$ that scale as a power of the Reynolds number, and (ii) bounded moments $\varphi_{\tau}$ for any finite $\tau>0$ when $\mathcal{R}_{e}\to\infty$.

2.2 Pseudo-dissipation as a causal and stationary multiplicative chaos

Multiplicative chaos has been introduced by Mandelbrot (1972) in the context of turbulence, following the work of Yaglom (1966), and formalized by Kahane (1985) in order to give a stochastic meaning to the dissipation field as depicted in the phenomenology of Kolmogorov (1962) and Obukhov (1962) (see recent mathematical developments in Rhodes & Vargas (2014)). In the Eulerian framework, for which we need to define spatial fields, dissipation is taken as the exponential of a isotropic Gaussian field logarithmically correlated over the integral length scale $L$. As briefly reviewed in Pereira et al. (2016), we can give a clear meaning to the dissipation field as depicted in KO62 while proposing a stochastic representation of it. At a given finite Reynolds number, implying a finite Kolmogorov length scale $\eta$, this representation reads

[TABLE]

where $X_{\eta}$ is a homogeneous and isotropic Gaussian field whose variance diverges logarithmically with $\eta$ and is logarithmically correlated over the integral length scale $L$. Being Gaussian, the field $X_{\eta}$ can be conveniently expressed as a linear operation (i.e. a convolution) on a white measure as

[TABLE]

where a regularized norm $|\boldsymbol{x}|_{\eta}$ is introduced, and whose exact expression has little importance at this stage, since the statistical properties of $X_{\eta}$ become independent on its precise form in the limit of infinite Reynolds number (i.e. $\eta\to 0$). We will take for instance, to illustrate our purpose, the homogeneous and isotropic small-scale cut-off $|\boldsymbol{x}|^{2}=|\boldsymbol{x}|^{2}+\eta^{2}$ for any $\boldsymbol{x}\in\mathbb{R}^{3}$. See Rhodes & Vargas (2014); Pereira et al. (2016) for a discussion on this subject.

It is then possible to show (see for example Pereira et al. 2016 for detailed computations and Rhodes & Vargas 2014 for a general approach) that

[TABLE]

Furthermore, locally, we have the convergence, for $|\boldsymbol{h}|>0$,

[TABLE]

where $g$ is a continuous and bounded function of its argument, and $\log^{+}(|\boldsymbol{x}|)=\log(|\boldsymbol{x}|)$ for $|\boldsymbol{x}|\geq 1$ and vanishes elsewhere. Doing so, we can derive the statistical properties of the modelled (Eulerian) pseudo-dissipation field (10) that fulfils the axiomatics of KO62, namely

[TABLE]

Recalling that the Kolmogorov scale $\eta$ behaves with the Reynolds number $\mathcal{R}_{e}$ as $\eta\propto L\mathcal{R}_{e}^{-3/4}$, we see that modelling the pseudo-dissipation as a multiplicative chaos (10) gives moments proportional to a power of the Reynolds number (14) as required by KO62. Furthermore, the locally averaged pseudo-dissipation

[TABLE]

is such that its moments remain bounded when $\eta\to 0$ and we have the following behaviour at small scales

[TABLE]

where $C_{q}$ is a scale and Reynolds number independent constant, that can be eventually calculated.

We can see from the proposition of Mandelbrot (1972) that taking the pseudo-dissipation (10) as the exponential of a log-correlated field (11) ensures that its moments behave as power laws of the Reynolds number (14), in particular its average is consistent with Kolmogorov’s prediction (7). Moreover, moments of its locally averaged version over a ball of radius $\ell$ (15) become independent of the Reynolds number and behave as power laws at small scales (16).

In the Lagrangian context of velocity gradient dynamics we are interested in (equation (2)), we would like to build up a unidimensional version of the multiplicative chaos used in the Eulerian framework (10) to model pseudo-dissipation as seen by a fluid particle along its trajectory. It has to be stationary, its variance particularly should be time independent, and, furthermore, it has to be causal. Causality is crucial here, not only as a physical requirement but, technically speaking, because the so-obtained multiplicative chaos has to be defined as a solution of a stochastic differential equation that will be coupled to the dynamics of $\mathsfbi{A}$ (3), as done by Girimaji & Pope (1990) while imposing that $\mbox{tr}(\mathsfbi{A}\mathsfbi{A}^{\top})$ follows the statistics of an exponentiated Ornstein-Uhlenbeck process, as already depicted in section 2.1. Thus, instead of the three-dimensional spatial field considered in (10), we will consider the following one-dimensional temporal stochastic field

[TABLE]

where now $X_{\tau_{\eta}}(t)$ is a zero-average causal and stationary Gaussian field whose variance $\mathbb{E}[X_{\tau_{\eta}}^{2}]$ blows up logarithmically with the Kolmogorov time scale $\tau_{\eta}$ (instead of (12)), and whose covariance is independent of the Reynolds number and decreases logarithmically over the integral time scale $T$ of turbulence (in contrast to (13)). As proved in appendix B, the unique solution $X_{\tau_{\eta}}(t)$ of the following stochastic dynamics fulfils these two requirements

[TABLE]

being $W$ a Gaussian white noise and $\beta_{\tau_{\eta}}(s)$ a random function known in the literature as a (${\tau_{\eta}}$-regularized) fractional Gaussian noise of Hurst exponent $H=0$ (Mandelbrot & Van Ness, 1968)

[TABLE]

We note that (18) is linear and so can be integrated in a simple fashion such that the solution $X_{\tau_{\eta}}(t)$ can be equivalently defined as

[TABLE]

where we have chosen formally the initial condition $X_{\tau_{\eta}}(-\infty)=0$. As we can see, the process $X_{\tau_{\eta}}$ is the sum of two Gaussian causal fields. The first process on the RHS of (20) is similar, in a certain sense, to the one in the formulation of the Eulerian modelling of pseudo-dissipation (11). It can be seen as an exponentially smoothed version of a fractional Brownian noise of vanishing Hurst exponent, and is known in the literature as a fractional Ornstein Uhlenbeck process (Cheridito et al., 2003). The second process on the RHS of (20) is a usual Ornstein-Uhlenbeck process whose variance blows up as $1/\tau_{\eta}$. Notice that the same white measure $W$ appears in both these Gaussian fields, making them correlated. This correlation will be responsible for a strong cancellation taking place in between these two processes. Actually, this cancellation, formalized in appendix B, will make the variance of $X_{\tau_{\eta}}\sim\log(T/\tau_{\eta})$ way smaller than the one of each of these formerly described Gaussian processes. The idea of adding a correlated process to a fractional Gaussian noise can be understood in practice with the heuristics developed in Schmitt (2003). Our contribution to this matter is to properly regularize the dynamics proposed by Schmitt (2003), injecting the large scale $T$ dependence not in the fractional Brownian noise, but through a linear damping, allowing us to rigorously define $X_{\tau_{\eta}}$ (appendix B).

It is then possible to show that the pseudo-dissipation defined as the exponential of $X_{\tau_{\eta}}$ fulfils the KO62 requirements, namely, high order moments behave as a power laws of the Reynolds number

[TABLE]

where the multiplicative constant (given explicitly in (84)) is independent of the Reynolds number, and moments of its coarse-grained version along the trajectory (9) behaving as power laws of the time scale $\tau$ in the asymptotic limit of infinite Reynolds number, according to

[TABLE]

where again the exact expression of the multiplicative constant which is both scale and Reynolds number independent is given in (B.4). The constant $\tilde{g}(0)$ has the same status as $g(0)$ in (14) and is precisely defined alongside detailed computations in appendix B.4.

2.3 Differences between stochastic models of pseudo-dissipation and comparison with numerical flows

As we have seen, modelling pseudo-dissipation as an exponentiated Ornstein-Uhlenbeck process, as done by Girimaji & Pope (1990) and discussed in section 2.1, or as a multifractal measure, like discussed in Borgas (1993) and considered by Schmitt (2003), who proposed in the context of multiplicative chaos a causal dynamics based on heuristics arguments that we defined rigorously in section 2.2, are very different matters. First of all, the first model is not consistent with the multifractal picture, which asks for (i) moments of $\varphi$ to behave as power laws of the Reynolds number, and (ii) moments of $\varphi_{\tau}$ to behave, in the asymptotic limit of infinite Reynolds number, as a power laws of the coarse-graining scale $\tau$. We have seen that indeed the second model, based on the multiplicative chaos, results in such statistical behaviour.

One may wonder whether these two different statistical behaviours could be tested against experiments and/or numerical simulations. Testing the behaviour of high order moments can be difficult since the different scalings are rather similar, and moreover, statistical convergence may be an issue. Instead, we could check the behaviour of the correlation function of the logarithm of pseudo-dissipation. Indeed, these two stochastic models have different predictions for the correlation of $\log\varphi$ along trajectories. On the one hand, the model of Girimaji & Pope (1990) predicts that $\log\varphi$ is an Ornstein-Uhlenbeck process and so $\mathbb{E}\left[\log\varphi(t)\log\varphi(t+\tau)\right]$ should decrease exponentially in time (over the integral time scale $T$). On the other hand, multifractal modelling asks for a logarithmic decrease over the time scale $T$. The analysis of Huang & Schmitt (2014), based on Lagrangian trajectories extracted from a numerical simulation among the highest Reynolds numbers available from modern computers (with Taylor based Reynolds number on the order of 400) shows that data, and more precisely the correlation function of $\log\varphi$, are consistent both with exponential and logarithmic decreases. In other words, even at the highest accessible Reynolds numbers it is very difficult to make a distinction between the two models, both being able to reproduce in a quantitative way numerical data. In the following, we will prefer to work with the multifractal model, because (i) it has a proper behaviour in the asymptotic limit of infinite Reynolds numbers and (ii) the free parameter of the models, $\mu^{l}$, which appears in both models (see (55) for the first model, and directly in (17) for the second model) has a clear meaning only in a multifractal framework, where it is denominated the intermittency coefficient.

3 Stochastic dynamics for the velocity gradient tensor with multifractal properties

Let us now turn back to the construction of a stochastic model of the velocity gradient tensor $\mathsfbi{A}$. To do so, let us first recall (3), the general definition of a stochastic dynamics for $\mathsfbi{A}$ along a Lagrangian trajectory, namely

[TABLE]

where $V_{ij}$ and $D_{ijkl}$ are respectively called the drift and diffusion terms, and $W_{ij}$ are the components of a tensorial white noise $\mathsfbi{W}$, i.e. uncorrelated in time and with independent elements (not to be confused with the scalar noise, denoted simply by a serifed letter $W$). At a fixed time $t$, they are a priori causal functionals of the tensorial functions $\mathsfbi{A}(s)_{s\in(-\infty,t]}$ and the tensorial Wiener process $\mathsfbi{W}(s)_{s\in(-\infty,t)}$. Following the approach of Girimaji & Pope (1990), we want to determine $V_{ij}$ and $D_{ijkl}$ in order to (i) include some of the crucial physics of the Navier-Stokes equations (2) that we are able to close in terms of $\mathsfbi{A}(t)$ and (ii) impose for the respective modelled pseudo-dissipation $\varphi=\mbox{tr}(\mathsfbi{A}\mathsfbi{A}^{\top})$ the multifractal structure described in section 2.2.

Let us first remark that to constrain the dynamics of $\mathsfbi{A}$ while imposing on one of its contractions $\varphi$ a certain precise statistical behaviour requires some approximations. Indeed, constraining $\varphi$, which is a scalar, does not determine uniquely $\mathsfbi{A}$, which is a tensor, since several degrees of freedom are lost in the way. Similarly, we cannot at this stage start from (22) and impose exactly for the respective pseudo-dissipation to be a multifractal process as described in section 2.2. The reason is that this multifractal process is defined with a simple scalar white noise $W$ (see (17) and (20)), whereas the dynamics of $\mathsfbi{A}$ is defined with a tensorial white noise $\mathsfbi{W}$ (22). In the sequel, we will follow the ideas developed in Girimaji & Pope (1990) to propose approximations that are consistent with underlying isotropic conditions and perform a mean-square estimation to determine the diffusion terms.

3.1 General formulation and isotropic conditions.

We consider (22) as a modelled stochastic evolution of the velocity gradient tensor $\mathsfbi{A}$. We are also interested in the respective dynamics of the pseudo-dissipation $\varphi=\mbox{tr}(\mathsfbi{A}\mathsfbi{A}^{\top})$ which is obtained from (22) after applying the chain rule (Ito’s lemma). Doing so, we obtain

[TABLE]

Assume now that the diffusion term $D_{ijkl}$ has a isotropic form, namely

[TABLE]

where $a$, $b$, $c$ are unknown and, in general, scalar functionals of $\mathsfbi{A}$ and time. The incompressibility condition imposes

[TABLE]

Computing the spurious drift term $D_{ijkl}D_{ijkl}$ entering the dynamics of $\varphi$ (23) and using the incompressibility condition $\mbox{tr}(\mathsfbi{A})=A_{ii}=0$ results in the following dynamics for $\varphi$

[TABLE]

3.2 Mean-square estimation of the diffusion term.

We want to impose a dynamics for $\varphi$ (26) as close as possible to the one of a multiplicative chaos, derived in appendix B.4 and expressed in (83), which we recall for convenience

[TABLE]

This can be done at the cost of an approximation since, as explained, the true dynamics of $\varphi$ (26), comprising a tensorial white noise $\mathsfbi{W}$, involves more degrees of freedom than the multiplicative chaos one (3.2), where only a scalar white noise $W$ appears.

We follow the approximation of Girimaji & Pope (1990) that consists in adopting a mean-square approximation of the diffusion term. We thus take the diffusion term from (23) such that its variance conditioned to $\mathsfbi{A}$ equals the variance of the diffusion term from (3.2) conditioned to $\varphi_{\tau_{\eta}}=\varphi=A_{ij}A_{ij}$. Since we are dealing with random distributions $\mathsfbi{W}$ and $W$, the meaning of their square is a priori not clear. Nonetheless, this can be clarified while considering the square of their integral over a finite range and then taking a limit as depicted in Girimaji & Pope (1990). Doing so, we obtain, for any time $t$,

[TABLE]

As in Girimaji & Pope (1990), we remark that the equality in (28) provides an underdetermined problem in relation to the free parameters $b$ and $c$ (one equation for two unknowns). To this regard, we assume that the dynamics of $\varphi=\mbox{tr}(\mathsfbi{A}\mathsfbi{A}^{\top})$ only depends on itself, and not on the contraction $\mbox{tr}(\mathsfbi{A}^{2})$. Using furthermore the incompressibility condition (25) we arrive at

[TABLE]

It would be interesting to study the dependence of the model on the choice $c\neq 0$, particularly if we want to build up separately dynamics for dissipation $\mbox{tr}(\mathsfbi{S}^{2})$ and enstrophy $-\mbox{tr}(\mathsfbi{\Omega}^{2})$. We keep these developments for further investigations. The dynamics of $\varphi$ (26) becomes

[TABLE]

We can see that this mean-square procedure suggests the replacement of the scalar white noise $W$ in the dynamics of $\varphi_{\tau_{\eta}}$ (3.2) by the scalar noise $A_{kl}W_{kl}/\sqrt{\varphi}$ that emerges from the $\varphi$ dynamics (30).

3.3 Determination of the drift term and introduction of the Recent Fluid Deformation closure

Analogously to Girimaji & Pope (1990), we would like to impose on $\varphi$ (30) some dynamics as close as possible to the one fulfilled by the multiplicative chaos (3.2). To do so, we need to provide an expression for the tensor $\mathsfbi{V}$ of the drift part of (30). This tensor is the only one able to give a flavour of the Navier-Stokes equation (2), and should be closed in terms of $\mathsfbi{A}$. The very first idea would be to take into account only the self-stretching term $\mathsfbi{V}=\mathsfbi{V}^{\mbox{\tiny{RE}}}$ (4). This would be done disregarding all the interesting physics of the joint action of the pressure Hessian and viscous term which compose the true dynamics of $\mathsfbi{A}$ (2). So instead, we will adopt the more sophisticated closure $\mathsfbi{V}=\mathsfbi{V}^{\mbox{\tiny{RFD}}}$ (5) given by the Recent Fluid Deformation (RFD) approximation and discussed in Chevillard & Meneveau (2006); Chevillard et al. (2008); Meneveau (2011). Although it reproduces some of the physics governed by the pressure and viscosity terms, this closure is not consistent with the constraint we want to impose for the pseudo-dissipation (3.2) when inserted in (30). To make it so, we will furthermore consider an additional damping term proportional to $\mathsfbi{A}$ and show that the associated multiplicative factor can be chosen such that the resulting dynamics for $\varphi$ is similar, in a certain sense, to the one fulfilled by the multiplicative chaos (3.2). This is similar to the procedure made by Girimaji & Pope (1990). Accordingly, we choose for the drift tensorial term $\mathsfbi{V}$ (22) the following model

[TABLE]

where the multiplicative factor $f(t)$ is chosen so that the respective drift term of the dynamics of $\varphi$ is as close as possible to the one of the multiplicative chaos (3.2). With this aim, and to find an expression for $f(t)$, we must identify the drift parts of the dynamics of $\varphi_{\tau_{\eta}}$ (3.2) with the those of the $\varphi$ dynamics (30).

As we have seen in section 3.2, the application of a mean-square procedure to the diffusion term entering in the dynamics of $\varphi$ (30) suggests to replace the scalar noise $W$ by $A_{kl}W_{kl}/\sqrt{\varphi}$. We will follow this suggestion and replace accordingly the fractional noise $\beta_{\tau_{\eta}}$ entering (3.2) and defined in (19), built from the scalar noise $W$, by its respective version $\hat{\beta}_{\tau_{\eta}}$ built on the noise $A_{kl}W_{kl}\sqrt{\varphi}$. This new noise reads

[TABLE]

Several reasons pushed us to consider this noise $\hat{\beta}_{\tau_{\eta}}$. First of all, as mentioned, the mean square procedure used in section 3.2 to determine the diffusion term $D_{ijkl}$ suggests the replacement of $W$ by $A_{kl}W_{kl}/\sqrt{\varphi}$. One could try to replace the scalar noise $W$ by a another scalar such as $\mbox{tr}(\mathsfbi{W})$, but in doing so one loses the correlation between $\hat{\beta}_{\tau_{\eta}}$ and the random term $A_{kl}W_{kl}$ from the dynamics of $\varphi$ (30). This correlation, as shown in the appendix B, is crucial in order to achieve multifractal statistics. One could then wonder that the statistics of $\hat{\beta}_{\tau_{\eta}}$ might be completely different from the initial Gaussian noise $\beta_{\tau_{\eta}}$ that forms the multiplicative chaos dynamics (3.2) and leads to exact multifractal properties of the positive field $\varphi_{\tau_{\eta}}$. It is to be checked numerically, a posteriori, that their imprints on the final dynamics have similar consequences. And indeed, as we verify in section 5, dedicated to numerical results, the model we are constructing with the noise $\hat{\beta}_{\tau_{\eta}}$ culminate in a pseudo-dissipation whose statistics is remarkably similar to the exact multifractal model $\varphi_{\tau_{\eta}}$ from (3.2). This being said, it is possible to verify numerically that the variances of $\hat{\beta}_{\tau_{\eta}}$ and $\beta_{\tau_{\eta}}$ are actually the same, i.e.

[TABLE]

Thus, as far as the variance is concerned, the noises $\hat{\beta}_{\tau_{\eta}}$ (32) and $\beta_{\tau_{\eta}}$ (19) behave similarly. Higher order moments will be different though, since $\beta_{\tau_{\eta}}$ is a Gaussian process and $\hat{\beta}_{\tau_{\eta}}$ is not.

Again analogously to Girimaji & Pope (1990), we equate the drift term of the multifractal process (3.2) with the drift term of pseudo-dissipation (30) using as a model for $\mathsfbi{V}$ the sum of the RFD closure’s drift term and a (non-linear) damping term (31), while replacing the noise $\beta_{\tau_{\eta}}$ by the new noise $\hat{\beta}_{\tau_{\eta}}$, to deduce an expression for the functional $f(t)$

[TABLE]

Let us note that the quantity $\mathbb{E}\left[\left(X_{\tau_{\eta}}\right)^{2}\right]$, a function of $\tau_{\eta}/T$ only (see (61)), already present in the dynamics of $\varphi_{\tau_{\eta}}$ (3.2), also appears in (34).

3.4 Assembling the model

By inserting the various terms discussed above in the dynamics of $\mathsfbi{A}$ (22), we forge a closed dynamics along a trajectory. For clarity and convenience, we recall here all the pieces gathered to build the model. For the drift term (31), $\mathsfbi{V}=\mathsfbi{V}^{\mbox{\tiny{RFD}}}+f(t)\mathsfbi{A}$, there are two contributions. The RFD part is given by (5), that is,

[TABLE]

being $\mathsfbi{C}_{\tau_{\eta}}$ the short-time Cauchy-Green tensor (6) and $\tau_{\eta}$ and $T$ the Kolmogorov and integral time scales respectively. The damping part is formed with the scalar function $f(t)$ derived just above in (34), which comprises the intermittency coefficient $\mu^{l}$, the quantity $\mathbb{E}\left[\left(X_{\tau_{\eta}}\right)^{2}\right]$, which is function of $\tau_{\eta}/T$ only, and justified in (61), that reads

[TABLE]

and can be eventually expressed with special functions, as given by a symbolic calculation software. A modification to the regularized fractional Gaussian noise of vanishing Hurst exponent (19) denoted as $\hat{\beta}_{\tau_{\eta}}$ (32), which reads

[TABLE]

enters also in the definition of the scalar function $f(t)$. Together with the isotropic diffusion term $D_{ijkl}$ (equations 24 and 29), constructed by a mean-square estimation procedure (section 3.2), one finally arrives at

[TABLE]

This proposed dynamics is closed. It yields, starting for instance at time $t=0$, the time evolution of $\mathsfbi{A}(t)_{t\geq 0}$ once some history $\mathsfbi{A}(t)_{t<0}$ is given. This boundary condition is required by the modification of the fractional Gaussian noise $\hat{\beta}_{\tau_{\eta}}$ (37) that seeks in the past to set its current value. In this sense the dynamics proposed in (38) is non-Markovian. We expect this noise to decorrelate as fast as $1/t^{2}$ because its Gaussian version $\beta_{\tau_{\eta}}$ does decorrelate this way (see developments in appendix B). So in practice, the noise $\hat{\beta}_{\tau_{\eta}}(t)$ (37) is well approximated by a truncation over a finite interval, say for instance $t\in[t-T,t]$, and this approximation is more accurate as the small time scale $\tau_{\eta}$ gets smaller compared to $T$.

4 Numerical procedure

This section is devoted to the numerical integration of the proposed new dynamics for the velocity gradient tensor $\mathsfbi{A}$ along a trajectory (38). For this purpose, we need a numerical approximation for the noise $\hat{\beta}_{\tau_{\eta}}[t]$ (32) which composes the function $f(t)$ (34), as well as its boundary condition. Denoting as $\Delta t$ a time marching step, the tensor $\mathsfbi{W}[dt]$ is such that each element $W_{ij}[dt]$ at time $t$ is an independent realization of a zero-mean Gaussian random variable with variance $\Delta t$.

Over a discrete set of time and starting, say, at time $t=0$, we will use henceforth the following numerical approximation of the noise $\hat{\beta}_{\tau_{\eta}}$ (32):

[TABLE]

This approximation, based on a truncation of the time integration, is, as mentioned in the end of section 3.3, more and more realistic as $\tau_{\eta}/T\to 0$. A numerical study devoted to the dependence of the approximation of the fractional Gaussian noise (32) by its truncated version (39) remains to be properly done, since this is the most demanding step of the numerical integration of the $\mathsfbi{A}$ dynamics. Nevertheless, we have performed simulations (data not shown) using a truncation over a shorter range $[t-T/2,t]$ and no relevant quantitative differences were observed on the statistical results, at least in the investigated range of $\tau_{\eta}$ values. To make sense of the discretization used in (39), we also need to choose $\Delta t$ small enough compared to $\tau_{\eta}$ to properly resolve the smooth kernel $(t-s+{\tau_{\eta}})^{-3/2}$ when $s\to t^{-}$.

Having set a numerical approximation of the noise $\hat{\beta}_{\tau_{\eta}}$ (39), we need now to give initial conditions. In the following numerical simulations, we will start at time $t=0$, with $\mathsfbi{A}[0]=\mathsfbi{\mathcal{N}}-\mbox{tr}(\mathsfbi{\mathcal{N}})\mathsfbi{I}/3$, where $\mathsfbi{\mathcal{N}}$ is a $3\times 3$ tensor such that elements are independent and normally distributed with zero average and unit variance. As boundary conditions, we simply take $\mathsfbi{A}[t]_{-T\leq t\leq 0}=\mathsfbi{A}[0]$. The precise forms of the initial and boundary conditions are not important (as long at they are trace-free) since we expect to reach a stationary regime independent of them at large time $t$. Finally, we use for time marching a simple Euler approximation, which is a delicate matter and deserve some comments. Even the multiplicative chaos dynamics (3.2) with the fractional Gaussian noise (19) exhibits a rather slow statistical convergence with respect to the ensemble size, as already observed for the Eulerian case in Pereira et al. (2016), where many realizations of high resolution 3d cubes were needed to assert their conclusions. Thus very large ensembles will be necessary to attain smooth statistics. With the Euler method, small time steps are required to reduce the numerical error associated with time integration, meaning that a fairly big amount of data will be needed, in particular for large Reynolds numbers, but this may translate into a substantial error accumulation towards the end of the integration. Unfortunately, for the proposed dynamics (38), the highly non-linear and non-Markovian character, together with the multiplicative nature of the noise and the strong correlations between the terms, make the implementation of higher order time integration methods a challenging problem. All things considered, we employ the Euler method but instead of producing very long trajectories we choose to generate a large ensemble of moderately shorter trajectories (always respecting an initial transient time). Also, for the larger Reynolds numbers, runs with varying $\Delta t$ revealed decisive to determine if discrepancies among the statistical behaviours of $\varphi$ built from (38) and $\varphi_{\tau}$ from (3.2) were caused by numerical errors or fundamental differences between the processes.

5 Numerical results

5.1 Parameters of the simulations

In the following simulations we rescale time by $T$ and take $T=1$, with no loss of generality. As for the intermittency coefficient in the Lagrangian framework $\mu^{l}$, we use the recently estimated value on high Reynolds number numerical turbulent data of Huang & Schmitt (2014), who found $\mu^{l}=0.3$. Notice that this is consistent with the commonly accepted value of the intermittency coefficient in the Eulerian framework $\mu^{e}=0.2$, as it may be checked with the phenomenological theory of Borgas (1993), which relates $\mu^{l}$ and $\mu^{e}$ through a non-linear relation (relating more precisely the respective singular spectra, using the vocabulary of the multifractal formalism). If one neglects the fluctuating nature of the spatial and temporal dissipative scales, an approximative linear relation can be drawn in a straightforward way: we observe that the moments of the pseudo-dissipation field $\varphi$ can be expressed in the Eulerian (14) and Lagrangian (21) frameworks, and, recalling that $\eta/L\propto\mathcal{R}_{e}^{-3/4}$ and $\tau_{\eta}/T\propto\mathcal{R}_{e}^{-1/2}$, one is led to $3\mu^{e}=2\mu^{l}$ by equating their Reynolds number dependence.

The positive quantity $\mathbb{E}[(X_{\tau_{\eta}})^{2}]$ contained in the dynamics of $\mathsfbi{A}$ (38) is given explicitly as an integral in (61). This integral could be numerically evaluated, but it may be analytically expressed in terms of error and generalized hypergeometric functions, which in turn must be numerically evaluated as well, but whose handling by popular math softwares can be quite efficient and controllable. We follow this path to compute $\mathbb{E}[\left(X_{\tau_{\eta}}\right)^{2}]$ for the various values of $\tau_{\eta}/T$ considered.

Using $\Delta t=2\times 10^{-3}\tau_{\eta}$, or $\Delta t=10^{-3}\tau_{\eta}$ for higher $\mathcal{R}_{e}$, we numerically integrate numerous realizations of the $\mathsfbi{A}$ dynamics over time intervals of dozens or hundreds of $T$, such that at least a total time of $3\times 10^{3}\,T$ is covered. We observe then, for the explored range $2.48\times 10^{-3}\leq\tau_{\eta}/T\leq 2.64\times 10^{-1}$, or equivalently $-6.0\leq\log(\tau_{\eta}/T)\leq-1.33$, the stationary statistics of $\mathsfbi{A}$ and its contraction $\varphi=A_{ij}A_{ij}$.

It is worth mentioning that extensive tests were successfully done in order to establish not only the stationarity of the generated processes but also its independence of the initial conditions. In this perspective, we checked (data not shown) that an ensemble of highly anisotropic initial conditions indeed relaxes to an isotropic state, at least with the same level of success as in subsection 5.3.

5.2 Variance and covariance of the velocity gradients, and the Reynolds number dependence

5.2.1 Variance

For illustrative purposes, let us first begin by showing some temporal profiles of a diagonal element $A_{11}$ (figure 1a), an off-diagonal element $A_{12}$ (figure 1b) and the related pseudo-dissipation $\varphi$ (figure 1c) along a Lagrangian trajectory, at a given Reynolds number quantified by $\log(\tau_{\eta}/T)=4.67$. We observe an erratic behaviour of gradients that is expected in a highly turbulent flow. We may also notice the non-symmetrical behaviour of the diagonal element $A_{11}$ fluctuations (figure 1a), reminiscent of the skewness phenomenon, and the highly intermittent $\varphi$ fluctuations. We invite the reader to take a look at section 5.6 for a more in-depth study of the non-Gaussian nature of these fluctuations.

The second behaviour we would like to check is whether the proposed stochastic model (38) is consistent with basic Kolmogorov phenomenology, such as the Reynolds number (or $\tau_{\eta}/T$) dependence of the variance of the elements of $\mathsfbi{A}$. In particular, we would like the average pseudo-dissipation $\mathbb{E}(\varphi)=\mathbb{E}[\mbox{tr}(\mathsfbi{A}\mathsfbi{A}^{\top})]$ to be proportional to $(T/\tau_{\eta})^{2}$, according to (7). Let us notice that we expect that since we are somehow putting this property by hands. Indeed, we use as a constraint for building the dynamics of $\mathsfbi{A}$ (38) that it fulfils as close as possible the dynamics of a multiplicative chaos $\varphi_{\eta}$ (3.2), which is such that $\mathbb{E}(\varphi_{\eta})=1/\tau_{\eta}^{2}$. It remains to check numerically that this behaviour is actually observed for $\varphi$. More precisely, we would like to check whether, as required by isotropy (see for instance Pope, 2000),

[TABLE]

We represent in figure 1(d) in a logarithmic scale the behaviour of the variances of components $A_{11}$ and $A_{12}$ and the average of $\varphi$ as a function of the ratio $\tau_{\eta}/T$. We see that the isotropic relations (40) are verified on our numerical simulations of the process $\mathsfbi{A}$ 38 for a wide interval of $\tau_{\eta}/T$. This shows that, in this sense, statistics of $\mathsfbi{A}$ are isotropic and behave in a consistent way with $\tau_{\eta}/T$.

5.2.2 Autocovariance

As far as second order statistics are concerned, another important property that should be checked in the simulated trajectories of $\mathsfbi{A}$ is the covariance of its elements. Kolmogorov’s phenomenology (Frisch, 1995) suggests that, along Lagrangian trajectories, elements of $\mathsfbi{A}$ (such as $A_{11}$ or $A_{12}$) are correlated over the small time scale $\tau_{\eta}$. We remark that if indeed this short-time correlation is observed in the present stochastic model, it is not clear to us how to derive it directly from the dynamics given in (38) (see devoted works on this subject by Afonso & Meneveau, 2010; Yu & Meneveau, 2010) since we are only imposing that pseudo-dissipation (which is a positive quantity) is correlated over the large time scale $T$. We will come back to this long-range correlation behaviour in the next section. As we are dealing with a highly nonlinear dynamics and not imposing in a rigorous and definitive way that the elements of $\mathsfbi{A}$ should be correlated over the short time scale $\tau_{\eta}$, it deserves to be precisely quantified in our simulations. We represent in figure 2 the estimation of the correlation function of the diagonal element $A_{11}$ as a function of (i) the rescaled time $\tau/T$ (figure 2(a)), and (ii) the rescaled time $\tau/\tau_{\eta}$ (figure 2(b)). Since the variance of $A_{11}$ increases as $\tau_{\eta}/T$ decreases, we focus on the normalized covariance, namely

[TABLE]

which is independent of $t$ by stationarity, and tends towards unity as $\tau\to 0$. Note also that we do not subtract the square of the average of the diagonal element since we have $\mathbb{E}(A_{11})=0$. A comparison between figures 2(a) and 2(c) shows clearly that the diagonal element $A_{11}$ is indeed correlated over the Kolmogorov’s time scale $\tau_{\eta}$. Similar observations can be made on the off-diagonal element $A_{12}$ (data not shown).

As we have seen, the diagonal element $A_{11}$ is correlated over the small time scale $\tau_{\eta}$, as it should be. One may wonder what is the characteristic correlation time scale of the amplitude of the diagonal element $A_{11}$, or equivalently the absolute value of the diagonal element $|A_{11}|$. To do so, we define the following (normalized and centred) autocorrelation function

[TABLE]

which, by stationarity, should be time independent.

We represent in figure 2 the behaviour of this correlation function (42) both as a function of $\log(\tau/T)$ (figure 2(b)) and $\log(\tau/\tau_{\eta})$ (figure 2(d)). We indeed observe two type of behaviours. When the lag $\tau$ is rescaled by the large scale $T$, figure 2(b) reveals that the short time portion of the correlations, typically when $\tau$ is of the order of $\tau_{\eta}$, does not coincide, whereas for large times (i.e. when $\tau\gg\tau_{\eta}$) the correlations vanish around the large scale $T$. This analysis is confirmed in figure 2(d) where it is seen that the correlations superimpose over the short time scale $\tau_{\eta}$ but not over larger time lags.

This study shows that the (signed) diagonal elements decorrelate over the short time scale $\tau_{\eta}$, whereas their absolute values decorrelate over the large time scale $T$. This is a hallmark of intermittency: the sign of fluctuations decorrelates very fast, whereas the amplitude remains correlated over the large time scale $T$. This is consistent with saying that intermittency and high values of gradients appear as bursts, that can be correlated over possibly large time scales (Frisch, 1995).

5.3 Comparisons between the dynamics of $\varphi$ and $\varphi_{\tau_{\eta}}$

In this section, we would like to study numerically the statistical properties of the pseudo-dissipation $\varphi=\mbox{tr}(\mathsfbi{A}\mathsfbi{A}^{\top})$ as obtained from the stochastic process (38), which are difficult to get analytically. We are also interested in comparing them to the ones given by the model $\varphi_{\tau_{\eta}}$ (3.2), that can be carried on analytically (appendix B). The statistics of $\varphi$ and $\varphi_{\tau_{\eta}}$ are indeed expected to be very similar, since their dynamics are themselves driven by similar equations. To justify this, compare the dynamics of $\varphi_{\tau_{\eta}}$ (3.2) to that of $\varphi$, which can be obtained while inserting (34) in the dynamics of $\mathsfbi{A}$ (38) and applying the chain rule, which reads

[TABLE]

We see that (3.2) and (5.3) are different for two reasons: (i) the noise $\beta_{\tau_{\eta}}$ (19) is replaced by the noise $\hat{\beta}_{\tau_{\eta}}$ (33) and (ii) the noise $\varphi_{\tau_{\eta}}W$ is replaced by $\sqrt{\varphi(t)}A_{kl}W_{kl}$.

Continuing our analysis, we represent in figure 3(a) the estimation of the PDFs of the realizations of the random variables $\log\varphi_{\tau_{\eta}}$ and $\log\varphi$, once centred (we subtract their respective means) and normalized (we divide them by their respective standard deviations) for different values of the ratio $\tau_{\eta}/T$. We recall that $\varphi_{\tau_{\eta}}$ is obtained while exponentiating $X_{\tau_{\eta}}$, accordingly with (17), which is itself given by a linear stochastic differential equation (18). We see that, as expected, the statistics of $\log\varphi_{\tau_{\eta}}$ is indeed normal for any $\tau_{\eta}/T$. As for $\log\varphi$, its PDFs closely follow those of $\log\varphi_{\tau_{\eta}}$, and hence their fluctuations are indeed observed very similar.

As also studied by Yeung & Pope (1989); Pope (1990); Huang & Schmitt (2014), we reproduce in figure 3(b) a similar estimation of the PDFs of two other quantities derived from $\mathsfbi{A}$, namely dissipation $\mbox{tr}(\mathsfbi{S}\mathsfbi{S}^{\top})$ and enstrophy $\mbox{tr}(\mathsfbi{\Omega}\mathsfbi{\Omega}^{\top})$, being $\mathsfbi{S}$ and $\mathsfbi{\Omega}$ respectively the symmetric and antisymmetric parts of the decomposition of $\mathsfbi{A}$. We see that whereas the statistics of $\varphi$ is very close to the one of a lognormal process, those of dissipation and enstrophy differ significantly from such a process. This fact holds also true in direct numerical simulations (DNS) of turbulence as shown in Yeung & Pope (1989); Pope (1990); Huang & Schmitt (2014).

Let us now focus on the averages of these fields. We represent in figure 3(c) the average values of the pseudo-dissipation fields $\varphi$ and $\varphi_{\tau_{\eta}}$. As expected, since this is imposed by hands, it is observed that $\tau_{\eta}^{2}\mathbb{E}[\varphi_{\tau_{\eta}}]\approx 1$ for all studied values of $\tau_{\eta}/T$. Up to statistical fluctuations, we also observe that $\mathbb{E}[\varphi]\approx\mathbb{E}[\varphi_{\tau_{\eta}}]$ for all $\tau_{\eta}/T$. Here, numerous tests with smaller time steps and bigger ensembles were necessary to establish this result, since systematic differences may appear at high Reynolds numbers if statistical convergence is not properly attained. In fact, the mild increases in the averages one observes in figure 3(c) for the smaller ratios $\tau/T$ (higher Reynolds) seem to be mainly due to numerical errors, since they are much more acute if greater time steps are used for integration or smaller ensembles are considered (data not shown). We also observe, and this is a limitation of the present model, that statistics are not fully consistent with those expected from isotropic turbulence. Indeed, for a homogeneous and isotropic velocity field, we expect that $\mathbb{E}[\varphi]=2\mathbb{E}[\mbox{tr}(\mathsfbi{S}\mathsfbi{S}^{\top})]=2\mathbb{E}[\mbox{tr}(\mathsfbi{\Omega}\mathsfbi{\Omega}^{\top})]$ (Pope, 2000), a relation which is not perfectly met by the model. Figure 3(c) shows that the average enstrophy is a bit too small, with $2\mathbb{E}[\mbox{tr}(\mathsfbi{\Omega}\mathsfbi{\Omega}^{\top})]\approx 0.88\mathbb{E}[\varphi]$, while the average dissipation a bit too high, $2\mathbb{E}[\mbox{tr}(\mathsfbi{S}\mathsfbi{S}^{\top})]\approx 1.17\mathbb{E}[\varphi]$, compared to the average pseudo-dissipation. We discuss possible improvements of the model in order to correct this bias in section 6.

Next, we study in figure 3(d) the correlation structure of the fields $\varphi(t)$ and $\varphi_{\tau_{\eta}}(t)$. With this aim, we estimate the autocovariance

[TABLE]

which, by stationarity, should be time independent. We know for sure, in the asymptotic regime $\tau_{\eta}\to 0$, that when it comes to the field $\varphi_{\tau_{\eta}}$ this correlation function behaves logarithmically with the time lag $\tau$, i.e.,

[TABLE]

We compare in figure 3(d) the estimation of the autocovariance $\mathbb{E}_{c}$ (44) for the field $\varphi_{\tau_{\eta}}$ with the asymptotic logarithmic behaviour (45). It is seen that, as $\tau_{\eta}/T$ gets smaller, the autocovariance becomes more consistent with the expected logarithmic behaviour. As for $\varphi$, it is found that its corresponding centred correlation function follow very closely those of $\varphi_{\tau_{\eta}}$, for all the ratios $\tau_{\eta}/T$ here considered. We may nonetheless remark that, once again, improved time resolutions and larger ensembles were decisive in establishing this numerical trend.

5.4 Joint PDF of the two invariants $R$ and $Q$

As an important characterization of the velocity gradient tensor $A_{ij}$, we study its two non vanishing invariants. The second invariant $Q$ is given by

[TABLE]

where $\boldsymbol{\omega}$ is the vorticity vector and $\mathsfbi{S}$ the symmetric part of $\mathsfbi{A}$, and can be interpreted as the competition between enstrophy and dissipation (per unit viscosity). Therefore, positive $Q$ represents rotation-dominated regions and negative $Q$ dissipation-dominated regions. Analogously, the third invariant $R$ is given by

[TABLE]

representing competition between enstrophy production and dissipation production (we refer to Tsinober, 2001; Wallace, 2009; Meneveau, 2011, for discussions on this topic). A numerical estimation of the PDF of the invariants $Q$ and $R$ is represented in figure 4(a). One may see that it is more elongated along the right tail of the Vieillefosse line and in the upper-left quadrant, precisely as observed in direct numerical simulations.

5.5 Alignments of vorticity with the eigenframe of the deformation rate

Another striking property of turbulence is the preferential alignment of vorticity with the strain eigendirection associated to the intermediate eigenvalue. We refer again to Tsinober (2001); Wallace (2009); Meneveau (2011) for further discussions. Figure 4(b) displays the probability density function of the cosine of the angle $\theta$ between vorticity and the strain eigenvectors. It indicates the preferential alignment of vorticity with the intermediate eigendirection of deformation, as it is observed in DNS. We can add that the alignments PDFs (figure 4b) do not depend significantly on the $\tau_{\eta}/T$ value in the studied range, as well as the $RQ$ plane from figure 4(a), where the contour lines just get a bit more spaced as $\mathcal{R}_{e}$ increases due to additional non-Gaussianity (data not shown).

5.6 Higher moments of the gradients and multifractal properties

We now characterize and quantify the intermittent properties of the statistics of $\mathsfbi{A}$ in a sense that will become more precise.

First we estimate the PDF of the diagonal and off-diagonal elements of $\mathsfbi{A}$. They are portrayed respectively in figures 5(a) and 5(b). All PDFs are normalized to unit variance and are arbitrarily vertically shifted for the sake of clarity. It is clear that in both cases, as the ratio $\tau_{\eta}$ decreases, the PDF undergoes a continuous shape deformation, from a Gaussian shape at large ratios towards large tails at small ratios. It is additionally noticeable that the cores of the diagonal elements PDFs (figure 5(a)) are not symmetrical. This continuous shape deformation of the PDFs as the Reynolds number increases is a hallmark of multifractal behaviour.

To further characterize the intermittency phenomenon and the non-Gaussian behaviour of fluctuations, we compute the skewness $\mathbb{E}[A_{ij}^{3}]/(\mathbb{E}[A_{ij}^{2}])^{3/2}$ of diagonal and off-diagonal elements. It is seen from figure 5(c) that, up to statistical uncertainties, only diagonal elements are skewed, as allowed by isotropy. Moreover, we contemplate that the skewness of diagonal elements exhibits a Reynolds number dependence (through the ratio $\tau_{\eta}/T$), as expected from multifractal phenomenologies. We superimpose in figure 5(c) the power law (dotted line) observed in experimental and numerical data (see equation 4.4 in Ishihara et al., 2007)

[TABLE]

which is very close to the prediction of the multifractal formalism (see for instance Frisch, 1995; Chevillard et al., 2012). It is worth noting that in figure 5(c) the power law (48) has been fitted beforehand on a log-log plot in order to fix the suitable multiplicative constant (not shown).

In the same spirit, we pursue the quantification of the level of intermittency through an estimate of the velocity gradients flatness, i.e. $\mathbb{E}[A^{4}_{ij}]/(\mathbb{E}[A^{2}_{ij}])^{2}$, again for both diagonal and off-diagonal elements of $\mathsfbi{A}$. The resulting behaviours are shown in figure 5(d). It is evident that the depicted flatnesses exceed the Gaussian value $3$. We superimpose in this representation (dotted lines) the power law observed in numerical and experimental flows (Ishihara et al., 2007)

[TABLE]

which is also very close to what is predicted by the multifractal formalism (Frisch, 1995; Chevillard et al., 2012). For our simulated process (38), the power law (49) seems to be representative of the statistics only for the smallest values of $\tau_{\eta}/T$ tested.

We emphasize that these numerically observed results for the velocity gradients higher-order moments illustrate the self-consistency of the present approach, once its single free parameter $\mu^{l}$ is given. Indeed, the power-law exponents as given in (48) and (49) can be readily related to the intermittency coefficient $\mu^{l}$ of the Lagrangian framework since it governs the statistics of the overall amplitude of the gradients, and more precisely pseudo-dissipation. They can also be derived, in the Eulerian framework, from the respective intermittency coefficient $\mu^{e}$ (see (10)) as it is traditionally done in the literature (Frisch, 1995). We see that arguments developed by Borgas (1993), as recalled in section 5.1, which relate both intermittency coefficients as $3\mu^{e}=2\mu^{l}$, give a consistent picture of the power-laws announced in (48) and (49).

5.7 Rotation rate of anisotropic particles

Many numerical and theoretical studies employ the Jeffery equation (Jeffery, 1922) to predict the time evolution of the orientation of axisymmetric ellipsoidal particles as they are advected and influenced by a turbulent velocity field. Specifically, Jeffery’s equation for the unit orientation vector $\boldsymbol{p}(t)$ in the ellipsoid major axis reads

[TABLE]

where $\mathsfbi{S}$ and $\mathsfbi{\Omega}$ are the strain and rotation rate tensors respectively, and $\lambda=(\alpha^{2}-1)/(\alpha^{2}+1)$, being $\alpha$ the particle’s aspect ratio. In this convention, the limit $\alpha\to\infty$ corresponds to rods (elongated particles along the axis of symmetry), whereas discs are given in the limit $\alpha\to 0$. See for instance Chevillard & Meneveau (2013) or Gustavsson et al. (2014) for details, and Voth & Soldati (2017) for a recent review on the subject. It was proposed by Parsa et al. (2012) to quantify, in experiments and numerical simulations, the variance of the rotation rate of anisotropic particles $\mathbb{E}[|\dot{\boldsymbol{p}}|^{2}]$ as a function of the aspect ratio $\alpha$, which turns out to be a precise and demanding way to test models for the velocity gradient tensor $\mathsfbi{A}$ (Chevillard & Meneveau, 2013).

We reproduce in figure 6 the behaviour of the normalized variance of the rotation rate $\mathbb{E}[|\dot{\boldsymbol{p}}|^{2}]/\mathbb{E}[\varphi]$ as estimated in DNS (Li et al., 2008) and in trajectories built from the RFD approximation (see Chevillard & Meneveau, 2013, for details), along with the outcome of the present stochastic model (38) for different $\tau_{\eta}/T$. Let us recall some of the interpretations that can be drawn from the behaviour of the rotation rate as a function of the aspect ratio $\alpha$. At high aspect ratios $\alpha\gg 1$, corresponding to elongated particles along the axis of symmetry (or rods), the dynamics of anisotropic particles resembles that of material lines (see for instance Guala et al., 2005, for a review on this matter). Therefore, a strong correlation between vorticity and the orientation vector is expected, leading in particular to a preferential alignment between them, as observed in numerical simulations (Pumir & Wilkinson, 2011). At low aspect ratios $\alpha\ll 1$, corresponding to disc shaped particles in the plane perpendicular to the orientation vector, it was observed that, contrary to rods, the orientation vector gets preferentially perpendicular to vorticity (Chevillard & Meneveau, 2013; Gustavsson et al., 2014). As a consequence of this preferential alignment it is expected that discs will tumble, implying a higher rotation rate than in the rod case. This is what is observed in DNS (Parsa et al., 2012): $\mathbb{E}[|\dot{\boldsymbol{p}}|^{2}]$ increases monotonically as the aspect ratio $\alpha$ decreases from rods ($\alpha\gg 1$) to discs ($\alpha\ll 1$).

It was shown in Chevillard & Meneveau (2013) that the stochastic model for velocity gradients coined with the RFD approximation is only able to satisfactorily reproduce the behaviour of the rotation rate observed in DNS for aspect ratios larger than unity. To this regard, the present model (38) exhibits an even more striking agreement for these aspect ratios. For particles with smaller aspect ratios ($\alpha<1$), it is seen that the RFD approximation predicts a decrease of the rotation rate, whereas in DNS it keeps increasing up to a saturation value. The present model (38) shows a similar saturation structure of the fluctuations with the aspect ratio $\alpha$, an improvement in the modelling of the rotation rate of anisotropic particles compared to the RFD approximation, but overall, even if the present model (38) exhibits many aspects of the statistical properties of the velocity gradient tensor $\mathsfbi{A}$, it still fails at providing a satisfying picture of the tumbling motions of discs in turbulence. Let us finally remark that the predictions of the present model on the rotation rate are almost independent of the Reynolds number (or equivalently $\tau_{\eta}/T$), whereas it is shown in former sections that the statistics of $\mathsfbi{A}$ do have a realistic dependence. This is consistent with the trends observed from DNS studies (Parsa et al., 2012).

6 Conclusion and perspectives

We have proposed a new stochastic time evolution of the velocity gradient tensor $\mathsfbi{A}$ along a Lagrangian trajectory (38), and studied numerically the statistical properties of its solution.

The main ingredients of this model are (i) a regularization by the joint action of closures for pressure Hessian and the viscous Laplacian of the exact self-stretching term entering the dynamics (5) and (ii) an additional stochastic closure that constrains the pseudo-dissipation to follow as close as possible the dynamics of a multifractal measure. Doing so, we end up with a non-Markovian dynamics that depends on the Reynolds number (through the ratio $\tau_{\eta}/T$) and on a free parameter $\mu^{l}$, the intermittency coefficient, taken to be equal to $0.3$ as observed in DNS of homogeneous and isotropic turbulence. To our knowledge, this is the first stochastic proposition for the dynamics of A which incorporates the intermittency parameter $\mu^{l}$ and is able to reproduce, in a realistic way, the statistics of the small scales of turbulence at any Reynolds numbers without relying on hierarchical structures for the couplings of multiscale velocity gradient tensors, as it was previously addressed by Biferale et al. (2007) and, more recently, by Johnson & Meneveau (2017).

To illustrate the realism of our model, we have numerically studied the statistical properties of its solution and showed that the variance and higher-order moments of the elements of $\mathsfbi{A}$ behave accurately with the Reynolds number as observed in experimental and numerical data, and as captured by the multifractal phenomenology of turbulence. Furthermore, geometrical statistical properties such as vorticity alignments with respect to the eigenframe of deformation and joint distribution of the invariants are also reproduced in a satisfactory manner.

Despite the realism of the so-obtained statistical behaviour of $\mathsfbi{A}$, and the theoretical success to gather in a certain sense the multifractal phenomenology and the velocity gradient dynamics, it is also shown that the tumbling phenomenon of anisotropic particles of small aspect ratios (i.e. discs) is still not accurately reproduced.

The present proposition (38) could be improved in several ways. First of all, the statistical properties of $\mathsfbi{A}$ from the model are not fully consistent with those of isotropic turbulence (see the behaviour of the averaged dissipation, enstrophy and pseudo-dissipation in figure 3(c), and the corresponding discussion provided in section 5.3. To this regard, a first step toward an improvement of the model would be to explore the dependence of the closure made in (29) where the choice $c=0$ is made. Another possible way to improve these statistics would be to consider more sophisticated closures of pressure Hessian and viscous term other than the RFD approximation, such as the propositions of Wilczek & Meneveau (2014) and, more generally, those of Johnson & Meneveau (2016). In a different line of development, it would be interesting to try to adjust the model to more complex anisotropic flows, in the spirit of Girimaji & Speziale (1995) for instance.

As for the statistics of rotation rate of anisotropic particles (section 5.7), an alternative way to explore the implied physics of the tumbling phenomenon and the implications on velocity gradient statistics would be to quantify the gradients’ third order correlation functions as suggested in Gustavsson et al. (2014). Such a study, supplemented by the analysis of homogeneous and isotropic turbulence DNS, would help us to better understand the underlying mechanisms leading to the tumbling phenomenon and to design accurate stochastic models able to reproduce them. We keep these perspectives for future investigations.

Acknowledgements

The authors thank the PSMN (Pôle Scientifique de Modélisation Numérique) computing centre of ENS de Lyon for numerical resources. R.M.P. acknowledges CAPES for financial support. L.M. thanks CNPq for partial support and the ENS-Lyon for its warm hospitality and great scientific environment during his collaboration visit. L.C. is supported by ANR grants Liouville ANR-15-CE40-0013 and MultiFracs ANR-16-CE33-0020.

Appendix A Exponentiation of an Ornstein-Uhlenbeck process: definition and consequences

A.1 Definition and moments of the pseudo-dissipation field

Pope (1990) has proposed to model the instance of the pseudo-dissipation field as seen by a Lagrangian particle as

[TABLE]

where $\hat{a}$ is the free parameter of the model. The dynamics defined by (51) has a unique statistically stationary solution. Without loss of generality, one may choose some arbitrary positive initial condition and take a shortcut to the stationary regime considering formally $t_{0}=-\infty$. This unique solution reads

[TABLE]

It is easily seen from (52) that moments of $\varphi$ are time-independent (stationary regime) and, considering for simplicity $q\in\mathbb{N}$, given by

[TABLE]

In particular, Kolmogorov’s prediction (7) is automatically fulfilled in this formalism, that is, $\mathbb{E}(\varphi)=1/\tau_{\eta}^{2}$.

Let us now compute the moments of the coarse-grained pseudo-dissipation $\varphi_{\tau}$ (9) along a Lagrangian trajectory. We have, for $q\in\mathbb{N}^{*}$,

[TABLE]

Notice that

[TABLE]

where each of the stochastic integrals in the argument of the exponential is a zero-average Gaussian process of variance

[TABLE]

and covariance

[TABLE]

Using now that $\mathbb{E}(e^{g})=e^{\frac{1}{2}\mathbb{E}(g^{2})}$ for any zero-average Gaussian variable $g$, we have

[TABLE]

Rescaling the dummy variables $s_{i}$ by $\tau$ and taking the limit $\tau\to 0$ then recovers the $q^{\mbox{th}}$-order moment of $\varphi$ expressed in (53).

A.2 Discussions on the Reynolds number dependence of the free parameter of the theory

To be consistent with both observed fluctuations of turbulence and the multifractal formalism (Borgas, 1993; Frisch, 1995), we expect at the very least a power law dependence of the $\varphi$ moments (52) with the Reynolds number $\mathcal{R}_{e}$. In this spirit, Girimaji & Pope (1990) suggest to set the free parameter $\hat{a}$, up to a Reynolds independent additive constant, equal to

[TABLE]

defining thus a dimensionless free parameter $\mu^{l}$ known in the language of the multifractal formalism as the intermittency coefficient. Hence, recalling that the ratio $T/\tau_{\eta}$ is proportional to $\sqrt{\mathcal{R}_{e}}$, we see that modelling a pseudo-dissipation trajectory as the exponential of an Ornstein-Uhlenbeck process (51) predicts in a consistent way the Reynolds number dependence of the higher-order moments (53) if the free parameter is chosen with an appropriate Reynolds number dependence (55).

Let us now discuss the applicability of the Refined Similarity Hypothesis (RSH) in this context. Recall that RSH bridges in a Eulerian context the statistics of the coarse-grained dissipation with the statistics of the spatial velocity increments (Frisch, 1995). In a Lagrangian context (Borgas, 1993), but using pseudo-dissipation instead of dissipation, the RSH implies in a similar way (see Chevillard et al., 2012, for details) that the statistics of temporal velocity increments is given by that of an appropriate power of $\varphi_{\tau}$. Among other features, if the statistics of the increments is Reynolds number independent in the inertial range, as observed in real flows, this also holds for $\tau_{\eta}^{2}\varphi_{\tau}$. It is then easy to see, since the parameter $\hat{a}$ diverges with the Reynolds number (55), that the high-order moments of $\tau_{\eta}^{2}\varphi_{\tau}$ (54) will also diverge for any $q>1$. This is not consistent with experimental observations showing that velocity fluctuations are independent of the viscosity in the inertial range. Thus, to this regard, the choice of an exponentiated Ornstein-Uhlenbeck process (51) as a model for pseudo-dissipation trajectories is not consistent with the standard phenomenology. We will see in the following that the multifractal model we consider indeed predicts simultaneously an appropriate Reynolds number dependence for the $\varphi$ moments (53) and an appropriate behaviour of the moments of $\tau_{\eta}^{2}\varphi_{\tau}$ (54) as the Reynolds number tends to infinity.

Appendix B A causal multiplicative chaos: Definition and statistical properties

B.1 Setup and notations

We consider the following linear stochastic differential equation

[TABLE]

with $W$ a Gaussian white noise and $\beta_{\tau_{\eta}}(t)$ a random function known in the literature as a (${\tau_{\eta}}$-regularized) fractional Gaussian noise of Hurst exponent $H=0$ (Mandelbrot & Van Ness, 1968):

[TABLE]

Equation 56 unique solution may be conveniently written as

[TABLE]

if we choose, formally, the initial condition $X_{\tau_{\eta}}(-\infty)=0$, taking a shortcut to define these Gaussian processes (58 and 57) starting from an infinitely ancient time. As we will see, this makes perfect sense and leads to finite variance processes (at a finite $\tau_{\eta}$). This enables us to work directly with a causal and stationary framework. From a numerical point of view, it is equivalent to start from any initial condition $X(t_{0})<+\infty$ or, similarly, truncate the stochastic integrals from (58) and (57) over $[t_{0},t]$ and propagate in time until reaching a stationary regime (see also the discussion provided in section 4). For convenience, we define the sub processes $X_{\tau_{\eta}}^{(1)}(t)$ and $X_{\tau_{\eta}}^{(2)}(t)$, such that $X_{\tau_{\eta}}(t)=X_{\tau_{\eta}}^{(1)}(t)+X_{\tau_{\eta}}^{(2)}(t)$, as

[TABLE]

which is known in the literature as a fractional Ornstein-Uhlenbeck (FOU) process (Cheridito et al., 2003) of Hurst exponent $H=0$, and

[TABLE]

a standard Ornstein-Uhlenbeck (OU) process.

A process similar to $X_{\tau_{\eta}}$, although different, has been considered heuristically by Schmitt (2003). In our formulation (56), the integral time scale $T$ comes into the dynamics through the damping term, which allows us to derive rigorously, in the stationary regime, its statistical properties. We preserve from the heuristics of Schmitt (2003) the correlated OU-process $X_{\tau_{\eta}}^{(2)}$ (60) which plays a crucial role in the determination of the statistics in the stationary regime.

B.2 Variance and covariance in the stationary regime

The Gaussian process $X_{\tau_{\eta}}(t)$ (58) reaches a stationary regime, with

[TABLE]

and

[TABLE]

Furthermore, in the asymptotic limit of infinite Reynolds number, i.e. ${\tau_{\eta}}\to 0$, we have the following logarithmic diverging behaviour

[TABLE]

In the asymptotic limit of infinite Reynolds number $\tau_{\eta}\to 0$, as opposed to the variance, the covariance remains bounded for $\tau>0$, and we note

[TABLE]

In the limit of vanishing time lag $\tau$, we have moreover the following logarithmic divergence of the correlation function

[TABLE]

B.3 Proofs

B.3.1 The Ornstein-Uhlenbeck process $X_{\tau_{\eta}}^{(2)}$

It is well known that

[TABLE]

as readily seen from (60). As for the covariance, we arrive at

[TABLE]

independently of $t$. Notice that by taking $\tau=0$ in (65) the variance itself may be written as

[TABLE]

We can see that both the covariance and the variance of $X_{\tau_{\eta}}^{(2)}$ diverge as the Reynolds number becomes infinite.

B.3.2 The fractional Ornstein-Uhlenbeck process $X_{\tau_{\eta}}^{(1)}$

Similarly, it is straightforward to verify that

[TABLE]

The process’ covariance is given by

[TABLE]

where $C_{\tau_{\eta}}(h)$ stands for the covariance of the fractional noise $\beta_{\tau}$ (57) and reads

[TABLE]

Therefore, it may be conveniently written as

[TABLE]

As for the variance, with the change of variables $t-s_{i}\rightarrow s_{i}$ for $i\in\{1,2\}$, we arrive at

[TABLE]

that can be worked out a bit to become

[TABLE]

Focusing on the integration of the dummy variable $h$, an integration by parts leads to

[TABLE]

B.3.3 Covariance of $X_{\tau_{\eta}}^{(1)}$ and $X_{\tau_{\eta}}^{(2)}$

Without loss of generality, consider for instance $\tau\geq 0$. One has

[TABLE]

where it is simple to give meaning to the following rule of calculation

[TABLE]

Here $1_{s_{1}\geq s_{2}}$ is the indicator function of the ensemble $\{(s_{1},s_{2})\in\mathbb{R}^{2},s_{1}\geq s_{2}\}$. Manipulating a bit, one gets

[TABLE]

where

[TABLE]

An integration by parts on the remaining integral finally leads to

[TABLE]

Similarly, recalling that $\tau\geq 0$, one may write

[TABLE]

B.3.4 Variance of $X_{\tau_{\eta}}$

Noticing that

[TABLE]

one obtains, using (66), (B.3.2) and (72),

[TABLE]

where

[TABLE]

Integrating (75) by parts results in

[TABLE]

with

[TABLE]

Now, inserting this into the variance (74), one finally derives

[TABLE]

which entails the proposition made in (61).

B.3.5 Asymptotics in the limit of infinite Reynolds numbers ${\tau_{\eta}}\to 0$

We rewrite the process’ variance rescaling the integration variable by ${\tau_{\eta}}$, obtaining

[TABLE]

with

[TABLE]

Notice that $r$ is a bounded function of its argument, in particular, $r(0)=1$ and, as $h\to\infty$, it behaves as

[TABLE]

To deduce the behaviour of the variance of the process $X_{\tau_{\eta}}$ (77) as $\tau_{\eta}\to 0$, we split the integral into two contributions, namely

[TABLE]

The first contribution is bounded with ${\tau_{\eta}}$ and tends to $\int_{0}^{1}r(h)dh$ when ${\tau_{\eta}}\to 0$. Nevertheless, the second one will diverge when ${\tau_{\eta}}\to 0$ and will thus dominate. To see how fast it diverges, we write it as

[TABLE]

Using the expansion (78) we see that $r(h)-\frac{1}{h}$ is integrable when $h\to\infty$, we can thus apply Lebesgue’s dominated convergence theorem and conclude that the first contribution is bounded with ${\tau_{\eta}}$ (and tends to $\int_{1}^{\infty}\left[r(h)-\frac{1}{h}\right]dh$ when ${\tau_{\eta}}\to 0$). Only the second contribution will diverge when ${\tau_{\eta}}\to 0$, so we can write

[TABLE]

To assess how this well behaved quantity diverges, we integrate by parts to observe that

[TABLE]

which concludes the proof of the proposition made in (62).

B.3.6 Asymptotics of the covariance of $X_{\tau_{\eta}}$

Starting from

[TABLE]

one has, using (65), (B.3.2), (71) and (B.3.3),

[TABLE]

Let’s denote the second term of the RHS of the former equation by $I_{\tau_{\eta}}(\tau)$. Notice first that one may integrate by parts the definition of $C_{\tau_{\eta}}$ (68) to stablish

[TABLE]

that may be used to obtain

[TABLE]

With the change of variables $h=s_{1}-s_{2}$ (while keeping the integration over $s_{1}$ and $v$), we perform a integration by parts over $h$, manipulate a bit, and arrive at

[TABLE]

The first term on the RHS of this equality cancels an identical term of opposite sign in (B.3.6). It is then not difficult to show that the two others terms on the RHS remain bounded, for a given $|\tau|>0$, when $\tau_{\eta}\to 0$. We have thus demonstrated that the covariance of $X_{\tau_{\eta}}$ remains bounded at infinite Reynolds number, and we designate

[TABLE]

The integral over $v$ may be carried on with the exact result provided in (76), which brings us to the simplified expression

[TABLE]

Performing an integration by parts over $h$, it is easy to justify that the second term in the RHS of the former equality remains bounded when $\tau\to 0$. As for the first term, it diverges when $\tau\to 0$ and will thus dominate the covariance of $X$ at small scales. For these reasons, we finally deduce the following equivalent of the covariance function at small scales:

[TABLE]

A similar integral has been encountered in section B.3.5, where it is shown that it diverges with $\tau$ according to $\log(T/|\tau|)$, which entails the proposition (64).

B.4 A causal multifractal process for pseudo-dissipation and its statistical properties

In the spirit of Schmitt (2003), we consider the process

[TABLE]

which contains the variance of $X_{\tau_{\eta}}$ in the stationary regime (77) and whose dynamics is given by

[TABLE]

The respective Lagrangian multiplicative chaos, which is causal and stationary, is readily obtained while exponentiating the Gaussian process $Y_{\tau_{\eta}}$:

[TABLE]

An application of Ito’s lemma leads to the stochastic dynamics of the pseudo-dissipation as seen by a Lagrangian particle along its trajectory, namely

[TABLE]

that will eventually reach a stationary regime for any bounded, non vanishing and positive initial condition.

From 82 we extract the statistical properties of $\varphi$. First, for positive integers $q$, its moments of order $q$ are given by

[TABLE]

As discussed in section B.3.5, the variance of $X_{\tau_{\eta}}$ diverges logarithmically with $\tau_{\eta}$. Assuming that the sub-leading term is constant, call it $\tilde{g}(0)$, we may write, when $\tau_{\eta}\to 0$,

[TABLE]

such that

[TABLE]

In a similar fashion, let us now compute the moments of the averaged pseudo-dissipation $\varphi_{\tau}(t)$ over a time interval $\tau$ (9). We get

[TABLE]

We have seen in section B.3.6 that the covariance of $X_{\tau_{\eta}}$ remains bounded for non vanishing time lag when $\tau_{\eta}\to 0$. Furthermore, in this limit of infinite Reynolds number, we have also seen that this asymptotical covariance diverges logarithmically at the origin. Assume now that the sub-leading terms remain bounded and call $\tilde{g}(t)$ such a bounded function. We can then write

[TABLE]

as it was found in the Eulerian context (13). It is then not difficult to show that $\tau_{\eta}^{2q}\mathbb{E}\left[\varphi_{\tau}^{q}\right]$ remains bounded when $\tau_{\eta}\to 0$ for $q<1+2/\mu^{l}$ to obtain the following behaviour at small scales:

[TABLE]

Note that the condition on $q$, that is $q<1+2/\mu^{l}$, ensures that the integrals in (B.4) exist.