3.1. Fractal dimension analysis

To scrutinise the aforementioned hypothesis, we adopt the fractal dimension analysis, as particle dispersion is multi-fractal in nature (Bec Reference Bec2005). Here, particle clustering shows distinct fractal features at different turbulent scales, or, in other words, particles with different inertia primarily respond to eddies of different sizes (Eaton & Fessler Reference Eaton and Fessler1994). Thus, we perform the fractal dimension analysis on the global particle dispersion, as well as on the individual particle clusters. First, we compute the correlation dimension ( $D_2$ ) of each particle species, which characterises the dispersion pattern (chaotic or organised) of the particles (Grassberger & Procaccia Reference Grassberger and Procaccia1983). The parameter $D_2$ is expressed as

(3.4) \begin{equation} D_2 = \lim _{l \to 0} \frac {1}{\log (l)} \log (C_l) ,\end{equation}

where $C_l = \sum p_i^2$ is the correlation integral, and $p_i$ is the probability of finding a particle pair at an inter-particle distance less than $l$ . By plotting the logarithmic correlation integral ( $\log (C_l)$ ) against the log of the normalised separation distance, $D_2$ can be approximated as the slope of the obtained curve. This is shown in figure 5 (a), and the corresponding $D_2$ value of each particle species is shown in the legend, marked with the particle species number as a subscript. The fundamental benefit of using $D_2$ to characterise the dispersion is that it is not dependent on the choice of the length scale. On the other hand, methods such as box counting depend on the box size (Borgani & Murante Reference Borgani and Murante1994) , and pair-correlation function is influenced by the radius (Monchaux, Bourgoin & Cartellier Reference Monchaux, Bourgoin and Cartellier2012) in which the number of particle pairs is counted.

For $D_2$ , a lower value depicts a more organised (clustered) dispersion of the particles; a higher value shows chaotic dispersion. As $D_2$ decreases with the increase in particle size in figure 5, it validates the qualitative observations of figure 4, where smaller particles are more dispersed as they start spreading laterally after a short axial distance from the inlet, while larger particles primarily travel along the jet axis. Thus, $\mathrm{P4}$ particles depict organised (clustered) dispersion.

By analysing $D_2$ of the individual particle cluster, we can probe if the fractal dimension changes based on the cluster size. If it does, it would indicate that the cluster size varies based on the turbulent scale (since particle clustering is governed by particles’ dominant correspondence to a turbulent scale), which would hint at the mechanism that causes particles to spread laterally at different downstream distances. Note that $D_2$ indicates the scale of clustering, where a lower $D_2$ depicts clustering due to the larger scales, while a higher $D_2$ relates clustering to the small-scale features (Tang et al. Reference Tang, Wen, Yang, Crowe, Chung and Troutt1992). Considering the visual observation in figure 4, the scale at which clustering occurs should increase with the increase in particle size. We witness this trend in figure 5 (b), where $D_2$ is plotted against the cluster volume ( $V_c$ ). The dispersion dimensionality of these curves also verifies the qualitative observation of figure 4. Here, $V_c$ is determined by using the density-based spatial clustering of application with noise (DBSCAN) method (Ester et al. Reference Ester1996), which segregates clustered particles from noise/unclustered particles. By computing the convex hull using the $x$ , $y$ and $z$ separations between particles in each cluster, $V_c$ is evaluated. The details of the DBSCAN method in the present context can be found in Saieed & Hickey (Reference Saieed and Hickey2024).

The downward curves in figure 5 (b) portray a critical observation: the $D_2$ rises quasi-linearly until it reaches an inflexion point. This inflexion point marks the boundary between chaotic (rising $D_2$ with respect to $V_c$ ) and organised (declining $D_2$ with respect to $V_c$ ) dispersion regimes; these inflexion points are approximately highlighted with a grey colour region. This trend, however, seems counterintuitive because the eddy size increases with the axial distance. But the rise in $D_2$ with $V_c$ on the left side of the grey region containing inflexion points means that particle clustering is governed by turbulence at smaller scales. This can be understood with the fact that, in the present $10 d_{\!j}$ domain length, which is essentially the near-field region of the jet, the turbulent fluctuations and thus the averaged $\mathrm{TKE}$ augments as the local eddy length scale ( $L_{\textit{eddy}}$ ) increases downstream, see figure 6. The rise in downstream turbulent fluctuations at the jet centreline can be witnessed in figure 3 (b). Additionally, the $\mathrm{TKE}$ in figure 6 is the bulk $\mathrm{TKE}$ of the flow, which is computed by averaging $\mathrm{TKE}$ in the radial direction for each grid point along the $x$ -axis. Therefore, because we are in the near-field region where shear is high at the inlet (shear layer is significant) and the jet is spreading, the bulk $\mathrm{TKE}$ increases as the velocity fluctuations are high in this region ( $x=10$ ).

Figure 6. Evolution of temporally and spatially averaged $\mathrm{TKE}$ and eddy length scale ( $L_{\textit{eddy}}$ ). Here, $L_{\textit{eddy}}$ is normalised with jet diameter $d_{\!j}$ . The dotted vertical lines mark the axial location ( $x_r$ ) of the respective particle initial radial dispersion, where ${\textit{St}}^*$ is the local Stokes number of the particles at that location (discussed in § 3.3).

In figure 6, the $L_{\textit{eddy}}$ is computed by identifying vortices with the $\mathcal{Q}$ -criterion and labelling dominant vortices on the $y$ - $z$ plane at each $x$ position by applying an approximate threshold ( $\mathcal{Q} \gt \overline {\mathcal{Q}}$ ). Each detected vortex has a peak value of $\mathcal{Q}$ . By identifying the location of this peak $\mathcal{Q}$ , and determining the average Euclidean distance between the locations of these peak $\mathcal{Q}$ values of vortex pairs on each $y$ - $z$ plane, the local separation of vortices is estimated. This local separation is temporally averaged for each particle pair at a given $x$ -position, which resulted in $L_{\textit{eddy}}$ . This provides a more accurate estimate of the local eddy length scale as compared with $r_{1/2}$ , which is often used. However, note that, although our later results converge based on this definition of $L_{\textit{eddy}}$ , the values of $L_{\textit{eddy}}$ are sensitive to the applied threshold. Also, note that $L_{\textit{eddy}}$ is qualitatively similar to $r_{1/2}$ of the jet; the change in threshold is expected to yield similar trends. This is because selecting a different threshold will alter the magnitude of $L_{\textit{eddy}}$ , but $L_{\textit{eddy}}$ will still increase axially, akin to the trend in figure 6; the formal assessment of this sensitivity is left for future consideration.

The discussion above infers that because particles are accelerated in the axial direction by the supersonic gas in the jet core, the natural reaction of the particles is to stay within the central region of the jet, as drag (governed by the fluid–particle slip velocity) is higher within this region. This limits $V_c$ at any axial location. Since the jet spreads downstream, particles get more space to form clusters, which raises $V_c$ . However, the rise in $V_c$ is limited by the jet boundary/shear layer. Further rise in $V_c$ can only occur if particles experience a significant radial dispersion and leave the central region of the jet. We infer that, at the inflexion point (inside the grey region in figure 5), the particles’ direction of travel develops a radial component. Beyond this, $V_c$ increases, but particles are now facing larger eddies with lower turbulent fluctuations (see figure 3). Additionally, now that the particles are out of the higher velocity central region of the jet that houses smaller eddies (recall figure 4) they lose their axial momentum, which causes particles to react and cluster at larger eddies. This is reflected by the drop in $D_2$ . Evidently, from $\mathrm{P1}$ to $\mathrm{P4}$ the $D_2$ value of the inflexion point drops and therefore the scales at which their corresponding clusters form increase. Based on these figures we hypothesise that each particle species leaves the dominantly longitudinal travel in the jet until it meets an eddy corresponding to its inertia (depicted by the inflexion point), which pushes it out in the lateral direction. Therefore, this inflexion point is a clear indication that particles start spreading radially based on the eddy they dominantly respond to, as they travel downstream.

3.2. Finite-time Lyapunov exponent

To quantify the local eddies that trigger the lateral particle dispersion, an approximation of the axial location ( $x_r$ ) where particles first leave the jet core is required. Although we could simply time average that radial particle velocity, this statistic smears out some of the important dynamics. To get a better sense for the radial dispersion, the rate of convergence and divergence of particle trajectories along the jet centreline is probed via the finite-time Lyapunov exponent (FTLE) (Haller Reference Haller2001). The FTLE determines the local attractor points by estimating the temporal evolution of the separation between particle pairs. In other words, FTLE characterises the evolution of particle separation distance over time. If particles only travel within the jet core, their radial velocity component will be substantially smaller than their axial velocity. This will result in small inter-particle separation in the radial direction. By examining the spike in the particle separation distance determined with FTLE, we can estimate the location where particles dominantly start travelling radially. Note that a major rise in FTLE will only occur if particles start dispersing radially because the axial fluid–particle slip velocity is very high; the particles injected in the jet at the inlet reach the outlet at approximately the same time if they stay within the central region of the jet.

We defined FTLE ( $\lambda (x)$ ) with the following mathematical expression:

(3.5) \begin{equation} \lambda (x ) = \frac {1}{\Delta t} \log \left | \frac {\delta r(x) }{ \delta r (x_0)} \right | ,\end{equation}

where $\delta r(x)$ is the separation between the trajectories of a particle pair at a given time, while $\delta r(x_0)$ depicts the initial separation between neighbouring particles. The value of $\delta r(x)$ is computed from the flow map, $\phi (x, u, t)$ , which is a function of particle position, velocity and time instance (Günther & Theisel Reference Günther and Theisel2016). The $\phi (x, u, t)$ maps a particle’s travel from the initial condition to the final destination after integrating along the particle trajectory for a given time duration. Thus, to find $\lambda (x)$ in the streamwise direction, we compute $\delta r(x)$ at each $x$ -location using the separation distance between two particle trajectories in the $y$ - $z$ plane at each $x$ -location. To better understand the evolution of the particle pair separation, in this case, $\delta r(x_0)$ is taken as the separation of the particles at the previous $y$ - $z$ plane. Thus, $\lambda (x)$ provides a local estimation of the particle separation at each $x$ -location. For this analysis, we temporally averaged $\lambda (x)$ of all the pairs of each particle species on a given $y$ - $z$ plane.

Figure 7. The FTLE ( $\lambda (x)$ ) computed along the streamwise direction, where $\lambda (x)$ is averaged over all pairs of a particle species on the given $y$ - $z$ plane along $x$ -axis to obtain a single curve. The red dotted line shows the axial location ( $x_r$ ) of peak particle dispersion.

The values of $\lambda (x) \gt 0$ and $\lambda (x) \lt 0$ exhibit divergence and convergence of inertial particle trajectories, respectively. From this, it can be inferred that the axial location of a positive peak corresponds to the maximum divergence of the particle trajectories, and it depicts the position of the first dominant radial dispersion ( $x_r$ ) of the particles. This trend can be observed in figure 7, where the peak $\lambda (x)$ is farther downstream with increasing particle size, such that the locations of the initial radial dispersion for $\mathrm{P1}$ , $\mathrm{P2}$ , $\mathrm{P3}$ and $\mathrm{P4}$ are $x_r = 3.3$ , $5.5$ , $6.8$ and $7.5$ , respectively. Here, the $\lambda (x)$ drops after $x_r$ . This is because, although particles are spreading, there are still particle pairs that are coming close to each other as they spread radially. Therefore, we see significant fluctuations in the present $\lambda (x)$ curves in figure 7.

Also note that similar to $D_2$ in § 3.1, $\lambda (x)$ values of these peaks drop from $\mathrm{P1}$ to $\mathrm{P4}$ , showing an overall drop in the chaotic dispersion (rise in clustering) of the particles. A similar trend can be seen in the two-dimensional (2-D) contour of $\lambda (x)$ computed within a thin slice at the centre of the domain, as depicted in figure 8. We only show the positive values of FLTE in figure 8 to reduce the noise that manifests from the fact that particles are mainly present at the centre of the domain, resulting in fewer particle trajectories, while most of the domain is empty. Moreover, we are only concerned with $\lambda (x) \gt 0$ , which quantifies radial dispersion. The depiction of $\lambda (x) \gt 0$ can also imply that in a particle pair, one particle is dispersing radially, while the other stays close to the jet centre. However, in the present discussion, we use $\lambda (x) \gt 0$ to represent dispersion because even in this case, at least one particle (of the particle pair) is still dispersing. Albeit the contours in figure 8 only show the central region of the domain, the location of the peak FTLE agrees reasonably well with figure 7.

Figure 8. Contour of FTLE ( $\lambda (x)$ ) in a thin slice at the centre of the domain. Here, only positive (dispersion) values of FTLE are shown. The dotted red colour line depicts the $x_r$ of the corresponding particle species.

3.3. Local particle characteristics

Now that $x_r$ has been determined, we can characterise the turbulent eddies that exist at these locations. But first, we have to characterise particle parameters at $x_r$ that cause particles to react to the local eddies. The most important of these parameters is the particle local Stokes number ( ${\textit{St}}^*$ ). The value of ${\textit{St}}^*$ is computed from the fluid variables at the particle position ( $x_{\!p}$ ) and is defined as

(3.6) \begin{equation} St^* = \frac {\tau _{\!p} (x_{\!p})}{\tau _{\textit{eddy}}} = \frac { \dfrac {\rho _{\!p} d_{\!p}^2} {18 \, \rho _g(x_{\!p}) \, \nu _g(x_{\!p})}} {\dfrac {L_{\textit{eddy}}} {u_{\textit{rms}}(x_{\!p})}} ,\end{equation}

where $\tau _{\!p} (x_{\!p})$ is the local inertial time scale of the particles, computed from the local fluid parameters (not the bulk fluid properties), as shown in (3.6), whereas, $\tau _{\textit{eddy}}$ is the characteristic time scale of the representative eddies at $x_{\!p}$ . Note that $L_{\textit{eddy}}$ is computed at each grid point along the $x$ -axis using the 2-D $y$ - $z$ slice; the details of the computation have been previously discussed in § 3.1. Here, $u_{\textit{rms}} (x_{\!p})$ is the root mean square of the local fluid velocity at the particle position ( $x$ , $y$ , $z$ ). Since ${\textit{St}}^*$ is defined with the local fluid parameters, the value of ${\textit{St}}^*$ will be different for each particle. To deal with this, we compute the average value of ${\textit{St}}^*$ at each $y$ - $z$ plane along the $x$ -direction for each particle species, and then we also temporally averaged ${\textit{St}}^*$ .

Starting from the initial local Stokes numbers of $1.5$ , $2.7$ , $4.2$ and $6.0$ for $\mathrm{P1}$ to $\mathrm{P4}$ at the inlet (see table 1), the averaged ${\textit{St}}^*$ of each particle species drops to ${\textit{St}}^* \approx 0.6$ at $x_r$ , as depicted in figure 6. The decline in ${\textit{St}}^*$ is expected as $L_{\textit{eddy}}$ and the corresponding eddy time scale ( $\tau _{\textit{eddy}}$ ) increase with the downstream distance. According to § 3.1, $L_{\textit{eddy}}$ depends on the threshold applied on $\mathcal{Q}$ -criterion to determine the dominant vortices. As a result, our ${\textit{St}}^*$ is also sensitive to this threshold. This value of ${\textit{St}}^* \approx 0.6$ is unexpected because, if ${\textit{St}}^* \approx 1.0$ were determined at $x_r$ , particle lateral diffusion could be attributed solely to PC, since in this case particle dispersion is characterised by the centrifugal mechanism (Uthuppan et al. Reference Uthuppan, Aggarwal, Grinstein and Kailasanath1994). This begs the question: Is it possible that because particles are in a supersonic flow facing a high drag force in the axial direction, they might experience a delay in responding to the local eddies even when their ${\textit{St}}^* \approx 1.0$ , and that this causes them to disperse at ${\textit{St}}^* \lt 1.0$ ? This explanation is plausible considering that the particle momentum is mainly in the axial direction as particles are injected at the inlet with random velocity in the $x$ -direction. But then why does their ${\textit{St}}^*$ converge to $0.6$ at $x_r$ of each particle species, in comparison with any value less than unity?

Figure 9. Particle axial acceleration ( $a_{p,x}$ ) normalised with the initial value ( $a_{p,x,0}$ ). The dotted lines represent the $x_r$ of the corresponding particle species.

The above-stated lag in the particles’ response to the local eddies can be characterised by the evolution of particle axial acceleration ( $a_{p,x}$ ). Since particles in the centre of the jet have a higher velocity, as these particles move radially, they travel into the lower momentum region of the flow, while possessing higher momentum. Thus, they have a higher slip, which drops their $a_{p,x}$ . If the downstream location of this peak occurs slightly before $x_r$ , then this potentially indicates that, based on particle momentum, particles take some time to correlate to the local eddies, and during this time their ${\textit{St}}^*$ drops to around $0.6$ . We can see this in figure 9, where we have marked these downstream peak acceleration values with circular markers. For each particle species, the circular marker is observed before $x_r$ , and the distance between the downstream $a_{p,x}$ peak and $x_r$ increases with the particle size.

Apart from the above discussion, the acceleration trends in figure 9 also agree with the past literature, where smaller particles initially accelerate rapidly and then experience a sharp drop in acceleration (Yu et al. Reference Yu, Yang, Hu and Wang2023). While larger particles accelerate at a slower rate (Sommerfeld Reference Sommerfeld1994; Sakakibara, Wicker & Eaton Reference Sakakibara, Wicker and Eaton1996), and their large inertia (Chen et al. Reference Chen, Zhang, Li and Pan2024) causes their acceleration to be prolonged. These effects also cause a gradual drop in the acceleration of larger particles as they travel downstream. This is a credible quantitative check, suggesting that our hypothesis regarding the lag experienced by the particles is reasonable.

Figure 10. Magnitude of the local gas velocity gradient at the particle position in the radial direction, $\Vert\partial u_{g,x} (x_{\!p})/\partial r\Vert = \sqrt {(\partial u_{g,x} (x_{\!p})/ \partial y)^2 + (\partial u_{g,x} (x_{\!p})/ \partial z)^2 }$ , non-dimensionalised with the respective particle aerodynamic time scale ( $\tau _{\!p}$ ). Here, $\Vert\boldsymbol{\cdot }\Vert$ represents the magnitude, and $\vartheta$ in the legend shows the exponential decay constant of the corresponding particle species.

3.4. Local eddy characteristics

Due to momentum diffusion, jets spread laterally as they evolve downstream; thus, particles located within the jet will experience a non-zero radial velocity. This results in a radial drag force on the particles, as it impacts the particle–fluid relative velocity ( $\boldsymbol{u}_{\!p} - \boldsymbol{u}_{\!f}(x_{\!p})$ ). Additionally, since the maximum jet velocity at any axial location is at the centre of the jet and decreases radially, there exists an axial velocity gradient in the radial direction. In this case, particles tend to migrate towards low shear regions, to sample the zero-acceleration points where the net force on them is zero (Chen, Goto & Vassilicos Reference Chen, Goto and Vassilicos2006). By examining the exponential decay constant ( $\vartheta$ ) of this temporally and azimuthally averaged local streamwise gas velocity gradient in the radial direction ( $\Vert\partial u_{g,x}(x_{\!p}) / \partial r\Vert$ ), we can estimate how slowly particles respond to this radial push. This is shown in figure 10, where $\vartheta$ decreases from $\mathrm{P1}$ to $\mathrm{P4}$ . In this case, $\vartheta$ is computed by curve fitting the local gas velocity gradient in the radial direction to an exponential decay function. For an arbitrary variable $\varPi$ , the exponential decay function is $\varPi (x) = \varPi _0 \, e^{-\vartheta x}$ .

The higher inertia of $\mathrm{P4}$ particles makes them less responsive to this local velocity gradient, causing the local velocity gradient to be high and $\vartheta$ to be low. Another interesting observation is that each particle curve plateaus roughly around $x_r$ of the corresponding particle. Apart from supporting the $x_r$ determined in § 3.2, $\vartheta$ and the location of the plateau of these curves can help predict the transition between the streamwise travel and radial particle dispersion in jet flows. This is because, unlike FTLE, which requires the computation of particle trajectories over time, $\Vert\partial u_g(x_{\!p}) / \partial r\Vert$ and $\vartheta$ are rather inexpensive and straightforward to compute. Thus, $x_r = x$ such that $ d/{\rm d}x \Vert\partial u_g(x_{\!p}) / \partial r\Vert$ approaches zero.

Figure 11. Contour of temporally averaged $\mathcal{Q}(y,z)$ on a 2-D thin slice in $y$ - $z$ planes at $x_r = 3.3$ , $5.5$ , $6.8$ and $7.5$ of the $\mathrm{P1}$ , $\mathrm{P2}$ , $\mathrm{P3}$ and $\mathrm{P4}$ particles, respectively. Where, $\mathcal{Q}_{max}(y,z)$ is the maximum value of $\mathcal{Q}(y, z)$ at each $x_r$ plane.

Figure 12. Mean 1-D energy spectrum computed at the $y$ - $z$ planes at $x_r=3.3$ , $5.5$ , $6.8$ and $7.5$ of the $\mathrm{P1}$ , $\mathrm{P2}$ , $\mathrm{P3}$ and $\mathrm{P4}$ particles, respectively.

To further understand the local eddy characteristics at $x_r$ of each particle species, we plot the 2-D contours of the $\mathcal{Q} (y,z)$ -criterion at $x_r$ of the considered particle sizes, as shown in figure 11. As expected, a noticeable expansion of the jet can be witnessed in this figure from $\mathrm{P1}$ to $\mathrm{P4}$ . Due to the shear layer at the outer edge of the jet, the vorticity of the larger eddies away from the jet centre is higher, creating a higher radial push on the particles. Hence, it is plausible that the decline in particle ${\textit{St}}^*$ and the expansion of $\mathcal{Q}(y,z)$ depicting higher energy eddies at the periphery of the jet create a correlation between the local eddies and the particles at their corresponding $x_r$ .

We can quantify this rise in the energy of larger eddies with the one-dimensional (1-D) energy spectrum computed from the two lateral fluctuation velocity components of the fluid at $x_r$ of each particle species on the $y$ - $z$ plane, as depicted in figure 12. It can be observed in this figure that at $x_r = 3.3$ smaller eddies contain higher energy, while at $x_r = 7.5$ larger eddies on the $y$ - $z$ plane possess higher energy. These local eddies at $x_r$ essentially expel particles from the jet via centrifugal force.

Figure 13. Second-order structure function of the local gas velocity ( $S^2_{u_g(x_{\!p})}(\Delta r)$ ) at the particle position, evaluated at the $y$ - $z$ planes at $x_r=3.3$ , $5.5$ , $6.8$ and $7.5$ of the $\mathrm{P1}$ , $\mathrm{P2}$ , $\mathrm{P3}$ and $\mathrm{P4}$ particles, respectively. The second-order structure function is normalised with its first non-zero value $S^2_{u_g (x_{\!p})} (\Delta r)_0$ .

Since larger eddies grow downstream, there is a steeper energy change between scales, which is crucial for radial particle dispersion. This energy distribution across scales can be assessed with the second-order structure function ( $S^2_{u_g(x_{\!p}) }(\Delta r)$ ) of local gas velocity (at the particle position ( $x_{\!p}$ )), which is the mean square of the velocity difference between two points with a separation distance $r$ . We evaluate $S^2_{u_g(x_{\!p})}(\Delta r)$ at the $y$ - $z$ plane at $x_r$ of each particle species, thus the $S^2_{u_g(x_{\!p})}(\Delta r)$ is mathematically given as

(3.7) \begin{equation} S^2_{u_g(x_{\!p})} (\Delta r) = \overline {(u_g(x_{\!p}, r + \Delta r) - u_g(x_{\!p}, r) )^2}, \ \ \mathrm{where} \ \ r=\sqrt {y^2 + z^2} .\end{equation}

In the above equation, the overline depicts the temporal average. Hence, if $S^2_{u_g(x_{\!p})}(\Delta r)$ is higher, then at the given $y$ - $z$ plane the velocity gradient across particle positions is higher. This will create a higher drag force on the particles in the radial direction. We can witness in figure 13 where the $S^2_{u_g(x_{\!p})}(\Delta r)$ curves increase with particle size, proving that particles are pulled out of the axial motion once they encounter eddies of higher energy. Note that $S^2_{u_g(x_{\!p})}(\Delta r) = 0$ at $r=0$ . Therefore in figure 13, $S^2_{u_g(x_{\!p})}(\Delta r)$ is normalised with its first non-zero value denoted as $S^2_{u_g(x_{\!p})}(\Delta r)_0$ .

3.5. Particle–jet interaction

Thus far, we have scrutinised the effect of gas on the particles as it directly governs particle dispersion. Now we focus our attention on the modulation of the jet by the inertial particles in the near field. This is crucial as the gas and particles are mutually coupled. Therefore, the local turbulence not only promotes particle dispersion, but the presence of particles also modulates the local eddy structures of the jet.

The local gas velocity experiences the immediate impact of particles, as particles lag behind the gas in the near field due to their inertia. On the one hand, a larger particle momentum means that the particles are less impacted by the small-scale turbulence fluctuations. On the other hand, the larger size of the particles promotes the correlation between particle and fluid motion via the increase in drag force on the particles. Note that we increase particle inertia by increasing $d_{\!p}$ instead of $\rho _{\!p}$ . Of course, the balance of these two effects will govern the particle acceleration as we witnessed in figure 9. In either case, particles will considerably decelerate the local gas. To analyse this, we start by assessing the joint probability distribution function ( $\mathrm{PDF}$ ) of $u_{\!p}$ and $u_g(x_{\!p})$ magnitude over the entire domain, as shown in figure 14. The correlation between $u_{\!p}$ and $u_g(x_{\!p})$ is anticipated to drop with the rise in particle size. However, we witness a bimodal distribution in this figure. From $\mathrm{P1}$ to $\mathrm{P4}$ , the $\mathrm{PDF}$ s transition from almost linear growth to where $u_g(x_{\!p})$ is high while $u_{\!p}$ is low. This highlights a significant lag between the particle and gas velocities, which is proportional to the drag force on the particle. This enhanced drag force on the larger particles significantly impacts the local gas (recall (2.10)).

Figure 14. Joint $\mathrm{PDF}$ of the particle ( $u_{\!p}$ ) and local gas $(u_g(x_{\!p}))$ velocities.

Figure 15. Turbulent kinetic energy ( $\mathrm{TKE} (x_{\!p})$ ) and dissipation rate ( $\varepsilon (x_{\!p})$ ) at the location of particles.

To probe the particle-based local modulation of the turbulence, we compute the turbulent kinetic energy ( $\mathrm{TKE} (x_{\!p})$ ) and dissipation rate ( $\varepsilon (x_{\!p})$ ) of the gas at the particle position, defined as

(3.8) \begin{align}&\quad\,\, \mathrm{TKE} (x_{\!p}) = \frac {1}{2} \overline {u^{\prime}_{g,i}(x_{\!p}) \ \ u^{\prime}_{g,i} (x_{\!p})} , \end{align}

(3.9) \begin{align}& \varepsilon (x_{\!p}) = 2 \overline { \nu _g (x_{\!p}) \left |\frac {\partial u^{\prime}_{g,i} (x_{\!p})}{\partial x_k} \ \frac {\partial u^{\prime}_{g,i} (x_{\!p})}{\partial x_k} \right |} , \end{align}

where $u^{\prime}_{g,i} =u_{g,i}- \overline {u_{g,i}}$ in the above equations represents the fluctuating component of the fluid velocity in the $i$ th direction. The averaging, $\overline {u_{g,i} (x,r)}$ , is done over time and in the azimuthal direction.

As the joint $\mathrm{PDF}$ s indicated (see figure 14), the $\mathrm{TKE} (x_{\!p})$ increases with the particle size, as shown in figure 15. Note that, unlike the overall $\mathrm{TKE}$ of the gas in figure 6, $\mathrm{TKE} (x_{\!p})$ dwindles downstream. This is reasonable as the slip velocity drops axially, and the back-reaction force (producing fluctuations in the local gas velocity) is the only force particles exert on the local gas, which depends on the particle–gas slip velocity (recall (2.10)). Thus, this drop in fluid–particle slip velocity downstream reduces the fluctuations particles induce in the local fluid velocity field, which results in $\mathrm{TKE} (x_{\!p})$ decline in the $x$ -direction. A more plausible explanation for this trend is that as the particles disperse radially downstream, they enter regions with low $\mathrm{TKE}(x_{\!p})$ . For example, according to figure 3, the flow $u_{\textit{rms}}$ increases downstream at the centreline, and eventually reaches a maximum value at $x\gt 6$ . On the other hand, the $\mathrm{TKE} (x_{\!p})$ reaches its maximum value around $x \approx 4$ and then decreases. Since most particles have started dispersing radially around $x\approx 6$ , they are not present in the higher $\mathrm{TKE}$ (central) region of the jet. Thus, $\mathrm{TKE}(x_{\!p})$ drops downstream. The increase in $\mathrm{TKE} (x_{\!p})$ until $x=4$ can be ascribed to the higher slip velocity and acceleration of particles close to the jet inlet, as they are introduced with a random velocity in the $x$ -direction at the inlet. We see a steeper $\varepsilon (x_{\!p})$ decay after $x=4$ similar to that of $\mathrm{TKE} (x_{\!p})$ . These trends do not resonate with the classical trends where smaller particles dissipate more energy, while larger particles enhance energy and consequently reduce the decay rate of the jet (Njue et al. Reference Njue, Salehi, Lau, Cleary, Nathan and Chen2021). This is because we are analysing the local near-field effects on the gas. In the far field, where gas velocity will experience a significant drop, smaller particles will decelerate with the gas, further taking energy from the flow. But the larger particles will retain their axial momentum. Therefore, beyond a certain axial position, their velocity will exceed the local gas velocity. In this case, they will start dragging the gas, imparting their kinetic energy onto the gas. Thus, it is clear that close to the inlet, the local turbulence rapidly enhances with particle size in the regime where particles experience higher acceleration via the drag force.

Figure 16. Relative kinetic energy of the particles, $\Delta \mathrm{KE}$ (left), and axial component of particle–gas slip velocity, $u_{\textit{pg}}$ (right).

Naturally, we expect that the rise in local turbulent fluctuations in the gas caused by the particles should also influence the particle dynamics. In this case, smaller particles should adapt to these fluctuations more rapidly than their larger counterparts. This can be scrutinised with the particle–gas relative kinetic energy ( $\Delta \mathrm{KE}$ ) and slip velocity in the $x$ -direction ( $u_{\textit{pg}}$ ), which are defined as

(3.10) \begin{align}& \Delta \mathrm{KE} = \frac {1}{2} m_{\!p} \overline {\big ( u^{\prime}_{p,i} - u^{\prime}_{g, i} (x_{\!p}) \big ) \big( u^{\prime}_{p,i} - u^{\prime}_{g, i} (x_{\!p}) \big)} , \end{align}

(3.11) \begin{align}&\qquad\qquad\quad u_{\textit{pg}} = \overline {u_{\!p, x} - u_{g, x} (x_{\!p})} . \end{align}

These computations are shown in figure 16. The negative sign of $u_{\textit{pg}}$ suggests that particles take energy from the local gas in the $x$ -direction and never drag the gas in the considered near-field region, as discussed earlier. We are only showing the axial component of $u_{\textit{pg}}$ to preserve its sign, while its radial components are negligibly small as compared with the axial component. Note that $\mathrm{P1}$ particles depict the highest slip velocity at $x\lt 2$ . This is because particles are introduced at the jet inlet with a random velocity in the $x$ -direction. This effect rapidly disappears as $\mathrm{P1}$ particles promptly adapt to the local gas velocity, which causes their slip velocity to drop.

Interestingly, the difference in $u^{\prime}_{p,i}$ and $u^{\prime}_{g,i} (x_{\!p})$ increases downstream, as indicated by $\Delta \mathrm{KE}$ , while the difference between their mean components drops. Moreover, this increase reaches its peak at $x_r$ . This suggests that, regardless of particle size (in the realm of finite particle inertia), particle velocity fluctuations do not adjust to the local fluctuations and therefore particles are quasi-ballistic to turbulent fluctuations. It is the particle mean velocity component that adapts to the local turbulent fluctuations. Despite the two-way coupling, comparing $\mathrm{TKE} (x_{\!p})$ and $\Delta \mathrm{KE}$ in figures 15 and 16, respectively, implies that the effect of the particles on the gas and that of gas on the particles is rather dissimilar. For the present ${\textit{St}}$ range, particles mainly alter the fluctuating component of the gas velocity, whereas the gas influences the mean velocity component of the particles. Since the fluid–particle interaction is a function of ${\textit{St}}$ (Ferrante & Elghobashi Reference Ferrante and Elghobashi2003), this effect should be tested on a broader range of ${\textit{St}}$ . It should be noted that we are referring to the dominant effect here. Of course, in a holistic view, both velocity components experience changes. The peak $\Delta \mathrm{KE}$ at $x_r$ is due to the rise in energy of the local eddies which pushes particles out of the jet’s central region. Once outside this region, $\Delta \mathrm{KE}$ drops because the local velocity fluctuations dwindle rapidly.

The smaller particles are more responsive to the local eddies (Squires & Eaton Reference Squires and Eaton1991), which promotes the gas–particle momentum exchange. The overall slip velocity in figure 16 demonstrates this trend. Furthermore, since the average particle velocity in the domain is inherently slow in responding to the local changes, we see that the inflexion point of $u_{\textit{pg}}$ occurs after the $x_r$ plane, unlike $\Delta \mathrm{KE}$ . These two tendencies affirm the validity of our previous observations.

The fluid–particle slip velocity can further validate the $x_r$ of each particle species. For this purpose, we define a second-order neighbouring particle slip velocity function ( $\chi ^2_{u_{\textit{pg}} } (x)$ ), which is given as

(3.12) \begin{equation} \chi ^2_{u_{\textit{pg}} } (x) = \overline { \left ( u_{\textit{pg}}(x + \varsigma ) - u_{\textit{pg}}(x) \right )^2 } ,\end{equation}

where $\varsigma$ represents the streamwise separation between two neighbouring particles. Note that the streamwise separation distance varies between particle pairs, as the inter-particle distance can differ at each instance and based on particle location. By averaging the slip velocity difference between neighbouring particles at each streamwise location, the relative change in slip velocity at each streamwise location in the jet is determined. The resulting $\chi ^2_{u_{\textit{pg}}} (x)$ is shown in figure 17, which depicts that the correlation of $u_{\textit{pg}}$ increases with the particle size. This is expected because larger particles retain their inertia for longer, and in the near-field region, the gas velocity is still very high. The interesting aspect is that, starting from $x = 0$ , each curve starts with a shallower slope, and then, at a certain $x$ , its slope becomes steeper. This shallower slope is also longer for larger particles. The change in slope of these curves occurs close to the $x_r$ of the corresponding particles. This analysis reinforces the credibility of $x_r$ .

Figure 17. Second-order neighbouring particle slip velocity function ( $\chi ^2_{u_{\textit{pg}}}(x)$ ). The dotted lines on the plot show the slope.

The pressing question at this point is, how do turbulent fluctuations of the local gas and particle velocity fluctuations in the axial and radial directions correlate? The correlation between particle and gas velocity fluctuations can be probed with the local gas mean Reynolds stress ( $\overline {u^{\prime}_i u^{\prime}_{\!j}}$ ) and the averaged product of particle velocity fluctuations. We acknowledge that, for particles, $\overline {u^{\prime}_i u^{\prime}_{\!j}}$ represents the correlation between particle fluctuating velocity components and cannot be called Reynolds stress. Here, $\overline {u^{\prime}_i u^{\prime}_{\!j}}$ defines the momentum transport of the gas and particles due to the local eddies. The normal stresses ( $\overline {u^{\prime}_i u^{\prime}_i}$ ) for both gas and particles should be positively correlated as the main thrust of the flow is in the $x$ -direction. Additionally, the jet is spreading laterally downstream, and particles are shadowing the gas with respect to their inertia. The top two plots of figure 18 illustrate this trend in the normal stress components. Here, particle fluctuations do not increase to the same level as the gas fluctuations, similar to the trend observed in figure 16, but they do show a tendency to escalate akin to the gas fluctuations. We also see this trend in $\overline {u^{\prime}_y u^{\prime}_z}$ (lateral components), where particle velocity fluctuations try to match the fluctuations in the gas velocity, albeit much more slowly, as particle axial inertia does not prefer the radial momentum transfer. It should be noted that the axial location of peak values of particle $\overline {u^{\prime}_y u_y'}$ and $\overline {u^{\prime}_z u_z'}$ correlates well with $x_r$ , which implies higher transport of particles in the radial direction around $x_r$ . This validates the robustness of $x_r$ for predicting particle radial dispersion in a jet.

In terms of the momentum transfer between the axial and lateral components ( $\overline {u^{\prime}_x u^{\prime}_y}$ and $\overline {u^{\prime}_x u^{\prime}_z}$ ), the downstream radial expansion of gas implies negative Reynolds stress, suggesting the growth in lateral fluctuations, as depicted in figure 18. However, particle axial to radial momentum transfer via velocity fluctuations should be minimal compared with the gas because of particle inertia in the axial direction. This momentum transfer should further drop with the increase in particle size, as shown in figure 18. These two phenomena are intuitive, but the compelling observation is that particles are almost entirely ballistic to the axial to lateral velocity fluctuations of the local gas. This disconnect between the gas Reynolds stresses and particle velocity fluctuations signifies that particle radial momentum is exclusively generated by the lateral combination of the gas velocity fluctuations. Thus, a turbulent eddy must align in a specific manner to exert centrifugal force dominantly in the radial direction for particles to be pulled out of the central region of the jet.

Figure 18. Comparison of the averaged product of particle fluctuating velocity components and mean Reynolds stresses ( $\overline {u^{\prime}_i u^{\prime}_{\!j}}$ ) of the local gas velocity fluctuations ( $u^{\prime}_{\!p}(x_{\!p})$ ). Note that, due to the symmetry of the jet, $\overline {u^{\prime}_y u^{\prime}_y} \approx \overline {u^{\prime}_z u^{\prime}_z}$ and $ \overline {u^{\prime}_x u^{\prime}_y} \approx \overline {u^{\prime}_x u^{\prime}_z}$ . Therefore, we only show $ \overline {u^{\prime}_y u^{\prime}_y}$ (top-right) and $ \overline {u^{\prime}_x u^{\prime}_y}$ (bottom-left) curves in this plot.

3.6. Effect of compressibility on dispersion

In a supersonic flow, the compressibility effects also govern the dispersion of particles (Pérez-Munuzuri Reference Pérez-Munuzuri2015). To assess the effect of compressibility on particle transport, we evaluate the magnitude of local density gradient and divergence of the gas velocity ( $u_g (x_{\!p})$ ) at the location of particles (dilatation of the local fluid). Due to their higher inertia, larger particles are unaffected by small variations in density. We see this trend in figure 19, where smaller particles prefer regions with the probability of minimal density change, while larger particles favour regions with high spatial change in the fluid density. Note that we also have negative values for our magnitude. This is because we preserved the sign of each component by multiplying it with the magnitude value. In this figure, the regions with negative density gradients also have significant dilatation. Hence, large particles accumulate in areas with a higher probability of dilatation, as shown in figure 19. This is potentially because the regions with strong dilatation have low vorticity, where particles tend to cluster due to PC. In other words, strong dilatation corresponds to stagnation points in the flow, where particles form clusters as the net force on them at these locations is zero (Chen et al. Reference Chen, Goto and Vassilicos2006).

Figure 19. The $\mathrm{PDF}$ of the magnitude of gas density gradient (left) and dilatation (right) computed at the position of particles. We preserved the sign of each component before taking the magnitude and multiplying these signs with the final magnitude value, which gave us a negative sign for the magnitude. The dotted red line represents the Gaussian curve.

Figure 20. Number of particles ( $n_s$ ) normalised with the total number of particles of a given species ( $n_{\!p}$ ) in the entire domain, sampled in the regions where $\rho _g (x_{\!p}) \gt \Vert\rho _g + \sigma (\rho _g)\Vert^*$ (left), $\boldsymbol{\nabla }\boldsymbol{\cdot }\rho _g (x_{\!p}) \gt \Vert \boldsymbol{\nabla }\boldsymbol{\cdot }\rho _g + \sigma (\boldsymbol{\nabla }\boldsymbol{\cdot }\rho _g)\Vert^*$ (middle), $\boldsymbol{\nabla }\boldsymbol{\cdot }u_g (x_{\!p}) \gt \Vert\boldsymbol{\nabla }\boldsymbol{\cdot }u_g + \sigma (\boldsymbol{\nabla }\boldsymbol{\cdot }u_g)\Vert^*$ (right). Here, $\sigma$ is the standard deviation, and the superscript $*$ represents that the curves are normalised with their maximum value.

Next, we quantify the number of particles of each species ( $n_s$ ) present in regions with density ( $\rho _g(x_{\!p}) \gt \Vert \overline {\rho _g} + \sigma (\rho _g)\Vert^*$ , here, the superscript $*$ depicts that the curve is normalised with its maximum value), density gradient ( $\boldsymbol{\nabla }\boldsymbol{\cdot }\rho _g(x_{\!p}) \gt \Vert \overline {\boldsymbol{\nabla }\boldsymbol{\cdot }\rho _g} + \sigma (\boldsymbol{\nabla }\boldsymbol{\cdot }\rho _g)\Vert^*$ ) and dilatation ( $\boldsymbol{\nabla }\boldsymbol{\cdot }u_g(x_{\!p}) \gt \Vert \overline {\boldsymbol{\nabla }\boldsymbol{\cdot }u_g} + \sigma (\boldsymbol{\nabla }\boldsymbol{\cdot }u_g)\Vert^*$ ) greater than spatial mean plus one standard deviation ( $\sigma$ ). The resulting number of particles in the regions that fulfil these conditions is presented in figure 20. Firstly, following the trends in literature, qualitatively larger particles reside in regions with elevated gas density (Zhang et al. Reference Zhang, Liu, Ma and Xiao2016; Xiao et al. Reference Xiao, Jin, Luo, Dai and Fan2020). The surprising aspect here is that $n_s$ of larger particles rises steeply as $\Vert\overline {\rho _g} + \sigma (\rho _g)\Vert^*$ increases, compared with the smaller particles. This is tied to the affinity of larger particles towards a higher-density gradient (see figure 20). Particles should naturally move to the location of the lower gradient. However, larger particles are slow in responding to the changes in the flow due to their higher inertia. On the other hand, smaller particles rapidly adapt to the variations in the flow and move to the low-density gradient and dilatation areas. Thus, larger particles reside in high-density gradient and dilatation regions. Note that regions with higher dilatation also exhibit low vorticity. Thus, larger particles inherently settle in the regions with high dilatation.