2. Theoretical model

Consider a concentrated suspension of rigid, non-colloidal, neutrally buoyant spheres flowing between two parallel, stationary plates under the influence of a pressure field $P^{\prime }$ (see figure 1). The suspending medium is Newtonian with a viscosity $\mu _0^{\prime }$. The Hele-Shaw geometry has a thickness of $2B^{\prime }$ in the depth direction (3), and a characteristic length scale of $L^{\prime }$ in the lateral directions, 1 and 2. Here $x_1^{\prime }$ and $x_2^{\prime }$ are the co-ordinates along the lateral directions, while $x_3^{\prime }$ is the co-ordinate in the depth direction, along which the equations and variables will be averaged. We denote the depth-averaged volume fraction by $\bar{\phi}$ and the depth-averaged velocity field by $\bar{{u_{i}}}^{\prime }$. Following the development of the suspension balance model in our previous work (Ramachandran & Leighton Reference Ramachandran and Leighton2008), and using $L^{\prime }$ as the length scale for the $x_1^{\prime }$ and $x_2^{\prime }$ directions, $B^{\prime }$ as the length scale for the $x_3^{\prime }$ direction, $U_c^{\prime }$ as the velocity scale and $\mu _0^{\prime } U_c^{\prime } L^{\prime }/B^{\prime 2}$ as the stress scale, the following dimensionless governing equations emerge:

(2.1a)\begin{gather} \boldsymbol{\nabla}_{2_j}u_j+\frac{\partial u_3}{\partial x_3}=0,\end{gather}

(2.1b)$$\begin{gather} \epsilon^2\boldsymbol{\nabla}_{2_j}\left[\mu_r\left(\boldsymbol{\nabla}_{2_j}u_i+\boldsymbol{\nabla}_{2_i}u_j\right)\right]+\epsilon^2\frac{\partial}{\partial x_3}\left(\mu_r\frac{\partial u_3}{\partial x_i}\right)+\frac{\partial}{\partial x_3}\left(\mu_r\frac{\partial u_i}{\partial x_3}\right)\nonumber\\ \hspace{-11pc}=\frac{\partial P}{\partial x_i}+\epsilon\boldsymbol{\nabla}_{2_j}\left(\alpha\tau {\mathsf{H}}_{ij}\right), \end{gather}$$

(2.1c)\begin{gather} \epsilon^4\boldsymbol{\nabla}_{2_j}\left[\mu_r\left(\frac{\partial u_3}{\partial x_j}+\mu_r\frac{\partial u_j}{\partial x_3}\right)\right] +\epsilon^2\frac{\partial}{\partial x_3}\left(2\mu_r\frac{\partial u_3}{\partial x_3}\right)=\frac{\partial P}{\partial x_3}+\epsilon\frac{\partial }{\partial x_3}\left(\alpha\tau\right), \end{gather}

(2.1d)\begin{gather} \chi \epsilon \left[\frac{\partial \phi}{\partial t}+\boldsymbol{\nabla}_{2_j}\left(u_j\phi\right)+\frac{\partial}{\partial x_3} \left(u_3\phi\right)\right] =\frac{2}{9}\,\epsilon^2\boldsymbol{\nabla}_{2_i}\left[f\boldsymbol{\nabla}_{2_j}\left(\alpha\tau {\mathsf{H}}_{ij}\right)\right] +\frac{2}{9}\,\frac{\partial}{\partial x_3}\left[f\frac{\partial }{\partial x_3}\left(\alpha\tau \right) \right]. \end{gather}

Equation (2.1a) is the continuity equation, (2.1b) represents the momentum equations in the lateral directions $(x_1\text { and } x_2)$, (2.1c) is the momentum equation in the shallow direction $(x_3)$ and (2.1d) is the particle conservation equation. Here, the fields $P$, $u_i$ and $\phi$ are the dimensionless pressure, velocity and volume fraction, respectively. The operator $\boldsymbol {\nabla }_2$ is the gradient operator defined in the 1–2 co-ordinates only. The subscripted indices $i$ and $j$ assume values of 1 and 2 only. The parameter $\epsilon$ is the aspect ratio,

(2.2)\begin{equation} \epsilon=\frac{B^{\prime }}{L^{\prime }}, \end{equation}

and is the small parameter that will be exploited in the asymptotic expansion to follow. The quantity $\alpha$ in (2.1) is the reduced normal stress (Zarraga et al. Reference Zarraga, Hill and Leighton2000), given by

(2.3)\begin{equation} \alpha(\phi)=2.17\phi^3\textrm{e}^{2.34\phi}\left(1-\frac{\phi}{\phi_m}\right)^{-0.01}. \end{equation}

By allowing the reduced normal stress to diverge at maximum packing fraction, we overcome the issue of singularity in the volume fraction at $x_3=0$ (e.g. see Miller & Morris Reference Miller and Morris2006; Ramachandran & Leighton Reference Ramachandran and Leighton2007b, Reference Ramachandran and Leighton2008). The results presented in this paper are relatively insensitive to values of the exponent less than $0.01$. In (2.1), $\tau$ is a local measure of the dimensionless shear stress (Ramachandran & Leighton Reference Ramachandran and Leighton2008), given by

(2.4)\begin{equation} \tau=\mu_r(\phi)\dot{\gamma}, \end{equation}

where the shear rate $\dot {\gamma }$ is defined as

(2.5)\begin{align} \dot{\gamma}=\sqrt{\epsilon^2\left(\boldsymbol{\nabla}_{2_j}u_i+\boldsymbol{\nabla}_{2_i}u_j\right)\boldsymbol{\nabla}_{2_j}u_i+\epsilon^4\left(\boldsymbol{\nabla}_{2_j}u_3+\frac{\partial u_j}{\partial x_3}\right)\boldsymbol{\nabla}_{2_j}u_3+ \left(\frac{\partial u_i}{\partial x_3}\right)^2+\epsilon^2\frac{\partial u_i}{\partial x_3}\boldsymbol{\nabla}_{2_i}u_3}. \end{align}

In (2.1), ${\mathsf{H}}_{ij}$ is the geometry tensor, defined as

(2.6)\begin{equation} {\mathsf{H}}_{ij}=\delta_{ij}-b\frac{m_im_j}{m_km_k}+d\frac{p_ip_j}{p_kp_k}, \end{equation}

where $m_i$ and $p_i$, which are the velocity and vorticity directions, are given by

(2.7a,b)\begin{equation} m_i=u_i \quad \mathrm{and} \quad p_i=\epsilon_{ij3}m_j. \end{equation}

Figure 1. Schematic of the flow geometry.

The boundary conditions for the velocity component in the depth direction, $u_3$, are

(2.8)\begin{equation} \left.u_3\right|_{x_3=0,1}=0, \end{equation}

due to the antisymmetry about the centerplane and the zero normal flow condition at the wall. The boundary conditions for $u_i$ are the symmetry of the in-plane components of the velocity field at the midplane,

(2.9)\begin{equation} \left.\frac{\partial u_i}{\partial x_3}\right|_{x_3=0}=0, \end{equation}

and the no-slip boundary condition at the wall,

(2.10)\begin{equation} \left.u_i\right|_{x_3=1}=0. \end{equation}

The boundary conditions for the particle concentration are the no-flux restriction at the wall, and symmetry at the midplane, i.e.

(2.11)\begin{equation} \left.\left[\chi \epsilon u_3\phi-\frac{2}{9}\,f(\phi)\frac{\partial }{\partial x_3}\left(\alpha\tau H_{33}\right)\right]\right|_{x_3=0,1}=0. \end{equation}

Equation (2.8) allows us to simplify the above equation to

(2.12)\begin{equation} \left.\left[\frac{2}{9}\,f(\phi)\frac{\partial }{\partial x_3}\left(\alpha\tau\right)\right]\right|_{x_3=0,1}. \end{equation}

3. Multiple time scale expansion

We follow the procedure in Ramachandran (Reference Ramachandran2013) and implement a multiple time scale expansion of the governing equations and boundary conditions (Pagitsas et al. Reference Pagitsas, Nadim and Brenner1986; Hinch Reference Hinch1991). To begin, we need to define a small parameter relevant to this problem. For this, consider the ratio of time scale for shear-induced migration over the depth and the convective time scale in the flow direction,

(3.1)\begin{equation} \frac{B^{\prime 2}/\left[(U_c^{\prime }/B^{\prime })a^{\prime 2}\right]}{L^{\prime }/U_c^{\prime }}=\chi \epsilon. \end{equation}

Here, the parameter $\chi$ is

(3.2)\begin{equation} \chi=\frac{B^{\prime 2}}{a^{\prime 2}}. \end{equation}

In the Taylor dispersion limit (Brenner & Edwards Reference Brenner and Edwards1993; Probstein Reference Probstein1994; Ramachandran Reference Ramachandran2013), $\chi \epsilon =B^{\prime 3}/a^{\prime 2}L^{\prime }$ is a small parameter, i.e. the length scale $L^{\prime }$ in the direction parallel to the Hele-Shaw plates needs to be significantly greater than the induction length scale $B^{\prime 3}/a^{\prime 2}$ (Chapman Reference Chapman1990; Nott & Brady Reference Nott and Brady1994; Hampton et al. Reference Hampton, Mammoli, Graham and Tetlow1997; Ramachandran & Leighton Reference Ramachandran and Leighton2007a; Lecampion & Garagash Reference Lecampion and Garagash2014).

In dimensionless terms the two time variables for the multiple time scale expansion are $t_f=t$ (the fast time scale) and $t_s=\chi \epsilon t$ (the slow time scale). The partial derivative $\partial /\partial t$ is thus transformed to

(3.3)\begin{equation} \partial/\partial t = \partial/\partial t_f+\chi \epsilon\partial/\partial t_s. \end{equation}

After implementing this transformation, the only governing equation in (2.1) that is modified is the particle balance equation (2.1d), and it becomes

(3.4)$$\begin{align} &\chi \epsilon \left[\frac{\partial \phi}{\partial t_f}+\chi \epsilon\frac{\partial \phi}{\partial t_s}+\boldsymbol{\nabla}_{2_j}\left(u_j\phi\right)+\frac{\partial}{\partial x_3} \left(u_3\phi\right)\right] \notag\\ &\quad =\frac{2}{9}\,\epsilon^2\boldsymbol{\nabla}_{2_i}\left[f\boldsymbol{\nabla}_{2_j}\left(\alpha\tau {\mathsf{H}}_{ij}\right)\right] +\frac{2}{9}\,\frac{\partial}{\partial x_3}\left[f\frac{\partial }{\partial x_3}\left(\alpha\tau \right) \right]. \end{align}$$

All the dependent variables in the problem are now functions of five independent variables: $t_f$, $t_s$, $x_1$, $x_2$ and $x_3$.

The dependent variables can now be written as the following regular perturbation expansions in the small parameter $\chi \epsilon$:

(3.5)\begin{equation} \left. \begin{gathered} \phi=\phi^{(0)}+(\chi \epsilon) \phi^{(1)}+(\chi \epsilon)^2\phi^{(2)}+ \cdots, \\ P=P^{(0)}+(\chi \epsilon) P^{(1)}+(\chi \epsilon)^2P^{(2)}+ \cdots,\\ u_i=u_i^{(0)}+(\chi \epsilon) u_i^{(1)}+(\chi \epsilon)^2u_i^{(2)}+ \cdots\,, \\ u_3=u_3^{(0)}+(\chi \epsilon) u_3^{(1)}+(\chi \epsilon)^2u_3^{(2)}+ \cdots. \end{gathered} \right\} \end{equation}

In the subsections that follow, we will implement the multiple time scale expansion procedure. As a note, the overline that will appear over several variables henceforth represents the operator for performing an average over the depth of the Hele-Shaw geometry, i.e.

(3.6)\begin{equation} \bar{\mathcal{F}}\left(x_1,x_2,t\right)=\int_0^1\mathcal{F}\left(x_1,x_2,x_3,t\right)\,\textrm{d}x_3. \end{equation}

Also, the expansion procedure produces a host of functions: $M$, $F$, $g$, $G$ and $Q_1$ through $Q_{12}$. These functions are summarized in Appendix B.

To conclude this subsection, we note that the depth averages of the volume fraction and the two-dimensional velocity field are given by the averages of their zeroth-order approximations, i.e.

(3.7a,b)\begin{equation} \bar{\phi}=\overline{\phi^{(0)}} \quad \text{and}\quad \overline{u_i}=\overline{u_i^{(0)}}, \end{equation}

implying that

(3.8a,b)\begin{equation} \overline{\phi^{(j)}}=0 \quad \text{and} \quad \overline{u_i^{(j)}}=0, \end{equation}

for $j\geq 1$.

3.1. The $O[(\chi \epsilon )^0]$ problem

To leading order in $\chi \epsilon$, the particle balance equation (3.4) is

(3.9)\begin{equation} \frac{2}{9}\,\frac{\partial}{\partial x_3}\left[f(\phi^{(0)})\frac{\partial }{\partial x_3}\left(\alpha\tau\right)^{(0)} \right]=0. \end{equation}

The no-flux boundary condition (2.12) at the walls at this order,

(3.10)\begin{equation} \left.\left[\frac{2}{9}\,f(\phi^{(0)})\frac{\partial }{\partial x_3}\left(\alpha\tau\right)^{(0)}\right]\right|_{x_3=0,1}, \end{equation}

combined with (3.9) yields the following expression for the particle pressure:

(3.11)\begin{equation} \left(\alpha\tau\right)^{(0)}=c'. \end{equation}

Here $c'$ is the constant of integration. To obtain the stress $\tau ^{(0)}$ required for computing the volume fraction profile, the continuity and momentum equations in (2.1a) to (2.1c) are examined at $O[(\chi \epsilon )^0]$ to yield

(3.12a)\begin{gather} \boldsymbol{\nabla}_{2_j}u_j^{(0)}+\frac{\partial u_3^{(0)}}{\partial x_3}=0, \end{gather}

(3.12b)\begin{gather} \frac{\partial}{\partial x_3}\left[\mu_r(\phi^{(0)})\frac{\partial u_i^{(0)}}{\partial x_3}\right]=\boldsymbol{\nabla}_{2_j} P^{(0)}, \end{gather}

(3.12c)\begin{gather} \frac{\partial P^{(0)}}{\partial x_3}=0. \end{gather}

As the pressure is invariant in the $x_3$ direction at this order, (3.12b) can be integrated with respect to $x_3$. Doing this once yields

(3.13)\begin{equation} \mu_r(\phi^{(0)})\frac{\partial u_i^{(0)}}{\partial x_3}=\boldsymbol{\nabla}_{2_j} P^{(0)}x_3, \end{equation}

where the symmetry condition (2.9) has been applied. The stress, $\tau$, at this order, defined using (2.4) and (2.5), can be deduced from the above result as

(3.14)\begin{equation} \tau^{(0)}=\left|{\boldsymbol{\boldsymbol{\nabla}}}_2 P^{(0)}\right|x_3. \end{equation}

Note that $\left |{\boldsymbol {\boldsymbol {\nabla }}}_2 P^{(0)}\right |$ is independent of $x_3$ due to (3.12c). Substitution of the above result in (3.11) yields

(3.15)\begin{equation} \alpha(\phi^{(0)})x_3=\alpha_w, \end{equation}

where $\alpha _w$ is the reduced normal stress evaluated at the wall. The above equation provides the volume fraction distribution for a given depth-averaged volume fraction $\bar {\phi }$,

(3.16)\begin{equation} \phi^{(0)}(\bar{\phi},x_3)=\alpha^{{-}1}\left[\frac{\alpha_w}{x_3}\right], \end{equation}

In figure 2(a) we have shown the volume fraction profile over the depth for different $\bar {\phi }$. Similar to tube flow of suspensions (Phillips et al. Reference Phillips, Armstrong, Brown, Graham and Abbott1992; Hampton et al. Reference Hampton, Mammoli, Graham and Tetlow1997; Fang et al. Reference Fang, Mammoli, Brady, Ingber, Mondy and Graham2002; Miller & Morris Reference Miller and Morris2006; Ramachandran & Leighton Reference Ramachandran and Leighton2008; Ramachandran Reference Ramachandran2013), there is a region of maximum packing fraction around the centreline due to the phenomenon of shear-induced migration (Koh, Hookham & Leal Reference Koh, Hookham and Leal1994; Lyon & Leal Reference Lyon and Leal1998; Fang et al. Reference Fang, Mammoli, Brady, Ingber, Mondy and Graham2002; Gao & Gilchrist Reference Gao and Gilchrist2008; Lecampion & Garagash Reference Lecampion and Garagash2014; Yadav, Reddy & Singh Reference Yadav, Reddy and Singh2016; Dontsov et al. Reference Dontsov, Boronin, Osiptsov and Derbyshev2019), which becomes thicker as the volume fraction increases.

Figure 2. Predictions of the zeroth-order profiles of (a) volume fraction ($\phi ^{(0)}$, (3.16)), (b) the function quantifying the in-plane velocity ($F$, (3.17)) and (c) the velocity along the thin direction ($\hat {Q}_1=Q_1//M$, (3.26)), as functions of the co-ordinate in the thin direction, $x_3$. The solid, dotted and dash–dotted curves are the results for depth-averaged volume fractions of 20 %, 30 % and 45 %. The inset in subfigure (a) provides a zoomed-in view of the volume fraction profiles near the centreline.

To get the velocity distribution, (3.13) may be divided by $\mu _r(\phi ^{(0)})$ and integrated in $x_3$ to produce the following expression for $u_i$ after the application of the no-slip boundary condition (2.10):

(3.17)\begin{equation} u_i^{(0)}={-}F(\bar{\phi},x_3)\boldsymbol{\nabla}_{2_i} P^{(0)}. \end{equation}

Here, the function $F$ is defined as

(3.18)\begin{equation} F(\bar{\phi},x_3)=\int_{x_3}^1\frac{\xi \,\textrm{d}\xi}{\mu_r[\phi^{(0)}(\bar{\phi},\xi)]}. \end{equation}

The velocity distribution is shown as a function of $x_3$ for different depth-averaged volume fractions in figure 2(b). Again, analogous to tube flow, the maximum packing region near the centreline leads to blunting of the velocity profile due to the high viscosity associated with this region (Koh et al. Reference Koh, Hookham and Leal1994; Lyon & Leal Reference Lyon and Leal1998; Fang et al. Reference Fang, Mammoli, Brady, Ingber, Mondy and Graham2002; Gao & Gilchrist Reference Gao and Gilchrist2008; Lecampion & Garagash Reference Lecampion and Garagash2014; Yadav et al. Reference Yadav, Reddy and Singh2016; Dontsov et al. Reference Dontsov, Boronin, Osiptsov and Derbyshev2019). As $\bar {\phi }$ increases, the average velocity for a given pressure gradient decreases, and the size of the blunted zone increases. Notably, the sizes of the maximum packing regions and the blunted zones are smaller than the corresponding ones for tube flow. This is because the area of a region of a given width near the centre represents a smaller fraction of the total area for tube Poiseuille flow than plane Poiseuille flow.

The depth-averaged velocity at this order is

(3.19)\begin{equation} \overline{u_i^{(0)}}=\overline{u_i}={-}M(\bar{\phi})\boldsymbol{\nabla}_{2_i} P^{(0)}. \end{equation}

Here, the function $M$ quantifies the mobility of the suspension, and is defined as

(3.20)\begin{equation} M(\bar{\phi})=\int_0^1 F(\bar{\phi},x_3)\,\textrm{d}x_3. \end{equation}

Integrating the continuity equation (3.12a), and applying the boundary conditions for $u_3$ (see (2.8)), we get the governing equation for the pressure distribution at this order,

(3.21)\begin{equation} \boldsymbol{\nabla}_{2_j}[M(\bar{\phi})\boldsymbol{\nabla}_{2_j} P^{(0)}]=0. \end{equation}

Finally, we derive an expression for $u_3^{(0)}$, which will be required in the $O(\chi \epsilon )$ analysis. Using the continuity equation (3.12a), we get

(3.22)\begin{equation} \frac{\partial u_3^{(0)}}{\partial x_3}=\boldsymbol{\nabla}_{2_j} [F(\bar{\phi},x_3)\boldsymbol {\nabla}_{2_j} P^{(0)}]. \end{equation}

Integrating the above equation with respect to $x_3$, we get

(3.23)\begin{equation} u_3^{(0)}=\boldsymbol{\nabla}_{2_j} [Q_1(\bar{\phi},x_3) \boldsymbol{\nabla}_{2_j} P^{(0)}], \end{equation}

where

(3.24)\begin{equation} Q_1(\bar{\phi},x_3)=\int_0^{x_3}F(\bar{\phi},\xi)\,\textrm{d}\xi. \end{equation}

Using the following modified version of the averaged continuity equation (3.21),

(3.25)\begin{equation} \boldsymbol{\nabla}_2^2 P^{(0)}={-}\frac{M^{\prime }(\bar{\phi})}{M(\bar{\phi})}\boldsymbol{\nabla}_{2_j}\bar{\phi}\boldsymbol{\nabla}_{2_j} P^{(0)}, \end{equation}

and the relationship between $\overline {u_j}$ and $\boldsymbol {\nabla }_2 P^{(0)}$ in (3.19), (3.23) can be simplified to

(3.26)\begin{equation} u_3^{(0)}=\hat{Q}_1(\bar{\phi},x_3)\boldsymbol{\nabla}_{2_j}\bar{\phi} \quad \overline{u_j}, \end{equation}

where

(3.27)\begin{equation} \hat{Q}_1(\bar{\phi},x_3)=\frac{M^{\prime }(\bar{\phi})}{M^2(\bar{\phi})}Q_1(\bar{\phi},x_3)-\frac{1}{M(\bar{\phi})}\frac{\partial Q_1}{\partial \bar{\phi}}(\bar{\phi},x_3). \end{equation}

Shown in figure 2(c) is the function $\hat {Q}_1$ as a function of $x_3$ for different $\bar {\phi }$. For positive gradients in the volume fraction along the flow direction $\left (\boldsymbol {\nabla }_{2_j}\bar {\phi }\ \overline {u_j}>0\right )$, the suspension viscosity increases in the flow direction, and this increase is stronger near the centreplane of the Hele-Shaw geometry than the walls. As a result, there is a decrease in the planar velocity field $u_i^{(0)}$ in the flow direction, which, by conservation of mass, leads to an expulsion of suspension from the centreplane to the walls, i.e. a positive $u_3^{(0)}$. Similarly, negative gradients in the volume fraction $(\boldsymbol {\nabla }_{2_j}\bar {\phi }\ \overline {u_j}<0)$ along the flow direction lead to negative $u_3^{(0)}$.

3.2. The $O[(\chi \epsilon )^1]$ problem

The particle balance equation (3.4) at $O[(\chi \epsilon )^1]$ is

(3.28)$$\begin{align} &\frac{\partial \phi^{(0)}}{\partial t_f}-\nabla_{2_j}\left[F\left(\overline{\phi},x_3\right)\phi^{(0)}\left(\overline{\phi},x_3\right)\nabla_{2_j} P^{(0)}\right]\notag\\ &\quad+ \frac{\partial}{\partial x_3} \left[u_3^{(0)}\phi^{(0)}\left(\overline{\phi},x_3\right)\right]=\frac{2}{9}\frac{\partial}{\partial x_3}\left[f\left(\phi\right)\frac{\partial }{\partial x_3}\left(\alpha\tau \right) \right]^{(1)}.\end{align}$$

In the above equation, we have substituted the expression for $u_i^{(0)}$ from (3.17). Integrating the above equation with respect to $x_3$ from 0 to 1, we get

(3.29)\begin{equation} \frac{\partial \bar{\phi}}{\partial t_f}-\boldsymbol{\nabla}_{2_j} [G(\bar{\phi}) \boldsymbol{\nabla}_{2_j} P^{(0)}]=0, \end{equation}

where

(3.30)\begin{equation} G(\bar{\phi})=\int_0^1F(\bar{\phi},\xi)\phi^{(0)}(\bar{\phi},\xi)\,\textrm{d}\xi. \end{equation}

In (3.29) we have applied the no-flux boundary condition in (2.12) at this order, and also used the result for the zeroth-order particle pressure from (3.11). Using the following result from the chain rule,

(3.31)\begin{equation} \frac{\partial \phi^{(0)}}{\partial t_f}=\frac{\partial \phi^{(0)}}{\partial \bar{\phi}}\frac{\partial \bar{\phi}}{\partial t_f}, \end{equation}

we can manipulate (3.28) and (3.29) to arrive at

(3.32) $$\begin{align} &\frac{\partial \phi^{(0)}}{\partial \overline{\phi}}\nabla_{2_j}\left[G\left(\overline{\phi}\right)\nabla_{2_j} P^{(0)}\right]-\nabla_{2_j}\left[F\left(\overline{\phi},x_{3}\right)\nabla_{2_j} P^{(0)}\phi^{(0)}\right]\notag\\ &\quad +\frac{\partial}{\partial x_{3}} \left(u_{3}^{(0)}\phi^{(0)}\right)=\frac{2}{9}\frac{\partial}{\partial x_3}\left[f\left(\phi^{(0)}\right)\frac{\partial }{\partial x_3}\left(\alpha\tau \right)^{(1)} \right]. \end{align}$$

To arrive at the right-hand side of the above equation from (3.28), we have used the constancy of $(\alpha \tau )^{(0)}$ over the depth (3.11). After the integration of the above equation in $x_3$ and the implementation of the algebraic steps discussed in Appendix C, we get the volume fraction perturbation $\phi ^{(1)}$ as

(3.33)\begin{equation} \phi^{(1)} =Q_7(\bar{\phi},x_3)\frac{\boldsymbol{\nabla}_{2_j}\bar{\phi} (-\boldsymbol{\nabla}_{2_j} P^{(0)})}{|{\boldsymbol{\boldsymbol{\nabla}}}_{2}P^{(0)}|}=Q_7(\bar{\phi},x_3)\left(\frac{\boldsymbol{\nabla}_{2_j}\bar{\phi}\hspace{0.1in}\overline{u_j}}{\left|\overline{\boldsymbol{u}}\right|}\right). \end{equation}

Note that $\phi ^{(1)}$ is proportional to the directional derivative of $\bar {\phi }$ along the depth-averaged velocity vector.

To derive the velocity perturbation, we begin with (2.1b) at $O[(\chi \epsilon )^1]$, which gives

(3.34)\begin{equation} \frac{\partial}{\partial x_3}\left(\mu_r\frac{\partial u_i}{\partial x_3}\right)^{(1)}=\boldsymbol{\nabla}_{2_i} P^{(1)}. \end{equation}

Integration of the above equation, followed by the application of the boundary condition at the wall and some simplifications (see Appendix D for details), yields

(3.35)\begin{equation} u_i^{(1)}=Q_{10}(\bar{\phi},x_3)\frac{\boldsymbol{\nabla}_{2_j}\bar{\phi}\boldsymbol{\nabla}_{2_j} P^{(0)}}{\left|\boldsymbol{\boldsymbol{\nabla}}_{2} P^{(0)}\right|}\boldsymbol{\nabla}_{2_i} P^{(0)}=\frac{Q_{10}(\bar{\phi},x_3)}{M(\bar{\phi})}\left(\frac{\boldsymbol{\nabla}_{2_j}\bar{\phi}\hspace{0.1in}\overline{u_j}}{\left|\overline{\boldsymbol{u}}\right|}\right)\overline{u_i}. \end{equation}

In figure 3 we have shown the functions $Q_7$ and $Q_{10}/M$ corresponding to the volume fraction and velocity distributions at $O\left (\chi \epsilon \right )$ for $\bar {\phi }=0.2, 0.3$ and 0.45. When $\boldsymbol {\nabla }_{2_j}\bar {\phi }\ \overline {u_j}>0$, the volume fraction perturbation is negative near the centre (see figure 3a). As explained in Ramachandran (Reference Ramachandran2013), this is due to two reasons: the displacement of a low volume fraction suspension into a region of high volume fraction, and the positive $u_3^{(0)}$ velocity which drives particles away from the centre and towards the walls when $\left (\boldsymbol {\nabla }_{2_j}\bar {\phi }\ \overline {u_j}>0\right )$. The greater the volume fraction, the smaller the perturbation due to the smaller contrast in the volume fraction between the centre and wall, and the weaker the magnitude of $u_3^{(0)}$. The negative value of $\phi ^{(1)}$ near the centre leads to a decrease in suspension viscosity there and, hence, a positive perturbation in the first-order velocity field (see figure 3b), the effect being stronger for higher volume fractions. Similar arguments can be used to deduce trends when $\boldsymbol {\nabla }_{2_j}\bar {\phi }\ \overline {u_j}<0$.

Figure 3. Predictions of the first-order profiles of (a) volume fraction ($Q_7$, (3.33)) and (b) the function quantifying the in-plane velocity, ($Q_{10}/M$, (3.35)), as functions of the co-ordinate in the thin direction, $x_3$. The solid, dotted and dash–dotted curves are the results for depth-averaged volume fractions of 20 %, 30 % and 45 %.

We end this subsection by commenting that it is possible to determine the pressure $P^{(1)}$ and the velocity component $u_3^{(1)}$ by employing the continuity equation, but these variables are not required to derive the macrotransport equation to be presented at the end of this section. Hence, the calculations of $P^{(1)}$ and $u_3^{(1)}$ will not be pursued.

3.3. The $O[(\chi \epsilon )^2]$ problem and the averaged equations

The particle balance equation (3.4) at $O[(\chi \epsilon )^2]$ is

(3.36)$$\begin{align} &\frac{\partial \phi^{(1)}}{\partial t_f}+\boldsymbol{\nabla}_{2_i}\left(u_i^{(0)}\phi^{(1)}+u_i^{(1)}\phi^{(0)}\right)+\frac{\partial}{\partial x_3} \left(u_3\phi\right)^{(1)}+\frac{\partial \phi^{(0)}}{\partial t_s}\nonumber\\ &\quad =\frac{\epsilon^2}{(\chi\epsilon)^2}\frac{2}{9}\,\boldsymbol{\nabla}_{2_i}\left[f\boldsymbol{\nabla}_{2_j}\left(\alpha\tau {\mathsf{H}}_{ij}\right)\right]^{(0)} +\frac{2}{9}\,\frac{\partial}{\partial x_3}\left[f\frac{\partial }{\partial x_3}\left(\alpha\tau \right) \right]^{(2)}. \end{align}$$

Note that we have preserved the leading order shear-induced migration term, $({2}/{9})\epsilon ^2\boldsymbol {\nabla }_{2_i}[f\boldsymbol {\nabla }_{2_j}(\alpha \tau {\mathsf{H}}_{ij})]^{(0)}$, even though it is of order $\chi ^{-2}$ relative to other terms at this order. As will become clear later, the key particle flux that emerges from the left-hand side of the equation, the Taylor dispersion flux, does not have a component in the direction normal to the streamlines defined by the velocity field $\overline {u_i}$ (see (3.33) and (3.35)). The leading order source of particle flux perpendicular to the flow streamline is the shear-induced migration flux and, hence, it has to be retained.

We can use (3.14), (3.15) and (3.19) to write

(3.37)\begin{equation} (\alpha\tau {\mathsf{H}}_{ij})^{(0)}=\frac{\alpha_w(\bar{\phi})}{M(\bar{\phi})}\left|\overline{\boldsymbol{u}}\right| {\mathsf{H}}_{ij}^{(0)}, \end{equation}

where

(3.38)\begin{equation} {\mathsf{H}}_{ij}^{(0)}=\delta_{ij}-b\frac{\overline{u_i}\ \overline{u_j}}{\overline{u_k}\ \overline{u_k}}+d\frac{\epsilon_{ik3}\overline{u_k}\ \epsilon_{jl3}\overline{u_l}}{\overline{u_p}\ \overline{u_p}}. \end{equation}

This modifies (3.36) to

(3.39)\begin{equation} \frac{\partial \bar{\phi}}{\partial t_s}+\boldsymbol{\nabla}_{2_i}\left(\overline{u_i^{(0)}\phi^{(1)}}+\overline{u_i^{(1)}\phi^{(0)}}\right)=\boldsymbol{\nabla}_{2_i}\left[\frac{2}{9}\,\overline{f^{(0)}}\boldsymbol{\nabla}_{2_j}\left(\frac{\alpha_w(\bar{\phi})}{M(\bar{\phi})}\left|\overline{\boldsymbol{u}}\right| {\mathsf{H}}_{ij}^{(0)}\right)\right]. \end{equation}

Let us now consider the two convective terms on the left-hand side of the above equation. The first one, which arises due to the volume fraction perturbation $\phi ^{(1)}$, can be written as

(3.40)\begin{equation} \overline{u_j^{(0)}\phi^{(1)}}={-}M(\bar{\phi})Q_{11}(\bar{\phi}) \left(\frac{\boldsymbol{\nabla}_{2_j}\bar{\phi}\ \overline{u_j}}{\left|\overline{\boldsymbol{u}}\right|}\right)\overline{u_i}, \end{equation}

where

(3.41)\begin{equation} Q_{11}(\bar{\phi})={-}\frac{1}{M(\bar{\phi})}\int_0^1F(\bar{\phi},\xi)Q_7(\bar{\phi},\xi)\,\textrm{d}\xi, \end{equation}

while the second one, which arises due to the perturbation, $u_i^{(1)}$, in the velocity field, simplifies to

(3.42)\begin{equation} \overline{u_j^{(1)}\phi^{(0)}}=M(\bar{\phi})Q_{12}(\bar{\phi})\left(\frac{\boldsymbol{\nabla}_{2_j}\bar{\phi}\hspace{0.1in}\overline{u_j}}{\left|\overline{\boldsymbol{u}}\right|}\right), \end{equation}

with

(3.43)\begin{equation} Q_{12}(\bar{\phi})=\frac{1}{M(\bar{\phi})}\int_0^1Q_{10}(\bar{\phi},\xi)\phi^{(0)}(\bar{\phi},\xi)\,\textrm{d}\xi. \end{equation}

A negative sign has been added to the definition of $Q_{11}$ in (3.41) to ensure that $Q_{11}$ is a positive function. Equation (3.39) thus becomes

(3.44)\begin{align} \frac{\partial \bar{\phi}}{\partial t_s}=\boldsymbol{\nabla}_{2_i}\left[M(\bar{\phi})h(\bar{\phi})\left(\frac{\boldsymbol{\nabla}_{2_j}\bar{\phi}\hspace{0.1in}\overline{u_j}}{\left|\overline{\boldsymbol{u}}\right|}\right)\overline{u_i}\right]+\frac{\epsilon^2}{(\chi\epsilon)^2}\boldsymbol{\nabla}_{2_i}\left[\frac{2}{9}\,\overline{f^{(0)}}\boldsymbol{\nabla}_{2_j}\left(\frac{\alpha_w(\bar{\phi})}{M(\bar{\phi})}\left|\overline{\boldsymbol{u}}\right| {\mathsf{H}}_{ij}^{(0)}\right)\right], \end{align}

where

(3.45)\begin{equation} h(\bar{\phi})=Q_{11}(\bar{\phi})-Q_{12}(\bar{\phi}). \end{equation}

After expanding out the shear-induced migration term in the above equation, we get

(3.46)$$\begin{align} \frac{\partial \bar{\phi}}{\partial t_s}&=\boldsymbol{\nabla}_{2_i}\left[h(\bar{\phi})\left(\frac{\boldsymbol{\nabla}_{2_j}\bar{\phi}\hspace{0.1in}\overline{u_j^{(0)}}}{\left|\overline{\boldsymbol{u}}\right|}\right)\overline{u_i}\right]+\frac{\epsilon^2}{(\chi\epsilon)^2}\boldsymbol{\nabla}_{2_i} (q_A(\bar{\phi})\left|\overline{\boldsymbol{u}}\right| {\mathsf{H}}_{ij}^{(0)}\boldsymbol{\nabla}_{2_j}\bar{\phi})\nonumber\\ &\quad +\frac{\epsilon^2}{(\chi\epsilon)^2}\boldsymbol{\nabla}_{2_i}\left[q_B(\bar{\phi})\boldsymbol{\nabla}_{2_j}\left(\left|\overline{\boldsymbol{u}}\right| {\mathsf{H}}_{ij}^{(0)}\right)\right], \end{align}$$

where

(3.47)\begin{equation} \left. \begin{gathered} q_A(\bar{\phi})=\frac{2}{9}\,\overline{f^{(0)}}\frac{\textrm{d}}{\textrm{d}\bar{\phi}}\left[\frac{\alpha_w(\bar{\phi})}{M(\bar{\phi})}\right],\\ q_B(\bar{\phi})=\frac{2}{9}\,\overline{f^{(0)}}\frac{\alpha_w(\bar{\phi})}{M(\bar{\phi})}. \end{gathered} \right\} \end{equation}

We now rewrite (3.29) with the pressure gradient $\boldsymbol {\nabla }_{2_j}P^{(0)}$ replaced by the depth-averaged velocity field $\overline {u_i}$ using (3.19),

(3.48)\begin{equation} \frac{\partial \bar{\phi}}{\partial t_f}+\boldsymbol{\nabla}_{2_i}\left[g(\bar{\phi})\overline{u_i}\right]=0, \end{equation}

where the function $g$ is the local flow-averaged volume fraction (Ramachandran & Leighton Reference Ramachandran and Leighton2007b; Ramachandran Reference Ramachandran2013), defined as

(3.49)\begin{equation} g(\bar{\phi})=\frac{G(\bar{\phi})}{M(\bar{\phi})}. \end{equation}

Multiplying (3.48) by $\chi \epsilon$ and (3.46) by $\left (\chi \epsilon \right )^2$, and then summing the resulting equations, yields

(3.50)$$\begin{align} \left(\chi\epsilon\right)\left[\frac{\partial \bar{\phi}}{\partial t}+\boldsymbol{\nabla}_{2_i}\left(g(\bar{\phi})\overline{u_i}\right)\right] &=\left(\chi\epsilon\right)^2\boldsymbol{\nabla}_{2_i}\left[h(\bar{\phi})\left(\frac{\boldsymbol{\nabla}_{2_j}\bar{\phi} u_j}{|\overline{\boldsymbol{u}}|}\right)\overline{u_i}\right]\nonumber\\ &\quad +\epsilon^2\boldsymbol{\nabla}_{2_i}(q_A(\bar{\phi})\left|\overline{\boldsymbol{u}}\right| {\mathsf{H}}_{ij}^{(0)}\boldsymbol{\nabla}_{2_j}\bar{\phi})\notag\\ &\quad +\epsilon^2\boldsymbol{\nabla}_{2_i} [q_B(\bar{\phi})\boldsymbol{\nabla}_{2_j} (|\overline{\boldsymbol{u}}| {\mathsf{H}}_{ij}^{(0)})]. \end{align}$$

Reverting to dimensional variables, we obtain

(3.51)$$\begin{align} \frac{\partial \bar{\phi}}{\partial t'}+\boldsymbol{\nabla}_{2_i}^{\prime }[g(\bar{\phi})\overline{u'_i}] &=\frac{B^{\prime 3}}{a^{\prime 2}}\boldsymbol{\nabla}_{2_i}^{\prime }\left[h(\bar{\phi})\left(\frac{\boldsymbol{\nabla}_{2_j}\bar{\phi} \ \overline{u_j^{\prime }}}{\left|\overline{\boldsymbol{u}}^{\prime }\right|}\right)\overline{u_i^{\prime }}\right]\nonumber\\ &\quad +\frac{a^{\prime 2}}{B^{\prime }}\boldsymbol{\nabla}_{2_i}^{\prime }(q_A(\bar{\phi})|\overline{\boldsymbol{u^{\prime }}}| {\mathsf{H}}_{ij}^{(0)}\nabla^{\prime }_{2_j}\bar{\phi})+\frac{a^{\prime 2}}{B^{\prime }}\boldsymbol{\nabla}_{2_i}^{\prime }[q_B(\bar{\phi})\boldsymbol{\nabla}_{2_j}^{\prime }(|\overline{\boldsymbol{u}^{\prime }}| {\mathsf{H}}_{ij}^{(0)})]. \end{align}$$

The above equation is the macrotransport equation for the depth-averaged volume fraction.

5. Viscous miscible fingering instability in Hele-Shaw suspension flows

As follows from § 4, we anticipate the case of a low viscosity suspension displacing a high viscosity suspension in a Hele-Shaw geometry to lead to a viscous miscible fingering instability. Here, we analyse the stability of a base state comprising a step increase in the volume fraction across $x_1=0$ along the flow direction (see schematic in figure 9). Although more sophisticated spectral approaches are possible (Ben et al. Reference Ben, Demekhin and Chang2002), we follow Tan & Homsy (Reference Tan and Homsy1986) by invoking the quasi-steady state approximation (QSSA) to deduce the relationship between the growth rate and wavenumber of a perturbation perpendicular to the flow in the $x_2$ direction. The QSSA assumes that the rate of change of the base state is small compared with the rate of growth of the perturbations. Thus, the base state is assumed to be ‘frozen’ in time, and the stability of perturbations about this state is examined. The advantage of this approach is that it permits a semi-analytical solution to the linear stability problem at $t=0$. Tan & Homsy (Reference Tan and Homsy1986) provide a comparison of growth rates from QSSA computations against rates from numerical solutions of the governing equations with a random noise as the initial condition. They show that except for a short initial period, there are negligible differences in the results of the growth rates from the two calculations, and that the QSSA treatment captures the essential features of the instability.

Figure 9. Definition sketch for the linear stability analysis of the viscous miscible fingering problem for the Hele-Shaw flow of suspensions. The initial configuration comprises a suspension at a low volume fraction $\bar {\phi }_-$ $(x_1<0)$ interfacing a suspension at a high volume fraction $\bar {\phi }_+$, $x_1>0$, in the form of a flat plane at $x_1=0$. The combination is displaced from left to right. Sinusoidal perturbations of the volume fraction, pressure and velocity fields are imposed at about the base configuration in the $x_2$ direction.

For the initial volume fraction distribution shown in figure 9, at time $t=0$, the base state solution profiles for the pressure, volume fraction and velocity obtained from (4.1) to (4.5) and written in the frame of reference of the shock velocity, $U_s^{\prime }$ (4.6) are

(5.1a–d)$$\begin{gather} \bar{\phi}^{(b)}=\bar{\phi}_-{+}\left(\bar{\phi}_+{-}\bar{\phi}_-\right)\mathcal{H}(x_1^{\prime }), \quad u_1^{\prime (b)}=U^{\prime }-U_s^{\prime }, \quad u_2^{\prime (b)}=0,\nonumber\\ \frac{\partial P^{\prime ^{(b)}}}{\partial x_1^{\prime }}={-}\frac{\mu_0^{\prime }U^{\prime }}{B^{\prime 2}M(\bar{\phi}^{(b)})}. \end{gather}$$

Here, the superscript $(b)$ denotes the base state, and $\mathcal {H}$ is the Heaviside step function. We then introduce normal mode perturbations to the base state volume fraction and pressure distributions,

(5.2a,b)\begin{align} \bar{\phi}=\bar{\phi}^{(b)}+\hat{\phi}({\hat{x}_1})\exp\left({({\sigma \hat{t}+\textrm{i}k{\hat{x}_2}})}\right), \quad P^{\prime }=P^{\prime (b)}+\frac{U^{\prime } \mu_{0} \chi}{B}\hat{P}({\hat{x}_1})\exp\left({({\sigma \hat{t}+\textrm{i}k{\hat{x}_2}})}\right). \end{align}

The perturbed equations obtained from (4.1)–(4.5) are rendered dimensionless, using the following scalings:

(5.3a–e)\begin{gather} \sigma =\dfrac{\sigma^{\prime }}{\left(\dfrac{U^{\prime }}{\chi B}\right)},\quad \hat{t}=\dfrac{t^{\prime }}{\left(\dfrac{\chi B^{\prime }}{U^{\prime }}\right)}, \quad k=k^{\prime }\chi B^{\prime }, \quad {\hat{x}_1}=\dfrac{x_1^{\prime }}{\chi B^{\prime }},\notag\\ \hat{P}=\dfrac{P^{\prime }}{\left(\dfrac{U^{\prime } \mu_{0} \chi}{B^{\prime }}\right)} \quad \text{with}\ \chi=\dfrac{B^{\prime 2}}{a^{\prime 2}}. \end{gather}

A linearization of the system of equations in (4.1)–(4.5) about the base state followed by the application of QSSA results in the following perturbed forms of the continuity and macro-transport equations:

(5.4a)\begin{align} &\dfrac{\textrm{d}^{2}\hat{P}}{\textrm{d}{\hat{x}_1}^{2}}-k^{2}\hat{P}=\dfrac{\dfrac{\textrm{d}M}{\textrm{d}\bar{\phi}}(\bar{\phi}^{(b)})}{M^{2}(\bar{\phi}^{(b)})}\dfrac{\textrm{d}\hat{\phi}}{\textrm{d}{\hat{x}_1}}, \end{align}

(5.4b)\begin{align} &\left[h(\bar{\phi}^{(b)})+\dfrac{\left(1-b\right)}{\chi^2}q_A(\bar{\phi}^{(b)})\right] \dfrac{\textrm{d}^2\hat{\phi}}{\textrm{d}\hat{x}_1^2}={-}\dfrac{\left(b+d\right)}{\chi^2}q_B(\bar{\phi}^{(b)})k^2\hat{u}_1\nonumber\\ &\quad +\left\{\dfrac{k^2}{\chi^2}\left[\left(1+d\right)q_A(\bar{\phi}^{(b)})+\left(1+d+\dfrac{(b+d)}{\chi}\right) q_B(\bar{\phi}^{(b)})\dfrac{\dfrac{\textrm{d}M}{\textrm{d}\bar{\phi}}(\bar{\phi}^{(b)})}{M(\bar{\phi}^{(b)})}\right]+\sigma\right\}\hat{\phi}\nonumber\\ &\quad +\left[\dfrac{\textrm{d}g}{\textrm{d}\bar{\phi}}(\bar{\phi}^{(b)})-\dfrac{U_s^{\prime }}{U^{\prime }}\right]\dfrac{\textrm{d}\hat{\phi}}{\textrm{d}\hat{x}_1}. \end{align}

The perturbed velocity $\hat {u}_1$ appearing in the above equation is given by

(5.5)\begin{equation} \hat{u}_1={-}M(\bar{\phi}^{(b)})\dfrac{\textrm{d}\hat{P}}{\textrm{d}\hat{x}_1} +\hat{\phi}\dfrac{\dfrac{\textrm{d}M}{\textrm{d}\bar{\phi}}(\bar{\phi}^{(b)})} {M(\bar{\phi}^{(b)})}. \end{equation}

The dimensionless ordinary differential equations (ODEs) in (5.4a), (5.4b) and (5.5) represent a system of four first-order ODEs in the variables $\hat {P}$, $\hat {u}_1$, $\hat {\phi }$ and $\textrm {d}\hat {\phi }/\textrm {d}\hat {x}_1$, and are applicable in each of the two domains, $\hat {x}_1<0$ and $\hat {x}_1>0$. The ODE system was solved using the eigenvector expansion technique on MATLAB^®. The boundary conditions employed were that the perturbed variables vanish as $\hat {x}_1\to -\infty$ and $\hat {x}_1\to +\infty$, and that $\hat {P}$, $\hat {u}_1$, $\hat {\phi }$ and the total particle flux in the $\hat {x}_1$ direction $(\hat {J}_1)$, derived separately for $\hat {x}_1<0$ and $\hat {x}_1>0$, match at $\hat {x}_1=0$. The expression for $\hat {J}_1$ is

(5.6)$$\begin{align} \hat{J}_1&=\dfrac{\left(1+d\right)}{\chi^2}q_B(\bar{\phi}^{(b)}) M(\bar{\phi}^{(b)})k^2\hat{P}+g(\bar{\phi}^{(b)})\hat{u}_1 +\left[\dfrac{\textrm{d}g}{\textrm{d}\bar{\phi}}(\bar{\phi}^{(b)})-\dfrac{U_s^{\prime }}{U^{\prime }}\right]\hat{\phi}\nonumber\\ &\quad -\left[h(\bar{\phi}^{(b)})+\dfrac{\left(1-b\right)}{\chi^2}q_A(\bar{\phi}^{(b)})\right] \dfrac{\textrm{d}\hat{\phi}}{\textrm{d}\hat{x}_1}. \end{align}$$

Figure 10 shows the relationship between $\sigma$ and $k$ for $B^{\prime }/a^{\prime }=20, 40$ and 100 for the case of $\hat {\phi }_-=0.4$ and $\hat {\phi }_{+}=0.5$. The behaviour is analogous to the observations of Tan & Homsy (Reference Tan and Homsy1986). For each aspect ratio $B/a$, $\sigma$ is positive in the long wave limit and increases as $k$ (see inset of figure 10). For short waves, however, the perturbations are stabilized by shear-induced migration normal to the flow direction. The dimensionless growth rate corresponding to the fastest growing mode, $\sigma _{m}$, is relatively unaffected by the aspect ratio (see figure 11a). However, the dimensionless wavenumbers $k_{m}$ and $k_{c}$ corresponding to $\sigma _{m}$ and the cut-off (see figure 11b,c), respectively, increase as the aspect ratio increases, i.e. a larger range of high frequency waves falls in the unstable region with an increase in $B^{\prime }/a^{\prime }$. This is due to a decrease in the shear-induced migration rate and an increase in the Taylor dispersivity with the increase in the aspect ratio. Figure 11 shows three key variables: $\sigma _{m}$, $k_{c}$ and $k_{m}$, over a range of $B^{\prime }/a^{\prime }$ values. We observe that when equivalent parameters are considered for the Hele-Shaw suspension flow problem and the solute problem (Tan & Homsy Reference Tan and Homsy1986), i.e. when the following substitutions are made,

(5.7a–c) \begin{equation} D_{{\parallel}}^{\prime }\sim\frac{U^{\prime }B^{\prime 3}}{a^{\prime 2}},\quad D_\perp^{\prime}{\sim}\frac{U^{\prime }a^{\prime 2}}{B^{\prime }} \quad \text{and}\quad \frac{D_\perp^{\prime }}{D_\parallel^{\prime }} \sim \frac{a^{\prime 4}}{B^{\prime 4}},\end{equation}

the results in figure 11 are in agreement with the theory of Tan & Homsy (Reference Tan and Homsy1986) in the limit $D_{\perp }^{\prime }/D_{\parallel }^{\prime }\ll 1$ (see table 1), which is valid for the suspension problem as $a^{\prime }\ll B^{\prime }$ in most suspension experiments. It follows that $k_{m}^{\prime }\propto a^{\prime 2/3}/B^{\prime 5/3}$, $k_c^{\prime }\propto {B^{\prime -1}}$ and $\sigma _{m}\propto U^{\prime }a^{\prime 2}/B{\prime ^3}$.

Figure 10. The dispersion relationship arising from a linear stability analysis of the macrotransport equation for a step increase in volume fraction in the flow direction from $\bar {\phi }_-=0.4$ to $\bar {\phi }_+=0.5$ (see figure 8b). The solid, dotted and dash–dotted curves are results for $B/a =$ 10, 40 and 100, respectively. The inset shows the same dispersion curves in the limit of small wavenumbers.

Figure 11. The maximum growth rate $\sigma _m^{\prime }$ (subfigure a), the wavenumber, $k_m^{\prime }$, corresponding to the maximum growth rate (subfigure b) and the cut-off wavenumber (subfigure c) as functions of the aspect ratio $B^{\prime }/a^{\prime }$ for $\bar {\phi }_-=0.4$ and $\bar {\phi }_+=0.5$.

Figure 12 examines the effect of volume fraction on the stability characteristics. It can be seen that for the same step change in volume fraction of 10 %, the growth rates and the range of wavenumbers for which $\sigma '$ is positive is greatly diminished for $\bar {\phi }_-=0.3$ and $\bar {\phi }_+=0.4$ as compared with $\bar {\phi }_-=0.4$ and $\bar {\phi }_+=0.5$. This is primarily due to the highly nonlinear influence of $\bar {\phi }$ on the Taylor dispersivity (see figure 4c), and to a lesser extent, on the suspension mobility (see figure 4a) and shear-induced migration rates (see figure 4e,f). It is also instructive to examine the role of the excess of the flow-averaged concentration over the depth-averaged concentration, i.e. $g (\bar {\phi })>\bar {\phi }$ on the stability curves, since this is a feature that is different from the work of Wooding (Reference Wooding1962) and Tan & Homsy (Reference Tan and Homsy1986) In figure 13(a) we have shown $\sigma '$ vs $k'$ for $B^{\prime }/a^{\prime }=40$, $\bar {\phi }_-=0.4$ and $\bar {\phi }_+=0.5$ with $g (\bar {\phi })$ set manually to $\bar {\phi }$ (dash–dotted curve). As can be seen, in the absence of a difference between $g (\bar {\phi })$ and $\bar {\phi }$, the growth rates are enhanced at all the shown wavenumbers, with the largest increase occurring near $k_m^{\prime }$. On the other hand, for $B^{\prime }/a^{\prime }=40$, $\bar {\phi }_-=0.3$ and $\bar {\phi }_+=0.4$, the opposite effect is observed (see figure 13b); the excess of $g (\bar {\phi })$ over $\bar {\phi }$ diminishes the growth rates of perturbations. This can be explained as follows. At the higher volume fractions, the shock velocity $U_s^{\prime }$ is $0.95U'$, i.e. less than $U^{\prime }$. This situation is similar to the meniscus accumulation phenomenon (Karnis & Mason Reference Karnis and Mason1967; Chapman Reference Chapman1990; Tang et al. Reference Tang, Grivas, Homentcovschi, Geer and Singler2000; Ramachandran & Leighton Reference Ramachandran and Leighton2007a; Ramachandran Reference Ramachandran2013; Luo, Chen & Lee Reference Luo, Chen and Lee2018), whereby particles accumulate at an advancing interface due to the difference in the particle average velocities behind the meniscus and at the meniscus. The greater the displacement, the thicker the accumulated layer at the shock front. As shown in figure 14, when a perturbed miscible interface is displaced in the flow direction, turning on the physical effect of $U_s^{\prime }< U^{\prime }$ would lead to an accumulation of particles at the displaced interface compared with the case when $U_s^{\prime }$ is equal to $U^{\prime }$. The accumulation manifests as a reduced growth rate of the perturbation. The situation would be exactly the opposite for the case when $U_s^{\prime }>U^{\prime }$, which occurs for $\bar {\phi }_-=0.3$ and $\bar {\phi }_+=0.4$. Consequently, there would be a depletion of particles at the perturbed interface compared with the case of $U_s^{\prime }=U^{\prime }$, leading to higher growth rates when $g (\bar {\phi })>\bar {\phi }$.

Figure 12. The effect of volume fraction on the stability curve, $\sigma '$ vs $k'$. For the solid curve, $\bar {\phi }_-=0.4$ and $\bar {\phi }_+=0.5$, while for the dash–dotted curve, $\bar {\phi }_-=0.3$ and $\bar {\phi }_+=0.4$. For both curves, $B^{\prime }/a^{\prime }=40$. Larger volume fractions lead to higher growth rates and a larger range of wavenumbers for which perturbations are unstable.

Figure 13. The effect of the excess of the flow-averaged concentration $g$ over the area-averaged concentration $\bar {\phi }$ on the relationship between $\sigma '$ and $k'$; (a) $\bar {\phi }_-=0.4$ and $\bar {\phi }_+=0.5$, (b) $\bar {\phi }_-=0.3$ and $\bar {\phi }_+=0.4$. For both subfigures, $B^{\prime }/a^{\prime }=40$. In each subfigure, the dash–dotted line is the result of a simulation where $g$ was manually set to be equal to $\bar {\phi }$.

Figure 14. A picture explaining why the case of $U_s^{\prime }< U^{\prime }$ would lead to lower growth rates compared with the case of $U_s^{\prime }=U^{\prime }$. A low volume fraction suspension displaces a high volume fraction suspension from left to right.

We conclude this section with two comments. First, the results above are obtained by invoking the QSSA for a step change in the volume fraction at $t=0$. At longer times, the base state volume fraction distribution will eventually relax to produce gradients over a length scale of $B^{\prime 3}/a^{\prime 2}$. As a consequence, the maximum growth rate and wavenumber will diminish, as observed by Tan & Homsy (Reference Tan and Homsy1986). An investigation of the stability characteristics of the base state solution for $t>0$ will be pursued in the future. Second, the critical viscosity predicted from our theory is unity, i.e. the suspension being displaced needs only to be more viscous than the upstream suspension to initiate the instability. However, it has been reported in several recent studies that a critical viscosity contrast less than unity has to be crossed to observe the viscous miscible fingering instability, and this has been attributed to the three-dimensional distribution of the two fluids at the mixing front (Lajeunesse et al. Reference Lajeunesse, Martin, Rakotomalala, Salin and Yortsos1999; Bischofberger, Ramachandran & Nagel Reference Bischofberger, Ramachandran and Nagel2014; Videbaek & Nagel Reference Videbaek and Nagel2019; Luo, Chen & Lee Reference Luo, Chen and Lee2020). Our stability analysis is implemented on a macrotransport equation obtained from an asymptotic expansion, and the inclusion of the Taylor dispersion term does account for the penetration of two fluids in the mixing zone (see the discussion around figure 3). But the observation of a critical viscosity contrast to trigger the instability would imply that we would need to extend the expansion to higher orders to capture the physics. This extension is also left to future work.