Cluster formation in particle-laden flows is a continuous phase transition

Numerical simulations for particle-laden 2D Taylor-Green flow are performed for various $St$ in the range $0.01\leq St\leq 3$. The particle dynamics is modeled using the Maxey-Riley equation (see materials and methods for the model).

To obtain the clustering characteristics from the spatial distribution of particles, we consider individual particles as nodes of the network. The adjacency matrix of the network has $N\times N$ elements, where $N$ is the number of particles in the flow domain. Particle positions at a time instant from the simulation form the basis for our network construction. The link between the particles is established based on the proximity determined from the radial distribution function ($g(r)$), which provides a statistical measure on the spatial distribution of particles in the neighbourhood of the particle in consideration  (?, ?) and is calculated using the following equation:

$$g(r)=\frac{N_{r}(r)}{\delta A}\bigg{/}{\frac{N}{A}}.$$

(1)

Here, $N_{r}$ is the number of particles found within a radial distance $r$ from the reference particle, $\delta A$ represents the differential area element of thickness $dr$ constructed at a distance $r$ from the reference particle. $N$ is the total number of particles in the flow domain and $A$ is the area of the flow domain. When particles are uniformly distributed around the reference particle, $g(r)=1$. However, if the particles are clustered in the vicinity of the reference particle, $g(r)>1$ and if the particles are sparsely distributed, $g(r)<1$. We calculate $g(r)$ for all the particles in the flow domain. From the calculated $g(r)$, we find the critical radius $r_{c}$, the first intersection point of $g(r)$ and $g(r)=1$, with $dg/dr\Big{|}_{r=r_{c}}<0$. Links are established between two particles if any one of the particle is within the radius $r_{c}$ of the other particle.

Let us consider two scenarios to illustrate the use of $g(r)$ in our network construction, a particle that is well entrenched in a cluster (Fig. 1a), and an isolated particle where particles are sparely distributed in its vicinity (Fig. 1b). The variation of $g(r)$ for both the particles is shown in Fig. 1c. Here, for the clustered particle, $g(r)$ is greater than one in its vicinity and decreases in the local neighborhood along the radial direction based on the spatial distribution of particles. Particles within the radius $r_{c}$ (see Fig. 1a ) are clustered relative to the reference particle and hence are connected to the reference particle in the network. For the isolated particle (see Fig. 1b), $g(r)<1$ in its vicinity and increases beyond one after a specific radial distance. From the variation of $g(r)$, we can infer that though there are particles in the neighborhood of the reference particle, they are not in the immediate vicinity. Hence, no connection is established with this isolated particle. Following this approach, the existence of a link between any two particles is encoded into a binary adjacency matrix ($A_{ij}$) as 1 if they are connected or 0 otherwise. The links established between the particles are unweighted and undirected. Since we do not consider any self-connections, $A_{ii}$ is set to 0.

Using network measures such as the degree of a node ($k$) and the local clustering coefficient ($C$), we characterize the clustering between particles. The degree of a node ($k$) quantifies the number of connections a node has with other nodes in the network. Higher $k$ indicates the presence of many particles in the vicinity of a node. To remove the dependence on the number of particles in the flow domain, we normalize $k$ with the total number of nodes in our system. The average normalized degree ($\langle k\rangle$), obtained by averaging the degree across all the nodes in the network (Eq. 7), provides a global estimate on the clustering amongst particles. For studying the changes in the network structure, we use degree distribution ($P(k)$ vs. $k$), where $P(k)$ is the probability of a randomly selected node having degree $k$.

The local clustering coefficient ($C$) captures the inter-connectivity between the neighbours of the reference node. $C$ is a local measure useful in capturing the cluster density between the neighbours of the reference node. A particle entrenched in a densely concentrated region of a cluster will have a higher value of $C$ than a particle in a moderately concentrated region. The value of $C$ lies between 0 and 1, where 0 corresponds to an absence of links between neighbors and 1 corresponds to all the neighbors being connected to each other. For obtaining a global estimate on how densely the particles are clustered, we average the clustering coefficient ($\langle C\rangle$) across all the nodes in the network.

From the spatial distribution of particles at each time instant, a proximity-based complex network is constructed and characterized using $\langle k\rangle$ and $\langle C\rangle$. As particles cluster in time, $\langle k\rangle$ and $\langle C\rangle$ of the networks increases. When particles asymptotically settle into the stagnation points, the values of $\langle k\rangle$ and $\langle C\rangle$ saturate at $0.25$ and $1$, respectively. Weakly inertial particles ($St<0.1$) cluster slowly which is evident from the gradual increase of $\langle k\rangle$ and $\langle C\rangle$ (Fig. 2 and movie S1).

As we increase $St$, the rate at which particles are centrifuged out increases  (?); hence, for $St\sim O(1)$, both $\langle k\rangle$ and $\langle C\rangle$ values increase rapidly. For $St\leq 0.7$, all particles in the flow domain asymptotically settle into the stagnation points (movies S1 and S2). Hence, as $t\rightarrow\infty$, both $\langle k\rangle$ and $\langle C\rangle$ saturate at $0.25$ and $1$ respectively. For $St>0.7$, a fraction of particles exhibit periodic motion without settling into the stagnation points  (?) (movie S3). As we increase the $St$, the fraction of particles exhibiting periodic motion increases. Hence, we observe a decrease in the saturated values of $\langle k\rangle$ and $\langle C\rangle$ for $St>0.7$. Furthermore, as a result of the periodic motion of particles, we observe undulations in the saturated $\langle k\rangle$ and $\langle C\rangle$ (Fig. 2). For $St>1$, all the particles exhibit periodic motion without clustering in the shear regions, which is reflected in the lower values of $\langle k\rangle$ and $\langle C\rangle$. Thus, the emergence of clusters for $0.01\leq St\leq 1$ saturates the average network measures at a high value reflecting the global clustering characteristics. These observations are illustrated and analyzed here in detail for the specific cases of $St=0.05$, $St=0.5$ and $St=1$ as shown in Fig. 3- 5.

We start with a random initial distribution of particles in the flow domain (Fig. 3-5(a) at $t=0$). In the corresponding network at $t=0$, most particles have a few connections ($k\approx 0$) while a few nodes have more than just a few connections ($k\approx 0.01$). Corresponding, we observe a power law behaviour in the degree distribution for randomly distributed particles (Fig. 3-5(b) at $t=0$). The inertia of particles comes into play for $t>0$, and the simulations are performed up to $t=100$.

As we progress in time, particles move out of their streamlines and form prominent cluster patterns. The time at which cluster patterns become apparent differ with $St$; for $St=0.05,0.5$ and $1$, patterns that resemble a cross emerge in the shear regions at $t=20,3$ and $2$, respectively (movies S1, S2 and S3). The scatter plots in Fig. 3-5(a) show the increase in the degree of particles that start clustering in the shear region. Correspondingly, we observe a distortion of the original power-law degree distribution in Fig. 3-5(b). As particles form dense clusters, the number of particles with higher $k$ increases significantly.

Thus, the corresponding network possesses a large fraction of nodes with higher degrees. In the degree distribution at $t=100$, the particles with high $k$ correspond to the particles near the stagnation points while particles with low $k$ correspond to the few particles commuting between these stagnation points in the shear regions. Finally, at very large $t$, for $St=0.05$ and $0.5$, all particles settle into the four stagnation points, thus exhibiting a steady state clustering (movies S1 and S2).

The trajectory of a particle initiated in the vortical region is governed by its $St$ and the underlying flow field. Particles with $St\leq 0.25$, do not cross the separatrices in the shear regions due to their weak inertia (?, ?). Instead, the particle gets slowly centrifuged outwards and asymptotically settles into a stagnation point in the edge of the vortex (see Fig.3c ). As this particle starts nearing the stagnation point, the number of particles in its neighbourhood start to increase and hence its $k$ value increases. The oscillations in $k$ arise as the particle approaches and moves farther from stagnation points due its the centrifugal motion. The peaks and troughs in $k$ occur when a particle is nearest and farthest from the stagnation points respectively. Moreover, the frequency of the oscillations in $k$ are directly related to the frequency of the spiral motion of the particle in the vortical region.

On the other hand, a particle with $St>0.25$ initialized inside a vortical region has sufficient inertia to cross the separatrices. Thus, particles originating in different regions can interact with each other  (?, ?) (evident in Fig. 4(c)). Since a particle in the vortical region is quickly centrifuged out and enters the shear region, $k$ of the particle increases sharply. The frequency of oscillation in $k$ is determined by the time the particle takes to move from one stagnation point to the other in the shear region.

Finally, for $St=1$, the particle response time matches with the flow time scale. Particles are ejected very rapidly from the vortex regions and form dense clusters in the shear regions (Fig. 5(a) and movie S3). As a result, $k$ of the particle increases sharply. These particles do not asymptotically settle into the stagnation point; instead, they move in specific periodic orbits in the shear region (Fig. 5(c)). The frequency of oscillation of $k$ is clearly a function of the frequency of particle motion in its periodic orbit. The patterns formed by a collection of such particles remain consistent with time (Fig. 5(a) at $t=100$).

The connections between the nodes in the network increase drastically as particles organize into cluster patterns that resemble a cross. The size of the largest component $(C_{M})$ is defined as the total number of nodes in the largest connected set of nodes in the network  (?). Here, we draw inferences about the process of clustering in particle-laden 2D Taylor-Green flow by studying the variation of $C_{M}$ of the network, which indicates the extent of clustering amongst particles at any given time. The normalized size of the largest component ($C_{M}/N$) changes drastically from $(C_{M}/N)\sim 0$ to $(C_{M}/N)\sim 1$ over a finite time interval ($\Delta t$) characterizing a phase transition. In this study, when $(C_{M}/N)\sim 1$ , the largest component is referred to as the giant component. The emergence of a giant component in the network often occurs via a phase transition  (?). Here, we use $C_{M}$ as the order parameter to characterize this phase transition.

For particles with $St=0.5$, the phase transition in the derived networks occurs around $t=2$ (see Fig. 6a-b) as $C_{M}$ increases rapidly in a short epoch owing to the rapid clustering of particles in the shear regions (see Fig.  4a). The results in Fig. 6 are obtained for an ensemble of $200$ realizations and for a time resolution of 0.02. The size of the giant component saturates and reaches a maximum for a significant epoch between $t=5$ to $10$ (see Fig. 6b). During this epoch, a significant proportion of particles commute in the shear regions between the stagnation points (see Fig. 4a). As particles asymptotically settle into the stagnation points, the number of particles commuting between stagnation points in the shear region reduces leading to the fragmentation of the network and a corresponding decrease in $C_{M}$. Note that, the number of particles that settle asymptotically into each stagnation point is not necessarily equal after complete clustering occurs in the flow field, and thus the size of the largest cluster not necessarily tends to $N/4$ as $t\to\infty$ (see Fig. 6b). The oscillations in $(C_{M}/N)$ occur due to the decrease in connections in the network as the number of particles commuting between stagnation points reduces.

The rate of phase transition in proximity networks is a strong function of $St$ (Fig. 6(c)). Significant variation in the $(C_{M}/N)$ observed in Fig. 6(c) occurs precisely when patterns that resemble a cross become apparent in the flow (Fig. (3-5)a). The time interval of phase transition ($\Delta t$) is determined by setting thresholds for the order parameter ($C_{M}$). We define $\Delta t=t_{1}-t_{0}$, where $t_{0}$ is the maximum time at which $C_{M}\leq N^{\gamma}$, i.e., when $C_{M}$ is a fractional power of the network size. Further, $t_{1}$ is the minimum time at which $C_{M}\geq AN$, i.e., when $C_{M}$ is a significant fraction of the network size. For the results presented in our analysis, we choose $A=4/5$ and $\gamma=3/5$. We observe that $\Delta t$ decreases as we increase $St$ from $0.01$ to $1$ (Fig. 7). Clearly, $\Delta t$ depends on the time taken for the particles to form clusters of discernible patterns in the flow.

A typical feature of phase transitions is the divergence of the variance of fluctuations of the order parameter, called susceptibility ($\chi$), around the critical time when the transition occurs. Susceptibility ($\chi$) is defined as $<C_{M}^{2}>-<C_{M}>^{2}$, where $<>$ is the ensemble average operator. From Fig. 6(d), we observe that susceptibility diverges at faster rate and at lower critical times ($t_{c}$) for higher values of $St$. Here, $t_{c}$ is identified as the time at which $\chi$ attains its peak value in Fig. 6(d). For lower values of $St$, this transition occurs over a longer time interval $\Delta t$ and at higher values of $t_{c}$ when compared to that during higher values of $St$. Shorter epochs of phase transition for particles with higher $St$ indicate the tendency of such particles to cluster rapidly. In Fig. 6(e), we plot the variation of $\chi$ as a function of $(1-t/t_{c})$ on a log-log scale. Firstly, we note in Fig. 6(e) that the critical behavior during phase transition for various $St$ exhibits a power-law scaling around the critical time as $\chi=(1-t/t_{c})^{\alpha}$, with critical exponent $\alpha=-2.8\pm 0.41$ with $90\%$ confidence.

For $St=0$, the particles are essentially tracer particles and do not form any spatial clusters as they follow the fluid streamlines indefinitely. Therefore, there is no phase transition observed for such particles, and thus $\Delta t\to\infty$ as $St\to 0$. For inertial particles with $St\sim O(1)$, the response time of particles is close to the flow timescale and the particles tend to cluster rapidly. Figure 7 shows the variation of $\Delta t$ with $St$ which clearly follows a power-law with a scaling exponent $\beta=-0.95\pm 0.01$ with $90\%$ confidence, initially for $St<0.25$. For particles with $0.25\leq St\leq 1$, the phase transition time varies in the form $\Delta t\sim e^{(C/St)}$ with the constant $C=0.18\pm 0.02$ with $90\%$ confidence. This deviation from power-law behaviour could possibly be due to the formation of caustics amongst the particles in this range $0.25\leq St\leq 1$ and the variation of $\Delta t$ is consistent with the rate of caustics formation  (?, ?).

In order to understand the clustering process using complex networks, we choose the steady 2-D Taylor-Green flow  (?, ?) to be our canonical flow (Fig. 8). This simple flow was chosen for our analysis as it has some of the features of turbulent flows such as coexisting regions of vorticity and shear (?, ?, ?).

The velocity profile equations for Taylor-Green flow are as follows:

$$u_{f}(x,y)=V_{0}\left[-cos\left({\frac{2\pi x}{L_{f}}}\right)sin\left({\frac{2\pi y}{L_{f}}}\right)\right],$$

(2)

$$v_{f}(x,y)=V_{0}\left[sin\left(\frac{2\pi x}{L_{f}}\right)cos\left(\frac{2\pi y}{L_{f}}\right)\right]$$

(3)

In the above equation, $u_{f}$ and $v_{f}$ are the two components of fluid velocity. Here, $L_{f}$ and $V_{0}$ are the length of the flow domain and the velocity of the fluid flow. The time scale, $t_{f}$, of the flow is calculated as $t_{f}=L_{f}/V_{0}$. These aforementioned parameters are used to non-dimensionalize the system of equations in order to minimize the number of input parameters for initiating the simulation. The particles are modeled as point particles having a significantly higher density compared to the surrounding fluid, i.e., $\rho_{p}/\rho_{f}=10^{3}$. Due to the high density ratio, the contributions of added mass effect, buoyancy and Basset history force were neglected (?). Therefore, following the Maxey-Riley equation  (?), only the effect of Stokes drag is considered in the particle dynamics.

We considered $5\times 10^{3}$ particles in our analysis to maintain a sufficiently low volume fraction, to ensure the validity of the unidirectional coupling and to keep the computational time within reasonable limits. In unidirectional coupling, the particle motion is affected by the local fluid structures but the collective motion of the particle does not affect the fluid structures. The dimensionless form of the equations are:

$$\frac{{d}\vec{x}_{p}(x,y)}{{d}t}=\vec{u}_{p},$$

(4)

$$\frac{{d}\vec{u}_{p}(x,y)}{{d}t}=\frac{1}{St}\left(\vec{u_{f}}-\vec{u}_{p}\right).$$

(5)

Here, $\vec{x}_{p}$, $\vec{u_{f}}$ and $\vec{u_{p}}$ are the non-dimensional particle position, flow velocity, and particle velocity vectors, respectively. For obtaining different spatial distributions of particles, we perform the analysis for various $St$ values from $0.01$ to $3$. We limit our study in the range $St$ less than $3$, as the assumption of neglecting the added mass effect and Basset history force becomes invalid for $St$ greater than $3$ (?).

At $t$ = 0, the particles are randomly distributed across the flow domain with initial velocities identical to that of the surrounding flow velocity. Then, the particle inertia is activated. We implement a periodic boundary condition across the edges of the flow domain. The simulations are performed up to $100$ non-dimensional units of time with a time step of $10^{-3}$. The simulation results are consistent with previously reported results from the literature  (?, ?, ?). Using the simulated data, we construct complex networks to study the particle clustering process.

Though $g(r)$ is defined as a continuous function, for estimating it, we need to discretize the space in the vicinity of the particle along the radial direction. Since, we construct our networks from $g(r)$, an optimal discretization resolution ($dr$) should be selected, depending upon the number of particles in the flow domain ($N$) and the size of the flow domain ($L_{f}$). Towards this purpose, we analyzed the variation of the normalized number of components ($N_{comp}/N$) of the network constructed from $N$ randomly distributed particles with respect to the variation of $dr$ for a fixed flow domain size $L_{f}$ (see Fig. 9). We fixed the size of the flow domain as $2\pi$ and varied $N$ from 100 to 50000.

For a specific value of $N$, for a small $dr$ values, $g(r)$ will have sharp peaks at the locations where particles are present and remain zero everywhere else, making it unsuitable for network construction. Hence, the conditions given in Eq. 3 will not be satisfied and the links between particles cannot be established. Thus, the network has a few connections ($N_{comp}/N\approx 1$) (see Fig. 9a). For very large $dr$, $N_{r}/\delta A\approx N/A$, hence $g(r)\approx 1$ becomes incapable of capturing the spatial inhomogeneities, leading to $N_{comp}/N\approx 1$ (see Fig. 9a). Between the two extremes, there are intermediate values of $dr$ where $N_{comp}/N$ reaches a minimum. We choose the $dr$ where $N_{comp}/N$ goes below unity and reaches a minimum in our study. At this resolution $dr$, the $g(r)$ calculated is considered optimal for constructing the network and capturing the spatial inhomogeneities. When the network is constructed using this optimal $dr$, the degree distribution of the network at $t$=0 corresponding to $N$ randomly distributed particles follows a power-law behavior. To verify this, we studied the variation of $\langle k\rangle$ with time for $St=0.5$ particles by varying $N$ from 1000 to 10000 (see Fig. 9b). The variation of $\langle k\rangle$ does not depend to the choice of $N$ and therefore indicative of the robustness of this method in analyzing the clustering characteristics in particle-laden flows.

The degree of a node $k_{i}$ can be expressed in terms of the adjacency matrix as follows,

$$k_{i}=\sum_{j=1}^{N}A_{ij}.$$

(6)

The average degree ($\langle k\rangle$) is calculate as

$$\langle k\rangle=\frac{1}{N}\sum_{i=1}^{N}k_{i}.$$

(7)

The local clustering coefficient ($C_{i}$)

$$C_{i}=\frac{2l_{i}}{k_{i}(k_{i}-1)}$$

(8)

where $l_{i}$ is the total number of links between the $k_{i}$ neighbors of node i. The average clustering ($\langle C\rangle$) is calculate as

$$\langle C\rangle=\frac{1}{N}\sum_{i=1}^{N}C_{i}.$$

(9)