1. Introduction
Particle-laden flows are ubiquitous in both nature and industrial applications. For example, in rivers, where the deposition of large grains can influence the solutal and nutrient exchange processes (Ferdowsi et al. Reference Ferdowsi, Ortiz, Houssais and Jerolmack2017). Or in oceans, where prediction for the accumulation of large plastic debris remains a topic of ongoing research (Cózar et al. Reference Cózar2014). In industrial applications the accumulation of particles in turbo-machineries can reduce the efficiency and even damage rotor or stator blades (Hamed, Tabakoff & Wenglarz Reference Hamed, Tabakoff and Wenglarz2006). Another example is in the paper-making industry, where the orientation of the fibres in the pulp suspension determines the mechanical strength of the final product (Lundell, Söderberg & Alfredsson Reference Lundell, Söderberg and Alfredsson2011). Given the importance of particle-laden flows, understanding phenomena such as transport and clustering is key for optimising engineering applications.
Studies of particle-laden flows (Salazar & Collins Reference Salazar and Collins2009; Toschi & Bodenschatz Reference Toschi and Bodenschatz2009; Balachandar & Eaton Reference Balachandar and Eaton2010; Mathai, Lohse & Sun Reference Mathai, Lohse and Sun2020), in general, have shown a rich phenomenology and can be broadly grouped into two categories. The first focuses on the flow response due to the presence of the particles. Most frequently, these particles are assumed to be spherical. Examples of these studies include investigations into the influence of particles on turbulent structures (Wang et al. Reference Wang, Abbas, Yu, Pedrono and Climent2018; Ardekani & Brandt Reference Ardekani and Brandt2019; Ramesh, Bharadwaj & Alam Reference Ramesh, Bharadwaj and Alam2019; Yu et al. Reference Yu, Xia, Guo and Lin2021), drag (Andersson, Zhao & Barri Reference Andersson, Zhao and Barri2012; Niazi Ardekani et al. Reference Niazi Ardekani, Costa, Breugem, Picano and Brandt2017; Bakhuis et al. Reference Bakhuis, Verschoof, Mathai, Huisman, Lohse and Sun2018; Spandan, Verzicco & Lohse Reference Spandan, Verzicco and Lohse2018; Wang, Xu & Zhao Reference Wang, Xu and Zhao2021) and the turbulent energy budget (Peng, Ayala & Wang Reference Peng, Ayala and Wang2019; Olivieri et al. Reference Olivieri, Brandt, Rosti and Mazzino2020). The second category focuses on explaining the dynamics of particles in these flows themselves. This becomes particularly interesting for non-spherical particles. For instance, on how shape influences particle rotation (Byron et al. Reference Byron, Einarsson, Gustavsson, Voth, Mehlig and Variano2015; Zhao et al. Reference Zhao, Challabotla, Andersson and Variano2015) or how particles cluster and preferentially align (Henderson, Gwynllyw & Barenghi Reference Henderson, Gwynllyw and Barenghi2007; Ni, Ouellette & Voth Reference Ni, Ouellette and Voth2014; Uhlmann & Chouippe Reference Uhlmann and Chouippe2017; Voth & Soldati Reference Voth and Soldati2017; Majji & Morris Reference Majji and Morris2018; Bakhuis et al. Reference Bakhuis, Mathai, Verschoof, Ezeta, Lohse, Huisman and Sun2019).
A review of the existing literature reveals that for spherical particles the particle dynamics in wall-bounded shear flows is reasonably well understood. However, apart from some recent work (Bakhuis et al. Reference Bakhuis, Mathai, Verschoof, Ezeta, Lohse, Huisman and Sun2019), little is yet known about the interactions of non-spherical particles in turbulent flows with curvature effects. The objective of this work is to fill this gap. Questions we ask are, for instance, how do non-spherical particles respond to shear flow with large-scale flow structures (due to different lift/drag forces), and do they exhibit preferential clustering or alignment? Another unresolved question regards the relation between the particle size compared with that of the flow features in the fluid phase. In particular, we want to work out the underlying physics. The answers to these questions can provide valuable insight into the underlying mechanism for particle-/bubble-induced drag reduction in wall-bounded turbulent flows, where particle geometry might affect momentum transport.
As a very well controlled shear flow, we choose Taylor–Couette (TC) flow and then study the two-way coupled dynamics of finite-size inertial anisotropic particles in this flow. TC flow is convenient for the following reasons. First, the flow regimes of TC flow are well understood and documented (Andereck, Liu & Swinney Reference Andereck, Liu and Swinney1986; Fardin, Perge & Taberlet Reference Fardin, Perge and Taberlet2014; Grossmann, Lohse & Sun Reference Grossmann, Lohse and Sun2016). Second, it is a closed system with exact balances (Eckhardt, Grossmann & Lohse Reference Eckhardt, Grossmann and Lohse2007) that is very accessible, both numerically and experimentally, due to its relatively simple geometry and symmetries.
To fully resolve the motion of the ellipsoidal particles and their two-way interaction with the surrounding fluid, we employ the immersed boundary method (IBM) using the moving least squares algorithm (Vanella & Balaras Reference Vanella and Balaras2009; de Tullio & Pascazio Reference de Tullio and Pascazio2016; Spandan et al. Reference Spandan, Meschini, Ostilla-Mónico, Lohse, Querzoli, de Tullio and Verzicco2017). The IBM is computationally advantageous for this application since the underlying grid is fixed and no computationally expensive remeshing is needed (Breugem Reference Breugem2012). Secondly, it is straightforward in the IBM to vary the size of the particles. The disadvantage of the IBM is that interparticle and particle–wall collisions need to be modelled. Here, we adopt the collision model of Costa et al. (Reference Costa, Boersma, Westerweel and Breugem2015), which has been widely tested and employed in studies for particle-laden flows (e.g. Niazi Ardekani et al. Reference Niazi Ardekani, Costa, Breugem, Picano and Brandt2017; Yousefi, Ardekani & Brandt Reference Yousefi, Ardekani and Brandt2020).
The present work is structured as follows. In § 2, we describe the TC set-up and give an overview of the investigated flow regimes. In § 3, we present the details of the numerical method describing the dynamics of the fluid and particles. In § 4, we present the spatial distributions of particles and categorise them into modes (i) to (iv). Then, we compare the modes for the two simulated particle sizes via the joint probability density function (p.d.f.) of the particle radial position versus the driving of the TC system. In § 5, we investigate the particle orientations with respect to the local cylinder tangent for the categorised modes. Here, we observe a strong particle alignment, which we correlate to the axial vorticity of the Taylor vortices. Finally, in § 7, we summarise our results.
2. Taylor–Couette set-up in the Taylor vortex flow regime
The TC set-up, as employed here, comprises a confined fluid layer between a coaxially rotating inner cylinder and a fixed outer cylinder (see figure 1a for a schematic). The dimensionless parameters characterising this system are the ratio of the inner radius 
$r_i$ and outer radius 
$r_o$ of the cylinders, i.e. 
$\eta \equiv r_i/r_o$, the aspect ratio of the domain, 
$\varGamma \equiv L/d$, and the Reynolds number of the inner cylinder, 
$Re_i \equiv r_i\omega _id/\nu$. Here, 
$L$ denotes the axial length of the cylinders, 
$d \equiv r_o - r_i$ the gap width, 
$\omega _i$ the angular velocity of the inner cylinder and 
$\nu$ the kinematic viscosity of the fluid. We set 
$\varGamma =2{\rm \pi} /3$ to allow one pair of Taylor vortices to fit within the domain (Ostilla-Mónico, Verzicco & Lohse Reference Ostilla-Mónico, Verzicco and Lohse2015) and 
$\eta =5/7$ to match the experimental T
$^3$C set-up (Bakhuis et al. Reference Bakhuis, Mathai, Verschoof, Ezeta, Lohse, Huisman and Sun2019). No-slip and impermeability boundary conditions are imposed on both cylinder walls. In the azimuthal and axial directions, periodic boundary conditions are used. We employ a rotational symmetry of order 6 to reduce computational cost such that the streamwise aspect ratio of our simulations 
$L_\varphi /d=(2{\rm \pi} r_i/6)/d=2.62$. The resulting streamwise domain length is sufficient to capture the mean flow statistics (Ostilla-Mónico et al. Reference Ostilla-Mónico, Verzicco and Lohse2015).
(a) Schematic of the TC configuration and geometrical definitions of the particle (not to scale). (b) The standard deviation of the normalised vertical velocity averaged over the domain and time versus 
$\textit {Ta}$.
The control parameter for the TC flow is 
$Re_i$ and we vary the values in the range 
$Re_i=[1.6\times 10^2,\ 8.0\times 10^3]$. The outer cylinder is fixed. For ease of comparison with existing numerical studies, we also define the Taylor number as
\begin{equation} \textit{Ta}\equiv\dfrac{(1+\eta)^6}{64\eta^4}Re_i^2, \end{equation}
with the corresponding values to the 
$Re_i$ range being 
$\textit {Ta}=[3.9\times 10^4,9.8\times 10^7]$. An overview of the cases is presented in table 1. This range of 
$\textit {Ta}$ covers the regimes known as Taylor vortex, wavy vortex and turbulent Taylor vortex flow (Grossmann et al. Reference Grossmann, Lohse and Sun2016). Within the chosen 
$\textit {Ta}$ range, the flow experiences a series of transitions. The lowest simulated 
$\textit {Ta}$ is chosen to lie slightly above the regime of circular Couette flow. The onset from circular Couette flow to Taylor vortex flow is estimated to occur at 
$\textit {Ta} \approx 1.0 \times 10^4$, which is determined from the critical Reynolds number 
$Re_{i,cr}(\eta )=(1+\eta ^2)/\{2\eta \alpha ^2[(1-\eta )(3+\eta )]^{1/2}\}$ with 
$\alpha = 0.1556$ for a resting outer cylinder (Esser & Grossmann Reference Esser and Grossmann1996). The transition point from Taylor vortex to wavy vortex flow has been investigated by numerous authors (Ahlers, Cannell & Lerma Reference Ahlers, Cannell and Lerma1983; Jones Reference Jones1985; Langford et al. Reference Langford, Tagg, Kostelich, Swinney and Golubitsky1988). Under current conditions it is empirically found to lie at around 
$\textit {Ta}=3\times 10^6$, by tracking the time- and volume-averaged standard deviation of the vertical velocity 
$u_z$ as a function of 
$\textit {Ta}$ (see figure 1
$b$). The visualisations of the aforementioned flow regimes are shown in figure 2.
(a–e) Instantaneous snapshots of the azimuthal flow field (arrows denote 
$u_r,u_z$) for various 
$\textit {Ta}$. (f–j) Corresponding time-averaged velocity fields. Here, the flow regimes are (a–c) Taylor vortex flow, (d) wavy vortex flow and (e) turbulent Taylor vortex flow.
Summary of simulation parameters. The first two columns denote the driving, expressed as either 
$\textit {Ta}$ or 
$Re_i$. The third column presents the grid resolution for 
$\ell /d=0.2$. The simulated cases corresponding to 
$\ell /d=0.1$ were performed on a 
$640\times 256 \times 480$ grid. The 
$\textit {Ta}$ 
$(Re_i)$ cases for which no 
$\ell /d=0.1$ were considered are indicated with (—). Here 
$0.1d/\eta _k$ and 
$0.2d/\eta _k$ denote the particle size to the Kolmogorov scale for 
$\ell /d=0.1$ and 
$\ell /d=0.2$, respectively.
3. Governing equations and numerical methods
3.1. Carrier phase
The velocity field is obtained by solving the incompressible Navier–Stokes equations in cylindrical coordinates. The continuity and momentum equations read (see e.g. Landau & Lifshitz Reference Landau and Lifshitz1987)
\begin{equation} \frac{1}{r} \partial_r(ru_r) + \frac{1}{r}\partial_\varphi u_\varphi +\partial_z u_z=0 \end{equation}and
$$\begin{gather} \partial_t u_\varphi + (\boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla}) u_\varphi + \frac{u_r u_\varphi}{r} ={-}\frac{1}{ r} \partial_{\varphi} p + \frac{1}{Re} \left(\boldsymbol{\Delta} u_\varphi + \frac{2}{r^2}\partial_\varphi u_r - \frac{u_\varphi}{r^2} \right) + f_\varphi, \end{gather}$$
$$\begin{gather}\partial_t u_r + (\boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla}) u_r - \frac{u_\varphi^2}{r} ={-}\partial_r p + \frac{1}{Re} \left( \boldsymbol{\Delta}u_r -\frac{2}{r^2}\partial_\varphi u_\varphi -\frac{u_r}{r^2} \right)+f_r, \end{gather}$$
$$\begin{gather}\partial_t u_z + (\boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla}) u_z ={-}\partial_z p + \frac{1}{Re} \boldsymbol{\Delta}u_z + f_z, \end{gather}$$
where 
$(\boldsymbol {u}\boldsymbol {\cdot }\boldsymbol {\nabla }) = u_r\partial _r + r^{-1}u_\varphi \partial _\varphi + u_z\partial _z$ and 
$\boldsymbol {\Delta } = r^{-1}\partial _r(r\partial _r) +r^{-2}\partial _{\varphi ^2} + \partial _{z^2}$. The last terms (
$f_\varphi$, 
$f_r$, 
$f_z$) on the right-hand side of (3.2a–c) denote the IBM forcing (see Appendix A for details).
We employ a fractional step method to numerically solve (3.1) and (3.2a–c). The velocity field is discretised using a conservative spatial, second-order, central finite-difference scheme and a temporal third-order Runge–Kutta scheme, except for the viscous terms that are treated implicitly with a Crank–Nicolson scheme. In the wall-normal direction, the grid is stretched with a clipped Chebychev type of stretching (cf. Ostilla-Monico et al. Reference Ostilla-Monico, Yang, van der Poel, Lohse and Verzicco2015). The grid is uniform in the azimuthal and axial directions. For more details, we refer the reader to Verzicco & Orlandi (Reference Verzicco and Orlandi1996). The grid resolution for the fluid phase is based on Ostilla et al. (Reference Ostilla, Stevens, Grossmann, Verzicco and Lohse2013), with the note that the grid aspect ratio here is 1.0 at mid-gap. This criterion is used to ensure sufficient nodes for the IBM, which is discussed in § 3.2. The time step is variable and satisfies the condition 
$\text {CFL}=0.3$; this restrictive limit is used owing to the explicit coupling of the particles to the fluid phase.
3.2. Particles
We use prolate ellipsoids as the dispersed phase in the TC flow. The control parameter for the particle is the ratio of particle size to the gap width, 
$\ell /d$, with 
$\ell$ the major axis of the ellipsoid (figure 1a). For our study, 
$\ell /d=0.1$ or 
$0.2$. The aspect ratio of the ellipsoid is 
$\varLambda \equiv \ell /b= 4$, with 
$b$ the minor axis of the particle. The exact value of 
$\varLambda =4$ was chosen to deviate significantly from spherical whilst maintaining a particle length that was small compared to the gap width as well as to maintain computational viability. Indeed a systematic study of 
$\varLambda$ may yield other interesting insights but this was outside the scope and computational costs of the present work. Sixteen particles are used in each simulation, yielding volume fractions of 
$0.01\,\%$ and 
$0.07\,\%$ for 
$\ell /d=0.1$ and 
$0.2$, respectively. The reported Stokes number is obtained via (cf. Voth & Soldati Reference Voth and Soldati2017)
\begin{equation} St\equiv\frac{\tau_p}{\tau_v},\quad \text{with}\ \tau_p= \frac{\varGamma}{18}\frac{b^{2}}{\nu}\frac{\varLambda\ln({\varLambda + \sqrt{\varLambda^2-1}})}{\sqrt{\varLambda^2-1}}\ \text{and}\ \tau_v=\dfrac{\nu}{u^2_\tau}. \end{equation}
Here, the density ratio is kept at 
$\varGamma =1$ and the friction velocity is defined as 
$u_\tau =\sqrt {\nu \partial _r \langle u_{\varphi }\rangle _{A,t}}|_{r_i}$ (average over time and inner cylinder).
The rigid particle dynamics is obtained by integrating the Newton–Euler equations. The governing equations and numerical methods regarding the particle translation, rotation and collisions are provided in Appendix A.
The ratio of the particle size to the Kolmogorov scale is estimated based on the global exact balance for TC flow (Eckhardt et al. Reference Eckhardt, Grossmann and Lohse2007)
\begin{equation} \eta_K/d = (\sigma^{{-}2} \textit{Ta} [Nu_\omega-1])^{{-}1/4},\quad \text{where}\ \sigma = (1+\eta)^4/(16\eta^2), \end{equation}
and the non-dimensional response parameter 
$Nu_\omega = J^\omega /J^\omega _{lam}$, the Nusselt number, which directly relates to the torque required to keep the angular velocity constant. The angular velocity flux is 
$J^\omega = r^3(\langle u_r \omega \rangle _{A,t} - \nu \partial _r \langle \omega \rangle _{A,t})$ and 
$J^\omega _{lam}= 2\nu \omega _i r_i^2r_o^2 /d^2$ its value for the laminar case (Eckhardt et al. Reference Eckhardt, Grossmann and Lohse2007). Here, we have used the Nusselt number as the dimensionless parameter to measure the transport of angular velocity flux because of the analogy between Rayleigh–Bénard and TC flow (see e.g. Grossmann et al. Reference Grossmann, Lohse and Sun2016). For all 
$\textit {Ta}$, we estimate 
$\ell /\eta _k$ based on (3.4) (see table 1).
To initialise the simulations, the single-phase flow is first advanced in time until a statistically stationary state is attained. Once this state has been achieved, the particles are inserted at random positions, with zero velocity within the domain, whilst ensuring that their initial distribution is spatially homogeneous. The initial orientations of the particles are also randomised. After inserting the particles, the simulations are run for at least 50 full inner cylinder rotations before collecting statistics. Between15 and 25 grid points per 
$\ell$ are used to ensure the boundary layers over the particles be sufficiently resolved. The particle Reynolds number is estimated as 
$Re_p\equiv \dot {\gamma } \ell ^2 /\nu$ and ranges from 
$O(0.1)$ to 
$O(60)$, with 
$\dot {\gamma }$ the average radial derivative of the azimuthal velocity in the bulk.
4. Spatial distributions of particles
4.1. Observed spatial modes
We examine the statistics of the particle positions. In particular, we select the cases 
$\textit {Ta}=3.9\times 10^4,\ 3.2\times 10^5,\ 1.9\times 10^6,\ 9.5\times 10^6$ and 
$9.8\times 10^7$ with particles (see also the single-phase flow fields in figure 2). In figure 3 we present the time-averaged particle distribution, after the initial transients, projected onto the 
$r$–
$z$ plane for both particle sizes: 
$\ell /d=0.1$ and 
$\ell /d=0.2$. Figure 3 reveals four distinct spatial patterns, which depend on both 
$\ell /d$ and 
$\textit {Ta}$. The characteristics of these flow different ‘modes’ are as follows.
a. Mode (i): steady large orbits (figure 3a,b,f).
b. Mode (ii): steady orbits, with particles spiralling closely around the vortex cores (figure 3d,h).
c. Mode (iii): a combination of modes (i) and (ii) (figure 3c,g).
d. Mode (iv): unsteady orbits, with particles distributed quite homogeneously throughout the domain (figure 3e,i,j).
Probability density function of the particle distribution. The average is taken over time and azimuthal direction. (a–e) Particles of size 
$\ell /d=0.1$ and (f–j) particles of 
$\ell /d=0.2$. The coloured circles on top of the contour plots denote distinguishable regimes of particle dynamics, which correspond to those in figures 4 and 6.
For mode (i), the rotational particle dynamics shows no stable alignment, but instead, a tumbling type of motion is observed. At this 
$\textit {Ta}$, the base flow is slightly above the transition point from the circular Couette flow to a Taylor vortex flow regime. We stress that the particles were released at random locations after the flow was fully developed – therefore, on release, each particle undergoes an inertial migration process before reaching its stable orbit.
For mode (ii), particle orbits are observed to cluster closely around the vortex cores and remain at a (small) finite distance away from the vortex core. The most pronounced example from figure 3 is at 
$\textit {Ta}=1.9\times 10^6$ for 
$\ell /d=0.2$ (figure 3h). Remarkably, the particle concentration at the core is much higher for the larger particles, thus indicating that this is definitely a finite-particle-size effect. Mode (ii) is accompanied also by a stable particle alignment, which is addressed in detail in § 5.
Mode (iii) consists of a combination of modes (i) and (ii). This regime is the transition between stable (limit-cycle-type) orbits and preferential clustering at the core. This mode is observed at 
$\textit {Ta}=1.0\times 10^6$ for 
$\ell /d=0.1$. Interestingly, for 
$\ell /d=0.2$, mode (iii) is observed at a lower value of 
$\textit {Ta}=3.2\times 10^5$. Hence, it appears that both particle clustering as well as the transitions to various particle orbit regimes are functions of 
$\ell /d$. In § 4.2, we study the regime transitions as a function of 
$\ell /d$.
One key feature of modes (i) to (iii) is the relative steadiness of the orbital regions as the particles trace their path through the domain. Note here that modes (i) to (iii) show orbits within a finite region of the domain and do not perfectly overlap each revolution (the pathline of a particle does not perfectly follow a streamline). The thickness of this stable region of the orbit is induced by multiple factors such as particle inertia and aperiodic angular dynamics (see § 5). Moreover, the orbital steady regions that a particle is observed to trace out are subject to the initial position of the particle. The steady orbital regions occur at a unique set of conditions and can be linked to two fundamental features of the Taylor vortices. Firstly, the background flow is completely steady and time-invariant. Secondly, the Taylor vortices are axisymmetric about the cylindrical axis (see § 2). The comparison with the background flow for modes (i) to (iii) is justified by noting that the flow state was not considerably altered by the particles. This observation is supported by the fact that 
$Nu_\omega$ was altered by just 
$0.02\,\%$ when ellipsoidal particles with 
$\ell /d=0.2$ are added to the flow at 
$\textit {Ta}=1.9\times 10^6$, which indicates that the overall transport properties and flow state are unaltered. This observation confirms that the influence of the particles on the flow is limited for the parameter range in which modes (i) to (iii) are observed.
For mode (iv), we now observe unsteady dynamics caused by unsteady Taylor vortices. For instance, as shown in figure 2(e), the flow corresponds to the wavy Taylor vortex regime. Due to the unsteadiness of the wavy Taylor vortices, particles experience spatial and temporal variations of the hydrodynamic loads. These variations prevent any stable orbits from happening. For 
$\textit {Ta}=9.5\times 10^6$ and 
$\ell /d=0.2$, particle trajectories tend to maintain some coherence (see figure 3i) and appear to trace out the complete vortex. However, for the same 
$\textit {Ta}$ and 
$\ell /d=0.1$ (figure 3e), as well as at larger 
$\textit {Ta}$ for 
$\ell /d=0.2$ (figure 3i), the coherence observed for mode (ii) is lost and the particles distribute nearly homogeneously throughout the domain. The number of collisions is very limited for mode (iv) due to the homogeneous distribution. The distributions for the latter combinations of 
$\textit {Ta}$ and 
$\ell /d$ are reminiscent of those for spheres and fibres in TC flow (Majji & Morris Reference Majji and Morris2018; Bakhuis et al. Reference Bakhuis, Mathai, Verschoof, Ezeta, Lohse, Huisman and Sun2019).
4.2. The transition from stable orbits to clustering at the core is size-dependent
In § 4.1, in some cases particles are observed to preferentially cluster at the vortex core. Now, we will examine the clustering behaviour in more detail. The 
$\textit {Ta}$ range for which clustering occurs is investigated via the p.d.f. of the particle radial distributions 
$P(r_p)$, with 
$r_p$ being the particle radial position. Distribution 
$P(r_p)$ is plotted versus 
$\textit {Ta}$ in figure 4. From this figure, two main observations can be made. (1) We observe that the modes, indicated by the colours at the top of the figures, do not occur exactly at the same 
$\textit {Ta}$ when comparing 
$\ell /d=0.1$ and 
$\ell /d=0.2$. (2) The magnitude of the peak is much more intense for 
$\ell /d=0.2$, suggesting that clustering is enhanced for the larger particles. These two effects are discussed below.
Joint p.d.f. of the particle normalised radial position versus 
$\textit {Ta}$. Particles with size ratio (a) 
$\ell /d=0.1$ and (b) 
$\ell /d=0.2$, and (c) spheres with identical volume to ellipsoids with 
$\ell /d=0.2$. For reference, the modes addressed in § 4.1 are indicated with (i)–(iv) in the coloured top bar and the colour coding corresponds to the one used in figures 3 and 6.
In § 4.1 we defined four modes characterised by the particle distributions in the flow field. We highlighted these regimes in figure 3 with four colours at the top; the colours correspond to those used in figure 4. Mode (i) corresponds to the stable particle orbits, resulting in helical particle trajectories. For this regime, we observe a preferential particle concentration away from the vortex centre, as is evident by the light-blue regions at 
$(r_p-r_i)/d$ around 
$0.2$ and 
$0.8$ in figure 4. In this regime, particles in the Taylor vortex are found to move further outwards. Beyond this regime for even larger 
$\textit {Ta}$, we observe mode (iii), which is the mixed regime in which both modes (i) and (ii) are observed simultaneously. Beyond this, mode (ii) is encountered; particles cluster at the central region of the Taylor vortices as evidenced by the peak in the p.d.f. at 
$(r_p-r_i)/d \approx 0.5$. On average we found for mode (ii) that a particle at 
$\ell /d=0.2$ translates around one 
$D_{\textit {eq}}$ per 10 full inner cylinder rotations for mode (ii). Here, 
$D_{\textit {eq}}$ refers to the equivalent diameter of a sphere. For the particles of size 
$\ell /d=0.1$ the translation time increased to one 
$D_{\textit {eq}}$ per 20 full inner cylinder rotations. Finally, beyond this regime, particles move outwards again and distribute more homogeneously when the Taylor vortex starts to destabilise into the wavy regime. When comparing figures 4(a) and 4(b), the transition between the regimes shifts to lower 
$\textit {Ta}$ for larger particles. Since there is this clear size dependence, it is, therefore, instructive to compare the corresponding particle time scale 
$\tau _p$, which is set by the particle size, with the fluid shear time scales, 
$\tau _\nu$: effectively, we compute the Stokes number 
$\textit {St}$ for the particles as defined earlier in (3.3). When considering 
$\textit {St}$ as the governing parameter, we observe that mode (iii) (transition to clustering) occurs at 
$\textit {St} \approx 0.7$ for 
$\ell /d=0.1$ and 
$\textit {St}\approx 1.0$ for 
$\ell /d=0.2$, suggesting that 
$\textit {St}$ is of a similar order of magnitude for mode (iii). However, the reason why we urge caution is that, aside from the small numerical parameter space, TC flows are inherently three-dimensional because of curvature effects. Therefore, particle dynamics becomes sensitive to spatially varying hydrodynamic forces (Trevelyan & Mason Reference Trevelyan and Mason1951). Curvature effects on particles in TC flow can be straightforwardly investigated by examining Jeffery solutions (Jeffery Reference Jeffery1922) for prolate ellipsoids in the Couette regime. In Appendix B evidence is provided that curvature effects on the particle rotation rate are already observed at Taylor numbers as low as 
$\textit {Ta}=1.0$.
Comparing clustering for the two different particle sizes, we observe a much higher magnitude of 
$P(r_p)$ in figure 4 for 
$\ell /d=0.2$ than for the case 
$\ell /d=0.1$. The clustering is weaker for smaller 
$\ell /d$ for the following reason: the clustering regime of 
$\ell /d=0.1$ falls together with the onset of the wavy Taylor vortex regime, while particles of 
$\ell /d=0.2$ start to cluster around 
$\textit {Ta}\approx 4\times 10^5$ (Taylor vortex regime).
4.3. A comparison with volume-equivalent spheres
In § 4.2, we found clustering to be much more pronounced for 
$\ell /d=0.2$ compared with 
$\ell /d=0.1$. Here, we ask whether the clustering is due to the larger size or to a specific feature of the ellipsoidal particle shape. To clarify this issue we simulate 16 spheres with the same volume as ellipsoids with 
$\ell /d=0.2$, using the same initialisation procedure and simulation settings. Interestingly, 
$P(r_p)$ does not indicate strong clustering for the spheres (cf. figure 4c). Note that the obtained results for spheres closely resemble the experimental work by Majji & Morris (Reference Majji and Morris2018). However, we emphasise that in the range of 
$\textit {Ta}=3.2\times 10^5$ to 
$\textit {Ta}=1.9\times 10^6$ the orbit of few spheres is close to the vortex core, which results in a p.d.f. that is comparable with the one observed in figure 3(b) where most particles spiral to the edge of the vortex. Since the spheres have comparable 
$\textit {St}$ with the ellipsoids with 
$\ell /d=0.1$, we conclude that the particle shape has a pronounced influence on the p.d.f. of the particle positions.
5. Statistics of the particle orientation
5.1. Angular dynamics
Up to this point, the spatial statistics of particles have been examined. A number of regimes were found, one of which is of specific interest since ellipsoidal particles were found to cluster at the central region of Taylor vortex cores. Additionally, clustering is observed to enhance when the particle size is increased from 
$\ell /d=0.1$ to 
$\ell /d=0.2$. As a follow-up, we examine the statistics of particle orientations, corresponding to the identified spatial modes. We examine the angle 
$\theta _z$ (see figure 5a) between the particle pointing vector 
$p_i$ and the local tangent of the cylinder (cf. Bakhuis et al. Reference Bakhuis, Mathai, Verschoof, Ezeta, Lohse, Huisman and Sun2019). Here, we make use of the symmetry of the particle and let 
$\theta _z \in [-{\rm \pi} /2,{\rm \pi} /2]$.
(a) Definition of the angle 
$\theta _z$ between the pointing vector 
$p_i$ and the tangent along the cylinder. By symmetry of the particle, 
$\theta _z \in [-{\rm \pi} /2,{\rm \pi} /2]$. (b) Angular time signal of a particle within a stable orbit (light-blue line, 
$\textit {Ta}=3.9\times 10^4$; cf. figure 3f) and for mode (ii) (yellow line, 
$\textit {Ta}=1.9\times 10^6$; cf. figure 3h). (c) Definition of the width of the p.d.f. 
$P(\theta _z)$. The width is measured for the highest peak of 
$P(\theta _z)$ at half-height.
Two typical time signals of 
$\theta _z$ are given in figure 5(b) for particles of size 
$\ell /d=0.2$. The signal for 
$\textit {Ta}=3.9\times 10^4$ belongs to a particle travelling along a stable orbit. Interestingly, the particle orientational dynamics in the steady Taylor vortex regime still shows a Jeffery-like character. These Jeffery-like features have been observed in a variety of cases that also do not strictly fulfil the criteria of Jeffery's equations (Wang et al. Reference Wang, Abbas, Yu, Pedrono and Climent2018; Kamal, Gravelle & Botto Reference Kamal, Gravelle and Botto2020). In contrast, the time signal of 
$\theta _z$ for 
$\textit {Ta}=1.9\times 10^6$ shows an interesting, nearly constant angle 
$\theta _z$. This may be visualised as if the axis of revolution always aligns with the local cylinder tangent. The particle does not flip but oscillates only relatively mildly. This case corresponds to preferential particle concentration at the vortex core (mode (ii); see figure 3h).
For illustration purposes eight particles of size 
$\ell /d=0.4$ at 
$\textit {Ta}=1.9\times 10^6$ are simulated (subject to same grid resolution for the corresponding case 
$\ell /d=0.2$). This larger particle displays mode (ii), too, with a slightly larger oscillatory character (see figure 5b). This indicates that finite-size effects remain visible for even larger particles, but nonlinear effects come into play, which are out of the scope of this study.
The statistics of 
$\theta _z$ are discussed in the following and linked to the spatial distribution regimes of § 4.1.
5.2. Angular statistics corresponding to the observed spatial distributions
Typical p.d.f.s of 
$\theta _z$ for various 
$\textit {Ta}$ and for 
$\ell /d=0.1$ and 
$\ell /d=0.2$ are given in figure 6(a,b). The shown cases correspond to the spatial distributions in figure 3. Several interesting features can be observed in the distribution of 
$P(\theta _z)$. In particular, we find that these features correlate to the different spatial particle distributions described earlier in § 4.1 as listed below.
a. Mode (i): steady large orbits. A positive preferred orientation (maximum of
$P(\theta _z)$ occurs at $\theta _z>0$).b. Mode (ii): orbits spiralling closely around the core. A sharp peak for
$P(\theta _z)$ located at $\theta _z\approx 0$.c. Mode (iii): a combination of modes (i) and (ii). Angular dynamics show modes (i) and (ii) on the border between clustering and stationary orbits.
d. Mode (iv): unsteady dynamics, particles distribute throughout the whole domain. A non-homogeneous distribution of
$\theta _z$ with the maximum of $P(\theta _z)$ occurring at negative $\theta _z$.
The p.d.f. of 
$\theta _z$ for various 
$\textit {Ta}$. (a) Particles with size 
$\ell /d=0.1$ and (b) particles with size 
$\ell /d=0.2$. (c) Decomposition of the orientational statistics for mode (iii), showing the origin of the two peaks. A fraction of the particles is close to the vortex core, whereas other particles are within a stable orbit. For reference, the experimental observations from Bakhuis et al. (Reference Bakhuis, Mathai, Verschoof, Ezeta, Lohse, Huisman and Sun2019) are added (cf. 
$\textit {Ta}=9.5\times 10^{10}$). The colour coding of the plots corresponds to that of figures 3 and 4.
For mode (i), the flow is in the steady Taylor vortex flow regime (
$\textit {Ta}=3.9\times 10^4$). From figure 6, for both 
$\ell /d=0.1$ and 
$\ell /d=0.2$ the angular statistics display a non-homogeneous distribution, with a positive preferred angle. For mode (ii), particles agglomerate at the Taylor vortex cores. Remarkably they show a strong preferential alignment at 
$\textit {Ta}=1.9\times 10^6$ as confirmed by the sharp peaks of the alignment probability of figure 6. The more defined peak of 
$P(\theta _z)$ observed for 
$\ell /d=0.2$, as compared with 
$\ell /d=0.1$, is related to the enhanced clustering (see § 4.2) and correlates with the result of stronger preferential alignment, thus indicating that the particle size plays a predominant role in the phenomenon. Tails of 
$P(\theta _z)$ can also be observed, for example at 
$\textit {Ta}=1.9\times 10^6$ in figure 6(a,b). These tails are the result of particles precessing in a stable orbit visible in the particle distribution in figure 3(d,h). Our explanation as to why particles may still intermittently exhibit behaviour akin to ‘stable orbits’ is because particles collide throughout these simulations at this 
$\textit {Ta}$ (
$O(5)$ for all particles per full inner cylinder rotation). Observations from such events indicate that collisions cause the particles to end up further away from the vortex core in a meta-stable orbit, which eventually decays to the stable preferential alignment at the vortex core.
For mode (iii), two peaks are observed for 
$P(\theta _z)$, as shown by the green curves in figure 6(a,b). These two peaks originate from the two spatial modes (i) and (ii), shown in figure 3(g). The contribution from the two modes can be explained by disentangling the two angular statistics, which we illustrate in figure 6(c). First, we separate the two spatial modes by taking a subsample of particles close to and far from the vortex cores. Next, the angular statistics of these subsamples are computed, resulting in two p.d.f.s of 
$\theta _z$ (black dashed lines in figure 6c). These separate p.d.f.s illustrate that particles close to the vortex core show preferential alignment at 
$\theta _z\approx 0$, whereas those in a stable orbit far from the core peak at 
$\theta _z\approx 0.09{\rm \pi}$. It is highlighted that this particle behaviour forms the transition point between stable orbits and clustering (mode (iii) in figure 4), and occurs within a very narrow 
$\textit {Ta}$ range.
For mode (iv), particles distribute homogeneously throughout the domain due to the instabilities and unsteadiness of the flow. The preferential alignment observed with axisymmetric Taylor vortices (mode (ii)) cannot exist when the flow undergoes a transition to wavy Taylor vortex flow. For 
$P(\theta _z)$ this results in a distribution that is flatter. However, some statistical preferential alignment persists. Bakhuis et al. (Reference Bakhuis, Mathai, Verschoof, Ezeta, Lohse, Huisman and Sun2019) reported a difference of 40 % between the lowest and highest values of 
$P(\theta _z)$ at 
$\textit {Ta}=1.0\times 10^{12}$ for cylinders of aspect ratio 
$\varLambda =5$. Within this work, at 
$\textit {Ta} = 9.8 \times 10^8$, the difference is about 67 % which suggests that the tendencies for particles to preferentially align are stronger at lower 
$\textit {Ta}$. We also highlight that, while the 
$\textit {Ta}$ values are much lower in our set-up than in Bakhuis et al. (Reference Bakhuis, Mathai, Verschoof, Ezeta, Lohse, Huisman and Sun2019), there is a general tendency for the peaks to shift for increasing 
$\textit {Ta}$ towards negative 
$\theta _z$ and lower 
$P(\theta _z)$. Indeed, the incipient trend is consistent with the distributions in Bakhuis et al. (Reference Bakhuis, Mathai, Verschoof, Ezeta, Lohse, Huisman and Sun2019). Further studies at larger 
$\textit {Ta}$ in the simulations will be necessary to verify this trend.
5.3. The most pronounced alignment of particles
In § 5.2, preferential alignment of the particle angle 
$\theta _z$ is observed in the case when particles agglomerate near the vortex core. Here, the objective is to determine the conditions for which alignment is strongest. For all 
$P(\theta _z)$, the width 
$w$ of the p.d.f. around the highest peak is calculated and taken at half-height (see figure 5c). The plot of 
$w$ versus 
$\textit {Ta}$ is given in figure 7. Remarkably, 
$w$ has a very pronounced minimum (note the double-logarithmic scale) with respect to 
$\textit {Ta}$, occurring at around 
$\textit {Ta}=7\times 10^5$ for 
$\ell /d=0.1$ and at 
$\textit {Ta}=4\times 10^5$ for 
$\ell /d=0.2$. This minimum corresponds to a particle that oscillates the least with respect to the local tangent of the cylinder. Intriguingly, the minimum is even more pronounced for 
$\ell /d=0.2$ compared with 
$\ell /d=0.1$. In the remainder of this work, we discuss the origin of this minimum in 
$w$ and offer an explanation for why this is observed in the investigated configuration.
Width, 
$w$, of the p.d.f. 
$P(\theta _z)$ versus 
$\textit {Ta}$. The definition of 
$w$ is sketched in figure 5.
6. The role of axial vorticity at the vortex core
6.1. The link between strong alignment and minimum axial vorticity
From our analysis in the preceding sections, we observed that at specific 
$\textit {Ta}$ values particles tend to get trapped within the vortex core and have a preferred orientation. Now, we aim to answer the question: Why do the particles align at this 
$\textit {Ta}$ value? In the following, we show that the preferential alignment is linked to the TC flow state which exhibits a minimum in the shear gradient at the Taylor vortex core. In particular, we base our analysis on the axial component of the vorticity of the flow, which is shown to be the key metric determining the preferential alignment.
In the spirit of Jeffery's equation for the rotation of an ellipsoid, we formulate an area average of the axial vorticity, 
$\omega _z=r^{-1} \partial _{r} ( r u_\varphi ) - r^{-1}\partial _\varphi u_r$ evaluated in the 
$r,z$ plane, based on the premise that a particle aligns with the local cylinder tangent. This vorticity component is held responsible for rotational flipping events of the ellipsoidal particle around the 
$z$ axis. To compute the axial vorticity, it is first instructive to identify a region of interest where particles would presumably cluster. Taking heed from the clusterings observed in figure 3, this region can be reasonably and safely assumed to coincide with the central regions of the Taylor vortices. Therefore, the average of 
$\omega _z$ is taken for each 
$r$–
$z$ plane over a circular patch positioned at the vortex cores. The patch has a radius assigned with cross-section dimensions similar to those of the particle. Here, we select a circle with radius 
$i\cdot b$, with 
$b$ the minor axis of the particle. The Taylor vortex cores are identified by a local minimum of the meridional velocity, 
$u_r^2+u_z^2$. Note that for the single-phase flow in the Taylor vortex regime, the vorticity in the 
$r,z$ plane is stationary and, due to the axial symmetry of the flow, independent of the coordinate 
$\varphi$. The averaged 
$\omega _z$ is normalised by the rotational velocity of the inner cylinder 
$\omega _i$ and plotted in figure 8 versus 
$\textit {Ta}$.
The average vorticity, 
$\langle | \omega _z | \rangle /\omega _i$, is computed for the single-phase flow situation. The area covered for the average vorticity, 
$\omega _z$, is a circle centred at the vortex core with radius 
$b$ and 
$2b$, where 
$b$ denotes the particle minor axis with 
$\ell /d=0.2$. The analysis for a patch with radius 
$3b$ is performed in the presence of particles of size 
$\ell /d=0.2$. The inset shows an instantaneous 
$\varphi$–
$r$ slice of 
$\omega _z$ for the single-phase flow and two-phase flow cases, highlighting the perturbed vorticity fields due to the presence of particles.
A minimum of 
$\langle |\omega _z|\rangle$ can be clearly observed at 
$\textit {Ta} \approx 6\times 10^5$. This minimum implies that, roughly close to this 
$\textit {Ta}$ value, the fluid torque applied to the body (following Jeffery's equations) will be the smallest. Indeed, comparing this minimum with the previously acquired most pronounced alignment in figure 7, a very good agreement is found: the preferential alignment of the particle at 
$\textit {Ta}=4.2\times 10^5$ and the minimum of average vorticity for the single-phase flow occurs at 
$\textit {Ta}=6\times 10^5$. Additional calculations with a larger nominal diameter equal to 
$2b$ (red circles in figure 8) also find the minimum to occur around this 
$\textit {Ta}$, which shows that this metric is quite robust.
For completeness, we further analyse the effect on 
$\langle |\omega _z|\rangle$ with particles for 
$\ell /d=0.2$. In contrast to the single-phase flow, the presence of particles introduces velocity gradients in their vicinity because of the formation of viscous boundary layers at the particle surface. Therefore, 
$\langle |\omega _z|\rangle$ is determined at a slightly larger patch with radius 
$3b$ in order to filter out spurious vorticity magnitudes with particles. Here, the volume occupied by the particle is excluded from the calculation. Patches smaller than 
$3b$ are possible although they result in larger variances due to the closer proximity to the wakes shed by the particles. This trend of 
$\langle |\omega _z|\rangle$ including the particles is shown in figure 8 with black symbols. As can be seen from the figure, the trend is closely followed, but when particles agglomerate near the vortex core, the metric becomes very sensitive to particle-induced gradients and skews the picture. We find the lower values of 
$\langle |\omega _z|\rangle$ to occur in parts of the domain in which the particle is absent. The pronounced results of 
$\langle |\omega _z|\rangle$ for the single-phase flow served as a good guideline in our understanding of strong alignment. However, distilling similar results with particles is challenging.
6.2. Lagrangian statistics of particle rotational energy
In this analysis we investigate the effect of 
$\omega _z$ on the particle dynamics in relation to the particle position. Here, we select a single representative particle (
$\ell /d=0.2$) from cases 
$\textit {Ta}=1.0\times 10^5$ and 
$\textit {Ta}=5.6\times 10^5$. For the purpose of our discussion we refer to these as case 1 and case 2, respectively. Note that each case is highlighted with a circle at the top of figure 8. Case 1 corresponds to a particle exhibiting a final spatial distribution with a stable orbit far from the vortex core, whereas case 2 converges to a strong alignment mode in the vicinity of the vortex core. For both cases the rotational kinetic energy is tracked over time (translational energy is left out of the analysis). Here, because of the tumbling motion of the particle, it is assumed that its rotational energy provides a reasonable metric of the particle Lagrangian dynamics. The particle selected for case 1 initially started out close to the vortex core, whereas that for case 2 initially started at the edge of the vortex and subsequently spiralled inwards.
The particle rotational energy, 
$E_r$, is given by
\begin{equation} E_r = \tfrac{1}{2}\hat{\boldsymbol{\omega}}^T \boldsymbol{\mathsf{I}}_p \hat{\boldsymbol{\omega}}, \end{equation}
with 
$\hat {\boldsymbol {\omega }}=\boldsymbol {\omega }_p/\omega _i$ the normalised rotational velocity and 
$\boldsymbol{\mathsf{I}}_p$ its moment of inertia tensor of the particle. The two cases are illustrated in figure 9. We observe local regions where the rotational kinetic energy is highest, which correspond to particle tumbling events. The distinct difference between cases 1 and 2 is when 
$E_r$ is examined at the vortex core. In fact, for case 2, the value of 
$E_r$ at the core is negligibly small and is well correlated with the low 
$\langle |\omega _z| \rangle$ value shown in figure 8. This local minimum, interestingly, only holds for a specific region close to the core, whereas outside of the core, the rotational kinetic energy is finite. This picture, therefore, illustrates that strong alignment occurs only when particles are within the local vorticity minimum of the core. In contrast, case 1 clearly shows tumbling events that can still be found at the vortex core. This phenomenon, therefore, illustrates the importance of axial vorticity in determining the preferential alignment of ellipsoids in TC flow.
The dimensionless space–time evolution of the rotational energy, 
$E_r$, of a particle (
$\ell /d=0.2$) for (a) 
$\textit {Ta}=1.0\times 10^5$ (case 1) and (b) 
$\textit {Ta}=5.6\times 10^5$ (case 2). In (a), the particle eventually spirals outwards and does not display mode (ii). In (b), the particle spirals inwards towards the core. The starting and ending positions of the particle are denoted by 
$\times$ and 
$\bigcirc$, respectively. Large magnitudes of 
$E_r$ correspond to tumbling events of the particle. The arrows denote the 
$(u_r,u_z)$ velocity field. (c) Accelerations 
$\ddot {\boldsymbol {x}}$ for case 2 relative to the single-phase case. The black arrows indicate the single-phase fluid accelerations 
$\dot {\boldsymbol {u}}$.
Based on our findings for the vorticity statistics, this secondary mechanism appears also to be a crucial factor that determines preferential alignment, in addition to the Stokes number effect discussed in § 4.2. Simply put, while clustering can be controlled by varying 
$\textit {St}$ (or 
$\ell /d$), the unique nonlinear axial vorticity fields at different 
$\textit {Ta}$ establish a nominal limit for which preferential alignment can occur for a given 
$\ell /d$. We emphasise that this mechanism is present only in the regime of Taylor vortex flow.
6.3. The clustering mechanism for ellipsoids
The observation that the particle alignment and clustering depend on the particle geometry (§ 4.3) raises the important question: what drives the migration of these ellipsoids? Since we observe the migration to the centre of the vortex core only for prolate ellipsoids, we know it is a geometrical effect. To further understand the underlying physics, we hypothesise that the particle acceleration 
$\ddot {\boldsymbol {x}}$ closely resembles the single-phase fluid acceleration 
$\dot {\boldsymbol {u}}$ (cf. (3.2a–c)). This assumptions is based on the observation that the particles have negligible influence on the flow (see § 4.1). In figure 9(c), we present the particle accelerations 
$\ddot {\boldsymbol {x}}$ and single-phase fluid accelerations 
$\dot {\boldsymbol {u}}$ (black arrows) for case 2. We find that 
$\langle (\ddot {\boldsymbol {x}} - \dot {\boldsymbol {u}})/\dot {\boldsymbol {u}} \rangle \approx 10\,\%$, confirming the strong correlation between 
$\ddot {\boldsymbol {x}}$ and 
$\dot {\boldsymbol {u}}$. Interestingly, from figure 9(c), the particle experiences a net acceleration towards the core in the plume ejecting area (around 
$z/d\approx (1.2,\,1.5)$, towards the outer cylinder). However, at the same time, the analysis also shows the particle to accelerate away from the core in the plume impacting area (
$z/d\approx (1.5,\,2)$, towards the inner cylinder). This seems to be the general picture for all steady vortex flow cases; a fine balance appears between the two respective effects, which depend on 
$\textit {Ta}$. The particle Reynolds number, 
$Re_p\equiv \ell \|\dot {\boldsymbol {x}}-\boldsymbol {u}\|/\nu$, depends on 
$\textit {Ta}$ and the location with respect to the vortex. For case 2 we find 
$Re_p$ of 
$O(0.1)$ at the vortex core and 
$O(1)$ at the edge of the vortex. Interestingly, the region where the particle shows a net acceleration towards the vortex core overlaps with the region of tumbling events induced by inertial effects (see § 6.2). Therefore, we note that the clustering may indeed emerge from tumbling events.
7. Conclusion
In this study, we investigated finite-size, neutrally buoyant, prolate ellipsoids (aspect ratio 
$\varLambda =4$) in TC flow. The explored flow regimes are governed by pure inner-cylinder driving and comprise the following regimes: Taylor vortex, wavy vortex and turbulent Taylor vortex flow. The fluid phase is simulated using direct numerical simulation, whereas particles are represented through an IBM approach. Two particle size ratios were considered: 
$\ell /d=0.1$ and 
$\ell /d=0.2$ for volume fractions 0.01 % and 0.07 %, respectively. Here, 
$\ell$ denotes the particle major axis and 
$d$ the gap width.
Upon releasing the particles at initially random location and orientation, we observe, after a transient, various distinctive particle distributions (figure 3). These distributions are according to the flow regime and particle spatial distributions categorised in modes (i) to (iv) (see § 4.1). Modes (i) to (iii) are observed in the Taylor vortex flow regime. Here, mode (i) corresponds to steady large orbits away from the core. Remarkably, for higher 
$\textit {Ta}$ in the Taylor vortex flow regime, particles get trapped in the vortex core. This particle distribution is denoted as mode (ii) and it is the focus of this work. Interestingly, the 
$\textit {Ta}$ range corresponding to mode (ii) is different for 
$\ell /d=0.1$ and 
$\ell /d=0.2$. Moreover, the particle concentration for mode (ii) was observed to be much higher for 
$\ell /d=0.2$ than for 
$\ell /d=0.1$ (figure 4). Mode (iii) is a transition in which mode (i) as well as mode (ii) are observed. Mode (iv) corresponds to particle distributions in the wavy vortex regime and turbulent Taylor vortex regime. Here, particles distribute throughout the domain due to the instabilities in the flow.
Furthermore, we find distinctive particle orientations for each mode. Let 
$\theta _z$ denote the angle between the particle axis of revolution and the local cylinder tangent. We find for mode (ii) a sharp peak around 
$\theta _z=0$ in the p.d.f. 
$P(\theta _z)$ (figure 7). The ability of particles to align is found to depend on three factors. Firstly, the gradient in the flow. We find the most pronounced alignment for particles with 
$\ell /d=0.2$ to occur at 
$\textit {Ta}=4.2\times 10^5$. This was observed to closely match the axial vorticity minimum at the vortex core (
$\textit {Ta}=6\times 10^5)$. In comparison, the axial gradient at 
$\textit {Ta}=1\times 10^5$ (mode (i)) is observed to be two orders of magnitude higher. From the particle Lagrangian dynamics a stable alignment was not observed (cf. figure 9a). Secondly, the ability of particles to cluster. Figure 9(b) indicates that a stable particle alignment (absence of rotational energy) only occurs once a particle is in the near vicinity of the vortex core. Lastly, the onset of instabilities in the flow. We observed in the transition from steady to wavy vortex flow that particles spread throughout the domain (cf. figure 3e,i,j). The corresponding spatial distributions (mode (iv)) become significantly flatter than the ones of mode (ii) (cf. figure 6).
Our results collectively indicate that shear, large-scale structures induced by the curvature of the domain, particle shape and particle size play an important interlinked role in the dynamics of the particles themselves. The TC flow is known for its rich variety of flow structures. By adding particles, we observed that the particle dynamics alters significantly when the driving parameter 
$\textit {Ta}$ is varied. Furthermore, we expect the effect of particle anisotropy, as expressed here by the parameter 
$\varLambda$, on particle dynamics to be of importance for the clustering behaviour and alignment statistics. It is suspected that oblate ellipsoids (
$\varLambda < 1$) will not experience strong alignment in the vortex cores, since in that case an alignment of the pointing vector with the direction of the vorticity will lead to enhanced tumbling behaviour. A study of the effects of anisotropy would be an interesting extension of the present work.
Acknowledgements
We wish to express our gratitude to H.W.M. Hoeijmakers for helpful comments on our manuscript. We acknowledge that the results of this research have been achieved using the DECI resource Kay based in Ireland at the Irish Center for High-End Computing (ICHEC) with support from the PRACE aisbl. This work was partly carried out on the national e-infrastructure of SURFsara, a subsidiary of SURF cooperation, the collaborative ICT organisation for Dutch education and research.
Funding
This work is part of the research programme of the Foundation for Fundamental Research on Matter with project number 16DDS001, which is financially supported by the Netherlands Organisation for Scientific Research (NWO). R.J.A.M.S. acknowledges the financial support from ERC (European Research Council) Starting Grant No. 804283 UltimateRB.
Declaration of interests
The authors report no conflict of interest.
Appendix A. Numerical method particles
A.1. Newton–Euler equations
The dynamics of rigid particles is governed by the Newton–Euler equations:
\begin{equation} \left.\begin{array}{c} \displaystyle\dfrac{1}{r_p}\dfrac{{\rm d}}{{\rm d}t}(r_p^2\dot{\varphi}_p)=\dfrac{6}{\rm \pi} \int (\boldsymbol{\tau}\boldsymbol{\cdot}\boldsymbol{n})_{\varphi} \, {\rm d}S + F_{\varphi},\\ \displaystyle\ddot{r}_p =\dfrac{6}{\rm \pi} \int (\boldsymbol{\tau}\boldsymbol{\cdot}\boldsymbol{n})_{r} \, {\rm d}S + r_p \dot{\varphi}^2_p + F_{r},\\ \displaystyle\ddot{z}_p =\dfrac{6}{\rm \pi} \int (\boldsymbol{\tau}\boldsymbol{\cdot}\boldsymbol{n})_{z} \, {\rm d} S + F_{z}. \end{array}\right\} \end{equation}
Equations (A1) are non-dimensionalised with the length scale 
$D_{\textit {eq}}$ (volume-equivalent diameter of a sphere) and velocity of the inner cylinder 
$\omega _ir_i$. The terms 
$F_\varphi$, 
$F_r$ and 
$F_z$ denote the collision force due to short-range particle–particle and particle–wall interactions (see Appendix A.2 for details).
The term 
$\int (\boldsymbol {\tau }\boldsymbol{\cdot }\boldsymbol {n})_{i} \, \textrm {d}S$ is computed as
\begin{equation} \int (\boldsymbol{\tau}\boldsymbol{\cdot}\boldsymbol{n})_{i} \, {\rm d}S\approx{-}\sum_{l=1}^{N_l} c_l V_l f^{l}_{i} + \dfrac{{\rm d}}{{\rm d}t} \int_V u_i \,{\rm d}V + \int_Vk_i\,{\rm d}V, \quad i= \varphi,r,z , \end{equation}
where 
$-\sum _{l=1}^{N_l}c_lV_lf^{l}_{i}$ is the force integrated over the shell of the particle. The ratio between the Lagrangian volume 
$V_l$ and Eulerian volume 
$V_e$ associated with a single Lagrangian marker with index 
$l$ is denoted as 
$c_l=V_l/V_e$. The hydrodynamic load in (A2) has in contrast to Cartesian coordinates (cf. Breugem Reference Breugem2012) an additional term 
$k_i$ stemming from the centrifugal component of (3.2a–c). The term 
$k_i$ is derived by integrating (3.2a–c) over the volume of the particle and observing that from (3.2a–c) the terms 
$u_ru_\varphi /r$ and 
$- u_r^2/r$ are non-vanishing. Here, 
$k_i\equiv (k_\varphi, k_r,k_z)^T$ with 
$k_\varphi = u_ru_\varphi /r$, 
$k_r =-u_\varphi ^2/r$ and 
$k_z =0$. A similar argument for the torque results in
\begin{equation} \int_{\partial V} \boldsymbol{r} \times (\boldsymbol{\tau}\boldsymbol{\cdot}\boldsymbol{n}) \,{\rm d}S = \dfrac{{\rm d}}{{\rm d}t}\int_V \boldsymbol{r} \times\boldsymbol{u} \,{\rm d}V - \int \boldsymbol{r} \times \boldsymbol{f} \,{\rm d}V + \int \boldsymbol{r} \times \boldsymbol{k} \,{\rm d}V. \end{equation}The orientation of the particle, described in a Cartesian frame, follows
\begin{equation} \boldsymbol{\mathsf{I}}_p\dfrac{{\rm d}\boldsymbol{\omega}^b_p}{{\rm d}t} + \boldsymbol{\omega}^b_p \times (\boldsymbol{\mathsf{I}}_p\boldsymbol{\cdot}\boldsymbol{\omega}^b_p) = \boldsymbol{\mathsf{R}}\cdot \int \boldsymbol{r}\times ( \boldsymbol{\tau}\boldsymbol{\cdot}\boldsymbol{n}) \,{\rm d}S + \boldsymbol{\mathsf{R}}\boldsymbol{\cdot} \boldsymbol{T} \end{equation}
and is solved in a local body frame (Allen & Tildesley Reference Allen and Tildesley1989) by employing a quaternion description of the orientation. Here 
$\boldsymbol{\mathsf{R}}$ represents the rotation matrix that converts the torque to the local body frame and is obtained via
\begin{equation} \boldsymbol{\mathsf{R}} =\begin{pmatrix} q_0^2+q_1^2-q_2^2-q_3^2 & 2(q_1q_2+q_0q_3) & 2(q_1q_3-q_0q_2) \\ 2(q_1q_2-q_0q_3) & q_0^2-q_1^2+q_2^2-q_3^2 & 2(q_2q_3-q_0q_2) \\ 2(q_1q_3+q_0q_2) & 2(q_2q_3-q_0q_1) & q_0^2-q_1^2-q_2^2+q_3^2 \end{pmatrix}, \end{equation}
with 
$q_i$ the four quaternion components describing the orientation of the particle. Then, after (A4) is solved for, the quaternions are updated for each particle with
\begin{equation} \dfrac{{\rm d}}{{\rm d} t} \begin{pmatrix} q_0 \\ q_1 \\ q_2 \\ q_3 \end{pmatrix} = \frac{1}{2}\begin{pmatrix} q_0 & -q_1 & -q_2 & -q_3 \\ q_1 & q_0 & -q_3 & q_2 \\ q_2 & q_3 & q_0 & -q_1 \\ q_3 & -q_2 & q_1 & q_0 \end{pmatrix} \begin{pmatrix} 0 \\ \omega_x^b\\ \omega_y^b\\ \omega_z^b \end{pmatrix}. \end{equation}The IBM, here, uses a moving least squares (MLS) approach to perform the interpolation and spreading of forces (Vanella & Balaras Reference Vanella and Balaras2009; Spandan et al. Reference Spandan, Meschini, Ostilla-Mónico, Lohse, Querzoli, de Tullio and Verzicco2017). The MLS approach is beneficial because the transfer functions formulated using MLS conserve momentum on both uniform and stretched grids. At the same time, MLS also retains reasonable accuracy for conserving the exchange of torque between the Eulerian and Lagrangian mesh on slightly stretched grids (Vanella & Balaras Reference Vanella and Balaras2009; de Tullio & Pascazio Reference de Tullio and Pascazio2016), which we employ in our study.
A.2. Short-range collisions: forces and torques
The collision force and torque, 
$\boldsymbol {F}=\{F_\varphi,\ F_r,\ F_z\}$ and 
$\boldsymbol {T}$, respectively, are obtained via a lubrication and soft sphere model. A brief overview of the contact model is presented here and for more details we refer the reader to Costa et al. (Reference Costa, Boersma, Westerweel and Breugem2015) and Ardekani et al. (Reference Ardekani, Costa, Breugem and Brandt2016). The lubrication force 
$F_l$ becomes active when the IBM fails to resolve the lubrication interaction (see also Appendix A.4). In that case the collision force is corrected by 
$\Delta F_l = -6{\rm \pi} \mu R_p u_{ij,n}(\lambda (\varepsilon ) -(\lambda (\varepsilon _{\Delta x}))$, with 
$\varepsilon \equiv \delta _{ij,n}/R_p$ the normalised distance to the radius, 
$\lambda$ the Stokes amplification factor, 
$\varepsilon _{\Delta x}$ the distance from which the model is applied and 
$\delta _{ij,n}$ the gap distance between particles 
$i$ and 
$j$. Parameter 
$R_p$ is locally calculated for the ellipsoid via a Gaussian radius of curvature yielding the best fitting sphere tangent to the surface. Here, 
$\varepsilon _{\Delta _x} =0.075$ for particle–wall collisions and 
$\varepsilon _{\Delta _x} =0.025$ for particle–particle collisions. Following Costa et al. (Reference Costa, Boersma, Westerweel and Breugem2015) roughness is taken into account for the lubrication model for distances 
$0\leq \varepsilon _\sigma <0.001$. The lubrication correction in this range reads 
$\Delta F_l = -6{\rm \pi} \mu R_p u_{ij,n}(\lambda (\varepsilon _\sigma ) -(\lambda (\varepsilon _{\Delta x}))$.
Once particles overlap, the contact model takes over. The contact model considers normal 
$\boldsymbol {F}_{ij,n}$ and tangential 
$\boldsymbol {F}_{ij,t}$ forces, defined as
\begin{equation} \boldsymbol{F}_{ij,n} ={-}k_n \boldsymbol{\delta}_{ij,n} -\eta_n {\boldsymbol{u}}_{ij,n}\quad \text{and} \quad \boldsymbol{F}_{ij,t} = \min(\| -k_t \boldsymbol{\delta}_{ij,t} -\eta_n \boldsymbol{u}_{ij,t}\|,\| -\mu \boldsymbol{F}_{ij,n} \|)\boldsymbol{t}_{ij}, \end{equation}
with 
$\boldsymbol {t}_{ij} = - (k_t\boldsymbol {\delta }_{ij,t}+\eta _t\boldsymbol {u}_{ij,t})/\| k_t\boldsymbol {\delta }_{ij,t}+\eta _t\boldsymbol {u}_{ij,t}\|$, and the coefficients
\begin{equation} k_{n} = \dfrac{m_e ({\rm \pi}^2 + \ln^2 e_{n,d})}{(N\Delta t)^2}\quad \text{and}\quad \eta_n ={-}\dfrac{2m_e \ln e_{n,d}}{N\Delta t}, \end{equation}
with 
$m_e=(m_i^{-1}+m_j^{-1})^{-1}$, and 
$m_i$ and 
$m_j$ the mass of particle 
$i$ and 
$j$, respectively. The coefficients for the tangential force read
\begin{equation} k_t=\dfrac{m_{e,t} ({\rm \pi}^2 + \ln^2 e_{t,d})}{(N\Delta t)^2}\quad \text{and}\quad \eta_t ={-}\dfrac{2m_{e,t} \ln e_{t,d}}{N\Delta t}, \end{equation}
with 
$m_{e,t}=2m_e/7$. The contact time scale is set to 
$N=8$ time steps of the Navier–Stokes solver. We set the coefficients 
$e_{n,d}=0.97$, 
$e_{t,d}=0.39$ and 
$\mu =0.15$ (cf. Costa et al. Reference Costa, Boersma, Westerweel and Breugem2015). The collision force is computed via substepping with an incremental time step of 
$\Delta t/50$. For each substep, an iterative scheme is used that converges once the criterion 
$|x_i^k - x_i^{k-1}|<\Delta z/100$ is met, with 
$x_i^k$ the rotated major axis of the ellipsoid in the inertial frame. Note that obtaining the shortest distance between an ellipsoid and a cylinder is a non-trivial problem. In Appendix A.3 we extend the iterative scheme from Lin & Han (Reference Lin and Han2002) to solve this problem. In Appendix A.4 we validate the collision model with data from available literature.
A.3. Collision of ellipsoids with the inner and outer cylinder
An important aspect of particle–particle or particle–wall collisions is the shortest distance between the two geometries. In the case of non-spherical ellipsoid–ellipsoid interaction, a successful approach has been demonstrated in Ardekani et al. (Reference Ardekani, Costa, Breugem and Brandt2016) by employing the iterative scheme from Lin & Han (Reference Lin and Han2002) to compute the shortest distance. Here, it is shown how by introducing a ghost particle the same algorithm can be used to obtain the shortest distance between an ellipsoid and a cylinder.
Given an ellipsoid 
$E$ for which we would like to calculate the shortest distance to the inner cylinder, we introduce ghost particle 
$E^\prime$. The shortest distance is obtained by positioning 
$E^\prime$ such that the line of shortest distance connecting 
$E$ and 
$E^\prime$: (i) passes through the centreline of the inner cylinder and (ii) is perpendicular to the revolution axis of the inner cylinder. This is achieved as follows. Let 
$\varphi _p$ denote the azimuthal position of 
$E$. Then, 
$E^\prime$ is located at 
$\varphi _p+{\rm \pi}$. Let 
$q =q_0 + q_1 i + q_2 j +q_3 k$ denote the quaternion of 
$E$. Then 
$E^\prime$ should be rotated such that 
$q^\prime = q_0 - q_1 i - q_2 j +q_3 k$. In figure 10 an overview is given for two different configurations, which illustrates the definitions for the distance computations.
Overview of a particle and its ghost particle for two different configurations (not to scale). Here 
$\delta _{\min }$ and 
$\delta _{\max }$ are used to calculate the shortest and longest distance between particle 
$E$ to the inner and outer cylinder, respectively. (a) Top view and (b) side view.
The shortest distance 
$\delta _{\min }$ between 
$E$ and 
$E'$ can now be found in the conventional way as described in Lin & Han (Reference Lin and Han2002). The shortest distance from 
$E$ to the inner cylinder is then obtained via 
$\frac {1}{2}(\delta _{\min }-2r_i)$.
The distance to the outer cylinder can be obtained by finding the largest distance 
$\delta _{\max }$ between 
$E$ and 
$E'$. One obtains 
$\delta _{\max }$ by picking the second root from the second-order polynomial in the algorithm. Effectively, 
$\delta _{\max }$ is obtained by interchanging the 
$\min$ and 
$\max$ functions when computing the step size (cf. Lin & Han Reference Lin and Han2002, (2.2)), i.e.
\begin{equation} \left.\begin{array}{c@{}} t_1 = \min \{ t \in [0,1]: (1 -t)c_1 + t c_2 \in E\},\\ t_2 = \max \{ t \in [0,1]: (1-t)c_1 + t c_2 \in E^\prime \}. \end{array}\right\} \end{equation}
The particle distance to the outer cylinder is obtained via 
$\frac {1}{2}(\delta _{\max }-2r_o)$.
A.4. Collision validations
Here, we validate our code similarly to that performed in Costa et al. (Reference Costa, Boersma, Westerweel and Breugem2015). The validation consists of two parts. First, we consider the force exerted on a single particle approaching the inner cylinder and compare it with a theoretical prediction developed for a sphere approaching a flat plane at creeping flow conditions (Brenner Reference Brenner1961; Jeffrey Reference Jeffrey1982). The sphere has a diameter 
$D$ such that 
$D/x=16$ (uniform grid), and is positioned in a TC configuration with particle diameter to gap width ratio 
$d/D=12$. The inner cylinder to particle radius complies with 
$r_i/r_p=200$ and the particle Reynolds number is 
$\textit {Re}_p=0.1$. We position the sphere in a quiescent flow (cylinders are fixed) at a fixed distance and impose a velocity on its surface in the radial direction (Reynolds number 
$\textit {Re}_p=0.1$). We then let the simulation run until a steady state is reached and report the final value of the force 
$f_r$ acting on the sphere.
In figure 11(a) we report the shortest distance 
$h$ versus normalised 
$f_r$. The results we observe are very similar to the results found in Costa et al. (Reference Costa, Boersma, Westerweel and Breugem2015), namely for distances of 
$h/r_p\gtrsim 0.03$ the numerical values of 
$f_r$ are in close agreement with the theoretical predictions. This figure shows that the simulation is not able to properly capture the asymptotic behaviour for 
$h/r_p<0.1$ due to the insufficient grid resolution between the sphere and cylinder.
(a) Normalised force exerted on a sphere approaching a cylinder versus the shortest distance between the geometries. The radius ratio is set to 
$r_i/r_p=200$. The theoretical prediction is from Brenner (Reference Brenner1961). (b) A spherical particle falling on a flat plate at 
$\textit {Ga}=130.4$ and density ratio 
$\rho _p/\rho _f=8.34$. The results are compared with experimental data reported in Gondret et al. (Reference Gondret, Lance and Petit2002) and numerical results from Costa et al. (Reference Costa, Boersma, Westerweel and Breugem2015).
The second validation is the replication of a sphere impacting on a flat wall, which bounces up and down until it is at rest. The non-dimensional parameters of the particle in relation to the fluid describing the problem are the Galileo number 
$\textit {Ga}=130.4$ and density ratio 
$\rho _p/\rho _f=8.34$, with 
$\textit {Ga}=U_gD/\nu$. The characteristic gravitational velocity is defined as 
$U_g=\sqrt {|\rho _p/\rho _f-1|gD}$. The computational domain is 
$30D$ in the vertical direction and 
$10D$ in the horizontal direction. The top and bottom boundaries are walls with zero velocity imposed on them. The side boundaries of the domain are periodic. A single sphere (uniform grid 
$D/\Delta x =16$) is then moved at the expected terminal velocity (
$u_z/U_g=1.25$) for four particle diameters starting from the initial position of 
$29D$ above the bottom surface. From then on, the particle is allowed to freely fall. For a close comparison, we set the collision parameters equal to those reported in table II of Costa et al. (Reference Costa, Boersma, Westerweel and Breugem2015). Our findings of the particle velocity, 
$u_z$, are in close agreement with those reported in Gondret, Lance & Petit (Reference Gondret, Lance and Petit2002) and Costa et al. (Reference Costa, Boersma, Westerweel and Breugem2015) (see figure 11b).
Appendix B. Curvature effects
The effect of the TC curvature on the particle motion is examined here under conditions where Jeffery's equations are approximately valid. This is performed by comparing the motion of a single particle in plane Couette flow with the motion of a particle in TC flow. The TC configuration is taken as in the main study with a radius ratio of 
$\eta =5/7$. The inner and outer Reynolds numbers of the cylinders are 
$Re_i=-0.48$ and 
$Re_o=0.48$, respectively (
$\textit {Ta}=1.0$). This flow is characterised as part of the circular Couette flow regime. Curvature effects are assessed by varying the particle size. The particle is fixed at the location where the flow satisfies 
$u_\varphi = 0$. In comparison, the plane Couette geometry is subject to the same conditions with 
$Re_l=0.48$ and 
$Re_u=-0.48$ for the lower and upper wall, respectively. The particle is pinned at mid-gap.
The particle dynamics is presented in figure 12. The velocity gradient at the particle centre is equal for both the TC system and plane Couette set-up. Two observations can be made. First, the particle in TC flow is observed to rotate at a frequency 
$\omega _b$ that is 
$22\,\%$ slower compared with plane Couette flow. Second, the difference in rotation rate differs more between particle sizes for TC flow (difference of 
$3.5\,\%$) in comparison with plane Couette flow (<1 %). The greater difference between different particle sizes for TC flow is understood as a curvature effect.
Curvature effects related to particle size. The rotational velocity 
$\omega _b$ of a single pinned (where 
$u_\varphi =0$) particle is tracked. The quantity is normalised with 
$\omega _d=u_w/d$, with 
$u_w$ the wall velocity and 
$d$ the gap width. (a) Plane Couette flow and (b) TC flow. For comparison, the result derived by Jeffery is included (Jeffery Reference Jeffery1922).
Open access
