Effect of gravity on clustering patterns and inertial particle attractors in Kinematic Simulation

TL;DR

This study investigates how gravity influences the clustering patterns and attractors of inertial particles in a kinematic simulation, revealing how attractor dimensions and positions vary with particle inertia and gravity effects.

Contribution

It introduces a systematic Lagrangian tracking method to identify and quantify particle clustering attractors under varying gravity and inertia conditions.

Findings

Attractor dimensions vary from 3D to 1D depending on gravity and inertia.

Gravity shifts attractors from horizontal to vertical positions.

Clustering patterns depend on Stokes and Froude numbers.

Abstract

In this paper, we study the clustering of inertial particles using a periodic kinematic simulation. The systematic Lagrangian tracking of particles makes it possible to identify the particles' clustering patterns for different values of particle inertia and drift velocity. The different cases are characterised by different pairs of Stokes number $St$ and Froude number $Fr$. For the present study $0\leq St \leq 1$ and {$0.4 \leq Fr \leq 1.4$}. The main focus is to identify and then quantify the clustering attractor - when it exists - that is the set of points in the physical space where the particles settle when time goes to infinity. Depending on gravity effect and inertia values, the Lagrangian attractor can have different dimensions varying from the initial three-dimensional space to two-dimensional layers and one-dimensional attractors that can be shifted from a horizontal to a vertical position.

Tables

Table 1. Table 1: Periodic KS parameters

$L_x=L_y=L_z$	$1$
$N$	$10$
$N_p$	$15625$
$u_{rms}$	$\mathrm{0\cdot 8703}$
$ℒ$	$\mathrm{0\cdot 2106}$
$L$	$1$
$\eta$	$\mathrm{0\cdot 0642}$
$𝒯$	$\mathrm{0\cdot 2420}$
$t_d$	$\mathrm{1\cdot 1491}$
$k_i^{}{}_{/}{}^{}k_i^{}{}_{min}{}^{}$	$15$
$k_{max}/k_{min}$	$\mathrm{15\cdot 5885}$
$R{e}_L$	$\mathrm{38\cdot 94}$

Table 2. Table 2: Different cases for studying the attractor topology for different ranges of S ​ t 𝑆 𝑡 St .

Case	$St$ range	$Fr$ range	observed patterns
			1D-H	1D-V	2D-L	3D
A	0.041-0.124	0.42-1.34	✓			✓
B	0.165-0.300	0.42-1.34	✓	✓		✓
C	0.331-0.413	0.42-1.34		✓	✓	✓
D	0.496-1.000	0.42-1.34			✓	✓

Table 3. Table 3: Different cases for studying the attractor topology for different ranges of F ​ r 𝐹 𝑟 Fr .

Case	$Fr$	$St$ range	observed patterns
			1D-H	1D-V	2D-L	3D
E	1.01	0-1	✓	✓		✓
F	0.717	0-1		✓	✓	✓
G	0.548	0-1		✓	✓

Equations

$$\frac{d\mathbf{X}}{dt}$$

$$\frac{d\mathbf{V}}{dt}$$

$$\mathbf{u}(\mathbf{x})=\sum _{i=1}^N\sum _{j=1}^N\sum _{K=1}^N{\mathbf{a}}_{\mathbf{ijl}}\cos ({\mathbf{k}}_{\mathbf{ijl}}.\mathbf{x})+{\mathbf{b}}_{\mathbf{ijl}}\sin ({\mathbf{k}}_{\mathbf{ijl}}.\mathbf{x})$$

$$k_i$$

$$k_j$$

$$k_l$$

$${\mathbf{a}}_{\mathbf{ijl}}\cdot {\mathbf{k}}_{\mathbf{ijl}}={\mathbf{b}}_{\mathbf{ijl}}\cdot {\mathbf{k}}_{\mathbf{ijl}}=0$$

$${\left|{\mathbf{a}}_{\mathbf{ijk}}\right|}^2={\left|{\mathbf{b}}_{\mathbf{ijk}}\right|}^2=2E(k)\Delta k_{ijl}/m_k$$

$$E(k)=A{k}^{-5/3}\text{mbox}for{k}_{min}\le k\le k_{max}$$

$$u_{rms}=\sqrt{\frac{2}{3}{\int }_{k_{min}}^{k_{max}}E(k)dk}$$

$$\mathcal{L}=\frac{3\pi }{4}\frac{{\int }_{k_{min}}^{k_{max}}{k}^{-1}E(k)dk}{{\int }_{k_{min}}^{k_{max}}E(k)dk}$$

$$\frac{d\mathbf{V}}{dt}=\frac{1}{{\tau }_a}\left(\mathbf{u}({\mathbf{x}}_p(t),t)-\mathbf{V}(t)+{\mathbf{V}}_d\right)$$

$$St={\tau }_a/\mathcal{T}={\tau }_a{u}_{rms}/\mathcal{L}$$

$$Fr=u_{rms}/\sqrt{g\mathcal{L}}$$

$$\gamma =V_d/u_{rms}$$

$$\gamma ={\tau }_{a}g/u_{rms}$$

$$Fr$$

$$\gamma$$

$$N_r\sim r^D$$

$$D=\frac{\Delta \ln N_r}{\Delta \ln r}$$

$${\Delta }_{mi}^2=(x_m-x_i{)}^2+(y_m-y_i{)}^2+(z_m-z_i{)}^2$$

$$\Delta =\frac{1}{N_p}\sqrt{\sum _{m=1}^{N_p}{\Delta }_{mn}^2}$$

$${\Delta }_H=\frac{1}{N_p}\sqrt{\sum _{m=1}^{N_p}(x_m-x_n{)}^2+(y_m-y_n{)}^2}$$

$${\Delta }_V=\frac{1}{N_p}\sqrt{\sum _{m=1}^{N_p}(z_m-z_n{)}^2}$$

Full text

Kinematic Simulation, Particle dispersion, inertia Particle, Multi-particle sets

††preprint: APS/123-QED

I Introduction:

Clustering could be defined as the propensity of an initially uniformly distributed cloud of particles to accumulate in some regions of physical space. This is an important phenomenon to understand in order to explore, identify and possibly monitor some natural or hand-made mixing processes such as those causing rain formation Falkovich et al. (2002), sediment transportation Pan et al. (2011), fuel mixing and combustion.

There are different ways to analyse particle clustering in turbulent flow and Direct Numerical Simulation (DNS) is the most widely used method (e.g. Cencini et al. (2006); Saw et al. (2006); Falkovich-Pumir (2004)). Particle clustering depends on both the flow conditions and the particle characteristics. Different flow conditions can lead to different clusters. The clustering mechanism would be different in the inertial or dissipation range of turbulent flow Bec et al. (2007). In our paper we only study the effect of the scales in the inertial range and this is possible by using a synthetic model where forcing and dissipation are not needed to develop the inertial range. While considering particle characteristics, most of the studies on particle clustering have been conducted in the absence of external forces on particles but the effect of gravity(external force) was discussed in relation to cloud physics and rain formation in Falkovich-Pumir (2004); Woittiez et al. (2008).

More recently, the effect of gravity on clustering mechanism has been further emphasized in Yongnam Park (2014); Gustavsson et al. (2014); Bec et al. (2014). In the present study, to observe the clustering pattern in the presence of gravity, the particles are initially uniformly distributed in the Kinematic Simulation (KS) flow. Though there is no particular difficulty in considering particles with different inertia in Kinematic Simulation, this study is limited to mono-dispersed seeding i.e. particles having the same inertia. Furthermore, the particles are considered small enough so that they neither affect the flow nor interact with each other (one-way coupling). The positions of particles are monitored as a function of time and a Lagrangian attractor is observed for some cases. That is, the initially distributed cloud of particles will end in a set of loci that does not evolve any further. The particles move within that set of loci which defines the structure of the Lagrangian attractor and its dependence on $St$ and $Fr$ numbers is studied here.

We only consider attractors with integer dimensions (one-dimensional and two-dimensional structures) which are easy to identify. Different types of methods can be found in the literature to identify and then quantify particle clustering patterns, e.g.: correlation dimension Bec et al. (2007), radial distribution function RDF Saw et al. (2006), average-distance-to-nearest-neighbour method Yongnam Park (2014). The selection of a method is mainly based on the objective of the study. For example, the RDF has the advantage of being directly related to the droplet collision rate. For the present work, the Box counting method BCM and the nearest-neighbour distance analysis are implemented to identify the integer dimensions of Lagrangian attractor in the presence of gravity.

The paper is organised as follows: in § II we introduce the KS model, its notations and its parameters. The different kinds of Lagrangian attractor are discussed and introduced in § III. A quantitative analysis is conducted in § IV. Section V summarises our main conclusions.

II Kinematic simulation technique:

Kinematic Simulation (KS) is a particular case of synthetic turbulence where the focus is on particle’s trajectory at the expense of solving the Navier-Stokes equation. An analytical formula ‘synthetic flow’ is used for the Eulerian flow field. The simplicity of the KS model excludes some features of real turbulent flow but capture the part of the physics which is required to perform Lagrangian particle tracking. Such is the idea with synthetic turbulence which retains less information than the whole flow, but tries to keep what is paramount for the Lagrangian story.

KS modelling has been successfully employed and validated Fung et al. (1992); Elliott and Majda (1996); Malik and Vassilicos (1999). This kind of simulation is much less computing-time consuming than DNS which is important for the present study where we need to run many cases (about 400 cases for up to 1200 turnover time). Each case corresponds to a given $St$, $Fr$ and time and involves 15625 particles.

With synthetic simulations, one can develop models where turbulence ingredients and complexity can be added step by step helping to understand their respective importance. These synthetic models can be a useful complement to Direct Numerical Simulation. In particular, KS was instrumental in discriminating between the role of Lagrangian and Eulerian correlations for vertical diffusion in stratified and rotating flows Nicolleau and Yu (2007). With KS it is also possible to play with the spectral law Nicolleau and Nowakowski (2011) and its consequences in terms of particle’s dispersion. We also refer to the work of Malik (2014a, b) for a discussion on how the work on KS can help to understand the sweeping effect on two-particle dispersion.

KS was first introduced as a way to understand particle dispersion rather than particles clustering but we propose here a work getting back to the main strength of KS. That is to provide a coherent Lagrangian framework where some parameters (e.g. spectra Nicolleau and Nowakowski (2011), waves Nicolleau and Yu (2007); Nicolleau et al. (2013), …) can be studied in details posing the basis for a comparison with experiments. Previous work Ijzermans et al. (2010); Meneguz and Reeks (2011) particularly supports the use of KS for studying the evolution of the particle cloud in the absence of gravity effect which made the study more about segregation than clustering.

As we are not interested in two-particle dispersion, we limit our study to small Reynolds numbers, more precisely to scale ratio ${k_{i}}_{max}/{k_{i}}_{min}=15$111$i=1$, 2 or 3.

II.1 Periodic KS method for isotropic turbulence

In Kinematic Simulation the underlying Eulerian velocity field is generated as a sum of random incompressible Fourier modes with a prescribed energy spectrum $E(k)$. Here we limit the study to a Kolmogorov type spectral law $E(k)\sim k^{-5/3}$. Using KS, the computational task reduces to calculate the trajectory of each particle placed in the turbulent field initially at $\mathbf{X}_{0}$. Each trajectory is, for a given initial condition, solution of the differential equation set:

[TABLE]

where $\mathbf{X}(t)$ is the particle’s position, $\mathbf{V}(t)$ its Lagrangian velocity and ${\bf u}_{E}$ the analytical Eulerian velocity used in KS. $\mathfrak{F}$ is a function relating the Lagrangian acceleration to the Eulerian and Lagrangian velocities.

In KS ${\bf u}_{E}$ takes the form of a truncated Fourier series, sum of $N_{k}=N^{3}$ Fourier modes:

[TABLE]

where $\mathbf{a_{ijl}}$ and $\mathbf{b_{ijl}}$ are the decomposition coefficients corresponding to the wavevector $\mathbf{k_{ijl}}$. In its general form the KS field can also be a function of time but we limit the study to a steady KS. The effect of introducing a time-dependence in the Fourier modes will be the objective of future study.

Unlike the classic KS decomposition (Fung and Vassilicos, 1998; Nicolleau and ElMaihy, 2006), here the wavevectors $\mathbf{k_{ijl}}=(k_{i},k_{j},k_{l})$ are implemented arithmetically to enforce a periodic condition for the flow field:

[TABLE]

where $(n_{i},n_{j},n_{k})$ are integers satisfying $1\leq n_{i}\leq N$. In practice, we choose $(L_{x}=L_{y}=L_{z})$ for creating isotropic turbulence and to ensure the flow incompressibility the Fourier coefficient vectors $\mathbf{a_{ijl}}$ and $\mathbf{b_{ijl}}$ are set orthogonal to the wavevector:

[TABLE]

Their magnitude is fixed by the energy spectrum, $E(k)$

[TABLE]

where $m_{k}$ is the number of wavevectors of wavenumber $k=\|\mathbf{k_{ijl}}\|$. The spectrum follows the Kolmogorov form in the inertial range,

[TABLE]

where $A$ is a constant. From the spectral law, the rms velocity and the integral length scale can be defined as follows:

[TABLE]

[TABLE]

The Kolmogorov length scale is defined as $\eta=2\pi/k_{max}$, whereas the largest physical scale is $L=2\pi/k_{min}$ which determines the inertial range $[\eta,L]$ over which (9) is observed. It is worth noting that $\mathcal{L}\simeq L$ for sufficiently large inertial ranges. However, here in contrast to other KS studies the inertial range is small and $L\simeq 5\mathcal{L}$. In this paper, non-dimensional numbers ($St$ and $Fr$) are based on the integral length scale $\mathcal{L}$ and for the sake of future comparisons both are reported in Table 1. The ratio between the largest length scale and the Kolmogorov length scale is $k_{max}/k_{min}$ and the associated Reynolds number is: $Re_{{L}}=(k_{max}/k_{min})^{4/3}$. This is the standard way to define a Reynolds number in KS and a DNS or an experiment yielding the same ratio $k_{max}/k_{min}$ would have a much larger Reynolds number. Finally, a characteristic time for normalisation can be $t_{d}={L}/u_{rms}$ or $\mathcal{T}=\mathcal{L}/u_{rms}$. All the periodic KS parameters are gathered in Table 1.

The particles are initially homogeneously distributed as shown in Fig. 1a and whenever a particle leaves the turbulence box domain (e.g. $\textbf{X}_{i}>L_{x}$), then it is re-injected from the opposite side as shown in Fig. 1b to keep the periodic condition.

II.2 Equation of motion

Following Abou-El-Azm and Nicolleau (2008) the equation of motion for the inertial particle is derived from Gatignol (1983); Maxey and Riley (1983) and consists of a drag force and drift acceleration (weight):

[TABLE]

where $\tau_{a}$ is the particle’s aerodynamic response time and $V_{d}=\tau_{a}\;{\bf g}$ the particle’s terminal fall velocity or drift velocity.

II.3 Non-dimensional parameters

Three non-dimensional parameters are introduced to make qualitative and quantitative analyses of the particle clustering.

• •

The Stokes number expresses the ratio between the particle’s response time (inertia effect) and the turbulence characteristic time

[TABLE]

It measures the relative importance of the particle inertia. In the limiting case $St=0$; the heavy particles recover the motion of the fluid tracers, whereas for $St\rightarrow\infty$ the heavy particles become less and less influenced by the surrounding velocity field.

• •

The Froude number is the ratio between inertial forces and gravitational forces.

[TABLE]

In our study the rms velocity $u_{rms}$ and inertial length scale are constant and $g$ is varied.

• •

The Drift parameter is the ratio between the particle’s drift velocity and the turbulence rms velocity:

[TABLE]

The Drift parameter can still be defined without gravity. Then $\gamma$ can be considered as measuring the effect of a mean velocity $V_{d}$.

If $V_{d}$ is caused by gravity:

[TABLE]

in this case the drift parameter is affected by both the gravity and the particle’s inertia.

• •

$\gamma$ can be expressed as a function of Stokes and Froude numbers so for a given turbulence the case corresponding to a constant gravity, that is varying $\tau_{a}$ only, is given by

[TABLE]

III Results and Discussion

III.1 Clustering Pattern variations in relation to time of evolution

The particles initially uniformly distributed in the flow field are allowed to evolve until a Lagrangian attractor is achieved. The shape of the attractor varies from clear one-dimensional structures (Figs 2 and 3) to three-dimension distributed structures (Fig. 4de) or two-dimensional curtain-like structures (Fig. 4dh). For a short time the attractor’s shape is time-dependent as shown in Figs 2 and 3. The time evolution of the cluster depends on non-dimensional parameters $St$ and $Fr$ as illustrated in Fig. 2 where it takes 4 times longer to reach the one-dimensional Lagrangian attractor than in the case of Fig. 3. In this paper, we do not intend to study the temporal evolution of cluster attractors and are only interested in attractor’s asymptotic form (i.e. for $t\to\infty$).

In this section, we focus on the qualitative measure of attractors and only a few cases are presented in Fig.4 and Fig5. A systematic quantification will be proposed in § IV which consists of a comprehensive set of data generated with small increments in $Fr$ and $St$ numbers. As the one-dimensional Lagrangian attractor is observed for various pairs of $St$ and $Fr$, it was difficult to suggest any definite relationship between one-dimensional clustering and $St$ and $Fr$ . Hence, a scheme of further classification is adopted to establish such a relationship. We use the following nomenclature:

• i)

1D-H : horizontal one-dimensional Lagrangian attractor as in Fig. 3bf,

• ii)

1D-V : vertical one-dimensional Lagrangian attractor as in Fig. 4bg,

• iii)

1D-HV : Intermediate one-dimensional Lagrangian attractor as in Fig. 4ag,cg,

• iv)

2D-L : two-dimensional vertical curtain-like layer as in Fig. 4dh (see also Woittiez et al. (2008)),

• v)

3D : any three-dimensional structure without any particular structure in the cloud as in Fig. 4de.

The qualitative results are split into three different categories which can take into account the effect of gravity and/or inertia:

• i)

keeping $St$ constant §III.2,

• ii)

keeping $Fr$ constant §III.3,

• iii)

keeping $\gamma$ constant §III.4.

III.2 Clustering Pattern variations in relation to Constant Stokes Number $St$

We can analyse the results by fixing the Stokes number $St$ and varying the Froude number $Fr$. Four representative cases with varying values of $Fr$ are shown in Fig. 4, in columns a, b, c and d, namely $St=0.165$, 0.249, 0.331 and 0.663. For each of the Stokes numbers, we explored within the range $0.548\leq Fr\leq 1.34$.

As the Froude number decreases, the particles first cluster on a one-dimensional Lagrangian attractor then they may rearrange on another 1D or 2D Lagrangian attractor (Fig. 4 bh, ch, dg and dh). The Lagrangian attractor also moves from a predominantly horizontal (Fig. 4 ae, af and bf) to a vertical direction (Fig. 4 ag, bg and cf) with decreasing value of $Fr$.

The qualitative shapes with varying values of $Fr$ are listed in Table 2. When one-dimensional structures are in-between vertical and horizontal (1D-HV) directions as in Fig. 4bf we ticked both 1D-H and 1D-V in the table.

Now we describe each case one by one according to observed patterns. For Case A, corresponding to low values of $St$, the decrease in $Fr$ (increase in gravity) forces the particles to cluster in a horizontal direction whereas the particles will disperse evenly as $Fr$ is further decreased. The increases in gravity destroy the particles’ clustering for low range of $St$. Similarly, for case B, initially distributed inertial particles cluster into a horizontal attractor (1D-H) Fig. 4bf and then an increase in the gravity effect (lower $Fr$) causes the particles to accumulate in the vertical direction (1D-V) as in Fig. 4bg. With further increases in the gravity they recover a 3D distribution as shown in Fig. 4bh.The appearance of 1D-V in case B shows the greater effectiveness of gravity at a relatively high range of the $St$ for a given value of $Fr$.

For further increases in the Stokes number $St$, case C, the horizontal structure (1D-H) is not observed. Rather a vertical 1D pattern (1D-V) is seen as in Fig. 4cf which can transform into a 1D-HV attractor as in Fig. 4cg with decreasing values of $Fr$. This implies that the particle inertia starts dominating over the flow Eulerian structure and allows the gravity to play a more important role. At higher values of $St$ (case D), there is no more one-dimensional clustering but some clustering can still be observed in the form of two-dimensional vertical curtain-like structures as shown in Fig. 4dh at low values of $Fr$.

III.3 Clustering Pattern variations in relation to Constant Froude Number $Fr$

A constant Froude number corresponds to the case of varying the particle’s property ($\tau_{a}$) for a given environment (turbulence and gravity) which exists in most of the experimental situations. The variations in clustering patterns are identified by keeping $Fr$ constant while varying $St$. For the purpose of qualitative measures, three different cases as shown in Table 3 are considered with small increments in the $St$ ranging [0-1].

In Fig. 4 cases of constant Fr correspond to the horizontal rows e, f, g, h. As $St$ increases, the particles’ one-dimensional clustering is first enhanced and then destroyed to eventually reappear in the form of a two-dimensional layer (2D-L).

For high values of $Fr$ (low gravity), case E corresponding to Fig. 4f, particles settled on horizontal one-dimensional structures (1D-H) for low values of $St$. The increase in $St$ values resulted into vertical one-dimensional structures (1D-V). For the mid-range values of $Fr$, case F corresponding to Fig. 4g, the clear one-dimensional horizontal structure (1D-H) is no longer observed but instead some intermediate (1D-HV) one-dimensional structures can be seen for low $St$ values which converge into a layered curtain-like (2D-L) structure as $St$ is increased. Finally, low values of $Fr$ (case G, Fig. 4h) allow the particles to accumulate predominantly in the direction of gravity, so vertical patterns are identified such as 1D-V and 2D-L structures.

III.4 Clustering Pattern variations in relation to constant drift parameter $\gamma$

It results from the previous discussion that the variations in inertial and gravity effects do not have a monotonic effect on the particle clustering. Physically gravity and inertia are combined effects but one can consider a particle subjected to a drift velocity without referring explicitly to gravity. This effect of drift can be assessed by identifying the patterns with the drift parameter $\gamma$ instead of $Fr$. So here we want to observe the variation in the particle attractor by keeping the drift parameter $\gamma$ constant. Fig. 5 shows the three cases $\gamma=0.135$, 0.689 and 2, as $St$ increases for the arbitrary time $t=300$. It is not possible to keep the same range for $St$ for all cases because $St$ and $\gamma$ are linked.

For low values of $0\leq\gamma\leq 0.2$ (Fig. 5 top), particles tend to accumulate as a horizontal attractor for low values of $St$. As observed previously, further increases in $St$ lead to particles scattering. Eventually particles disperse evenly in the flow for high values of $St$.

As $\gamma$ increases $0.2\leq\gamma\leq 0.8$ (Fig. 5 middle), horizontal attractors are not observed any more for low values of $St$, instead some intermediate 1D-HV and vertical 1D-V attractors are observed with increasing values of $St$. Further increases in inertia disperse the particles evenly in the flow field. Third case $0.8<\gamma<2$ (Fig. 5 bottom) corresponds to relatively high values of $\gamma$. The particles are trapped in a 2D-L structure. The increase in $St$ results into the particles dispersing more homogeneously on this 2D-L attractor.

To summarise: an increase in $\gamma$ can lead to a 1D-V or 2D-L rather than a 1D-H attractor and an increase in $St$ destroys the one dimensional attractor leading to the particles reorganising on a 2D-L.

This is in agreement with Yongnam Park (2014) who mention gravity-driven clustering of inertial particles in turbulence and report a different kind of particle clustering caused purely by gravity, that is, clustering in a vertical stripe pattern formed when strong gravity acts on heavy particles.

IV Quantification of Clustering Patterns

Visualizations of the particle cloud for small discrete increments of the non dimensional parameters $St$ and $Fr$ can be tedious. It means looking at about 400 cases in this study in a systematic order. Beyond the simple visualisation, it is important to quantify the Lagrangian attractors using an appropriate method for spatial clustering. Two different methods are considers here: the box-counting method (BCM) and the average-distance-to-nearest-neighbour $\Delta$. The average-distance-to-nearest-neighbour is eventually chosen for the final quantitative analysis.

IV.1 Box Counting Method

The box counting method (BCM) is a commonly used method to determine the fractal dimension of an object. Though in our simulation the range of scales is too short to observe the fractal patterns described in Nicolleau and ElMaihy (2004), BCM remains a useful tool to discriminate between one-dimension, two-dimension and three-dimension clustering patterns. The fractal dimension, $D$, represents the relation between the box size, $r$ and $N_{r}$ the number of boxes needed to cover the cloud of particles, that is:

[TABLE]

It is straightforward to obtain the fractal dimension from a log-log plot:

[TABLE]

A validation of the method is made on three clear identified shapes, namely the one-dimensional Lagrangian attractor, the two-dimensional curtain-like layered pattern and a three dimensional distribution. As shown in Fig. 6, the difference between these three patterns is clearly captured.

The BCM is sensitive to the achievement of the attractor, that is if few particles have not settled on the attractor they can alter the box counting results. So, with this method it is necessary to make sure that the cloud shape is the asymptotic final one which requires a very long time. As illustrated in Fig. 2: though at a very short time $t=10$ the position and shape of the 1D-V is obvious, it is necessary to wait up to $t=1200$ to get the final cluster position that will allow a correct measure for the BCM. All the cases are reported in Fig. 7 for clusters having reached their final shapes (attractors). Iso-contours of $D$ as a function of $(St,Fr)$ are plotted in Fig. 7 and it appears that for $St>0.45$ strong one-dimensional clustering has vanished.

A major problem with the box counting method is to discriminate between the dimension $D$ of very similar patterns as shown in Fig. 8.Therefore, we cannot be sure that the box-counting result is meaningful in a case with no clear structures (i.e. without an integer dimension. As mentioned earlier, in order to be accurate, the BCM must be applied to the Lagrangian attractor. If the particle cloud has not settled on the attractor as at time $t=10$ or $t=600$ in Fig. 2, the BCM will not educe the 1D patterns. So, in practice, it means running the cases for long times until the particles have all settled on the Lagrangian attractor which can be prohibitive.

IV.2 Nearest-neighbour distance analysis

The advantage of using this approach is that it is not necessary to reach the final cluster at $t\to\infty$, a snapshot at earlier times gives us a clear idea of the kind of Lagrangian attractor to expect. For example in Fig. 2 the kind of 1D-V Lagrangian attractor is clear at $t=10$ there is no need to wait for the final asymptotic shape at $t=1200$. Though other methods may struggle to pick up the structure at $t=10$ and will only give the correct diagnosis when all the points have settle on the attractor, that is for $t=1200$, the nearest-neighbour distance analysis will work for intermediary times. The average-distance-to-the-nearest-neighbour $\Delta$ Yongnam Park (2014) is introduced to systematically quantify the clustering patterns. At a given time for each particle $\mathbf{X}_{m}$ its nearest neighbour $\mathbf{X}_{n}$ is such that

[TABLE]

is a minimum for $i=n$. Then we define the average-distance-to-the-nearest-neighbour as

[TABLE]

Where $\mathbf{X}_{n}=(x_{n},y_{n},z_{n})$ is the nearest particle’s neighbour.

We get three obvious benchmark values for this method.

• i)

If the particles are homogeneously distributed as at time $t=0$, (Fig. 1a) then $\Delta\simeq L_{x}/N=1/25=0.04$.

• ii)

If the particles are distributed on a surface-like attractor 2D-L as in Fig. 5) then $\Delta\simeq L_{x}/N^{3/2}=1/25^{3/2}=0.008$.

• iii)

If the particles are distributed on a line-like attractor as in Fig. 2) then $\Delta\simeq L_{x}/N^{3}=1/25^{3}=6.4\,10^{-5}$.

In practice, the method will detect a one-dimensional structure (1D-H or 1-DV) for $\Delta\leq 0.008$ while 2D layered structures are observed for $0.01\leq\Delta\leq 0.014$.

We applied the average-distance-to-nearest-neighbour method to all run cases to see the variations in the attractor patterns for the same time $t=300$. Colour-wise blue (darkest spots for $St\leq 0.4$) corresponds to the 1D Lagrangian attractor, yellow-green (light grey) to the 2D-L and dark red (darkest spots for $St\geq 0.4$) to 3D structures. Iso-contours of $\Delta$ are plotted as functions of ($St$, $Fr$) in Fig. 9a and ($St$, $\gamma$) in Fig. 9b to see the effect of varying gravity and inertia on the Lagrangian attractors. Fig. 9 confirms that the clustering of inertial particles is not a monotonic function of either $St$ or $Fr$ number. However, it is possible to identify regions in the plane $(St,Fr)$:

• i)

In agreement with Fig. 7, there is no one-dimensional structures for $St\geq 0.5$.

• ii)

We can refine the analysis in term of $\gamma$: In Fig. 9b, for large values of the Stokes number ($St>0.3$) and gravity effects ($\gamma>0.8$) the 2D-L structures are predominant. This is also in agreement with Gustavsson et al. (2014) whose calculations show that for large values of $St$ particles may cluster strongly.

Another advantage of the average-distance-to-nearest-neighbour method is that variations in horizontal $\Delta_{H}$ and vertical $\Delta_{V}$ directions can be identified separately which help to monitor anisotropic patterns. In practice, $\Delta_{H}$ and $\Delta_{V}$ are defined as follows:

[TABLE]

[TABLE]

Fig. 10a shows the iso-contours of $\Delta$ as a function of ($St$, $Fr$) in the region where one-dimensional Lagrangian attractors are observed, i.e. $0.05\leq St\leq 0.4$. The two different types of one-dimensional attractor either horizontal (1D-H) or vertical (1D-V) can be further analysed in Figs 10b and c. Fig. 10c shows the ratio $\Delta_{H}/\Delta$ where one-dimensional attractors exist, that is when $\Delta\leq 0.008$. $\Delta_{H}/\Delta\leq 0.5$, that is blueish points (light grey area), indicates 1D-V structures; whereas $\Delta_{H}/\Delta\geq 0.75$, that is redish points (dark spots for $Fr\geq 0.7$), indicates 1D-H structures. Fig. 10b describes a similar relationship based on $\Delta_{V}$. So it is clear from the points’ colours distribution that horizontal attractors are predominant for large $Fr$ while vertical attractors are prevalent as $Fr$ decreased.

V Conclusion

We used Kinematic Simulation (KS) to study the clustering pattern of particles with inertia subjected to gravity effects. For some combined inertia and gravity effects ($St$, $Fr$), the particles cluster on a fixed space subset. That subset can be one-dimensional or two-dimensional. In most of the cases, the particles did not cluster and disperse occupying the most of the periodic box.

Using KS, it became possible to investigate many combinations of ($St$, $Fr$) and educe and classify those one-dimensional or two-dimensional subsets. Though KS retains only part of the turbulence physics, it helps to understand the clustering patterns and the effect of gravity on these patterns.

The main results can be summarised as follows (and in a more synthetic presentation in Figure 11):

• •

The effect of gravity may reduce or enhance inertial particles clustering (as noticed in Gustavsson et al. (2014); Bec et al. (2014)) depending on the Stokes number. This effect can lead to strongly anisotropic clusterings (1D or 2D-L) very clearly evidenced by the KS model.

• •

The 1D structure is better observed with the synthetic flow, as in real flows unsteadiness may prevent the particles from reaching that asymptotic state. These 1D attractors move from the horizontal to the vertical direction as the $Fr$ number decreases.

• •

For our range of Froude numbers, we found two critical Stokes numbers: for $St>St_{cr1}=0.3$ there is never occurrence of a horizontal (1D-H) type attractor and no 1D-type attractor is found for $St>St_{cr2}=0.5$.

• •

For low values of $Fr$, curtain-like two-dimensional layered structures similar to the ‘curtain-like manifolds’. already observed in Woittiez et al. (2008) are recovered as the high gravity prevents the inertial particles from settling uniformly in the turbulent flow.

VI Acknowledgement

M. Farhan gratefully acknowledges support funding from the Department of Mechanical Engineering of the University of Engineering and Technology Lahore, Pakistan.

This work was supported by EPSRC grant EP/L000261/1