Data-driven model for Lagrangian evolution of velocity gradients in incompressible turbulent flows

Abstract Velocity gradient tensor, $A_{ij}\equiv \partial u_i/\partial x_j$, in a turbulence flow field is modelled by separating the treatment of intermittent magnitude ($A = \sqrt {A_{ij}A_{ij}}$) from that of the more universal normalised velocity gradient tensor, $b_{ij} \equiv A_{ij}/A$. The boundedness and compactness of the $b_{ij}$-space along with its universal dynamics allow for the development of models that are reasonably insensitive to Reynolds number. The near-lognormality of the magnitude $A$ is then exploited to derive a model based on a modified Ornstein–Uhlenbeck process. These models are developed using data-driven strategies employing high-fidelity forced isotropic turbulence data sets. A posteriori model results agree well with direct numerical simulation data over a wide range of velocity-gradient features and Reynolds numbers.


Introduction
The velocity gradient (VG) evolution in incompressible turbulent flows depends upon four processes -inertial, pressure, viscous, and large-scale forcing if present.The inertial contribution arises from the momentum of the fluid and is local.The pressure contribution is given by the pressure Hessian, which is highly nonlocal (as pressure is governed by the elliptic partial differential equation -Poisson equation) and can be divided into an isotropic part which is local, and an anisotropic part which is nonlocal.The viscous forces on the evolution of velocity gradients at any point in the flow depend on its immediate neighborhood, thus also representing a nonlocal contribution.The large-scale forcing effect stems from the application of external forces that sustain the turbulence in the flow and is a function of the nonlocal forces in the immediate neighborhood.The interaction between all these effects leads to the complicated behavior of small-scale intermittency and multifractality in turbulence (Yakhot & Donzis 2017;Sreenivasan & Antonia 2017;Meneveau & Sreenivasan 1991).Despite these features, turbulence small scales exhibit a certain degree of universality (Kolmogorov 1941;Sreenivasan 1998;Schumacher et al. 2014).
Due to its theoretical significance and its practical utility in a variety of applications, there have been numerous attempts at modeling the Lagrangian evolution of the velocity gradient tensor, which is a difficult pursuit due to the inherent challenges of modeling nonlocality and intermittency.As discussed above, in velocity-gradient evolution, the local inertial and isotropic pressure Hessian contributions are closed, while the remaining non-local processes require closure.The earliest attempts (Cantwell 1992;Girimaji & Speziale 1995;Martın et al. 1998) at modeling VG dynamics neglected the non-local terms and modeled only the closed restricted Euler (RE) equations (Vieillefosse 1982;Cantwell 1992) which led to finite-time singularity.Beginning with the work of Girimaji & Pope (1990), a series of stochastic velocity gradient models followed over the years that used diverse closure techniques for modeling the effects of the non-local pressure and viscous contributions leading to statistically stationary solutions (Chertkov et al. 1999;Jeong & Girimaji 2003;Chevillard & Meneveau 2006;Chevillard et al. 2008;Wilczek & Meneveau 2014).Some of the recent modeling efforts, such as the recent deformation of Gaussian field (RDGF) model proposed by Johnson & Meneveau (2016), the multifractal process-based stochastic model by Pereira et al. (2018), the model based on temporal correlation of strain and rotation rate by Leppin & Wilczek (2020), and the data-driven VG model employing tensor-basis neural network for pressure term's closure (Tian et al. 2021), have shown improvements over previous models.However, all these models are unable to simultaneously capture both the statistical properties of intermittency and the geometric features of small-scale turbulence with high accuracy.There is a need for a robust VG model that accurately predicts both these essential features of small-scale dynamics and is generalizable to a wide variety of turbulent flows.
Modeling the nonlinear dynamics of the VG tensor (   ≡   /  ) is challenging due to its multifractal and intermittent nature.To overcome this complexity, in this work we adopt a different approach: we isolate the universal features of turbulence from intermittency.Kolmogorov's refined similarity hypothesis (Kolmogorov 1962) suggests that considering intermittency, using the local average of dissipation rate better encapsulates the universality of turbulence than its global average.Extending the same principle, we use the absolute local pseudo-dissipation rate or VG magnitude ( = √︁       ) to normalize the VG tensor as follows (Girimaji & Speziale 1995): where  ≡ || ||  = √︁     . (1.1) The normalized velocity gradient tensor,    , is a mathematically bounded tensor and is statistically nearly universal across different types of turbulent flows at different Reynolds numbers (Das & Girimaji 2019, 2022).The geometric shape features of the turbulence small scales are encoded in the    tensor, while its scale and intermittency lie in the VG magnitude  (Das & Girimaji 2020).The idea here is to develop separate models for    and  tailored for capturing the dynamical behavior of each uniquely, resulting in an overall improved prediction of    evolution.Within the    framework, the effects of different turbulence processes including the nonlocal pressure, viscous, and forcing, are nearly universal, wellbehaved, and more amenable to modeling (Das & Girimaji 2020, 2022) than the nonlocal terms in    space modeled in previous studies.Taking advantage of these properties, we model the unclosed pressure and viscous contributions in the mathematically bounded statespace of    .On the other hand, the intermittency is captured in a separate model of magnitude .
To model the nonlocal terms, we employ a simple data-driven approach based on our physical understanding of the small-scale dynamics.For this, we utilize highly resolved direct numerical simulations (DNS) data and rely on a lookup table approach in a four-dimensional compact space.It provides a more exact representation of the flow physics as compared to other data-driven modeling methods such as neural networks, due to the compactness of the    space.This results in a generalizable closure of the non-local pressure and viscous processes at high enough Reynolds numbers within the    framework.
Modeling the VG magnitude  does not require any additional closures.We model the evolution of magnitude  (pseudodissipation rate ∼  2 ) within the framework of the Ornstein-Uhlenbeck (OU) process (Uhlenbeck & Ornstein 1930), due to its near lognormal probability distribution and exponential decay of its auto-correlation (Kolmogorov 1962;Oboukhov 1962;Yeung & Pope 1989).Although multifractal formalism (Mandelbrot 1974) suggests that pseudodissipation rate is not precisely lognormal, studies have shown support for the lognormal framework of modeling the temporal dynamics of VG magnitude (Pope & Chen 1990;Girimaji & Pope 1990;Huang & Schmitt 2014).Therefore, we model the VG magnitude as a Reynolds number-dependent modified lognormal process.We further incorporate DNS data-based physical modifications within the OU process, with the expectation of capturing the intermittent nature of small-scale turbulence more accurately than a simple lognormal process.
Overall, this work presents a data-driven Lagrangian model to accurately reproduce all the essential characteristics of VG dynamics in turbulent flows for a broad range of Reynolds numbers with minimal computational effort.The novelty of our model lies in two primary features of our approach: (1) we model the geometric aspect of VG dynamics separately from its intermittent magnitude, and (2) within the compact space of VG geometry, we use a cleverly articulated, simple but exact lookup table approach for the data-driven closure of the non-local pressure and viscous processes.The remaining sections of the paper are arranged as follows.In section 2 we discuss the properties and present the governing differential equations for the normalized VG tensor and VG magnitude in an incompressible turbulent flow.The entire modeling methodology is described in section 3, including the philosophy of the modeling approach and its generalizability, formulation of the model equations and closures, and a complete model summary.The numerical procedure of the simulations performed using the model is outlined in section 4. Finally the results of the model are compared with that of DNS and previous models in section 5 and the conclusions are presented in section 6.

Governing equations
The Navier-Stokes and continuity equations for velocity fluctuations,   , in an incompressible turbulent flow can be written as where,  is the pressure fluctuation,  is the kinematic viscosity, and   represents the largescale forcing.Pressure and viscous effects are the key nonlocal processes in turbulence.Forcing, which causes the energy production at large scales to compensate for the viscous dissipation at small scales, can be expressed in the following general form for most commonly encountered flows with a mean flow: where,   =   +   is the total velocity and indicates ensemble averaging or spatial averaging in homogeneous directions.The effective forcing   varies from one turbulent flow to another depending on the mean flow field as well as the inhomogeneity and anisotropic nature of the flow geometry (Rogallo 1981).In homogeneous isotropic turbulence with no mean flow, forcing simply entails injecting energy at the lowest wavenumbers (Eswaran & Pope 1988;Donzis & Yeung 2010).
From equation (2.1), the governing equation for the velocity gradient tensor can be derived as: where Here, / = / +   /  is the material or substantial derivative.The first two terms on the right-hand-side (RHS) of equation ( 2.3a) represent the non-linear effects, but they are local in space.The tensor    is the anisotropic pressure Hessian tensor,    is the viscous Laplacian tensor, and    is the anisotropic forcing tensor.   and    tensors represent the non-local effects in velocity gradient dynamics and    depends on the nature of forcing.

Normalized VG tensor
The evolution equation for    in the flow frame of reference, derived from equation (2.3a), is where  =  is the time increment normalized by local VG magnitude, and the timescale  is referred to as the local timescale.Here, .5)are the normalized anisotropic pressure Hessian, viscous Laplacian, and anisotropic forcing tensors, respectively.In the    equation (2.4), the first three terms on the RHS are closed and represent the nonlinear () -inertial and isotropic pressure Hessian -effects.The next three terms constitute the non-local pressure (), viscous (), and forcing () effects on    evolution that require closure.Consideration of the    evolution in local timescale is not only consistent with Kolmogorov's refined similarity hypothesis but also leads to the added practical advantage that all the terms on the RHS of equation 2.4 are normalized tensors.These normalized pressure Hessian and viscous Laplacian tensors are not necessarily bounded, but are well-behaved in the phase plane of    invariants and are considerably more amenable to modeling than the unnormalized tensors in the    equation (Das & Girimaji 2022).They further exhibit a nearly universal behavior across turbulent flows of different Reynolds numbers (Das & Girimaji 2020, 2022).Therefore, the Lagrangian evolution of    can be modeled in the local timescale  without any explicit dependence on the VG magnitude.The magnitude dependence comes in only when determining the    evolution in real time.
In order to most effectively employ data-driven techniques, we will now establish the minimum number of free parameters required to completely describe the    tensor.A detailed derivation can be found in Das & Girimaji (2020); only the main outcomes are summarized below.Without any loss of generality, we can express    in the principal Focus on Fluids articles must not exceed this page length (eigen) reference frame of normalized strain-rate tensor,    , as follows where  1  2  3 (2.6)Here, ( ) represents tensors in the principal reference frame of    , and   are the eigenvalues of    in decreasing order, such that  1 (> 0) is the most expansive strain-rate,  3 (< 0) is the most compressive strain-rate, and the intermediate strain-rate  2 can be positive, negative or zero, in incompressible flows.Further,   are the components of the normalized vorticity vector ( ì ) along the strain-rate eigendirections.Since the signs of the strain-rate eigenvectors are not uniquely determined by the eigendecomposition, we consider the eigendirections that provide all vorticity components to be of the same sign (either all positive or all negative).
Applying the constraints of incompressibility ( b = 0) and normalization ( b  b  = 1), the b  state-space can be reduced to a four-dimensional space of only four independent variables (shape-parameters) -, ,  2 ,  2 .Here,  and  are the second and third invariants of the tensor, respectively: All the remaining elements of b  can be determined uniquely once these four variables are known, as shown below: ) ) .8c)This shows that , ,  2 , and  2 completely define the tensor b  and thence the geometricshape of the local flow streamlines.These four variables are also mathematically bounded as follows: , (2.9) . (2.10)

VG magnitude
The velocity gradient magnitude or pseudodissipation rate has been shown to have a nearly lognormal distribution (Kolmogorov 1962;Oboukhov 1962;Yeung & Pope 1989;Monin & Yaglom 2013).For this reason, we consider the dynamics of the logarithm of VG magnitude: (2.11) which is expected to exhibit a near-normal distribution in a turbulent flow field.We introduce the scalar variable -standardized VG magnitude: which exhibits a nearly standard normal distribution, N (0, 1), in a variety of turbulent flows (Das & Girimaji 2022).The evolution equation for  * , derived from equation (2.3a) is given by: where  * =  .Here,  * is referred to as the global timescale and it represents the timescale normalized by the global mean of VG magnitude.This normalization is in essence similar to normalization by the Kolmogorov timescale (  ∼ 1/  2 1/2 ) and is found to be more appropriate for examining VG magnitude than the local timescale used for    .The four terms on the RHS of the above equation represent the nonlinear, pressure, viscous, and forcing effects, respectively, on the VG magnitude evolution.

Model Formulation
Instead of considering the evolution of    directly, this work isolates the modeling of the normalized velocity gradient tensor (   ) from that of the magnitude () as shown in the schematic of figure 1.In this section, the universal features of the dynamics that need to be captured by a velocity-gradient model and that can be used to our advantage in the modeling approach, are first discussed.This is followed by the main modeling strategies and a detailed description of the complete model.Finally, all the model equations and parameters are summarized.

Generalizability of modeling VG dynamics
As outlined in figure 1, the large scales of motion in a turbulent flow depend upon the flow geometry and driving mechanism of the flow.It is therefore difficult to develop generalizable models for the large scales that will apply to different turbulent flows.Models of smallscale dynamics are likely to be more generalizable in comparison since the small scales in In this model, we separate    into normalized velocity gradient tensor (   ) and velocity gradient magnitude (), such that the tensor    is nearly universal across different turbulent flows while the scalar  reflects all the Reynolds number dependence.The universality is evident in the PDF and higher order moments of    (Das & Girimaji 2019) as well as in the mean evolution of    invariants (Das & Girimaji 2020, 2022), that are insensitive to the variation of Taylor Reynolds number (  ) across different turbulent flows.Therefore,    evolution modeled using DNS data of only one turbulent flow at a given Reynolds number may be considered universal (up to a modeling approximation) and can be applied to reproduce the    -dynamics of different turbulent flows at different Reynolds numbers.
The magnitude , on the other hand, exhibits a strong dependence on the Reynolds number of the flow.The mean and variance of its logarithm ( = ln ), plotted in figure 2, clearly increase with increasing   .Preliminary results suggest that  and   follow approximate scaling laws with   , as indicated in the figures.The scaling law obtained for  2   is in close agreement with the   -scaling of the logarithm of pseudodissipation rate reported by Yeung & Pope (1989).But further advanced simulations and analyses are required to develop universal scaling laws for  and  2  .In fact, these two quantities are the input parameters of our model for VG magnitude (section 3.4), representing the characteristic Reynolds number dependence of velocity gradients.
As summarized above, the advantage of this modeling framework is that the ninecomponents tensorial variable    is nearly universal, and one can develop a potentially generalizable    -model applicable to different types of turbulent flows at wide-ranging Reynolds numbers.Only the scalar -model is Reynolds number dependent, which can be represented by scaling laws that are likely generalizable across different types of turbulent flows.

Modeling strategy
The modeling approach of this work constitutes the following: (i)    -model: The    dynamics in a local timescale (equation 2.4) is a function of    and other normalized non-local tensors, and it does not explicitly depend on magnitude .Therefore, we formulate a stochastic model for the Lagrangian evolution of    in the local timescale ( ) without any explicit dependence on  * .
(ii) Closure of nonlocal processes: As inferred from DNS data in our previous analysis (Das & Girimaji 2020, 2022), the conditional statistics of the normalized non-local tensors can be reasonably approximated as exclusive functions of    .Thus, we develop DNS datadriven closure models (generalizable at high   ) for capturing the conditional mean nonlocal effects of normalized pressure and viscous processes within the four-dimensional bounded state-space of b  .The fluctuations of these nonlocal effects as well as the effect of large-scale forcing are modeled in the stochastic diffusion term using moment constraints.
(iii)  * -model: We model the evolution of VG magnitude in global timescale ( * ) within the framework of Ornstein-Uhlenbeck (OU) process (Pope & Chen 1990) in three different ways.The first model is a simple OU model for  * decoupled from    dynamics.The second and third models are modified OU models with    -dependence incorporated into the  * evolution using a DNS data-based diffusion process.
(iv) Timescale: In addition, an ordinary differential equation provides the relation between the local and the global timescales.Finally, the    and  * models are combined to form an integrated system of model equations representing the Lagrangian evolution of    in global time.

Model for normalized VG tensor
The Lagrangian dynamics of normalized velocity gradient tensor,    , is modeled here as a diffusion process (Karlin & Taylor 1981).A diffusion process is a continuous-time Markov process and is represented by a stochastic differential equation (SDE).The SDE for    commonly used in previously developed models (Girimaji & Pope 1990;Chevillard & Meneveau 2006;Chevillard et al. 2008;Johnson & Meneveau 2016) is of the form where    is a tensor-valued isotropic Wiener process such that The    tensor represents the drift coefficient tensor and     constitutes the diffusion coefficient tensor of the model.Taking the trace of equation ( 3.1), one can show that   =   = 0 satisfies the incompressibility condition   = 0. Starting from the above equation and using the properties of an Itô process (Kloeden et al. 1992), one can derive the following SDE for    in local timescale  (see appendix B for derivation): where, and the Wiener process satisfies It is important to note that all the drift and diffusion coefficient tensors of this system of SDEs are dimensionless.The drift tensor of the    equation has two parts due to the normalization: (i)    is obtained from the drift tensor of the    equation,    , and (ii)    is obtained from the diffusion tensor of the    equation,     .The diffusion tensor of the    equation,     , is also obtained from     .The tensor    relates the drift and diffusion processes in the dynamics such that despite the random stochastic forcing term,    remains mathematically bounded.All the coefficient tensors are modeled in the specific functional forms given above.
It can be proved that any system of SDEs for    , that complies with the above forms of drift and diffusion terms, clearly satisfies the incompressibility constraint: Equation ( 3.3) further satisfies the mathematical constraint of normalization: which ensures that the Frobenius norm of the tensor  is equal to unity at all times.The proofs are presented in appendix C and D. Equation (3.3) leads to a Fokker Planck equation (Pope 1985) for the joint PDF F() of the tensor    : Now, the exact differential equation for the joint PDF of    , F(), in a turbulent flow can be derived from the    governing equation (2.4) as The drift and diffusion coefficient tensors need to be modeled in a way that F() ≈ F().In data-driven modeling, F() and its moments are taken from high-fidelity DNS data.However, requiring the PDFs to be identical is very challenging.In this work, we constrain the equations of    -moments up to third order to obtain the parameters of diffusion coefficient tensor along the lines of Girimaji & Pope (1990).This modeled diffusion process is, therefore, consistent up to order three, although the numerical results of the model show reasonable agreement of much higher-order moments.

Drift coefficient tensor
Comparing the terms of equations (3.8) and (3.9), the inertial and isotropic pressure Hessian terms are exact, and considering that the role of     is to model the large-scale forcing effect and that of    is to maintain the unit Frobenius norm of    , the drift coefficient tensor    takes the form: 10) The term    represents the dynamics of the inertial and isotropic pressure Hessian contributions as well as the conditional mean of the anisotropic pressure Hessian and viscous contributions.
The conditional mean normalized anisotropic pressure Hessian and viscous Laplacian tensors, ℎ   | and    | , require closure modeling.As discussed in section 2.1, the tensor b in the principal reference frame of the strain-rate tensor can be expressed as a function of only four bounded variables.Therefore, in order to take advantage of this four-dimensional bounded state-space of b, the conditional averaging of the normalized pressure Hessian and viscous Laplacian tensors is performed in the    principal reference frame.For a rotation tensor, , with columns constituted by the right eigenvectors of , the required tensors in the principal reference frame are b Then, the conditional mean pressure Hessian and viscous Laplacian tensors in the flow reference frame can be recovered as follows: since  is a function of .Therefore, the conditional mean pressure Hessian and viscous Laplacian tensors in the flow reference frame can be obtained if the conditional mean pressure Hessian and viscous Laplacian tensors in the principal frame are known and the local strain-rate eigenvectors are known.
As shown in section 2.1, in the principal reference frame, the tensor b  can be uniquely defined by a set of only four bounded variables -, ,  2 and  2 .Therefore, the conditional mean pressure Hessian and viscous Laplacian tensors in the principal frame can be modeled as a function of only four bounded variables, as follows: The goal is to develop a data-driven model for the above tensors in terms of a four-dimensional input.Das & Girimaji (2020, 2022).Therefore, the simpler and more accurate data-driven approach of direct tabulation based on DNS data is used in this work.
The approach can be summarized as follows: (i) The finite space of , ,  2 and  2 , as given in equation ( 2.10), is discretized into (60, 60, 30, 30) uniform bins.This discretization strikes the appropriate balance between sampling accuracy in the bins and the desired details of nonlocal flow physics to be captured.Other discretizations are tested to show convergence to this combination for the most accurate results.
(ii) The conditional expectations of the tensors, h  |, ,  2 ,  2 and τ  |, ,  2 ,  2 , are computed in this discrete phase-space, using DNS data set of forced isotropic turbulence (see appendix G for details of the dataset).Note that only one lookup table is required to model the mean nonlocal dynamics for turbulent flows of different   .Further, this datadriven closure is rotationally invariant since the input and output variables do not depend on the flow reference frame.
(iii) This lookup table can then be accessed by an inexpensive array-indexing operation, to determine the conditional mean pressure and viscous dynamics in the principal frame for a given (, ,  2 ,  2 ) at any point in the flow field or following a fluid particle.This is then transformed to the flow reference frame (equation 3.12) based on the local eigendirections of strain-rate tensor, to be used in    for computations.This completes the modeling of the mean drift tensor    as a function of the local    and it is straightforward to show that our data-driven model for    is Galilean invariant.The proof for the same is presented in appendix E.
Our previous work has shown the universality of    statistics and associated nonlocal processes across different types of turbulent flows (Das & Girimaji 2022).Therefore, in this work, we restrict ourselves to forced isotropic turbulence.We use DNS data of isotropic turbulence at different Reynolds numbers (appendix G) to illustrate the universality of different modeling components.Although the data of different   are considered, since the nonlocal terms requiring closure are insensitive to   variation, only the   = 427 case is used for final comparison.

Diffusion coefficient tensor
As discussed at the beginning of this section, the interrelationship between the tensors     and    is important in guaranteeing that the mathematical and physical constraints of    are satisfied.This relationship holds if we use their functional forms as given in equation (3.4), in terms of     from the    SDE.For this, we assume a general isotropic form of the four-dimensional tensor,     , following previous models (Girimaji & Pope 1990;Chevillard et al. 2008;Johnson & Meneveau 2016): where  1 ,  2 ,  3 are constant dimensionless diffusion coefficients of the model.Only two of these three coefficients are independent subject to the incompressibility condition: From equations (3.4), (3.14) and (3.15), the diffusion coefficient tensor of the    equation is and (3.17) To determine the constant coefficients,  2 and  3 , we use a moments constraint method similar to Girimaji & Pope (1990).In this method, the equations of second and third order moments of    are constrained to the values obtained from DNS.First, the SDEs for the second () and third () invariants are derived from the    -SDE (3.3) using Itô's lemma (appendix D): Taking mean of the above equations and substituting the expressions for     and    from equations (3.16) and (3.17), yields the following differential equations of the moments - and  : To model a statistically stationary solution of turbulence, the rate-of-change of moments must be driven to zero while ensuring that the moment values converge to that of DNS.For this, we equate the RHS to negative of the error term: where ,  are the global mean of ,  obtained from DNS data.Here,  represents the rate of convergence of these moments and is set to unity.An a priori simulation of the    model equations is run in the normalized timescale  , with an ensemble of 40000 particles.At each time step, the above system of nonlinear equations is solved using Newton's method to determine the values of the coefficients  2 ,  3 .In this a priori run, the coefficients converge to the following values:  2 = 0.0099 ,  3 = −0.064(3.22) as the model's moments,  and  , converge very close to the DNS values of  and .These optimized diffusion coefficient values are used in the stochastic model for    and are insensitive to the Reynolds number.

Model for VG magnitude
The Lagrangian evolution of the scalar  * (equation 2.12) is modeled using a modified lognormal approach.The magnitude  has a nearly lognormal probability distribution and exponential decay of autocorrelation in time (Kolmogorov 1962;Oboukhov 1962;Yeung & Pope 1989).The exponentiated Ornstein-Uhlenbeck (OU) process is a statistically stationary process that satisfies both these properties (Uhlenbeck & Ornstein 1930;Pope & Chen 1990) and is therefore ideal for modeling  * .While it has been pointed out that pseudodissipation rate ( 2 ) cannot be precisely lognormal in the context of multifractal formalism (Mandelbrot 1974;Meneveau & Sreenivasan 1991), the OU process models the overall dynamics of  quite accurately (Pope & Chen 1990; Girimaji & Pope 1990).In fact, a recent analysis of Lagrangian trajectories in high Reynolds number turbulence (Huang & Schmitt 2014) has shown evidence that the autocorrelation function of  is consistent with both the exponential decay prescribed by the OU process as well as the logarithmic decay suggested by the multifractal framework, and the two are nearly indistinguishable at such high Reynolds numbers (Pereira et al. 2018).Since the focus of this work is to accurately reproduce the overall Lagrangian dynamics of the velocity gradients in turbulence, we model the velocity gradient magnitude as a Reynolds number-dependent modified lognormal process, without explicitly accounting for multifractal behavior.
The OU process is a stationary continuous Gaussian Markov process that is often used in modeling systems of finance, mathematics, and physical and biological sciences (Pope & Chen 1990; Klebaner 2012).It further shows the property of mean-reversion.The SDE for a general OU process  * evolving in time  * is given by (Girimaji & Pope 1990) where, ,  > 0 are parameters of the model,  * =   is the non-dimensional global timescale and  * is the increment of a Wiener process or a Gaussian random variable with zero mean and variance  * .The parameter  represents the rate of mean-reversion, and without loss of generality it is set to unity since the model propagates in timescale  * which is already normalized.The expected value  * = 0, by construction, in DNS data.Therefore, the general form of  * -SDE used in this work is The diffusion coefficient, , is modeled in three different ways as described below.

Model 1 -simple OU process
Here, we disregard the dependence of  * on    and consider the simple OU dynamics that satisfies the global mean and global variance of  * .In this case, the diffusion coefficient  is taken to be a constant value, which is calculated as follows.The equation for the global mean is obtained from equation (3.24), Thus, the model maintains a stationary mean value of  * once the solution is driven to the zero mean value by the mean-reversion property.Next, the equation for the global variance is obtained from equation (3.24) using Itô's product rule, For a statistically stationary solution, we must have Therefore, the final form of the  * -SDE for Model 1 is given by Here, the value of the variance  * 2 is obtained from DNS data and it can be a function of   .(3.30) The conditional variance remains statistically stationary only if Here, the value of  0 is calculated based on the solutions of models 1 and 2. Model 3 also depends upon    through its invariants (, ).

Model summary
The resulting model for the Lagrangian evolution of the complete velocity gradient tensor in a turbulent flow is described by a system of stochastic differential equations for the normalized velocity gradient tensor,    , and a separate stochastic differential equation for the standardized VG magnitude,  * .
The final system of equations for    in local timescale ( ): = (   +    ) +       (3.34) Here, the diffusion coefficient values are The conditional mean normalized pressure Hessian and viscous Laplacian tensors are obtained from the data-driven closure in the strain-rate () eigen reference frame as a function of the current (, ,  2 ,  2 ), followed by a rotation to the flow reference frame using the local eigenvectors of : • Model 1 - where,  0 = 0.103.The diffusion coefficients of models 2 and 3 are obtained from the tabulated conditional variance of  * invariant with   (figure 3).The VG magnitude and the VG tensor are then given by  =  The numerical solution of the model equations involves numerically propagating the velocity gradient tensor, in terms of the variables    and  * , of an ensemble of 40000 particles.As the initial conditions for the simulations, the particles are picked at random from a randomly generated incompressible isotropic velocity field.The trajectories are advanced for a total time period of approximately 1200  following these steps at each update: (i) The    SDEs (equation 3.34) are numerically propagated in the normalized local timescale  , using a second-order weak predictor-corrector scheme (see appendix D) with a constant time increment  .At each step, the conditional mean nonlocal pressure and viscous contributions are calculated based on the current (, ,  2 ,  2 ) values, using the (60, 60, 30, 30) sized lookup-table.
(ii) The  * SDE (equation 3.38, 3.39 or 3.40) is advanced using the second-order weak predictor-corrector numerical scheme (appendix F) in the global timescale  * , using a firstorder approximation of the increment  * =   / for a fixed value of  .
(iii) The global timescale,  * , is obtained for every particle at each  by numerically solving the ordinary differential equation 3.43 using the implicit second-order Trapezium rule method.The model solution propagates all particles at a uniform local time increment of  = 0.01, but the global timescale  * varies from one particle to the other depending upon its current velocity gradient magnitude.A particle with a smaller magnitude requires fewer steps in  to reach a certain  * , than a particle with a larger magnitude.The VG magnitude  * evolves in global time  * , which approximately scales with Kolmogorov timescale.On the other hand,    evolves in local timescale  which varies depending on the local value of .Issues may arise when  <  , i.e. when    evolves faster than  * , and appropriate measures should be taken to ensure that the  is suitable to propagate the    equations.However, particles with such low  values do not contribute significantly toward the overall velocity gradient statistics.Convergence of the model's results for  = 0.05, 0.01 and 0.002 suggest that  = 0.01 is sufficient here for accurate statistical modeling.
The incompressibility and normalization constraints are automatically upheld by the model, but are only valid up to the order of numerical error.Therefore, to avoid the accretion of numerical errors over large periods of time, hard constraints of   = 0 and ||||  = 1 are enforced after every update.The computation time is approximately 1.5-2 hours on a single processor for the model's simulations to achieve statistically stationary solutions.The results of the model's simulations for the three different  * -models are illustrated in the next section as model 1 if equation (3.38) is used, model 2 if equation (3.39) is used, and model 3 if equation (3.40) is used, each along with the    equation (3.34).The convergence of all the major results have been tested for these models by performing the simulations with 40000 and 100000 particles.

Results and comparison with DNS data
This section presents a statistical analysis of the solutions of the three models and a comparison with the statistics of the corresponding DNS data and some previous models.First, the statistics of  * are illustrated, followed by the statistics of    .Finally the complete velocity gradient tensor    -statistics are shown.The time evolution of the model's statistics are illustrated as a function of the global normalized time  * .The time-converged statistical results are plotted by averaging over multiple time realizations of the model's solution, separated by at least 5  , well after statistical stationarity has been achieved.The time axis is in logscale.

VG magnitude
First, the time evolution of the mean and standard deviation of  * are plotted for all the three  * -model equations in figure 4. Note that the time axes are plotted in logscale to display the transients clearly.The numerical simulation starts from an initially random field, which is inconsistent with the DNS values of  and   , and therefore the initial values of  * and  *  are different from zero and unity, respectively.Over time, the model's solution evolves toward the DNS value achieving statistically stationary state at about  ≈ 7  ( * ≈ 6), where  is the real time.As expected, the global mean of  * is captured equally well by all three models due to the mean-reverting property of the OU process.The global standard deviation of  * is reproduced accurately by both models 1 and 3.However, it is worse in model 2 compared to the simple OU model (model 1).This indicates that imposing  * to satisfy the variance conditioned on local streamline geometry (, ) in model 2 does not necessarily guarantee that the global variance of  * is automatically satisfied.This justifies the need for the third model to satisfy both the global and conditional standard deviation values.
The converged probability density function (PDF) of the standardized VG magnitude,  * , for all the models and DNS data are plotted in figure 5.It is clear that models 1 and 3 are able to reproduce the  * PDF very well, while model 2 shows deviation from the desired DNS result.The plot in the log-linear scale confirms that the converged PDFs of models 1 and 3 agree well with that of DNS even near the extreme tails of the PDFs.Next, the converged PDFs of the VG magnitude (/  ) in each of the three models and DNS are plotted in figure 6.All models capture the peak of the PDF reasonably well but model 2 deviates at higher values of magnitude while models 1 and 3 perform better.It further shows that model 3 is able to reproduce the tails of the PDF slightly better than model 1.

Normalized VG tensor
The conditional mean trajectories (CMTs) in the phase plane of normalized velocity gradient invariants (, ) are examined as an a priori test of the data-driven closure used to capture the conditional mean nonlocal effects of pressure and viscosity (section 3.3) on the    -dynamics.The - CMTs are obtained by integrating the vector field of conditional mean velocity (ṽ) in the - plane: due to the inertial, pressure and viscous processes in the turbulent flow.Note that the effect of the large-scale forcing is not included here because it is not accounted for in the data-driven closure but rather in the stochastic forcing (diffusion) term in the    -SDE, which can not be tested a priori.The - CMTs obtained directly from DNS data are plotted in figure 7(a).
As discussed in Das & Girimaji (2022), trajectories closer to the origin converge toward the attractor near the origin (represents pure-shear streamlines) while trajectories that are outside the separatrix loop are attracted toward the bottom line attractor (represents purestrain streamlines).This behavior is almost exactly replicated by the - CMTs computed using the model's data-driven closure for the conditional mean pressure Hessian and viscous Laplacian tensors (section 3.3.1) in the equation 5.1 (figure 7b).The close resemblance between the two is somewhat expected given the very nature of the lookup table approach for closure.Now, we compare the a posteriori results of the normalized VG tensor of the model with that of DNS.First, we study the moments of the second () and third () invariants of the tensor, which are important quantities as they determine the geometric shape of the local flow streamlines.The evolution of up to fourth-order moments of  and  are plotted for each model in figure 8.It is first evident that all three models with the same    -SDE but different  * -SDEs produce nearly identical ,  moment values.Thus, it appears that the variation in the  * model does not have a discernible impact on the    statistics of the models.Starting from a randomly generated set of initial conditions, the    -SDE drives the solution toward convergence to a statistically stationary state at  ≈ 72  ( * ≈ 60).Therefore, the    model takes approximately 10 times as long as the  * model to reach stationarity.Up to at least fourth-order converged moments of both  and  are reasonably close to the DNS values.Further, the time evolution of moments of correlation between  and , i.e.  , and  2  2 , are plotted in figure 9.These moments show a slightly larger deviation from the DNS values as compared to all the other moments.
The evolution of - joint PDF is now investigated for the propagation of model 3.The solutions of the other two models show similar trends and are, therefore, not presented separately.The - joint PDF is plotted with ensembles of only 40000 particles at different The alignment of vorticity with strain-rate eigendirections is another key feature of smallscale turbulence.From equation (2.6) in the eigen reference frame of the strain-rate tensor, one can show that the cosine of the angles of alignment between the vorticity vector and the three strain-rate eigenvectors are given by, The angles  1 ,  2 ,  3 represent the angles of alignment of vorticity with the most expansive, intermediate, and most compressive strain-rate eigenvectors, respectively.In figure 12, the ).Yet, the    statistics of the model appear to be unchanged with the variation of  * model equation.

VG tensor
After examining the statistical results of    and  individually, we now test the models' performance in capturing the overall velocity gradient (   ) statistics.First, the PDFs of the longitudinal ( 11 ) and transverse ( 12 ) velocity gradients, normalized by their global root-mean-square values, are examined for all the three models in figure 13.For comparison, we have also plotted the corresponding PDFs obtained from the DNS data and two of the recent velocity gradient models that have shown improved results compared to the other models in the literature -(i) recent deformation of Gaussian field (RDGF) model by Johnson & Meneveau (2016) and (ii) physics-informed machine learning (PIML) model by Tian et al. (2021).Our models are able to reproduce the skewed  11 -PDF and the symmetric  12 -PDF, as observed in DNS, with reasonable accuracy.They show significant improvement  in capturing the PDFs of both  11 and  12 compared to both RDGF and PIML models.
On closer observation, it is evident that while the PDFs are captured nearly perfectly in the densely populated part by all three models, there are smaller differences near the tails of the PDFs.Models 1 and 3 predict a slightly heavier-tailed distribution of  11 than DNS, while model 2 produces a more accurate  11 -PDF.On the other hand, model 3 appears to capture the  12 -PDF tails slightly more accurately than the other two.
In order to determine the finer differences in the PDFs, we examine the higher order moments of the velocity gradient magnitude, , as well as the velocity gradient components,  11 and  12 .These moment values for all three models of this work, DNS data, and the RDGF model of Johnson & Meneveau (2016) are presented in table 2. The moment value produced by a model that is closest to the DNS is marked in bold type font.It is evident that the moments of magnitude  are best captured by model 3.The skewness, kurtosis, and 6  ℎ order moment of the longitudinal component  11 are reproduced best in model 2, although model 3 is not far behind and is slightly better than model 1.The skewness of the transverse component  12 is correctly captured as zero by all the models, maintaining a symmetrical probability distribution in each case.The kurtosis and 6  ℎ order moment of  12 are also captured most closely by model 3. Overall, model 3 provides the most accurate representation of the probability distributions and moments of the velocity gradient tensor.
Finally, the PDFs of the dissipation rate (      ), enstrophy (      ), and pseudodissipation rate (       ), normalized by their global means, are computed from the converged stationary state solution of all the three models and plotted in figure 14.The PDFs obtained from the DNS data and those available from the RDGF model (Johnson & Meneveau 2016) are also illustrated for comparison, along with the PDFs for the initial Gaussian field used in our model's simulations.It is interesting to note that the model is able to start from this Gaussian field and develop a turbulent flow field solution closely resembling that of DNS with the characteristic PDF-tails at extreme values.It is clear that all three models reproduce the heavy-tailed probability distributions of both dissipation and enstrophy more accurately than the RDGF model.Model 2 provides the most accurate representation of the dissipation PDF while models 1 and 3 over-predict the probability of occurrence of large dissipation rates near the tails of the PDF.Enstrophy, which is more intermittent in nature than dissipation rate (Yeung et al. 2018;Buaria et al. 2019), is captured best by model 3. Models 1 and 2 under-predict the probability density of enstrophy near the extreme tails.Taking the sum of the dissipation rate and enstrophy results in the pseudodissipation rate, which is reproduced quite accurately by model 3, even near the extreme tails.Overall, the results of model 3 constitute the closest representation of the velocity gradient statistics in turbulent flows.

Conclusion
A stochastic model for the Lagrangian evolution of velocity gradient (VG) tensor in an incompressible turbulent flow is presented.The bounded and well-behaved dynamics of the normalized velocity gradient tensor (   ) is modeled separately from the intermittent velocity gradient magnitude ().The main nonlocal flow physics of pressure and viscous processes are well-behaved and amenable to modeling in the bounded framework of    .Additionally, we can reduce the    space into a four-dimensional compact space.Thus, Itô's product rule: For SDEs of two scalar variables,  1 and  2 , given by the SDE of the product of the two variables is where    is the mean drift coefficient tensory,    is the additional drift coefficient tensor and     is the diffusion coefficient tensor.

Appendix C. Incompressibility constraint
To prove that the system of SDEs of    in equation (3.4) satisfies the incompressibility constraint, we take the trace on both sides of the    -SDE: Now, since   = 0, we have Further, since   = 0 by construction, it can be easily showed that

Appendix E. Galilean invariance
Now we demonstrate that the approach of closure modeling of the normalized anisotropic pressure Hessian () and viscous Laplacian () tensors satisfies Galilean invariance.The tensor  is modeled as  =  h   (E 1) where h is the pressure Hessian tensor in the principal frame of strain-rate tensor ().This h is obtained from data-driven closure as a function of b, also in principal reference frame.Note that h = h(, ,  2 ,  2 ), all four of which are either frame invariant or specifically defined in the principal reference frame and therefore h is unaltered by frame rotation.It is evident from equation (E 6) that the new tensor  also rotates by the same angles with

Figure 1 :
Figure 1: Flowchart to explain the behavior of velocity gradient tensor and its constituents in turbulence.

Figure 2 :
Figure 2: Statistics of  from DNS datasets of forced isotropic turbulent flows at different   : (a) global mean  as a function of   (in natural log scale); dashed line represents a linear least-squares fit of the data (  = −0.2+ 0.6 ln   ); and (b) variance  2  =  2 −  2 as a function of   (in natural log scale); dashed line represents a linear least-squares fit of the data ( 2  = −0.074+ 0.07 ln   ).

Figure 4 :
Figure 4: Evolution of  * statistics: (a) mean,  * and (b) standard deviation,   * , for the three models with different  * equations.The DNS statistics are marked by dashed lines.The time axis is in logscale.

Figure 5 :
Figure 5: PDF of standardized VG magnitude  * in: (a) linear-linear scale and (b) linear-log scale, for the three models.The black solid line with symbols represent the  * -PDF from DNS data.

Figure 7 :
Figure 7: Conditional mean trajectories in the - plane due to the inertial, pressure and viscous effects obtained using (a) DNS data and (b)    data-driven model.Background contours represent the speed of the trajectory at each point, given by the magnitude of the conditional mean velocity vector, | ṽ|.

Figure 8 :Figure 9 :Figure 10 :Figure 11 :
Figure 8: Evolution of  and  statistics in global normalized time * .Means: (a)  and (b)  ; second order moments: (c)  2 and (d)  2 ; third order moments: (e)  3 and (f)  3 ; fourth order moments: (g)  4 and (h)  4 for the three models with different  * equations.The dashed lines represent the DNS statistics.The time axis is in log-scale.

Figure 12 :
Figure 12: PDFs of absolute values of cosine of angles between vorticity vector and strain-rate eigenvectors (1 -most expansive, 2 -intermediate, 3 -most compressive).The solid lines are the PDFs obtained from DNS data.

Figure 13 :
Figure 13: PDFs of: (a) longitudinal component of velocity gradient tensor,  11 / √︃  2 11 , and (b) transverse component of velocity gradient tensor,  12 / √︃  2 12 , in log-linear scale obtained from the solutions of the three models.The solid line marked with symbols represent the PDFs obtained from DNS data.The dashed and dash-dotted lines represent the PDFs obtained from previous models -RDGF (Johnson & Meneveau 2016) and PIML (Tian et al. 2021), respectively.

Figure 14 :
Figure 14: PDFs of: (a) dissipation rate,       / √︁       , (b) enstrophy,       / √︁       , and (c) pseudodissipation rate,       / √︁       , in log-linear scale obtained from the solutions of the three models.The black solid line marked with symbols represent the PDFs obtained from DNS data; black dash-dotted line marks the PDFs for the initial field used in the model's simulations; dashed line represent the PDFs from the RDGF model of Johnson & Meneveau (2016).
Thus, = [ 1  2  3 ] (E 2) where   are the right eigenvectors of  corresponding to its eigenvalues   and   constitute the columns of the rotation matrix .Let us rotate the coordinate frame of the observer by certain angles, using a rotation matrix .Let the tensors and vectors in new reference frame be marked by .Then the tensor  becomes =     (E 3)and its eigenvectors also rotate by the same angles since  =     =⇒      =     =⇒    =     =⇒    =     where   =   (E 4)Since   constitute the columns of the rotated tensor  , we can say =   (E 5)Therefore, using equations (E 1) and (E 5), the pressure Hessian tensor in the new reference frame becomes,  =  h   =   h     =     (E 6) Recent studies(Parashar et al. 2020; Tian et al. 2021) have used tensor-basis neural network to model the unnormalized pressure Hessian (   ) and viscous Laplacian (   ) tensors in the    -equation as a function of    .Since    constitutes an unbounded space and the behavior of the tensors    and    is not necessarily invariant across turbulent flows of different Reynolds numbers, network-based modeling becomes essential.However, in our case (, ,  2 ,  2 ) form a bounded state-space and the conditional mean dynamics of h  and τ  in the bounded b  space is nearly unaltered with Reynolds number variation for a broad range of Therefore, the final SDE of  * for Model 2 is given by  * 2|,  −  * |,  2 )  * .(3.32) Here, the conditional variance values are obtained from DNS data of one   by discretizing the ,  space into 30 × 30 bins.The same conditional variance table is applicable when modeling turbulent flows of different Reynolds numbers, as evident from figure 3. Tis  * -model is weakly coupled with the    dynamics since it depends on , .Model 1 ensures that the constant diffusion coefficient captures the accurate global variance of  * , while model 2 enforces the accurate modeling of conditional variance of  * for a given , . Fnally, in Model 3 we propose an adjustment to model 2 such that both the conditional and global variances are satisfied.For this, we propose to shift the conditional variance based diffusion coefficient of model 2 by a constant value,  0 , as follows: * = − *  * + √︃ 2(

Table 1 :
Components of model based on DNS data.
2 ,  2   and    | =   τ |, ,  2 ,  2   .(3.37)The above coefficient values and the data-driven closure can be applied to model velocity gradient dynamics of incompressible turbulent flows irrespective of the Taylor Reynolds number.The final stochastic differential equation for  * in global timescale ( * =  ): models can be used to simulate the dynamics of VG magnitude for any Reynolds number, provided the corresponding parameter values are known.To reconcile between the two different timescales - used for    evolution and  * used for  * evolution -a simple, closed ordinary differential equation, for a given  .No closure is required for this equation.Finally, the Lagrangian evolution of the velocity gradient tensor,    , is obtained by multiplying  and    at different global time *   +  ,    =     (3.41)Ingeneral, the parameters in the above equation are Reynolds number dependent as shown in figure 2. For the case of   = 427, the parameter values are:  = 2.7493 ,   = 0.589655 ,  = 18.64173.(3.42)The above  * * * .Overall, the velocity gradient model presented here consists of three types of data-based components that are listed in table 1.The four-dimensional lookup tables for normalized pressure Hessian and viscous Laplacian tensors and the twodimensional table for conditional variance of  * can be used in modeling velocity gradient dynamics independent of Reynolds numbers; only the three scalar parameters of the model which are statistics of the VG magnitude are Reynolds number dependent and potentially generalizable in the future with universal scaling laws.4. Numerical procedure

Table 2 :
Third, fourth and sixth order moments of VG magnitude ( = √︁       ), longitudinal VG component ( 11 ), transverse VG component ( 12 ) from DNS data, model 1, model 2, model 3 and RDGF model of Johnson & Meneveau (2016).For each moment, the DNS value and the model's value closest to DNS are written in bold type font.

Table 3 :
Grid points    Details of forced isotropic incompressible turbulence data.