ABSTRACT

INTRODUCTION
Driven by numerous applications of particle laden turbu- lence, e.g.  dust in stars or the rain formation in atmospheric clouds, there are many numerical and theoretical studies in the literature, see e.g. (Brownlee, 1985; Shaw, 2003). Under- standing the motion of inertial particles is therefore essential to get insight into these phenomena and to develop sound mod- els.  The aim here is to determine the dynamical properties of inertial particles in isotropic turbulence using finite time tessellation-based measures.  The particle velocity is crucial to understand particle dynamics, i.e. convergence, divergence of particles and the vortical or swirling motion of the particle clouds. The challenging task is that we only know the velocity of the particles at discrete points in space, i.e. at the particle position. Voronoi tessellation of the particle position has previously been used, e.g. in Monchaux et al. (2010), for analyzing pref- erential concentration of inertial particles.  In our recent work, Oujia et al. (2020), we proposed a Lagrangian approach to de- termine the divergence of particle velocity using the volume change rate of the Voronoi cells.  Here we propose a further new Lagrangian approach using Delaunay tessellation, to de- termine the spatial velocity derivatives numerically by consid- ering subsequent time instants. In particular the trace of the ve- locity gradient tensor, corresponding to the divergence, which yields information about sources and sinks and its antisym- metric part, corresponding to the curl, which characterizes the rotation can thus be computed.  The helicity, i.e.  the scalar product of vorticity and velocity, can be likewise calculated and swirling motion of particle clouds can consequently be quantified.  The dynamical properties of the particle velocity can be assessed and detailed statistics of the divergence and of the rotation can be performed. Moreover, we will address the presence of extreme val- ues in these quantities which is reflected in heavy tails in the probability distribution functions (PDFs). Extreme events are a generic feature of turbulence, see ee.g., Moffatt (2021), and make its prediction particularly difficult, e.g., for tornadoes, large floods and other extreme events with devastating impact. In the present work we will show that the divergence and ro- tation values of the particle velocity likewise have extreme in- tensities, much larger than their standard deviation. The remainder of the manuscript summarizes first the DNS flow data we analyze.  Then we describe the method to determine the divergence and rotation based on Delaunay tessellation of the particle positions. Numerical results and ex- treme divergence and rotation values are presented afterwards and finally some conclusions are drawn.
FLOW DATA AND METHODS
Particle position and velocity data are generated by DNS of particle-laden homogeneous isotropic turbulence, presented in Matsuda et al. (2014).  The incompressible Navier-Stokes equations are solved in a 2π–periodic cube with a fourth or- der finite difference scheme.  Statistically stationary flow is obtained by forcing at large scales. Uniformly distributed dis- 1 crete particles are then injected into the fully developed flow and are tracked in the Lagrangian framework. Maxey's model (Maxey, 1987) for inertial heavy point particles with Stokes drag is used and the inertial dynamics is controlled by the Stokes number, St = τρ/τη, where τρ is the particle relaxation time and τ n the Kolmogorov time. The equation of particle velocity v p j is given by
where u p j is the fluid velocity at particle position x p j.  The subscript p denotes the quantity at the position of a particle, and the subscript j denotes the particle identification number. State of the art high resolution DNS with N 3 = 512 3 grid points is performed for the Taylor-microscale Reynolds num- ber Re λ = 204, where Re λ ≡ ι ′ λ / ν, ν is the kinematic viscos- ity, and λ  is the Taylor microscale. The number of particles N is 1.5 × 10 7 and the considered Stokes numbers are St = 0.5, 1, 2 and 5. Particles with different Stokes numbers, including St = 0 were tracked in an identical turbulent flow. We apply 3D Delaunay tessellation (Delaunay, 1934) us- ing the Quickhull algorithm (Barber et al. , 1996) to the par- ticle positions at two consecutive time instants.  To define a cell corresponding to a particle, we use the dual graph of the Delaunay tessellation.  For stability reasons, instead of using the circumcenter of the Delaunay cell, as done for the Voronoi tessellation, we use the center of gravity to define the vertices of the cell. To compute the divergence of the particle velocity B ( V p ), we consider the local number density averaged over a cell, which is the inverse of the corresponding volume.  Using the fact that particles satisfy the conservation equation of the den- sity n p , D i n p = -n p V · V p where D i = d + v · V is the La- grangian derivative and considering two time instants t k and t k+1 = t k + Δt of the Delaunay tessellation with time step Δt, we can determine the volume change. Thus we obtain (Oujia et al., 2020),
The curl of the particle velocity, is defined by computing the circulation of the velocity field of particles over a cell. This can also be expressed as the divergence of the velocity of the fluid which has been rotated in a direction π/2 with respect to the direction of the curl. We define ν + = A x ν, ν + = A x ν and v ± = A z v where A x , A y and A z are rotation matrices around the different axes.  We obtain that the curl of the particle velocity C ( v p ) is given by
The accuracy and reliability of the method to compute the divergence, curl and furthermore the velocity gradient tensor of the particle velocity are studied in Oujia et al. (2022).  We find that the spatial order of the method is the inverse of the dimension of the space.
RESULTS
Figure 1 shows two-dimensional cuts of the particles colored with the divergence B(v p ) (top) and the enstrophy ||& (Vp)||3 (bottom) for St = 0 (left) and 1 (right). The particles for St = 0, i.e. fluid particles, are randomly distributed because the flow is incompressible. We can deduce that for St = 0 the divergence values different from zero are numerical errors due to the sampling.  Their distribution is homogeneous and we find mostly values close to zero, except in some regions, which are characterized by strong enstrophy values.  The enstrophy of the fluid particles (bottom, left) shows a non homogeneous distribution and is spotty, characteristic for intermittent quan- tities. In contrast to fluid particles, inertial particles (Figure 1, right), here for St = 1, are clustered. We can observe the pref- erential concentration of particles, corresponding to clusters in regions of low fluid vorticity and voids in regions of large fluid vorticity. This observation is explained by the fact that inertial particles are ejected from high fluid vorticity regions due to the centrifugal force.  We also find the high convergence and divergence value are observed close to each other.  We con- jecture that the particles are crossing.  In a sufficiently large region of low fluid vorticity, we can observe that divergence and enstrophy values of inertial particles are close to zero. PDFs of divergence and rotation are shown in Figure 2. As the turbulence is isotropic, we do not distinguish between the different vector components of C ( V p ), and we compute the histograms of all these components together. For fluid particles in the continuous setting, the divergence of the fluid velocity vanishes exactly, whereas in the discrete setting D ( V p ) differs from zero due to the sampling and the computation of the dis- crete divergence. For inertial particles, the PDFs of the diver- gence and of the rotation deviate from those for fluid particles. This nicely illustrates that for large Stokes numbers heavy tails and extreme values can be observed in both quantities, going up to more than 50 and 40 times of the standard deviation for the divergence and the rotation, respectively. The PDFs of the divergence and of the rotation have much heavier tails com- pared to the fluid vorticity, and the extreme values increase significantly with the Stokes number.  This is quantified fur- ther in Table 1, where each of the variance and flatness of the divergence S ( V p ) and curl C ( V p ) is given as a function of the Stokes number. We observe that the variance of the diver- gence increases as the Stokes number increases.  The flatness of the divergence first increases with St with a maximum for St = 1 and then decreases. The large flatness value is a sign of strong non-Gaussianity and the spatially intermittent behav- ior of the divergence for St = 1. These observations suggest that the global number of particles with large divergence val- ues increases with the Stokes number. Note that for St > 1 tails become shallower with increasing St when the PDFs are nor- malized by the standard deviation, not shown here. We recall that the divergence of the particles is induced by the vorticity of the fluid flow via the centrifugal force.  The decay of the flatness for St ≥ 2 can be explained by the particle's inertia. In fact, heavy particles are less subjected to the small scale fluc- tuations of the fluid vorticity, so that the particles follow less the fluid flow and thus have less tendency to accelerate sud- denly.  The variance of the curl decreases from St = 0 to 0.5 and then increases as the Stokes number increases.  The flat- ness decreases from St = 0 to 0.5 and then increases with St with a maximum for St = 2 and then decreases.  The reduc- tion of the variance and flatness of the curl between St = 0 and 0.5 can be explained by the fact that the particles in the high velocity regions in the fluid domain are ejected due to the centrifugal force. However, for a small St the particle velocity at a given position is still close to that of the fluid.  This ex- plains the reduction of the curl value between St = 0 and 0.5 as well as that of the variance. The increase of the variance of the curl value for larger Stokes numbers could be explained by the higher inertia. Due to the sling effect, the particle trajecto- ries may cross each other, or the particles may enter high fluid vorticity regions and experience strong vortical motion. Figure 3 shows the PDF of the relative helicity of parti- cle velocity for different Stokes numbers and randomly dis- tributed particles. The relative helicity of the particle velocity is defined as
We can observe that for St < 1 there are two maxima at :H p = ±1, which means that we have a higher probability that the vorticity is oriented in the same or opposite direction as the particle velocity. This corresponds to helical motion, which is a signature of coherent structures in the particle laden flow. In contrast we find for St = 5 two minima at H p = ±1 and a maximum close to 0 which means that the vorticity has higher probability to be orthogonal to the particle velocity, which in- dicates incoherent particle flow motion. We can observe a tran- sition from a convex to a concave shape at St = 2. Note that the asymmetry of the distribution can be explain by the presence of helicity asymmetry in the fluid flow due to initial forcing.
CONCLUSIONS
We analyzed high resolution DNS data of particle laden turbulence considering different Stokes numbers.  We pro- posed a semi-Lagrangian method to determine the divergence and rotation of the particle velocity.  To this end we consid- ered the time change of the Delaunay tessellation of the parti- cle positions and computed the spatial derivatives in the La- grangian frame.  Consequently the divergence, the curl and likewise the helicity of the particle motion could be quantified. We showed that the divergence and curl PDFs exhibit heavy
tails and we found that for large Stokes numbers extreme val- ues are above 40 times the corresponding standard deviations. In future work, we plan to perform multiscale analysis using multiresolution tessellations. ACKNOWLEDGEMENTS: T.O. and K.S. acknowledge partial fund- ing from the Agence Nationale de la Recherche (ANR), grant ANR- 20-CE46-0010-01. Centre de Calcul Intensif d'Aix-Marseille is ac- knowledged for granting access to its high performance computing resources.  K.M. acknowledges financial support from JSPS KAK- ENHI Grant Number JP20K04298.  The DNS data analysed in this project were obtained using the Earth Simulator supercomputer sys- tem of JAMSTEC.