1. Introduction
The dispersion of heavy inertial particles in columnar vortices is a fundamental and widely occurring process. Several manufactured devices depend critically on the interaction between inertial particles and vortical flow structures, such as centrifugal separation, swirling biomass combustors and aircraft trailing vortices (Paoli et al. Reference Paoli, Vancassel, Garnier and Mirabel2008; Paoli & Shariff Reference Paoli and Shariff2016). The interaction between suspended inertial particles and vortical flow structures also plays a significant role in geophysical flows. Examples include dust devils, formed when a concentrated vortex tube lofts sand and debris (Zhao et al. Reference Zhao, Gu, Yu, Ge, Li and Feng2004; Fenton et al. Reference Fenton, Reiss, Lemmon, Marticorena, Lewis and Cantor2016; Simpson & Glezer Reference Simpson and Glezer2016), and dune morphology, which is impacted by the interaction between suspended sand grains and coherent vortical structures above dune beds (Müller & Gyr Reference Müller and Gyr1986; Best Reference Best2005; Bristow et al. Reference Bristow, Blois, Best and Christensen2020). In flows where the concentration of the disperse phase is very low, the fluid and particle phases may be partially decoupled: whereas the dispersion of particles is determined by the carrier flow, the latter is insensitive to the particle phase, thus evolving as in particle-free conditions (Elghobashi & Truesdell Reference Elghobashi and Truesdell1992). This flow regime is known as the one-way coupling regime in the classification of Elgobashi (Reference Elgobashi2006). However, when the concentration of the disperse phase is large enough, the feedback force from the particles on the fluid may result in significant modulation of the suspending flow (Elghobashi & Truesdell Reference Elghobashi and Truesdell1993; Ferrante & Elghobashi Reference Ferrante and Elghobashi2003; Elgobashi Reference Elgobashi2006). This regime, known as the two-way coupling regime, occurs in semidilute suspensions, characterized by large particle-to-fluid density ratio 
$\rho _p/\rho _f=O(10^3)$, particle volume fraction in the range 
$10^{-6} \leqslant \langle \phi _p \rangle \leqslant 10^{-3}$ and order-unity mass loading 
$M=O(1)$. Although it is well established that particles heavier than the suspending fluid tend to be ejected out of vortical regions by a mechanism known as preferential concentration (Squires & Eaton Reference Squires and Eaton1990, Reference Squires and Eaton2014; Wang & Maxey Reference Wang and Maxey1993), the joint evolution of the particle phase and carrier vortical flow under two-way coupling is poorly understood. To that end, we investigate the modulating effect of inertial particles on a representative vortex tube, namely the Lamb–Oseen vortex (Lamb Reference Lamb1945), in the semidilute regime. In particular, we analyse the effect of particle inertia on the decay of this elementary flow structure.
To help with the subsequent discussion, we adopt the Lamb–Oseen vortex as a prototypical isolated vortex tube with Gaussian distribution of vorticity (Lamb Reference Lamb1945). The axial and radial velocity components are zero, while the azimuthal velocity reads
\begin{equation} u_{f,\theta}(r, \theta)=\frac{\varGamma}{2{\rm \pi} r} \left( 1-\mathrm{exp}\left(-\frac{r^2}{r_0^2+4\nu t}\right) \right), \end{equation}
where 
$(r, \theta )$ are polar coordinates, 
$\varGamma$ is the circulation, 
$r_0$ is the initial core scale and 
$\nu$ is the kinematic viscosity (Lamb Reference Lamb1945). The velocity profile (1.1) is a self-similar solution to the Navier–Stokes equations. The strength of the vortex depends on the circulation Reynolds number 
$\textit {Re}_{\varGamma }=\varGamma /2{\rm \pi} \nu$. Experimental observation of trailing wing-tip vortices (Bhagwat & Leishman Reference Bhagwat and Leishman2002; Jacquin et al. Reference Jacquin, Fabre, Sipp, Theofilis and Vollmers2003; Paoli & Shariff Reference Paoli and Shariff2016) show that these vortex tubes are well approximated by the Lamb–Oseen profile (1.1), and that the flow remains largely laminar even at high Reynolds numbers due to the stabilizing effect of rotation (Spalart Reference Spalart1998). Consequently, the most significant characteristic of the Lamb–Oseen vortex is its remarkable stability.
The property is substantiated by several studies. Using linear stability analysis, Michalke & Timme (Reference Michalke and Timme1967) showed that a general vortex tube whose only vorticity extremum lies on the axis, such as the Lamb–Oseen vortex (1.1), is stable with respect to two-dimensional (2-D) perturbations. For a three-dimensional (3-D) vortex tube, Lessen, Singh & Paillet (Reference Lessen, Singh and Paillet1974) and Lessen & Paillet (Reference Lessen and Paillet1974) showed that the flow remains stable to axial perturbations. Moreover, in an attempt to map all possible normal modes for a 3-D Lamb–Oseen vortex, Fabre, Sipp & Jacquin (Reference Fabre, Sipp and Jacquin2006) showed that all these modes are stable, further emphasizing the extreme resilience of this vortex in constant-density single-phase flows. For an isolated vortex, the existence of instabilities requires flow stratification. Sipp et al. (Reference Sipp, Fabre, Michelin and Jacquin2005) showed that a 2-D instability of the Rayleigh–Taylor type could grow in a heavy-core Lamb–Oseen vortex. In the case of non-isolated vortex tubes, external strain fields applied on the Lamb–Oseen vortex may lead to axial linear instabilities (Moore, Saffman & Stuart Reference Moore, Saffman and Stuart1975; Sipp & Jacquin Reference Sipp and Jacquin2003).
For particle-laden vortex columns, several authors addressed the role of inertia on the dispersion of particles, under one-way coupling. Using direct numerical simulations and Lagrangian tracking, Marshall (Reference Marshall2005) studied the one-way coupled dispersion of inertial particles in a columnar vortex enveloped by turbulence. He found that inertial particles, initially dispersed within the vortex core, quickly form annular structures that swell radially outwards. Once these clusters expand beyond the protective rotational core, turbulence breaks them apart into smaller particle clusters. Ravichandran & Govindarajan (Reference Ravichandran and Govindarajan2015) analysed the clustering behaviour of inertial particles around point vortices, but did not consider particle feedback on the carrier fluid either. They showed that the particles accumulate into ring-shaped clusters (2-D analogues of the annular clusters in Marshall (Reference Marshall2005)) around point vortices, which increases the odds of particle–particle collision or coalescence. This process is speculated to be a driver of water droplet coalescence in warm clouds (Deepu, Ravichandran & Govindarajan Reference Deepu, Ravichandran and Govindarajan2017).
Despite being overlooked in most prior studies of particle-laden vortices, it is well established that two-way coupling may result in significant modulation of the carrier flow. The extent of the modulation is in part controlled by particle inertia, quantified by the Stokes number 
${St}=\tau _p/\tau _f$ and mass loading 
$M=\rho _p/\rho _f\langle \phi _p\rangle$ (Balachandar & Eaton Reference Balachandar and Eaton2010; Elghobashi Reference Elghobashi2019), where 
$\tau _p$, 
$\tau _f$, 
$\rho _f$, 
$\rho _p$ and 
$\langle \phi _p\rangle$ are the particle response time, the fluid characteristic time scale, fluid density, particle density and average particle volume fraction, respectively. Non-inertial particles, which are merely advected by the flow and follow Taylor's particle dispersion theory (Taylor Reference Taylor1922), do not have a significant modulating effect.
Conversely, particles with non-zero inertia deviate from the fluid streamlines, causing micro-stresses on the fluid (Maxey Reference Maxey1991). Collectively, the particles may exert a strong forcing capable of modulating flow structures on scales much larger than the particles (Kasbaoui Reference Kasbaoui2019). This mechanism is known as preferential concentration, which results in a bias in particle distribution towards straining regions of the flow over rotational ones (Squires & Eaton Reference Squires and Eaton1990, Reference Squires and Eaton2014; Wang & Maxey Reference Wang and Maxey1993). Particle clusters, several orders of magnitude larger than the particle diameter, may form as a result of preferential concentration (Gualtieri et al. Reference Gualtieri, Picano, Sardina and Casciola2013; Baker et al. Reference Baker, Frankel, Mani and Coletti2017; Kasbaoui, Koch & Desjardins Reference Kasbaoui, Koch and Desjardins2019a). Kasbaoui et al. (Reference Kasbaoui, Koch, Subramanian and Desjardins2015) showed that the collective effects of dispersed particles can activate non-modal instabilities in a linear shear flow, one that shares the same remarkable stability as the Lamb–Oseen vortex. The coupling between the two phases can also result in significant modulation of turbulent flows, including the modulation of turbulent coherent structures (Richter & Sullivan Reference Richter and Sullivan2014) and augmentation or attenuation of turbulence intensity (Ahmed & Elghobashi Reference Ahmed and Elghobashi2000; Gualtieri, Battista & Casciola Reference Gualtieri, Battista and Casciola2017; Kasbaoui Reference Kasbaoui2019; Kasbaoui, Koch & Desjardins Reference Kasbaoui, Koch and Desjardins2019b).
To the best of our knowledge, only Druzhinin (Reference Druzhinin1994, Reference Druzhinin1995) has considered the flow modifications induced by particles on circular flows. He employed the two-fluid model and solved the inviscid Euler equations with the assumption that the vortex flow retained an axial symmetry. The study was restricted to weakly inertial particles in 2-D columnar vortices with algebraic vorticity characterized by a ‘solid’ (constant-vorticity) rotation in the vicinity of the vortex centre and a potential (zero-vorticity) rotation at the periphery, which is a steady solution to the 2-D Euler equations. The assumption of weak particle inertia enables semi-analytical treatment in which the particle velocity field and concentration are determined from the carrier fluid velocity with a small inertial correction. Similarly to other authors, Druzhinin (Reference Druzhinin1994) showed that the particles form a ring-shaped cluster that expands outwardly from the core at a rate controlled by their inertia and leaves the core depleted of particles. For these weakly inertial particles, Druzhinin (Reference Druzhinin1994) suggests that particle dispersion should follow a spiralling motion caused by an inward fluid momentum flux resulting from the outward particle momentum flux.
This paper is organized as follows. In § 2, we present the governing equations and Eulerian–Lagrangian simulation method used in this study. The configuration and flow parameters are presented in § 3. In § 4, we review the characteristics of the Lamb–Oseen vortex in single-phase flows and use analytical expressions to validate our solution procedure. In § 5, we analyse particle-laden Lamb–Oseen vortices. We show that the presence of particles results in dramatically shortened life span of the vortex and provide scaling laws for the vortex core radius, vorticity, enstrophy as well as the growth rate of the void fraction bubble resulting from the outward migration of the particles. Concluding remarks are given in § 6.
2. Governing equations
We describe the dynamics of the carrier phase and suspended particles in the Eulerian–Lagrangian framework following the volume-filtered formulation of Anderson & Jackson (Reference Anderson and Jackson1967). The conservation equations for the carrier phase are given by the incompressible Navier–Stokes equations,
$$\begin{gather} \frac{\partial \phi_f}{\partial t}+\boldsymbol{\nabla}\boldsymbol{\cdot} (\phi_f\boldsymbol{u}_f) =0, \end{gather}$$
$$\begin{gather}\rho_f\left( \frac{\partial \phi_f\boldsymbol{u}_f}{\partial t}+\boldsymbol{\nabla}\boldsymbol{\cdot} (\phi_f\boldsymbol{u}_f \boldsymbol{u}_f)\right) ={-}\boldsymbol{\nabla} p+\mu_f\nabla ^2 \boldsymbol{u}_f+\boldsymbol{F}_p +\boldsymbol{\nabla} \boldsymbol{\cdot}\boldsymbol{R}_\mu, \end{gather}$$
where 
$\boldsymbol {u}_f$ is the fluid velocity, 
$\phi _f=1-\phi _p$ is the local fluid volume fraction, 
$\phi _p$ is the local particle volume fraction, 
$p$ is the pressure, and 
$\rho _f$ and 
$\mu _f$ are the fluid density and viscosity, respectively. The term 
$\boldsymbol {R}_\mu$ is a byproduct of filtering the fluid stress tensor, and is closed using an effective viscosity model, 
$\boldsymbol {R}_\mu =\mu _f(\phi _f^{-2.8}-1)[\boldsymbol {\nabla } \boldsymbol {u}_f + \boldsymbol {\nabla } \boldsymbol {u}_f^{\rm T}-\frac {2}{3}(\boldsymbol {\nabla }\boldsymbol {\cdot } \boldsymbol {u}_f) I]$ (Gibilaro et al. Reference Gibilaro, Gallucci, Di Felice and Pagliai2007). This term vanishes in the semidilute regime.
The coupling between phases occurs at two levels. A first pathway for flow modulation involves a reduction of volume available to the fluid due to the presence of particles.
This effect, which is captured by the volume fraction factor 
$\phi _f$ in the continuity and momentum equations, is negligible in the semidilute regime. Pakseresht & Apte (Reference Pakseresht and Apte2019) showed that effects due to volumetric displacement play a role when the particle volume fraction exceeds 5 %. In the present regime considered, the peaks of local particle volume fraction are of the order of 
$10^{-3}$, resulting in 
$\phi _f\simeq 1$ even within particle clusters. The dominant pathway for flow modulation occurs through momentum exchange between particles and gas, which is captured by the term 
$\boldsymbol {F}_p$. In the present work, momentum exchange is expressed as
\begin{equation} \boldsymbol{F}_p= {-\phi_p \boldsymbol{\nabla}\boldsymbol{\cdot} \tau|_p} +\rho_p\phi_p f_d\frac{(\boldsymbol{u}_p-\boldsymbol{u}_{f|p})}{\tau_p},\end{equation}
where 
$\boldsymbol{\tau} =-p\mathcal {I}+\mu (\boldsymbol {\nabla } \boldsymbol{u}_f+\boldsymbol {\nabla } \boldsymbol{u}_f^{\rm T})/2$ is the total fluid stresses, 
$\boldsymbol {u}_p$ is the Eulerian particle velocity, 
$\boldsymbol {u}_{f|p}$ is the fluid velocity at the particle location, and 
$f_d$ is an inertial drag correction. The first term in (2.3) represents the effect of the undisturbed fluid stresses at the particle locations. The second term in (2.3) represents the drag force exerted by the particle phase back onto the fluid phase. In this study, we use the inertial drag correction given by Tenneti, Garg & Subramaniam (Reference Tenneti, Garg and Subramaniam2011),
\begin{equation} f_d=\frac{1+0.15{Re}_p^{0.687}}{(1-\phi_p)^3}+\frac{5.81\phi_p}{(1-\phi_p)^3}+\frac{0.48\phi_p^{1/3}}{(1-\phi_p)^4}+\phi_p^3 {Re}_p \left(0.95+\frac{0.61\phi_p^3}{(1-\phi_p)^2}\right),\end{equation}
where 
${Re}_p$ is the particle Reynolds number. In the limit, where 
$\phi _p\ll 1$, this correlation collapses onto the well-known Schiller–Naumann drag correlation for a sphere.
For particle-laden flows with high density ratio (
$\rho _p/\rho _f\gg 1$), the interphase coupling is dominated by the drag force, whereas the effect of undisturbed fluid stresses becomes negligible. From a scaling analysis of the drag force in (2.3), it is clear that the mass loading 
$M=\rho _p\langle \phi _p\rangle /\rho _f$ determines the strength of the coupling. Thus, if the mass loading is vanishingly small, the disperse phase has no effect on the stability or decay of the Lamb–Oseen vortex. Conversely, if the mass loading is of order unity, the term (2.3) represents a significant source or sink of energy that may enhance or attenuate perturbations traversing the flow (Kasbaoui et al. Reference Kasbaoui, Koch, Subramanian and Desjardins2015).
The particle phase is treated in the Lagrangian reference frame. The evolution of the 
$i$th particle is given by the equation of motions for a spherical particle (Maxey & Riley Reference Maxey and Riley1983),
$$\begin{gather} \frac{{\rm d}\boldsymbol{x}_p^i}{{\rm d}t}(t)={\boldsymbol{u}_p^i}(t), \end{gather}$$
$$\begin{gather}\frac{{\rm d}\boldsymbol{u}_p^i}{{\rm d}t}(t)= \frac{1}{\rho_p}\boldsymbol{\nabla}\boldsymbol{\cdot} \boldsymbol{\tau}(\boldsymbol{x}_p^i,t)+ f_d\frac{\boldsymbol{u}_f(\boldsymbol{x}_p^i,t)-\boldsymbol{u}^i_p(t)}{\tau_p}, \end{gather}$$
where 
$\boldsymbol {x}_p^i$ and 
$\boldsymbol {u}_p^i$ are the position and velocity of the 
$i$th particle, respectively, 
$\tau _p=\rho _pd_p^2/(18\mu )$ is the particle response time and 
$d_p$ is the particle diameter. In order to focus on inertial effects, gravity is ignored. Owing to the large density ratio and low volume fraction in the semidilute regime, other hydrodynamic and particle–particle interactions may be neglected.
Note that, in the equations above, the instantaneous particle volume fraction and Eulerian particle velocity field are obtained from Lagrangian quantities using
$$\begin{gather} \phi_p(\boldsymbol{x},t) =\sum_{i=1}^N V_p g(\|\boldsymbol{x}-\boldsymbol{x}_p^i\|), \end{gather}$$
$$\begin{gather}\phi_p \boldsymbol{u}_p(\boldsymbol{x},t) =\sum_{i=1}^N \boldsymbol{u}_p^i(t) V_p g(\|\boldsymbol{x}-\boldsymbol{x}_p^i\|), \end{gather}$$
where 
$V_p={\rm \pi} d_p^3/6$ is the particle volume. In these expressions, 
$g$ represents a Gaussian filter kernel of size 
$\delta _f=3\Delta x$, where 
$\Delta x$ is the grid spacing. Further details on the computation of these terms, as well as other aspects of the solver, can be found in Capecelatro & Desjardins (Reference Capecelatro and Desjardins2013b).
3. Configuration and flow parameters
The impact of inertial particles on the Lamb–Oseen vortex is investigated with the configuration schematized in figure 1. Following the approach of Takahashi, Ishii & Miyazaki (Reference Takahashi, Ishii and Miyazaki2005), four counter-rotating vortices are initiated in a square domain with periodic boundary conditions applied in both directions. Each quadrant containing one of the four vortices has a size 
$30r_0$ by 
$30r_0$, where 
$r_0$ is the initial vortex core scale. The spacing between two neighbouring vortices is 
$30r_0$, such that vortex–vortex interactions are negligible for the length of the simulations considered, as shown in § 4. Monodisperse particles are seeded randomly in the domain with velocities equal to the fluid velocity at their locations. In the subsequent analysis, data collected from all four quadrants are used in the computation of statistics after application of appropriate symmetries.
Four counter-rotating vortices are arranged in the computational domain with sufficient spacing to allow for independent evolutions. Initially, the inertial particles are seeded randomly. Radial profiles are computed by averaging statistics from all four quadrants with appropriate sign changes.
Table 1 lists a summary of non-dimensional parameters for the nine cases considered in this study. In all these simulations, the carrier fluid has density 
$\rho _f=1.2\ \mathrm {kg\ m^{-3}}$ and viscosity 
$\mu _f = 1.8 \times 10^{-5}\ \mathrm {kg}\ {\rm m}^{-1}\ {\rm s}^{{-1}}$. The initial vortex core scale is 
$r_0=0.05\ \mathrm {m}$. Case SP is a reference unladen Lamb–Oseen vortex at the circulation Reynolds number 
$\textit {Re}_{\varGamma }=1700$, given by the vortex circulation 
$\varGamma = 1.57\ \mathrm {m}^{2}\ {\rm s}^{{-1}}$. In the remaining cases, monodisperse particles are seeded with density 
$\rho _p = 1000\ \mathrm {kg}\ \mathrm {m}^{-3}$, average volume fraction 
$\langle \phi _p\rangle =1.2\times 10^{-3}$ and mass loading 
$M=1$. In cases A–C, the Stokes number is varied while the Reynolds number is fixed at 
$\textit {Re}_{\varGamma }=1700$. The Stokes number is varied by changing the particle diameter to be 
$d_p = 57\ \mathrm {\mu } \mathrm {m}$ for case A, 
$d_p = 80\ \mathrm {\mu } \mathrm {m}$ for case B and 
$d_p = 113\ \mathrm {\mu } \mathrm {m}$ for case C. These values yield circulation-based Stokes numbers equal to 
$\textit {St}_{\varGamma }=\tau _p/\tau _f=0.1$, 0.2 and 0.4, respectively, where 
$\tau _f=2{\rm \pi} r_0^2/\varGamma$ is the characteristic vortex time scale. In cases D–F, the Reynolds number is varied from 
$800$ to 5000 by varying the vortex circulation 
$\varGamma$. For these cases, the Stokes number is maintained constant at 
$\textit {St}_{\varGamma }=0.4$ by changing the particle diameter accordingly.
Non-dimensional parameters considered in the simulations. Case SP corresponds to an unladen (single-phase) Lamb–Oseen vortex. Cases A–C correspond to particle-laden Lamb–Oseen vortices where the Stokes number is varied by varying the particle diameter. In cases D–F, the Reynolds number is varied by changing the vortex circulation, while the particle diameter is varied to maintain constant Stokes number. In all these runs, a Cartesian grid is employed with uniform resolution corresponding to about 102 grid points across the initial vortex core scale.
The parameters chosen here are such that collective effects of the particle phase dominate over discrete effects. Effects due to the discrete forcing exerted by the Lagrangian particles have a scale comparable to the interparticle distance 
$L_p$. When the particles are first seeded, this distance may be estimated as 
$L_p\sim (V_p/ \langle \phi _p\rangle )^{1/3}=({\rm \pi} /6)^{1/3}\langle \phi _p\rangle ^{-1/3}d_p$. For a semidilute flow with 
$\langle \phi _p\rangle \sim 10^{-3}$, the interparticle distance is thus comparable to the particle diameter, i.e. 
$L_p\sim d_p$. Kasbaoui et al. (Reference Kasbaoui, Koch, Subramanian and Desjardins2015) argue that, in shear-dominated flows, particle-induced fluctuations on scales smaller than the viscous dissipation scale are quickly suppressed by viscosity at a super-exponential rate. For the Lamb–Oseen vortex, the vortex core scale represents the viscous scale of the flow. In all cases chosen, the ratio 
$r_0/d_p$ is large (see table 1), such that fluctuations due to the discrete particle forcing fall well below the viscous dissipation scale. In this way, the coupling between the fluid and particles occurs primarily through collective particle dynamics at scales comparable with 
$r_0$, rather than discrete effects at scales comparable with the interparticle distance.
Owing to their finite inertia, the initial particle velocity may affect the response of the fluid. Particles that are seeded with velocities far away from their equilibrium velocity must be accelerated or decelerated, which results in a premature modification of the carrier flow. To avoid such effects, the particles are initialized with velocities that match the fluid velocity at the particle locations, i.e. with zero slip velocity. The reader is referred to Appendix B for a discussion of additional simulations conducted with non-zero initial slip velocity.
The equations of motion are integrated using the finite-volume flow solver NGA (Desjardins et al. Reference Desjardins, Blanquart, Balarac and Pitsch2008) with the Eulerian–Lagrangian strategy of Capecelatro & Desjardins (Reference Capecelatro and Desjardins2013b) on a Cartesian mesh. This computational infrastructure has been applied and validated in several prior studies from dilute to dense particle-laden flows (Capecelatro & Desjardins Reference Capecelatro and Desjardins2013a; Arolla & Desjardins Reference Arolla and Desjardins2015; Kasbaoui Reference Kasbaoui2019; Kasbaoui et al. Reference Kasbaoui, Koch and Desjardins2019b,Reference Kasbaoui, Koch and Desjardinsa; Rao & Capecelatro Reference Rao and Capecelatro2019; Wang et al. Reference Wang, Fong, Coletti, Capecelatro and Richter2019). For further details on how this method compares with the two-fluid method and other kinetic-based Eulerian–Eulerian methods, the reader is referred to Patel et al. (Reference Patel, Desjardins, Kong, Capecelatro and Fox2017).
In order to provide good resolution of the vortex core, we use the spatial resolution 
$3072\times 3072$, which yields the ratio 
$r_0/\Delta x\simeq 51$, i.e. 102 grid points lie across the initial vortex core scale. This resolution yields adequate grid convergence as is shown in Appendix A. It is also noteworthy that this resolution is significantly higher than what is commonly employed in studies of single-phase vortex dynamics. For example, Orlandi (Reference Orlandi2007) found that 30 grid points across the core (
$r_0/\Delta x= 15$) were sufficient to capture the viscous and convective dynamics of merging vortex pairs. Despite the relatively low volume fraction considered in this study, a great number of Lagrangian particles must be tracked owing to the large scale separation between the particle scales at 
$O(d_p)$ and the modulation scales at 
$O(r_0)$. The total number of particles tracked, 
$N$, is 
$2.7\times 10^7$, 
$9.6 \times 10^6$ and 
$3.4 \times 10^6$ for cases A, B and C, respectively. For cases D, E and F, 
$N$ is 
$1.2\times 10^6$, 
$9.7\times 10^6$ and 
$1.8\times 10^7$, respectively.
For the present simulations of particle-laden vortices, we shall note that the axial direction has a finite thickness 
$L_z=3d_p$, is resolved with one grid point (
$n_z=1$) and is taken to be periodic. This numerical set-up, referred to as ‘pseudo-2-D’, serves to represent the spherical geometry of the particles and compute quantities such as volume fraction and number density. Further details can be found in Kasbaoui et al. (Reference Kasbaoui, Koch and Desjardins2019a), where the authors study the onset of particle clustering in ‘pseudo-2-D’ simulations of semidilute particle-laden linear shear flow.
4. Fluid motion of the single-phase Lamb–Oseen vortex
In this section, we present numerical results pertaining to the single-phase flow Lamb–Oseen vortex. The goal is twofold. First, this section intends to be a short review of the main characteristics of a tubular Lamb–Oseen vortex in a single-phase flow, i.e. prior to adding particles. Second, as a full solution of the Navier–Stokes equations, the theoretical expressions of a Lamb–Oseen vortex offer the opportunity to validate our numerical approach.
Figure 2 shows the isocontours of the velocity and vorticity field from the single-phase flow simulation (case SP) at 
$\textit {Re}_{\varGamma }=1700$. The snapshots are given between 
$t/\tau _f=0$ and 30 in a view centred on one of the four vortices in the domain. Here, 
$\tau _f=2{\rm \pi} r_0^2/\varGamma$ represents the time it takes to complete a full rotation for a fluid particle located at a radial distance where the azimuthal velocity reaches its peak. This time scale is indicative of inertial effects associated with a Lamb–Oseen vortex. As shown in figure 2, the velocity field has an annular structure whereby it is made of concentric circles. The azimuthal velocity reaches a peak velocity at the location 
$r_c$, which denotes the core radius. The vorticity has a Gaussian distribution. Viscous effects induce the dampening of the vorticity at the core and the growth of the core radius 
$r_c$ as the vortex tube spreads outwards. This process occurs at the slow viscous time scale 
$\tau _\mu =\rho _f r_0^2/\mu$. In figure 2, no significant viscous effects can be seen by the naked eye within the time range 
$t/\tau _f=0$ to 30. This is due to the large Reynolds number considered, which relates to the inertial and viscous time scales according to 
$\textit {Re}_{\varGamma }=\tau _\mu /\tau _f$. Further, as expected from the stable nature of the Lamb–Oseen vortex, the flow retains its cylindrical symmetry throughout the simulation.
Isocontours of the normalized velocity magnitude 
$\|\boldsymbol {u}_f\|/(\varGamma /2{\rm \pi} r_0)$ (first row) and axial vorticity 
$\omega _z/\omega _{c,0}$ (second row) at successive instants for the single-phase Lamb–Oseen vortex at 
$\textit {Re}_{\varGamma }=1700$.
We now show that the numerical approach is capable of reproducing the evolution of a Lamb–Oseen vortex. We compare the simulation data with the analytical expressions for vortex core radius, vorticity at the centre and total enstrophy. These expressions are derived from the analytical form given in (1.1) and read
$$\begin{gather} r_c(t) ={\rm arg\,max}_r u_{f,\theta}= \sqrt{\alpha (r_0^2+4 \nu t)} , \end{gather}$$
$$\begin{gather}\omega_c(t) = \omega_z(r=0,t)=\frac{\varGamma}{{\rm \pi} r_0^2+4{\rm \pi} \nu t}, \end{gather}$$
$$\begin{gather}\epsilon(t) = \iint_{\mathbb{R}^2} \omega_z^2 \,{\rm d}S =\frac{\varGamma^{2}}{2{\rm \pi} r_c^2}, \end{gather}$$
where 
$\alpha =1.2564$ is the Oseen parameter.
Figure 3 shows a comparison of the numerical results and analytical expressions in (4.1), (4.2) and (4.3). The figure shows results from case SP at 
$\textit {Re}_{\varGamma }=1700$ and two auxiliary simulations at 
$\textit {Re}_{\varGamma }=830$ and 560 obtained by changing the viscosity 
$\mu$ while keeping 
$\varGamma$ and 
$r_0$ identical. Time is normalized by the vortex characteristic time scale 
$\tau _f$, which is the same in all three runs. We observe that larger dynamic viscosity leads to faster decay of the vortex characterized by growing core radius and decaying centre vorticity. Enstrophy, which quantifies the vortical energy contained in this tubular vortex, decays at a rate that approaches 
$(t/\tau _\mu )^{-1}$. In all of figures 3(a–c), we observe that the numerical results match closely with the analytical expressions for the entire duration of simulation considered (
$0\leqslant t/\tau _f \leqslant 120$). This excellent agreement with the expressions for an isolated vortex tube validates our numerical approach and confirms that the interaction between the four vortices in the domain is negligible during this time frame.
Evolution of the normalized (a) vortex core radius, (b) centre vorticity and (c) enstrophy for an unladen Lamb–Oseen vortex at 
$\textit {Re}_{\varGamma }=1700$ (black 
$\circ$ /solid), 
$\textit {Re}_{\varGamma }=830$ (red 
$\square$ /dashed) and 
$\textit {Re}_{\varGamma }=560$ (blue 
$\diamond$ /dash-dotted). Symbols correspond to simulation data. Lines correspond to analytical expressions.
5. Lamb–Oseen vortex with suspended particles
5.1. Overview
The Lamb–Oseen vortex from case SP is now laden with inertial particles. We investigate the effects of particle inertia on this prototypical vortex tube under strong interphase coupling due to mass loading 
$M=1$. We note that in cases A, B and C, the Stokes number 
$\textit {St}_{\varGamma }$ is varied by varying the particle diameter without changing the parameters of the carrier Lamb–Oseen flow.
Figure 4 shows the isocontours of particle volume fraction for the three Stokes numbers 
$\textit {St}_{\varGamma }=0.1$, 0.2 and 0.4 at successive instants. Time is non-dimensionalized by the vortex characteristic time scale 
$\tau _f$ and varies between 0 and 12. Similar to figure 2, the view is centred on one of the four vortices in the domain. In all cases, it is seen that the initial random particle distribution becomes progressively highly inhomogeneous. As particles are expelled from the vortex core, a void fraction bubble forms. This region is eventually fully depleted of all particles. The latter cluster into a ring structure surrounding the void fraction bubble and vortex core. The clustering pattern is similar to those reported by Druzhinin (Reference Druzhinin1994) in two dimensions, and analogous to the annular clusters in the 3-D simulations of Marshall (Reference Marshall2005).
Isocontours of particle volume fraction normalized by the average volume fraction for cases A (top), B (middle) and C (bottom) at five non-dimensional times 
$t/\tau _f$. As time progresses, the inertial particles are expelled out of the vortex cores, forming a ring cluster around a void fraction bubble. The rate at which these structures form increases with particle inertia.
The evolution of the azimuthally averaged volume fraction distribution for case A is shown in figure 5. The particle volume fraction within the initial vortex core region (
$r< r_0$) decreases with time until fully depleted, consistent with the formation of a void fraction bubble. The particles start clustering at a radial distance 
${\sim }1.5r_0$, which then shifts outwards to 
${\sim }4.5r_0$. The expansion of the particle cluster slows down with time, showing the transient nature of this process. We also note that the volume fraction within the cluster caps at approximately 1.6 times higher than the average volume fraction. Thus, these clusters remain dilute, with negligible chance for particle–particle interaction.
Evolution of the azimuthally averaged volume fraction for a particle-laden Lamb–Oseen vortex with 
$\textit {St}_{\varGamma }=0.4$. Owing to their inertia, the particles are progressively expelled from the vortex core, leading to the formation of a ring particle cluster.
In figure 6(a), we show the time evolution of the radius of the void fraction bubble 
$r_b$, normalized by the initial vortex core scale 
$r_0$, for the three cases at 
$\textit {St}_{\varGamma }=0.1$, 0.2 and 0.4, and 
$\textit {Re}_{\varGamma }=1700$. Here, we define 
$r_b$ as the location of the peak particle volume fraction. We note that increasing particle inertia results in faster formation of a void bubble due to the faster outward migration of higher inertia particles. In figure 6(b), we fit the evolution of 
$r_b$ using least-squares power-law regression, 
$C_b(t/\tau _f)^m+r_{b,0}/r_0$. The coefficient 
$r_{b,0}$ represents the location where the particles start accumulating. In figure 6(b), 
$r_{b,0}/r_0\sim 1.4$ in agreement with naked-eye observations from figures 4 and 5. Note that the particle accumulation for 
$\textit {St}_{\varGamma }=0.1 - 0.4$ starts at approximately the same location since 
$r_{b,0}$ does not vary considerably. Similarly, the power-law exponent 
$m\sim 0.5$ does not vary considerably among the three cases A–C, and indicates a long-term evolution that scales with 
$\sqrt {t}$. Based on these fits, we see that the speed at which the ring cluster expands, captured by 
${\rm d}r_b/{\rm d}t$, is largest for 
$t\ll 1$ and drops as 
$1/\sqrt {t}$ for 
$t\gg 1$. The power-law coefficient 
$C_b$ increases with Stokes number, which denotes a larger growth of the bubble radius with 
$\textit {St}_{\varGamma }$.
Normalized radius of the void fraction bubble for 
$\textit {St}_{\varGamma }=0.1$ (black 
$\circ$), 0.2 (red 
$\square$) and 0.4 (blue 
$\diamond$). Increasing particle inertia leads to faster formation of a void bubble.
In figure 7, we compare the isocontours of the axial vorticity from the single-phase and particle-laden Lamb–Oseen vortices in cases A, B and C. Contrary to the single-phase case, where the vorticity field retains a cylindrical symmetry at all times, the presence of inertial particles induces a quick breakdown of the vortex flow structure and the loss of symmetry. Owing to the modulating effect of the disperse phase, the fluid vorticity at the vortex core decays faster than that of the single-phase flow. The rate of decay increases with increasing particle inertia. Moreover, the smearing of the vortex over time is further enhanced with increasing Stokes number.
Isocontours of normalized vorticity (
$\omega _z/\omega _c$) for single-phase and different 
$\textit {St}_{\varGamma }$ numbers at five non-dimensional times 
$t/\tau _f$.
5.2. Time scale of particle-laden vortex and preferential concentration
As shown in § 5.1, the dynamics of the particle-laden Lamb–Oseen vortex depend greatly on particle inertia. It is clear that, besides the vortex characteristic time scale 
$\tau _f$, there is another time scale related to the particle phase that controls the evolution of this particle-laden vortex. Here, we show that the appropriate time scale is given by the clustering time scale 
$\tau _c=\tau _f/\textit {St}_{\varGamma }$ imposed by the preferential concentration mechanism.
To understand why a particle migrates away from the vortex cores, we consider the continuity equation for the particle phase given by
\begin{equation} \frac{\partial n}{\partial t}+\boldsymbol{\nabla}\boldsymbol{\cdot} (n\boldsymbol{u}_p) =0, \end{equation}
where 
$n$ represents the particle number density and 
$\boldsymbol {u}_p$ the particle velocity field. This equation can be rearranged as
\begin{equation} \frac{\partial n}{\partial t} +\boldsymbol{u}_p\boldsymbol{\cdot}\boldsymbol{\nabla}n={-}n\boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{u}_p , \end{equation}
where it becomes clear that the right-hand side of (5.2) controls the accumulation or depletion of particles. The reason is, unlike the divergence-free fluid phase, the particle phase is highly compressible. Following Maxey (Reference Maxey1987) and Ferry & Balachandar (Reference Ferry and Balachandar2001), the right-hand side of (5.2) can be approximated as 
$-n\boldsymbol {\nabla }\boldsymbol {\cdot } \boldsymbol {u}_p \simeq \tau _p n \boldsymbol {\nabla } \boldsymbol {u}_f\boldsymbol {:}\boldsymbol {\nabla } \boldsymbol {u}_f$ in the limit of small Stokes numbers 
$St \ll 1$. This term can be further decomposed as 
$\tau _p n \boldsymbol {\nabla } \boldsymbol {u}_f\boldsymbol {:}\boldsymbol {\nabla } \boldsymbol {u}_f=\tau _p n (\textit {\mathsf {S}}^2-\textit {\mathsf {R}}^2)$, where 
$\textit {\mathsf {S}}^2$ and 
$\textit {\mathsf {R}}^2$ represent the second invariants of the fluid strain-rate tensor 
$\boldsymbol {\textit {\mathsf {S}}}=(\boldsymbol {\nabla } \boldsymbol {u} +\boldsymbol {\nabla } \boldsymbol {u}^{\rm T})/2$ and rotation tensor 
$\boldsymbol {\textit {\mathsf {R}}}=(\boldsymbol {\nabla }\boldsymbol {u} -\boldsymbol {\nabla } \boldsymbol {u}^{\rm T})/2$ (Squires & Eaton Reference Squires and Eaton2014). When the local strain exceeds the local rotation, this term is positive, leading to particle accumulation. Conversely, when this term is negative, the particles are expelled. This phenomenon is known as preferential concentration. For a particle-laden Lamb–Oseen vortex, the preferential mechanism causes inertial particles to be expelled from the highly rotational vortex core, resulting in the generation of ring-shaped particle cluster surrounding the void bubble.
Since the time scale of the particle expulsion from the vortex cores is controlled by (5.2), we analyse the scaling of each term,
\begin{equation} \underbrace{\vphantom{\frac{1}{t}}\frac{\partial n}{\partial t}}_{n_0/\tau_c} +\underbrace{\vphantom{\frac{1}{t}}\boldsymbol{u}_p\boldsymbol{\cdot}\boldsymbol{\nabla} n}_{u_{p,0} n_0/r_0} =\underbrace{\vphantom{\frac{1}{t}}\tau_p n (\textit{\mathsf{S}}^2-\textit{\mathsf{R}}^2)}_{n_0\tau_p u_{f,0}^2/{r^2_0}}, \end{equation}
where 
$\tau _c$ is the characteristic clustering time scale, 
$n_0=\langle \phi _p\rangle /({\rm \pi} d_p^3/6)$ is the average number density, 
$u_{f,0}=\varGamma /2{\rm \pi} r_0$ is the characteristic fluid velocity and 
$u_{p,0}$ is the characteristic particle velocity. By balancing the rate of change of number density with the preferential concentration term, i.e. the right-hand side of (5.3), we see that the clustering time scale is given by 
$\tau _{c}={r^2_0}/(\tau _p u_{f,0}^2)=\tau _f/\textit {St}_{\varGamma }$.
5.3. Formation of void fraction bubble and particle clustering
The simulation results support the argument that the time scale 
$\tau _{c}=\tau _f/\textit {St}_{\varGamma }$ controls the large-scale dynamics of the particle phase. To illustrate this aspect, we show, in figure 8, the isocontours of the normalized volume fraction field for cases A, B and C at 
$t/\tau _c=0$, 3, 6, 9 and 12. When considered at the same non-dimensional time 
$t/\tau _c$, the particle distributions present similar large-scale patterns regardless of the particle Stokes number. Specifically, it can be visually observed that the void fraction bubbles have similar sizes, which shows that the evolution of these structures is dominated by preferential concentration.
Isocontours of normalized particle volume fraction for three cases A, B and C (top to bottom) at non-dimensional times 
$t/\tau _c=$0, 3, 6, 9 and 12. The similar evolution of the void fraction bubble shows that the clustering time scale 
$\tau _c$ controls the large-scale evolution of the particle phase.
We should note that there is additional small-scale clustering seen in figure 8 that does not evolve on the time scale 
$\tau _c$. These small structures manifest as void fraction regions of size 
$\ll r_b$. The small-scale clustering is weak for case 
$\textit {St}_{\varGamma }=0.1$, and strongest for 
$\textit {St}_{\varGamma }=0.4$. The origin of these structures could be tied to the formation of caustics (Wilkinson & Mehlig Reference Wilkinson and Mehlig2005; Gustavsson et al. Reference Gustavsson, Meneguz, Reeks and Mehlig2012; Ravichandran & Govindarajan Reference Ravichandran and Govindarajan2015) resulting from particle trajectory crossing (Kasbaoui et al. Reference Kasbaoui, Koch and Desjardins2019a). These effects arise from the discrete nature of the particle phase, and are generally stronger with increasing Stokes number.
The radial profiles of the azimuthally averaged normalized volume fraction are shown in figure 9. The averaging procedure reduces the effect of small-scale clustering, allowing us to focus on the large scales. In figure 9(a), the radial profiles are shown at non-dimensional times between 0 and 
$64$ where time is normalized by the vortex characteristic time 
$\tau _f$. Owing to the different rates at which inertial particles are expelled out of the vortex cores, the radial profiles of volume fraction are further shifted outwards with increasing Stokes number. These profiles collapse narrowly when reported at the same non-dimensional times normalized by 
$\tau _c$ as seen in figure 9(b). Nevertheless, small differences in the particle distribution increase with time. At 
$t/\tau _c=12$, the particle distribution is slightly narrower with increasing Stokes number, resulting in higher peak volume fraction within the ring cluster. This effect is due to nonlinear wave steepening similar to that observed in linear shear by Kasbaoui et al. (Reference Kasbaoui, Koch and Desjardins2019a). Trajectories of the peak volume fraction plotted with lines in figure 9(a) show that the volume fraction within the ring cluster saturates. As argued in Kasbaoui et al. (Reference Kasbaoui, Koch and Desjardins2019a), this saturation is a result of particle trajectory crossing in semidilute flows. The effect of this mechanism is to smooth the nonlinear wave steepening, which limits the ability to concentrate volumetric regions of inertial particles into one-dimensional (lines) or 2-D (sheets) structures.
Profiles of azimuthally averaged particle volume fraction for 
$\textit {St}_{\varGamma }=0.1$ (black 
$\circ$), 0.2 (red 
$\square$) and 0.4 (blue 
$\diamond$) at (a) 
$t/\tau _f=0$, 8 and 64, and (b) 
$t/\tau _c= 0$, 3 and 12. The thick dashed lines in (b) represent trajectories of peak volume fraction. The collapse of the profiles in (b) shows that the clustering time scale 
$\tau _c$ controls the evolution of the ring cluster.
5.4. Modulation of the fluid phase
We now focus on the modulation induced by the inertial particles on the carrier phase. We discuss the change to the structure of the flow field in terms of the size of the vortex core, vorticity distribution and enstrophy.
Figure 10(a) shows the evolution of the vortex core radius 
$r_c$ with the non-dimensional time 
$t/\tau _c=St\times t/\tau _f$ for cases A–C. We observe that the vortex core grows at a rate that is significantly faster than in the unladen case. As seen from (4.1) and figure 3(a), viscosity controls the evolution of 
$r_c$ in the single-phase case, where the structure of the Lamb–Oseen vortex is made out of successive annuli of fluid at different azimuthal velocities. Momentum diffusion across these layers causes the vortex core to grow at the rate 
$\sqrt {t/\tau _\mu }$. In the two-phase setting, the interaction between particles and fluid phase dominates over momentum diffusion. The outward flux of particles induced by preferential concentration imposes the clustering time scale 
$\tau _c$ on the carrier phase.
Evolution of the vortex core radius with non-dimensional time 
$t/\tau _c$. (a) Cases A–C: 
$\textit {Re}_{\varGamma }=17\,000$ and 
$\textit {St}_{\varGamma }=0.1$ (black 
$\circ$), 0.2 (red 
$\square$) and 0.4 (blue 
$\diamond$). (b) Cases C–F: 
$\textit {St}_{\varGamma }=0.4$ and 
$\textit {Re}_{\varGamma }=800 - 5000$. The solid line in (b) represents equation (5.5).
As shown in figure 10(a), the initial evolution of 
$r_c$ nearly collapses for all three cases with 
$\textit {St}_{\varGamma }=0.1$, 0.2 and 0.4 particles. This trend holds for as long as there remain particles within the core, until about 
$t/\tau _c\sim 3$. The data in figure 10(a) suggest that a power-law fit,
\begin{equation} \frac{r_c^2(t)}{\alpha r_0^2}-1=\beta \left(\frac{t}{\tau_c}\right)^m \quad \text{for}\ t/\tau_c<3, \end{equation}
can capture the early growth of the vortex core (
$t/\tau _c<3$). The 
$\tau _c$ scaling in (5.4) embeds the idea that the growth of the core is driven by preferential concentration. For all three Stokes numbers considered, the power-law exponent 
$m$ is approximately equal to unity (see figure 10a). The coefficient 
$\beta$ has a weak non-monotonic variation with 
$\textit {St}_{\varGamma }$ around 
${\sim }0.5$.
Note that, in the case of vanishingly small particle inertia (
$\textit {St}_{\varGamma }\rightarrow 0$), the core radius evolution is similar to the single-phase case, i.e. driven by viscosity such that 
$r_c^2(t)/(\alpha r_0^2)-1= 4t/\tau _\mu$. Thus, a relationship that takes into account both viscous and preferential concentration effects for 
$\textit {St}_{\varGamma }=O(1)$ is
\begin{equation} \frac{r_c^2(t)}{\alpha r_0^2}-1=\frac{4t}{\tau_\mu} + \beta \frac{t}{\tau_c} = \frac{4t}{\tau_\mu}\left(1+\frac{\beta}{4}\textit{St}_{\varGamma}\textit{Re}_{\varGamma}\right) \quad \text{for}\ t/\tau_c<3, \end{equation}
where the relationships 
$\textit {Re}_{\varGamma }=\tau _\mu /\tau _f$ and 
$\textit {St}_{\varGamma }=\tau _f/\tau _c$ were used. In (5.5), the viscous growth is recovered when 
$\textit {St}_{\varGamma }\rightarrow 0$. For simulations cases A–C (
$\textit {St}_{\varGamma }=0.1 - 0.4$ and 
$\textit {Re}_{\varGamma }=1700$), growth due to preferential concentration dominates over the viscous growth since 
$(\beta /4)\textit {St}_{\varGamma }\textit {Re}_{\varGamma }\gg 1$.
In addition to the Stokes number 
$\textit {St}_{\varGamma }$, (5.5) shows a striking dependence on the Reynolds number 
$\textit {Re}_{\varGamma }$. The latter has no effect on 
$r_c$ in single-phase vortex tubes, which is evidence of the resilience of these structures. Here, it is suggested that dispersed inertial particles help break down these flow structures at high Reynolds numbers by enhancing the spreading of the vortex tube. To show the validity of the correlation (5.5), we report the evolution of the core radius from cases C–F in figure 10(b). In these runs, the Reynolds number is varied from 800 to 5000 while the Stokes number is held constant at 
$\textit {St}_{\varGamma }=0.4$. For all these cases, the data collapse onto the solid black lines representing (5.5).
The expulsion of the particles away from the core eventually decouples the evolution of the core radius 
$r_c$ from the particle phase. Beyond 
$t/\tau _c\sim 3$, the core is fully depleted. The inertial particles accumulate at the radial distance 
$r_b$ corresponding to the void fraction bubble radius. Note that the correlations in figure 6(b) yield coefficients that are independent of 
$\textit {St}_{\varGamma }$ when scaled with 
$\tau _c$, namely,
\begin{equation} {r_b}/{r_0}\simeq 0.4 \sqrt{t/\tau_c}+1.4 . \end{equation}
The void fraction bubble continues growing at a rate faster than that of the vortex core 
$r_c$, as shown by the growing ratio 
$r_b/r_c$ in figure 11. The decoupling of the two phases results in the pullback of 
$r_c$ beyond 
$t/\tau _c\sim 3$. The absence of particles in the core region causes the flow to rearrange itself in a way more consistent with that of a single-phase vortex tube.
Ratio between void bubble radial location 
$r_b$ and fluid core radius 
$r_c$ for 
$\textit {St}_{\varGamma }=0.1$ (black 
$\circ$), 0.2 (red 
$\square$) and 0.4 (blue 
$\diamond$).
It is noteworthy that Druzhinin (Reference Druzhinin1995) found a similar expression to (5.6) for the growth of the void bubble. Assuming weakly inertial particles (
$\textit {St}_{\varGamma }\ll 1$) and employing the method of characteristics, Druzhinin (Reference Druzhinin1995) provides
\begin{equation} \ln\left(\frac{r_b}{r_{b,0}}\right)+\frac{r_b^2}{r_{b,0}^2} +\frac{r_b^4}{4r_{b,0}^4}=\frac{1}{2}\frac{t}{\tau_c}+\frac{5}{4} , \end{equation}
where 
$r_{b,0}=\sqrt {2}r_0$ is the starting position of the volume fraction crest. Expression (5.7) captures the same dependence on the clustering time scale 
$\tau _c$. Further, at sufficiently short times 
$t/\tau _c$, the logarithmic and quartic terms are negligible, making expression (5.7) in qualitative agreement with (5.6).
The interaction between particles and fluid results in a weakened vortex characterized by wider vorticity distribution and diminished magnitude. Figure 12 shows the evolution of the isocontours of axial vorticity with non-dimensional time 
$t/\tau _c$. We observe that the large-scale features of the vorticity field, such as the spreading of the vortical core, unfold at the same time scale 
$\tau _c$ for all three Stokes number cases. By comparing figure 12 with figure 8, we also observe that a second vorticity peak is generated at the location of the concentration peak due to the preferential accumulation of particles. A similar vorticity peak is also noted in Druzhinin (Reference Druzhinin1995).
Isocontours of vorticity for the three cases A, B and C (top to bottom) at non-dimensional times 
$t/\tau _c=0$, 3, 6, 9 and 12.
The spreading of the vortical core is well visualized in figure 13, which shows radial profiles of the azimuthally averaged fluid vorticity normalized by the initial core vorticity 
$\omega _{c,0}$. The tails of the vorticity distributions (
$r>1.5r_0$) for different Stokes number cases collapse when the time scale 
$\tau _c$ is employed for the entire range 
$t/\tau _c=0 - 12$. However, the vorticity in the core region displays a more complicated time dependence. During the early stages 
$0< t\lesssim 4\tau _c$, the vorticity in the core regions drops considerably due to the outward forcing exerted by the ejection of the inertial particles from the core. The decay of the core vorticity ceases once all particles have been ejected from the core. The vorticity at the centreline stabilizes at about 
$0.50\omega _{c,0}$, 
$0.56\omega _{c,0}$ and 
$0.60\omega _{c,0}$, for 
$\textit {St}_{\varGamma }=0.1$, 0.2 and 0.4, respectively. Once the core is fully depleted, slow viscous dissipation becomes the dominant mechanism for vorticity decay in the vortex core. The numerical results of vorticity reduction for 
$\textit {St}_{\varGamma }=0.1$ and 
$M=1$ at times 
$t/\tau _c>3$ are in excellent quantitative agreement with the analytical solution obtained in Druzhinin (Reference Druzhinin1995):
\begin{equation} \frac{\omega_c(t)-\omega_{c,0}}{\omega_{c,0}}={-}\frac{M}{1+M}={-}0.5. \end{equation}
Although not pursued in this study, expression (5.8) suggests that the decay of the centre vorticity increases with mass loading. The agreement between the present simulations and expression (5.8) deteriorates with increasing 
$\textit {St}_{\varGamma }$ due to the finite accuracy of the asymptotic assumption in Druzhinin (Reference Druzhinin1995).
Radial profiles of fluid vorticity normalized by the initial core vorticity for 
$\textit {St}_{\varGamma }=0.1$ (black solid), 0.2 (red dashed) and 0.4 (blue dash-dotted). (a–d) Profiles at various 
$t/\tau _f$ instants, showing considerable variation of the profile width with Stokes number. (e–h) Profiles at various 
$t/\tau _c$ instants, showing better collapse of the tails of the vorticity distribution across 
$\textit {St}_{\varGamma }$.
Besides preferential concentration, small-scale clustering induced by particle trajectory crossing further accelerates the decay of particle-laden vortex tubes. Figure 12 shows the presence of small-scale perturbations with characteristic size much smaller than the core and bubble radii. These localized vorticity gradients are stronger with increasing Stokes number and correlate with the small-scale clustering observed in figure 8. The effect of these clusters on the carrier flow is captured by the flow enstrophy,
\begin{equation} \epsilon(t)=\iint |\boldsymbol{\omega}|^2\,{\rm d}^2\boldsymbol{x}. \end{equation}This quantity is sensitive to the sharp vorticity gradients resulting from the small-scale clustering, hence, providing an indirect measure of particle–fluid interactions at those scales.
In figure 14, we report the enstrophy, normalized by its initial value 
$\epsilon _0$, for the three different 
$\textit {St}_{\varGamma }$ numbers and in the single-phase case. During the ejection phase (
$0\leqslant t\leqslant 3\tau _c$), enstrophy grows considerably. The amplitude of the growth increases with increasing Stokes number. This behaviour plays a significant role in the accelerated decay of the particle-laden Lamb–Oseen vortex, since enstrophy controls the rate at which the total kinetic energy 
$E(t)=\iint (|\boldsymbol {u}_f|^2/2)\,{\rm d}^2\boldsymbol {x}$ is dissipated. This is clearly seen from the equation of conservation of energy,
\begin{equation} \frac{{\rm d}E}{{\rm d}t}={-}\nu \epsilon+ \iint \boldsymbol{u}_f\boldsymbol{\cdot} \boldsymbol{F}_p\,{\rm d}^2\boldsymbol{x}, \end{equation}where a 2-D flow is assumed. One can see that inertial particles may accelerate the decay by two mechanisms: (i) through the modulation of enstrophy, and (ii) negative production due to the drag force (second term in (5.10)).
Decaying of the enstrophy at the vortex centre (normalized by its initial value) for various Stokes number. (a) Enstrophy decays for 
$\textit {St}_{\varGamma }=0.1$ (black 
$\circ$), 0.2 (red 
$\square$) and 0.4 (blue 
$\diamond$) based on time scale 
$\tau _f$. (b) Enstrophy decays with a non-dimensional time based on time scale 
$\tau _f$ in log–log scale.
Furthermore, the growth of enstrophy observed in figure 14 is a remarkable feature of the two-way coupled particle-laden vortex, in contrast with single-phase 2-D vortices whose enstrophy decays monotonically. Similar to (5.10), one may also derive an equation for enstrophy for a 2-D semidilute particle-laden vortex,
\begin{equation} \frac{{\rm d}\epsilon}{{\rm d}t}={-}\nu \iint \boldsymbol{\nabla} \boldsymbol{\omega}\boldsymbol{:}\boldsymbol{\nabla} \boldsymbol{\omega}\,{\rm d}^2\boldsymbol{x}+ \iint \boldsymbol{\omega} \boldsymbol{\cdot}\boldsymbol{\nabla}\times\boldsymbol{F}_p\, {\rm d}^2\boldsymbol{x} . \end{equation}
Given that the viscous term in (5.11) is strictly negative, it is the momentum exchange (second term on the right-hand side of (5.11)) that drives the growth of 
$\epsilon$. The 
$\tau _c$ scaling does not collapse the different Stokes number curves in figure 14(a), suggesting that it is the small-scale clustering rather than preferential concentration that is the main momentum exchange mechanism driving the growth of 
$\epsilon$. Once all particles have been driven out of the vortex core (
$t\gtrsim 3\tau _c$), enstrophy decays gradually. The power-law fits in figure 14(b) show that the decay of enstrophy is faster with increasing Stokes number. This is expected since the viscous term in (5.11) depends on the magnitude of the local vorticity gradients, which are larger with increasing particle inertia.
6. Conclusion
In this work, we investigate the flow modulation induced by inertial particles suspended in a prototypical vortex tube. The base flow is a 2-D Lamb–Oseen vortex seeded with randomly dispersed inertial particles. Six cases are considered where the Stokes number 
$\textit {St}_{\varGamma }$ is varied in the range 0.1–0.4 and the Reynolds number 
$\textit {Re}_{\varGamma }$ is varied in the range 800–5000. The mass loading 
$M=\rho _p/\rho _f\langle \phi \rangle$ being equal to unity ensures that the two phases are tightly coupled. The simulations show that inertial particles accelerate the decay of the vortex tube by orders of magnitude compared to an unladen vortex tube. The cause of the faster decay is strong interphase momentum exchange, whereby particles ejected from the vortex core cause an outward dispersion of the fluid momentum and enhanced dissipation due to small-scale clustering. Within the range of Stokes numbers considered, 
$0.1\leqslant \textit {St}_{\varGamma } \leqslant 0.4$, increasing particle inertia results in faster decay of the vortex and stronger inhomogeneity in the particle phase.
The interaction of inertial particles with a Lamb–Oseen vortex can be decomposed into two stages: an inertial clustering stage 
$0\leqslant t \lesssim 3\tau _c$, followed by a viscous stage 
$t >3\tau _c$. During the first stage, preferential concentration is the main driver of clustering. Owing to this mechanism, the particles, which are initially randomly distributed, are expelled from the vortex centre. This process causes the growth of a void fraction bubble until the core is fully depleted. The inertial particles cluster around the void fraction bubble into a ring of high particle volume fraction. This void bubble and ring have faster formation speed for particles with larger inertia and grow continuously at a rate that scales as 
$t^{-1/2}$. Within the concentrated clusters, the particle volume fraction increases with time initially until it saturates at 
${\sim }1.6$ times the average volume fraction for 
$\textit {St}_{\varGamma }=0.4$. The clustering time scale imposed by preferential concentration, 
$\tau _c=\tau _f/\textit {St}_{\varGamma }$, controls the large-scale dynamics. This results in a particle distribution displaying similar patterns for all 
$\textit {St}_{\varGamma }$ cases when considered at the same 
$t/\tau _c$. We also noticed that the growth of the void fraction bubble drives the growth rate of the vortex core radius 
$r_c$ to be much larger than that in the single-phase case. The axial vorticity, which initially has a Gaussian shape, is progressively flattened as the particles cause the expansion of the vortex. The dynamics of the fluid phase follow the clustering time scale 
$\tau _c$, until 
$t\sim 3\tau _c$. At the end of the clustering stage, the vortex core is depleted from all particles, and a viscous stage flows. The decay of the vortex significantly slows down due to an effective decoupling of the particle phase and the rotational fluid core. The evolution of the rotational core becomes similar to that of a single-phase vortex.
The simulations also show the presence of quick breakdown of the vortex flow symmetry and small-scale clusters, which have not been considered in previous studies. These inhomogeneities in the particle distribution have a scale significantly smaller than the void fraction bubble or fluid core radius. These small-scale clusters are likely to be due to particle trajectory crossing, which occurs when the local particle velocity field becomes multivalued. This effect relating to the discrete nature of the particle phase is well evidenced in the Eulerian–Lagrangian simulations used in this study. Owing to the coupling with the fluid, these small-scale clusters induce large gradients in the fluid vorticity field and cause a significant increase of enstrophy during the clustering stage. This is in contrast to single-phase 2-D vortices where enstrophy decreases monotonically. The accelerated decay of the particle-laden vortex is partly due to preferential concentration on the scale of the vortex core, and partly due to the modulation of enstrophy by the small-scale clusters. This work shows that the injection of inertial particles could provide an effective strategy for the control and suppression of resilient vortex tubes.
Declaration of interests
The authors report no conflict of interest.
Appendix A. Effect of varying grid resolution on the azimuthally averaged profiles
In this appendix, we investigate the effect of grid resolution on the Eulerian fluid and particle data obtained by our numerical method presented in § 2. For this grid dependence study, we focus on the case with 
$\textit {St}_{\varGamma }=0.4$ particles in a Lamb–Oseen vortex at Reynolds number 
$\textit {Re}_{\varGamma }=1700$.
The grid spacing is varied to yield three different resolutions. We denote as ‘fine’ the resolution obtained using a Cartesian grid of size 
$3072\times 3072$. This resolution is the one described in § 3, and used elsewhere in this study. Two additional simulations are conducted with ‘medium’ and ‘coarse’ resolutions corresponding to grid sizes 
$1536\times 1536$ and 
$768\times 768$, respectively. In these two simulations, the number of grid points across the initial vortex core scale is 51 and 25, respectively.
Figure 15(a) shows radial profiles of the fluid vorticity normalized by the initial centre vorticity at instants 
$t/\tau _f=0$, 4 and 60. As the grid resolution increases, the curves collapse onto the data with ‘fine’ resolution. This collapse holds true throughout the time integration.
Effect of varying the grid resolution from fine (solid), medium (dash-dotted), to coarse (dashed) in Eulerian–Lagrangian simulations of particle-laden Lamb–Oseen vortex with 
$\textit {St}_{\varGamma }=0.4$ and Reynolds number 
$\textit {Re}_{\varGamma }=1700$. (a) Radial profile of fluid vorticity at 
$t/\tau _f=0$, 4 and 60. (b) Azimuthally averaged volume fraction distribution at 
$t/\tau _f=0$, 15 and 30.
Figure 15(b) shows the azimuthally averaged volume fraction field at 
$t/\tau _f=0$, 15 and 90. The figure shows that changes in the volume fraction profiles are small and diminish with increasing resolution. The reader interested in a detailed study of the grid convergence behaviour of the method, and how it compares to other Eulerian–Eulerian methods, is referred to the study of Patel et al. (Reference Patel, Desjardins, Kong, Capecelatro and Fox2017).
Appendix B. Effect of particle initial velocity
In simulations of particle-laden flows, the choice of initial particle velocity has an impact on the evolution of the two-phase flow, particularly when the particles have finite inertia and two-way coupling is significant. Owing to the drag term in (2.3), the slip velocity 
$\boldsymbol {u}_s=\boldsymbol {u}_p-\boldsymbol {u}_{f|p}$ controls the magnitude of the force exerted by the particle phase onto the fluid phase. When the initial slip velocity 
$\boldsymbol {u}_s$ is zero, the initial drag force is also zero. This allows the fluid and particle phases to adapt to each other gradually. In contrast, if the initial slip velocity is 
$\boldsymbol {u}_s\neq 0$, the non-zero initial drag force may cause a delayed response of the particles and/or premature modification of the fluid phase.
To illustrate this aspect, we consider three different initializations of the particle-laden Lamb–Oseen vortex at 
$\textit {St}_{\varGamma }=0.4$ and 
$\textit {Re}_{\varGamma }=1700$. In the first initialization, the particles are seeded with zero slip velocity, i.e. they have velocities equal to the fluid velocity at the particle locations. In the second initialization, the particles are initialized with zero velocities. In this case, the non-dimensional slip velocity 
$\boldsymbol {u}_s/(\varGamma /2{\rm \pi} r_0)$ varies between 0 and 1 depending on the distance from the vortex core. In the third initialization, the particles are seeded with random velocities in the range 
$-\varGamma /2{\rm \pi} r_0$ to 
$\varGamma /2{\rm \pi} r_0$ corresponding to non-dimensional slip velocity 
$\boldsymbol {u}_s/(\varGamma /2{\rm \pi} r_0)$ between 
$-2$ and 2.
The evolution of the radial profiles of volume fraction are shown in figures 16(a–c). Compared to the case where the particles are initialized with zero slip velocity, the radial ejection and clustering of the particles are delayed when the initial slip velocity is non-zero. In these cases, it takes longer for the particles to reach their equilibrium velocity with the fluid, and to start being centrifuged out. Figure 16(d) shows that the vortex core grows slower when the initial slip velocity is non-zero. This behaviour is expected since the vortex core radius grows as a result of the outward forcing exerted by the particles. In the simulations where the initial slip velocity is non-zero, the particles take longer to start centrifuging out. Thus, the delayed formation of the annular particle cluster results in the delayed growth of the vortex core radius.
Comparison of the evolution of a particle-laden Lamb–Oseen vortex at 
$\textit {St}_{\varGamma }=0.4$ and 
$\textit {Re}_{\varGamma }=1700$ with different initial particle velocities: (
$\circ$) matching the fluid velocity at the particle locations; (
$\square$) zero velocity; and (
$\diamond$) random velocity. (a–c) Normalized radial profile of volume fraction at successive instants. (d) Vortex core radius. The lines in (d) correspond to power-law fits of the type 
$A(t/\tau _f)^m$.
Open access
