2. Mathematical description of the problem

We assume that the suspended particles are sufficiently small that the buoyancy force is in balance with the Stokes drag and the particle inertia is negligible. Thus, the Eulerian description for the volume fraction of particles, $\phi ^*$, can be expressed as

(2.1)\begin{equation} \dfrac{\textrm{D}^*}{\textrm{D}^* t^*}{\phi^*}+V_0^*\left(\boldsymbol{\nabla}^*\phi^*\boldsymbol{\cdot}\hat{\boldsymbol{e}}_g\right) = \kappa^* \nabla^{*2} \phi^*, \end{equation}

where $\textrm {D}^*/\textrm {D}^* t^*$ is the material derivative, $V_0^*$ is the Stokes settling speed obtained by the balance between the particle drag and buoyant force, $\hat {\boldsymbol {e}}_g$ is the unit direction vector of the gravitational acceleration $g^*$, $\kappa ^*$ is the diffusivity and the superscript $*$ indicates the dimensional quantity of the physical variables. Following Herbolzheimer & Acrivos (Reference Herbolzheimer and Acrivos1981), we assume a dilute suspension, such that the viscosity of the liquid phase is not altered by the presence of the suspended particles. As shown in the schematic in figure 2(b), for a water column with an inclined angle $\theta$, under the Boussinesq approximation, the mass and momentum equations are as follows:

(2.2)$$\begin{gather} \boldsymbol{\nabla}^*\boldsymbol{\cdot}\boldsymbol{u}^* = 0, \end{gather}$$

(2.3)$$\begin{gather}\frac{\textrm{D}^*}{\textrm{D}^* t^*}\boldsymbol{u}^* ={-}\frac{1}{\rho^*}\boldsymbol{\nabla}^* p^* + \boldsymbol{\nabla}^* \mathcal{B}^* + \nu^*\nabla^{*2}\boldsymbol{u}^*, \end{gather}$$

where $p^*$ is the modified pressure term obtained by subtracting the hydrostatic pressure of the clear-fluid column from the total pressure field and $\boldsymbol {\nabla }^*\mathcal {B}^*$ is the body-force term written in a conservative form. After subtracting the hydrostatic pressure of the clear-fluid column, $\boldsymbol {\nabla }^*\mathcal {B}^*$ becomes the buoyant forcing due to $\phi ^*$ and can be written as

(2.4)\begin{equation} \mathcal{B}^* = \int^{\boldsymbol{x}^*}_{\boldsymbol{x}^*_0}\phi^* g^{\prime *}(\check{\boldsymbol{x}}^*) \hat{\boldsymbol{e}}_g\boldsymbol{\cdot} \textrm{d}\check{\boldsymbol{x}}^{*}, \end{equation}

where $\boldsymbol {x}^* = (x^*, y^*, z^*)$, $\boldsymbol {x}_0^*$ represents the reference position, $g^{\prime *} = (\rho _s^*/\rho ^*-1)g^*$ is the reduced gravity and $\hat {\boldsymbol {e}}_g = (-\cos {\theta }\hat {\boldsymbol {e}}_x, -\sin {\theta }\hat {\boldsymbol {e}}_y)$ represents the unit direction vector of $g^{\prime *}$ in the present configuration.

For the high-aspect-ratio case in which the liquid column is assumed to be fairly long in the $x$-direction, we consider the fully developed flow section where the variations in the unperturbed momentum and concentration in the streamwise direction ($x$) are zero. Without variation in both $x$- and $z$-directions (see figure 2b), the transport equation, (2.1), of $\phi ^*$ can be reduced to

(2.5)\begin{equation} \dfrac{\partial\phi^*}{\partial t^*}-V_0^*\sin\theta\dfrac{\partial \phi^*}{\partial y^*} = \kappa^*\dfrac{\partial^2\phi^*}{\partial y^{* 2}}. \end{equation}

It should be noted that in (2.5), if one considers only the molecular diffusivity, the value of $\kappa ^*$ can be infinitesimally small according to the Stokes–Einstein equation. Examples are Herbolzheimer & Acrivos (Reference Herbolzheimer and Acrivos1981) and Davis et al. (Reference Davis, Herbolzheimer and Acrivos1983), in which the problem configurations become two fluid flows separated by a density jump. More recent studies have shown that additional ‘hydrodynamic diffusion’ exists owing to particle–particle and particle–hydrodynamics interactions, which can dominate the diffusion of the suspended particles over the molecular effect (e.g. Davis Reference Davis1996; Segre et al. Reference Segre, Liu, Umbanhowar and Weitz2001). However, hydrodynamic diffusion is still not a well-defined physical quantity and its magnitude can vary among different flow configurations. In this regard, different values of $\kappa ^*$ are studied in § 6, where we show that the instability characteristics asymptotically converge when $\kappa ^*$ becomes infinitesimally small. Typically, when the volume fraction of suspended particles is modelled using (2.5), the bottom boundary condition is specified as $V_0^* \phi ^* + \kappa ^* \partial \phi ^*/\partial y^* = 0$. However, this is usually the case where flow turbulence exists to support resuspension (e.g. Chou & Fringer Reference Chou and Fringer2008), and the diffusivity ($\kappa ^*$) should be replaced by the eddy diffusivity. In this study, turbulence is not considered and, because $\kappa ^*$ is typically very small, the bottom wall is a deposition-dominant regime such that at the bottom boundary $V_0^* \phi ^* + \kappa ^* \phi ^*/\partial y^* > 0$ becomes deposition.

If the problem starts from a fully suspended water column with a volume fraction $c_0^*$, the layer of deposition is usually thin and its thickness is negligible relative to the domain width (i.e. $\delta _d = 0$ in figure 2b). Therefore, we only consider the layer in which $\phi ^* \leqslant c_0^*$ (the grey areas in figure 2) as the particle-laden liquid column, which eliminates the complexities while dealing with the transition to the dense deposition zone on the wall. Under this condition, because $\phi ^*$ is uniform within the suspension layer in our domain of interest except at the diffusive interface where $0 < \phi ^* \leqslant c_0^*$, $\kappa ^* \partial \phi ^*/\partial y^*$ becomes zero on the bottom boundary. Thus, the bottom boundary condition for solving (2.5) in this study is to specify a flux $-V_0^*\phi ^*$, such that $\phi ^*$ at the bottom remains as $c^*_0$, except when the diffusive interface (i.e. $\phi ^* < c^*_0$) approaches the bottom boundary near the end of the calculation. In fact, starting from a fully suspended water column with a volume fraction $c_0^*$, using the boundary conditions $\phi ^* = 0$ at the top ($y^* = b^*$) and specifying a flux $-V_0^*\phi ^*$ at the bottom, solving (2.5) (Ogata & Banks Reference Ogata and Banks1961) gives a descending diffusive interface specified by a time-varying error function given by

(2.6)\begin{equation} \phi^*(\zeta^*) = \dfrac{1 - \mathrm{erf}\left(\dfrac{\zeta^*}{2}\right)}{2}c_0^*, \end{equation}

where

(2.7)\begin{equation} \zeta^*\left(y^*, t^*\right) = \dfrac{y^*-b^*+V_0^*\sin\theta t^*}{\sqrt{\kappa^* t^*}}. \end{equation}

2.1. Scaling

The velocity scale is obtained by calculating the upper limit of the velocity magnitude for the undisturbed flow field in the steady state under the condition that both $\kappa ^*$ and $V_0^*$ are zero and that the interface is located at the middle of the vertical dimension (i.e. $y^* = 0.5 b^*$). Assuming a fully developed, steady undisturbed flow (i.e. $\partial \boldsymbol {u}^*/\partial x^* = 0$) without any variation in the $z$-direction, as shown in figure 2(c), the governing equations can be reduced to a momentum equation in the $x$-direction of a planar flow written as

(2.8)$$\begin{gather} 0 ={-}\dfrac{\partial p^*}{\partial x^*}-\phi^*g^{\prime *}\cos\theta+\nu^*\dfrac{\partial^2 u^*}{\partial y^{* 2}}, \end{gather}$$

(2.9)$$\begin{gather}0 ={-}\dfrac{\partial p^*}{\partial y^*}-\phi^*g^{\prime *}\sin\theta, \end{gather}$$

where

(2.10) \begin{equation} \phi^* = \left\{\begin{array}{@{}ll} c^*_0 & \mbox{if}\ 0\leqslant y^*<\dfrac{b^*}{2}\\ 0 & \mbox{if}\ \dfrac{b^*}{2}\leqslant y^*< b^* \end{array},\right. \end{equation}

That is, the flow, governed by (2.8) and (2.9), becomes a two-layered flow separated by the interface of $\phi ^*$ (see figure 2c).

For a water column with a long but finite length in the $x$-direction, owing to the existence of end walls in the $x$-direction, the net flow at each cross-section must be zero. By imposing this zero-net-flow condition at the cross-section, along with the no-slip boundary condition at the walls and continuity for $u^*$ and stress ($\nu ^* \partial u^*/\partial x^*)$ at the interface ($y^* = b^*/2$), the solution can be obtained as

(2.11) \begin{equation} u^* = \left\{ \begin{array}{@{}ll} \dfrac{1}{4\nu^*}c_0^*g^{\prime *}\cos\theta\left(y^{*2}-\dfrac{b^*}{2}y^{*}\right) & \mbox{if}\ 0 \leqslant y^* < \dfrac{b^*}{2}\\ \dfrac{1}{4\nu^*}c_0^*g^{\prime *}\cos\theta\left({-}y^{*2}+\dfrac{3b^*}{2}y^{*}-\dfrac{b^{* 2}}{2}\right) & \mbox{if}\ \dfrac{b^*}{2} \leqslant y^* < b^* \end{array},\right. \end{equation}

which gives a two-layered flow in counter directions. Substituting $y^* = b^*/4$ into the top identity or $y^* = 3b^*/4$ into the bottom identity of (2.11), the maximum magnitude of the velocity can be found as $c_0^*g^{\prime *} b^{* 2}\cos \theta / (64\nu ^*)$. Thus, we use

(2.12)\begin{equation} U_c^* = \frac{c_0^*g^{\prime *} b^{* 2}}{64\nu^*} \end{equation}

as the velocity scale of the present study.

For an inclined water column associated with a long dimension in the $x$-direction, we use $b^*$, the domain length in the $y$-direction (see figure 2b), as our characteristic length scale, $U_c^*$ (2.12) as the velocity scale, $\mu ^* U_c^*/ b^*$ as the pressure scale and $c_0^*$ as the scale for the volume fraction. Substituting these dimensional parameters into (2.1)–(2.4) gives the dimensionless transport equation of $\phi$, written as

(2.13)\begin{equation} \frac{\partial \phi}{\partial t} + \left(\boldsymbol{u}-R_v\sin{\theta}\hat{\boldsymbol{e}}_y\right)\boldsymbol{\cdot}\boldsymbol{\nabla}\phi = \frac{1}{RePr}\nabla^2\phi, \end{equation}

and the dimensionless mass and momentum equations for the fluid, written as

(2.14)$$\begin{gather} \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u} = 0, \end{gather}$$

(2.15)$$\begin{gather}\frac{\partial}{\partial t}\boldsymbol{u} + \boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla}\boldsymbol{u} ={-}\frac{1}{Re}\boldsymbol{\nabla} p + \boldsymbol{\nabla}\int^{\boldsymbol{x}}_{\boldsymbol{x}_0} \frac{64}{Re}\phi\hat{\boldsymbol{e}}_g\boldsymbol{\cdot} \textrm{d} \check{\boldsymbol{x}} + \frac{1}{Re}\nabla^2\boldsymbol{u}, \end{gather}$$

where

(2.16a–c)\begin{equation} Re = \frac{U_c^* b^*}{\nu^*},\quad R_v = \frac{V^*_0}{U_c^*}, \quad \mbox{and}\quad Pr = \frac{\nu^*}{\kappa^*}, \end{equation}

are the Reynolds number, velocity ratio and Prandtl number, respectively.

2.2. Specification of time-dependent base flow

Each physical variable is expressed as a summation of its base state and perturbation, which can be written as

(2.17) \begin{equation} \left.\begin{array}{rclcr@{}} \displaystyle u(x,y,z,t) & = & U(y,t) & + & \tilde{u}(x,y,z,t), \\ v(x,y,z,t) & = & & & \tilde{v}(x,y,z,t), \\ w(x,y,z,t) & = & & & \tilde{w}(x,y,z,t), \\ p(x,y,z,t) & = & P(x,y,t) & + & \tilde{p}(x,y,z,t), \\ \phi(x,y,z,t) & = & \varPhi(y,t) & + & \tilde{\phi}(x,y,z,t). \end{array}\right\} \end{equation}

Substitution of (2.17) into (2.13)–(2.15) yields equations describing the base flow and the temporal evolution of disturbances. The base-flow equation is as follows:

(2.18)$$\begin{gather} \dfrac{\textrm{d} U}{\textrm{d} x} = 0\,, \end{gather}$$

(2.19)$$\begin{gather}\frac{\partial U}{\partial t} ={-}\frac{1}{Re}\frac{\partial P}{\partial x} - \frac{64}{Re}\cos{\theta}\varPhi + \frac{1}{Re}\frac{\partial^2 U}{\partial y^2} , \end{gather}$$

(2.20)$$\begin{gather}0 ={-}\frac{1}{Re}\frac{\partial P}{\partial y} - \frac{64}{Re} \sin{\theta}\varPhi , \end{gather}$$

(2.21)$$\begin{gather}\frac{\partial \varPhi}{\partial t} - R_v\sin{\theta}\frac{\partial\varPhi}{\partial y} = \frac{1}{RePr}\frac{\partial^2\varPhi}{\partial y^2}, \end{gather}$$

where $y \in [0, 1]$. Starting from a water column fully occupied by particle-laden fluid (i.e. $\varPhi (y, 0) = 1$), using the normalised boundary conditions, $\varPhi = 0$ at $y = 1$ and the mass flux $= -\varPhi R_v\sin \theta$ at the bottom wall, the solution of (2.21) can be simply obtained as

(2.22)\begin{equation} \varPhi(\zeta) = \frac{1-\mathrm{erf}\left(\dfrac{\zeta}{2}\right)}{2}, \end{equation}

where

(2.23)\begin{equation} \zeta(y, t) = \frac{\sqrt{PrRe} \left(y - 1 + R_v \sin{\theta} t\right)}{\sqrt{t}}. \end{equation}

To obtain the base flow $U$ by solving (2.19), we note the net flow at each cross-section must be zero, as previously mentioned in § 2.1. That is, using the forcing owing to the time-varying $\varPhi$ given by (2.22), with the no-slip boundary condition at both the bottom and top walls and the condition of zero net flow in the cross-section, one can obtain the time-varying base-flow solution ($U$) by solving (2.19) numerically. Here, we use the second-order central difference for spatial discretisation and second-order Crank–Nicolson method for time discretisation to obtain $U(y, t)$.

2.3. Temporal evolution of disturbances

In this study, we examine the transient instability by measuring the growth of the total energy, including the kinetic and potential energies. The kinetic energy can be obtained from the wall-normal velocity and the wall-normal vorticity (Schmid & Henningson Reference Schmid and Henningson2001; Schmid Reference Schmid2007), the potential energy in a generalised form (South & Hooper Reference South and Hooper1999; Sameen & Govindarajan Reference Sameen and Govindarajan2007; Jerome et al. Reference Jerome, Chomaz and Huerre2012) can be converted from the volume fraction $\phi$. Let a state vector of disturbances, $\tilde {\boldsymbol {q}}$, be composed of disturbances of the wall-normal velocity $\tilde {v}$, wall-normal vorticity $\tilde {\eta }$ and volume fraction $\tilde {\phi }$. The energy of these disturbances can be defined as a weighted energy norm written as

(2.24)\begin{equation} \langle \tilde{\boldsymbol{q}}, \tilde{\boldsymbol{q}} \rangle_E = \tilde{\boldsymbol{q}}^H\boldsymbol{\mathsf{E}}\tilde{\boldsymbol{q}}, \end{equation}

where

(2.25a–c) \begin{equation} \tilde{\boldsymbol{q}} = \left[\begin{array}{@{}c@{}} \tilde{v}\\ \tilde{\eta}\\ \tilde{\phi} \end{array} \right],\quad \boldsymbol{\mathsf{E}} = \left[ \begin{array}{@{}ccc@{}} \dfrac{1}{2k^2}\left(k^2 - \mathcal{D}^2\right) & 0 & 0 \\ 0 & \dfrac{1}{2k^2} & 0\\ 0 & 0 & \dfrac{1}{2}RePr \end{array}\right],\quad k = \sqrt{\alpha^2 + \beta^2} \end{equation}

the superscript $H$ denotes the Hermitian transpose, $\boldsymbol{\mathsf{E}}$ is the energy conversion matrix, $\alpha$ is the streamwise wavenumber, $\beta$ is the spanwise wavenumber and $\mathcal {D}^2 = \textrm {d}^2/{\textrm {d}y}^2$ ($\mathcal {D} = \textrm {d}/{\textrm {d} y}$). Given that any disturbance input can be amplified differently with time, we are interested in initial perturbations that lead to maximum amplification. Therefore, we are interested in finding the initial disturbances that generate the optimal energy gain, $G_n(t; t_0)$, at each time point. The optimal energy gain within a time interval $[t_0, t]$ can be described as

(2.26)\begin{equation} G_n(t; t_0) = \mathop{\arg\max}\limits_{\tilde{\boldsymbol{q}}(t_0)}\frac{\langle \tilde{\boldsymbol{q}}(t), \tilde{\boldsymbol{q}}(t) \rangle_E}{\langle \tilde{\boldsymbol{q}}(t_0), \tilde{\boldsymbol{q}}(t_0)\rangle_E} = \mathop{\arg\max}\limits_{\tilde{\boldsymbol{q}}_0}\frac{\langle \tilde{\boldsymbol{q}}, \tilde{\boldsymbol{q}} \rangle_E}{\langle \tilde{\boldsymbol{q}}_0, \tilde{\boldsymbol{q}}_0\rangle_E}. \end{equation}

The resulting optimisation problem can be solved by propagating the disturbances forward and backward in time using the governing and corresponding adjoint equations (Farrell & Moore Reference Farrell and Moore1992; Schmid Reference Schmid2007; Luchini & Bottaro Reference Luchini and Bottaro2014).

The normal mode formulation for the components of $\tilde {\boldsymbol {q}}$ (2.24) is written as

(2.27) \begin{equation} \begin{bmatrix} \tilde{v}\\ \tilde{\eta}\\ \tilde{\phi} \end{bmatrix} =\begin{bmatrix} \hat{v}(y,t)\\ \hat{\eta}(y,t)\\ \hat{\phi}(y,t) \end{bmatrix} \exp{\mathrm{i}(\alpha x + \beta z)}. \end{equation}

Based on Equations (A9), (A12), and (A13) derived in Appendix A, the evolution of the disturbance field can be obtained by solving the Orr–Sommerfeld equation with baroclinic forcing, the Squire equation and the transport equation in a compact notation, written as follows:

(2.28)\begin{equation} \left(\boldsymbol{\mathsf{M}}\frac{\partial}{\partial t} - \boldsymbol{\mathsf{L}}\right)\hat{\boldsymbol{q}} = 0, \end{equation}

with

(2.29a,b) \begin{equation} \boldsymbol{\mathsf{M}} = \begin{bmatrix} \mathcal{D} - k^2 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{bmatrix} \quad\mbox{and}\quad \boldsymbol{\mathsf{L}} = \begin{bmatrix} \mathcal{L}_{OS} & 0 & \mathcal{L}_{Cv} \\ \mathcal{L}_{C\eta} & \mathcal{L}_{SQ} & 0 \\ \mathcal{L}_{C\phi} & 0 & \mathcal{L}_{\phi} \end{bmatrix}, \end{equation}

where

(2.30a)$$\begin{gather} \mathcal{L}_{OS} = \frac{1}{Re}(\mathcal{D}^2-k^2)^2 - \mathrm{i}\alpha U(\mathcal{D}^2-k^2) + \mathrm{i}\alpha \mathcal{D}^2U, \end{gather}$$

(2.30b)$$\begin{gather}\mathcal{L}_{Cv} ={-}\frac{64}{Re}k^2\sin{\theta}, \end{gather}$$

(2.30c)$$\begin{gather}\mathcal{L}_{C\eta} ={-}\mathrm{i}\beta\mathcal{D}U, \end{gather}$$

(2.30d)$$\begin{gather}\mathcal{L}_{SQ} = \frac{1}{Re}(\mathcal{D}^2-k^2) - \mathrm{i}\alpha U, \end{gather}$$

(2.30e)$$\begin{gather}\mathcal{L}_{C\phi} ={-}\mathcal{D}\varPhi, \end{gather}$$

(2.30f)$$\begin{gather}\mathcal{L}_{\phi} = \frac{1}{RePr} (\mathcal{D}^2 - k^2) - \mathrm{i}\alpha U + R_v\sin\theta \mathcal{D}. \end{gather}$$

Because all disturbance quantities vanish on the walls, boundary conditions are given with

(2.31a–d)\begin{equation} \hat{v} = 0, \quad \mathcal{D}\hat{v} = 0, \quad \hat{\eta} = 0, \quad \hat{\phi} = 0 \quad\mbox{on}\ \partial\varOmega_y, \end{equation}

where $\partial \varOmega _y$ denotes the boundaries in the $y$-direction. It should be noted that the boundary condition $\hat {\phi } = 0$ on $\partial \varOmega _y$ is not consistent with the condition for the base flow, where we specify a flux for transport of the volume fraction at the bottom boundary. However, as we also solve for a set of adjoint equations in what follows, a flux-type boundary condition makes it very difficult to obtain its adjoint counterpart. In addition, because the instabilities of the momentum and volume fraction of particles are fully coupled, we set $\hat {\phi } = 0$ at the bottom boundary because of vanishing momentum disturbance in (2.31). Here, $\hat {\boldsymbol {q}}$ at time $t$ can be linked to the initial disturbance $\hat {\boldsymbol {q}}_0$ at time $t_0$ by an evolution matrix $\boldsymbol{\mathsf{A}}(t)$,

(2.32)\begin{equation} \hat{\boldsymbol{q}} = \boldsymbol{\mathsf{A}}(t)\hat{\boldsymbol{q}}_0, \end{equation}

which represents the time integration of (2.28) from the initial time $t_0$ to a later time $t$. Equation (2.26) can then be rewritten as follows:

(2.33) \begin{equation} G_n(t; t_0) = \mathop{\arg\max}\limits_{\tilde{\boldsymbol{q}}_0}\frac{\langle \tilde{\boldsymbol{q}}, \tilde{\boldsymbol{q}} \rangle_E}{\langle \tilde{\boldsymbol{q}}_0, \tilde{\boldsymbol{q}}_0\rangle_E} = \mathop{\arg\max}\limits_{\hat{\boldsymbol{q}}_0}\frac{\langle\boldsymbol{\mathsf{A}}(t)\hat{\boldsymbol{q}}_0, \boldsymbol{\mathsf{A}}(t)\hat{\boldsymbol{q}}_0 \rangle_E}{\langle \hat{\boldsymbol{q}}_0, \hat{\boldsymbol{q}}_0\rangle_E} = \mathop{\arg\max}\limits_{\hat{\boldsymbol{q}}_0}\frac{\langle\boldsymbol{\mathsf{A}}^{{\dagger}}(t)\boldsymbol{\mathsf{A}}(t)\hat{\boldsymbol{q}}_0, \hat{\boldsymbol{q}}_0 \rangle_E}{\langle \hat{\boldsymbol{q}}_0, \hat{\boldsymbol{q}}_0\rangle_E}, \end{equation}

where the superscript ${{\dagger}}$ denotes the adjoint operator and $\boldsymbol{\mathsf{A}}^{{\dagger}} (t)$ is the evolution matrix of the associated adjoint equation:

(2.34)\begin{equation} \left(-\boldsymbol{\mathsf{M}}^{{\dagger}}\frac{\partial}{\partial t} - \boldsymbol{\mathsf{L}}^{{\dagger}}\right)\hat{\boldsymbol{q}}^+= 0, \end{equation}

where the superscript $+$ denotes the complex conjugate. We note that in (2.28), $\boldsymbol{\mathsf{M}}$ is a self-adjoint operator, such that $\boldsymbol{\mathsf{M}}^{{\dagger}} = \boldsymbol{\mathsf{M}}$ in (2.34). The boundary conditions for the solutions of (2.34) are

(2.35a–d)\begin{equation} \hat{v}^+= 0, \quad \mathcal{D}\hat{v}^+= 0, \quad \hat{\eta}^+= 0, \quad \hat{\phi}^+= 0 \quad\mbox{on}\ \partial\varOmega_y. \end{equation}

The last identity in (2.33) implies that finding an initial disturbance that causes the maximum gain is equivalent to seeking an eigenvector that corresponds to the maximum eigenvalue of $\boldsymbol{\mathsf{A}}^{{\dagger}} (t)\boldsymbol{\mathsf{A}}(t)$. This can be accomplished by solving the ordinary differential equation system composed of (2.28) and (2.34). This process begins from guessing an initial, random disturbance field $\hat {\boldsymbol {q}}_0$ at the initial time $t_0$. The initial disturbance is propagated to a later time $t$ using (2.28) with an appropriate time-marching scheme. The resulting disturbance at $t$ is then propagated back to $t_0$ using (2.34). Variations in the resulting disturbance and energy gain between two consecutive loops decrease as the number of loops increases. For details of the computational procedure used, readers are referred to Schmid & Brandt (Reference Schmid and Brandt2014). In our case, the termination criterion for the iterations is that the relative variation between two consecutive loops needs to be less than $10^{-5}$. Equation (2.28) and its adjoint form (2.34) are evolved using the Chebyshev collocation method (Weideman & Reddy Reference Weideman and Reddy2000; Trefethen Reference Trefethen2019) in space and the second-order marching scheme in time accompanied by the implicit Euler method as a starter following Hack & Zaki (Reference Hack and Zaki2015). The $(n+1)$th step of the state vector is thus calculated using

(2.36)\begin{equation} \hat{\boldsymbol{q}}^{(n+1)} = (3\boldsymbol{\mathsf{M}}-2\Delta t\boldsymbol{\mathsf{L}}^{(n+1)})^{{-}1}\left(\boldsymbol{\mathsf{M}}\left(4\hat{\boldsymbol{q}}^{(n)}-\hat{\boldsymbol{q}}^{(n-1)}\right)\right). \end{equation}

We adopted Chebyshev points of the second kind to discretise the physical domain in the wall-normal direction. As the spectral grids are defined within the interval $[-1,1]$, the derivative with respect to $y$ must be stretched by a factor of two. To impose $m$ boundary constraints, an $N$-point second-kind Chebyshev polynomial interpolant of order $N$ is downsampled on to an $N$-point first-kind Chebyshev polynomial interpolant of order $(N-m)$, which excludes the boundary points at $y = \pm 1$. In this way, an $N$-by-$N$ differentiation matrix becomes an $(N-m)$-by-$N$ matrix, whose vanishing rows are then replenished by the boundary information, following the work of Driscoll & Hale (Reference Driscoll and Hale2016) and Trefethen (Reference Trefethen2019).

Computations were performed using the Julia programming language (Bezanson et al. Reference Bezanson, Edelman, Karpinski and Shah2017), and the packages of Cairo.jl (Bezanson et al. Reference Bezanson, Fischer, Holy, Nolta and Shah2012), Dierckx.jl (Barbary et al. Reference Barbary, Sundar and Ortner2020), FFTW.jl (Frigo & Johnson Reference Frigo and Johnson2005), Makie.jl (Danisch et al. Reference Danisch and Krumbiegel2021), JLD.jl (Holy et al. Reference Holy, Hinsen, Short and Kornblith2020), Roots.jl (Verzani et al. Reference Verzani, Ignatiev and Rackauckas2021), SpecialFunctions.jl (Johnson et al. 2021), and ToeplitzMatrices.jl (Noack et al. Reference Noack, Olver and Xu2021).

3. Effect of particle settling on the base flow

In this section, we investigate changes in the base flow in response to a fast-settling, particle-laden layer. To examine this clearly, (2.19) is rearranged to become

(3.1)\begin{equation} \frac{\partial U}{\partial t} = \frac{1}{Re_s\sin{\theta}}\left(-\frac{\partial P}{\partial x} - 64\cos{\theta}\varPhi + \frac{\partial^2 U}{\partial y^2}\right), \end{equation}

where

(3.2)\begin{equation} Re_s = R_v Re = \frac{V_0^* b^*}{\nu^*} \end{equation}

is the sedimentation Reynolds number, which indicates the balance between particle settling and the viscous effect. In (3.1), we can see that, given a fixed inclined angle, $\theta$, so that the buoyant forcing, $64\cos {\theta }\varPhi$, is fixed, $Re_s$ appears to be the only parameter that controls the base flow. Figure 3 presents profiles of the streamwise velocity for various $Re_s\sin {\theta }$ with $\theta = 30^\circ$, $Re = 215$, and $PrRe = 1000$ at ten time points, which show that typically, base flows of various $Re_s\sin {\theta }$ attain maximum speed at the time point at which the particle-laden interface reaches the middle of the domain in the vertical direction. However, as $Re_s\sin \theta$ increases (i.e. fast settling), the development of the base flow is somewhat delayed relative to the settlement of the particle-laden interface, as indicated by the grey lines in the figure. In figure 3, we can also see that the maximum speed the base flow can reach decreases as $Re_s\sin {\theta }$ increases. Moreover, the temporal development of the velocity profile converges with decreasing $Re_s\sin {\theta }$. In (3.1), $Re_s\sin \theta$ indicates how fast the flow can adapt to any change owing to the settlement of the particle-laden interface. If $Re_s\sin {\theta } \ll 1$, one anticipates a quasi-steady state in which the flow can adjust itself rapidly in response to the combined forcing of the pressure and buoyancy. In this case, the flow can attain its maximum strength at each time step. In contrast, when $Re_s\sin {\theta } \gg 1$, the flow cannot react instantly to the particle settling, so the flow develops after the descent of the particle-laden interface.

Figure 3. Time series of the streamwise velocity profiles for various $Re_s\sin {\theta }$. The grey horizontal lines indicate positions of the particle-laden interface.

4. Energy growth of disturbances

For the sake of clarity, here, we first focus on a representative case with $R_v = 0.01$, $Re = 215$, $PrRe = 1000$, $\theta = 30^\circ$ and $(\alpha, \beta ) = ({\rm \pi} /2, {\rm \pi}/2)$, which is denoted as the base case hereafter. In this case, values for $\alpha$ and $\beta$ are chosen such that both the algebraic and exponential instabilities can be clearly seen, as discussed in the following. Figure 4 presents time histories of various $G_n(t;t_0)$. In this figure, each curve represents different $G_n(t;t_0)$ resulting from the disturbances initialised at different time instants, $t_0$, during the early time stage (see inset (a)). Typically, the energy gain $G_n(t;t_0)$ undergoes algebraic growth and decay in the early stage (see inset (a) in figure 4), followed by a duration of stable flow development (i.e. $G_n(t;t_0) = 1$), and then exponential growth to a peak magnitude significantly greater than that of the early algebraic case. In the present cases, owing to the descent of the particle-laden interface, the buoyant forcing resulting from the particle-laden layer eventually vanishes, such that the flow energy always decays at the end. In figure 4, the time step $t = 215$ is the time point at which particles completely settle on to the upper-facing wall. Hereafter, we also consider the maximum amplification of the energy gain at time $T$ acquired from perturbations introduced at all time instants $t$ during the time interval $[t_0, T]$. Our focus is on the maximum energy gain, $G(t)$, caused by any disturbances, which is defined as follows (Blesbois et al. Reference Blesbois, Chernyshenko, Touber and Leschziner2013):

(4.1)\begin{equation} G(t) = \mathop{\arg\max}\limits_{t_0}{G_n}(t; t_0). \end{equation}

An example of $G$ for the base case is highlighted by the black dash-dotted line in figure 4.

Figure 4. Time histories of the optimal energy gain $G_n(t;t_0)$ in the base case with the disturbance initialised at different time points during the beginning stage (see inset (a)). The dash-dotted line indicates the corresponding maximum energy gain $G$ (4.1).

4.1. Algebraic growth

Figure 5 presents contours of the peak value of $G_n(t;t_0)$ during the algebraic instability as a function of the streamwise and spanwise wavenumbers in cases wherein all other parameters are the same as those in the base case. In general, we can see that with a fixed $\beta$, the energy growth during the algebraic instability diminishes with increasing $\alpha$. The maximum growth is found at $\beta = 0.68$ when $\alpha = 0$. To examine the mechanism causing the algebraic instability, we set $\alpha = 0$, such that the governing differential equations, (2.27) and (2.28), become

(4.2)$$\begin{gather} (\mathcal{D}^2 - \beta^2) \displaystyle\frac{\partial}{\partial t} \hat{v} ={-}\frac{64}{Re}\beta^2\sin{\theta}\hat{\phi}+ \frac{1}{Re}(\mathcal{D}^2-k^2)^2\hat{v}, \end{gather}$$

(4.3)$$\begin{gather}\displaystyle\frac{\partial}{\partial t} \hat{\eta} ={-}\textrm{i}\beta\mathcal{D}U\hat{v} + \frac{1}{Re}(\mathcal{D}^2-k^2)\hat{\eta}, \end{gather}$$

(4.4)$$\begin{gather}\displaystyle\frac{\partial}{\partial t} \hat{\phi} ={-} R_v\sin{\theta}\mathcal{D}\hat{\phi} - \mathcal{D}\varPhi\hat{v} + \frac{1}{PrRe}(\mathcal{D}^2-k^2)\hat{\phi}. \end{gather}$$

At this point, one can see that the algebraic instability can be attributed to two potential contributions. The first possible contribution is the buoyant forcing (the first term in the right-hand side of (4.2), which can increase the magnitude of $\hat {v}$. The second contribution is the possible growth in $\eta$ due to vortex tilting, also known as the lift-up mechanism (4.3) (Brandt Reference Brandt2014; Schmid & Brandt Reference Schmid and Brandt2014). To identify the dominant mechanism, figure 6 presents $G_n(t;t_0)$ and representative snapshots of streamlines on the $yz$-plane using $t_0 = 0.0$ in the case with $(\alpha, \beta ) = (0, {\rm \pi}/2)$, with all other parameters the same as those in the base case. In figure 6(b), we can see that the instability begins with perturbations that have a value of zero in the streamwise direction ($\tilde {u} \approx 0$), but vary in the spanwise direction, which results in a pair of streamwise vortices. As shown in figure 6(c), when the peak of $G_n(t;t_0)$ is attained during algebraic instability (see figure 6a), the transverse perturbation components (shown by the streamlines) do not change appreciably, but $\tilde {u}$ significantly increases, which, with its spanwise variation, causes growth in $\tilde {\eta }$. Therefore, as consistently reported in the context of transient behaviours in shear instabilities (e.g. Brandt Reference Brandt2014; Schmid & Brandt Reference Schmid and Brandt2014), the energy growth during algebraic instability can be attributed to the lift-up mechanism.

Figure 5. Contours of the peak values of the maximum energy gain, $G_n(t;t_0)$, during algebraic instability as a function of $\alpha$ and $\beta$ in cases where all other parameters are the same as those in the base case. The solid lines are contour lines starting at $100$ with increments of $100$.

Figure 6. (a) Time history of the development of the energy gain that yields the maximum energy growth ($t_0 = 0$) in the case with $(\alpha, \beta ) = (0, {\rm \pi}/2)$, with all other parameters the same as in the base case during the algebraic instability and the corresponding disturbance flow field at (b) $t = 0$ and (c) $= 37$, the time point at which the energy growth reaches its local maximum. In (b,c), the flow field is presented by the colour contours of the streamwise component $\tilde {u}$ and the streamlines on the $yz$-plane, for which the solid and dashed lines indicate clockwise and counterclockwise rotation, respectively.

4.2. Exponential energy growth

As shown in figure 4, $G_n(t;t_0)$ begins to increase when $t \approx 100$. The greatest exponential growth (see the pink line in figure 4) corresponds to the evolution of the perturbation initialised at $t_0 = 67.5$, which is approximately the time when the descending particle-laden interface reaches one third of the domain. That is, unlike the algebraic instability regime in which earlier $t_0$ leads to the greater and earlier growth of $G_n(t;t_0)$ (see the inset in figure 4), in the exponential growth regime, flow perturbations at the time point at which the interface reaches approximately one third of the domain produces optimal energy growth. This can be further examined from the energy equation, which is obtained by substituting (2.17) into (2.15), multiplying with $\hat {\boldsymbol {u}}$, and integrated over a control volume within a wave period, denoted as $\varOmega$. The resulting energy equation is expressed as

(4.5)\begin{equation} \displaystyle\frac{\textrm{d} }{\textrm{d} t} E ={-}\int_\varOmega \hat{u}\hat{v}\mathcal{D}U\,\textrm{d}\boldsymbol{x} - \frac{1}{Re}\int_\varOmega \left(\boldsymbol{\nabla} \times {\hat{\boldsymbol{u}}}\right)^2\,\textrm{d}\boldsymbol{x}, \end{equation}

where

(4.6)\begin{equation} E = \tfrac{1}{2}\int_\varOmega \hat{\boldsymbol{u}}\boldsymbol{\cdot} \hat{\boldsymbol{u}}\,\textrm{d} \boldsymbol{x} \end{equation}

is the kinetic energy of the disturbance. Moreover, substituting (2.17) into (2.13), multiplying with $\hat {\phi }$, and integrating over a control volume within a wave period gives

(4.7)\begin{equation} \dfrac{\textrm{d}}{\textrm{d}t} E_{\hat{\phi}} ={-}\int_\varOmega \hat{\phi}\hat{v}\mathcal{D}\varPhi\,\textrm{d}\boldsymbol{x} - \frac{1}{PrRe} \int_\varOmega \boldsymbol{\nabla}\hat{\phi}\boldsymbol{\cdot}\boldsymbol{\nabla}\hat{\phi}\,\textrm{d} \boldsymbol{x}, \end{equation}

where

(4.8)\begin{equation} E_{\hat{\phi}} = \tfrac{1}{2}\int_\varOmega \hat{\phi}^2 \,\textrm{d}\boldsymbol{x} \end{equation}

is the generalised potential energy of the disturbances (South & Hooper Reference South and Hooper1999; Sameen & Govindarajan Reference Sameen and Govindarajan2007; Jerome et al. Reference Jerome, Chomaz and Huerre2012). In both (4.5) and (4.7), the first terms on the right-hand sides represent the energy production that transfers energy from the shear of the base flow, while the second terms on the right-hand sides of these two equations represent dissipation. Equations (4.5) and (4.7) clearly show that the energy of disturbances can grow due to the interaction between the base-flow shear and disturbance fields of momentum and concentration, respectively. The exponential growth of the flow energy is much stronger than that in the algebraic case and would dominate the flow instability throughout the descent of the particle-laden interface. Figure 7 shows the contours of the peak values of the maximum energy gain ($G$) at various $\alpha$ and $\beta$ values in cases wherein all other parameters are the same as those in the base case. It can be seen that for a fixed $\alpha$, $G$ monotonically increases with decreasing $\beta$, which indicates that perturbations with a streamwise variation are most susceptible to the shear of the base flow, and the spanwise variation can stabilise the flow. Moreover, figure 7 shows that the most unstable case (i.e. greatest $G$) occurs at $\alpha \approx {\rm \pi}$, which means that the perturbation leading to the greatest amplification of energy has a wavelength which is comparable to twice the domain width, $b$, i.e. a vortex with the size of the domain width.

Figure 7. Contours of peak values of the maximum energy gain, $G$, as a function of $\alpha$ and $\beta$ in cases wherein all other parameters are the same as those in the base case.

At this point, it is also interesting to investigate the evolution of the disturbance field in response to the descent of the particle-laden interface. We examine this in a case wherein all parameters are the same as those in the base case but for the setting $(\alpha, \beta ) = ({\rm \pi} /2, 0)$. The panels in the upper row of figure 8 show five consecutive snapshots of the streamlines in the $yx$-plane at $z = 0.5$ at representative time points for the disturbance field with $t_0 = 0$. Initially, streamwise perturbations form counter-rotating pairs of spanwise vortices (i.e. $\tilde {\omega }_z$) along the $x$-direction. During the descent of the particle-laden interface, the upper panels in figure 8 show that the vortices in the upper layer (i.e. the clear-fluid layer) are distorted toward the left and those in the lower layer (i.e. the particle-laden layer) are distorted toward the right. This vortex distortion can be explained by the equation of the spanwise vorticity $\tilde {\omega }_z$, which is obtained by applying the curl operator to the momentum equations for disturbances (i.e. (A2)–(A4)). Neglecting the diffusion term, the temporal change of $\tilde {\omega }_z$ can be expressed as follows:

(4.9)\begin{equation} \frac{\partial \tilde{\omega}_z}{\partial t} + U\frac{\partial \tilde{\omega}_z}{\partial x} = \tilde{v}\mathcal{D}^2U. \end{equation}

Equation (4.9) shows that the temporal change of $\tilde {\omega }_z$ at fixed spatial points is caused by transport by the base flow and the transfer from the vorticity of the base flow through the perturbation component $\tilde {v}$, i.e. $\tilde {v}\mathcal {D}^2 U$. Here, the second derivative of the base-flow velocity is typically one order of magnitude greater than the velocity itself. Therefore, unless the streamwise wavenumber $\alpha$ is high (e.g. ${>}O(10)$), the temporal change in $\tilde {\omega }_z$ is dominated by $\tilde {v}\mathcal {D}^2 U$. The panels in the lower row of figure 8 show contours of $\tilde {v}\mathcal {D}^2 U$ that correspond to the upper panels. Because $\mathcal {D}^2 U$ has different signs in the clear-fluid and particle-laden layers (see figure 3) and $\tilde {v}$ has the same sign throughout the $y$-direction, the $\tilde {u}\mathcal {D}^2 U$ field usually jumps to a different sign across the descending particle-laden interface, thus distorting the streamline patterns along the $y$-direction.

Figure 8. Consecutive snapshots of the streamlines in the $yx$-plane (upper panels) of the disturbance flow field that leads to the maximum energy gain during exponential instability in the case wherein all of the parameters are the same as those in the base case except that $(\alpha, \beta ) = ({\rm \pi} /2, 0)$, and the corresponding contours of $\tilde {v}\mathcal {D}^2 U$ (lower panels).

5. Effect of ${R}_{v}$ and $Re$ on disturbances

As presented in § 3, with fixed $\theta$ and $Pr$ values, the combined effect of $R_v$ and $Re$ dominates the base flow, whereas the disturbance fields are dominated by $R_v$ and $Re$ separately (see (2.28)–(2.30)). In this section, we examine the energy gain as a function of $R_v$ and $Re$ in the exponential growth regime. As shown in the inset of figure 4, in the exponential regime of the present cases, once the energy begins to grow, $G_n(t; t_0)$ can easily reach $10^4$ and typically soon increases to $O(10^{6})$ in the case with the highest energy. Here, we set $G = 10^4$ as a criterion to distinguish between two regimes (i.e. $G > 10^4$ or $< 10^4$). Although, this criterion is somewhat arbitrary, we note that the overall trend that we present subsequently would not change if other values were adopted. Figure 9 presents the contours of $G = 10^4$ as a function of $\alpha$, $\beta$ and $R_v$, with fixed $Re$ (${=}215$) and $\theta$ (${=}30^\circ$ values). In the figure, the area enclosed by the contour is the regime in which $G > 10^4$, whereas outside the contour, $G < 10^4$. In figure 9, we can see that the area enclosed by the contour of $G = 10^4$ shrinks as $R_v$ increases. This enclosed area becomes zero when $R_v = 0.5$, which indicates the stabilising effect of the descent of the particle-laden interface, mainly because of the weakened base-flow field as $R_v$ increases, as discussed in § 3. The contours tend to remain unchanged as $R_v$ decreases (e.g. $R_v < 0.001$), which means that for given $Re$ and $\theta$ values, a maximum instability regime can be obtained with $R_v \rightarrow 0$.

Figure 9. Contour lines of $G = 10^4$ as a function of $\alpha$ and $\beta$ for various $R_v$ in the cases wherein all other parameters are the same as those in the base case.

Figure 10 presents the contour lines of $G = 10^4$ at $\alpha = 3.0$ as a function of $\beta$, $Re$ and $R_v$. Again, in this figure, the area beneath the contour line in each case represents the regime in which $G > 10^4$. Therefore, the smaller enclosed area with increasing $R_v$ indicates that the flow becomes more stable as particles settle faster, the same trend can also be seen in figure 9. Generally, with increasing $Re$, flows become unstable when certain threshold values of $Re$ are reached. However, in the cases with greater $R_v$ (e.g. $R_v \geqslant 0.028$), the contour lines in figure 10 imply that the flows become stabilised when $Re$ further increases (e.g. $Re > 500$). This is because during the transient stage, with the greater $Re$, it requires more time for the flow to be fully developed. In this case, the presence of the descending particle-laden interface hinders the development of the flow, such that the shear and resulting instability are suppressed.

Figure 10. Contour lines of $G = 10^4$ as a function of $\beta$ and $Re$ with a fixed $\alpha = 3.0$ for various $R_v$ in cases wherein all other parameters are the same as those in the base case.

6. Effect of the Prandtl number on disturbances

In this section, we investigate the influence of the Prandtl number, $Pr$. Figure 11 shows the time evolution of the maximum energy gain, $G(t)$, in the base case with various $Pr$. In this figure, each curve of $G(t)$ is obtained based on the collection of $G_n(t;t_0)$ values with disturbances initialised at $t_0$ ranging from $0$ to $204$ with an equal time spacing $\Delta t_0 = 10.75$. It can be seen that increasing $Pr$ results in an increase of $G(t)$. This is due to the fact that $PrRe$ controls the diffusion of mass transport, (2.30f). With a fixed $Re$, a higher $Pr$ leads to a sharper particle-laden interface in the base flow and higher energy dissipation in disturbances. Moreover, figure 11 shows that $G(t)$ converges with increasing $Pr$, which implies that when there are large particles in association with infinitesimal diffusivities, asymptotic results can be anticipated with increases in $Pr$.

Figure 11. Time histories of the maximum energy gain $G$ in cases with $PrRe$ varying from $1000$ to $9000$ and all other parameters being the same as those in the base case. Here, each curve of $G(t)$ is obtained based on the collection of $G_n(t;t_0)$ (see (4.1)) in which $t_0 \in [0, 204]$ with an increment $\Delta t_0 = 10.75$.

7. Influence of the inclined angle, ${\theta }$, on disturbances

In this section, we investigate the influence of the inclined angle, $\theta$, on the energy growth of disturbances. In figure 12, we present the temporal evolution of the maximum energy gain, $G(t)$, in cases with various $\theta$. Because $\theta$ varies between cases, the total durations of all cases are different depending on $V_0^*\sin {\theta }$. Therefore, for easier comparison, we rescale the time by multiplying by $R_v \sin {\theta }$. In this way, the inclined angle is taken into account in the new time scale, and the new dimensionless time $tR_v \sin {\theta }$ represents the time relative to particle settling, i.e. $tR_v \sin {\theta } = 1$ when the particle laden layer ($\varPhi = 1$) fully settles on to the bottom. In the figure, each curve of $G(t)$ is obtained based on the collection of $G_n(t;t_0)$ values with disturbances initialised at $t_0 R_v\sin {\theta }$, ranging from $0.0$ to $0.6$ with an equal time spacing $\Delta t_0 R_v\sin {\theta } = 0.05$. Figure 12 shows that increasing $\theta$ delays the development of energy growth, which is evident in aspects of both the algebraic and exponential growth. Moreover, the peak of $G(t)$ decreases with increases in $\theta$, except that when $\theta = 50^\circ$, no exponential growth occurs. These results are mainly due to the decrease in buoyant forcing in the momentum of the base flow when $\theta$ increases (see (2.19)), which reduces the generated velocity shear, thereby suppressing instability. Another important effect of $\theta$ is on the descending speed of the particle-laden interface. As shown in the second term on the left-hand side of (2.21), increasing $\theta$ yields higher descending speed of the particle-laden interface. As discussed in § 3, a high descending speed of the interface may restrict the base flow from reaching its quasi-steady state. As a result, the strength of the base flow decreases and the flow, thus, becomes more stable. Moreover, figure 12 shows that $G(t)$ dramatically decreases when $\theta$ becomes greater than $40^\circ$, indicating that flow instabilities are significantly suppressed. When $\theta$ further increases to $55^\circ$, it can be seen that values of $G(t)$ at $t R_v\sin {\theta } > 0.50$ are much reduced compared with all the other cases with smaller $\theta$. In fact, when $\theta = 55^\circ$, increase of $G(t)$ is only attributed to algebraic growth. This can be seen from figure 13, where we plot all curves of $G_n(t; t_0)$ that make $G(t)$ with $\theta = 55^\circ$ in figure 12. Each curve of $G_n(t; t_0)$ has either one or two peaks, in which the first peak corresponds to the algebraic growth and the second is due to the exponential growth. However, unlike the base case as shown in figure 4 where the exponential growth is much stronger, even the greatest value of the second peak of $G_n(t; t_0)$ in figure 13 (note the dash-dotted curve started at $t = 52$) is found to be much less than the first peak values of $G_n(t; t_0)$. Therefore, one may conclude that the flow can be significantly stabilised in its exponential regime when $\theta > 40^\circ$. At this point, it is worthwhile to mention that the previous numerical study by Chang et al. (Reference Chang, Chiu, Hung and Chou2019) suggested an optimal inclined angle $\theta _{op} = 40^\circ \text {--}50^\circ$ such that the sedimentation process can be accelerated to a maximum extent without causing flow instability. Moreover, an empirical value of ${\approx }45^\circ$ was suggested by the wastewater engineering practice (Kao Reference Kao1978) to reach the maximum separation efficiency. Our finding here seems to support these previous findings from different study approaches.

Figure 12. Time histories of the maximum energy growth with different $\theta$ in cases wherein all other parameters are the same as those in the base case. Here, each curve of $G(t)$ is obtained based on the collection of $G_n(t;t_0)$ (see (4.1)) in which $t_0 R_v\sin {\theta } \in [0, 0.6]$ with an increment $\Delta t_0 R_v\sin {\theta } = 0.05$.

Figure 13. Time histories of all $G_n(t; t_0)$ that make the curve of $G(t)$ with $\theta = 55^\circ$ in figure 12. The dash-dotted line denotes the curve associated with $t_0 = 52$, which shows two peaks of $G_n$. The thick grey line represents the curve of $G$ obtained based on the collection of the local maximum of all $G_n$ curves (4.1).

8. Applicability

It should be noted that in this study, an important assumption we have made is the small particle size, such that particles can soon reach their terminal speed and our scalar transport equation used to describe the volume fraction of particles is valid. This is the case when the particle relaxation time, $\tau _p$, is small. For spherical particles, this means that $\tau _p^* = \rho _s^*d_0^{*2}/(18\mu ^*) \ll 1$, in which $d^*_0$ is the particle diameter. In liquid flows for which the $\mu \approx 10^{-3}\ \textrm {kg}\ \textrm {m}^{-1}\ \textrm {s}^{-1}$ (e.g. water) or greater, this can be easily met if the particle diameter $d_0 < O(1\ \mbox {mm})$, which is very common in suspension flows in various engineering and industrial applications (Galvin et al. Reference Galvin, Callen, Zhou and Doroodchi2005; Madge, Romero & Strand Reference Madge, Romero and Strand2005). Another assumption is that the particle concentration is dilute such that the viscosity cannot be significantly altered in the presence of suspended particles (see § 2). According to the theoretical work of Batchelor & Green (Reference Batchelor and Green1972), the flow viscosity can be significantly altered when the volume fraction $\phi ^*$ is $O(0.1)$. Therefore, the applicability of the present findings can be extended to suspension liquids in which $\phi ^* < 0.1$. However, even when $\phi ^*$ is beyond this criterion in the highly concentrated case, the procedure presented here can be modified by the inclusion of the existing models that describe the flow viscosity as empirical functions of $\phi ^*$.

The major contribution of the present study is the investigation of the case in which particles settle fast, such that the unsteadiness of the base flow is not negligible. According to our findings, this is the case when the parameter $R_v = V_0^*/U^*_c > 0.001$ (see figures 9 and 10). To see this more clearly, using the formula of $V_0^*$ based on Stokes’ law, i.e.

(8.1)\begin{equation} V_0^* = \dfrac{d_0^{*2}}{18\mu^*}\left(\rho^*_s-\rho^*_0\right)g^*, \end{equation}

and the definition of $U_*$ (2.12), $R_v$ can be reduced to

(8.2)\begin{equation} R_v = \dfrac{32}{9}\dfrac{d_0^{*2}}{c_0^*b^{*2}}. \end{equation}

That is, $R_v$ appears to be a parameter indicating the competitive influences among the width of the liquid column ($b^*$), particle size ($d_0^*$) and volume fraction ($c_0^*$). As seen from (8.2), $R_v > 0.001$ requires that $d^*_0 > 0.00028 b^* \sqrt {c_0^*}$, which can be easily met in realistic conditions. For example, in the laboratory, the magnitude of $b^*$ is typically $O(1\ \mbox {cm})$. Considering a case in which $c_0^* \sim O(0.01)$, $d^*_0$ needs to be at least $O(1\ \mathrm {\mu }\mbox {m})$ to satisfy $R_v > 0.001$, which is very common for suspended materials. In such cases, the findings of the present study can provide useful guidance in optimising the design for the separation process of suspension liquid.

9. Concluding remarks

In conclusion, we have analysed the stability of particle-laden flows in long, tilted water columns in batch settling mode beyond the quasi-steady approximation, which is essential for the fast settling of particles in the problem of sedimentation. By introducing a settling time scale in the momentum and transport equations, the transient behaviour of base flow was resolved in terms of Reynolds number ($Re$), Prandtl number ($Pr$) and the settling speed ratio ($R_v$). The magnitude of base flow increases as the settling speed is reduced, which attains its maximum value when the settling speed becomes infinitesimal. Based on the unsteady base-flow solutions, the non-modal analysis was utilised to analyse the temporal evolution of the disturbance flow field. The time history of the disturbance flow energy gain experiences an algebraic and an exponential growth stage owing to the lift-up mechanism of wall-normal disturbance and the shear instability of base flow, respectively. The flow instability is enhanced by the increase of Prandtl number owing to the sharpening of the particle-laden interface, and suppressed by the increase of settling speed because of less disturbance energy being extracted from the base flow. Considering the efficiency of sedimentation, for example, in the removal of suspended particles, there exists an optimal inclined angle $\theta \approx 50^\circ$, which provides the maximum capacity for sedimentation enhancement while maintaining flow stability to prevent the resuspension of the sedimenting particles.