1. Introduction
Taylor–Couette flow (TCF) between two coaxial cylinders has been one of the paradigmatic flows for a long time, since it allows the study of various physical phenomena in fluid dynamics, such as flow instability and transition (Coles Reference Coles1965; DiPrima, Eagles & Ng Reference Diprima, Eagles and Ng1984; Marques & Lopez Reference Marques and Lopez1997; Caton, Janiaud & Hopfinger Reference Caton, Janiaud and Hopfinger1999; Hristova et al. Reference Hristova, Roch, Schmid and Tuckerman2002), nonlinear dynamics (Pfister & Rehberg Reference Pfister and Rehberg1981; Mullin, Cliffe & Pfister Reference Mullin, Cliffe and Pfister1987), pattern formation (Andereck, Liu & Swinney Reference Andereck, Liu and Swinney1986; Boubnov, Gledzer & Hopfinger Reference Boubnov, Gledzer and Hopfinger1995; Lim, Chew & Xiao Reference Lim, Chew and Xiao1998; Guillerm et al. Reference Guillerm, Kang, Savaro, Lepiller, Prigent, Yang and Mutabazi2015) and turbulence (Lathrop, Fineberg & Swinney Reference Lathrop, Fineberg and Swinney1992a,Reference Lathrop, Fineberg and Swinneyb; Huisman et al. Reference Huisman, Van Gils, Grossmann, Sun and Lohse2012; Grossmann, Lohse & Sun Reference Grossmann, Lohse and Sun2016). Moreover, flows of a dispersed phase (bubbles, drops or particles) and the migration behaviour of dispersions have been widely examined in the Taylor–Couette system due to numerous practical applications in chemical engineering, such as the oil industry, filtration and transformation of thermal energy (Chan & Leal Reference Chan and Leal1981; Wereley & Lueptow Reference Wereley and Lueptow1999; Resende et al. Reference Resende, Tardioli, Fernandez, Ferreira, Giordano and Giordano2001; Rudman Reference Rudman2004; Climent, Simonnet & Magnaudet Reference Climent, Simonnet and Magnaudet2007). In this work, we explore a semidilute suspension in a Taylor–Couette flow using rheological laws such as the suspension balance model (SBM): (1) to understand if such a model could be used as a tool to examine suspensions in a TCF, (2) to characterize the particle migration and (3) to examine parameters that might be challenging to be analysed using experiments including friction and torque coefficients. Although several experiments were performed to examine the behaviour of particles in a TCF (Ali et al. Reference Ali, Mitra, Schwille and Lueptow2002; Majji, Banerjee & Morris Reference Majji, Banerjee and Morris2018; Majji & Morris Reference Majji and Morris2018), to the best of the authors’ knowledge, little work has been performed to predict the flow transitions in a TCF using a continuum model.
At low Reynolds numbers, neutrally buoyant particles in dilute suspensions follow the streamlines. However, micron-size particles undergo a self-diffusion phenomenon arising from hydrodynamic interactions between particles at higher concentrations. It triggers the shear-induced particle migration in non-uniform shear flows (Eckstein, Bailey & Shapiro Reference Eckstein, Bailey and Shapiro1977; Leighton & Acrivos Reference Leighton and Acrivos1987; Phillips, Armstrong & Brown Reference Phillips, Armstrong and Brown1992; Sierou & Brady Reference Sierou and Brady2004). To understand the shear-induced diffusion, a shear flow of particulate suspensions in the Taylor–Couette geometry has been utilized in various studies (Phillips, Armstrong & Brown Reference Phillips, Armstrong and Brown1992; Tetlow et al. Reference Tetlow, Graham, Ingber, Subia, Mondy and Altobelli1998; Morris & Boulay Reference Morris and Boulay1999; Fang et al. Reference Fang, Mammoli, Brady, Ingber, Mondy and Graham2002). At low Reynolds regime, particles migrate from regions of a higher shear rate to those of a lower shear rate. As a result, a non-uniform distribution of the particle concentration is caused and the velocity profile is altered.
In Taylor–Couette flows, the flow becomes unstable and a transition occurs at sufficient ratio rates of the inner and outer cylinders which induce the centrifugal instability. For the flow of a pure Newtonian fluid, the flow instability and transitions have been well-documented by numerous studies as mentioned earlier. For example, Andereck, Liu & Swinney (Reference Andereck, Liu and Swinney1986) determined transitions as functions of Reynolds numbers based on the rotating angular velocity of the inner and outer cylinders, and mapped out the observed flow patterns. By contrast, a few studies have examined the instability and transition of finite-size particles contained in the Taylor–Couette system with the rotating inner cylinder and stationary outer one. Ali et al. (Reference Ali, Mitra, Schwille and Lueptow2002) carried out a linear stability analysis (LSA) for dilute suspensions 
$({\phi _b} \le 0.05)$ of rigid spherical particles in cylindrical Couette flow, where 
${\phi _b}$ is the bulk particle volume fraction. Their result of LSA showed that dispersed particles destabilize the flow and the degree of destabilization depends on the density ratio between the particle and the suspending fluid. However, their experiments with neutrally buoyant particles of 
${\phi _b} \le 0.005$ present a stabilizing effect of particles by increasing the effective viscosity. Dherbécourt et al. (Reference Dherbécourt, Charton, Lamadie, Cazin and Climent2016) and Rida et al. (Reference Rida, Cazin, Lamadie, Dherbécourt, Charton and Climent2019) performed a series of experiments for neutrally buoyant particles of size dp = 800–1500 μm (
$d_p$ is the diameter of particles), gap-to-particle ratio of 7.3–13.75 and volume fractions ranging from 0.01 to 0.08, and observed an improved mixing in TVF. Majji & Morris (Reference Majji and Morris2018) conducted experiments for the flow of a dilute suspension 
$({\phi _b} = 0.001)$ confined in a finite cylindrical annulus. Similar to the pure Newtonian fluid flow, they observed that very dilute suspension undergoes transitions from the circular Couette flow (CCF) via Taylor vortex flow (TVF) to wavy vortex flow (WVF). They also showed that, in the CCF, particles migrate to an equilibrium location near the middle of the annulus with an offset toward the inner cylinder because of the competition between the shear gradient in the flow and the wall interactions. For the TVF state, however, particles were trapped in a circular equilibrium region inside each vortex. Unlike in the CCF and TVF, particles in the WVF did not reach a steady equilibrium location and they were uniformly distributed in the annular region.
In another study, Majji, Banerjee & Morris (Reference Majji, Banerjee and Morris2018) performed a set of flow-visualization experiments on inertial transitions of neutrally buoyant suspensions for various Reynolds numbers (Re) and particle volume fractions 
$(0 \le {\phi _b} \le 0.3)$ to investigate the influence of particle loading and size on flow transitions. For 
$0.05 \le {\phi _b} \le 0.15$ and 
$\mathrm{\epsilon }\;( = d\textrm{/}a) = 60$, where d is the gap width between two cylinders and a is the radius of particles, when the Reynolds number of the inner cylinder was slowly decreased (the decreasing-Re protocol), they discovered additional non-axisymmetric flow states, namely spiral vortex flow (SVF) and ribbons (RIB). For the reduced particle size (i.e. 
$\mathrm{\epsilon } = 200$) at 
${\phi _b} = 0.1$, only the RIB was observed between TVF and CCF. While these are absent in a similar TCF with a pure Newtonian fluid, both spiral vortices and ribbons are known to appear in counter-rotating TCF (Coles Reference Coles1965; Andereck, Liu & Swinney Reference Andereck, Liu and Swinney1986), as primary bifurcations from the CCF in the narrow-gap limit. Their experiments with 
${\phi _b} = 0.1$ also showed a hysteresis during the transitions of WVF ↔ TVF and TVF ↔ SVF. By further increasing the particle volume fractions (i.e. 
${\phi _b} = 0.2,\;0.3$) at 
$\mathrm{\epsilon } = 60$, another non-axisymmetric flow pattern, the so-called WSVF, was observed. In addition, by decreasing Re, the suspensions transitioned following the sequence of WSVF → WVF → WSVF → TVF → SVF → RIB → CCF for 
${\phi _b} = 0.2$, whereas, for 
${\phi _b} = 0.3$, only non-axisymmetric structures appeared and the sequence was simplified as WSVF → SVF → CCF. Recently, Gillissen & Wilson (Reference Gillissen and Wilson2019) employed the LSA of axisymmetric perturbations of the two-fluid theory to examine the stability analysis for the Taylor–Couette flow of suspensions from dilute to highly concentrated suspensions 
$(0 \le {\phi _b} \le 0.5)$. They predicted the critical, effective Taylor number based on the effective suspension viscosity as a function of particle volume fraction 
$({\phi _b})$ and revealed the destabilization effect of the flow due to non-Brownian spherical particles. Recently, Ramesh, Bharadwaj & Alam (Reference Ramesh, Bharadwaj and Alam2019) visualized flows and measured the velocity field using particle image velocimetry (PIV) for the flow of non-Brownian suspension in a narrow-gap Taylor–Couette cell. They found that the bifurcation sequence for dilute suspensions 
$({\phi _b} \lt 0.05)$ remained the same as in the flow of a pure Newtonian fluid (i.e. CCF → TVF → WVF with increasing Re). However, for 
${\phi _b} \ge 0.05$, the non-axisymmetric SVF appeared as the primary bifurcation when Re was slowly reduced and the transition route of WVF → TVF → SVF → CCF was presented by decreasing Re. On the other hand, new coexisting states, namely TVF + WVF and TVF + SVF, were uncovered when the rotation rate of the inner cylinder was gradually ramped up. The TVF + WVF occurred with the onset of the primary instability from the CCF for 
$0.05 \le {\phi _b} \le 0.12$, whereas the TVF + SVF appeared between CCF and TVF as the secondary bifurcation for 
${\phi _b} \ge 0.13$. They also found that the secondary bifurcation from TVF to WVF is subcritical (or hysteretic) for 
${\phi _b} \ge 0.05$ and the value of the critical Reynolds number is reduced for both primary and secondary transitions. Lately, Ramesh & Alam (Reference Ramesh and Alam2020) presented a new pattern for 
${\phi _b} \ge 0.1$, the so-called interpenetrating spiral vortices (ISVs) which is a coexisting state of upward and downward propagating spiral (or helical) vortices, in the suspension Taylor–Couette flow of a higher radius ratio (η) and aspect ratio (
$\varGamma$).
The focus of this study is to conduct, for the first time, numerical simulations for non-Brownian and non-colloidal suspension flows employing rheological constitutive models known as the SBM. In particular, we assume that the suspension undergoes the shear-induced migration while the inertial migration of particles is neglected. To compare our results with the experimental data reported by Majji, Banerjee & Morris (Reference Majji, Banerjee and Morris2018), herein, we consider a suspension flow of 
${\phi _b}\; = 0.1$ with various particle sizes of 
$\mathrm{\epsilon } = d\textrm{/}a = 60,\;200$. The primary and secondary bifurcations are investigated with varying Reynolds number based on the rotating angular velocity of the inner cylinder, and the transition states are presented by determining the flow patterns. We characterize the flow structures, particle concentration fields and the dimensionless coefficients related to the torque exerted on the inner cylinder. This paper is organized as follows. In § 2, mathematical formulations for the flow of suspension are described following with the rheological constitutive laws, numerical procedure and parameters. The results and discussion are presented in § 3. Section 4 addresses the summary and conclusion of this study.
2. Problem formulation
A flow of neutrally buoyant, non-colloidal, rigid, spherical particles suspended in a viscous fluid of viscosity 
$\mu $ and density 
$\rho $ is considered. The suspension is confined between two coaxial cylinders of the gap width d with the inner one rotating at a constant angular velocity 
$({\varOmega _i})$, while the outer one is fixed (figure 1). The radii of the inner and outer cylinders are 
${r_i}$ and 
${r_o}( = {r_i} + d)$, respectively. We assume that the suspensions are in the Stokes flow regime, i.e. 
$R{e_p} \ll 1$, where 
$R{e_p} = \rho \dot{\gamma }{a^2}\textrm{/}\mu $ is the particle Reynolds number based on the fluid shear rate 
$\dot{\gamma }$ and the radius of the suspended particles a. We also consider the suspensions in the limit of 
$Pe \to \infty$ where 
$Pe = 6{\rm \pi} \mu \dot{\gamma }{a^3}\textrm{/}kT$ is the Péclet number defined with the thermal energy kT.
Figure 1. Schematic diagram of a flow of suspensions in Taylor–Couette geometry.
2.1. Balance equations and constitutive laws
To model the flow of suspension, we employ the SBM introduced by Nott & Brady (Reference Nott and Brady1994) and later by Morris & Brady (Reference Morris and Brady1998), in which mass and momentum balances are written for the bulk suspension and they are solved simultaneously for the particle concentration and bulk velocity.
Using the SBM for non-Brownian particles suspended in an incompressible fluid, the conservation equations for particle mass and momentum are obtained simply by ensemble and volume averages over the particles as (Drew & Lahey Reference Drew, Lahey and Roco1993)
\begin{gather}\; \frac{{\partial \phi }}{{\partial t}} + \boldsymbol{\nabla }\boldsymbol{\cdot }\phi {\langle \boldsymbol{u}\rangle _p} = 0,\end{gather}
\begin{gather}{\rho _p}\phi \frac{{{\textrm{D}_p}{{\langle \boldsymbol{u}\rangle }_p}}}{{\textrm{D}t}} = {\langle \boldsymbol{b}\rangle _p} + {\langle \boldsymbol{F}\rangle _p} + \boldsymbol{\nabla }\boldsymbol{\cdot }{\langle \boldsymbol{\varSigma }\rangle _p},\end{gather}
where the substantial derivative, 
${\textrm{D}_p}/\textrm{D}t = \partial /\partial t\; + {\langle \boldsymbol{u}\rangle _p}\boldsymbol{\cdot }\boldsymbol{\nabla }$, follows the average particle motion (Nott & Brady Reference Nott and Brady1994; Morris & Brady Reference Morris and Brady1998). It has been assumed that the density of a particle 
${\rho _p}$ is constant. 
${\langle \;\rangle _p}$ represents a particle-phase average and thus 
${\langle \boldsymbol{u}\rangle _p}$ is the average particle velocity. In (2.2), 
${\langle \boldsymbol{b}\rangle _p}$, 
${\langle \boldsymbol{F}\rangle _p}$ and 
${\langle \boldsymbol{\varSigma }\rangle _p}$ are the average body force, hydrodynamic force and stress of the particles, respectively (Nott & Brady Reference Nott and Brady1994; Morris & Brady Reference Morris and Brady1998).
The balance equations for the suspension are derived from the mass and momentum conservation equations by taking the averages over the bulk suspensions (Drew & Lahey Reference Drew, Lahey and Roco1993), and can be stated as
\begin{gather}\boldsymbol{\nabla }\boldsymbol{\cdot }\langle \boldsymbol{u}\rangle = 0,\end{gather}
\begin{gather}\frac{{\textrm{D}\langle \rho \boldsymbol{u}\rangle }}{{\textrm{D}t}} = \langle \boldsymbol{b}\rangle + \boldsymbol{\nabla }\boldsymbol{\cdot }\langle \boldsymbol{\varSigma }\rangle ,\end{gather}
where the material derivative is that following the suspension-averaged, 
$\textrm{D}/\textrm{D}t = \partial /\partial t\; + \langle \boldsymbol{u}\rangle \boldsymbol{\cdot }\boldsymbol{\nabla }$ (Nott & Brady Reference Nott and Brady1994; Morris & Brady Reference Morris and Brady1998). 
$\langle \;\rangle $ denotes a suspension average and then 
$\langle \boldsymbol{u}\rangle $ is the suspension-averaged velocity. Here, 
$\langle \boldsymbol{b}\rangle $ and 
$\langle \boldsymbol{\varSigma }\rangle $ are the average body force acting on the suspension and the average suspension stress, respectively.
For neutrally buoyant particles, the particle mass and momentum equations ((2.1) and (2.2)) can then be coupled by vanishing Reynolds number and taking infinite Péclet number (Morris & Brady Reference Morris and Brady1998; Morris & Boulay Reference Morris and Boulay1999). Therefore, the particle mass conservation (2.1) reduces to
\begin{equation}\frac{{\partial \phi }}{{\partial t}} + \langle \boldsymbol{u}\rangle \boldsymbol{\cdot }\boldsymbol{\nabla }\phi ={-} \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{j},\quad \boldsymbol{j} = \frac{{2{a^2}}}{{9\mu }}f(\phi )\boldsymbol{\nabla }\boldsymbol{\cdot }{\langle \boldsymbol{\varSigma }\rangle _p}.\end{equation}
Here, 
$\boldsymbol{j}$ implies the particle migration flux relative to the mean motion of the suspension; 
$f(\phi )$ is the sedimentation hindrance function that indicates the sedimentation rate of a homogeneous suspension of spheres at a volume fraction (ϕ) to the isolated Stoke settling velocity. We take the form reported by Miller & Morris (Reference Miller and Morris2006) for the sedimentation hinderance function as 
$f(\phi ) = (1 - \phi /{\phi _m}){(1 - \phi )^{\alpha - 1}}$ with 
$\alpha = 4$. This formula is similar to the one proposed by Richardson & Zaki (Reference Richardson and Zaki1954) as 
$f(\phi ) = {(1 - \phi )^\alpha }$, but it has been altered to ensure that particle migration ceases when the volume fraction approaches to the maximum packing 
${\phi _m}$ (Miller & Morris Reference Miller and Morris2006). We also ignore the impact of the particle lift force causing the inertial migration in the flow.
For non-Brownian suspensions, the constitutive law for the suspension stress stated in (2.4) can be given by (Nott & Brady Reference Nott and Brady1994; Morris & Brady Reference Morris and Brady1998)
\begin{equation}\langle \boldsymbol{\varSigma }\rangle ={-} {\langle p\rangle _f}\boldsymbol{\mathsf{I}} + 2\mu \langle \boldsymbol{S}\rangle + {\langle \boldsymbol{\varSigma }\rangle _p},\end{equation}
where 
${\langle p\rangle _f}$ is the average pressure in the fluid, I is the identity tensor and 
$\langle \boldsymbol{S}\rangle $ is the bulk suspension rate of strain defined as 
$\boldsymbol{S} = [\boldsymbol{\nabla }\boldsymbol{u}\; + \; {(\boldsymbol{\nabla }\boldsymbol{u})^\textrm{T}}]/2$. The particle stress 
${\langle \boldsymbol{\varSigma }\rangle _\textrm{p}}$ is defined as (Morris & Boulay Reference Morris and Boulay1999)
\begin{equation}{\langle \boldsymbol{\varSigma }\rangle _p} ={-} \mu {\mu _n}(\phi )\dot{\gamma }\boldsymbol{\mathsf{Q}} + \; 2\mu {\mu _p}(\phi )\langle \boldsymbol{S}\rangle ,\end{equation}
where 
$\dot{\gamma } = {(2\boldsymbol{S}:\boldsymbol{S})^{1/2}}$ is the local shear rate. The first term in (2.7) represents the particle contribution to the normal stress. The function 
${\mu _n}(\phi )$ is the ‘normal stress viscosity’ that can be expressed as 
${\mu _n}(\phi ) = {K_n}{(\phi /{\phi _m})^2}{(1 - \phi /{\phi _m})^{ - 2}}$ where 
${K_n} = 0.75$ (Morris & Boulay Reference Morris and Boulay1999; Miller & Morris Reference Miller and Morris2006); 
$\mu_p$ is the particle contribution to the shear viscosity. The constant tensor Q describes the anisotropy of the normal stresses as
\begin{equation}\boldsymbol{\mathsf{Q}} = \left( \begin{array}{*{20}{c}} \lambda_{\textrm{1}}&\textrm{0}&\textrm{0}\\ \textrm{0}&\lambda_{\textrm{2}}&\textrm{0}\\ \textrm{0}&\textrm{0}&\lambda_{\textrm{3}} \end{array} \right).\end{equation}
The values of 
${\lambda _1} = 0.8$, 
${\lambda _2} = 1$ and 
${\lambda _3} = 0.5$ have provided a good fit to the experimental data in a wide-gap Couette flow of concentrated suspensions (Morris & Boulay Reference Morris and Boulay1999), while a set of 
${\lambda _1} = 0.6$, 
${\lambda _2} = 1.2$ and 
${\lambda _3} = 1.2$ has been also adopted to the flow of highly concentrated suspensions in the same geometry (Fang et al. Reference Fang, Mammoli, Brady, Ingber, Mondy and Graham2002). In this study, in a narrow-gap Couette flow, we choose 
${\lambda _1} = {\lambda _2} = 1$ and 
${\lambda _3} = 0.8$ using our numerical tests. The last term of (2.7) is the particle contribution to the shear stress, and it can be combined with the shear stress of the fluid 
$2{\mu}\langle \boldsymbol{S}\rangle $. Finally, the suspension stress can be written as
\begin{equation}\langle \boldsymbol{\varSigma }\rangle ={-} {\langle p\rangle _f}\boldsymbol{\mathsf{I}} - {\mu}{\mu _n}(\phi )\dot{\gamma }\boldsymbol{\mathsf{Q}} + 2\mu {\mu _s}(\phi )\langle \boldsymbol{S}\rangle .\end{equation}
Here, 
${\mu _s}(\phi )$ is the effective shear viscosity of the bulk suspension, where 
${\mu _s} = 1 + {\mu _p}$. To define the effective shear viscosity of the suspension, we employ the Krieger's empirical correlation 
${\mu _s}(\phi ) = {(1 - \phi /{\phi _m})^{ - 1.82}}$ where ϕm = 0.68 is the maximum packing particle volume fraction (Krieger Reference Krieger1972). Hereafter, the suspension average 
$\langle \;\rangle $ is omitted for simplicity.
2.2. Control parameters
We have defined the suspension Reynolds number as 
$R{e_s} = \rho {r_i}{\varOmega _i}d/\mu {\mu _s}({\phi _b})$, then used it as a control parameter based on the bulk particle volume fraction 
${\phi _b}$ (Majji, Banerjee & Morris Reference Majji, Banerjee and Morris2018). The radius ratio 
$(\eta = {r_i}/{r_o})$ of concentric cylinders is 
$\eta = 0.877$, and the aspect ratio 
$(\varGamma = {L_z}/d)$ is 
$\varGamma = 4\;({\cong} 2{\lambda _c}/d)$ where 
${\lambda _c}$ is the critical wavelength for the flow of a pure Newtonian fluid. Note that Ali et al. (Reference Ali, Mitra, Schwille and Lueptow2002) and Majji, Banerjee & Morris (Reference Majji, Banerjee and Morris2018) found the same critical wavenumber for suspensions from the LSA as it is for a pure Newtonian fluid. To compare our results with Majji, Banerjee & Morris (Reference Majji, Banerjee and Morris2018), we have considered suspensions of 
${\phi _b} = 0.1$ with the particle size 
$\mathrm{\epsilon } = d/a = 60,\;200$. For non-dimensionalization, we have chosen the gap width d as the scale for length, the rotating velocity of the inner cylinder 
${r_i}{\varOmega _i}$ as the scale for velocity and 
$\rho {d^2}/\mu $ as the time scale. While the suspension Reynolds number has been varied up to 
$R{e_s} = 180$, based on the relation between the particle and suspension Reynolds number as 
$R{e_p} = R{e_s}{\mu _s}/{\mathrm{\epsilon }^2}$, the maximum particle Reynold number was computed as 
$R{e_p} \lt \; {10^{ - 1}}$, i.e. negligible particles inertia. For convenience, we also use a reduced radial coordinate as 
$x = (r - {r_i})/d\; \in [0,1]$.
2.3. Numerical schemes and boundary conditions
The governing equations (2.3)‒(2.5) were discretized using a finite volume method in a cylindrical coordinate system (r, φ, z) shown in figure 1. A second-order central difference scheme was used for spatial discretization of derivatives except for the convective term 
$(\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{\nabla }\phi )$ of the particle mass conservation (2.5) where we employed the QUICK (quadratic upstream interpolation for convective kinematics) scheme for the discretization (Kang & Mirbod Reference Kang and Mirbod2020). A hybrid scheme was utilized for the time advancement: nonlinear terms and cross-diffusion terms were explicitly advanced by a third-order Runge–Kutta scheme, and the diffusion terms were implicitly advanced by the Crank–Nicolson method (Kang et al. Reference Kang, Meyer, Mutabazi and Yoshikawa2017a,Reference Kang, Meyer, Yoshikawa and Mutabazib). A fractional step method was used for time integration, and the Poisson equation that resulted from the second stage of the fractional step method was solved by a fast Fourier transform (FFT) (Kim & Moin Reference Kim and Moin1985).
The no-slip boundary condition was employed at the cylindrical surfaces, and the migration flux of particles was set to be zero at the walls. The flow and particle concentrations were also assumed to be periodic in the axial direction (z). These can be expressed as
\begin{equation}\left. {\begin{array}{c@{}} {\boldsymbol{u} = {r_i}{\Omega_i},\quad \boldsymbol{j}\boldsymbol{\cdot }\boldsymbol{n} = 0\quad at\;\; r = {r_i},}\\ {\boldsymbol{u} = 0,\quad \boldsymbol{j}\boldsymbol{\cdot }\boldsymbol{n} = 0\quad at\;\; r = {r_o},}\\ {\boldsymbol{u}(r,\varphi ,z) = \boldsymbol{u}(r,\; \varphi ,\; z + \varGamma ),}\\ {\phi (r,\varphi ,z) = \phi (r,\varphi ,z + \varGamma ).} \end{array}} \right\}\end{equation} Computations were performed in a grid system with 
$96(r) \times 512(\varphi ) \times 64(z)$ determined by grid refinement study. The radius ratio of 
$\eta = 0.877$ leads to a long domain length in the azimuthal direction. Moreover, for the WVF, flow and particle concentration fields vary with the azimuthal direction. For these reasons, we used large enough grid points in the azimuthal direction to resolve the variations of both velocity and concentration fields. More grid points were allocated near the cylinder walls in the radial direction (r) with 
$\Delta {r_{min }} = \; 0.005d$, where 
$\Delta{r_{min}}$ is the minimum grid size, whereas the grid cells in the azimuthal (φ) and axial (z) directions were uniform.
2.4. Validation
To verify our code, we performed a direct comparison with previous studies for the wide-gap steady Couette flow shown in figure 2 that presents the variations of particle volume fraction and azimuthal velocity for a wide-gap steady Couette flow of 
${\phi _b} = 0.5$ explored by Phillips, Armstrong & Brown (Reference Phillips, Armstrong and Brown1992). Our computed distribution of particle volume fraction, shown in figure 2(a), is consistent with the measurements of Phillips, Armstrong & Brown (Reference Phillips, Armstrong and Brown1992) and also previous numerical predictions that employed the SBM (Morris & Boulay Reference Morris and Boulay1999; Fang et al. Reference Fang, Mammoli, Brady, Ingber, Mondy and Graham2002). The profile of azimuthal velocity also shows a good agreement with that of Morris & Boulay (Reference Morris and Boulay1999) as displayed in figure 2(b).
Figure 2. Comparisons of (a) particle volume fraction (ϕ) and (b) azimuthal velocity (
${u_\varphi }$) normalized by 
${r_i}{\varOmega _i}$ for 
${\phi _b} = \; 0.5$ in a wide-gap Couette flow with η = 0.269.
3. Results
3.1. Base flow
To analyse in detail the suspension flow transition and structure, we performed numerical simulations in our proposed geometry. In the classical Taylor–Couette flow, the transition from CCF to TVF occurs at a primary bifurcation when the Reynolds number (Re) is larger than its critical value 
$(Re \ge R{e_c})$. Below the critical Reynolds number, the slowly rotating inner cylinder causes a simple shear flow where the base state is stationary, an axisymmetric CCF of suspension. In the CCF, the velocity and particle concentration depend only on the radial coordinate, i.e. 
$\boldsymbol{u} = {u_\varphi }(r){\boldsymbol{e}_\varphi }$ and 
$\phi = \phi (r)$, see figure 3. As can be observed, the distributions of particle concentration and azimuthal velocity are linear over the gap, and the velocity profile is almost in-line with that of a pure fluid flow. Therefore, one can expect that the profiles of a semidilute suspension 
$({\phi _b} = 0.1)$ in a narrow gap display dissimilar properties with those of the highly concentrated suspension flow in a wide-gap Couette flow presented in figure 2.
Figure 3. Profiles of the base flow; (a) particle volume fraction (ϕ), (b) azimuthal velocity (
${u_\varphi }$) for 
$R{e_s} = 110$ and for 
${\phi _b} = \; 0.1$ of both particle sizes.
It should be noted that Majji & Morris (Reference Majji and Morris2018) and Baroudi, Majji & Morris (Reference Baroudi, Majji and Morris2020) experimentally observed the inertial migration that particles migrate to an equilibrium location near the middle of the annulus with an offset toward the inner cylinder for very dilute suspensions (ϕb = 0.001 for d/a = 40 and 0.01 for d/a = 60.8); however, none of the previous studies reported inertial migration and the corresponding concentration profiles for semidilute and concentrated suspensions. One expects that for very dilute suspensions the lift force acting on the particles is dominant, leading to particle migration towards the middle between walls as theoretically verified by Cox & Brenner (Reference Cox and Brenner1968), Ho & Leal (Reference Ho and Leal1974), Vasseur & Cox (Reference Vasseur and Cox1976), McLaughlin (Reference Mclaughlin1993), Hogg (Reference Hogg1994) and Asmolov (Reference Asmolov1999). Indeed, the collision of particles with each other and the cylindrical walls might not be negligible for semidilute and concentrated suspensions; however, the related lift force acting on particles due to particle inertia, that needs to be added as an additional term in SBM, has not yet been identified in any literature. Because of these reasons, we assumed that the particle inertia and the corresponding lift force are negligible, while the shear-induced diffusion is dominant in the flow which results in the linear profile of the particle concentration over the annulus gap as shown in figure 3(a).
3.2. Flow instability
As the rotation rate of the inner cylinder and consequently 
$R{e_s}$ increases, the centrifugal force exceeds the viscous force and the flow becomes unstable in the annulus region. To characterize the critical states of the flow, we have employed the Landau model which describes the evolution of the flow perturbation in its weakly nonlinear regime and it can be determined as (Landau & Lifshitz Reference Landau and Lifshitz1976; Guckenheimer & Holmes Reference Guckenheimer and Holmes1983; Kang et al. Reference Kang, Meyer, Mutabazi and Yoshikawa2017a,Reference Kang, Meyer, Yoshikawa and Mutabazib, Reference Kang, Meyer, Yoshikawa and Mutabazi2019)
\begin{equation}\frac{{\textrm{d}A}}{{\textrm{d}t}} = \sigma (1 + i{c_1})A - l(1 + i{c_2})|A{|^2}A + \cdots .\end{equation}
The parameters σ and l represent the linear growth rate of the perturbation and the Landau constant where the sign of l indicates the nature of the bifurcation (i.e. supercritical vs subcritical), and the constants c 1 and c 2 are the linear and nonlinear dispersion coefficients. We have introduced the norm of the radial velocity component at the central surface in the annulus to define the amplitude of the perturbation 
$|A|$, as
\begin{equation}|A |= \frac{1}{{2{\rm \pi} \varGamma }}\int_0^\varGamma {\int_0^{2{\rm \pi} } {|{{u_r}({r_m},\; \varphi ,\; z)} |\,\textrm{d}\varphi \,\textrm{d}z} } ,\quad \textrm{where}\;\; {r_m} = ({r_i} + {r_o})/2.\end{equation}
The growth rates (σ) obtained from the evolution of the amplitude 
$|A|$ for several 
$R{e_s}$ near the threshold are plotted in figure 4(a) where the linear extrapolation gives the critical value of suspension Reynolds number 
$(R{e_{s,c}})$. The predicted critical value is 
$R{e_{s,c}} = 120$ for the suspension of 
$\mathrm{\epsilon } = 60$ and 
${\phi _b} = 0.1$. We also found the same critical Reynolds number 
$\textrm{(}R{e_{s,c}})$ for 
$\mathrm{\epsilon } = 200$, and it appears that the particle size has no influence on the critical Reynolds number. In contrast, Majji, Banerjee & Morris (Reference Majji, Banerjee and Morris2018) found a reduction in critical suspension Reynolds number by increasing the particle size. This difference might be because, in the experiments, particles are impacted by finite inertia along with the existence of flow curvature in rotating flow, thus undergo inertial migration leading to non-uniform particle distribution that cannot be predicted by the SBM.
Figure 4. (a) Growth rates of perturbations near the critical Reynolds number for 
$\mathrm{\epsilon } = 60$. (b) The derivative of the amplitude logarithm plotted against the square of the amplitude for 
$R{e_s}\; = 122$, 
$\mathrm{\epsilon } = 60$ and 
${\phi _b} = 0.1$.
Note that the critical Reynolds number 
$\textrm{(}R{e_c})$, using LSA, for a pure Newtonian fluid with 
$\eta = 0.877$ is found to be 119.1 (DiPrima, Eagles & Ng Reference Diprima, Eagles and Ng1984). Majji, Banerjee & Morris (Reference Majji, Banerjee and Morris2018) also experimentally observed the primary bifurcation at Re = 120 for a pure fluid flow, where 
$Re = \rho {r_i}{\varOmega _i}d/\mu $. By linear extrapolation of the growth rates, we also obtained Rec = 119.4 for the flow of a pure Newtonian fluid (i.e. 
${\phi _b} = 0$). However, when Re is based on the viscosity of suspending fluid, the suspension critical Re can be computed as 
$R{e_c} = 160.3$ for 
${\phi _b} = 0.1$, meaning that, by adding particles in a viscous fluid, the flow is blocked and it is stabilized. Also, adding particles increases the effective viscosity of suspensions resulting in stronger viscous force leading to the higher shear rate, thereby a larger rotation rate of the inner cylinder is required to induce more intense centrifugal force to reach the critical value (Ali et al. Reference Ali, Mitra, Schwille and Lueptow2002).
The critical Reynolds number 
$\textrm{(}R{e_{s,c}})$ is compared with the previous theoretical and experimental results in table 1. Note that the critical values reported in previous works (Majji, Banerjee & Morris Reference Majji, Banerjee and Morris2018; Gillissen & Wilson Reference Gillissen and Wilson2019; Ramesh, Bharadwaj & Alam Reference Ramesh, Bharadwaj and Alam2019) were converted to our definition for Res accordingly. We employed the maximum packing particle volume fraction as of ϕm = 0.68 in the Krieger's empirical correlation while they considered ϕm = 0.55 in their studies. We have obtained the critical value lower than that reported by Gillissen & Wilson (Reference Gillissen and Wilson2019) who predicted 
$R{e_{s,c}} = 127$ using the LSA of the two-fluid theory. This might be because the interactions of particles, which is significant for concentrated suspensions, have been ignored in their analysis. However, our prediction for 
$R{e_{s,c}}$ at 
$\mathrm{\epsilon } = 60$ shows a good agreement with the experiments of Majji, Banerjee & Morris (Reference Majji, Banerjee and Morris2018). The 
$R{e_{s,c}}$ is also consistent with the value reported by Ramesh, Bharadwaj & Alam (Reference Ramesh, Bharadwaj and Alam2019) for a radius ratio 
$\eta = 0.914$, but it reveals a distinct disagreement with the experimental result performed by Ramesh & Alam (Reference Ramesh and Alam2020) for 
$\eta = 0.94$. This could be due to the curvature effect arising from different radius ratios.
Table 1. Comparison of critical Reynolds numbers (
$R{e_{s,c}}$) in the Taylor–Couette flow of non-Brownian suspension for 
${\phi _b} = 0.1$. Note that, for comparison, we converted the 
$R{e_{s,c}}$ reported in previous literatures by the definition of 
$R{e_s}$ used in this study.
In (3.1), the sign of the Landau constant l determines the type of transition. The transition from the base flow is supercritical if 
$l \gt 0$; otherwise, it is subcritical 
$(l \lt 0)$ (Kang et al. Reference Kang, Meyer, Mutabazi and Yoshikawa2017a, Reference Kang, Meyer, Yoshikawa and Mutabazi2019). The sign of l can be identified from the behaviour of the instantaneous growth rate 
$\textrm{d}\,\textrm{ln}\,|A|/\textrm{d}t$ as a function of 
$|A{|^2}$ at a vanishing 
$|A{|^2}$. We plotted 
$\textrm{d}\,\textrm{ln}\,|A|/\textrm{d}t$ versus 
$|A{|^2}$ at the critical Reynolds threshold 
$R{e_{s,c}}$, in figure 4(b). The intersection with the vertical axis provides the linear growth rate 
$(\sigma )$ of the amplitude 
$|A|$, where the slope at the origin (i.e. 
$|A{|^2} = 0$) determines the nonlinear bifurcation characteristics. This reveals that the instability appears through a supercritical bifurcation 
$(l \gt 0)$. We also found that the transition is supercritical (non-hysteretic) for both particle sizes.
Majji, Banerjee & Morris (Reference Majji, Banerjee and Morris2018) have also studied the hysteresis in flow transitions. They have observed flow transitions and structures by slowly increasing and decreasing Reynolds number (i.e. increasing-Re ramp and decreasing-Re ramp). They showed the hysteresis for the secondary bifurcation, during the increasing-Re ramp, occur at higher Reynolds number (TVF → WVF); however, the transition for the primary instability has not been captured in their experiment.
Ramesh, Bharadwaj & Alam (Reference Ramesh, Bharadwaj and Alam2019) have also investigated the presence of hysteresis by up-sweep (increasing-Re) and down-sweep (decreasing-Re) experiments. They indicated that the secondary bifurcation from TVF to WVF is subcritical (or hysteretic), while the critical Re for the up-sweep is higher than that for the down-sweep. Nonetheless, the primary bifurcation leading to the onset of SVF has not shown the hysteresis for 
${\phi _b} = 0.2$ and the transition was supercritical similar to our results (in particular for 
$\mathrm{\epsilon } = 60$). They also reported that they could not determine the type of the transition for 
${\phi _b} = 0.1$ because different flow patterns have been uncovered during the up-sweep and down-sweep behaviours. Therefore, detail examination of the nature of the transition for different values of 
${\phi _b}$ remains still unknown.
3.3. Flow patterns for small particles 
$(\mathrm{\epsilon } = 200)$
With the onset of the instability, axisymmetric counter-rotating vortices are formed, namely the TVF, in a suspension with 
${\phi _b} = 0.1$ and 
$\mathrm{\epsilon } = 200$ confined between two concentric cylinders (see figure 5). In figure 5(a–c), the spatial and temporal characteristics of the suspension flow are clearly depicted. These vortices are stationary in time and symmetric in the azimuthal direction. As can be seen in figure 5(d), the particles are migrated inside the vortices and they are distributed uniformly in the azimuthal direction. The flow and particle concentration fields at the cross-section of the r–z plane are also presented in figure 5(e–h). The vortices exhibit mirror symmetry that is characterized with a pair of counter-rotating vortices (figure 5e,f). These vortices result in the convective flow and the gradient in the shear, causing the shear-induced particle migration towards the core of the vortices (figure 5g,h). When neutrally buoyant particles are transferred by the convective flow, the particles are simultaneously migrated to the region of the lower shear rates, which is the centre of the vortices.
Figure 5. Flow and particle concentration fields for 
$R{e_s} = 125$, 
$\mathrm{\epsilon } = 200$ and 
${\phi _b} = 0.1$: (a) space–time diagram of radial velocity 
$(u_r)$ at the mid-gap (x = 0.5) and a given 
$\varphi $, (b) iso-surface of 
$- {\lambda _2} = 0.01$ (Jeong & Hussain Reference Jeong and Hussain1995), contours of (c) azimuthal vorticity (
${\omega _\varphi }$) and (d) particle volume fraction versus 
$\phi$, where 
$\phi$ indicates the nondimensionalized length in the azimuthal direction, (
$\phi $) at the central surface (x = 0.5), contours of (e) azimuthal vorticity (
${\omega _\varphi }$) with velocity vectors, (f) azimuthal velocity (
${u_\varphi }$) with streamlines, (g) particle volume fraction (
$\phi $) and (h) local shear rate (
$\dot{\gamma }$) with velocity vectors in an r–z plane (
$\varphi = {\rm \pi}$). Velocity vectors were plotted for every four and two points in the radial (x) and axial (z) directions, respectively.
As 
$R{e_s}$ increases, the WVF is set in the annulus region at 
$R{e_s} = 135$ (figure 6) where the space–time diagram clearly shows an oscillation of vortices in time (figure 6a). The oscillation propagates in both axial and azimuthal directions resulting from the rotating wave of the wavy vortices (figure 6b,c). In contrast to the particle concentrations of the TVF state shown in figure 5(d), here, fewer particles migrate towards the vortices and they are distributed non-uniformly. Figure 6(e–h) display the flow and particle concentration fields in the r–z plane. Counter-rotating vortices, established in the annulus gap, are connected by the convective flow (figure 6e,f). Consequently, fewer particles flow into clockwise-rotating vortices with negative 
${\omega _\varphi }$, migrated by the shear and accumulated in the core of the vortices, where the local shear rate 
$(\dot{\gamma })$ is lower (figure 6g,h).
Figure 6. Flow and particle concentration fields for 
$R{e_s} = 135$, 
$\mathrm{\epsilon } = 200$ and 
${\phi _b} = 0.1$: (a) space–time diagram of radial velocity 
$(u_r)$ at the mid-gap (x = 0.5) and a given 
$\varphi $, (b) iso-surface of 
$- {\lambda _2} = 0.02$, contours of (c) azimuthal vorticity (
${\omega _\varphi }$) and (d) particle volume fraction (
$\phi $) at the central surface (x = 0.5), contours of (e) azimuthal vorticity (
${\omega _\varphi }$) with velocity vectors, (f) azimuthal velocity (
${u_\varphi }$) with streamlines, (g) particle volume fraction (
$\phi $) and (h) local shear rate (
$\dot{\gamma }$) with velocity vectors in an r–z plane (
$\varphi = {\rm \pi}$). Velocity vectors were plotted for every four and two points in the radial (x) and axial (z) directions, respectively.
When 
$R{e_s}$ increases further, at moderately high 
$R{e_s}$, the waviness in vortices becomes stronger where the distinct accumulation of particles does not appear any longer in the particle concentration field of the WVF state (figure 7). The distribution of particles becomes roughly uniform in the bulk flow (figure 7d,h) that is due to the axial oscillation of the WVF state (Majji & Morris Reference Majji and Morris2018). The shear-induced particle migration appears to occur slowly due to the diffusion in the flow of non-colloidal suspensions. For the TVF state, the suspension flow is symmetric in the azimuthal direction and therefore particles have enough time to reach the equilibrium state due to the shear-induced diffusion. On the other hand, the WVF oscillates more rapidly in the axial direction against the diffusion time; thus, particles do not have enough time to undergo the shear-induced diffusion. Consequently, particles are almost uniformly distributed within the annulus gap since the suspension flow pulsates more strongly with the increase of 
$R{e_s}$.
Figure 7. Contours of azimuthal vorticity (
${\omega _\varphi }$) with velocity vectors and particle volume fraction (
$\phi $) in an r–z plane (
$\varphi = {\rm \pi}$) and at the central surface (x = 0.5) for 
$R{e_s} = 150$ (a–d) and 
$R{e_s} = 170$ (e–h), 
$\mathrm{\epsilon } = 200$ and 
${\phi _b} = 0.1$. Velocity vectors were plotted for every four and two points in the radial (x) and axial (z) directions, respectively.
To further analyse the flow structure, we represent mean concentration profiles for several 
$R{e_s}$ plotted in figure 8(a). Here, 
${\langle X\rangle _A}$ denotes the averaged value over a cylindrical area 
$A(r)= 2{\rm \pi} r{L_z}$ at a given r, i.e. 
${\langle X\rangle _A} = (1/A)\int\!\!\!\int {Xr\,\textrm{d}\varphi \,\textrm{d}z} $. As shown earlier, the profile for 
$R{e_s} = 125$ (TVF) has higher values of particle concentration in the middle of the annulus gap where particles accumulate in the core of the vortices. At higher 
$R{e_s}$ where the WVF state appears, the concentration profile gradually becomes flat in the bulk flow. This verifies that particles remain almost uniformly dispersed in the bulk flow. The profiles of the mean azimuthal velocity are also plotted in figure 8(b) for several 
$R{e_s}$. The presence of counter-rotating vortices in the annulus gap leads to the distortion in the mean velocity profile. With the increasing 
$R{e_s}$, the velocity is more distorted because the vortices are intensified in the suspension flow and the slope of the velocity at both walls becomes sharper.
Figure 8. Profiles of (a) mean concentration 
$({\langle \phi \rangle _A})$ and (b) mean azimuthal velocity 
$({\langle {u_\varphi }\rangle _A})$ for 
$\mathrm{\epsilon } = 200$ and 
${\phi _b} = \; 0.1$.
3.4. Flow patterns for large particles $(\mathrm{\epsilon } = 60)$
Above the threshold 
$R{e_{s,c}}$, unlike the suspension of small particles 
$\textrm{(}\mathrm{\epsilon } = 200)$, non-axisymmetric counter-rotating vortices are developed in the annulus of large particle suspensions 
$\textrm{(}\mathrm{\epsilon } = 60)$. Figure 9(a–c) clearly describe the spatial and temporal characteristics of the suspension flow of non-axisymmetric vortices. The SVF characterized by the inclination of vortices and the travelling to the axial direction (Andereck, Liu & Swinney Reference Andereck, Liu and Swinney1986) is definitely verified. Interestingly, bands of higher particle concentration are presented only in vortices of positive 
${\omega _\varphi }$ (figure 9d), whereas they appeared in both vortices for the TVF state of small particle suspensions 
$\textrm{(}\mathrm{\epsilon } = 200)$ (figure 5d). This results from the interaction between the convective flow and the shear-induced migration (figure 9e–h). Non-axisymmetric counter-rotating vortices built in the annulus gap are linked with each other (figure 9e,f), while the axisymmetric ones are separated (figure 5f). As expressed by the streamlines shown in figure 9(f), the suspensions flow from counter-clockwise vortices (with negative 
${\omega _\varphi }$) to the clockwise vortices (with positive 
${\omega _\varphi }$). Particles are then transported into the clockwise-rotating vortices and they are migrated by the shear gradient caused by the vortices (figure 9h). Eventually, particles are concentrated in the core of clockwise-rotating vortices (figure 9g) resulting in the bands of higher concentration observed in figure 9(d).
Figure 9. Flow and particle concentration fields for 
$R{e_s} = 125$, 
$\mathrm{\epsilon } = 60$ and 
${\phi _b} = 0.1$: (a) space–time diagram of radial velocity at the mid-gap (x = 0.5) and a given 
$\varphi $, (b) iso-surface of 
$- {\lambda _2} = \; 0.01$, contours of (c) azimuthal vorticity (
${\omega _\varphi }$) and (d) particle volume fraction (
$\phi $) at the central surface (x = 0.5), contours of (e) azimuthal vorticity (
${\omega _\varphi }$) with velocity vectors, (f) azimuthal velocity (
${u_\varphi }$) with streamlines, (g) particle volume fraction (
$\phi $) and (h) local shear rate (
$\dot{\gamma }$) with velocity vectors in an r–z plane (
$\varphi = {\rm \pi}$). Velocity vectors were plotted for every four and two points in the radial (x) and axial (z) directions, respectively.
By increasing 
$R{e_s}$ to 130 and further to 135, the intensity of spiral vortices becomes stronger (figure 10a,d). As a consequence, more particles accumulate in the vortex core (figure 10b,c,e,f), which is because the intensified vortices cause stronger shear-induced diffusion by generating larger shear gradient in the flow. The distributions of the mean concentration profile 
$({\langle \phi \rangle _A})$ verify the enhancement of accumulation quantitatively (figure 10g), in which the particle concentrations are enhanced by increasing 
$R{e_s}$ in the mid-gap.
Figure 10. Contours of azimuthal vorticity (
${\omega _\varphi }$) with velocity vectors, particle volume fraction (
$\phi $) in an r–z plane (
$\varphi = {\rm \pi}$) and at the central surface (x = 0.5) for 
$R{e_s} = 130$ (a–c) and 
$R{e_s} = 135$ (d–f) with 
$\mathrm{\epsilon } = 60$. (g) Profiles of mean concentration 
$({\langle \phi \rangle _A})$ for 
$\mathrm{\epsilon } = 60$. Velocity vectors were plotted for every four and two points in the radial (x) and axial (z) directions, respectively.
With increase in the Reynolds number to 
$R{e_s} = 140$, an oscillation occurs in the axial travelling wave (figure 11a). These oscillations arise from the appearance of a wavy pattern propagating to the azimuthal direction in spiral vortices (figure 11b,c). We refer to these non-axisymmetric flow structures as WSVF. Figure 11(d) shows the structure of particle distribution in the central surface (x = 0.5) between two cylinders that reveals combined features of SVF and WVF (figures 6d and 9d). As can be seen, more particles are gathered in non-wavy spiral vortices of positive 
${\omega _\varphi }$, displaying distinct bands of higher particle concentration. On the other hand, few particles are collected in oscillating wavy vortices of negative 
${\omega _\varphi }$ where the marked band does not appear in the oscillating vortices. As shown for the SVF state in figure 9, the connected vortices transfer particles to the vortex core of positive 
${\omega _\varphi }$ (figure 11e,f). The transferred particles are concentrated in the vortex core by the shear-induced diffusion resulting in migration of particles to the centre of the vortices with lower shear rates 
$(\dot{\gamma })$ (figure 11g,h). It appears that the particles do not undergo sufficient shear-induced particle migration in wavy vortices due to the axial oscillation as mentioned earlier; therefore, particles are mixed and dispersed more broadly in the region of wavy vortices.
Figure 11. Flow and particle concentration fields for 
$R{e_s} = 140$, 
$\mathrm{\epsilon } = 60$ and 
${\phi _b} = 0.1$: (a) space–time diagram of radial velocity at the mid-gap (x = 0.5) and a given 
$\varphi $, (b) iso-surface of 
$- {\lambda _2} = 0.025$, contours of (c) azimuthal vorticity (
${\omega _\varphi }$) and (d) particle volume fraction (
$\phi $) at the central surface (x = 0.5), contours of (e) azimuthal vorticity (
${\omega _\varphi }$) with velocity vectors, (f) azimuthal velocity (
${u_\varphi }$) with streamlines, (g) particle volume fraction (
$\phi $) and (h) local shear rate (
$\dot{\gamma }$) with velocity vectors in an r–z plane (
$\varphi = {\rm \pi}$). Velocity vectors were plotted for every four and two points in the radial (x) and axial (z) directions, respectively.
A transition to WVF occurs at higher 
$R{e_s} = 165$ (figure 12). As we observed earlier, the axial oscillation of vortices in the WVF limits the accumulation of particles caused by the shear-induced diffusion. Although stronger vortices induced by higher rotation rate of the inner cylinder trigger large shear rates and the gradient in the flow, few particles are sparsely gathered that are mixed and more evenly distributed in the annulus.
Figure 12. Flow and particle concentration fields for 
$R{e_s} = 165$, 
$\mathrm{\epsilon } = 60$ and 
${\phi _b} = 0.1$: (a) space–time diagram of radial velocity at the mid-gap (x = 0.5) and a given 
$\varphi $, (b) iso-surface of 
$- {\lambda _2} = 0.03$, contours of (c) azimuthal vorticity (
${\omega _\varphi }$) and (d) particle volume fraction (
$\phi $) at the central surface (x = 0.5), contours of (e) azimuthal vorticity (
${\omega _\varphi }$) with velocity vectors, (f) azimuthal velocity (
${u_\varphi }$) with streamlines, (g) particle volume fraction (
$\phi $) and (h) local shear rate (
$\dot{\gamma }$) with velocity vectors in an r–z plane (
$\varphi = {\rm \pi}$). Velocity vectors were plotted for every four and two points in the radial (x) and axial (z) directions, respectively.
While particles were more uniformly dispersed in the annulus for the WVF of suspension (figure 12), accumulating particles to the vortex core has been clearly detected in the SVF (figures 9 and 10) and partially in the WSVF (figure 11). For the SVF and WSVF, the suspension flow travels to the axial direction with a constant rotating wave. As a result, the suspended particles have sufficient time to reach their equilibrium states. However, as stated earlier, the particles in the WVF do not reach the steady state situation since the flow oscillates in the axial direction (Majji & Morris Reference Majji and Morris2018). In fact, the time frequency of the oscillations are very small compared with the time scale of the inertial migration of particles.
Figure 13 displays distributions of mean volume fraction and azimuthal velocity between two cylinders for various 
$R{e_s}$ representing different flow structures ranging from CCF to WVF. For the SVF and WSVF, a gentle peak appears in the mean concentration profiles due to the particle accumulation at the middle of the gap (figure 13a). By contrast, the distribution becomes flatter in the WVF because of the roughly uniform dispersion of particles. Similar to the suspension flows of 
$\mathrm{\epsilon } = 200$, the mean velocity profile is more distorted by increasing 
$R{e_s}$ and the slope of the velocity at both walls becomes steeper (figure 13b).
Figure 13. Profiles of (a) mean concentration 
$({\langle \phi \rangle _A})$ and (b) mean azimuthal velocity 
$({\langle {u_\varphi }\rangle _A})$ for 
${\phi _b} = 0.1\; $ suspensions of particle size 
$\mathrm{\epsilon } = 60$.
3.5. Transition scenario
In figure 14, flow patterns for 
${\phi _b} = 0.1$ suspension of both particle sizes (i.e. 
$\mathrm{\epsilon } = 60$ and 200) are summarized and compared with the pure Newtonian flow. Beyond 
$R{e_{s,c}}$, from CCF where the primary bifurcation occurs, various transitions appear for both particle sizes. This is because axisymmetric counter-rotating vortices (TVF) are formed for 
$\mathrm{\epsilon } = 200$ similar to the flow of a pure Newtonian fluid, while non-axisymmetric flow structures (SVF) appear for 
$\mathrm{\epsilon } = 60$. Although we did not discover the ribbons (RIB), which is a stationary pattern arising from the nonlinear interactions between up- and down-propagating spirals different from the experiment of Majji, Banerjee & Morris (Reference Majji, Banerjee and Morris2018), the onsets of TVF and SVF for the corresponding particle sizes 
$(\mathrm{\epsilon })$ are in good agreement with their experiment. It appears that the onset of the suspension primary instability 
$\textrm{(}R{e_{s,c}} = 120)$ is almost similar to that of a pure fluid flow, indicating no particle effects. In Taylor–Couette flow, the destabilization from CCF is induced by the centrifugal force resulting from the curvature of the cylinders and the rotating velocity. Therefore, we can infer that migration of particles has a very weak influence on destabilizing the flow. We also obtained the same 
$R{e_{s,c}}$ for suspensions of both particle sizes, while Majji, Banerjee & Morris (Reference Majji, Banerjee and Morris2018) found a smaller value for the flow of larger particles. As mentioned earlier, this difference might be due to the finite inertial force acting on particles, which is neglected in this study. On the other hand, Ramesh, Bharadwaj & Alam (Reference Ramesh, Bharadwaj and Alam2019), for radius ratio 
$(\eta = 0.914)$ and particle size 
$(\mathrm{\epsilon } = 70)$, found the transition from CCF to SVF without going through RIB under the ramp-down protocol. Note that we did not find the coexisting pattern (WVF + TVF) which was detected by Ramesh, Bharadwaj & Alam (Reference Ramesh, Bharadwaj and Alam2019) during the ramp-up experiment. For suspension of small particles 
$(\mathrm{\epsilon } = 200)$, the stationary TVF transitions to the travelling WVF occur at 
$R{e_s}$ (
$\ge 130$) that is lower than 
$R{e_s} \ge 135$ for a pure Newtonian fluid flow and also 
$R{e_s} \ge 135$ for 
$\mathrm{\epsilon } = 200$ reported by Majji, Banerjee & Morris (Reference Majji, Banerjee and Morris2018). In our study, for suspension of 
$\mathrm{\epsilon } = 60$, an additional non-axisymmetric flow pattern, WSVF, appears between SVF and WVF that is not consistent with the experimental observation. In the experiment reported by Majji, Banerjee & Morris (Reference Majji, Banerjee and Morris2018), an axisymmetric TVF has been detected between SVF and WVF and the WSVF has been found at higher 
${\phi _b}$ (= 0.2, 0.3). These might be because they performed their experiments in a finite cylindrical annulus with a relatively short length (
$\varGamma$ = 20.5), while we have used the periodic boundary conditions in the axial direction. Moreover, flows of suspension could be affected by the inertia of particles that has been neglected in our computations.
Figure 14. Phase diagrams of flow patterns for 
${\phi _b} = 0.1$ of both particle sizes and their comparison with a pure Newtonian flow.
In short, transitions from CCF via TVF to WVF occur for the suspension flow of small particles 
$(\mathrm{\epsilon } = 200)$, similar to the flow of a pure Newtonian fluid. For the suspension of large particles 
$(\mathrm{\epsilon } = 60)$, CCF transitions appear via SVF and WSVF to WVF.
3.6. Friction coefficient
The friction coefficient is a dimensionless measure of the torque that the suspension exerts on the inner cylinder. For the current geometry, it can be given by (Kang, Yang & Mutabazi Reference Kang, Yang and Mutabazi2015)
\begin{equation}{C_{{M_z}}} = \frac{{{M_z}}}{{\rho {\rm \pi}r_i^2{L_z}{{({r_i}{\varOmega _i})}^2}/2}},\end{equation}
where 
${M_z}$ is the axial momentum over the cylindrical surface on the inner cylinder. Figure 15 presents the variations of the friction coefficient 
${C_{{M_z}}}$ for both pure Newtonian fluid and suspensions versus 
$R{e_s}$ together with the values for the CCF with a uniform particle concentration over the annulus given by
\begin{equation}{C_{{M_z},lam}} = \frac{8}{{\eta (1\; + \; \eta )}}Re_s^{ - 1}.\end{equation}
For the CCF, the values are in-line with 
${C_{{M_z},lam}}$, although the distribution of particle concentration is non-uniform over the annulus as shown in figure 3(a). With the development of counter-rotating vortices (i.e. the TVF and SVF states), the friction coefficient 
${C_{{M_z}}}$ increases sharply. The values for suspension of 
$\mathrm{\epsilon } = 200$ (at TVF state) are larger than those of 
$\mathrm{\epsilon } = 60$ (at SVF state). However, the coefficients of suspensions are reduced compared with those of a pure fluid flow. When oscillating waves (i.e. the WSVF and WVF states) occur, the friction coefficient 
${C_{{M_z}}}$ gradually decreases with the increase of 
$R{e_s}$. In particular, for the WSVF state, it decreases more steeply in comparison with that of the WVF state. However, it appears that the friction coefficients have practically the same values for suspensions of different particle sizes and the flow of a pure fluid in the WVF state. Accordingly, we could deduce that the spiral pattern weakens the friction acting on the inner cylinder and the particle size has little influence on the friction coefficient when particles are completely mixed by the oscillating wave in the WVF state.
Figure 15. Friction coefficients (
$C_{M_{z}}$) versus 
$R{e_s}$ for 
${\phi _b} = 0.1$ suspensions of both particle sizes and its comparison with pure fluid. The solid line corresponds to the value for CCF given by (3.4).
In addition, the non-axisymmetric vortices have the azimuthal momentum 
${M_\varphi }$ and the azimuthal friction coefficient 
${C_{{M_\varphi }}}$ (Kang, Yang & Mutabazi Reference Kang, Yang and Mutabazi2015), but we found that the computed values of 
${C_{{M_\varphi }}}$ are very small compared with 
${C_{{M_z}}}$ 
$({C_{{M_\varphi }}}\sim {10^{ - 3}}{C_{{M_z}}})$; therefore, the values of 
${C_{{M_\varphi }}}$ are not shown in the graph.
3.7. Torque and pseudo-Nusselt number
The friction coefficient is related to the dimensionless torque G acting on the inner cylinder that can be stated as (Guillerm et al. Reference Guillerm, Kang, Savaro, Lepiller, Prigent, Yang and Mutabazi2015)
\begin{equation}\; G = \frac{{{C_{{M_z}}}Re_s^2}}{4}{\left( {\frac{\eta }{{1 - \eta }}} \right)^2}.\end{equation}For infinite cylindrical annulus, the dimensionless torque of laminar flow with a uniform particle concentration over the annulus is stated as
\begin{equation}\; {G_{lam}} = \frac{{2\eta \,R{e_s}}}{{(1\; + \; \eta ){{(1 - \eta )}^2}}}.\end{equation} In analogy with thermal convection, the pseudo-Nusselt number (
${N^\omega }$) introduced by Eckhardt, Grossmann & Lohse (Reference Eckhardt, Grossmann and Lohse2007) can be defined as
\begin{equation}{N^\omega } = G/{G_{lam}}.\end{equation}
The pseudo-Nusselt number 
${N^\omega }$ is the analogue of the Nusselt number Nu of the thermal convection. The quantity 
${N^\omega }$ measures how effective the transverse convective angular velocity transport is in terms of the purely molecular transverse transport (Eckhardt, Grossmann & Lohse Reference Eckhardt, Grossmann and Lohse2007); it also represents a direct measure of the dissipation rate 
$(\varepsilon )$ of the kinetic energy through the relation 
${N^\omega } = {(1\; + \; \eta )^2}\varepsilon /4$ (Guillerm et al. Reference Guillerm, Kang, Savaro, Lepiller, Prigent, Yang and Mutabazi2015).
Figure 16(a) displays the variation of the dimensionless torque G with 
$R{e_s}$. Here, for comparison, the values for a pure Newtonian fluid have been computed by
\begin{equation}G = \frac{{{C_{{M_z}}}\,R{e^2}}}{{4{{[{{(1 - 0.1/{\phi _m})}^{ - 1.82}}]}^2}}}{\left( {\frac{\eta }{{1 - \eta }}} \right)^2},\end{equation}
where 
${\phi _m} = 0.68$. The dimensionless torque G increases steadily with 
$R{e_s}$. However, the slope is changed when transitions occur. The values of G for the flow of suspensions are certainly larger than those for a pure fluid flow. This means that suspended particles enhance the torque exerted on the inner cylinder; in other words, more power is required to rotate the cylinder with a given angular velocity. Above the critical values, the torque coefficient is significantly enhanced. As depicted in the inset of figure 16(a), the suspension of 
$\mathrm{\epsilon } = 200$ at the TVF state requires a slightly larger torque compared with the suspension of 
$\mathrm{\epsilon } = 60$ at the SVF state. When the WVF and WSVF states appear, the slope (or increasing rate) of the dimensionless torque G is slightly reduced. Also, for suspension of 
$\mathrm{\epsilon } = 60$ at the WSVF state 
$(140 \le R{e_s} \le 160)$, the values are a little smaller than those of 
$\mathrm{\epsilon } = 200$ at the WVF state; therefore, we could conclude that the spiral vortices travelling to the axial direction cause the reduction of torque. However, the particle size does not impact the torque in the WVF state 
$(R{e_s} \ge 165)$ where particles are almost uniformly dispersed in the annulus region.
Figure 16. Variations of (a) dimensionless torque (G) and (b) pseudo-Nusselt number (
${N^\omega }$) with 
$R{e_s}$ for 
${\phi _b} = 0.1$ suspensions of both particle sizes and their comparison with pure fluid.
Figure 16(b) shows the variation of the pseudo-Nusselt number 
$({N^\omega })$ with 
$R{e_s}$. As mentioned earlier, 
${N^\omega }$ indicates the efficiency of the transfer of the angular velocity in the radial direction by the perturbed flow. With the onset of vortices 
$( \ge R{e_{s,c}})$, the pseudo-Nusselt number increases rapidly with 
$R{e_s}$; however, the values of suspensions are lower than those of a pure fluid flow as plotted in the inset of figure 16(b). Furthermore, the values of 
${N^\omega }$ for the suspension of 
$\mathrm{\epsilon } = 60$ at SVF state are slightly smaller than those of 
$\mathrm{\epsilon } = 200$ at TVF state similar to figure 16(a). Apparently, it could be concluded that particles suspended in a viscous fluid reduce the transverse momentum transfer; moreover, the axial travelling wave of spiral vortices decays the momentum transport in the radial direction. When the transition appears in the flow, the slope of the pseudo-Nusselt number is changed similar to the torque coefficient. Nevertheless, the pseudo-Nusselt numbers of a pure fluid and suspensions are practically the same in the WVF state where particles are almost uniformly distributed in the annulus region of the cylinder.
3.8. Current density of the angular velocity
The pseudo-Nusselt number 
$({N^\omega })$ can also be defined by the current density of the angular velocity 
$(\omega = {u_\varphi }/r)$ (Eckhardt, Grossmann & Lohse Reference Eckhardt, Grossmann and Lohse2007). The angular velocity current 
$({J^\omega })$ is derived by averaging the 
${u_\varphi }$-momentum equation over a cylindrical surface 
$A(r) = 2{\rm \pi} r{L_z}$, defined as
\begin{equation}{J^\omega } = J_c^\omega + J_v^\omega = {r^3}\left[ {\langle {u_r}{\omega}\rangle_A - {{\left\langle {\frac{{\mu {\mu_s}(\phi )}}{\rho }\frac{{\partial \omega }}{{\partial r}}} \right\rangle }_A}} \right].\end{equation}
The transverse transport of angular velocity 
${J^\omega }$ is the sum of the convective transport 
$J_c^\omega $ and of the viscous diffusion transport 
$J_v^\omega $. The quantity 
${J^\omega }$ indicates the transverse current of the azimuthal flow and can be related to the pseudo-Nusselt number as well as the dimensionless torque. Then, the pseudo-Nusselt number can be stated as 
${N^\omega } = {J^\omega }/J_{lam}^\omega $ where 
$J_{lam}^\omega $ is the current density of angular velocity in the laminar flow with a uniform particle concentration over the annulus, which is given by 
$J_{lam}^\omega = 2R{e_s}/\eta (1\; + \; \eta )$.
In figure 17, the variations of the angular velocity current 
${J^\omega }$ and its constitutive parts (
$J_c^\omega $ and 
$J_v^\omega $) are represented for various flow states. The current density is conserved in the radial direction 
$\textrm{(d}{J^\omega }\textrm{/d}r = 0)$. The contribution of 
$J_v^\omega $ on the total transport of angular velocity is dominant at low 
$R{e_s}$, while the portion of 
$J_c^\omega $ is enhanced by increasing 
$R{e_s}$, especially in the middle of the annular gap where the cores of counter-rotating vortices exist. As mentioned earlier, the axial travelling wave of spiral patterns reduces the transverse momentum transfer, leading to the decrease of 
${N^\omega }$. This could also be interpreted by the profiles of 
$J_c^\omega $. As can be seen in figure 17, the quantity of the convective transport 
$J_c^\omega $ is slightly reduced when the spiral pattern appears, i.e. for the SVF and WSVF states. Therefore, we could conclude that the axial travelling wave decays the convective momentum transfer in the radial direction; in consequence, it results in the reduction of the pseudo-Nusselt number and torque acting on the inner cylinder.
Figure 17. Radial distributions of angular velocity current (
${J^\omega }$) and the constitutive parts 
$(J_c^\omega ,J_v^\omega )$ for 
$R{e_s} = 125$; (a) 
$\mathrm{\epsilon } = 200$ (TVF) and (b) 
$\mathrm{\epsilon } = 60$ (SVF) and for 
$R{e_s} = 155$; (c) 
$\mathrm{\epsilon } = 200$ (WVF) and (d) 
$\mathrm{\epsilon } = 60$ (WSVF).
4. Conclusion
The Taylor–Couette flow of neutrally buoyant, non-colloidal and non-Brownian rigid particles has been numerically examined in this study. Numerical simulations for 
${\phi _b} = 0.1$ suspension of both particle sizes (i.e. 
$\mathrm{\epsilon }\;( = \,d\textrm{/}a) = 60$ and 
$\mathrm{\epsilon }\;( = \,d\textrm{/}a) = 200$) undergoing the shear-induced particle migration have been conducted using rheological constitutive models known as SBM. Same as the experiments by Majji, Banerjee & Morris (Reference Majji, Banerjee and Morris2018) and Ramesh, Bharadwaj & Alam (Reference Ramesh, Bharadwaj and Alam2019), we also considered that the flow of suspensions is in the Stokes flow regime; thus, the inertia of particles is neglected 
$(R{e_p} \ll 1)$, and the Péclet number (Pe) is very large 
$(Pe = \infty )$. However, we assumed that the shear-induced migration is dominant in the flow of suspensions while the inertial migration is prevalent in the experiments reported by Majji, Banerjee & Morris (Reference Majji, Banerjee and Morris2018) and Ramesh, Bharadwaj & Alam (Reference Ramesh, Bharadwaj and Alam2019). We have varied the Reynolds number of suspensions 
$(R{e_s})$ based on the rotating angular velocity of the inner cylinder and the effective viscosity of suspension. The flow and particle concentration fields have been investigated to clarify the effect of non-colloidal particles on the flow instability and transitions in suspension flows.
At low Reynolds number, below the critical value, the CCF of suspension occurs where the velocity and particle concentrations are basically linear depending only on the radial coordinate (r). It also reveals that the semidilute suspension 
$({\phi _b} = 0.1)$ has little impact on the azimuthal velocity over the narrow-gap where the inertial migration of particles is neglected.
The onset of the primary bifurcation has been estimated by employing the Landau model. We found that the particle size does not affect the critical values of suspension Reynolds number where the predicted critical Reynolds number was found to be 
$R{e_{s,c}} = 120$ for 
${\phi _b} = 0.1$ suspensions of both particle sizes (i.e. 
$\mathrm{\epsilon } = 60, 200$). This implies that particles dispersed in a viscous fluid lead to stabilizing the flow, meaning that suspensions require a larger rotation rate of the cylinder to cause instability in the flow. The obtained critical values show a good agreement with the previous experiments reported by Majji, Banerjee & Morris (Reference Majji, Banerjee and Morris2018) in a finite cylindrical annulus gap; therefore, it is expected that neither the length of the cylinder nor the inertial migration might impact the onset of the primary instability. We have also evaluated the type of the primary instability and we found that the transition appears through a supercritical (or non-hysteretic) bifurcation for both suspensions’ particle sizes (i.e. 
$\mathrm{\epsilon } = 60, 200$) similar to the experiments.
By varying Reynolds number 
$(R{e_s})$, different transition scenarios have been observed for 
${\phi _b} = 0.1\;$suspensions of both particle sizes. While we confirmed the non-axisymmetric patterns, such as the wavy vortices and spiral vortices to occur for 
${\phi _b} = 0.1$ suspensions, in contrast to Majji, Banerjee & Morris (Reference Majji, Banerjee and Morris2018), we did not observe the RIB state. This discrepancy might be due to the impact of the finite length of the cylinders and/or the influence of the microstructure of the particles that was not considered in this analysis. Also, in the experiment, particles undergo the inertial migration leading to non-uniform particle distribution, while the SBM models the shear-induced migration that occurs at very low Reynolds number. In general, the present work and the experiments reported by Majji, Banerjee & Morris (Reference Majji, Banerjee and Morris2018) confirm that different states can appear in the Taylor–Couette suspension flow even if the aspect ratio and the radius ratio are 
$\varGamma = 20.5$ and 
$\eta = 0.877$, respectively.
For suspensions of small particles 
$(\mathrm{\epsilon } = 200)$, the suspension flow transitions from the CCF via TVF to WVF, i.e. CCF → TVF → WVF, similar to the flow of a pure Newtonian fluid. In the TVF state, particles migrate towards the centre of counter-rotating vortices due to the shear-induced diffusion and accumulate in the core of vortices. However, for the WVF, particle migration is not significant due to the axial oscillation of the WVF. Although few particles are collected in vortex cores when weak oscillating waves appear, particles are almost uniformly dispersed in the annulus gap at higher Reynolds numbers. By contrast, for suspensions of large particles 
$(\mathrm{\epsilon } = 60)$, as the Reynolds number increases, transitions occur following the sequence of CCF → SVF → WSVF → WVF. Unlike the suspension flow of small particles, additional non-axisymmetric counter-rotating vortices, i.e. the SVF and WSVF states, appear between the CCF and WVF states. In the SVF state, particles migrated by the shear are transferred by the convective flow into vortex cores of positive azimuthal vorticity 
$({\omega _\varphi })$, which links counter-rotating vortices resulting in particle accumulation only in vortices of positive 
${\omega _\varphi }$. For the WSVF state, two features of the SVF and WVF states are presented in the particle concentration field. More particles are migrated in non-wavy spiral vortices of positive 
${\omega _\varphi }$, while they are mixed and dispersed in wavy vortices because of the axial oscillation.
We also evaluated the friction and torque coefficients for 
${\phi _b} = 0.1$ suspensions of both particle sizes and compared our results with the values for a pure Newtonian fluid. When counter-rotating vortices occur, the coefficients are considerably enhanced. Beyond the critical Reynolds number, the friction coefficient 
${C_{{M_z}}}$ rises with the Reynolds number in the flow regimes of non-oscillating waves (i.e. the TVF and SVF states). However, it gradually decreases with the increase of the Reynolds number after the appearance of oscillating waves (i.e. the WSVF and WVF states). By contrast, the dimensionless torque G and the pseudo-Nusselt number 
${N^\omega }$ steadily increase as the Reynolds number rises. It appears that suspended particles significantly enhance the torque, which is the required power to rotate the inner cylinder in comparison with the values of a pure fluid flow, while the particles cause the reduction of 
${C_{{M_z}}}$ and 
${N^\omega }$ above the threshold of 
$R{e_{s,c}}$. Moreover, the axial travelling wave of spiral vortices reduces both the friction coefficient 
$({C_{{M_z}}})$ and the torque (G) acting on the inner cylinder. It also weakens the transverse momentum transfer 
$({N^\omega })$ in the radial direction by decreasing the convective transport of angular velocity. However, the coefficients are practically the same in the WVF regime where particles are almost uniformly dispersed in the annulus gap by the axial oscillation of WVF. It appears that the axial oscillating vortices do not allow the particles to undergo the shear-induced diffusion.
Our study represents the very first numerical analysis of suspensions in a Taylor–Couette flow using continuum models of suspensions including SBM. It also represents a new approach to numerically deal with the flow instability and transitions that take place in the rotating flow of non-Brownian, non-colloidal suspensions undergoing shear-induced diffusion with various aspect ratio, radius ratio, gap width and particle size to understand the role of particles in Taylor vortices (Ramesh & Alam Reference Ramesh and Alam2020). The present work can be served as a first step towards extending the analysis of dense suspensions in the Taylor–Couette geometry using SBM where the inertial migration of particles is neglected. Further experiments should be performed to verify the inhomogeneity in concentration distribution observed herein. Examining whether the observed patterns during both ramp-up and ramp-down are associated with hysteresis and/or new patterns might occur for various concentrations is the subject of our current investigation.
Acknowledgements
This research was supported by National Science Foundation award no. 1854376 and Army Research Office award no. W911NF-18-1-0356.
Declaration of interests
The authors report no conflict of interest.
