Instability of a dusty shear flow

Particle-laden flows are prevalent in various environmental, astrophysical, and industrial settings. In the environmental context, these flows play pivotal roles in shaping the landscape around us through sediment transport (burns2015sediment), influencing weather patterns via clouds (pruppacher1998microphysics; shaw2003particle), sea-spray (veron2015ocean), volcanic eruptions (bercovici2010two), and gravity currents (necker2005mixing). Astrophysical scenarios include particle-laden flows in cosmic dust clouds, stellar winds, and the transport of interstellar dust particles in accretion disks and proto-planetary disks, which are crucial for forming and evolving planetary systems (youdin2005streaming; fu2014effects; homann2016effect). In industrial and engineering applications, particle-laden flows arise in processes like spray drying (straatsma1999spray; birchal2006spray), powder handling (baxter2000simulation), and fluidized bed reactors (bi2000state). Particle-laden flows typically involve multiple components, with at least a carrier phase and a dispersed phase, leading to physics occurring at multiple scales. In natural scenarios, the carrier flows are often turbulent, and the suspended particles are advected and sheared by this turbulence. Modelling particle-laden flows within a sufficiently dilute limit often involves considering the momentum exchange from the fluid to the particles while neglecting feedback from particles to the fluid (one-way coupling). However, this feedback is significant, especially when the mass fraction of particles to fluid is on the order of unity (e.g., dusty gas flows or water droplets in the air), as it can introduce new dynamics into the system. In this study, we investigate a novel instability in a particle-laden simple shear flow arising from two-way coupling, where the feedback force from particles to fluid is considered.

The non-uniform distribution of particles, also known as particle segregation or particle banding, is observed in particle-laden flows across diverse engineering and environmental contexts. It can occur due to various segregation mechanisms such as preferential clustering, differential settling velocities, turbulent dispersion, and particle-particle interactions. For example, in turbulent flows, inertial particles can cluster preferentially in regions of high strain and lower vorticity, forming regions with high and low particle concentration (maxey1987gravitational; bec2005multifractal; fiabane2012clustering). This phenomenon can arise solely from one-way coupling, and no feedback force is required. In pneumatic conveying systems transporting granular materials, such as powders or grains, bands of particle-rich and particle-deficient regions can form (known as rope formation) due to agglomeration and flow dynamics (klinzing2011pneumatic; huber1994characterization; lain2013characterisation; zhang2023coarse). In turbulent wall-bounded flows carrying suspended particles, such as industrial pipelines transporting slurries, particle segregation can occur due to differential settling velocities, turbulent dispersion and thermophoresis. This segregation can lead to the formation of particle-rich layers near the bottom or walls of the flow channel, with particle-deficient regions in the core of the flow (marchioli2003direct; sardina2012wall; picano2009spatial). The particles form extremely long clusters, called ropes, and align preferentially with the low-speed turbulent streaks, contributing to their stabilization and suppression of bursting (dave2023mechanisms). Despite the additional stresses resulting from particles, the alteration of near-wall coherent structures results in a notable decrease in Reynolds shear stresses and partial relaminarization of the near-wall flow. Fluidized bed reactors used in chemical engineering often exhibit particle banding phenomena (harris1994solitons; liu2016cfd; gilbertson2001segregation). Observations of sediment-laden river flows and estuarine environments have revealed the formation of bands of sediment deposition influenced by flow dynamics, sediment transport mechanisms, coagulation and channel morphology (gibbs1986segregation; sondi1995sedimentation; ogami2015dynamic). These sediment bands play a crucial role in shaping riverbeds, deltas, and coastal environments.

In astrophysical accretion disks, such as those around young stars or black holes, there are regions where gas and dust particles orbit around a central object. The dust-gas system exhibits Keplerian motion, alongside radial and azimuthal drifts between the dust and gas. When the gas and dust move at slightly different velocities, this disparity can result in a relative drift between the two components. This relative motion creates a shearing force that can amplify small perturbations/disturbances in the dust distribution - known as streaming instability (youdin2005streaming; chiang2010forming). The name “streaming instability” reflects the differential streaming motion between gas and dust particles within the disk. The instability arises from the relative drift between the gas and dust phases, a universal consequence of radial pressure gradients. Growth occurs even though the two components interact only via dissipative drag forces. Streaming instability exhibits growth rates that are slower than dynamical time scales but faster than drift time scales. As a result, dust particles start clumping together, even without self-gravity, forming dense structures or bands. Thus, the dust particles can localise, leading to additional dynamics or instabilities due to the Keplerian shear (as we see in this paper) and self gravity. Thus, streaming instability is crucial in various astrophysical processes, including planetesimals’ formation and dust grains’ growth in protoplanetary disks.

The hydrodynamic stability characteristics of particle-laden flows can be altered by modifying existing instabilities or generating new types of instabilities, as one considers the feedback from particles. Early studies by kazakevich1958investigations; sproull1961viscosity observed a notable reduction in the resistance coefficient when dust was added to the air flowing turbulently through a pipe. It was thought that adding particles altering the effective viscosity led to this modification; however, this contradicts Einstein’s formula for the suspension viscosity. saffman1962stability was among the first to provide an analytical model for studying the stability of a dusty planar flow. He proposed that inertial particles extract energy from turbulent fluctuations, thereby damping them. He modelled both particle and fluid phases as a continuum using a two-fluid model. The momentum exchange between both phases is accounted for using a linear Stokes drag. The modal analysis ultimately resulted in a modified Orr-Sommerfeld equation for uniformly distributed particles. saffman1962stability deduced that finer particles with low inertia could induce destabilization due to a reduction in effective kinematic viscosity. Although the effect of particles on the viscosity of dusty gas is negligible, it effectively increases the gas density, thereby reducing the kinematic viscosity. Conversely, coarser particles with high inertia could lead to stabilization through dissipation by Stokes drag. He concluded that dust merely alters the waves present in a clean gas and may not introduce any additional instabilities. Subsequent studies have numerically solved the modified Orr-Sommerfeld equation for various base flows, confirming Saffman’s conclusions (michael1964stability; asmolov1998stability; tong1999two; klinkenberg2011modal; sozza2022instability). sozza2020drag; sozza2022instability investigated a dusty Kolmogorov flow and demonstrated that increasing the particle mass loading reduces the amplitude of the mean flow and turbulence intensity. They observed that turbulence suppression is more pronounced for particles with smaller inertia. The study concluded that while inertia significantly influences particle dynamics, its impact on flow properties is negligible compared to mass loading.

Notably, Saffman’s analysis does not account for gravitational effects. Including gravity can introduce buoyancy effects that may destabilize the flow more easily (herbolzheimer1983stability; shaqfeh1986effects; borhan1988sedimentation). For example, magnani2021inertial investigated dusty Rayleigh-Taylor turbulence and found that the interface between the two phases becomes unstable in the presence of gravity forces, evolving into a turbulent mixing layer that broadens over time. However, in a particle-laden Rayleigh–Bénard system (prakhar2021linear), it was observed that particles tend to inhibit the onset of natural convection, thereby stabilizing the system. This is because particles act as a distributed drag force and heat source in the fluid, similar to a porous medium. Saffman’s analysis, also, focusing solely on uniform particle distribution, overlooks the potential effects of non-uniform distributions. A study by narayanan2002temporal has demonstrated the presence of additional unstable modes under large Stokes numbers and high mass loading conditions in a particle-laden mixing layer where the particle distribution is localised. Additionally, investigations utilising non-uniform particle distributions have highlighted the emergence of novel instabilities (senatore2015effect; warrier2023stability).

The non-uniform distribution of inertial particles arises naturally in vortical flows due to their preferential accumulation. As noted earlier, it has long been recognized that inertial particles tend to be centrifuged from vortical regions and cluster in regions of high strain (maxey1987gravitational), a phenomenon known as preferential accumulation. A numerical investigation by shuai2022accelerated demonstrated that in a Lamb-Oseen vortex, inertial particles are expelled from the vortex core, forming a ring-shaped cluster and a void fraction bubble that expands outward. However, without accounting for the two-way interaction, the vortex would decay slowly due to viscous dissipation. When the two-way coupling is considered, it is observed that the feedback from clustered particles flattens the vorticity distribution and leads to an accelerated vortex decay. It is noted that as the inertia of the particles increases, the vorticity decays even more rapidly. A follow-up numerical study by shuai2022instability on a particle-laden Rankine vortex revealed that the system becomes unstable due to the two-way interaction, which is also validated by analytical linear stability analysis. The feedback force from the particles triggers a novel instability, which can cause the breakdown of an otherwise resilient vortical structure. In the context of the merging of a pair of co-rotating vortices laden with inertial particles, it has been shown that (shuai2024merger) the feedback force from particles significantly alters the monotonic merging behaviour as observed without feedback. The vortices push apart for a while due to a net repulsive force from particles ejected from the vortex cores, thus delaying their merging.

Apart from modal instability, the interplay between particle inertia and shear from the base flow can lead to transient growth of perturbations in particle-laden flows via non-modal growth mechanisms such as the Orr mechanism or the ‘lift-up’ effect. Performing a non-modal analysis, klinkenberg2011modal showed that transient growth in a particle-laden channel flow increases proportionally with the particle mass fraction. Similar non-modal instabilities have been observed and studied in dusty gas flows, such as the Blasius boundary layer with a localized dust layer (boronin2014modal), plane channel suspension flow with a Gaussian layer of particles (boronin2016nonmodal), and stably stratified Blasius boundary layer flow (parente2020modal). Additionally, the inclusion of gravitational effects in a simple shear flow has been shown to alter the uniformity of particle distribution, leading to the formation of local particle clusters and, subsequently, to transient growth (kasbaoui2015preferential).

The previously mentioned studies mostly used an Eulerian-Eulerian model for particle-laden flows. However, there are three primary methods for modelling particle-laden flows: (i) Eulerian-Eulerian modelling, (ii) Eulerian-Lagrangian modelling, and (iii) fully resolved simulations. The Eulerian-Eulerian method (see jackson2000dynamics; drew2006theory) treats both the particle and fluid phases as interpenetrating continua in an Eulerian framework, assuming particles are small ($\Delta x\gg d_{p}$) and sufficiently densely distributed to be treated as a fluid continuum. While computationally more affordable, this method may introduce errors due to the difficulty of modelling terms in Eulerian-Eulerian methods that require closure. In contrast, the Eulerian-Lagrangian method treats the fluid phase as a continuum within an Eulerian framework while representing the particle phase as discrete entities tracked individually in a Lagrangian manner. Here, the fluid phase is solved at a coarser scale, typically with $\Delta x$ a few times larger than $d_{p}$, resulting in a scalable and cost-effective approach. However, empirical coupling between the particle and fluid phases introduces some approximation errors. Few studies that considered Eulerian-Lagrangian modelling and stability of particle-laden flows are meiburg2000vorticity; richter2013momentum; senatore2015effect; kasbaoui2015preferential; wang2019modulation; kasbaoui2019turbulence; pandey2019clustering; shuai2022instability. The fully resolved method, for example, realised by the immersed boundary method (see uhlmann2005immersed; breugem2012second; kempe2012improved; dave2023volume; Kasbaoui2024high-fidelity), is employed when the grid spacing $\Delta x$ is significantly smaller than the particle diameter $d_{p}$. While offering high accuracy, this method requires extensive resolution, making it costly and impractical for large-scale applications. Thus, in this study, we only employ the Eulerian-Lagrangian and Eulerian-Eulerian methods, the details of which are discussed respectively in §2.1 and §4.1.

Here, we investigate the stability of a dusty simple shear flow with non-uniformly distributed particles. In the absence of particles, a simple shear flow is modally stable to infinitesimal perturbations. Similarly, a non-uniform particle distribution without any background flow remains unaffected. Thus, each system under consideration is linearly stable when inspected in isolation. However, when considering the combined system (see the schematic in figure 1), where a simple shear flow is superimposed on a band of particles, the background flow advects the particles, and the feedback from the particles induces an instability in the system. This paper reveals that this seemingly simple yet crucial particle-laden system exhibits a novel type of instability when incorporating two-way coupling. In §2, we demonstrate the presence of this novel instability through Eulerian-Lagrangian numerical simulations. The mechanism underlying the genesis of this new instability is described in §3. An analytical study of the system is carried out in §4 using an Eulerian-Eulerian framework. Following the approach of saffman1962stability, linear stability analysis is employed to establish the existence and modal nature of the instability. §5 offers a comparative analysis between the Eulerian-Lagrangian and Eulerian-Eulerian results. Finally, we conclude in §6.

In this section, we present a novel instability occurring in a particle-laden simple shear flow. Along with the momentum transport from fluid to particle phase, we consider the feedback from the particle to fluid phase (two-way coupling), which is crucial for the instability to occur. We employ an Eulerian-Lagrangian (EL) method to illustrate this instability. Below, we describe the method used:

In the previous section, we observed that instability arises from the two-way coupling between the particle and fluid phases, driven by the finite inertia of the particles. Surprisingly, even weakly inertial particles ($St=10^{-3}$) triggered the instability. In this section, we delve into the instability mechanism in the small particle inertia limit while still considering the particle-fluid coupling. We employ the concept of wave interaction to elucidate this instability mechanism. In the weak particle inertia limit, we will demonstrate that our particle-laden system resembles a stratified fluid system. The fluid exhibits an effective modified density, which is stratified due to the particle concentration gradient. This density stratification creates edge waves, and their interaction leads to instability, similar to the case of the Taylor-Caulfield instability (taylor1931effect; caulfield1994multiple; balmforth2012dynamics). However, there is a caveat: the wave generation mechanism differs slightly here, which we will address as we proceed.

In the limit of weak particle inertia, following ferry2001fast; rani2003evaluation, the feedback force $\textbf{F}_{p}$ in equation (1) for dusty-gas system can be approximated as $\textbf{F}_{p}=-\rho_{p}\,\phi_{p}\,\textbf{a}+\textit{O}(\tau_{p})$, where $\textbf{a}=\left(\partial\textbf{u}/\partial t+\textbf{u}\cdot\boldsymbol{\nabla}\textbf{u}\right)$ represents the fluid acceleration term. Substituting it back into equations (2.1b) yields a simplified form representing a fluid with modified density (in the inviscid limit) as

where $\rho=\rho_{f}+\rho_{p}\,\phi_{p}$ represents an effective fluid density due to the suspended particles. Thus, a spatial inhomogeneity in the particle concentration ($\phi_{p}$) can result in a variation in the effective density of this composite fluid even if $\rho_{f}$ is a constant. Additionally, the conservation of the number of particles yields

Together, equations (2.1a), (3) and (4) resemble a stratified incompressible fluid system and is known as the single-fluid continuum model describing a particle-laden system. Taking the curl of equation (3) after scaling it with $\rho$ gives the evolution equation for the flow vorticity ($\textbf{q}=\boldsymbol{\nabla}\times\textbf{u}$) as (for a two-dimensional case)

where, we have used the relation between a and $\boldsymbol{\nabla}p$ from equation (3) here for further simplification. Equation (5) indicates that vorticity can arise in the system due to the misalignment between fluid acceleration and the gradient in particle concentration due to the baroclinic source term $\textbf{a}\times\boldsymbol{\nabla}\rho$. In the subsequent paragraphs, we will demonstrate the precise way by which vorticity is generated and how it gives rise to propagating waves that can interact and lead to instability.

The base state flow may inherently possess a vorticity field. However, the additional vorticity disturbance created would be responsible for generating waves. To analyze this, we need to consider the evolution equation for the disturbance vorticity field. Without loss of generality, let’s consider a general two-dimensional system in the $x$–$y$ plane with uni-directional flow in the horizontal direction and vertically stratified particle distribution. We can decompose the relevant quantities into base state and disturbance (denoted by $\tilde{()}$) parts as $\textbf{u}=U(y)\,\hat{\textbf{x}}+\tilde{\textbf{u}}(x,y)$ and $\rho=\rho_{b}(y)+\tilde{\rho}(x,y)$. The vorticity field $\textbf{q}=q_{z}\,\hat{\textbf{z}}$ will be along the $\hat{\textbf{z}}$ direction and can be decomposed as $q_{z}=Q_{z}(y)+\tilde{q}_{z}(x,y)$, where the base state vorticity is $Q_{z}(y)=-U^{\prime}(y)$. Substituting these expressions into equation (5) and considering the leading disturbance terms, we obtain the linearized evolution equation for the disturbance vorticity field as

where the operator $\mathcal{D}/\mathcal{D}t=\partial/\partial t+U(y)\,\partial/\partial x$ represents the linearized material derivative, and $()^{\prime}$ represents the operation $d/dy$. Utilizing the relationship between the vertical displacement field $\tilde{\eta}$ and the vertical velocity field $\tilde{u}_{y}=\mathcal{D}\tilde{\eta}/\mathcal{D}t$, we can rewrite equation (6) as

i.e., the quantity inside the bracket is materially conserved. In other words, the disturbance vorticity $\tilde{q}_{z}$ is related to the horizontal disturbance velocity $\tilde{u}_{x}$ and the vertical displacement $\tilde{\eta}$ as

The constant term represents a bulk vorticity in the background flow. Since we are primarily interested in the relative vorticity generation $\tilde{q}_{z}$ with respect to this bulk vorticity, we can set the constant term to zero without loss of generality. From this general formulation, let us consider our special case where the base state flow is a simple shear flow, i.e., $U=\Gamma\,y$, and the base state particle number density has a top-hat distribution. Without loss of generality, we assume $\Gamma>0$. Since the particle concentration has sharp variations at locations $y=h/2$ and $y=-h/2$, the effective density $\rho$ of the fluid also varies sharply at these locations. The base state density $\rho_{b}(y)$ takes a constant value $\rho_{1}=\rho_{f}$ outside the particle band and $\rho_{2}=(\rho_{f}+\langle\phi_{p}\rangle\,\rho_{p})>\rho_{1}$ within the particle band, as shown in schematic figure 4. These density jumps cause the locations of the jumps ($y=\pm h/2$) to act like interfaces separating different density fluids, marked as I and II in the figure 4. To begin with, we consider each of these interfaces in isolation and demonstrate that they support propagating waves due to disturbances.

Consider the interface-I at $y=h/2$, where the effective density drops from $\rho_{2}$ to $\rho_{1}$ (as we move along the positive $y$ direction). As a result, $\rho_{b}^{\prime}(y)=-\Delta\rho\,\delta(y-h/2)<0$, where $\Delta\rho=(\rho_{2}-\rho_{1})>0$ and $\delta(\cdot)$ represents a Dirac delta function. Let’s assume the interface is perturbed sinusoidally ($\tilde{\eta}$) as shown in the figure 4. From equation (8), we can deduce that a disturbance vorticity field is generated, given by $\tilde{q}_{z}=(\tilde{u}_{x}+\Gamma\,\tilde{\eta})\,\rho_{b}^{\prime}/\rho_{b}$. Generally, the associated horizontal perturbation velocity $\tilde{u}_{x}$ will have a discontinuity at the interface, which changes sign once we cross the interface. Physically, for a wave to be supported at the interface, there can be no self-induced $\tilde{u}_{x}$ for a wave-like solution at the interface. Thus, $\tilde{u}_{x}=0$ at the interface. Then the generated vorticity disturbance is directly proportional to the interface perturbation as $\tilde{q}_{z}=\Gamma\,\tilde{\eta}\,\rho_{b}^{\prime}/\rho_{b}$. Since $\Gamma>0$, $\rho_{b}>0$, and as we saw earlier $\rho_{b}^{\prime}<0$ at the interface-I, the generated vorticity disturbance will be out of phase with the interface displacement field. The maximum $\tilde{q}_{z}$ (counter-clockwise) occurs at the troughs of the $\tilde{\eta}$ field, and the minimum $\tilde{q}_{z}$ (clockwise) occurs at the crests of the interface perturbation, as can be seen in the figure 4 (the crests and troughs of the generated vorticity disturbances are shown as circular arrows). Since the density gradient is localised at the interface, the resulting vorticity disturbance would also be localised at the interface as a Dirac delta function with sinusoidal variation in the flow direction. The generated disturbance vorticity field induces a vertical flow velocity field $\tilde{u}_{y}$, which has a maximum at the positive sloping node and minimum at the negative sloping node of the $\tilde{\eta}$.The vertical velocity field lags behind the interface displacement by $90^{\circ}$. The crests and troughs of the $\tilde{u}_{y}$ field are represented as upward and downward vertical arrows in the figure 4. The $90^{\circ}$ phase difference serves as the perfect criterion for the propagation of the interface displacement as a wave. Upon visual inspection of figure 4, one can intuitively deduce that the upward $\tilde{u}_{y}$ at the positively sloping $\tilde{\eta}$ node and similarly the downward $\tilde{u}_{y}$ at the negatively sloping $\tilde{\eta}$ node would result in an intrinsic propagation of the wave leftward relative to the base state flow at the interface location. Also, the $90^{\circ}$ phase difference ensures that the wave would not experience any amplification or damping but only undergo propagation (see carpenter2011instability).

This is similar to the generation of a propagating edge wave in a background shear flow where there is a jump in vorticity/shear rate present (see valis2006atmospheric). A generalised edge wave in a stratified fluid, which has both jumps in shear rate (from $\Gamma_{2}$ to $\Gamma_{1}$) and density (from $\rho_{2}$ to $\rho_{1}$) across an interface, should have a phase speed $c$ of the form

Here, we don’t have a jump in shear rate (i.e., $\Gamma_{1}=\Gamma_{2}=\Gamma$), but only a jump in density. Thus, the resulting wave at the interface-I would propagate with a phase speed of $c_{\textrm{I}}=\Gamma\,h/2-\Gamma\,At/k$, where $k$ is the spatial wavenumber in the $x$ direction of the disturbance, and we used the definition of the Atwood number $At=(\rho_{2}-\rho_{1})/(\rho_{2}+\rho_{1})$, a non-dimensional number which frequently appears in the instability of density stratified flows. In non-dimensional terms (as one chooses the length scale to be $h$ and the time scale to be $\Gamma^{-1}$), the phase speed is $c_{\textrm{I}}^{*}=c_{\textrm{I}}/(\Gamma\,h)=1/2-At/k^{*}$, or the dispersion relation $c_{\textrm{I}}^{*}\,k^{*}=k^{*}/2-At$, which we will formally obtain as a large $k^{*}=k\,h$ (i.e., large separation $h$) asymptotic expression of the dispersion relation of the full system in §4.2. Here, $()^{*}$ represents non-dimensional quantities.

Similarly, another propagating wave would be generated at interface-II in isolation. At this interface, since the particle concentration varies such that $\rho_{b}^{\prime}=\Delta\rho\,\delta(y+h/2)>0$, the vorticity disturbance generated would be localised at the interface and in phase with $\tilde{\eta}$. As a result, the vertical velocity disturbance field thus produced would lead $\tilde{\eta}$ by $90^{\circ}$, and thus, the intrinsic propagation of the wave would be rightward relative to the base state flow at the interface location, as shown in figure 4. The propagation speed would be $c_{\textrm{II}}=-\Gamma\,h/2+\Gamma\,At/k$ or $c_{\textrm{II}}^{*}=-1/2+At/k^{*}$.

Now that we have established that there would be a leftward propagating wave at interface-I and a rightward propagating wave at interface-II when they are in isolation or equivalently when the interfaces are far separated (i.e., $h\rightarrow\infty$ or $k^{*}\gg 1$), let us consider the scenario as the separation $h$ between the interfaces is reduced. One can expect the interfacial waves to perceive other’s presence and interact as they come closer. This is because, although the interfacial vorticity fields are localised Dirac delta functions, the vertical velocity fields typically exhibit exponential decay (wave evanescence) away from the interface location (e.g., $\exp\left({-k\lvert y-h/2\rvert}\right)$ as we see in §4.2 as the eigenfunction associated with $\tilde{u}_{y}$). As the waves approach each other, they could interact. Still, for modal instability to occur, two essential conditions must be satisfied (see carpenter2011instability): (i) the phase speeds of the waves should be stationary with respect to each other (‘phase-locking’) and (ii) the relative phase of the waves should lead to the mutual growth of the interface disturbances. As the interfaces approach each other, the individual wave speeds also approach one another and are expected to become equal (and equal to zero, i.e., $c_{\textrm{I}}=c_{\textrm{II}}=0$) for $k\,h=k^{*}=2\,At$. However, along with this natural convergence, interactions play a role such that the waves become stationary relative to each other at a slightly larger separation (or non-dimensional wavenumber $k^{*}_{\textrm{cutoff}}>2\,At$). The leftward propagating wave at interface-I can slow down the rightward propagating wave at interface-II and vice versa (for $\Gamma>0$), thus achieving phase-locking at a larger separation. Once the phase-locking is achieved, it can be maintained even though the isolated phase speeds are generally unequal (except at $k^{*}=2\,At$) due to mutual interaction. Adjustments in the phase difference between the waves maintain an induced phase speed that exactly cancels the intrinsic propagation speed. Therefore, condition (i) would be satisfied for $0<k^{*}<k^{*}_{\textrm{cutoff}}$ through adjustments in the relative phase difference of the waves.

Similar to the propagation criterion, upon visual inspection of figure 4, one can deduce that an upward $\tilde{u}_{y}$ component at the crests of $\tilde{\eta}$ and equivalently a downward $\tilde{u}_{y}$ component at the troughs of $\tilde{\eta}$ will lead to the growth of the wave. As the isolated waves become phase-locked, the interaction occurs in such a way that the $\tilde{u}_{y}$ field of one wave at an interface and the $\tilde{\eta}$ field of the other interface satisfy this criterion, and vice versa. Thus, by condition (ii), the phase-locked waves mutually trigger amplification and arrest propagation, leading to instability through wave interaction for any smaller separation of the interfaces (or non-dimensional wavenumbers in the range $0<k^{*}<k^{*}_{\textrm{cutoff}}$). This intricate interaction between the phase-locked waves is fundamental to understanding the mechanism behind the generation of modal instability in such dust-stratified fluid systems.

We will now address the earlier caveat, highlighting an important distinction from the classical Taylor-Caulfield instability. In the Taylor-Caulfield instability, interfacial waves are generated due to a Boussinesq-type baroclinic torque driven by buoyancy/gravity effects. This occurs in a system with an overall stable stratification ($\rho_{1}<\rho_{2}<\rho_{3}$), where $\rho_{1}$ is at the top, $\rho_{2}$ is at the middle and $\rho_{3}$ is at the bottom layers. A pair of stable surface gravity waves are generated at each interface, separating these layers due to the Boussinesq baroclinic effects. One of them propagates leftward while the other propagates rightward; together, they interact to produce a stationary wave which is dampening at each interface. However, when the interfaces are brought closer, in a background shear flow ($\Gamma>0$), the interaction between the leftward propagating wave at the upper interface and the rightward propagating wave at the lower interface leads to instability.

However, in our system, no gravitational/buoyancy effects are present. Therefore, no restoring mechanism could generate propagating waves from Boussinesq baroclinic effects. Instead, the oft-neglected non-Boussinesq baroclinic effects are responsible for wave generation, acting as a different type of restoring mechanism. The waves generated closely resemble edge waves or vorticity waves rather than surface gravity waves, as only one stable propagating wave is generated at each interface. Similar to the Taylor-Caulfield instability, instability in our system requires a leftward propagating wave at the top interface and a rightward propagating wave at the bottom interface (noting that this preference is due to $\Gamma>0$ and would be reversed if $\Gamma<0$). Thus, for instability to occur, $\rho_{3}$ must be smaller than $\rho_{2}$, which we have assumed to be the same as $\rho_{1}$. However, if the system were configured with $\rho_{3}=\rho_{1}>\rho_{2}$, the direction of waves at each interface would be reversed. In such a scenario, the interaction would only amplify their propagation speed, suppressing phase-locking and avoiding any chance of instability. Similarly, if $\rho_{3}>\rho_{2}>\rho_{1}$ ($\rho_{3}<\rho_{2}<\rho_{1}$), then there would be leftward (rightward) propagation of waves at both the interfaces, whose interaction also can not produce instability. With the evidence of instability from EL simulations and an explanation of the mechanism involving interacting edge waves, we now proceed towards a quantitative evaluation of the growth rate of the instability using a linear stability analysis.

In this section, we employ an analytical approach, complemented by numerical assistance, to investigate the system more formally, illustrating the existence and criteria and quantifying the growth rate of the instability. To achieve this, we incorporate a continuum description of the particle phase, outlined below.

In the previous sections, we demonstrated the existence and mechanism of the instability using a linearized continuum model for the system. What remains is a quantitative comparison of the predictions from linear stability analysis with fully non-linear EL simulations.

We conducted a series of EL simulations across various parameter ranges to obtain the dispersion curve, maintaining most settings as described in §2.2. In these simulations, the Stokes number ($St$) was kept constant at $10^{-3}$. To ensure a fair comparison with the inviscid LSA, the EL simulations were carried out with a fixed $Re_{L_{x}}=5000$, which is sufficiently high to make viscous effects negligible. The box domain size was set as $L_{x}=2\,\pi\,\textrm{m}$, $L_{y}=6\,\pi\,\textrm{m}$ and $L_{z}=3\,d_{p}$. The average volume fraction $\langle\phi_{p}\rangle$ of particles within the band was adjusted to achieve various mass loading values. To confine only one wave within the domain in the $x$ direction, we fixed the dimensional wavenumber to $k=1\,\textrm{m}^{-1}$. By adjusting the bandwidth $h$, we obtained different non-dimensional $k$ values. In all selected cases, the ratio $h/d_{p}$ was set to be large (on the order of $10^{2}$ to $10^{3}$), ensuring that fluctuations caused by discrete particle forcing remained well below the viscous dissipation scale. This configuration enabled fluid-particle coupling primarily through collective particle dynamics at scales comparable to $h$, rather than discrete effects at inter-particle distances. Since $L_{x}$ is fixed, varying $h$ leads to a change in the ratio $L_{x}/h$ across the range from $2.5\,\pi$ to $32\,\pi$. However, $L_{y}$ remains sufficiently large in all these cases to avoid any possible interactions with other periodic simulation boxes in the $y$ direction.

The initial conditions of the numerical simulations consist of a combination of the base state (a simple shear flow with a top-hat number density profile) and the perturbation eigenmodes obtained from the corresponding LSA. To perturb our system, we utilize the analytic expressions for the eigenmodes of the velocity field obtained for the corresponding $St=0$ case. The initial particle number density field is left unperturbed since the eigenmodes are theoretically singular Dirac delta functions. Although this initial arrangement may not correspond precisely to an eigenmode of the system under consideration (due to the omission of finite $St$ aspects and perturbations in the number density field), it will closely resemble the eigenmode. As time progresses, the applied perturbations are expected to excite the most unstable mode, i.e., the respective eigenmode of the system. The initial perturbation fields used for the flow field can be visualized in the representative figure 10(a, f, k), and the corresponding initial particle distribution can be seen in figure 11(a, f, k). To ensure that the simulations capture the linear regime described by LSA and to delay any nonlinear effects, the perturbations are initialized with a small relative amplitude $\epsilon=10^{-2}$. The location of the particles are randomly initialised within each cell inside the particle band, defined by $\lvert y\rvert\leq h/2$.

Figure 10 illustrates the time evolution of the perturbation vorticity field (normalised by its initial maximum value) for a case with $M=1$ and $St=10^{-3}$, depicted for three different non-dimensional wavenumbers: $k=0.25$ (top panel), $k=0.5$ (middle panel), and $k=1.2$ (bottom panel). As earlier, the plots are presented in a dimensional spatial domain of $2\pi\times 2\pi$ for better visualisation. The variation in wavenumbers is noticeable from the variation in $h$ in the initial disturbance fields (panels a,f,k). The initial evolution of the perturbation fields confirms the excitation of the respective eigenmodes. For the smaller wavenumber (top panel), it can be seen that the perturbations grow over time, and the vorticity field amplifies and spreads across a larger region, indicating the presence of instability. At $k=0.5$ (middle panel), the instability becomes more pronounced and transitions to the nonlinear regime rapidly, as expected theoretically, since the maximum growth corresponds to a non-dimensional wavenumber closer to $k=0.5$. The theory predicts no linear instability for the larger wavenumber (bottom panel), which is also evident from the snapshots where there is minimal growth of the vorticity field. Figure 11 shows the corresponding evolution of the total particle number density field. It can be observed that although the number density field is initially not perturbed, eventually, the eigenmodes are getting excited. Additionally, it can be noticed that the evolution of the number density field is coherent with the corresponding perturbation vorticity field at all times for all wavenumbers. This correspondence also holds theoretically for a top-hat profile, as we have already seen using LSA that the vorticity perturbation and number density perturbation fields both are Dirac delta functions. The snapshots show that the case with wave number $k=0.5$ transitions to the nonlinear regime much faster than the others.

In order to determine the total perturbation growth rates from the present simulations, we calculate the disturbance kinetic energy associated with the perturbations as

where $\textbf{U}=\Gamma\,y\,\hat{\textbf{x}}$ is the base state simple shear flow, is subtracted from u - the total fluid velocity obtained from EL simulations, to obtain the disturbance fluid velocity. The integration is performed over the entire volumetric domain of the periodic simulation box.

Figure 12(a) shows the evolution of total perturbation energy ($E$), scaled with the initial value ($E(t=0)=E_{0}$), computed from EL simulations, presented in a semi-logarithmic scale for various non-dimensional wavenumbers $k$, corresponding to a mass loading of $M=1$. The plots exhibit a linear trend at initial times, indicating an exponential dependence on time. If the growth rate of velocity perturbations is denoted as $\sigma$, then the growth rate of energy perturbations would be $2\,\sigma$. Based on this argument, one may derive an expression to assess the growth rate $\sigma$ from EL simulations as

Figure 12(b) shows the time evolution of the growth rate estimated for EL simulation results using the expression given in equation (18). Unlike the prediction from LSA, which suggests a constant value, the instability growth rate exhibits variation over time in EL simulations. The figure indicates that the growth rate peaks between $t=0$ and $t=2$, depending on the seeded mode. The initial transient behaviour may attributed to imperfections in the initial conditions. The time window surrounding the peak growth rate in figure 12(b) corresponds to the exponential growth of kinetic energy observed in figure 12(a), thereby closely aligning with the dynamics predicted by LSA. Eventually, the growth rate decreases as the perturbation kinetic energy saturates, with the possibility of later increase due to nonlinear effects. An estimate for comparison with LSA would be the peak value of the growth rate for EL simulations obtained from equation (18). The non-dimensional growth rates $\sigma$ are thus obtained and plotted against the corresponding non-dimensional wavenumber $k$ for two different mass loading values in figure 13(a), represented by circle markers. The growth rate numerically obtained from LSA for a smooth number density profile ($m=4$) is also compared. The results show good agreement, although some discrepancies are observed in the numerical values. Considering that we are comparing completely different simulation classes: fully nonlinear EL simulations with LSA of the Euler-Euler model for the system under consideration, these differences are expected. Furthermore, the growth rate obtained from EL simulations indicates that energy perturbations decay beyond the critical wavenumber, contrary to what linear theory predicts. This could be attributed to the viscous effects in EL simulations.

In order to quantify the contribution of a specific wave mode $k$ to the overall dynamics, we compute the kinetic energy associated with this mode using Fourier decomposition as

where $\overline{(\cdot)}$ denotes complex conjugation, and $\hat{\textbf{u}}_{k}$ is the Fourier amplitude of the perturbation velocity field, given by

The normalization for the Fourier amplitude is chosen such that Plancherel’s identity yields $\int E_{k}\,dk=E$. The peak growth rate corresponding to the seeded eigenmode can thus be extracted using equations (19) with (18), and are also shown in figure 13 using diamond markers. These growth rates are also adequately consistent, qualitatively and quantitatively, with the LSA results. Additionally, the figure shows an average estimate of the growth rate over a non-dimensional time window $\Delta t=0.5$ centred around the peak value of total perturbation energy, depicted using square markers, which also compares well with the LSA results. The differences between these various growth rate estimates can also be observed to be very small.

Figure 13(b) displays the growth rate estimates from EL simulations for the $k=0.5$ case plotted for various mass loading values. As the mass loading increases, the growth rate also increases, indicating that the strength of the two-way coupling determines the strength of the instability. The trend in growth rate aligns well with the LSA results (shown as a dotted curve). Below a critical mass loading, however, the EL simulations show damped modes instead of the neutral modes observed in LSA, likely due to the inherent viscous nature of EL simulations, as mentioned earlier.

In this study, we have investigated a novel type of instability arising from the interaction between a simple shear flow and a non-uniform distribution of particles. Individually, each component is stable; a simple shear flow remains stable against infinitesimal perturbations, and a localised band of particles is stable without any gravitational or buoyancy effects. However, when considering the system as a whole, the momentum feedback from the particles to the fluid phase renders the system unstable. We utilised EL simulations to show the existence of instability, treating the fluid phase as a continuum and the particles as discrete using Eulerian and Lagrangian descriptions, respectively. The two-way coupling between the fluid and particle phases is crucial for the observed instability. The strength of the momentum feedback from particles to the fluid is directly proportional to the particle mass fraction and inversely proportional to particle inertia. Consequently, increasing the particle mass fraction enhances the instability while increasing particle inertia suppresses it. The simulations were conducted at high Reynolds numbers, indicating that the instability originates from inviscid effects. Moreover, we conducted simulations where we disabled the momentum feedback from particles to the fluid, while retaining the momentum exchange from the fluid to the particle phase. As anticipated, the system remained stable under these conditions, with particles solely advected by the shear flow.

Furthermore, a continuum model of the particle phase is employed to illustrate this instability’s modal nature and explain its generation mechanism. In the semi-dilute limit, the particle phase can be modelled as a separate pressureless fluid. The model adopts an Eulerian approach for particle and fluid phases, known as the two-fluid model. Linear stability analysis (LSA) of this model was conducted under the $St=0$ and $St\ll 1$ limits to quantify the growth rate of instability, utilising analytical and numerical techniques. The analytical investigation primarily focused on a simple top-hat distribution of particles, while a numerical treatment was reserved for smooth, localised particle distributions. In both scenarios, the analysis demonstrated that instability emerges as the momentum feedback from particles to the fluid phase is considered. Surprisingly, the instability persists even in the $St\rightarrow 0$ limit, suggesting that although it originates from particle inertia, its existence does not solely depend on inertia when the fluid-particle coupling is accounted for. However, increasing $St$ results in a damping effect on the growth rate of instability. The growth rate estimates from EL simulations matched the LSA results reasonably well.

To gain a mechanistic view of the instability, we examine the regime of small particle inertia ($St\rightarrow 0$). In this limit, the two-fluid model simplifies into a single-fluid model that describes the continuum evolution of a fluid with an effective density due to the suspended particles. The non-uniformity in particle distribution is thus reflected in the effective density of this fluid. The system now resembles a density-stratified fluid capable of supporting propagating edge waves at interfaces where vorticity or density exhibits discontinuities. Our system has two density interfaces, which individually support propagating edge waves generated due to non-Boussinesq baroclinic effects in isolation in a background shear flow. Unlike in the case of classical edge waves, the baroclinic torque arises from the misalignment between base state density stratification and disturbances in fluid acceleration. The propagating waves at both interfaces feel each other’s presence once brought closer, and they interact through phase locking and mutual amplification, culminating in the instability observed.

The present analysis has made several assumptions, and the relaxation of many of those would augment the instability found in the present study. The background flow is a simple shear flow, thus disallowing any shear instabilities that might arise (drazin2004hydrodynamic). We have also neglected the effects of gravitational settling, non-zero fluid inertia, concentration-dependent viscosity and hydrodynamic interactions. As was discussed earlier, in a simple shear flow of uniformly distributed negatively buoyant inertial particles, velocity perturbations can lead to preferential concentration of particles in regions of low vorticity. kasbaoui2015preferential showed that the variation of gravitational force exerted on the suspension thus induced would result in an algebraic instability. Ignoring gravitational effects, if the effects of fluid inertia were included instead, the particles would experience lift forces. drew1975lift showed that the lift force can lead to an instability in a bounded simple shear flow. A departure from the diluteness approximation would manifest first as a concentration-dependent viscosity and then as non-zero particle pressure. With the former physics, non-uniformity in particle concentration can lead to viscosity contrasts in the suspension, potentially causing short-wave viscous instabilities (see, e.g., hooper1983shear; hinch1984note; sahu2011linear). The latter, particle pressure due to hydrodynamic interactions, is responsible for shear-induced migration in dense suspensions (leighton1987shear). Shear-induced migration would be absent in a simple shear flow. However, the perturbation-induced migration could lead to a short-wave instability depending on the importance of concentration-driven diffusion relative to shear-induced migration (goddard1996migrational). Coupling even some of the physics mentioned above with the novel instability found in the present study would pave the way for a better understanding of instabilities in dusty shear flows and how they can have fundamentally different transition routes than their particle-free counterparts.

[Funding]A.R. acknowledge SERB project CRG/2023/008504 for funding. A.V.S.N. thanks the Prime Minister’s Research Fellows (PMRF) scheme, Ministry of Education, Government of India. A.R. and A.V.S.N. acknowledge the support of the Geophysical Flows Lab (GFL) at IIT Madras. M.H.K acknowledges support from the US National Science Foundation (award #2148710, CBET-PMP).

[Declaration of interests]The authors report no conflict of interest.

[Author ORCIDs]A. V. S. Nath, https://orcid.org/0000-0003-2144-2978; A. Roy, https://orcid.org/0000-0002-0049-2653; M. H. Kasbaoui, https://orcid.org/0000-0002-3102-0624.