2. Direct numerical simulation data

We analyse particle position and velocity data obtained by 3-D DNS of particle-laden homogeneous isotropic turbulence presented in Matsuda et al. (Reference Matsuda, Onishi, Hirahara, Kurose, Takahashi and Komori2014). The DNS was performed for a cubic computational domain with a side length of $2{\rm \pi}$. The incompressible Navier–Stokes equation was solved using a fourth-order finite-difference scheme. Statistically steady turbulence was obtained by forcing at large scales, i.e. for $k < 2.5$, where $k$ is the wavenumber. Discrete particles were tracked in the Lagrangian framework. The equation of particle motion is given by

(2.1)\begin{equation} d_t {\boldsymbol v_p}_j = - \frac{{\boldsymbol v_p}_j - {\boldsymbol u_p}_j}{\tau_p} , \end{equation}

where ${\boldsymbol v_p}_j$ is the particle velocity vector and ${\boldsymbol u_p}_j$ is the fluid velocity vector at particle position ${\boldsymbol x_p}_j$. Note that subscript $p$ denotes the quantity at the position of a particle (e.g. ${\boldsymbol u_p}_j \equiv {\boldsymbol u}({\boldsymbol x_p}_j)$), and the subscript $j$ denotes the particle identification number ($\,j=1,\ldots ,N$, where $N$ is the total number of particles). The Stokes drag was assumed for the drag force. The relaxation time $\tau _p$ is independent of the particle Reynolds number. See Onishi, Baba & Takahashi (Reference Onishi, Baba and Takahashi2011), Matsuda et al. (Reference Matsuda, Onishi, Hirahara, Kurose, Takahashi and Komori2014) and Matsuda & Onishi (Reference Matsuda and Onishi2019) for details on the computational method.

The number of grid points for the flow field is $N_g^{3} = 512^{3}$. The kinematic viscosity is $\nu =1.10\times 10^{-3}$. The r.m.s. of velocity fluctuation $u'$, energy dissipation rate $\epsilon$ and Taylor-microscale-based Reynolds number $Re_\lambda$ ($\equiv u' \lambda /\nu$, where $\lambda$ is the Taylor microscale) of obtained turbulence were $u'=1.01$, $\epsilon =0.344$ and $Re_\lambda = 204$. Note that Matsuda et al. (Reference Matsuda, Onishi, Hirahara, Kurose, Takahashi and Komori2014) confirmed that $Re_\lambda = 204$ is large enough to be representative of high-Reynolds-number turbulence at $k\eta > 0.05$, where $\eta \equiv \nu ^{3/4} \epsilon ^{-1/4}$ is the Kolmogorov scale (i.e. $\eta =7.90\times 10^{-3}$). It should be also noted that the number of grid points is sufficiently large for resolving the turbulent flow so that $k_{max} \eta$ reached 2.03, where $k_{max} \equiv N_g/2$. The number of particles $N$ is $1.5 \times 10^{7}$ and the Stokes number $St$ ($\equiv \tau _p/\tau _\eta$, where $\tau _\eta \equiv \nu ^{1/2}\epsilon ^{-1/2}$) is $St = 0.05$, 0.1, 0.2, 0.5, 1, 2 and 5. Particles with different Stokes numbers were tracked in an identical turbulent flow. After the turbulent flow had reached a statistically steady state, particles were randomly seeded satisfying a Poisson distribution. Effects of the reaction of particle motion to fluid flow and the interactions between particles were neglected because these effects are typically small in the time scale of $\tau _\eta$ for sufficiently dilute particles such as cloud droplets in atmospheric turbulence (Matsuda & Onishi Reference Matsuda and Onishi2019). In the DNS data the particle number in a volume of $\eta ^{3}$ is about 0.03 and close to the conditions in atmospheric clouds. For such particle densities, the point-particle approximation and the one-way coupling are justified for most of the particles. For this study, we additionally consider data of randomly distributed particles with fluid velocity at the particle positions. These data are analysed as the fluid particle case; i.e. $St=0$. Note that all statistical results are averaged over 10 snapshots at time $21 T_0 \le t \le 30 T_0$, where $T_0 \approx L_0/u'$ and $L_0$ is the representative length which gives the domain length of $2{\rm \pi} L_0$ (i.e. $L_0=1$).

Figure 1 shows two-dimensional (2-D) cuts of the particle number density for different Stokes numbers. We can clearly observe void areas for $St=1$. For $St < 1$, the void areas are less clear. For $St > 1$, they are larger but less pronounced than for $St=1$.

Figure 1.

Spatial distribution of the particles coloured with the divergence $\mathcal {D}_p$ for (a) $St=0.05$, (b) 0.2, (c) 1 and (d) 5 at time $t = 24 T_0$ for a slice of thickness $4\eta$.

3. Voronoi tessellation to compute the particle velocity divergence

To understand the dynamics of inertial particles, in particular the clustering, we consider the particle density $n$ in the continuous setting, as shown in many studies (e.g. Maxey Reference Maxey1987). It satisfies the conservation equation which can be written with the Lagrangian derivative, $D_t = \partial _t + {\boldsymbol v}\boldsymbol {\cdot } \boldsymbol {\nabla }$, as

(3.1)\begin{equation} D_t n = -n \boldsymbol{\nabla} \boldsymbol{\cdot} {\boldsymbol v}. \end{equation}

This illustrates that we need to know the divergence of the particle velocity, which is a source term on the right-hand side of (3.1). The problem to determine $\boldsymbol {\nabla } \boldsymbol {\cdot } {\boldsymbol v}$ is that we only know the discrete particle distribution, and, hence, we only know the particle velocity at the particle positions, i.e. ${\boldsymbol v_p}_j = {\boldsymbol v}({\boldsymbol x_p}_j)$, but not anywhere else. In this study we apply the Voronoi tessellation to overcome this problem.

The Voronoi tessellation (or diagram) is a technique to construct a decomposition of the space, i.e. the fluid domain, into a finite number of Voronoi cells $C_i$. When a finite number of particles $p_i$ are dispersed in space, a Voronoi cell $C_i$ is defined as a region closer to a particle than other particles. The volume of a Voronoi cell is referred to as Voronoi volume and denoted by $V_{p_i}$. The cell $C_i$ can be interpreted as the zone of influence of the particle $p_i$. The larger the number of particles in a given volume, the smaller the Voronoi volume. The diagram will allow us to identify particles inside clusters (corresponding to small cells) and particles inside void regions (corresponding to large cells). A survey on Voronoi diagrams, a classical technique in computational geometry, can be found in Aurenhammer (Reference Aurenhammer1991).

We apply 3-D Voronoi tessellation to the DNS data using the Quickhull algorithm provided by the Qhull library in python (Barber, Dobkin & Huhdanpaa Reference Barber, Dobkin and Huhdanpaa1996), which has a computational complexity of $O(N \log (N))$.

Here we propose a method to compute the divergence of the particle velocity ${\mathcal {D}}\equiv \boldsymbol {\nabla } \boldsymbol {\cdot } {\boldsymbol v}$ in a discrete manner. Dividing (3.1) by the particle density $n$, we obtain

(3.2)\begin{equation} {\mathcal{D}} = -\frac{1}{n} D_t n. \end{equation}

To calculate the Lagrangian derivative of $n$, we define the local number density $n_p$ as the number density averaged over a Voronoi cell, which is given by the inverse of the Voronoi volume $V_p$; i.e. $n_p = 1/V_p$. We consider two time instances, $t^{k}$ and $t^{k+1} = t^{k} + {\rm \Delta} t$, where ${\rm \Delta} t$ is the time step and the superscript $k$ denotes the discrete time index.

The time-averaged number density change in the period of ${\rm \Delta} t$ is given by $\overline {D_t n}^{{\rm \Delta} t} = (n_p^{k+1}-n_p^{k})/{\rm \Delta} t$, where $\bar {n}^{{\rm \Delta} t} = (n_p^{k+1}+n_p^{k})/2 + O({\rm \Delta} t)$. Thus, we obtain a finite-time discrete divergence of the particle velocity in the period of ${\rm \Delta} t$:

(3.3)\begin{equation} {\mathcal{D}}_p = -\frac{2}{{\rm \Delta} t} \frac{n^{k+1}_p-n^{k}_p}{n^{k+1}_p+n^{k}_p} + O({\rm \Delta} t) = \frac{2}{{\rm \Delta} t} \frac{V^{k+1}_p-V^{k}_p}{V^{k+1}_p+V^{k}_p} + O({\rm \Delta} t) \end{equation}

which depends on the choice of ${\rm \Delta} t$ and $N$. This shows that the divergence of the particle velocity can be estimated from subsequent Voronoi volumes, given that the time step is sufficiently small and the number of particles sufficiently large. To obtain the subsequent Voronoi volumes $V_p^{k+1}$, the particle positions were linearly advanced by ${\boldsymbol v_p}$; i.e. ${\boldsymbol x_p}^{k+1} = {\boldsymbol x_p}^{k} + {\boldsymbol v_p}{\rm \Delta} t$. The time step was set to ${\rm \Delta} t = 10^{-3}$, a value which is sufficiently small. A study of the influence of the step size and the number of particles can be found in appendix A. The results in the appendix show that neither the time step nor the particle number have a significant influence on the statistics of fluid particles when dividing or multiplying the time step by a factor two, or dividing the number of particles by two. However, for inertial particles the situation changes, except for the time step reduction by a factor two, while keeping $N$ fixed, for which we obtain an almost identical PDF of $\mathcal {D}_p$. This confirms that ${\rm \Delta} t$ is sufficiently small. Doubling the time step with fixed $N$ leads to a sudden drop of the tails in the PDF, which can be explained with the Courant–Friedrichs–Lewy (CFL) condition. The particle number dependence for fixed ${\rm \Delta} t$, which is related to the $N$ dependence of the mean separation length, shows that for $N/2$, the tails in the PDF become lighter, but the extreme values are almost unchanged. For details, we refer to figure 7 in appendix A.

In order to consider the difference with the fluid motion, we also compute the discrete divergence of the fluid velocity at Voronoi cells using (3.3). In this case, $V_p^{k+1}$ is obtained from fluid particle positions ${\boldsymbol x_p}^{k+1} = {\boldsymbol x_p}^{k} + {\boldsymbol u_p}{\rm \Delta} t$ instead.

Figure 2 shows a 2-D particle distribution identical to figure 1(c) and a magnified view with a corresponding Voronoi diagram. We can observe that Voronoi cells of particles in clusters are relatively small, while those corresponding to particles outside clusters are large.

Figure 2.

Voronoi tessellation generated by particles for $St = 1$. Two-dimensional particle distribution in a slice of thickness $4\eta$ (a), corresponding to a zoom of figure 1(c). A magnified view with Voronoi tessellation (b).

4. Numerical results

In the following we present numerical results for inertial particles in isotropic turbulence considering seven different Stokes numbers. Fluid particles are likewise analysed, which allow us to understand the statistical properties and to assess the numerical precision of the divergence approximation in (3.3).

We compute the Voronoi volume for different Stokes numbers and we compare statistical properties with those obtained for random particles. Figure 3(a) shows the PDF of the Voronoi volumes $V_p$, which is normalised by the mean volume $\overline {V_p}=(2{\rm \pi} )^{3}/N$. For randomly distributed particles, the PDF of the Voronoi volume becomes a gamma distribution (Ferenc & Néda Reference Ferenc and Néda2007). For the 3-D case, the PDF of the Voronoi volume is given by $\varGamma (5,1/5)$, where $\varGamma (k,\theta )$ corresponds to the PDF $f_{Vp}(x)=\varGamma (k)^{-1}\theta ^{-k}x^{k-1}\exp (-x/\theta )$ and where $k$ and $\theta$ are shape and scale parameters, respectively. In figure 3(a) we have confirmed that the PDF for randomly distributed particles in our results agrees well with the gamma distribution of $\varGamma (5,1/5)$. We can see that, as the Stokes number increases and is getting closer to $1$, the number of small Voronoi cells for $V_p/\overline {V_p} \lesssim 0.5$ increases, and then decreases after exceeding $St = 1$. The number of the large Voronoi cells increases as the Stokes number increases and it stabilises. The PDF of the Voronoi volume is often used to determine ‘cluster cells’ and ‘void cells’. Monchaux et al. (Reference Monchaux, Bourgoin and Cartellier2010) adopted the intersection of the PDF of the Voronoi volume and the gamma distribution as thresholds. A Voronoi cell smaller than the smaller threshold (the left intersection) is defined as a cluster cell, and, similarly, a Voronoi cell larger than the larger threshold (the right intersection) is defined as a void cell. In our results, the threshold to determine cluster cells is $V_p/\overline {V_p} \sim 0.5$.

Figure 3.

Probability distribution functions of (a) Voronoi volume normalised by the mean and (b) divergence for particle velocity for different Stokes numbers and for randomly distributed particles ($St=0$). The inset of (a) compares the PDF for randomly distributed particles with the gamma distribution. The inset of (b) shows the probability distribution functions for positive divergence values in log-log representation.

The divergence of the particle velocity is then computed using (3.3). Figure 3(b) shows the PDF of the divergence of the particle velocity. The probability distribution functions are almost symmetric, centred around 0 with stretched exponential tails.

Table 1 shows the variance and flatness of the divergence ${\mathcal {D}}_p$ as a function of the Stokes number. The variance of the divergence increases as the Stokes number increases. For small Stokes numbers ($St \le 0.2$), the variance is strongly reduced. The flatness of the divergence first increases with $St$, reaches its maximum value around $St=0.5$ and then decreases again. This observation suggests that the tails decay faster for larger $St$. To understand which part of the divergence of the particle velocity is due to a geometrical error, we added the PDF of the divergence of the fluid velocity ($St=0$). Note that for fluid particles in the continuous setting, the divergence of the fluid velocity vanishes exactly, whereas in the discrete setting ${\mathcal {D}}_p$ differs from zero, because the deformation of a Voronoi cell is not exactly the same as the deformation of a fluid volume in the continuous setting. For the divergence of the fluid velocity, we find indeed values ranging between $-150$ and $+150$. Hence, we can deduce that for Stokes numbers less than $0.2$, the divergence is mostly due to a geometrical effect, but for larger Stokes numbers, physical effects predominate. The numerical precision of the divergence computation has been assessed in appendix A. In appendix B the geometrical effect on the divergence is further discussed for the case of randomly distributed particles.

Table 1.

Variance and flatness of divergence ${\mathcal {D}}_p$ as a function of the Stokes number.

To get further insight into the cluster formation, we plot the joint PDF of the divergence ${\mathcal {D}}_p$ and the Voronoi volume $V_p$ normalised by its mean. Figure 4 shows the joint PDF for the cases of $St=1$ and 5 where the variance of ${\mathcal {D}}_p$ for inertial particle velocity is larger than that for fluid velocity. For these Stokes numbers, high probability is found along the line of ${\mathcal {D}}_p=0$. This indicates that a large number of the particles has quite small divergence, e.g. the divergence of 48.5 % of particles for $St=1$ is smaller than the standard deviation for $St=0$. This means that a group of local particles moves together like an incompressible flow. We can find large variance of the divergence that corresponds to the variance shown in figure 3(b). Such high probability of large positive/negative divergence is observed for cluster cells ($V_p/\overline {V_p} \lesssim 0.5$). A possible reason is that the probability of finding particles in cluster regions is higher than in void regions, as shown in figure 3(a). Another possible explanation is that the ‘caustics’ (Wilkinson & Mehlig Reference Wilkinson and Mehlig2005) of the particle density, where the particle velocity at a single position is multi-valued, result in a large velocity difference of two particles with a small separation distance and cause the extreme divergence values observed in the joint probability distribution functions.

Figure 4.

Joint PDF of the volume of Voronoi cells in log scale and divergence in linear scale for (a) $St=1$ and (b) $St =5$.

One would expect that particles in clusters are more likely to experience negative divergence. However, to maintain the statistically stationary state of the PDF of Voronoi volumes, positive divergence values must be comparable to the negative divergence values in clusters. If cluster volumes showed only negative divergence (cluster formation) and void volumes showed only positive divergence (void formation), particles would concentrate at infinitesimal volumes at the end, i.e. cluster destruction and void destruction are necessary for a statistically dynamic equilibrium.

In figure 4 it is also observed that large positive divergence values are mostly found at smaller volumes compared to large negative divergence values. This trend becomes more pronounced as the Stokes number increases and suggests asymmetry of the probability of positive and negative $\mathcal {D}_p$. The larger probability of positive values of $\mathcal {D}_p$ for small volumes (particles in clusters) is attributed to the crossing motion of particles in clusters. When two particles are approaching each other the divergence is negative, but if two particles have large relative approaching velocities the particles start leaving from each other within the finite time ${\rm \Delta} t$ and the divergence changes to positive. This would happen prominently in caustic regions. As the Stokes number increases, the deviation of particle motion from fluid flow becomes larger and the relative particle velocity increases. Thus, the influence of the crossing motion of particles is more pronounced. These results further suggest that cluster formation is driven by the underlying flow field, but the cluster formation coincides with cluster destruction in a fine cluster with the dissipation scale. Thus, the divergence and convergence of particles in clusters are not determined only by the underlying flow field. To confirm that, we coloured in figure 1 each particle by its divergence value $\mathcal {D}_p$. It is observed that particles with large positive and negative values are concentrated in fine clusters, and, thus, the cluster formation coincides with cluster destruction, as expected. Particles with large positive/negative divergence are found not in all clusters, but some portion of clusters have only moderate divergence. That is, there are two types of clusters, those associated with strong cluster formation and cluster destruction and those which are formed and destroyed by modest change. Figure 1 further illustrates that the clusters with strong formation and destruction are intermittently distributed for all Stokes numbers, even for the case of $St \le 0.2$, and the portion of such clusters seems to increase as the Stokes number increases.

We also compute the mean of the divergence as a function of the volume which corresponds to the conditional average of the divergence, defined as

(4.1)\begin{equation} \langle {\mathcal{D}}_p \rangle_{Vp} = \frac{1}{P(V_p/\overline{V_p})} \int_{-\infty}^{+\infty} {\mathcal{D}}_p P({\mathcal{D}}_p,V_p/\overline{V_p}) \, \textrm{d}{\mathcal{D}}_p , \end{equation}

and shown in figure 5(a). Negative values correspond to convergence of the particles and positive values to divergence. We can see that the mean of the divergence $\langle {\mathcal {D}}_p \rangle _{Vp}$ is negative for large volumes and turns to positive as the volume becomes smaller. This indicates that the cluster formation is active for large volumes, while cluster destruction/void formation is active for small volumes. For the considered Stokes numbers, the zero-crossing points are $V_p/\overline {V_p} \approx 0.1\text {--}0.2$. This can be expressed in terms of the dissipation scale $\eta$ of the flow which is $\eta = 7.9 \times 10^{-3}$. Using the relation $\overline {V_p}/\eta ^{3} = (2{\rm \pi} )^{3}/(N_p \eta ^{3}) \approx 33.54$, we can conclude that the scales of the zero-crossing point are ranging in $V_p^{1/3}/\eta \approx 1.49\text {--}1.88$. Moreover, we observe that when the Stokes number increases, the divergence takes larger values for both negative and positive sides. This means that the cluster formation becomes intensified as the Stokes number increases, and the cluster destruction for small volumes is also amplified. That is, when the Stokes number increases, for $V_p/\overline {V_p} \gtrsim 0.2$, cluster formation is intensified, whereas for $V_p/\overline {V_p} \lesssim 0.2$, cluster destruction is intensified. Strong cluster destruction suggests that clusters are destroyed before they become finer clusters. This is consistent with the fact that fine clusters are less observed as the Stokes number becomes larger than one.

Figure 5.

Mean divergence $\langle {\mathcal {D}}_p \rangle _{Vp}$ (a) and skewness (b) as a function of the Voronoi volume for different Stokes numbers. A moving average filter has been applied for the skewness. The inset shows a zoom to illustrate the presence of negative values.

The difference in magnitude between positive and negative values of the conditional divergence in figure 5(a) can be related to the second-order velocity structure function $S_2(r) = \langle |u(x+r) - u(x)|^{2} \rangle$. For isotropic turbulence, we have $S_2(r) \propto C_2 r^{2/3}$ for $\eta \ll r \ll l _E$ (where $l_E$ is the energy scale) and $S_2(r) \propto r^{2}$ for $r \ll \eta$. The structure function of the inertial particle velocity is similar to the one of the fluid velocity for $\eta \ll r \ll l_E$, for $r \ll \eta$ it deviates from $r^{2}$-scaling and does not approach zero for finite Stokes numbers (Salazar & Collins Reference Salazar and Collins2012; Ireland, Bragg & Collins Reference Ireland, Bragg and Collins2016). This means that for $r \gg \eta$, the r.m.s. particle velocity difference is $\sqrt {\langle |\delta v|^{2} \rangle } \propto r^{1/3}$. Thus, the increase of the velocity difference with $r$ is statistically slower than the increase of the separation distance $r$, i.e. $\mathcal {D}_p \sim \sqrt {\langle |\delta v|^{2} \rangle }/r$ decreases with $r$. For $r \ll \eta$, $\mathcal {D}_p \sim \sqrt {\langle |\delta v|^{2} \rangle }/r$ increases with decreasing $r$ because the r.m.s. particle velocity difference at zero separation distance has a finite value for inertial particles. Thus, the absolute values of divergence decrease as the volume increases.

To quantify the asymmetry of the joint probability distribution functions in figure 4, we plot in figure 5(b) the skewness of the divergence as a function of the Voronoi volume. We observe a similar trend, positive values for small volumes and negative values for large volumes, for the shown Stokes numbers. The negative values of the skewness are more significant than those of $\langle {\mathcal {D}}_p \rangle _{Vp}$, and, thus, we can confirm the asymmetry of the joint PDF at intermediate volumes in the range more clearly, e.g. for $St=1$ from $V_p/\overline {V_p} = 10^{-1}$ to almost $10^{1}$. For all Stokes numbers, $\langle {\mathcal {D}}_p \rangle _{Vp} \approx 0$ and the skewness is close to zero for $V_p/\overline {V_p} > 10^{2}$. This indicates that the joint PDF is symmetric, which suggests that void formation almost balances void destruction.