2.1. Description of the system

We investigate here a dilute suspension of particles in a pressure-driven channel flow, a schematic of which is shown in figure 1. The impact of this suspension on the flow is characterized by a two-way coupling, modelled using the formulation of Saffman (Reference Saffman1962), specifically through the application of Stokes drag, with the addition of viscosity variations due to particle concentration. The viscosity variation terms are derived from Govindarajan (Reference Govindarajan2004). The particulate suspension is treated as a continuous medium whose dynamics is describable by a field equation. The momentum balance and continuity equations for the fluid respectively are

(2.1)\begin{gather} \rho_{f}\left(\frac{\partial \boldsymbol{u_d}}{\partial t_d}+ \boldsymbol{u_d} \boldsymbol{\cdot} \boldsymbol{\nabla_d}\boldsymbol{u_d}\right)=-\boldsymbol{\nabla_d} p_d + \boldsymbol{\nabla_d}[\mu^{tot}_d\boldsymbol{\cdot} (\boldsymbol{\nabla_{d}}\boldsymbol{u_{d}}+ (\boldsymbol{\nabla_{d}}\boldsymbol{u_{d}})^T)] +K N(\boldsymbol{v_d}-\boldsymbol{u_d}), \end{gather}

(2.2)\begin{gather}\boldsymbol{\nabla_d}\boldsymbol{\cdot}\boldsymbol{u_d}=0, \end{gather}

while the particle suspension satisfies momentum balance and continuity respectively given by

(2.3)\begin{gather} m N\left(\frac{\partial \boldsymbol{v_d}}{\partial t_d}+ \boldsymbol{v_d} \boldsymbol{\cdot} \boldsymbol{\nabla_d}\boldsymbol{v_d}\right)=-KN (\boldsymbol{v_d}-\boldsymbol{u_d}), \end{gather}

(2.4)\begin{gather}\frac{\partial N}{\partial t_d}+\boldsymbol{\nabla_d}\boldsymbol{\cdot}(N \boldsymbol{v_d})=0. \end{gather}

Here, the subscript $d$ represents a dimensional variable and $\rho _{f}$ is the dimensional density of the fluid. The total dimensional viscosity $\mu ^{tot}_d=\mu _f+\mu _p$, where $\mu _f$ is the dimensional viscosity of the fluid and $\mu _p$ is the contribution to viscosity due to the particles. Also, $m$ and $\tau =m/K$ are the mass and relaxation time of a spherical dust particle, $N$ is their number density per unit volume. The quantity $K$ is the drag coefficient given by $6{\rm \pi} r\mu _f$ for a sphere of radius $r$. We establish the $x$-axis to be aligned with the channel centreline, the $y$-axis to be oriented in the wall-normal direction and the $z$-axis to be oriented perpendicular to the plane of the figure. The fluid velocity is given by $\boldsymbol {u}=(u_x,u_y,u_z)$. The particles are assumed to be continuously distributed in the flow, and to have a continuous velocity variation in space and time, and their dynamics may thus be described as a field, with $\boldsymbol {v}=(v_x,v_y,v_z)$. Our analysis thus precludes situations of diverging number density and of particle collisions. The mass fraction of particles in the suspension is

(2.5)\begin{equation} f=mN/\rho_{f}=\frac{4{\rm \pi} Nr^3}{3}\frac{\rho_p}{\rho_f}. \end{equation}

The density of the solid making up the particles is $\rho _p$. Unless otherwise specified, we work in the limit of $\rho _p/\rho _f \to \infty$, so the volume fraction occupied by the particles is negligible. We non-dimensionalize equations (2.1)–(2.4) using the channel-centreline mean velocity $U_m$ (i.e. the steady-state fluid velocity at the centre of the channel), the half-channel width $H$ and the viscosity $\mu _f$ of the fluid as scales. In dimensionless coordinates, the channel walls are positioned at $y=\pm 1$, with the suspended particles concentrated around $y=\pm a_p$. We prescribe a mean dust mass-fraction profile

(2.6)\begin{equation} \bar{f}(y)=f_{max}\left[\exp\left\{\frac{-(y-a_p)^2}{2\sigma^2}\right\}+ \exp\left\{\frac{-(y+a_p)^2}{2\sigma^2}\right\}\right],\end{equation}

with particles concentrated in two layers of thickness $\sigma$, but our numerical method is general and suitable for any desired particle concentration profile. The location $a_p$ of the maximum in particle concentration is an important parameter. If the same particles were to be uniformly distributed across the channel, the loading would be given by

(2.7)\begin{equation} f_{ave}=\frac{\displaystyle \int_{-1}^{{+}1}\bar{f}(y;\sigma,a_p){{\rm d} y}}{\displaystyle \int_{-1}^{{+}1}{{\rm d} y}} \simeq \sqrt{2{\rm \pi}}f_{max} \sigma. \end{equation}

The Reynolds and Stokes numbers, which will emerge out of the non-dimensionalization, are given respectively by

(2.8)\begin{equation} R \equiv \frac{HU_m}{\mu_f/\rho_{f}} \quad {\rm and} \quad S \equiv \frac{\tau}{\rho_{f} H^{2}/\mu_f} = \frac{2}{9}\frac{r^2}{H^2} \frac{\rho_p}{\rho_f}.\end{equation}

In terms of the average density in the flow, we may also define an effective Reynolds number $R_{eff}=(1+\bar {f}_{ave})R$. These two quantities, the Reynolds number $R$ and the Stokes number $S$, along with the thickness $\sigma$ of the particle-laden layer, the mass loading for a given $\sigma$ as measured by $f_{max}$ and the location $a_p$ of the maximum in particle concentration, are the parameters which determine this problem.

Figure 1. Schematic of the flow under consideration. The walls are situated at $y=\pm 1$, with the green curve and red arrows representing the mean velocity profile, $U(y) = 1 - y^{2}$. The particles are concentrated around $y = \pm a_p$ within a band of size $\sigma$. The mean particle mass fraction, $\bar f$, given by (2.6), is depicted on the right. Note that the volume fraction that heavy particles occupy will be much smaller.

2.2. Linear stability equations

After non-dimensionalizing, we split all quantities in (2.1)–(2.4) into their basic and fluctuating parts, as $\boldsymbol {u} = \boldsymbol {U}+\hat {\boldsymbol {u}}$, $\boldsymbol {v} = \boldsymbol {U}+\hat {\boldsymbol {v}}$, $p = P+\hat {p}$, $f=\bar {f}+\hat {f}$ and $\mu ^{tot}=\bar \mu + \hat \mu$. Here, a hat represents a perturbation quantity, while an upper case or overbar denotes a mean quantity. In parallel shear flows, we have $\boldsymbol {U}=U(y)\boldsymbol {e_x}$, where $\boldsymbol {e_x}$ is a unit vector in the streamwise direction, and $\bar f=\bar f(y)$. For small particulate volume fraction, the local viscosity is linearly related to the local particle concentration, as

(2.9)\begin{equation} \mu^{tot}_d=\mu_f\left[1+\dfrac{f}{\gamma}\right],\end{equation}

where $\gamma \propto \rho _{p}/\rho _{f}$. In accordance with Einstein's law, we take the proportionality constant to be $0.4$. In the limit of infinite $\gamma$, the dimensional viscosity remains at $\mu _f$ everywhere. The viscosity is non-dimensionalized by the viscosity of the pure fluid, $\mu _{f}$. The mean viscosity is described as $\bar {\mu }=1+\bar {f}/\gamma$, and the perturbation viscosity is $\hat {\mu }=\hat {f}/\gamma$, obtained from (2.9).

Upon linearization of (2.1)–(2.4) we have

(2.10)\begin{align} \left(\frac{\partial \hat{\boldsymbol{u}}}{\partial t}+ \hat{\boldsymbol{u}} \boldsymbol{\cdot} \boldsymbol{\nabla}\boldsymbol{U}+\boldsymbol{U} \boldsymbol{\cdot} \boldsymbol{\nabla}\hat{\boldsymbol{u}}\right) & =-\boldsymbol{\nabla}\hat{p}+\frac{1}{R}[\boldsymbol{\nabla} \bar{\mu}\boldsymbol{\cdot}(\boldsymbol{\nabla}\hat{\boldsymbol{u}}+( \boldsymbol{\nabla}\hat{\boldsymbol{u}} )^T)+ \boldsymbol{\nabla}\hat{\mu}\boldsymbol{\cdot} (\boldsymbol{\nabla}\boldsymbol{U}+ (\boldsymbol{\nabla}\boldsymbol{U} )^T) \nonumber\\ &\quad +\hat{\mu} \nabla^2 \boldsymbol{U} +\bar{\mu}\nabla^2 \hat{\boldsymbol{u}} ] +\frac{\bar{f}}{SR}(\hat{\boldsymbol{v}}-\hat{\boldsymbol{u}}), \end{align}

(2.11)\begin{gather} \boldsymbol{\nabla}\boldsymbol{\cdot}\hat{\boldsymbol{u}}=0, \end{gather}

(2.12)\begin{gather} \left(\frac{\partial \hat{\boldsymbol{v}}}{\partial t}+ \hat{\boldsymbol{v}} \boldsymbol{\cdot} \boldsymbol{\nabla}\boldsymbol{U}+\boldsymbol{U} \boldsymbol{\cdot} \boldsymbol{\nabla}\hat{\boldsymbol{v}}\right)={-} \frac{1}{SR}(\hat{\boldsymbol{v}}-\hat{\boldsymbol{u}}), \end{gather}

(2.13)\begin{gather} \frac{\partial \hat{f}}{\partial t}+\boldsymbol{\nabla}\boldsymbol{\cdot} (\,\bar{f} \hat{\boldsymbol{v}})+\boldsymbol{\nabla}\boldsymbol{\cdot} (\,\hat{f} \boldsymbol{U})=0. \end{gather}

We start by performing a normal-mode analysis, considering single Fourier modes for the perturbation quantities ($\hat {\boldsymbol {u}}$, $\hat {\boldsymbol {v}}$, $\hat {p}$, $\hat {f}$, $\hat {\mu }$) in both the $x$-direction and $z$-direction, as well as in time. In other words, the perturbation quantities are written in normal-mode form, with each mode given by

(2.14)\begin{equation} (\hat{\boldsymbol{u}}, \hat{\boldsymbol{v}}, \hat{f},\hat{\mu})=\tfrac{1}{2}[(\boldsymbol{u}(y),\boldsymbol{v}(y), f(y),\mu(y))\exp\{{\rm i}\alpha(x-ct) + {\rm i}\beta z\}+{\rm c.c.} ].\end{equation}

We may now express the particle velocity in terms of the flow velocity using (2.12), to get

(2.15)\begin{equation} (v_{x},v_{y},v_{z})=({\mathcal{M}}u_{x}-SR{\mathcal{M}}^{2}U^{\prime}u_{y},{\mathcal{M}}u_{y},{\mathcal{M}}u_{z}), \end{equation}

where

(2.16)\begin{equation} {\mathcal{M}}=\frac{1}{1+{\rm i}\alpha SR(U-c)}. \end{equation}

Additionally, by taking the divergence of (2.10), we write the pressure Laplacian in terms of the velocity field as

(2.17)\begin{align} -\nabla^2 p&=2{\rm i}\alpha U^{\prime}u_{y}-[ 2\bar{\mu}^{\prime}\nabla^{2}u_{y}+2\bar{\mu}^{\prime\prime}Du_{y}+2{\rm i}\alpha U^{\prime}D\mu+2{\rm i}\alpha U^{\prime\prime}\mu]\nonumber\\ &\quad -\frac{1}{SR}[\,\overline{f^{\prime}}({\mathcal{M}}-1)u_{y}+\bar{f}{\mathcal{M}}^{\prime}u_{y}-{\rm i}\alpha SR{\mathcal{M}}^{2}U^{\prime}\bar{f}u_{y}]. \end{align}

We apply the operator $\nabla ^2$ to the $y$ component of the vector equation (2.10), use (2.15) and (2.17) and divide throughout by $-{\rm i}\alpha ^2$ to express the resulting equation in terms of the variables $u_y$ and $\mu$

(2.18)\begin{align} -{\rm i}c \nabla^{2}\dfrac{u_{y}}{-{\rm i}\alpha } &={\rm i}[U^{\prime\prime} -U\nabla^{2}]\dfrac{u_{y}}{-{\rm i}\alpha }+\frac{1}{\alpha R}[\bar{\mu}'' ( -\nabla^2 + 2 D^2 ) + 2 \bar{\mu}' D \nabla^2 + \bar{\mu} \nabla^4 ] \dfrac{u_{y}}{-{\rm i}\alpha } \nonumber\\ &\quad +\frac{1}{\alpha R} [U^{\prime\prime\prime} + 2 U^{\prime\prime} D - U^{\prime} \nabla^2 + 2 U^{\prime} D^2] \mu \nonumber\\ &\quad +{\rm i}[({\mathcal{M}}^{2}U^{\prime}\bar{f})^{\prime}-(U-c){\mathcal{M}}\bar{f}^{\prime}D-(U-c){\mathcal{M}}\bar{f}\nabla^{2} ]\dfrac{u_{y}}{-{\rm i}\alpha}. \end{align}

Using (2.15) and the continuity equation (2.11), we can convert equation (2.4) into an expression that includes the variables $u_{y}$ and $\mu$, leading, after dividing throughout by $-{\rm i}\alpha$, to

(2.19)\begin{equation} (U-c)\gamma\mu+\left[- RS{\mathcal{M}}^{2}U^{\prime}\bar{f}^{\prime}+ \frac{({\mathcal{M}}\bar{f})^{\prime}}{{\rm i}\alpha}\right] u_{y}=0. \end{equation}

Equations (2.18) and (2.19) represent the linear stability equations for three-dimensional perturbations. For a clean parallel shear flow without particles, Squire (Reference Squire1933) had shown that, for every three-dimensional perturbation mode satisfying the stability equations, there exists a corresponding two-dimensional perturbation mode at a lower Reynolds number, displaying the same growth rate. Saffman (Reference Saffman1962) had shown that Squire's theorem applies in the case of a dusty channel with uniform particle loading. If we substitute $(\alpha ^2+\beta ^2)=\alpha _{2D}^2$, $\alpha R=\alpha _{2D} R_{2D}$ and $u_{y}/\alpha =u_{y,2D}/\alpha _{2D}$ into the aforementioned equations (2.18) and (2.19), these equations become equivalent to those of a two-dimensional system with the wavenumber denoted as $\alpha _{2D}$, the Reynolds number as $R_{2D}$ and the velocity eigenfunction $u_{y,2D}$. We thus show that Squire's theorem may be extended to dusty channels with inhomogeneous loading, including viscosity variation as well. The last, for viscosity stratified shear flow without particles, was shown by Jose (2024, private communication). Thus, for two-dimensional perturbations, these equations become (2.20) and (2.21). Therefore, while a non-modal study would require us to study three-dimensional perturbations, since our purpose is to obtain linear instability at low Reynolds number, it is sufficient to perform a two-dimensional calculation, by setting $\beta =0$ in (2.14). In fact, results for any given three-dimensional single mode may be obtained directly from an equivalent two-dimensional one by simple rescaling.

The two-dimensional equations for linear perturbations, after appropriate elimination and reduction, can be written in terms of the perturbation streamfunction $\psi (y)$ and the perturbation viscosity $\mu (y)$ as

(2.20) \begin{align} &[(U_{*} - c)(D^2 - \alpha^2) - U_{*}'']\psi + D(J\bar{f}' \psi) = \frac{1}{{\rm i}\alpha R}[\bar{\mu}(D^2 - \alpha^2)^2 + 2\bar{\mu}' D^3 + \bar{\mu}'' D^2\nonumber\\ &\quad - 2\alpha^2 \bar{\mu}^{\prime} D + \alpha^2 \bar{\mu}'']\psi + \frac{1}{R}[U' D^2 + 2U'' D + U''' + \alpha^2 U']\mu, \end{align}

and

(2.21)\begin{equation} -(U-c)\gamma\mu+[-{\rm i}\alpha RS{\mathcal{M}}^{2}U^{\prime}\bar{f}^{\prime}+ ({\mathcal{M}}\bar{f})^{\prime}]\psi=0,\end{equation}

where

(2.22)\begin{gather} U_{*}\equiv U+J\bar{f}, \quad {\mathcal{M}}= \frac{1}{1+{\rm i}\alpha (U-c)SR},\quad J=(U-c){\mathcal{M}}, \end{gather}

(2.23a)\begin{gather} u_y=-{\rm i}\alpha\psi, \end{gather}

(2.23b)\begin{gather}u_{x}=D\psi , \end{gather}

(2.23c)\begin{gather}v_{x}={\mathcal{M}} u_{x}-({\mathcal{M}}^{2}SRU^{\prime})u_{y} , \end{gather}

(2.23d)\begin{gather}v_y={\mathcal{M}} u_y . \end{gather}

The operator $D$ is defined as $D = {\rm d}/{\rm d}y$, and a prime denotes a derivative in $y$ of a mean quantity. Note that $v_x$ and $v_y$ in (2.23c) and (2.23d) can be expressed in terms of $\psi$ using (2.23a) and (2.23b). The boundary conditions are

(2.24)\begin{equation} \psi(y={\pm} 1)=D\psi(y={\pm} 1)=\mu(y=\pm1)=0.\end{equation}

For given mean flow $U(y)$, streamwise wavenumber $\alpha$, base particle loading $\bar f(y)$, particle to fluid density ratio and fixed Reynolds and Stokes numbers, (2.20)–(2.24) define an eigenvalue problem, which yields a spectrum of eigenvalues $c$ and corresponding eigenfunctions, $(\psi (y),\mu (y))$. If even one eigenvalue has a positive imaginary part, i.e. $c_{im}> 0$, we have an exponential growing mode.

In the limit of $\gamma \to \infty$, we have $\bar \mu =1$ and $\mu =0$, so (2.20) becomes

(2.25) \begin{equation} [(U_{*}-c)(D^2-\alpha^2)-U_{*}^{\prime\prime}]\psi+(J \bar{f}')^{\prime}\psi+(J\bar{f}')D\psi=\frac{1}{{\rm i}\alpha R}(D^2-\alpha^2)^2 \psi,\end{equation}

and is now decoupled from (2.21). When there is no particulate suspension, we have $\bar {f} = 0$, and the system (2.25) reduces to the well-known Orr–Sommerfeld equation

(2.26)\begin{equation} \left[U(D^2 - \alpha^2) - U'' + \frac{i}{\alpha R} (D^2 - \alpha^2)^2 \right] \psi = c (D^{2}-\alpha^{2})\psi.\end{equation}

In the case of a homogeneous suspension ($\,\bar {f} = \text {const}.$) along with $\gamma \to \infty$, the system reduces to that of Saffman (Reference Saffman1962).

2.3. Balance of perturbation kinetic energy

Whenever the flow is unstable, there is an exponential increase in perturbation kinetic energy. It is useful to derive the positive and negative contributors to this quantity. To do this, we multiply the linear equations for the fluid flow (in $\hat {\boldsymbol {u}}$) and for the particulate flow (in $\hat {\boldsymbol {v}}$), given as (2.10) and (2.12), by the respective complex conjugates $\hat {\boldsymbol {u}}^{*}$ and $\hat {\boldsymbol {v}}^{*}$. Upon averaging over a wavelength in the streamwise direction, we derive the evolution of perturbation kinetic energy $\hat {E}$ to be described by

(2.27) \begin{align} \partial_t \int \hat{E} \, {\rm d} V &=- \int \frac{\partial U_i}{\partial x_j} \hat{u}_i \hat{u}_j\, {\rm d} V -\frac{1}{R} \int \bar{\mu} |\partial_i \hat{u}_j|^2 \, {\rm d} V \nonumber\\ &\quad -\int \bar{f}\frac{\partial U_i}{\partial x_j} \hat{v}_i \hat{v}_j\, {\rm d} V -\frac{1}{S R} \int \bar{f} |\hat{u}_i - \hat{v}_i|^2 \,{\rm d} V \nonumber\\ &\quad -\dfrac{1}{R}\int\frac{\partial^2 \bar{\mu}}{\partial x_j \partial x_i} \hat{u}_i \hat{u}_j\, {\rm d} V -\dfrac{1}{R}\int\frac{\partial U_{j}}{\partial x_{i}} \hat{\mu} (\partial_j \hat{u}_i+\partial_i \hat{u}_j)\, {\rm d} V, \end{align}

where

(2.28)\begin{equation} \hat{E}=\tfrac{1}{2}(\hat{u}_{i}^{2}+\bar{f}\hat{v}_{i}^{2}), \end{equation}

and $V$ indicates a volume of fluid extending from wall to wall and over one perturbation wavelength in the streamwise direction. We then introduce the normal-mode forms of the perturbations, given by (2.14), into (2.27), and average over the streamwise direction $x$, to get

(2.29)\begin{align} 2\alpha c_{im} \int E \, {{\rm d} y} &={-}\frac{1}{4} \int \frac{\partial U_i}{\partial x_j} (u_i u_j^* + u_i^* u_j) \, {{\rm d} y} -\frac{1}{2R} \int \bar{\mu}|\partial_i u_j|^2 \, {{\rm d} y} \nonumber\\ &\quad -\frac{1}{4} \int \bar{f}\frac{\partial U_i}{\partial x_j} (v_i v_j^* + v_i^* v_j) \, {{\rm d} y} -\frac{1}{2S R} \int \bar{f} |u_i - v_i|^2 \, {{\rm d} y} \nonumber\\ &\quad -\frac{1}{4R}\int \frac{\partial^2 \bar{\mu}}{\partial x_j \partial x_i} ( u_i^* u_j+ u_i u_j^*){{\rm d} y} \nonumber\\ &\quad -\frac{1}{4R}\int \frac{\partial U_{j}}{\partial x_{i}} \{\mu (\partial_j u_i^*+\partial_i u_j^*)+\mu^{*} (\partial_j u_i+\partial_i u_j)\}{{\rm d} y} \nonumber\\ & \equiv \int (W_+-W_-+W_{p+}-W_{p-}+W_{\mu,1}+W_{\mu,2}){{\rm d} y}, \end{align}

where

(2.30)\begin{equation} E(y)=\tfrac{1}{4}(u_{i}(y)^{2}+\bar{f}v_{i}(y)^{2}), \end{equation}

and $W_+(y)$ and $W_-(y)$, respectively, are the production and dissipation of perturbation kinetic by the fluid, while $W_{p+}(y)$ and $W_{p-}(y)$, respectively, give the production and dissipation of perturbation kinetic energy of the particles. The last two terms $W_{\mu,1}$ and $W_{\mu,2}$ arise due to viscosity stratification.

2.4. Numerical method

We employ the Chebyshev spectral collocation method to discretize the system given by (2.20)–(2.24) at $n$ discrete points in the domain. The Chebyshev collocation points, defined as $y_{Cheb, j}=\cos [{\rm \pi} j/(n-1)], j=0,1,2,3,\ldots,n-1$, are naturally clustered close to the walls. Such a discretization would resolve the near-wall region well, where variations are large. But, for a small number of collocation points it would leave the particle layer, where variations are also large, not well resolved. In order to get results insensitive to the number of collocation points, we employ a stretching function to cluster a sufficient number of grid points into the particle-laden layer. Such a stretching function was used in Govindarajan (Reference Govindarajan2004) in a different context, and works well in the present situation as well. It is given by

(2.31)\begin{equation} y_j=\frac{a}{\sinh (b y_b)}[\sinh \{(y_{Cheb,j}-y_b) b\}+\sinh (b y_b)],\end{equation}

where

(2.32)\begin{equation} y_b=\frac{1}{2b} \log \left[\frac{1+(e^b-1) a}{1+(e^{-b}-1) a}\right], \end{equation}

is a constant, $a$ signifies the location around which clustering is desired and $b$ serves to determine the level of clustering. Once written in discrete form, each boundary condition may be applied by replacing one row of the discrete system appropriately. To solve (2.25) after discretization, we utilize the LAPACK FORTRAN package. For our all simulations, we use $n=81$, and verify our answers with $n=121$. At this resolution, the results are insensitive to the number of grid points, as well as to whether stretching is employed or not. But stretching improves the physical appearance of the eigenfunctions.

Since the chosen mass-fraction profile corresponds to the clustering of particles in the vicinity of $y=\pm a_p$, we select $a$ to be equal to $\pm a_p$ and set $b$ to the value of $2$ or $4$ in (2.31). We obtain eigenvalues correct to five decimal places for the most part, and at least to four places everywhere. At Stokes numbers of $10^{-2}$ or higher, however, the accuracy drops to three decimal places, and we do not venture into this regime to make our conclusions.

To validate our approach, we first perform computations using a uniform particle profile across the channel. Figure 2 shows neutral stability boundaries provided by Klinkenberg et al. (Reference Klinkenberg, De Lange and Brandt2011), compared with the present computations. The agreement is excellent for two different particle Stokes numbers as well as for the clean channel. The mode of instability which appears in all these cases is the traditional Tollmien–Schlichting (hereafter TS) instability, which is modified by the introduction of particles.

Figure 2. Validation for the case of uniform particle loading, in the form of neutral stability curves, with $\bar f = 0.05$, $S=5\times 10^{-5}$ and $S=2.5\times 10^{-4}$. Symbols correspond to Klinkenberg et al. (Reference Klinkenberg, De Lange and Brandt2011), while solid lines are from present computations. The region within the curves is unstable. The black dotted vertical line marks the minimum Reynolds number for $S=5 \times 10^{-5}$ at which instability is seen, termed the critical Reynolds number $R_{crit}$.

We are now in a position to study the instability mechanism. In the following section we derive a minimal equation set which allows us to highlight the basic physics.

3.1. Motivation

We saw in figure 2 that, in the case of constant particle loading, the TS mode of instability is modified by particles. Even with non-uniform particle loading, under non-overlap conditions, the same is observed. On the other hand, under overlap conditions, the picture is very different, and an example is shown in figure 4. Here, the TS mode is seen as a minor blip on the right of the figure. Two other modes of instability are now seen, which occur at much lower Reynolds number. The fact that these modes are distinct from the TS mode is evident from the separate regions in $\alpha \unicode{x2013}R$ space they occupy. To distinguish between them, the two lower-Reynolds-number modes of instability will be termed short wave and long wave, respectively, while remembering that the so-called short-wave mode actually has perturbation wavelengths of $O(1)$, i.e. comparable to the channel width (a wavelength of $\alpha =1$ is $2{\rm \pi}$ times the half-width). The long-wave modes extend from $O(1)$ to far lower wavenumbers. The short-wave mode of overlap instability occurs over the smallest Reynolds numbers, ranging from a few hundreds to a few thousands, while the long-wave mode spans decades in the Reynolds number, with the instability Reynolds number and the typical wavelength increasing together. In the following section we show that both of these are overlap modes of instability, caused by the variation of particle concentration within the critical layer. In § 3.3 we perform a similar analysis for the wall layer.

Figure 4. The three distinct modes of instability, shown by the shaded regions. A specific choice of parameters is made here, where overlap conditions prevail: the peak of the mean particle concentration profile has an amplitude $f_{max}=0.70$, and is positioned at $a_p=0.75$. The thickness of the particle layer is $\sigma =0.1$ and the Stokes number is $S=8\times 10^{-4}$. This figure is representative of a wide range of parameters under overlap conditions. The points marked $S$, $L$ and $T$ are representative of short-wave, long-wave and TS modes, respectively, and will be elaborated on.

In explaining the mechanism for the low-Reynolds-number instabilities, we may pursue one of two approaches. For both of them, we must begin by deriving the dominant balance in the critical layer. Once we have the lowest-order equations in the critical layer, we could solve the equations in the inner (critical) layer, and perform a matching with the outer layer (inviscid) solutions to obtain the full solutions. But this would yield no extra information, since we can already solve the full solutions. We therefore follow a second approach: of writing down a minimal composite theory for particulate shear flow. This theory (Narasimha & Govindarajan Reference Narasimha and Govindarajan2000; Govindarajan & Narasimha Reference Govindarajan and Narasimha2001; Bhattacharya et al. Reference Bhattacharya, Manoharan, Govindarajan and Narasimha2006) will obtain a reduced set of equations describing the stability problem. The reduced equations will contain all terms in the complete stability equations (3.23) which participate in the dominant balance somewhere in the flow, and none of the terms which do not participate in this anywhere.

3.2. Dominant balances in the critical layer

We first summarize existing knowledge in the context of a clean fluid, and then derive dominant balances within the critical layer in particulate shear flow.

In the Orr–Sommerfeld equation (2.26) for a clean fluid, it is seen that the highest, i.e. fourth-order, derivative term in $y$ is scaled by the inverse of the Reynolds number. Now, even if the Reynolds number approaches infinity, this term may not be dropped, because if it is, we will not be able to satisfy all four boundary conditions associated with (2.26). This is thus a classical singular perturbation problem (Van Dyke Reference Van Dyke1964), where the highest derivative term becomes as big as the terms on the left-hand side in some portions of the flow. There are two layers (Lin Reference Lin1945a,Reference Linb, Reference Lin1946) where viscous effects are important and gradients are large: the wall layer, of thickness $\epsilon _w \sim R^{-1/2}$, and the critical layer, of thickness $\epsilon \sim R^{-1/3}$ where, as defined above, $U \sim c$. It is the latter which is of primary interest to us to explain the mechanism of the overlap instabilities. To perform a similar analysis for particulate shear flow, we limit ourselves here to the regime where $RS \sim O(1)$, which is reasonable for dilute particle suspensions at high Reynolds number. A similar analysis may be carried out for any order of magnitude of this quantity. There are three layers we pay attention to on each side of the centreline, and these are shown in figure 3 for one half of the channel. There is also a wall layer shown, which will be discussed separately in § 3.3. There are now two critical layers: for the fluid and for the particle flow, of thickness $\epsilon$ and $\delta$, respectively, and the layer where the particles are concentrated, of thickness $\sigma$, which is pre-specified. The scales $\epsilon \ll 1$ and $\delta \ll 1$ are as yet unknown, and will be determined below. Figure 3(a) is a schematic for conditions where the particle layer and the critical layer are distinct, which we shall refer to as the non-overlap condition, and figure 3(b) depicts the overlap condition.

We derive equations within the critical (inner) layer in the inner variables $\xi$ and $\lambda$, defined as

(3.1a,b) \begin{equation} \xi=\frac{y-y_c}{\epsilon}\quad {\rm and} \quad \lambda=\frac{y-y_c}{\delta} , \end{equation}

and will select $\epsilon$ and $\delta$ to ensure that the derivatives of the fluid velocity components in $\xi$ and the particle velocity components in $\lambda$ are $O(1)$. In addition, it is useful to define

(3.2)\begin{equation} \chi=\frac{y-a_p}{\sigma}. \end{equation}

To derive the dominant balances we write the relevant variables in the form of series expansions within the critical layer as

(3.3a––c)\begin{equation} u_y=\sum_{n=0}^{\infty} \epsilon^n u_{y, n}(\xi), \quad v_y =\sum_{n=0}^{\infty} \delta^n v_{y,n}(\lambda) \quad {\rm and} \quad v_x =\sum_{n=0}^{\infty} \delta^n v_{x, n}(\lambda).\end{equation}

In this layer the mean flow may be written in the following expansion:

(3.4)\begin{equation} U(y)-c=(y-y_c)U_{c}^{'}+\frac{(y-y_c)^2}{2\!} U_{c}^{''}+\cdots.\end{equation}

The relative magnitudes within the critical layer of the two components of flow can be established from the continuity equation. We have

(3.5)\begin{equation} \sum_{n=0}^{\infty}\epsilon^n \left[{\rm i} \alpha u_{x,n} + \frac{1}{\epsilon} D_\chi u_{y,n}\right] = 0. \end{equation}

Constructing hierarchies of equations of different powers of $\epsilon$ yields

(3.6)\begin{equation} u_{y,0}=0,\end{equation}

and shows that the coefficient of a particular power of $\epsilon$ in $u_x$ is related to that which is one order higher in $u_y$. This is in fact a natural consequence of incompressibility. For the particle field, from (2.23), and using (3.3a–c) and (3.4), we get

(3.7)\begin{equation} \left[ U-c -\frac{i}{\alpha SR}\right]v_x={\rm i}\frac{U^{\prime}}{\alpha}v_y+\frac{1}{\alpha^2 SR}Du_{y} . \end{equation}

At the next two orders, using (3.7) as well as the incompressibility condition $D_\xi u_{y, 1} =- {\rm i} \alpha u_{x,0}$ in the critical layer, we get

(3.8a,b)\begin{equation} v_{x, 0} = u_{x, 0}, \quad v_{x,1}=\frac{\epsilon}{\delta} \left\{\frac{i}{\alpha}D_{\xi}u_{y,2} -{\rm i}\alpha SRU_c^{'}\xi v_{x,0} \right\}-SRU_{c}^{'}v_{y,1}. \end{equation}

This yields $\delta \sim \epsilon$, and without loss of generality, we choose $\delta = \epsilon$. The critical-layer thickness as perceived by the fluid and the particles is thus identical.

Using the third row of the matrix equation (2.25) along with (3.3a–c) and (3.4), and collecting terms at the lowest order in the expansion, we obtain

(3.9)\begin{equation} v_{y, 0} = u_{y, 0}, \end{equation}

which we know to be $0$ from (3.6). In other words, the expansions for the normal velocity components for the particles also begin one order higher than the streamwise component. At the next two orders, from (2.23)

(3.10)\begin{equation} \left[ U-c -\frac{i}{\alpha SR}\right]v_y=-\frac{i}{\alpha SR}u_{y},\end{equation}

we get

(3.11a,b)\begin{equation} v_{y,1}= u_{y,1} \quad {\rm and} \quad v_{y,2}=u_{y,2} -{\rm i}\alpha SR U_c^{'}\xi v_{y,1}. \end{equation}

From (3.3a–c), (3.8a,b) and (3.11a,b), we can see that the components of $\boldsymbol {v}$ and $\boldsymbol {u}$ differ from each other only at order $\epsilon$ relative to their largest value in the critical layer. This is consistent with the expectation that the particle velocity field must closely follow the fluid velocity field for low Stokes numbers. This analysis yields a measure of the difference between the two.

Finally, we may derive the dominant balance for fluid velocity in the critical layer from (2.23), along with (3.3a–c) and (3.4), and the Taylor expansion

(3.12)\begin{equation} \bar{f}= \begin{cases} \bar{f}_{c}+(y-y_{c})\bar{f}_{c}^{\prime}+ \dfrac{(y-y_c)^2}{2}\bar{f}_{c}^{\prime\prime}+\cdots & \text{overlap case} \\ \bar{f}_c & \text{non-overlap case} \end{cases}, \end{equation}

and under overlap conditions we may rewrite this in the relevant variable $\chi$. Now, there are different choices possible for the small parameters. We therefore a priori retain all terms which may participate in the dominant balance, and after some algebra obtain the following composite lowest-order equation:

(3.13)\begin{equation} \left[(1+\bar{f}) \xi U_{c}^{\prime}D_{\xi}^{2}+{\rm i}\frac{1}{\alpha R \epsilon^3}D_{\xi}^{4}-\frac{\epsilon}{\sigma}U_{c}^{\prime}(D_{\chi}\bar{f})(I-\xi D_{\xi})\right]u_{y,1}=0, \end{equation}

where $I$ is the identity operator. We will have one of the following four distinct cases arising. Case 1 includes the non-overlap case, while the others are for overlap conditions.

Case 1: either the particle layer and critical layer are well separated, or they overlap, but $1\sim {1}/{\alpha R\epsilon ^3}\gg {\epsilon }/{\sigma }$, i.e. the size of the particle layer significantly exceeds that of the critical layer. The third term in (3.13) now becomes negligible and in both cases this equation simplifies to

(3.14)\begin{equation} [-{\rm i}U_{c}^{\prime} \xi D_{\xi}^{2}+D_{\xi}^{4}]u_{y,1}=0. \end{equation}

The balance, with $\epsilon ={\alpha R}^{-1/3}$, is identical to the case with no particles, for the following reasons. When the particle layer and critical layer are well separated, the mass fraction $\bar {f}$ and its derivative ($D_\chi \bar {f}$) / $\sigma$ are practically negligible or much less than $O(1)$ in the critical layer, regardless of how small the particle layer $\sigma$ is.

Case 2: $1\ll {1}/{\alpha R\epsilon ^3}\sim {\epsilon }/{\sigma }$, i.e. the particle-laden layer is markedly smaller than the critical layer. The scaling that emerges is $\epsilon \sim ({\sigma }/{\alpha R})^{1/4}$ and, upon replacing $\sim$ by equality, (3.13) simplifies to

(3.15)\begin{equation} [ D_{\xi}^{4}+{\rm i}\{(D_{\chi}\bar{f})_{c} U_{c}^{\prime}\}(I-\xi D_{\xi})]u_{y,1}=0. \end{equation}

The variation in particle concentration is as important as the largest viscous effects in the critical layer.

Case 3: $1\sim {\epsilon }/{\sigma } \gg {1}/{\alpha R\epsilon ^3}$. The size of the particle layer is comparable to that of the critical layer, and the critical layer is much wider than that dictated by the scaling on the inverse Reynolds number. This case is a mathematical possibility, but is unlikely to occur physically, since the critical layer at large Reynolds will normally be influenced by the Reynolds number, and become thinner as Reynolds number increases. Under these hypothetical conditions, viscous effects appear only at higher order, since (3.13) simplifies to

(3.16)\begin{equation} [(1+\bar{f}) \xi U_{c}^{\prime}D_{\xi}^{2}-(D_{\chi}\bar{f}) U_{c}^{\prime}(I-\xi D_{\xi})]u_{y,1}=0,\end{equation}

where we have set $\epsilon =\sigma$.

Case 4: $1\sim {\epsilon }/{\sigma } \sim {1}/{\alpha R\epsilon ^3}$. The size of the particle layer is comparable to that of the critical layer, and viscosity plays a significant role as well. Equation (3.13) now becomes

(3.17)\begin{equation} [U_{c}^{\prime}(1+\bar{f}) \xi D_{\xi}^{2}+{\rm i} D_{\xi}^{4}-\kappa(D_{\chi}\bar{f}) U_{c}^{\prime}(I-\xi D_{\xi})]u_{y,1}=0.\end{equation}

Here, we define $\epsilon =(\alpha R)^{-1/3}$ and ${\epsilon }/{\sigma }\equiv \kappa$.

Equation (3.17) completely describes the critical layer for the thickness $\sigma$ we consider. Cases 2 to 4 correspond to overlap conditions, and the variation of the particle concentration within the critical layer, estimated by $D_{\chi }\bar {f}$, is an important player in the critical-layer balance, altering it fundamentally. We have thus established that the concentration profile will alter the fundamental nature of shear flow instability only under overlap conditions. This effect would be absent with uniform particle loading, where only under case 4 will we have a factor $(1+\bar {f})$ which merely rescales the Reynolds number.

3.3. Dominant balance in the wall layer

In shear flows, the critical and wall layers are often not well separated, and the overlap mode of instability presents such a case. We conduct the exercise below only to confirm that wall effects are not bringing in new physics into the instability. As before we define inner variables

(3.18a,b)\begin{equation} \xi=\frac{y-y_w}{\epsilon_w}\quad {\rm and} \quad \lambda=\frac{y-y_w}{\delta_w}, \end{equation}

where $y_w=\pm 1$ are the wall locations. Also, we have $U_w\equiv U(y_w)=0$ and $U_{w}^{'}\equiv U^{\prime }(y_w)$. We expand the variables in the form

(3.19a––c)\begin{equation} u_y=\sum_{n=0}^{\infty} \epsilon_{w}^n u_{y, n}(\xi), \quad v_y =\sum_{n=0}^{\infty} \delta_{w}^n v_{y,n}(\lambda) \quad {\rm and} \quad v_x =\sum_{n=0}^{\infty} \delta_{w}^n v_{x, n}(\lambda). \end{equation}

At the lowest order, $u_{y,0}=0$, and $v_{y,0}$ is proportional to this quantity and thus vanishes. At the next order, after some algebra, we obtain the scaling $\delta _w=\epsilon _w$, and we can use the incompressibility condition $D_\xi u_{y, 1} =- {\rm i} \alpha u_{x,0}$. Additionally, we have the following equations:

(3.20)\begin{gather} v_{x,0}=\dfrac{u_{x,0}}{1-{\rm i}\alpha cSR}, \end{gather}

(3.21a,b)\begin{gather}v_{x,1}=\frac{u_{x,1}}{1-{\rm i}\alpha cSR}-\left(\frac{SRU_{w}^{\prime}}{1-{\rm i}\alpha cSR}\right)v_{y,1} \quad {\rm and} \quad v_{y,1}=\frac{u_{y,1}}{1-{\rm i}\alpha cSR}. \end{gather}

We can express the dominant-balance composite equation for the flow as follows:

(3.22)\begin{equation} \left[-c\left(1+ \frac{\bar{f}}{1-{\rm i}\alpha cSR}\right) D_{\xi}^{2} + i\frac{1}{\alpha R\epsilon_w^2}D_{\xi}^4 - \frac{\epsilon_w}{\sigma}\left(\frac{cD_{\chi}\bar{f}}{1-{\rm i}\alpha cSR}\right)D_{\xi}\right]u_{y,1} = 0.\end{equation}

The structure of (3.22) in the wall layer is the same as (3.13) in the critical layer. There are changes in the coefficients, which change the scaling of $\epsilon _w$ to be $(\alpha R)^{-1/2}$. Again, when $\epsilon _w \sim \sigma$, along with significant overlap of the wall layer and the particle layer, the first derivative of the particle concentration profile is among the biggest terms at the lowest order. In our range of study, the wall layer is always very thin and well separated from the critical layer. All three terms are important when $\epsilon _w \sim (\alpha R)^{-1/2}$. In our investigation, we always positioned the particle layer at a considerable distance from the wall layer. Consequently, the mass fraction $\bar {f}$ and its derivative become negligible in the wall layer, and the dominant balance (3.22) is unaffected by the presence of particles.

3.4. Construction of the minimal composite theory

In the rest of the channel, a particle-laden counterpart of the Rayleigh equation, where all viscous effects are neglected in the stability operator, is valid. We now construct a reduced equation which includes every term in the complete (2.25) from which any of the dominant terms in the three layers originates, and neglect all other terms.

The final minimal composite equation is

(3.23) \begin{equation} [(U_{*}-c)(D^2-\alpha^2)-U'']\psi-(J'\bar{f}')\psi+(J\bar{f}')D\psi=\frac{1}{{\rm i}\alpha R}D^4 \psi,\end{equation}

with, as before,

(3.24)\begin{equation} U_{*} \equiv U+J\bar{f} \quad J\equiv \frac{U-c}{1+{\rm i}\alpha (U-c)SR}. \end{equation}

The terms that are not included in the above equation compared with the full (2.25) are: $-J^{''}\bar {f}\psi -(1/{\rm i}\alpha R)(-2\alpha ^2D^2+\alpha ^4)\psi$, apart from all the viscosity stratification effects which vanish from the minimal physics in a dilute suspension. This is because the derivative of the mean viscosity is $O(f_{max}/[\gamma \sigma ])$, which is small for a dilute suspension. Moreover, since $J$ is a function of $y$, we see that (3.23) represents a significant reduction of the complete stability operator. It comprises the inviscid stability operator of Rayleigh and the highest-order derivative in the viscous operator. For a constant particle loading, the only effect due to particles would come from the modified effective mean-flow profile $U_*$. Besides these, terms appear that are due to variations in the particle concentration profile which are critical. Note importantly that the minimal equation does not reduce to the Orr–Sommerfeld equation in any limit.

It remains to be seen whether the essential physics is contained in the minimal composite equation. The best parameter to make this explicit is the location, $a_p$, of the maximum in particle concentration. In figure 5 we show how the critical Reynolds number (the lowest Reynolds number for instability), $R_{crit}$, changes with $a_p$. The purple line with particles represents the full solution to (2.25), whereas the black line is the solution of the minimal composite equation (3.23). As $a_p$ increases, the particle-laden layer shifts towards the wall. We see that, below $a_p \sim 0.6$, there is a continuous increase in the critical Reynolds number. But beyond this, we see a sudden and large drop in $R_{crit}$, going down to values less than half of that a clean channel ($R_{crit}=5772.2$). At large $a_p$, however, i.e. when the particle-laden layer is very close to the wall, the trend is reversed again and a large stabilization is seen. The complete trend is captured by the minimal composite equation, although, as is to be expected, the agreement with the full solution is only qualitative.

Figure 5. With a specified amplitude of $f_{max}=0.1$, a peak width of $\sigma =0.1$ for the mass-fraction profile and a Stokes number of $S=2.5\times 10^{-4}$, the purple curve and the black curve in (a) illustrate the critical Reynolds number as a function of the position of the peak, $a_p$ of the mass fraction profile for the full (2.25) and the minimal equation (3.23), respectively. The green, black and red bands in (b) illustrate the notional critical layer, the wall layer and the particle-laden layer, respectively, as they vary with $a_p$.

The sensitivity of the critical Reynolds number to the location of the particle-laden layer is now seen to have its root in the dynamics within the critical layer of the dominant disturbance. The critical layer is shown in the inset of figure 5, as a function of $a_0$. A notional thickness of $R^{-1/3}$ is shown in this sketch. Shown in the same figure is the linear movement of the particle-laden layer with $a_p$. The large changes in stability occur in the regime of $a_p$ when the two layers overlap. A major portion of the disturbance kinetic energy is known to be produced within the critical layer in clean channel flow. We shall see that this is true of particulate flow too. Note that the wall layer (also shown in the figure) is unimportant in the dominant balance.

3.5. Summary of instability features in the overlap and non-overlap contitions

Before we discuss the energy budget, we summarize some additional features of the modes of instability. The maximum, $f_{max}$, in the particle concentration and the thickness, $\sigma$, of the particle-laden layer, have a quantitative, rather than qualitative, effect. We fix $\sigma$ at $0.1$.

Figures 6(a) and 6(b) summarize the variation of the critical Reynolds number (given in colour with contour lines) with the Stokes number and the particle loading. We examine two situations: at $a_p=0.4$, where the overlap mechanism is not in operation, and at $a_p=0.75$, where it is. It is evident that both in quality and quantity the two situations are very different. Under non-overlap conditions, i.e. where the particle-laden layer lies in a different part of the channel from the critical layer (figure 6a), we see stabilization as particle loading is increased. The stabilization is enormous in some portions of the regime, with the critical Reynolds number being extremely sensitive to either the particle loading or the Stokes number or both. The effect is largest at moderate particle Stokes number and is non-monotonic in the Stokes number. We now turn to figure 6(b), where completely the opposite trend is seen in response to increase in loading. For small to moderate Stokes number, with increase in loading, the flow is highly destabilized, with a sharp drop in critical Reynolds number. At high Stokes number, we see a reduction in the effect of particle loading, and a stabilization this time. The reason for this non-monotonicity could be as follows. Holding fixed the mass fraction $f_{max}$ of particles, as we increase the Stokes number, we are increasing the size of individual particles (from (2.8) we see that the particle radius $r$ scales as the square root of the Stokes number) and therefore reducing the number of particles. Thus, beyond a certain Stokes number, the forcing of the fluid by the particles comes down, and so does the effect of particle loading. We can check the consistency of this argument by examining the effect of Stokes number and mass loading on stability while holding the number density $N$ in (2.5) constant. Representative lines of constant $N$ are shown in figure 6. Following these lines from low $S$ and $f_{max}$ upwards, we again see different trends under overlap and non-overlap conditions. When the critical and particle-laden layers are distinct (figure 6a), and the particle number density is low (the two blue lines to the right of the figure), under non-overlap conditions, the critical Reynolds number is practically insensitive to changes in Stokes number and loading, i.e. the lines of constant $N$ appear parallel to lines of constant $R_{crit}$. This indicates that, at small particle number density, the critical Reynolds number is practically a function of $N$. At higher number densities (the two blue lines to the left of the figure) we see a strong stabilizing effect with increase in $S$ while $N$ is held constant. This is in sharp contrast to the stability response under overlap conditions, where we see that the effect of increasing Stokes number (and at the same time, mass loading) and constant $N$ is monotonically and strongly destabilizing.

Figure 6. Phase plot of the critical Reynolds number, shown in colour, in (a) and (b) as a function of the Stokes number $S$ and the particle loading strength $f_{max}$, and in (c) as a function of the particle loading location $a_p$ and $f_{max}$. (a) A case where there is no overlap mechanism in operation, with $a_p=0.40$, and (b) where it is in operation, with $a_p=0.75$. Note the difference in the colour bars in the two figures. In both figures $\sigma =0.1$, and the blue lines represent curves of constant particle number density $N$. The value of $N$ decreases from left to right, with the non-dimensional quantity $9\sqrt {2}H^3N\sqrt {\rho _{f}/\rho _p}$ being $[1.00\times 10^7, 4.42\times 10^5, 1.40\times 10^4, \ {\rm and } \ 1.25\times 10^3]$ for plot (a), and $[1.00\times 10^7, 8.94\times 10^5, 1.12\times 10^5, \ {\rm and }\ 3.95\times 10^4]$ for plot (b). The other parameters in (c) are the same as in figure (5). The green curve serves as a visual guide for observing the sharp change in the critical Reynolds number due to a change in the mode of instability. Some jitter in the colours is visible, and is due to the interpolation of results on a finitely spaced grid.

Figure 6(c) shows the variation of the critical Reynolds number with $a_p$ and $f_{max}$. A sharp contrast between the high critical Reynolds numbers of the TS mode and the far lower ones of the overlap mode is evident. At higher levels of $f_{max}$ the critical layer moves away from the wall, and $R_{crit}$ attains values as low as $200$ in this range.

Figure 7 shows the dependence of the neutral boundaries on $f_{max}$ and $a_p$. A similar figure appears in Rudyak & Isakov (Reference Rudyak and Isakov1996) as well. The long-wave mode is absent for smaller particle loading. The long-wave mode is odd in $u_y$ whereas the short-wave and TS modes are even. The two even modes can go through merging bifurcations, as seen in figure 7(b) at $a_p \sim 0.695$. After merger it is not easy to distinguish the boundary between the TS and the short-wave modes. The entire range of $a_p$ shown in this figure is small, underlining the sensitivity of the stability boundaries to this parameter. The even and odd overlap modes do not merge, but instead show a region of intersection, while each mode retains its character. They may always be distinguished by stipulating for the desired centreline conditions in the numerics. The long-wave instability is particularly interesting because, in channel flows, the most unstable perturbations are widely believed to be those which are even, with a maximum at the centreline, in the normal perturbation $u_y$. This assumption is so widespread that instability computations are often performed in the half-channel by imposing this symmetry at the centreline. Note that we use the terminology ‘odd’ mode going by the normal perturbation $u_y$.

Figure 7. (a) Stability boundaries for different amplitude of particle loading, with $a_p=0.75$ and $S= 8\times 10^{-4}$. (b) Sensitive dependence of the stability boundaries on the location $a_p$ of the particle concentration peak. Here, $f_{max}=0.4$, $S=2.5\times 10^{-4}$ and $\sigma =0.1$. The short-wave overlap mode undergoes a merger with the TS mode just past $a_p=0.695$. Also, there is significant intersection between the long-wave and short-wave modes.

The very fact that the long-wave mode is odd and the short-wave even means that their eigenfunctions are completely different in character, as seen in figures 8–12. The eigenfunctions have been normalized to set the maximum value of the streamfunction $\psi$ to unity. The short-wave instability in figure 9 exhibits much stronger streamwise velocity fluid perturbations compared with the TS mode throughout most of the channel, except in a narrow region where the TS streamwise velocity perturbations are pronounced. Also the near-wall structure of the streamwise velocity presents a distinct, wider and more symmetric reverse arrowhead shape than the TS. The eigenfunctions in the particle perturbation velocity components are strikingly different in the short-wave and the TS modes. In the short-wave mode, the particle dynamics is seen to follow the dynamics of the fluid, especially as seen in the streamwise velocity components. The normal velocity component is more peaky for the particles and more rounded for the flow. In contrast, the particles in the TS mode show a thinner region of strong streamwise velocity, but their wall-normal velocity is everywhere weak. The fact that the relevant portions of the eigenfunction profiles are thicker for the short-wave than for the TS mode is a consequence of the lower Reynolds numbers in the former, and we expect this from our critical-layer analysis above. For the short-wave mode, we compare the eigenfunctions from the full solution in figure 9 with those from the minimal composite equation in figure 10. Although the arrowhead shape is now distorted, the overall similarity in the eigenfunction structure between the two is striking. This is strong visual evidence that the dominant physics is contained in the minimal composite theory.

Figure 8. Typical eigenfunctions of the TS mode in the $x$–$y$ plane. This mode depicts point ‘$T$’ in figure 4 at $R=14\,000$, $\alpha =1.0$. A streamwise extent of two wavelengths is shown here and for all following eigenfunctions: (a) $\hat u_{x}(x,y)$, (b) $\hat u_{y}(x,y)$, (c) $\hat v_{x}(x,y)$ and (d) $\hat v_{y} (x,y)$.

Figure 9. Typical eigenfunctions of the short-wave mode shown at point ‘$S$’ in figure 4, where ${R}=1000$, $\alpha =1.6$: (a) $\hat u_{x} (x,y)$, (b) $\hat u_{y} (x,y)$, (c) $\hat v_{x} (x,y)$ and (d) $\hat v_{y} (x,y)$.

Figure 10. Characteristic eigenfunctions of the short-wave mode, obtained from the minimal composite equation (3.23) using identical parameters to those employed for figure 9: (a) $\hat u_{x} (x,y)$, (b) $\hat u_{y} (x,y)$, (c) $\hat v_{x} (x,y)$ and (d) $\hat v_{y} (x,y)$.

The eigenstructure of the long-wave instability is shown in figure 11 from the full equation. Again, the eigenfunctions from minimal composite theory, given in figure 12, are strikingly similar to those of the complete solution. The energy budgets in the following section will be rendered very surprising, given how different the odd and even modes are in their eigenstructure.

Figure 11. Typical eigenfunctions of the long-wave mode, at point ‘L’ in figure 4, where ${R}=3000$, $\alpha =0.6$: (a) $\hat u_{x}(x,y)$, (b) $\hat u_{y}(x,y)$, (c) $\hat v_{x}(x,y)$ and (d) $\hat v_{y} (x,y)$.

Figure 12. Characteristic eigenfunctions of the long-wave mode obtained from the minimal composite equation (3.23) using identical parameters to those employed for figure 11: (a) $\hat u_{x}(x,y)$, (b) $\hat u_{y}(x,y)$, (c) $\hat v_{x}(x,y)$ and (d) $\hat v_{y} (x,y)$.