Capturing velocity gradients and particle rotation rates in turbulence

Turbulent fluid flows exhibit a complex small-scale structure with frequently occurring extreme velocity gradients. Particles probing such swirling and straining regions respond with an intricate, shape-dependent orientational dynamics, which sensitively depends on the particle history. Here, we systematically develop a reduced-order model for the small-scale dynamics of turbulence, which captures the velocity gradient statistics along particle paths. An analysis of the resulting stochastic dynamical system allows pinpointing the emergence of non-Gaussian statistics and non-trivial temporal correlations of vorticity and strain, as previously reported from experiments and simulations. Based on these insights, we use our model to predict the orientational statistics of anisotropic particles in turbulence, enabling a host of modeling applications for complex particulate flows.

Turbulent fluid flows exhibit a complex small-scale structure with frequently occurring extreme velocity gradients. Particles probing such swirling and straining regions respond with an intricate, shape-dependent orientational dynamics, which sensitively depends on the particle history. Here, we systematically develop a reduced-order model for the small-scale dynamics of turbulence, which captures the velocity gradient statistics along particle paths. An analysis of the resulting stochastic dynamical system allows pinpointing the emergence of non-Gaussian statistics and non-trivial temporal correlations of vorticity and strain, as previously reported from experiments and simulations. Based on these insights, we use our model to predict the orientational statistics of anisotropic particles in turbulence, enabling a host of modeling applications for complex particulate flows.
Turbulent flows show complex dynamics with a wide range of dynamically active scales [1][2][3], which play an important role for the dispersal of pollutants and aerosols in the atmosphere [4,5], the transport of microorganisms in the ocean [6][7][8][9], as well as the mixing of reactants in turbulent combustion [10,11]. The smallest turbulent scales, which are essentially independent of the boundaries and anisotropies of the large-scale flow [1], have a profound impact on the dynamics and collision rates of small suspended particles, like plankton in the ocean [6][7][8][9] as well as droplets and ice crystals in clouds [12][13][14][15][16]. Even in the simplest case of very small, neutrally buoyant particles, which passively follow the velocity field, highly nontrivial, shape-dependent rotational motion has been observed [17][18][19][20][21][22]. Theoretically, this intricate dynamics is not well understood.
The spinning and tumbling of particles immersed in a turbulent flow are determined by the complex interplay of particle shape and the small-scale structure of the turbulent flow field, as encoded in the gradients of the velocity field A ij = ∂u i /∂x j [23]. Because particle rotations are very sensitive to various small-scale features of turbulence such as non-Gaussian fluctuations, the local flow topology and, most importantly, the temporal correlation of strain and vorticity along Lagrangian trajectories, capturing this complex motion with theoretically insightful reduced-order models for turbulence so far remained elusive. The challenges in predicting these aspects of turbulent velocity gradients ultimately arise from the nonlinear, nonlocal, and dissipative dynamics of the governing Navier-Stokes equations.
Over the past years, a variety of reduced-order models for the velocity gradient statistics based on stochastic differential equations (SDEs) has been developed [23][24][25][26][27]. In these models, the effects of nonlocal pressure and viscous diffusion result in unclosed terms, to which diverse closure techniques have been applied. Closure theories range from models based on the coarse-grained velocity gradient as perceived by a tetrad of tracer particles * Michael Wilczek: michael.wilczek@ds.mpg.de [24,28], or the deformation of fluid elements [25] to functional closures based on Gaussian random fields [26], as well as combinations of these approaches [27]. The most advanced reduced-order SDE models successfully reproduce many of the characteristic geometric and statistical properties of the turbulent small scales [18,26,27,[29][30][31]. However, all current models struggle to capture important aspects of the temporal correlation of strain rate and rotation rate, which in particular leads to poor predictions for the orientational dynamics of particles immersed in turbulent flows.
Here, we develop a minimal model for the velocity gradients in turbulence, which enables profound theoretical insights. Starting from an exact statistical evolution equation, we systematically constrain its structure based on tensor function representation theory. By using an ensemble approach, we construct a physically consistent model which complies with important homogeneity constraints of turbulent fields. Based on an analysis of the associated Fokker-Planck equation, we establish a clear interpretation of its nonlinear dynamics. Specifically, we identify the dynamical mechanisms which control the degree of non-Gaussianity and temporal correlations of vorticity and strain. We test our predictions against high-resolution simulation results of fully developed turbulence, and show that our model successfully captures the temporal auto-correlations of rotation rate and strain rate. Coupled to the equations for the orientation dynamics of ellipsoidal particles, our model, furthermore, accurately reproduces the tumbling and spinning rates of particles in turbulent flows.
The foundation of our model is an exact, unclosed SDE for the one-point statistics of homogeneous isotropic turbulence, which is derived from the Navier-Stokes equations using stochastic calculus [26]. The resulting SDE takes the form Here A is the stochastic process corresponding to the velocity gradient field at the position of a fluid particle, and the conditionally averaged fields are evaluated at the same position. The tilde denotes the traceless part of arXiv:2005.02682v1 [physics.flu-dyn] 6 May 2020 the tensor, e.g. A 2 = A 2 − 1 3 Tr(A 2 )I . The first term on the right-hand side captures the non-linear local selfamplification of the velocity gradient, and includes the local isotropic part of the pressure Hessian H ij = ∂p ∂xi∂xj , which is obtained from the pressure Poisson equation ∆p = Tr(H) = −Tr(A 2 ) [32]. The conditional average H|A contains information about the mean nonlocal, deviatoric part of the pressure Hessian given a velocity gradient configuration A, which a priori depends on the full flow field due to the pressure Poisson equation. The conditional Laplacian ν∆A|A encodes viscous effects in the velocity gradient tensor evolution. The term dF is a Gaussian, temporally delta-correlated tensorial forcing, that is consistent with isotropy, homogeneity and incompressibility Tr( To close the conditional mean pressure Hessian and the conditional mean Laplacian terms, we express them as isotropic, tensor-valued functions of the symmetric and anti-symmetric part of the velocity gradient, the strain rate S and rotation rate W , respectively. Using tensor function representation theory, one can derive a complete and irreducible representation in terms of a small number of tensorial terms [33][34][35][36][37][38]. The individual tensorial terms are comprised of combinations of S and W , with coefficient functions that depend on isotropic invariants of S and W . The conditional mean traceless, symmetric pressure Hessian, for example, can be expressed as a linear combination of seven tensorial terms with appropriate coefficient functions [cf. SM]. Previous studies [26,39] showed that the most important features of the dynamics can already be captured by retaining terms up to the lowest possible order (i.e. second order in the pressure Hessian and up to first order in the viscous term) with constant coefficients. To enable analytical insights, we therefore truncate the general tensorial expansion and consider the closure This expression still contains five scalar parameters, which need to be further constrained. A general limitation of single-point closures is that they lack the possibility to include physical constraints which depend on information from the full field. For homogeneous turbulence, for example, the velocity gradient field fulfills the Betchov constraints [40] Tr(A 2 ) = 0 and Tr(A 3 ) = 0, which encode the balance of enstrophy and dissipation, as well as of their production. So far, velocity gradient models need careful calibration to fulfill these constraints [27]. An intriguing alternative to achieve a model which is physically consistent with homogeneous turbulence is to consider an ensemble of Lagrangian fluid elements that sample the full velocity gradient field. We then achieve consistency with the Betchov constraints by identifying spatial averages over the field with ensemble averages over the Lagrangian fluid elements. This can be used to derive analytical expressions from (1) for two of the parameters, which allows us to constrain our closure (2)-(3) [cf. SM]. One additional parameter can be fixed by nondimensionalizing the velocity gradient model with the Kolmogorov time scale τ η , which implies Tr(S 2 ) = 1/2. Thereby, the parameter space is reduced by three dimensions, and Betchov's homogeneity constraints are fulfilled by design. Besides the forcing amplitude, which we fix for the following considerations (see SM for the impact of the forcing amplitude), this leaves two free parameters: α, γ.
The impact of these free parameters on the nonlinear dynamics of the velocity gradient model can be revealed from the Fokker-Planck equation (FPE) corresponding to (1), which governs the evolution of the full probability density function (PDF) f (A; t) of the velocity gradient tensor (implied summation): (4) Here, Q ijkl (0) denotes the forcing covariance [cf. SM], and the nonlinear and linear drift terms are given by The parameter α controls the strength of the strain selfamplification in the velocity gradient dynamics. For a vanishing self-amplification (α = −1) and parameters determined such that the Betchov constraints are fulfilled, we find that (4) has an exact Gaussian solution [cf. SM].
Remarkably, even in this case, the FPE contains nonlinear drift terms. We demonstrate below that the strength of strain self-amplification controls the Gaussianity as well as important features of the small-scale topology of the predicted velocity gradient statistics (see Fig. 1).
Further analysis of the FPE shows that the single-time statistics is independent of the parameter γ for isotropic turbulence. In this case, the velocity gradient PDF is a function of the tensor invariants only, and one can read- . This result is related to a recently reported gauge symmetry of the pressure Hessian [41]. However, we show below that for the two-time statistics, and in particular the autocorrelations of vorticity and strain, the γ-term turns out to be crucial (see Fig. 2).
To determine appropriate values for the free parameters, we perform parameter scans and compare our model results with velocity gradient statistics obtained from direct numerical simulations (DNS) of the Navier-Stokes equation. We conducted simulations of homogeneous isotropic turbulence with 2048 3 grid points at a Taylorscale Reynolds number of R λ ≈ 509. For the subsequent data analysis, 25 statistically independent snapshots are taken into account. For the parameter scans, we numerically solve (1) using the Euler-Maruyama method [42] with a time step of ∆t = 0.0002. For all simulations shown here, we have integrated an ensemble of 10 5 Gaussian initial conditions for 5 × 10 6 time steps, which corresponds to 1000τ η in physical time, after an initial transient of 100τ η . Initial simulations of (1) revealed the occurrence of rare rogue trajectories exploring far-out regions of the phase space, which leads to non-convergent statistics and may introduce numerical instabilities in the determination of our parameters. We identified the second-order truncation of the unclosed terms as the origin of this shortcoming, which can be remedied by including a nonlinear damping term. The auxiliary term A, which is added to (3), is constructed to damp trajectories that diverge far from the ensemble mean and is negligibly small for the major, dynamically most relevant part of phase space; we set = −10 −8 ((Tr(W 2 ) + 1/2) 4 + (Tr(S 2 ) − 1/2) 4 ). Figure 1(a) illustrates how the strength of strain selfamplification controls the Gaussianity of the predicted velocity gradient statistics. As α deviates from −1, the single-component PDFs become non-Gaussian with increasingly heavy tails. For α = −0.6, the standardized PDFs of the velocity gradient components of our model agree well with DNS results within seven standard deviations (cf. Fig. 1(a)). The match is particularly good for the transverse components, for which our model captures the vanishing skewness and closely matches the DNS kurtosis of  Fig. 1(b) shows the joint PDF of the standardized invariantsQ andR for different values of α and from DNS. As α is tuned from −1.3 to −0.6, the joint PDF first extends along the left part of the Vieillefosse line, becomes symmetric for Gaussian statistics (α = −1), and finally extends along the right part of the Vieillefosse line. This corresponds to a shift of probability from flow regions with two compressive principal strain directions to regions with two extensional principal strain directions [43]. For α = −0.6, our model qualitatively captures the shape of the PDF as observed in DNS (lower right panel) and experiments [23], although the probability of velocity gradient configurations along the right part of the Vieillefosse line is underestimated. Nonetheless, since our model inherently fulfills the Betchov constraints, the mean of our modelR-Q PDF lies accurately at R = Q = 0 for all values of α.
Strain self-amplification also impacts another important aspect of the the small-scale topology, the alignment between the vorticity vector and the principal strain rate axes. Fig. 1(c) shows the PDFs of the cosine of the angle between the vorticity vector and the three eigenvectors of the strain-rate tensor. Our model (with α = −0.6) accurately captures the alignment of the vorticity with all three eigenvectors. In particular, it captures the wellknown preferential alignment of the vorticity with the eigenvector to the intermediate eigenvalue [23,44,45]. The alignment strength decreases with decreasing selfamplification and for α = −1, when the strain selfamplification vanishes and the model statistics are Gaussian, as expected, no preferential alignment is observed.
While we showed analytically that the γ-term has no effect on the single-time statistics, it determines tempo- ral correlations of the velocity gradients. This can be rationalized from the fact that γ(SW − W S) essentially rotates the principal strain vectors about the axis given by the vorticity vector with a rotation rate proportional to the vorticity magnitude [26]. This directly impacts the temporal correlation of velocity gradients and allows to precisely control them. Figure 2(a) compares temporal correlations of velocity gradients C ij (t 0 )C ij (t 0 + τ ) / C mn (t 0 ) 2 C pq (t 0 + τ ) 2 (implied summation) of our model to DNS results, where C ij = S ij or W ij . For the DNS results, we continued our simulations with 10 6 Lagrangian tracer particles and collected data from the statistically stationary state. For γ = −1.1, our model matches the vorticity autocorrelation very well. Importantly, it also captures the previously observed [46] shorter correlation time of the rate of strain compared to the rate of rotation, although differences occur in the shape of the correlation function. These results show in particular that the rotation of the strain-rate tensor as encoded by the γ-term is responsible for a decrease of the correlation time of the strain rate and an increase of the rotation-rate correlation time. Having established a model which captures the different temporal correlations of strain and vorticity along with the central non-Gaussian features of small-scale turbulence, we can use it to predict the tumbling and spinning rates of Lagrangian particles. To this end, we couple our model to Jeffery's equation [47], which describes the orientational dynamics of ellipsoidal particles: Here, p denotes the particles' symmetry axis, and λ is the particles' aspect ratio, i.e. the ratio of the length along the symmetry axis to the length perpendicular to it. The rotation of an axisymmetric particle can be decomposed into spinning (rotation around the symmetry axis) and tumbling (rotation around an axis perpendicular to the symmetry axis) [17], with the squared spinning rate Ω 2 p = 1 2 ω · p 2 , where ω is the vorticity, and the squared tumbling rateṗ iṗi . In Fig. 2(b) the nondimensionalized mean square tumbling and spinning rates as predicted by our model are shown for different values of γ as a function of the particles' aspect ratio. Fig. 2(b) shows that especially the tumbling rates of disk-like particles increase when the temporal correlations are modified by increasing the magnitude of the coefficient γ. When our model exhibits the most realistic correlation times, i.e. for γ = −1.1, the rotation rates predicted by our model agree very well with the ones observed in our DNS and literature [17,20] for the full range of particle shapes.
In particular, our model predicts the high tumbling rates of disk-like particles observed in DNS and experiments. The comparison of the results for different values of γ in Fig. 2(a) and Fig. 2(b) indicates that the realistic autocorrelation times of our model are crucial for an accurate prediction of tumbling rates of suspended particles.
In summary, we have analyzed the dynamics and statistics of velocity gradients in turbulence in the framework of a minimal, physically consistent reduced-order model. Our combined analytical and computational analysis showed that strain self-amplification controls the non-Gaussianity as well as the small-scale topology of the velocity gradient dynamics, and identified the rotation of the strain eigenvectors by the vorticity as the major factor in determining the temporal correlations of velocity gradients. As a result, we obtained a reduced-order model for the small scales of turbulence that captures the different correlation times of strain and vorticity in turbulence. We showed that the reduced-order model can be used to accurately predict the orientational statistics of suspended anisotropic particles, enabling a host of modeling applications for complex particulate flows.
Based on tensor function representations theory and the systematic implementation of physical constraints, our closure approach explicitly uncovers the general tensorial structure of the unclosed terms, which also provides a firm foundation for future advancements. For example, we expect that the inclusion of higher-order terms and coefficient functions which depend on velocity gradient tensor function invariants will lead to further quantitative improvement. Machine learning approaches [48] could turn out to be instrumental in achieving such improved parameterizations.
Supplemental material: Capturing velocity gradients and particle rotation rates in turbulence To obtain a closed statistical evolution equation for the velocity gradient tensor (Eq. (1) in the main text), the conditional mean pressure Hessian and the conditional mean viscous Laplacian terms need to be specified. Both terms are tensor-valued functions of the velocity gradient tensor, whose structure is constrained by statistical isotropy. For a second-order tensorial function M which depends on the velocity gradient tensor A, statistical isotropy implies where Q is an arbitrary orthogonal matrix. Additional constraints arise when M is symmetric or anti-symmetric and depends on symmetric or anti-symmetric tensorial arguments. The general structure of the tensorial function can then be determined using tensor function representation theory [33][34][35]37], which expresses the function as a linear combination of tensorial terms (also called generators or form-invariants) with scalar coefficient functions that may depend on scalar tensor invariants of A. In the mathematical literature representations for two kinds of functions exist: For polynomial [49][50][51] and for general tensor functions [33][34][35]37]. For the closure, we use representations for general tensor functions, as these do not require assumptions about the functional form and also contain in general fewer terms than the representations of polynomial functions [37,38]. Based on this, the conditional traceless and symmetric pressure Hessian takes the form: where the B (n) are: and the coefficients b (n) are functions of the scalar invariants Tr(S 2 ), Tr(S 3 ), Tr(W 2 ), Tr(SW 2 ), Tr(S 2 W 2 ) and Tr(S 2 W 2 SW ). Note that there are six independent scalar invariants of S and W (and hence A), i.e. one independent scalar invariant more than typically discussed in the turbulence literature [1,23,43]. This difference is rooted in the fact that the number of invariants in an irreducible functional basis of a tensor is not necessarily equal to the number of independent entries of the tensor, which is typically considered. Furthermore, an irreducible functional basis is not necessarily minimal, meaning that a different representation with fewer elements might exist [36,37]. The irreducibility of the functional basis presented here has been explicitly proven in [36]. Similar to the conditional pressure Hessian, one can construct the most general tensor function representation for the viscous Laplacian term. Since the function representations only hold for symmetric or anti-symmetric tensor functions, one first has to decompose the Laplacian of A into the Laplacian of its symmetric part S and anti-symmetric part W . We find:

II. DERIVATION OF THE ADAPTIVE COEFFICIENTS
The closure for the velocity gradient model, as given by eqs. (2) and (3) in the main text, depends on five coefficients. The number of independent coefficients can be reduced by non-dimensionalizing the equation of motion (Eq. (1) in the main text) by the Kolmogorov time scale τ η , which implies Tr(S 2 ) = 1/2. Here and in the following · denotes an ensemble average. The Betchov constraints for homogeneous turbulence, Tr(A 2 ) = 0 and Tr(A 3 ) = 0, can be used to further reduce the number of parameters to a total of two independent parameters, α and γ. For the three model coefficients fixed in this way, analytical expressions can be found. They are constructed by deriving stochastic differential equations from the model SDE for each constrained quantity ( Tr(S 2 ) , Tr(W 2 ) , Tr(A 3 ) ) using Itô's formula. After averaging, one finds for example the following equation for Tr(S 2 ) : and accordingly for the other constraints. Here, N ij and L ij are the non-linear and linear terms of the FPE (eqs. (5), (6) in the main text), respectively, and Q ikjl (0) denotes the forcing covariance, detailed in section IV. Combined with Eq. (1) in the main text, the averages can be evaluated. As for constant constraints the left-hand sides vanish, we can solve for three of the parameters and obtain where a = α Tr(A 2 S 2 ) + 6 Tr(SW 2 ) Tr(A 2 S) + 2 Tr(A 2 ) Tr(A 2 S) and is the auxiliary damping term, discussed in the main text. These expressions are evaluated in every simulation timestep to update the coefficient values. Fig. S1(a) shows an example of the parameter values which are dynamically obtained throughout a simulation. Fig. S1(b) demonstrates that the constraints are indeed fulfilled to very good precision throughout the simulation.

III. SIMULATION DETAILS
The integration of the SDE (Eq. (1) in the main text) has been implemented in Python using the Euler-Maruyama method with a time step size ∆t = 0.0002. For all simulations shown here, we have integrated an ensemble of 10 5 Gaussian initial conditions for 5 × 10 6 time steps, which corresponds to 1000τ η in physical time, after an initial transient of 100τ η .
Furthermore, an integration of Jeffery's equation (Eq. (7) in the main text) has been implemented to obtain the rotation rates of Lagrangian particles. Initially, 10 4 randomly oriented symmetry axes of the particles were created for each particle aspect ratio and then evolved with Jeffery's equation using the Euler method, where S and W were taken from our model simulation. The simulation of the particles was done for logarithmically spaced aspect ratios ranging from 0.02 to 40.

IV. IMPACT OF THE FORCING AMPLITUDE
The forcing term dF ij = f D ijab dW ab is based on tensorial Gaussian white noise dW ab with statistically independent components, and has a covariance which is designed to be consistent with isotropy, homogeneity and tracelessness of A [26,29] in the model SDE (Eq. (1) in the main text). It can be motivated from a stochastic Gaussian force term in the Navier-Stokes equation [26] and assures stationary statistics by counteracting dissipative terms of the model. For small values of f (f ≈ 0.05) the deterministic dynamics are only little perturbed, leading to more excursive trajectories. Larger values of f (f ≈ 0.2) regularize the dynamics by increasing the influence of the Gaussian forcing. For large forcing amplitude f , the forcing dominates over the deterministic dynamics of the model, leading to Gaussian statistics. Therefore, as f increases, the tails of the component PDFs become less pronounced and the preferential alignments become weaker (see Fig. S2(a) and (c). Furthermore, the probability density along the Vieillefosse line in the lower left quadrant of the R-Q PDF decreases (see Fig. S2(b)). Based on a parameter scan of the forcing amplitude f and comparison with direct numerical simulations, we set f = 0.08. The strain self-amplification, controlled by the parameter α, determines the degree of Gaussianity of the modeled velocity gradient statistics. In fact, for α = −1, the Gaussian velocity gradient PDF is an exact solution to the Fokker-Planck equation (Eq. (4) in the main text), which is consistent with the Betchov constraints. Here, R −1 ijkl = 1 Tr(S 2 ) [4δ ij δ kl + δ il δ jk ] is the pseudo-inverse of the velocity gradient covariance tensor R ijkl = A ik A jl . The delta function stems from the incompressibility condition of the velocity field, Tr(A) = 0.

V.2. Independence of the single-time statistics of γ
Here we show that the single-time statistics is independent of γ under the assumption of statistical isotropy. This allows us to control the temporal correlations in our model through γ without changing the single-time statistics. For isotropic statistics, the velocity gradient PDF is a scalar function of the six isotropic invariants of A [37]: f (A) = F Tr(S 2 ), Tr(S 3 ), Tr(W 2 ), Tr(SW 2 ), Tr(S 2 W 2 ), Tr(S 2 W 2 SW ) . (S11) For the γ-term in the Fokker-Planck equation (Eq. (4) in the main text), we therefore obtain where the X n are the six scalar invariants listed above. One can now straightforwardly compute that the tensor contraction γ(S ik W kj − W ik S kj ) ∂Xn ∂Aij = 0 term by term, which shows that the γ-term vanishes in the Fokker-Planck equation.