1. Introduction
Particle-laden turbulent flows are prevalent in both natural phenomena and industrial applications, such as sandstorms (Di Renzo & Urzay Reference Di Renzo and Urzay2018; Zhang & Zhou Reference Zhang and Zhou2020), combustors (Apte et al. Reference Apte, Mahesh, Moin and Oefelein2003) and fluidized beds (Rokkam, Fox & Muhle Reference Rokkam, Fox and Muhle2010; Kolehmainen et al. Reference Kolehmainen, Ozel, Boyce and Sundaresan2016). These flows are inherently complicated due to the existence of turbulent fluctuation and also the boundary layer in shaping the intricate movement of particles. Despite its fundamental significance in multiphase flows, understanding the motion and distribution of particles within these systems is also important for efficiently manipulating the momentum or heat transports in industrial processes. There are also a couple of reviews regarding multiphase turbulent flow (Balachandar & Eaton Reference Balachandar and Eaton2010; Voth & Soldati Reference Voth and Soldati2017; Elghobashi Reference Elghobashi2019; Mathai, Lohse & Sun Reference Mathai, Lohse and Sun2020; Brandt & Coletti Reference Brandt and Coletti2022; Ni Reference Ni2024).
Previous studies investigating the interaction between turbulent flows and particles have often focused on homogeneous isotropic turbulence (HIT) as an idealized system. In HIT, it has been observed that inertial particles exhibit a non-uniform spatial distribution characterized by clusters and voids, a phenomenon known as preferential concentration (Squires & Eaton Reference Squires and Eaton1991; Saw et al. Reference Saw, Shaw, Ayyalasomayajula, Chuang and Gylfason2008; Monchaux, Bourgoin & Cartellier Reference Monchaux, Bourgoin and Cartellier2012). Through calculating the divergence of the particle velocity, Maxey (Reference Maxey1987) found that particles tend to accumulate in regions of high strain rate and low vorticity in the small Stokes number limit (
$St \ll 1$). At first, this phenomenon was explained by the centrifugal effect caused by vortices in the turbulent flow (Abrahamson Reference Abrahamson1975; Reade & Collins Reference Reade and Collins2000). As 
$St$ increases, the coupling between the particle dynamics and the local fluid velocity field decreases, while particle path-history interactions with turbulence become more significant (Bragg & Collins Reference Bragg and Collins2014; Bragg, Ireland & Collins Reference Bragg, Ireland and Collins2015b). In high Reynolds number flows, particles with large 
$St$ can cluster due to acceleration by high-velocity gradients, known as the sling effect (Falkovich, Fouxon & Stepanov Reference Falkovich, Fouxon and Stepanov2002; Gustavsson & Mehlig Reference Gustavsson and Mehlig2016). Goto & Vassilicos (Reference Goto and Vassilicos2008) found that point particles accumulate in regions with positive flow acceleration divergence, termed the sweep-stick mechanism. Bragg, Ireland & Collins (Reference Bragg, Ireland and Collins2015a) introduced a Stokes number (
$St_r$) based on the eddy turnover time scale, suggesting drift mechanisms for particle clustering, and noted that the sweep-stick mechanism is valid only for 
$St_r \ll 1$ in the inertial range.
Turbulent flows encountered in industrial settings are often bounded by solid walls. In such wall-bounded flows, the presence of no-slip and no-penetration boundary conditions leads to a sharp decrease in turbulence intensity near the wall. Consequently, particles disperse more rapidly in regions of higher turbulence intensity (Caporaloni et al. Reference Caporaloni, Tampieri, Trombetti and Vittori1975; Reeks Reference Reeks1983), resulting in a tendency for particles to accumulate in the viscous sublayer at higher concentrations compared with other regions (Marchioli & Soldati Reference Marchioli and Soldati2002; Bernardini Reference Bernardini2014). This phenomenon is known as turbophoresis, where several studies have indeed demonstrated the high particle concentration in the near-wall regions compared with the bulk (Marchioli & Soldati Reference Marchioli and Soldati2002; Sardina et al. Reference Sardina, Schlatter, Brandt, Picano and Casciola2012). It has been observed that particle accumulation at the wall becomes pronounced when the particle relaxation time matches the local turbulence time scale. Additionally, numerical simulations have shown that particles moving away from the wall are associated with the ejection events (Vinkovic et al. Reference Vinkovic, Doppler, Lelouvetel and Buffat2011). Guha (Reference Guha1997, Reference Guha2008) discussed the transport and deposition of particles and theoretical modelling and this theoretical model can quantify the physical mechanisms for particle transport, such as Brownian motion, turbulent diffusion, turbophoresis, etc. Instead of a smooth wall, Zhao & Wu (Reference Zhao and Wu2006) have examined the particle disposition over rough walls in a ventilation duct. Zahtila et al. (Reference Zahtila, Chan, Ooi and Philip2023) have appraised various Eulerian-based models of turbophoresis and turbulent diffusive coefficients over a wide range of particle Stokes numbers in turbulent pipe flow. Johnson, Bassenne & Moin (Reference Johnson, Bassenne and Moin2020) have recently examined the effects of biased sampling and turbophoresis on particle concentration profiles in turbulent channel flows, applying their model to wall-modelled large-eddy simulations. Recently, Zhang, Cui & Zheng (Reference Zhang, Cui and Zheng2023) used a theoretical model based on Johnson et al. (Reference Johnson, Bassenne and Moin2020) to explain the physical mechanisms of the wall-normal concentration profiles for different charged particles in channel flow.
In certain natural and industrial environments, additional forces come into play, further complicating the dynamical system. In systems subject to rotation, such as Taylor vortex reactors, centrifugal forces significantly alter the spatial distribution of particles by driving them away from the centre of rotation. To investigate the combined effects of shear turbulence and rotation, the Taylor–Couette (TC) flow (Couette Reference Couette1890; Taylor Reference Taylor1923b; Grossmann, Lohse & Sun Reference Grossmann, Lohse and Sun2016) serves as a paradigmatic system. The TC flow consists of fluid confined in the gap between two rotating cylinders. The dimensionless parameters regarding the geometry can be described by the radius ratio 
$\eta = r_i/r_o$ and the aspect ratio 
$\varGamma = L/d$, where the 
$r_i$ (
$r_o$) is the inner (outer) cylinder radius, 
$L$ and 
$d$ are the axial length and the gap between two cylinders, respectively. As the rotation of the inner and outer cylinders drive the flow, the two governing dimensionless parameters are the two Reynolds numbers (Grossmann et al. Reference Grossmann, Lohse and Sun2016), namely
\begin{equation} R e_{i,o}=\frac{{r_{i,o}} \omega_{i,o} d}{\nu}, \end{equation}
where 
$Re_i$ and 
$Re_o$ are the Reynolds numbers of the inner and outer cylinders, respectively, 
$\omega _i$ and 
$\omega _o$ are the angular velocities of inner and outer cylinders, respectively, 
$\nu$ is the kinematic viscosity of fluid. When the outer cylinder is fixed, the relevant control parameter is solely given by the inner cylinder Reynolds number 
$Re_i = r_i \omega _i d / \nu$. Alternatively, one can characterize the driving of the flow by Taylor number, which is
\begin{equation} T a=\frac{(1+\eta)^4}{64 \eta^2} \frac{\left(r_{o}-r_{i}\right)^2\left(r_{i}+r_{o}\right)^2\left(\omega_{i}-\omega_{o}\right)^2}{\nu^2}, \end{equation}
and regarding the rotation ratio, one can use the inverse Rossby number 
$Ro^{-1}$.
The pioneering work of flow structure in TC flow was done by Taylor (Reference Taylor1923b), who revealed the existence of the Taylor vortex. Since then, four distinct flow patterns have been found at different control parameters, namely circular Couette flow (CCF), Taylor vortex flow (TVF), wavy vortex flow (WVF) and turbulent TC flow. Due to the rich variety of flow characteristics in TC flow, this system has gained significant attention. The foci of the studies are on the flow instabilities (Taylor Reference Taylor1923b; Rüdiger et al. Reference Rüdiger, Gellert, Hollerbach, Schultz and Stefani2018), ultimate regimes (Froitzheim et al. Reference Froitzheim, Ezeta, Huisman, Merbold, Sun, Lohse and Egbers2019; Hamede, Merbold & Egbers Reference Hamede, Merbold and Egbers2023), properties of inner and outer boundary layers (Brauckmann & Eckhardt Reference Brauckmann and Eckhardt2017) and pattern formation (Koschmieder Reference Koschmieder1993). There are also studies focusing on multiphase TC flow, where the important issue is the drag modulation led by dispersed phases, such as bubbles (van den Berg et al. Reference van den Berg, Luther, Lathrop and Lohse2005; Spandan et al. Reference Spandan, Ostilla-Mónico, Verzicco and Lohse2016b; Spandan, Verzicco & Lohse Reference Spandan, Verzicco and Lohse2018), droplets (Spandan, Lohse & Verzicco Reference Spandan, Lohse and Verzicco2016a; Su et al. Reference Su, Zhang, Wang, Yi, Xu, Fan, Wang and Sun2025), polymers (Zhang et al. Reference Zhang, Fan, Su, Xi and Sun2025) and finite-sized particles (Wang et al. Reference Wang, Yi, Jiang and Sun2022). Besides, Bakhuis et al. (Reference Bakhuis, Verschoof, Mathai, Huisman, Lohse and Sun2018) have examined how the finite size spherical particles influence the drag in turbulent TC flow. Instead of focusing on spherical objects, cases with rigid fibres and ellipsoidal particles in TC flow have also been studied (Bakhuis et al. Reference Bakhuis, Mathai, Verschoof, Ezeta, Lohse, Huisman and Sun2019; Assen et al. Reference Assen, Ng, Will, Stevens, Lohse and Verzicco2022). In the meantime, there are also several studies examining the mixing problems in chemical processing applications (Schrimpf et al. Reference Schrimpf, Esteban, Warmeling, Färber, Behr and Vorholt2021) and the biomedical field (Curran & Black Reference Curran and Black2005).
Taylor–Couette flow is an ideal system for investigating dispersed multiphase flows with the interplay of shear turbulence and rotation. Studies on dispersed multiphase TC flow date back to the 1920s (Taylor Reference Taylor1923a). Subsequently, many studies found the drift of neutrally buoyant particles from high to low shear regions in TC flows (Tetlow et al. Reference Tetlow, Graham, Ingber, Subia, Mondy and Altobelli1998; Fang et al. Reference Fang, Mammoli, Brady, Ingber, Mondy and Graham2002; Qiao, Deng & Wang Reference Qiao, Deng and Wang2015; Li, Ku & Lin Reference Li, Ku and Lin2020; Kang & Mirbod Reference Kang and Mirbod2021). A comprehensive experiment conducted by Majji & Morris (Reference Majji and Morris2018) observed the inertial migration of neutral particles in CCF, TVF and WVF with Reynolds numbers ranging from 83 to 151.4. They reported the effect of the neutral particle inertial migration on the transition of TC flows (Baroudi, Majji & Morris Reference Baroudi, Majji and Morris2020). Li et al. (Reference Li, Ku and Lin2020) adopted Eulerian–Lagrangian simulation to examine the point-particle equilibrium position in CCF with Reynolds numbers ranging from 60 to 90. Dash, Anantharaman & Poelma (Reference Dash, Anantharaman and Poelma2020) considered cases with larger parameter range with the Reynolds number extended to 
$O(10^3)$. They found the distinct role played by the particles where, instead of simply destabilizing the flow, neutral particles may also decrease the growth of flow instabilities.
However, existing studies predominantly focus on particle densities close to that of the carrier phase, also often neglecting turbulent flow regimes in multiphase TC flow. The question of how particle inertia plays a role in shaping their dispersion behaviour and distribution in multiphase TC flow remains poorly addressed. Therefore, the primary objective of this study is to explore the influence of particle inertia on the spatial distribution of particles in turbulent TC flow.
In the present study, we analyse the effect of particle inertia on particle preferential concentration. We further examine how varying densities of particles influence the radial distribution of particles. Our objective is to uncover the mechanisms governing particle radial concentration profiles at various 
$St$. The remainder of the paper is organized as follows. The governing equations and problem set-up are introduced in § 2. Then, the results are shown in § 3. Concluding remarks are given in § 4.
2. Governing equations and problem set-up
In the present investigation, we carry out direct numerical simulations (DNS) of three-dimensional (3-D) particle-laden turbulent TC flow, as shown in figure 1.
Figure 1. Sketch of the set-up of particle-laden TC flow.
2.1. Carrier phase
For the carrier phase, the governing equations are non-dimensionalized using the inner cylinder rotation speed, 
$\omega _i r_i$, and the gap between the cylinders, 
$d$, which are given by
\begin{gather} \frac{\partial \pmb{u}}{\partial t} + (\pmb{u} \boldsymbol{\cdot}\boldsymbol{\nabla}) \pmb{u}={-}\boldsymbol{\nabla} p+\frac{1}{Re_i} \nabla^{2} \pmb{u}, \end{gather}
\begin{gather}\boldsymbol{\nabla}\boldsymbol{\cdot} \pmb{u}=0. \end{gather}
In (2.1), u and p represent the velocity and pressure of the fluid, respectively, and 
$Re_i = d \omega _i r_i/\nu$ is the Reynolds number defined by the inner cylinder rotation.
Direct numerical simulation of the carrier phase is performed using a second-order accurate finite-difference method in cylindrical coordinates 
$( \pmb {e}_r, \pmb {e}_{\theta }, \pmb {e}_z)$ (Verzicco & Orlandi Reference Verzicco and Orlandi1996; Ostilla et al. Reference Ostilla, Stevens, Grossmann, Verzicco and Lohse2013) with uniform grid spacing in the azimuthal and axial directions and a non-uniform grid spacing using a clipped Chebyshev-type clustering method in the radial direction. The computational domain is 
$r_i \leqslant r \leqslant r_o$, 
$0 \leqslant \theta \leqslant 2{\rm \pi}$ and 
$0 \leqslant z \leqslant L_z$, where 
$r_i$ and 
$r_o$ are the inner and outer cylinder radii, respectively, and 
$L_z$ is the axial length of the domain. The no-slip boundary condition is applied at the inner and outer cylinder walls, and a periodic boundary condition is applied in the azimuthal and axial directions.
The numerical approach has been adopted in many of our previous studies to simulate the problem of single phase and multiphase turbulent flows under a Cartesian geometry (Zhao et al. Reference Zhao, Zhang, Wang, Wu, Chong and Zhou2022Reference Zhao, Wang, Wu, Chong and Zhoua,Reference Zhao, Zhang, Wang, Wu, Chong and Zhoub; Guo et al. Reference Guo, Wu, Wang, Zhou and Chong2023; Meng et al. Reference Meng, Zhao, Wu, Wang, Zhou and Chong2024; Zhang & Zhou Reference Zhang and Zhou2024). In order to further validate our numerical code within a cylindrical geometry, we performed simulations of TC flow with the same grid resolutions and control parameters as in Ostilla et al. (Reference Ostilla, Stevens, Grossmann, Verzicco and Lohse2013), and compared our results for 
$Nu_\omega$ in table 1. Here, 
$Nu_\omega$ is defined as
\begin{equation} N u_\omega=\frac{J}{J_{l a m}}, \quad \text{with } J=r^3\left( \left\langle u_r \omega\right\rangle_{A, t}-v \partial_r\langle\omega\rangle_{A, t}\right), \end{equation}
where 
$\langle {\cdot } \rangle _{A, t}$ denotes the azimuthal, axial and time average, and 
$J_{lam}$ is the angular velocity flux for laminar flow. The results are in good agreement with the results of Ostilla et al. (Reference Ostilla, Stevens, Grossmann, Verzicco and Lohse2013). In this study, we used the control parameters in the second row of table 1, and the flow state of the TC system under this parameter state is a classical regime with boundary layers of the Prandtl–Blasius type and a turbulent bulk with Taylor rolls (Grossmann et al. Reference Grossmann, Lohse and Sun2016). Figure 2 presents the profiles of the mean azimuthal velocity and root-mean-square (r.m.s.) of the fluid velocity fluctuation. The mean azimuthal velocity profile is characterized by the presence of Taylor rolls, where the fluid rotates more rapidly near the inner cylinder and more slowly near the outer cylinder. The azimuthal and radial components of velocity fluctuations rise sharply along the wall-normal direction, exhibiting distinct peaks near the walls, and decrease towards the channel centre. This behaviour is similar to the streamwise and wall-normal r.m.s. fluctuations observed in channel flow. The axial r.m.s. increases from the wall, reaching its maximum value at the centre of the gap.
Table 1. Validation of TC simulations at 
$Re_i=160,2500$. Here, 
$Re_i$ is the Reynolds number of the inner cylinder, 
$N_{\theta } \times N_r \times N_z$ the grid resolution, 
$\varGamma$ the aspect ratio, 
$\eta$ the radius ratio, 
$Nu_{\omega }$ the normalized angular velocity flux and 
$\epsilon$ is the relative error of 
$Nu_\omega$ comparing with Ostilla et al. (Reference Ostilla, Stevens, Grossmann, Verzicco and Lohse2013).
Figure 2. Radial profiles of the mean azimuthal velocity (
$a$) and the r.m.s. of the fluid velocity fluctuation (
$b$) for the Reynolds number 
$Re_i=2500$.
2.2. Lagrangian particle tracking
In this study, we consider the inertial particles using the Lagrangian point-particle approach (Gatignol Reference Gatignol1983; Maxey & Riley Reference Maxey and Riley1983; Maxey Reference Maxey1987; Tsai Reference Tsai2022). Particles are released into the flow when the turbulent flow is fully developed. Low volume fractions of particles have been considered, and thus the feedback of particles to the flow field and the particle–particle collisions have been neglected (Elghobashi Reference Elghobashi1991, Reference Elghobashi1994; Elghobashi, Balachandar & Prosperetti Reference Elghobashi, Balachandar and Prosperetti2006; Balachandar & Eaton Reference Balachandar and Eaton2010). A small diameter 
$d_p$ and large density ratio 
$\rho ^* = \rho _p/\rho _f$ are considered in this study. The Lagrangian particle tracking considers Stokes drag in the governing equation of the particle motion. The dimensionless equations of particle translational motion in cylindrical coordinates read as
\begin{gather} \frac{\mathrm{d} v_r}{\mathrm{d} t}=\frac{C_D}{St}\left(u_r-v_r\right)+\frac{v_\theta^2}{r}, \end{gather}
\begin{gather}\frac{\mathrm{d} v_\theta}{\mathrm{d} t}=\frac{C_D}{St}\left(u_\theta-v_\theta\right)-\frac{v_\theta v_{\mathrm{r}}}{r}, \end{gather}
\begin{gather}\frac{\mathrm{d} v_z}{\mathrm{d} t}=\frac{C_D}{St}\left(u_z-v_z\right), \end{gather}
where 
$C_D = 1+ 0.15Re_p^{0.687}$ is the drag coefficient, 
$Re_p = d_p^* |\pmb {u} - \pmb {v}|Re_i$ is the particle Reynolds number and 
$\pmb {v}$ and 
$\pmb {u}$ are the particle and fluid dimensionless velocity, respectively. The particle diameter 
$d_p^*$ is normalized by the gap width 
$d$; note that the particles are treated as point-like particles, where the particle diameter is a virtual size for 
$St$ and 
$Re_p$ calculations only. Also, 
$v_r$, 
$v_\theta$ and 
$v_z$ are the particle velocities in the radial, azimuthal and axial directions, respectively. The particle Stokes number is defined as 
$St = \rho ^* d_p^{\ast 2} Re_i/18$, where 
$\rho ^* = \rho _p/\rho _f$ is the density ratio between particle and fluid. The particles experience periodic boundary conditions in the azimuthal and axial directions. For radial boundary conditions, complete elastic collisions occur at the smooth wall (Marchioli et al. Reference Marchioli, Soldati, Kuerten, Arcen, Taniere, Goldensoph, Squires, Cargnelutti and Portela2008). To compute the value of physical quantities at the particle location, we use a tri-linear interpolation method for the interpolation and compute the values at the particle centre position.
In the wall-bounded flow, there are three ways to determine the Stokes number, as listed in (2.7)–(2.8)
\begin{gather} St = \frac{\tau_p }{H/U_{bulk}}, \end{gather}
\begin{gather}St_{\eta} = \frac{\tau_p}{\tau_{\eta}}, \end{gather}
\begin{gather}St^+ = \frac{\tau_p }{\nu / u_{\tau}^2}. \end{gather}
In (2.7), the Stokes number is based on the bulk velocity 
$U_{bulk}$ and the distance of the walls 
$H$. For the TC system, the bulk Stokes number can be defined based on the inner cylinder rotation speed 
$\omega _i r_i$ and the cylinder gap 
$d$, which is equivalent to the 
$St$ in (2.4)–(2.6). The Kolmogorov Stokes number based on the Kolmogorov time scale is defined in (2.8), where 
$\tau _{\eta } = (\nu /\epsilon )^{1/2}$, where 
$\epsilon$ is the energy dissipation rate. The viscous Stokes number (based on friction velocity) is defined in (2.9), where 
$u_{\tau } = \sqrt {\tau _{w}/\rho }$ represents the friction velocity. In TC flow, the mean wall shear stress, 
$\tau _{w}$, is defined based on torque (Huisman et al. Reference Huisman, Scharnowski, Cierpka, Kähler, Lohse and Sun2013). Since the angular velocity flux 
$J$ and torque 
$\mathcal {T}$ are the same for both cylinders, 
$\tau _{w}$, 
$u_{\tau }$ and 
$\delta _{\nu }$ at the inner and outer cylinders follow these ratios
\begin{gather} \tau_{w,i}/\tau_{w,o} = 1/{\eta}^2, \end{gather}
\begin{gather}u_{\tau,i}/u_{\tau,o} = 1/\eta, \end{gather}
\begin{gather}\delta_{\nu,i}/\delta_{\nu,o} = \eta. \end{gather}
Figure 3 shows the azimuthal velocity profiles normalized by the friction velocity. In this study, the ratio of particle diameter to the gap of cylinder is 
$d_p^* = 0.005$, and the density ratio ranges from 10 to 288 for particles with different inertia. The corresponding Stokes numbers are listed in table 2. The normalized particle diameters are based on the inner and outer viscous length scales as 
$d_p/\delta _{\nu,i} = 0.8$ and 
$d_p/\delta _{\nu,o} = 0.57$, and the normalized particle diameter is based on the Kolmogorov length scale with 
$d_p/\eta = 0.33$. The volume fraction of particles is fixed at 
$\phi _v = 5.5 \times 10^{-5}$, and at this small volume fraction and Stokes number, one-way coupling is sufficient to describe the particle motion.
Figure 3. Azimuthal velocity profiles near the inner (
$(r-r_i)/d \in (0,1/2)$) and outer (
$(r_o-r)/d \in (1/2,1)$) cylinders are presented for the Reynolds number 
$Re_i=2500$. The distance from the wall, 
$y^+$, is normalized by the viscous length scale 
$\delta _{\nu }$, and the azimuthal velocity, 
$u^+$, is scaled by the friction velocity 
$u_{\tau }$. For the inner boundary layer, 
$y^+$ is defined as 
$(r-r_i)/\delta _{\nu,i}$, and 
$u^+$ is defined as 
$(\omega _i r_i - u_{\theta }(r)) / u_{\tau,i}$. For the outer boundary layer, 
$y^+$ is defined as 
$(r_o-r)/\delta _{\nu,o}$, and 
$u^+$ is defined as 
$u_{\theta }(r)/u_{\tau,o}$.
Table 2. Details of the inertial particle in the numerical simulations. The value of 
$d_p^* = d_p / d$ is the particle diameter normalized by the gap width, 
$\phi _v$ is the volume fraction, 
$\rho ^* = \rho _p/\rho _f$ is the density ratio between the particle and fluid, 
$St=\rho ^\ast d_p^{\ast 2}Re_i/18$ is the particle bulk Stokes number and 
$St^+_{i,o}=\tau _p/\nu u_{\tau,i,o}^2$ are the viscous Stokes numbers based on the inner and outer boundary layers. The normalized particle diameters are based on inner and outer viscous length scales as 
$d_p/\delta _{\nu,i} = 0.8$ and 
$d_p/\delta _{\nu,o} = 0.57$. The Kolmogorov Stokes number is 
$St_{\eta }= \tau _p/\tau _{\eta }$, and the normalized particle diameter is based on the Kolmogorov length scale with 
$d_p/\eta = 0.33$.
3. Results and discussion
3.1. Particles preferential concentration
To study how the particle inertia influences its distribution in TC flow, we first present the 3-D visualization of the instantaneous particle distributions with varying Stokes numbers, as illustrated in figure 4. Different patterns of particle distribution emerge depending on 
$St$. When the particle inertia is minimal (
$St=0.034$), the particles are dispersed throughout the entire domain. With an increase in particle inertia (i.e. with higher 
$St$), particles tend to cluster, forming distinct parallel stripe-like structures in the axial direction. For example, four clearly defined strips become evident at 
$St=1$, as seen in figure 4(
$c$). We also draw the iso-surface of the radial velocity in the 3-D visualization, which shows the footprint of Taylor rolls – a hallmark feature of TC flow. One sees that these stripe-like particle clusters coincide with the inward velocity region of Taylor rolls. From the top view of the particle distribution shown in figures 4(
$d$)–4(
$\,f$), it is observed that the striped particle clusters are located near the outer walls as the particle inertia becomes more pronounced with increasing 
$St$. To more clearly illustrate the spatial distribution of the particles, we plot their distribution in the 
$(r, t)$ plane, as shown in figure 5. Particles with small inertia are uniformly distributed in space. As particle inertia increases, they tend to accumulate at the edges of the Taylor rolls. In the case of the largest Stokes number, most particles concentrate between the Taylor roll pairs near the outer cylinder.
Figure 4. Three-dimensional visualization (a–c) and top view (d–f) instantaneous snapshots of the particle distribution for the (a,c) 
$St=0.034$, (b,e) 
$St=0.34$, (c,f) 
$St=1$; The iso-surface of radial velocity has also been drawn in the 3-D visualization (a–c) where the reddish surface represents the outward velocity and the bluish one represents the inward velocity, particularly at the dimensionless radial velocities 
$-0.1$ and 
$0.1$. The size of the particles has been magnified approximately 20 times for better visualization.
Figure 5. Instantaneous snapshots of particle distribution in the 
$(r,z)$ plane for (
$a$) 
$St=0.034$, (
$b$) 
$St=0.34$ and (
$c$) 
$St=1$. The contour plots are the angular velocity of fluid at 
$\theta = {\rm \pi}$, and the particles are selected near the 
$\theta = {\rm \pi}$ region.
In numerous industrial applications, the accumulation of particles near the wall region often leads to collisions or adhesions between particles, making this phenomenon a key aspect to examine. To characterize the particle distribution in the vicinity of the outer wall, we utilize a two-dimensional (2-D) Voronoï diagram analysis to quantify the preferential concentration of particles. The Voronoï diagram is a spatial tessellation in which each Voronoï cell is defined by the particle location based on the distance to adjacent particles. Consequently, in regions where particles cluster, the area of the Voronoï cells are smaller compared with those in adjacent regions. For the analysis of the particle distribution in the near outer-wall region, we consider particles in the region of 
$r_0-1.5d_p^* \leqslant r \leqslant r_0-0.5d_p^*$. As an example, figure 6 illustrates the 2-D Voronoï diagrams depicting the instantaneous distribution of particles near the outer wall for 
$St = 0.034$.
Figure 6. (
$a$) Instantaneous snapshot of particle distribution at near the outer-wall region for the 
$St = 0.034$, (
$b$) the 2-D Voronoï diagrams corresponding to (
$a$).
As the area of the Voronoï cells can represent particle clustering, we examine the probability density functions (p.d.f.s) of Voronoï cell areas to quantify the preferential concentration of particles. The p.d.f.s of particle Voronoï area in the near outer-wall region are presented in figure 7(
$a$), with the solid lines representing the p.d.f.s of simulation data. The dashed line corresponds to the p.d.f. of Voronoï cell areas (normalized by the mean area) for randomly distributed particles, as described by a 
$\varGamma$-distribution (Ferenc & Néda Reference Ferenc and Néda2007).
Figure 7. Statistics of 2-D Voronoï diagrams for particles near the outer wall at various Stokes numbers: (
$a$) the normalized Voronoï area p.d.f.s, where the grey dashed curve shows the 
$\varGamma$-distribution for a 2-D random distribution; (
$b$) the normalized standard deviation 
$\sigma / \sigma _{\varGamma }$ of the Voronoï area distribution; (
$c$) the relative p.d.f. defined by the ratio of the p.d.f.s to the 2-D 
$\varGamma$-distribution counterpart. Here, the reddish (greenish) area denotes the regime having cluster (void) of particles; (
$d$) the p.d.f.s of the Voronoï cells’ aspect ratio, defined by the ratio of azimuthal length to axial length, where the grey dashed curve shows the case for 2-D random distribution.
The 
$\varGamma$-distribution of 2-D Voronoï area is approximated by
\begin{equation} f(A^+)=\frac{343}{15} \sqrt{\frac{7}{2 {\rm \pi}}} (A^+)^{5 / 2} \exp \left(-\frac{7}{2} A^+\right). \end{equation}
Here, 
$A^+ = A/\bar {A}$ represents the Voronoï cell area normalized by the mean value. At 
$St=0.034$, the p.d.f. of the Voronoï cell area closely aligns with the 
$\varGamma$-distribution, with only slight deviations in the regime where 
$A/\bar {A}<0.6$. However, as the particle 
$St$ increases, the deviation from the 
$\varGamma$-distribution becomes more pronounced. Moreover, the p.d.f. values are not only larger than the 
$\varGamma$-distribution in the small area regime but also in the large area regime. Our results suggest that the preferential concentration of particles becomes increasingly significant with an increase in 
$St$. Next, we investigate the variation of the standard deviation of the Voronoï cell area normalized by the standard deviation of the 
$\varGamma$-distribution (Ferenc & Néda Reference Ferenc and Néda2007), with respect to the particle Stokes number (
$St$), as illustrated in figure 7(
$b$). By comparing 
$\sigma /\sigma _{\varGamma }$ across different 
$St$ values, we observe a clear trend of increasing particle preferential concentration with higher Stokes numbers.
To closely inspect the particle preferential distribution, we examine the ratio of respective p.d.f.s to the 
$\varGamma$-distribution p.d.f. counterpart, namely the relative p.d.f. (Monchaux, Bourgoin & Cartellier Reference Monchaux, Bourgoin and Cartellier2010) in figure 7(
$c$). For each relative p.d.f. there are two intersection points with the value 
$1$. According to Monchaux et al. (Reference Monchaux, Bourgoin and Cartellier2010), these intersections present the formation of clusters (red) and voids (green). With the increase in 
$St$, particle clusters become more pronounced, which is consistent with the results shown in figures 7(
$a$) and 7(
$b$).
In addition, to further quantify the characteristics of the particle distribution near the outer wall, we examine the p.d.f. of the aspect ratio of the Voronoï cell: 
$l_{\theta,A}/l_{z,A}$. Here, 
$l_{\theta,A}$ is the maximum length in the azimuthal direction of a Voronoï cell and 
$l_{z,A}$ is the maximum length in the axial direction. In figure 7(
$d$), it is evident that the p.d.f. curves of 
$l_{\theta,A}/l_{z,A}$ on the left are significantly higher than that of random distribution case. This indicates that particles prefer to cluster in the axial direction, which can also be seen in figures 8(
$a$)–8(
$c$). In addition, as 
$St$ increases, the left tail of the p.d.f. gets higher and the peak shifts to the left. This suggests that the greater the inertia of the particles, the more pronounced the aggregation of particles near the outer wall and the formation of an azimuthal stripe structure.
Figure 8. (
$a$–
$c$) Snapshots depicting instantaneous distributions of particles at different Stokes numbers in the near outer-wall (
$\theta \unicode{x2013}z$) plane, with the colour representing the fluctuation of the azimuthal velocity 
$u^{\prime }_{\theta }$ of the fluid. (
$d$) The p.d.f.s of the fluid azimuthal velocity fluctuation. (
$e$) Radial profiles of the particle numbers with 
$u^{\prime }_{\theta }>0$ compared with those with 
$u^{\prime }_{\theta }<0$. The grey dashed curve represents the result for a random distribution of particles.
Subsequently, we analyse the radial profile of the ratio of particle numbers with 
$u^{\prime }_{\theta }>0$ to those with 
$u^{\prime }_{\theta }<0$, as depicted in figure 8(
$e$). We can see from figure 8(
$e$) that, for particles with 
$St>0.034$, the number of particles with 
$u^{\prime }_{\theta }>0$ is lower than that with 
$u^{\prime }_{\theta }<0$, implying that particles prefer to accumulate in low-speed streaks, and note that this conclusion is only valid if the Eulerian p.d.f. of 
$u_{\theta }^{\prime }$ is symmetric, as shown in figure 8(
$d$). For particles with 
$St=0.034$, their trajectory is very close to that of the fluid particle (i.e. the tracer particle) due to their small inertia, and their velocity characteristics are close to those of randomly distributed particles.
3.2. Profile of particle radial concentration
In wall-bounded particle-laden turbulent flow, the wall-normal particle distribution is also worth studying. The presence of the wall boundary affects the movement and distribution of particles along the wall-normal direction. To examine the underlying mechanism governing the radial distribution of particle concentration in particle-laden TC flow and the dominant forces, we employ a theoretical model inspired by Johnson et al. (Reference Johnson, Bassenne and Moin2020), where their work considered the situation of plane Couette flow and disentangled the dominant effects into biased sampling and turbophoresis.
Following the concept of Johnson et al. (Reference Johnson, Bassenne and Moin2020), we derived the corresponding disentangled relation for the particle concentration in particle-laden TC flow with an additional centrifugal term
\begin{equation} C^*\left(r\right)=\alpha \exp \left(\underbrace{\frac{1}{St} \int^{r} \frac{\left\langle u_{r} |\eta\right\rangle}{\left\langle v_r^2 |\eta\right\rangle} \mathrm{d} \eta}_{\textit{biased sampling }}-\underbrace{\int^{r} \frac{\mathrm{d} \ln \left\langle v_{r}^{2} | \eta\right\rangle}{\mathrm{d} \eta} \mathrm{d} \eta}_{\textit{turbophoresis }}+\underbrace{\int^{r} \frac{\left\langle \dfrac{v_{\theta}^2}{r} | \eta\right\rangle}{\left\langle v_{r}^{2} |\eta\right\rangle} \mathrm{d} \eta}_{\textit{centrifugal effect }}\right), \end{equation}
which is derived from the radial direction single-particle position–velocity p.d.f. Detail of the derivation can be found in Appendix A. The numerical expression for the dimensionless concentration 
$C^*(r)$ is
\begin{equation} C^*\left(r\right)= \frac{n_rV}{n_0V_{slice}}. \end{equation}
Here, 
$n_r$ and 
$V_{slice}$ represent the particle numbers in each slice and slice volume, 
$n_0$ and 
$V$ are the total particle number and the volume of the entire domain, 
$u_r$ is the radial component of the fluid velocity sensed by the particles. 
$v_r$ and 
$v_{\theta }$ are the radial and angular components of the particle velocity, respectively, 
$\alpha$ is an integration constant and 
$\langle {\cdot } \rangle$ is the ensemble average done by averaging over all particles in the slice at the corresponding radial position. There are three terms known as phoresis integrals on the right-hand side of (3.2). The first term is the effect of biased sampling on the particle concentration profile
\begin{equation} I_b= \frac{1}{St}\int^{r} \frac{\left\langle u_{r} |\eta\right\rangle}{\left\langle v_r^2 |\eta\right\rangle} \mathrm{d} \eta. \end{equation}The second term represents the effect of turbophoresis on the particle concentration profile
\begin{equation} I_t={-} \int^{r} \frac{\mathrm{d} \ln \left\langle v_{r}^{2} | \eta\right\rangle}{\mathrm{d} \eta} \mathrm{d} \eta. \end{equation}The last term is given by
\begin{equation} I_c= \int^{r} \frac{\left\langle \left.\!\dfrac{v_{\theta}^2}{r} \right| \eta\right\rangle}{\left\langle v_{r}^{2} |\eta\right\rangle} \mathrm{d} \eta, \end{equation}which quantifies the centrifugal effect on the particle radial distribution.
Figure 9 shows the particle concentration profiles with various Stokes numbers (
$St$) at 
$Re_i = 2500$. The figure displays simulation results (solid curves) and the theoretical points (open symbols) from (3.2). It shows that the equation can nicely describe the variation of the particle concentration. The concentration of particles decreases significantly near the inner wall and in the bulk with increasing 
$St$, whereas there is a substantial increase near the outer wall.
Figure 9. Radial profiles of particle concentration at various Stokes numbers. The curves are the results from DNS and the symbols represent the results based on (3.2).
The benefit of using (3.2) is to disentangle the contributions from biased sampling, turbophoresis and the centrifugal effect. In such a way, one can examine how these three effects come into play in different regions of turbulent TC flow (near inner wall, outer wall and in the bulk). Thus, we examine the radial profiles of the sampling bias 
$I_b$ (3.4), turbophoresis 
$I_T$ (3.5) and centrifugal effect 
$I_c$ (3.6) integrals, as shown in figures 10, 11 and 12 respectively.
Figure 10. (
$a$) Radial profiles of negative of biased sampling integral (3.4); (
$b$) radial profiles of fluid radial velocity at particle location.
Figure 11. (
$a$) Radial profiles of turbophoresis integral (3.5); (
$b$) radial profiles of r.m.s. radial velocity of particles.
Figure 12. (
$a$) Radial profiles of centrifugal integral (3.6); (
$b$) radial profiles of azimuthal velocity of particles. As the reference, the black curve shows the radial profile of the fluid azimuthal velocity.
The first phoresis integral accounts for the average drag force on the ensemble of particles at a given radial location. As depicted in (3.4), it is proportional to 
$\langle u_{r} |\eta \rangle$, which represents the fluid radial velocity sampled at the particle position. For small 
$St$, a relatively random distribution of particles can be observed, as shown in figure 4(
$a$), resulting in a vanished drag term. Conversely, with sufficiently large 
$St$, where particle inertia becomes dominant, heavy particles tend to accumulate in regions with a radially inward velocity. It is notable that, for plane Couette flow, this biased sampling also occurs for large 
$St$, causing heavy particles to gather in low-speed streaks in the ejecting region (Rashidi, Hetsroni & Banerjee Reference Rashidi, Hetsroni and Banerjee1990; Eaton & Fessler Reference Eaton and Fessler1994; Marchioli & Soldati Reference Marchioli and Soldati2002).
In TC flow, as heavy particles cluster around regions with inward radial velocity (
$u_r<0$) more than those with outward velocity (
$u_r>0$), this results in a negative value of the biased sampling term. In figure 10(
$a$), we plot 
$-I_b$ in logarithmic scale vs radial position, also including the averaged fluid velocity vs radial position for reference in figure 10(
$b$). Indeed, in the bulk, the particles experience progressively stronger inward radial velocity for larger 
$St$, causing 
$I_b$ to monotonically decrease in the radial direction. This indicates that the biased sampling provides a net force that pushes the particles towards the inner walls.
The second phoresis integral addresses the turbophoresis pseudo-force, accounting for the migration of particles to regions with smaller radial velocity variance. With the no-slip and no-penetration boundary conditions imposed at the inner and outer walls, the radial velocity variance vanishes at both walls, leading to a tendency for particle migration toward the walls. In figure 11(
$a$), the turbophoresis integrals are plotted for different 
$St$, while the r.m.s. radial particle velocity (
$v_{r,rms}^{\prime }$) is also depicted in figure 11(
$b$). Notably, due to the vanished velocity variance at the two walls, a minimum 
$I_t$ is indeed observed in the bulk. It is apparent that the velocity variance declines at a greater rate at the outer wall than at the inner wall, resulting in a larger 
$I_t$ at the outer wall compared with the inner wall. This indicates a greater tendency to attract particles to the outer wall relative to the inner wall.
Finally, the radial dependency of the centrifugal effect (3.6) can be observed in figure 12(
$a$), showing a monotonic increase of 
$I_c$ along the radial direction, irrespective of 
$St$. The rotation of the inner wall causes particles to gain angular velocity, as depicted in figure 12(
$b$), illustrating the radial profiles of the azimuthal velocity of the particles. Consequently, particles exhibit a tendency to migrate outward in the domain due to the centrifugal effect.
We now comprehend the distinct roles of biased sampling, turbophoresis and the centrifugal effect in driving particles in TC flow: biased sampling leads to the migration of particles toward the inner wall, in contrast to the roles played by centrifugal and turbophoresis, where particles move toward the outer wall and both walls, respectively. To grasp the overall effect, we compare the relative dominance of the three mechanisms in particle concentration profiles, by examining the ratios of the intensities 
$I_b/I_c$ and 
$I_t/I_c$, as depicted in figure 13. From the figure, it is evident that the magnitude of both quantities is consistently less than one, signifying that the centrifugal effect holds the greatest strength at any radial position. For the smallest explored 
$St$ (
$=0.034$), the strength of biased sampling is nearly equivalent to that of the centrifugal effect, while the effect of turbophoresis is negligibly small. As a result, with two competing effects of similar strength, the particle concentration remains relatively uniform throughout the domain. However, the situation of relative dominance changes for heavy particles (large 
$St$). The relative strength of biased sampling progressively weakens with increasing 
$St$. Concurrently, turbophoresis becomes significant, with the magnitude of 
$I_t/I_c$ exceeding 60 % in the bulk and reaching almost 100 % near the outer wall. Indeed, a noticeable accumulation of heavy particles near the outer wall is observed due to the combined effect of the centrifugal effect and turbophoresis, pushing heavy particles toward that wall.
Figure 13. (
$a$) Radial profiles of the ratio of biased sampling integral to centrifugal integral; (
$b$) radial profiles of the ratio of turbophoresis integral and centrifugal integral.
4. Concluding remarks and outlook
To investigate the impact of particle inertia on the particle distribution in wall-bounded turbulent flows, we conducted DNS of particle-laden turbulent TC flow at 
$Re_i=2500$, introducing particles with various Stokes numbers. The simulations employed one-way coupled Lagrangian tracking of particles, with particle motions solely governed by the Stokes drag. Given the very low solid volume fraction considered, both two-way coupling and inter-particle collisions were disregarded. The focus of the study is on the spatial distribution of particles in TC flow. We have first examined the particle distribution near the outer cylindrical wall via a 2-D Voronoï diagram analysis given that particles tend to cluster around that wall for large enough 
$St$. Then, an approach utilizing the conservation equation for particle concentration has been employed, which was originally proposed by Johnson et al. (Reference Johnson, Bassenne and Moin2020). The benefit of this approach is to disentangle the effects caused by biased sampling, turbophoresis and centrifugal forces. The highlights of findings are below.
It is evident that, as particle inertia increases, a larger number of particles tend to accumulate in the outer-wall region, leading to an increase in preferential concentration near the outer wall. To quantify the extent of particle accumulation near the outer wall, we employed a 2-D Voronoï diagram analysis. While particles with low inertia exhibit a distribution close to random, an increase in particle inertia causes particles to cluster in the low-speed streaks at the outer wall, resulting in the formation of more clusters and voids. We further elucidated the emergence of stripe-like structures by showing the ratios of azimuthal length to axial length for the Voronoï cells. Additionally, by counting the number of particles in regions of low and high azimuthal velocity, we observed a preference for particles to accumulate in low-speed fluid structures.
Then, we differentiated the effects of biased sampling, turbophoresis and centrifugal forces by separately considering their relative strengths. By comparing the relative strength of these forces, we summarized the mechanisms of particle transport in TC flow in a schematic diagram depicted in figure 14. In terms of particle transport in the radial direction, centrifugal effects are predominant, with biased sampling effects being the second most influential for small inertial particles, while turbophoresis effects become second most important for large inertial particles. Each of these effects leads to distinct processes: centrifugal forces cause particles to migrate toward the outer wall; biased sampling forces cause particles to move toward the inner wall due to the fact that particles accumulate in the Taylor rolls with inward radial velocity; turbophoresis transports particles from the bulk region to the two walls due to the high turbulence intensity in the bulk region compared with that near the walls. However, it should be noted that the dominance of the centrifugal effect due to rotation is only valid on sufficiently large scales (at least within the range of parameters we have explored). This is because, on sufficiently large scales, particles behave more like tracers, where the inertia of the particles becomes negligible.
Figure 14. Schematic diagram showing the dominant terms in shaping the radial distribution of inertial particles at various 
$St$ in multiphase TC flow. The bluish and pinkish zones present the region of boundary layers and bulk, respectively. Three arrows of distinct colours represent the relative strength of the biased sampling, turbophoresis and the centrifugal effect. Here, the length of the arrow presents the intensity of the effect, and the directions of arrows represent the direction of particle migration led by the corresponding effect. The colour of the spheres represents different particle inertia, with the darker colours indicating particles with larger 
$St$ and lighter colours representing particles with smaller 
$St$.
Finally, there are several crucial considerations that were not addressed in this study, which merit further investigation. For instance, at higher particle mass loading, turbulence is notably influenced by particle feedback forces, thus affecting both the turbulent flow and the particle distribution. Besides, the range of density ratio of particles studied in this work is relatively narrow. For the case of low density ratio, the equation of particles with Stokes drag alone is not enough to describe the flow physics of particle movement. In the future, the range of particle density ratio could be expanded and more relevant forces should be considered, such as the added mass force and history force. In addition, only the leading-order translational dynamics of particles has been considered in this work. One can further examine the rotational dynamics of particles. The findings in this work provide groundwork for future study on the impact of inertial particles on the particle distribution in rotating wall-bounded turbulent flows.
Funding
This work was supported by the Natural Science Foundation of China under grant nos 11988102, 92052201, 12372219, 12421002 and 12072185.
Declaration of interests
The authors report no conflict of interest.
Appendix A. Derivation of the radial particle concentration profile
The wall-normal particle concentration profile in turbulent channel flow has been discussed by Johnson et al. (Reference Johnson, Bassenne and Moin2020). Here, we extend their proposed model to the radial particle concentration in TC flow with cylindrical coordinates. The derivation is based on the single-particle position–velocity p.d.f. along the radial direction. An additional centrifugal term emerges to account for the effect of rotation.
The radial direction single-particle position–velocity p.d.f. is denoted as 
$f(r, v_r, ; t)$
\begin{equation} f\left(r, v_r ; t\right)=\left\langle\delta(r-\hat{r}(t)) \delta\left(v_r-\hat{v}_r(t)\right)\right\rangle, \end{equation}
where 
$\delta (x)$ is the Dirac delta function and 
$\langle {\cdot } \rangle$ represents ensemble averaging over a disperse phase. Differentiating (A1) in time
\begin{equation} \frac{\partial f}{\partial t} = \left\langle \frac{\partial}{\partial t} \delta(r - \hat{r}(t)) \delta(v_r - \hat{v}_r(t)) \right\rangle. \end{equation}Using the chain rule for differentiation inside the average
\begin{equation} \frac{\partial f}{\partial t} = \left\langle -\dot{\hat{r}} \frac{\partial}{\partial r} \delta(r - \hat{r}(t)) \delta(v_r- \hat{v}_r(t)) - \dot{\hat{v}}_r \frac{\partial}{\partial v_r} \delta(r - \hat{r}(t)) \delta(v_r - \hat{v}_r(t)) \right\rangle, \end{equation}
and substituting the particle radial dynamic equation 
$\dot {\hat {r}}=\hat {v}_r; \dot {\hat {v}}_r=\hat {a}_r$,
\begin{equation} \frac{\partial f}{\partial t} ={-} \frac{\partial}{\partial r} \left\langle \hat{v}_r \delta(r - \hat{r}(t)) \delta(v_r - \hat{v}_r(t)) \right\rangle - \frac{\partial}{\partial v_r} \left\langle \hat{a}_r \delta(r - \hat{r}(t)) \delta(v_r - \hat{v}_r(t)) \right\rangle, \end{equation}
then we have the evolution of 
$f(r, v_r, ; t)$
\begin{equation} \frac{\partial f}{\partial t}+\frac{\partial\left(v_r f\right)}{\partial r}+\frac{\partial\left(\left\langle a_r | r, v_r \right\rangle f\right)}{\partial v_r}=0, \end{equation}
where 
$v_r$ is taken as the value of 
$\hat {v}_r$ at the point 
$r=\hat {r}$ and 
$v_r=\hat {v}_r$, 
$\langle a_r | r, v_r \rangle$ represents the conditional average acceleration given 
$r=\hat {r}$ and 
$v_r=\hat {v}_r$. Note that, for a complete expression, the collision term should also be included in (A5) to account for the change of velocities due to particle–particle collision. However, the contribution of this term to particle concentration vanishes due to momentum conservation of inter-particle collisions.
Since the particle phase in the turbulent TC flow is statistically homogeneous in the azimuthal and axial directions, the normalized particle concentration 
$C^*(r; t)$ can be obtained by integrating 
$f$ over the particle velocity 
$v_r$
\begin{equation} C^*(r;t) = C_0 \int_{-\infty }^{\infty} f(r, v_r ; t)\,\mathrm{d}v_r. \end{equation}
Here, 
$C_0$ is the bulk particle concentration. The particle radial direction momentum conservation is obtained as a first-order moment of 
$f$, i.e. multiplying (A5) by 
$v_r$ and 
$C_0$ and integrating over the particle radial velocity 
$v_r$
\begin{equation} \frac{\partial\left(\left\langle v_r | r\right\rangle C^*\right)}{\partial t}+\frac{\partial\left(\left\langle v_r^2 | r\right\rangle C^*\right)}{\partial r}-\left\langle a_r |r\right\rangle C^*=0. \end{equation}
The particle collisional term can be neglected because of particle momentum conservation by each collision event. In the statistically homogeneous state, the first term 
${\partial (\langle v_r | r\rangle C^*)}/{\partial t}=0$, and the particle radial momentum balance equation can be simplified as
\begin{equation} \frac{\mathrm{d}}{\mathrm{d} r}\left(\left\langle v_r^2 | r\right\rangle C^*\right)=\left\langle a_r | r\right\rangle C^*. \end{equation}
Substituting the particle radial acceleration 
$a_r = (u_r - v_r)/St +{v_\theta ^2}/{r}$ into the expression (A8), the particle radial momentum balance equation at the statistically homogeneous state can be rewritten as
\begin{equation} \left\langle v_r^2 | r\right\rangle \frac{\mathrm{d} C^*}{\mathrm{~d} r}=\left(\frac{\left\langle u_r | r\right\rangle}{St}-\frac{\mathrm{d}\left\langle v_r^2 | r\right\rangle}{\mathrm{d} r} + \frac{\left\langle v_\theta^2|r\right\rangle}{r} \right) C^*. \end{equation}
