2.1. Transient problem

To write the problem in non-dimensional form, we use the channel height $l_c=h$, the average velocity $u_c=Q/h$, the residence time $t_c=h^2/ Q$, and $p_c=\mu _2 Q/h^2$, as the characteristic length, velocity, time and pressure. Hereafter, all new variables refer to non-dimensional variables, and those introduced previously refer to their non-dimensional counterparts scaled with their characteristic values defined above.

Using this scaling, the dimensionless density and surface tension are the Reynolds number $Re=\rho Q/\mu _2$ and the inverse of the capillary number $Ca = \mu _2 Q / \gamma h$. With these non-dimensional variables, the continuity and momentum equations yield

(2.1)\begin{gather} \boldsymbol{\nabla}\boldsymbol{\cdot} \boldsymbol{v}_{i} =0 , \end{gather}

(2.2)\begin{gather}Re \left(\dfrac{\partial \boldsymbol{v}_{i}}{\partial t}+ \boldsymbol{v}_{i} \boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{v}_{i} \right) = \boldsymbol{\nabla}\boldsymbol{\cdot} \hat{\boldsymbol{T}}_i , \end{gather}

with the reduced stress tensor given by $\hat {\boldsymbol {T}}_i=-\hat {p}_i \boldsymbol {I} + \mu ( \boldsymbol {\nabla } \boldsymbol {v}_i + \boldsymbol {\nabla } \boldsymbol {v}_i^{\rm T} )$, and $\boldsymbol {v}_i$ and $\hat {p}_i$ the velocity and reduced pressure of fluid $i$, respectively. The viscosity ratio is $\mu =\mu _1/\mu _2$ if $y < \varGamma$, and $\mu =1$ if $y > \varGamma$.

To simplify the formulation and computation of the problem, (2.1)–(2.2) are written in terms of the reduced stress tensor $\hat {\boldsymbol {T}}_i$, with the reduced pressure field defined in terms of the pressure $p_i$ and the external force, $\hat {p}_i = p_i -\boldsymbol {f}\boldsymbol {\cdot } (\boldsymbol {x}-\boldsymbol {x}_p)$. The uniform volumetric force $\boldsymbol {f}$ that is acting on both liquids is assumed to have only non-zero vertical component $\boldsymbol {f}=f \boldsymbol {e_y}$. Introducing the stress tensor $\boldsymbol {T}_i = - p_i \boldsymbol {I} + \mu ( \boldsymbol {\nabla } \boldsymbol {v}_i + \boldsymbol {\nabla } \boldsymbol {v}_i^{\rm T} )$ in (2.2), and taking into account that $\widehat {\boldsymbol {T}_i} = \boldsymbol {T}_i + \boldsymbol {f} \boldsymbol {\cdot } ( \boldsymbol {x}-\boldsymbol {x}_p ) \boldsymbol {I}$, we would recover the momentum equation written in the traditional way.

To simplify the characterization of the position of the interface and the particle in order to carry out the computation and interpretation of the results, we define the normalized positions of the particle and the interface, $\xi = [ y_p-(\varGamma _{L/2}+d/2) ] / [ 1-(d+\varGamma _{L/2}) ]$ and ${\eta =\varGamma _{L/2}/(1-d)}$, with $\varGamma _{L/2}=\varGamma (L/2)$. The variable $\xi$ is in the range $0 \leq \xi \leq 1$, and the two extreme values of $\xi$ correspond to the particle touching the interface ${y_p=\varGamma _{L/2}+d/2}$ when $\xi =0$, or the channel wall $y_p=1-d/2$ when $\xi =1$, respectively (see figure 2). The normalized position of the interface $0\leq \eta \leq 1$ defines the extreme cases in which only fluid 2 runs through the channel ($\eta =0$) and in which the particle is squeezed between the interface and the upper wall ($\eta =1$), as illustrated in figure 2.

Figure 2. Definitions of the normalized positions (a,b) of the particle $\xi$ with arbitrary $\eta$, and (c,d) of the interface $\eta$.

2.2. Boundary conditions of the transient problem

The system of equations given by (2.1)–(2.2) is subject to the following boundary conditions. First, we enforce equilibrium of forces on the particle and assume that the particle inertia is negligible, $\int _{\varSigma _p} {\boldsymbol {T}}_2 \boldsymbol {\cdot } \boldsymbol {n}_p \,\mathrm {d} \varSigma = \boldsymbol {0}$. After substituting the stress tensor by its reduced counterpart and using Gauss's theorem to write $\int _{\varSigma _p} \boldsymbol {f} \boldsymbol {\cdot } (\boldsymbol {x}-\boldsymbol {x}_p) \boldsymbol {n}_p \,\mathrm {d} \varSigma = \int _{\mathcal {V}_p} \boldsymbol {\nabla }\{\boldsymbol {f} \boldsymbol {\cdot }(\boldsymbol {x}-\boldsymbol {x}_p)\} \,\mathrm {d}\mathcal {V} = \mathcal {V}_p \boldsymbol {f}$, this condition reads

(2.3)\begin{equation} \int_{\varSigma_p} \hat{\boldsymbol{T}}_2 \boldsymbol{\cdot} \boldsymbol{n}_p \,\mathrm{d} \varSigma = \mathcal{V}_p \boldsymbol{f} , \end{equation}

with $\boldsymbol {n}_{p}$ the unit vector normal to the particle pointing towards the fluid, and $\mathcal {V}_p$ the volume occupied by the particle. Thus the external force $\boldsymbol {f}$ acting on the liquids can also be understood as the volumetric migration force acting on the particle with hydrodynamics origin. Note that the vector equation (2.3) determines the terminal velocity $V$ and the body force $f$.

Furthermore, on the surface of the particle, we also impose the torque balance and the no-slip condition, considering that the particle moves as a rigid solid:

(2.4a)\begin{gather} 0 = \int_{\varSigma_{p}} \boldsymbol{e}_z \boldsymbol{\cdot} \left\{ \hat{\boldsymbol{T}}_2 \times (\boldsymbol{x}-\boldsymbol{x}_p)\right\} \boldsymbol{\cdot} \boldsymbol{n}_p \,\mathrm{d}\varSigma , \end{gather}

(2.4b)\begin{gather}\boldsymbol{v}_{2} = \boldsymbol{\varOmega} \times (\boldsymbol{x}-\boldsymbol{x}_p) + \frac{{\rm d}\kern0.7pt\boldsymbol{x}_p}{{\rm d}t} \quad \hbox{at}\ \varSigma_{p} , \end{gather}

where $\boldsymbol {e}_z=\boldsymbol {e}_x \times \boldsymbol {e}_y$, and the rotation occurs within the plane with rotational velocity $\boldsymbol {\varOmega }= \varOmega \boldsymbol {e}_z$.

No-slip boundary conditions are imposed at the channels walls $y=0$ and $y=1$, i.e.

(2.5)\begin{equation} \boldsymbol{v}={-}V \boldsymbol{e}_x , \end{equation}

whereas at the borders of the computational domain $x=\pm L/2$, we assume periodic boundary conditions

(2.6a–c)\begin{equation} \left.\boldsymbol{v}\right|_{L/2} = \left.\boldsymbol{v}\right|_{{-}L/2} , \quad \left.\partial_x \boldsymbol{v}\right|_{L/2} = \left.\partial_x \boldsymbol{v}\right|_{{-}L/2} ,\quad \left.p\right|_{L/2}= \left.p\right|_{{-}L/2} - \Delta p . \end{equation}

The problem under study is a pressure-driven flow, with the pressure drop $\Delta p$ being equal to the corresponding value for the Poiseuille flow plus a correction term due to the presence of the particle. This pressure drop is unknown initially, and it is calculated as part of the solution after the position of the particle is imposed. To calculate $\Delta p$ and the location of the interface position $\varGamma$, we enforce the total flow rate $1=Q_1+Q_2$ and the flow rate ratio $Q_1/Q_2$ such that

(2.7)\begin{gather} \int_{0}^{\varGamma} (u_1+V) \,\mathrm{d} y + \int_{\varGamma}^{1} (u_2+V) \,\mathrm{d} y = 1 , \end{gather}

(2.8)\begin{gather}\int_{0}^{\varGamma} (u_1+V) \,\mathrm{d} y = \frac{Q_1}{Q_2} \int_{\varGamma}^{1} (u_2+V) \,\mathrm{d} y , \end{gather}

where it has been taken into consideration that the average velocity has been chosen as the characteristic velocity to render the problem non-dimensional. Once $\mu _1/\mu _2$ and ${\eta }$ are chosen, there is only a value of $Q_1/Q_2$ that satisfies mass conservation $Q_{1}+Q_{2}=1$.

The implicit equation $q(x,y,t) \equiv \varGamma (x,t)-y=0$ is used to describe the evolution of the interface. Because $q=0$ defines the interface at all times, the material derivative must satisfy

(2.9)\begin{equation} \dfrac{\partial q}{\partial t} + \boldsymbol{v} \boldsymbol{\cdot} \boldsymbol{\nabla} q = 0 \quad \hbox{at} \ \varGamma . \end{equation}

Then the unit vector normal to the surface $\varGamma$ pointing from fluid 2 towards fluid 1 is defined using the function $q$ as $\boldsymbol {n}=\boldsymbol {\nabla } q/|\boldsymbol {\nabla } q|=(\varGamma _x \boldsymbol {e}_x-\boldsymbol {e}_y)/\sqrt {1+\varGamma _x^2}$.

Additionally, at the fluid–fluid interface, we also impose the continuity of velocities and the jump condition on the stress tensor,

(2.10a)\begin{gather} {}[\boldsymbol{v}] = \boldsymbol{0}, \end{gather}

(2.10b)\begin{gather}\boldsymbol{n} \boldsymbol{\cdot} [ \hat{\boldsymbol{T}}] = Ca^{{-}1} \left( -\boldsymbol{n}\,\boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{n} \right) \quad \text{at} \ \varGamma, \end{gather}

with the bracket operator indicating the jump of the variable included between them, e.g. ${[{\mathcal {L}}]={\mathcal {L}}_2-{\mathcal {L}}_1}$. Note that $[\hat {\boldsymbol {T}}]$ has been used, since $[\boldsymbol {T}]= [\hat {\boldsymbol {T}}]$.

2.3. Global linear stability problem

To initiate the analysis, we calculate the steady-state solution of the problem's variables ${\hat {p}_{i,0}, \boldsymbol {v}_{i,0}, \varGamma _0, V_0, \varOmega _0, f, \Delta p}$ for a given position of the particle $\boldsymbol {x}_{p,0}$ by integrating the system of equations described in the previous subsection without the unsteady terms in (2.2), (2.4b) and (2.9). This steady-state solution, denoted by the subscript $0$, was studied by Ruiz-Martín et al. (Reference Ruiz-Martín, Rivero-Rodriguez and Sánchez-Sanz2022b). To perform a global stability analysis, we superimpose a time-periodic small-amplitude disturbance upon this steady two-dimensional solution. The flow will be said to be globally stable if it is stable to all finite-amplitude perturbations (Theofilis Reference Theofilis2003).

To carry out the stability analysis, we perturb both the interface and the position of the particle. To do so, we start by introducing a small normal disturbance $\delta \boldsymbol {n}_{0}$ to the unperturbed fluid–fluid interface $\varGamma _0$, such that $\boldsymbol {x} + \delta \boldsymbol {n}_{0} \in \varGamma$ for $\boldsymbol {x} \in \varGamma _0$, where $\varGamma$ represents the perturbed interface. Similarly, the position of the centre of mass of the particle is perturbed as ${\boldsymbol {x}_{p} = \boldsymbol {x}_{p,0} + \Delta \boldsymbol {x}_{p}}$, such that $\boldsymbol {x}_p$ is the centre of mass of $\varSigma _{p}$, and $\boldsymbol {x}_{p,0}$ is the centre of mass of $\varSigma _{p,0}$, where $\Delta \boldsymbol {x}_p$ is a small vertical displacement of the particle such that $\Delta \boldsymbol {x}_p \boldsymbol {\cdot } \boldsymbol {e}_x =0$. Alternatively, it is sometimes more convenient to consider the normal displacement of the surface $\delta _p$, such that $\boldsymbol {x} + \delta _p \boldsymbol {n}_0 \in \varSigma _p$ for $\boldsymbol {x} \in \varSigma _{p,0}$, which also lies on the perturbed surface as long as $\delta _p = \Delta \boldsymbol {x}_p \boldsymbol {\cdot } \boldsymbol {n}_0$, where $\boldsymbol {n}_0=-\boldsymbol {n}_p$ is the normal of the unperturbed surface.

Notice that $\delta$ and $\Delta \boldsymbol {x}_{p}$ are a priori unknown and will be determined as part of the calculation process. For this reason, to impose the boundary conditions at the perturbed geometry, our resolution strategy starts by writing (2.3), (2.4a) and (2.4b) for the particle, and (2.9) and (2.10) for the interface on the unperturbed geometry (see Appendix B for a detailed explanation of this step). After doing this, in a second step, we perform the stability analysis for the variables by introducing the expansion

(2.11a)\begin{gather} \psi(\boldsymbol{x}, t) = \psi_0(\boldsymbol{x}) +{\epsilon\,{\rm e}^{\lambda t}}\,\psi_1 (\boldsymbol{x}) , \end{gather}

(2.11b)\begin{gather}\delta(\boldsymbol{x}, t) = {\epsilon\,{\rm e}^{\lambda t}}\,\delta_1 (\boldsymbol{x}) , \end{gather}

(2.11c)\begin{gather}\delta_p(\boldsymbol{x}, t) = {\epsilon\,{\rm e}^{\lambda t}}\,\delta_{p,1} (\boldsymbol{x}) , \end{gather}

(2.11d)\begin{gather}\boldsymbol{x}_p(t) = \boldsymbol{x}_{p,0} + {\epsilon\,{\rm e}^{\lambda t}}\,\boldsymbol{x}_{p,1} , \end{gather}

where $\psi =\{\hat {p}_i, \boldsymbol {v}_i, \varGamma, V, \varOmega$}, and the subscript $0$ indicates the steady-state values of the corresponding variables. The small vertical displacement of the particle is ${\Delta \boldsymbol {x}_p={\epsilon \,{\rm e}^{\lambda t}}\,\boldsymbol {x}_{p,1}}$ with $\boldsymbol {x}_{p,1} \boldsymbol {\cdot } \boldsymbol {e}_x=0$, and the amplitude of the perturbation is small, ${{\epsilon \,{\rm e}^{\lambda t}} \ll 1}$. Additionally, $\delta _{p,1} = \boldsymbol {x}_{p,1} \boldsymbol {\cdot } \boldsymbol {n}_0$ is the first-order displacement of the particle in the direction normal to $\varSigma _{p,0}$. Notice that neither $\Delta p$ nor $\boldsymbol {f}$ is expanded.

The eigenvalue $\lambda =\lambda _r + {\rm i} \lambda _i$ is a complex number, with its real $\lambda _r$ and imaginary $\lambda _i$ parts representing the perturbation's growth rate and its oscillation frequency, respectively. The base flow is then unstable when the real part of at least one eigenvalue is positive, $\lambda _{r}>0$. Unlike Salgado Sánchez (Reference Salgado Sánchez2020), we do not consider perturbations in the value of the volumetric force $\boldsymbol {f}$, which is why this variable is not in the list of expanded variables in (2.11). Note that for a generic variable $\psi _{i,j}$, the subscript $j=0$ represents the base state, while $j=1$ represents the linearized variable. The subscript $i$ has been defined previously to refer to the lower layer when $i=1$ and the upper layer when $i=2$. Thus $\boldsymbol {v}_{1,0}$, $\boldsymbol {v}_{1,1}$ are the velocity fields of the lower liquid for the unperturbed and linearized problems, and $\boldsymbol {v}_{2,0}$, $\boldsymbol {v}_{2,1}$ represent the base-state and perturbed velocities corresponding to the upper liquid. Subsequently, the jump of the variable across the interface is denoted as ${[{\mathcal {L}}]_j={\mathcal {L}}_{2,j}-{\mathcal {L}}_{1,j}}$.

Substituting (2.11a) in the dimensionless Navier–Stokes (NS) equations, we obtain for order $O(\epsilon )$ the following system for the linearized problem:

(2.12)\begin{gather} \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{v}_{i,1} = 0 , \end{gather}

(2.13)\begin{gather}Re\, (\alpha_{NS} \lambda \boldsymbol{v}_{i,1} + \boldsymbol{v}_{i,0} \boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{v}_{i,1} + \boldsymbol{v}_{i,1} \boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{v}_{i,0}) =\boldsymbol{\nabla} \boldsymbol{\cdot} \hat{\boldsymbol{T}}_{i,1} , \end{gather}

with $\hat {\boldsymbol {T}}_{i,1}=-\hat {p}_{i,1}\boldsymbol {I} + \mu ( \boldsymbol {\nabla } \boldsymbol {v}_{i,1} + \boldsymbol {\nabla } \boldsymbol {v}_{i,1}^{\rm T} )$ the perturbed reduced stress tensor at first order. The artificial parameter $\alpha _{NS}$ is introduced in the equation to study the contribution of the shear mode on the interfacial stability (see Appendix C). In all calculations shown below, $\alpha _{NS}=1$, with $\alpha _{NS}=0$ only in Appendix C to identify the influence of this term in the development of the instability.

2.4. Boundary conditions of the stability problem

This subsection describes the boundary conditions to be imposed in the perturbed problem. To facilitate the reading of this subsection, the algebra has been reduced to a minimum. The interested reader can find all mathematical details in the Appendix B.

As explained above, we perturb the position of the particle by introducing a small displacement $\Delta \boldsymbol {x}_p$. The perturbed position of the particle $\boldsymbol {x}_p$ is unknown and will be obtained as part of the calculation process. Consequently, the boundary conditions (2.3), (2.4a) and (2.4b) cannot be imposed explicitly. To circumvent this difficulty, we followed the methodology developed by Rivero-Rodriguez, Perez-Saborid & Scheid (Reference Rivero-Rodriguez, Perez-Saborid and Scheid2018) to write these boundary conditions referred to the unperturbed particle's position $\boldsymbol {x}_{p,0}$. This process resulted in (B3), (B5) and (B6), respectively, as outlined in detail in Appendix B. Once this is done, we introduce the expansion (2.11) in (B3), (B5) and (B6) to give the boundary conditions at the surface of the particle to first order in $\epsilon$:

(2.14a)\begin{gather} \int_{\varSigma_{p,0}} \,\left(\boldsymbol{n}_{0} \boldsymbol{\cdot} \hat{\boldsymbol{T}}_{2,1} + Re \,\delta_{p,1} \boldsymbol{v}_{2,0} \boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{v}_{2,0}\right)\, {\rm d}\varSigma =\boldsymbol{0} , \end{gather}

(2.14b)\begin{gather}\int_{\varSigma_{p,0}} \boldsymbol{e}_z \boldsymbol{\cdot} \,\left\{\left(\boldsymbol{n}_{0} \boldsymbol{\cdot} \hat{\boldsymbol{T}}_{2,1} + \delta_{p,1}\,Re \,\boldsymbol{v}_{2,0} \boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{v}_{2,0}\right) \times (\boldsymbol{x}-\boldsymbol{x}_{p,0} ) - \boldsymbol{n}_{0} \boldsymbol{\cdot} \hat{\boldsymbol{T}}_{2,0} \times \boldsymbol{x}_{p,1}\,\right\} {\rm d}\varSigma =0 , \end{gather}

(2.14c)\begin{gather}\boldsymbol{v}_{2,1} + \delta_{p,1} \boldsymbol{n}_0 \boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{v}_{2,0} = \boldsymbol{\varOmega}_1 \times (\boldsymbol{x} -\boldsymbol{x}_{p,0})+ \boldsymbol{\varOmega}_0 \times (\delta_{p,1} \boldsymbol{n}_0 - \boldsymbol{x}_{p,1}) +\alpha_\epsilon \lambda \boldsymbol{x}_{p,1} \quad \text{at} \ \varSigma_{p,0} . \end{gather}

The artificial parameter $\alpha _{\epsilon }$ takes the value $\alpha _{\epsilon }=1$ in all cases except in Appendix C, when we evaluate the effect of the migration mode in the solution of the linearized system and we impose $\alpha _{\epsilon }=0$.

Next, the velocities at the upper and lower walls, $y=1$ and $y=0$, are also expanded:

(2.15)\begin{equation} \boldsymbol{v}_{i,1} ={-} V_1 \boldsymbol{e}_x , \end{equation}

where $i=1$ for the lower wall, and $i=2$ for the upper one. The perturbation in the particle's position in the horizontal direction is considered through a perturbation of the terminal velocity that keeps the inter-particle distance $L$ constant to conserve the periodic character of the problem. Periodic boundary conditions are expanded as well:

(2.16a–c)\begin{equation} \left.\boldsymbol{v}_{i,1}\right|_{L/2} = \left.\boldsymbol{v}_{i,1}\right|_{{-}L/2} , \quad \left.\partial_x \boldsymbol{v}_{i,1}\right|_{L/2} = \left.\partial_x \boldsymbol{v}_{i,1}\right|_{{-}L/2} ,\quad \left.p_{i,1}\right|_{L/2}= \left.p_{i,1} \right|_{{-}L/2} . \end{equation}

Additionally, it is worth mentioning that the pressure drop $\Delta p$ has not been expanded.

Following the same procedure as with the particle, the interface separating the fluids is perturbed by applying a normal differential displacement $\delta \boldsymbol {n}_0$ that creates a volume as a result of the displacement of the interface from $\varGamma _0$ to its perturbed position $\varGamma (\boldsymbol {x},t)$, where the boundary conditions of the problem must be enforced. To do so, we write the boundary conditions referred to their unperturbed counterpart $\varGamma _0$, and we substitute the perturbation (2.11) in the resulting equations (B7), (B10) and (B11). Collecting terms of first order, the kinematic condition (2.9), the continuity of velocities (2.10a) and the stress balance (2.10b) can be written as

(2.17a)\begin{gather} \boldsymbol{v}_{2,1} \boldsymbol{\cdot} \boldsymbol{n}_0 - {\boldsymbol{D}}_{S} \boldsymbol{\cdot} ( \delta_{1} \boldsymbol{v}_{2,0} ) = \alpha_\delta \lambda \delta_1 \quad \mbox{at}\ \varGamma_0 , \end{gather}

(2.17b)\begin{gather}{}[\boldsymbol{v}]_1 + \delta_1 \boldsymbol{n}_0 \boldsymbol{\cdot} \boldsymbol{\nabla} [\boldsymbol{v}]_0= \boldsymbol{0} \quad \mbox{at}\ \varGamma_0 , \end{gather}

(2.17c)\begin{gather}{}[\hat{\boldsymbol{T}}]_1 \boldsymbol{\cdot} \boldsymbol{n}_0 - {\boldsymbol{D}}_{S} \boldsymbol{\cdot} \left(\delta_1 [\hat{\boldsymbol{T}}]_0\right) + Re \,\delta_1 [\boldsymbol{v} \boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{v}]_0 = Ca^{{-}1}\,{\boldsymbol{D}}_{S} \boldsymbol{\cdot} \boldsymbol{\mathcal{B}}_1 \quad\mbox{at}\ \varGamma_0 , \end{gather}

where ${\boldsymbol {D}}_{S} \boldsymbol {\cdot } \boldsymbol {b} = \boldsymbol {\nabla }_S \boldsymbol {\cdot } \boldsymbol {b} - (\boldsymbol {\nabla }_S \boldsymbol {\cdot } \boldsymbol {n} )\boldsymbol {n} \boldsymbol {\cdot } \boldsymbol {b}$ and $\boldsymbol {\nabla }_S \boldsymbol {b}=(\mathcal {I}-\boldsymbol {n}\boldsymbol {n})\boldsymbol {\cdot }\boldsymbol {\nabla } \boldsymbol {b}$ for any arbitrary vectorial or tensorial field $\boldsymbol {b}$. The tensor ${\boldsymbol {\mathcal {B}}_1}$ is the first term of the perturbation of the rotation tensor $\boldsymbol {\mathcal {B}}= \boldsymbol {n}_{S_0} \boldsymbol {n}_{S}$, with $\boldsymbol {n}_{S_0}$ and $\boldsymbol {n}_{S}$ the tangential vectors to the unperturbed $\varGamma _0$ and perturbed interface $\varGamma$. As particularized from the three-dimensional case given in Rivero-Rodriguez et al. (Reference Rivero-Rodriguez, Perez-Saborid and Scheid2018), the tangential vector reads $\boldsymbol {n}_{S}=\boldsymbol {n}_{S_0} + \boldsymbol {n}_{S_0} \boldsymbol {\cdot } (\boldsymbol {\nabla }_S \delta ) \boldsymbol {n}_{0}$, to give, after using the expansion (2.11b), the term of order $O(\epsilon )$ of the rotation tensor ${\boldsymbol {\mathcal {B}}_1 = (\boldsymbol {\nabla }_S\delta _1) \boldsymbol {n}_{0}}$ (Ruiz-Martín et al. Reference Ruiz-Martín, Rivero-Rodriguez and Sánchez-Sanz2022b). The contribution of the interfacial mode on the stability of the flow is evaluated by introducing the artificial parameter $\alpha _{\delta }$, which takes the value $\alpha _{\delta }=0$ only when we study the importance of this term in the stability of the interface (Appendix C).

2.5. Numerical method

Once the problem is formulated in non-dimensional form, and taking into account that the size of the particle is set constant and equal to $d=0.2$, it is possible to identify the parametric dependence of the variables of the base-state problem $\psi _0$, $f$, $\Delta p$ with $\mu _1/\mu _2$, $Re$, $Ca$, $\eta$, $\xi$ and $L$.

The functional dependence will be obtained by integrating numerically the system of equations detailed above in (2.1)–(2.2). Table 1 provides the range of variation of the values of the non-dimensional parameters. After imposing the parameters of the problem, the flow variables $\psi$ are calculated simultaneously using a Newton's iterative method that continues until the error is below $10^{-6}$ (see Ruiz-Martín et al. Reference Ruiz-Martín, Rivero-Rodriguez and Sánchez-Sanz2022b).

Table 1. Relevant non-dimensional parameters.

Once the base flow has been computed, we can solve the linearized system of equations formed by (2.12) and (2.13) with the boundary conditions (2.14), (2.15), (2.16a–c) and (2.17), such that one obtains the eigenfunctions of the velocity $\boldsymbol {v}_{i,1}$ and pressure fields $\hat {p}_{i,1}$, the correction for the terminal $V_1$ and rotational $\varOmega _1$ velocities, particle displacement $\boldsymbol {x}_{p,1}$ and interface displacements $\delta _1$ along with the eigenvalues $\lambda$.

The system of equations is solved using the finite element method implemented in Comsol Multiphysics. Linear Lagrangian elements have been used for pressure, and quadratic Lagrangian elements for the rest of the variables. The liquid–liquid interface at zeroth order has been solved with the help of the arbitrary Lagrangian Eulerian method also implemented in the software.

4.1. Particle position effect

To examine the influence of the particle position on the stability of the interface, we show in figure 4 the neutral stability diagrams in the $L\unicode{x2013}\xi$ parametric space for different viscosity and flow rate ratios. In this figure, vertical solid lines represent the value of $L$ beyond which the interface is unstable in the case with no particle. In particular, figure 4(b) shows that when the particle is immersed in the most viscous liquid ($\mu _{1}/\mu _{2}=0.10$), the presence of the particle has a clear destabilizing effect, narrowing the combination of parameters for which the interface is stable.

Figure 4. Neutral stability diagram in the $\xi \unicode{x2013}L$ plane (dashed lines) and comparison with the particle-free value of $L$ beyond which the interface is unstable (vertical solid lines) for $Re = 0.05$, $Ca^{-1}=1$ for different $\eta$ values, or equivalently flow rate ratios (see table 2), with (a) $\mu _{1}/\mu _{2}=10$, (b) ${\mu _1/\mu _2=0.10}$. In (a), the dashed orange line corresponding to $\eta =0.10$, $Q_{1}/Q_{2} = 2.40 \times 10^{-3}$ is not depicted because it is always unstable. In (b), the solutions for $\eta =0.10$, $Q_{1}/Q_{2} = 5.50 \times 10^{-2}$ (solid orange line) and $\eta =0.25$, $Q_{1}/Q_{2} = 2.25 \times 10^{-1}$ (solid brown line) are always stable in the particle-free case.

Interestingly, this figure illustrates how the particle can stabilize the interface when $\mu _1/\mu _2=10$. As shown in Yiantsios & Higgins (Reference Yiantsios and Higgins1988), the two-phase Poiseuille flow becomes unstable to long-wavelength perturbations when the thickness of the most viscous fluid is sufficiently low (thin layer effect). Remarkably, when the particle travels in the less viscous fluid ($\mu _{1}/\mu _{2}=10$), an increment of the volumetric flow rate shifted the unstable modes towards the cases in which the particle is close to either the interface $\xi < 0.2$ or the upper wall $\xi >0.6$ (see figure 4). In this case, intermediate values $0.2<\xi <0.6$ seem to be stable, $\lambda _r<0$, with the presence of the particle stabilizing the interface for $L>26$. For instance, the interface is stable for $Q_1/Q_2=1.68$ and $L>26$ (blue line in figure 4a) when the particle is located at $0.30<\xi <0.60$, with this region being unstable in the particle-free problem. It is worth mentioning that the interaction between the particles induced by interface perturbations with a region of influence higher than the inter-particle distance is not considered in this work, and the disturbance always fits within the inter-particle distance $L$. Hence as we increase $L$, new unstable modes that could not develop at smaller values might emerge.

Table 2. Values of the flow rate ratio $Q_{1}/Q_{2}$ for $\mu _{1}/\mu _{2}=10$ and $\mu _{1}/\mu _{2}=0.10$ for different values of $\eta$, corresponding to the lines plotted in figure 4.

4.2. Inertial and capillary effects

In this subsection, we investigate the influence of inertial and capillary effects in the two-dimensional problem. In the frame of the local stability analysis carried out by Yiantsios & Higgins (Reference Yiantsios and Higgins1988) at small or moderate Reynolds numbers, the imaginary wave speed $c_i$ of the two-layer Poiseuille flow depends linearly on the Reynolds number. Using our global stability analysis, we illustrate in figure 5 the dependence on the Reynolds number of the most dangerous eigenvalue $\lambda _r|_{max}$ for the interfacial mode. (Note that the most dangerous interfacial eigenvalue is that with the largest real part $\lambda _r|_{max}$ and $\lambda _i \ne 0$.) For the flow parameters considered here, $\lambda _r|_{max}$ is independent of the particle position $\xi$, and essentially follows the trend defined by the interface in the absence of the particle. For sufficiently low values of the Reynolds number, capillary effects become dominant and the particle changes the stability of the flow. Indeed, for $Re=0.05$ and the flow parameters considered in figure 5, the flow is stable only when the particle is sufficiently far from the interface $\xi >0.50$.

Figure 5. Most dangerous eigenvalue as a function of the Reynolds number for $Ca^{-1}=1$, with dotted black line representing the level $\lambda _{r}|_{max}=0$. The inset in (b) shows the solution in the zoomed region $Re<1$. The red line represents the particle-free configuration solution, the violet line is for the particle located at $\xi =0.90$, the green line is for $\xi =0.50$, and the blue line is for $\xi =0.10$. Flow parameters are $\mu _{1}/\mu _{2}=0.1$, $L=20$ and (a) $\eta =0.07$, $Q_{1}/Q_{2}=3.31 \times 10^{-2}$, (b) $\eta =0.40$, $Q_{1}/Q_{2}=0.59$.

In figure 6, we show the evolution of the absolute real part of the most dangerous eigenvalue with the inverse of the capillary number for these flow parameters. Figure 6 illustrates that the presence of the particle becomes more noticeable in the limit of small surface tension, whilst as $Ca^{-1} \gg 1$, the evolution of $\lambda _r|_{max}$ is independent of the particle position, and identical to the result obtained with the particle-free configuration. This is consistent with the fact that at high values of the surface tension, the interface deformation due to the presence of the particle becomes negligible, and consequently the influence of the particle on the interfacial mode decreases. The size of the perturbations is limited by the inter-particle distance, so that beyond a certain value of the surface tension, the eigenvalue is seen to increase linearly with $Ca^{-1}$. In figure 6(a), the greatest eigenvalue $\lambda _r|_{max}$ computed for small values of the volumetric flow ratio is negative for all values of $Ca^{-1}$ when the particle is far from the interface. For $\xi <0.6$, approximately, $\lambda _{r}|_{max}$ becomes positive at low surface tension, with small stability regions that emerged for intermediate values of $Ca^{-1}$ before $\lambda _{r}|_{max}$ finally achieves a positive constant value that is independent of the surface tension. As the flow ratio $Q_1/Q_2$ increases, the unstable behaviour of the interface emerges approximately at the same capillary number for all values of $\xi$, as illustrated in figure 6(b). This figure shows the variations on the evolution of the growth rate in the limit of small surface tension $Ca^{-1}$ due only to a change in the frequency of the wave.

Figure 6. Absolute real part of the most dangerous eigenvalue as a function of the inverse of the capillary number for $Re=0.05$ and the same flow parameters as considered in figure 5. Solid lines for stable (negative) eigenvalues and dashed lines for unstable (positive) values.

To study this in more detail, we represent in figure 7 the evolution of the most dangerous eigenvalue $\lambda _r|_{max}$ for the interfacial mode using $L=22$ and $Q_1/Q_2=2.05$. In this figure, we consider both the global stability solution as well as the results obtained for perturbations of wavenumber $k_L=2{\rm \pi} /L$. Note that the swift variations close to $Ca^{-1} = 10^{-3}$ and $Ca^{-1} =10^{-1}$ for the particle-free configuration are related to a change in the wavenumber of the perturbation with the maximum growth rate. In essence, the wavenumber of the most dangerous perturbation increases as $Ca^{-1} \to 0$, as the stabilizing effect of surface tension decreases at small values of the wavenumber. On the other hand, in the limit $Ca^{-1} \to \infty$, the wavelength is limited by the inter-particle length $L$, and beyond a certain value, the most dangerous eigenvalue is seen to depend linearly on $Ca^{-1}$.

Figure 7. Most dangerous eigenvalue as a function of the inverse of the capillary number. (a) The result obtained in the particle-free configuration is compared with the solution for perturbations of wavenumber $k = 2 {\rm \pi}n / L = n k_{L}$. (b) The evolution of the most dangerous eigenvalue in the two-dimensional problem for different particle positions is compared with the particle-free global stability analysis (red line). Solid lines for stable (negative) eigenvalues and dashed lines for unstable (positive) values. Flow parameters are ${Re=0.05}$, ${\mu _1/\mu _2=0.1}$, $\eta =0.60$, $Q_{1}/Q_{2}=2.04$, $L=22$.

Neutral stability diagrams for $\mu _1/\mu _2=0.10$ and different particle positions ($\xi =0.10$, $0.50$, ${0.90}$) are shown in figure 8 for $Re=0.05$ (figure 8(a) and solid lines in figure 8(b)) and for $Re=5$ (dashed lines in figure 8b). Figure 8 illustrates that at lower values of the Reynolds number, capillary migration becomes dominant and the interface deformation induced by the particle more significant. Consequently, the particle position has a higher impact on the stability of the flow, whereas for dominant inertial effects, the neutral lines collapse in the $\eta$–$L$ plane for the different particle positions. In addition, the unstable region for higher values of the Reynolds number widens towards lower values of the flow rate ratio and inter-particle distance (lower values of the perturbation wavelength). On the other hand, comparing the neutral stability diagram presented in figure 8(a) and that obtained for $Re=0.05$ and $Ca^{-1}=1$ with the results plotted in figure 8(b) with solid lines – which have also been obtained for $Re=0.05$ but for a higher surface tension $Ca^{-1}=10$ – we can see the stabilizing effect of the surface tension for lower values of the inter-particle distance. This is in agreement with the results obtained for the stability of the two-phase Poiseuille flow (Yiantsios & Higgins Reference Yiantsios and Higgins1988).

Figure 8. Neutral stability diagram in the $Q_{1}/Q_{2}\unicode{x2013}L$ plane for $\mu _{1}/\mu _{2}=0.10$ and (a) $Ca^{-1}=1$, $Re=0.05$, and (b) $Ca^{-1}=10$, $Re=0.05$ for solid lines and $Re=5$ for dashed lines. The particle position is $\xi =0.10$ for the blue line, $\xi =0.50$ for the green line, and $\xi =0.90$ for the violet line. The red line corresponds to the particle-free configuration.

If the particle is immersed in the most viscous fluid ($\mu _{1}/\mu _{2}\leq 1$), then the liquid–liquid interface becomes unstable at low values of the flow rate ratio when the particle is sufficiently close to the interface, independently of the inter-particle distance. Also, particle positions closer to the interface widen the unstable region towards lower values of the particle distance $L$. As mentioned before, the increase of the surface tension stabilizes the flow for relatively low values of $L$. Remarkably, the presence of the particle is stabilizing when the particle is in the less viscous fluid ($\mu _{1}/\mu _{2}=10$) in a narrow region at relatively high flow rate ratios ($Q_{1}/Q_{2} \sim 1$), regardless of the inter-particle distance considered (see figure 9). This stabilizing effect is more pronounced when the particle is far from the interface and appears only at low values of $Ca$ or large values of $Re$, as shown in figure 10.

Figure 9. Neutral stability diagram in the $Q_{1}/Q_{2}\unicode{x2013}L$ plane for different values of the particle position, $Re = 0.05$, $\mu _{1}/\mu _{2}=10$ and $Ca^{-1}=1$. The particle position is (a) $\xi =0.50$ and (b) $\xi =0.90$. Red lines correspond to the particle-free configuration. In (b), the shaded region represents the parameter range for which we were unable to compute a solution. The particular case highlighted by the star marker ($\bigstar$) in the figures is analysed in detail in figure 14.

Figure 10. Neutral stability diagram in the $Q_{1}/Q_{2}$–$L$ plane for $\mu _{1}/\mu _{2}=10$, $Ca^{-1}=10$ and (a,c,e) $Re=0.05$, (b,d,f) $Re=5$. Plots in each row are dedicated to a different particle position $\xi$. Red lines correspond to the particle-free configuration. Letters S and U denote the stable and unstable regions, respectively.

To understand further the effect of the particle on the interfacial instability, the eigenfunctions for the most unstable interfacial mode are plotted in figures 11 and 12, considering two sets of parameters. With the first set, the particle stabilizes the interface, whilst with the second set, the particle destabilizes the interface, which is stable in the particle-free problem. Due to the slenderness of the computational domain, the vertical axis has been stretched in figures 11 and 12 to facilitate the reading of the figures. The ellipsoidal geometry of the particle is only a visual consequence of that change of scale.

Figure 11. Contour plots of the perturbed streamlines (light grey lines). The positions of the unperturbed and perturbed interfaces are depicted with solid and dashed lines, respectively, for (a) a particle located at $\xi =0.90$ with $\lambda = -3.4486 \times 10^{-5}+0.19373{\rm i}$, and (b) the particle-free configuration with $\lambda = 4.4319 \times 10^{-6} +0.38565{\rm i}$. The rest of the flow parameters are $Re=0.05$, $\mu _{1}/\mu _{2}=10$, $\eta =0.69$ ($Q_{1}/Q_{2}=0.629$), $L=30$, $Ca^{-1}=1$.

Figure 12. Contour plots of the perturbed streamlines (light grey lines). The position of the unperturbed and perturbed interfaces are depicted with solid and dashed lines, respectively, for (a) a particle located at $\xi =0.90$ with $\lambda = 5.0833 \times 10^{-5}+0.20982{\rm i}$, and (b) the particle-free configuration with $\lambda = -1.6913 \times 10^{-6}+0.3139{\rm i}$. The rest of the flow parameters are $Re=0.05$, $\mu _{1}/\mu _{2}=10$, $\eta =0.95$ ($Q_{1}/Q_{2} = 3.17$), $L=30$, $Ca^{-1}=1$.

From figure 9, we identified a set of parameters in which the particle stabilizes the interface – $Re=0.05$, $\mu _{1}/\mu _{2}=10$, $\eta =0.69$ ($Q_{1}/Q_{2} = 0.629$), $\xi =0.90$, $L=30$, $Ca^{-1}=1$ – before plotting the real part of the eigenfunctions in figure 11. To facilitate the comparison, we included in figure 11(b) the eigenfunction in the case without particles.

On the other hand, figure 12 includes the eigenfunctions of the most unstable mode for the same flow parameters considered previously, but with the interface closer to the upper wall: $Re=0.05$, $\mu _{1}/\mu _{2}=10$, $\eta =0.95$ ($Q_{1}/Q_{2}=3.17$), $\xi =0.90$, $L=30$, $Ca^{-1}=1$. In this case, the particle induces the formation of a recirculating region near the particle that extends to affect the lower fluid, destabilizing the interface. The eigenfunction of the particle-free configuration is included in figure 12(b), illustrating the effect of the particle on the flow pattern. A more detailed plot of the region surrounding the particle is included in figure 13, comparing the cases $\eta =0.65$ and $0.95$. In this figure, we can identify the recirculating regions formed near the particle surface and at the lower and more viscous liquid.

Figure 13. Contour plots of the perturbed streamlines in the fluid near the particle for (a,c) $\eta =0.65$ ($Q_{1}/Q_{2} = 0.629$) with $\lambda = -3.4486 \times 10^{-5}+0.19373{\rm i}$, and (b,d) $\eta =0.95$ ($Q_{1}/Q_{2} = 3.17$) with $\lambda = 5.0833 \times 10^{-5}+0.20982{\rm i}$. The real part of the velocity field is at the top, and the imaginary part is at the bottom. The unperturbed interface position $\varGamma _0$ is shown as a solid red line.The rest of the flow parameters are $Re=0.05$, $\mu _{1}/\mu _{2}=10$, $\xi =0.90$, $L=30$, $Ca^{-1}=1$.

The analysis of the eigenfunctions suggests that these regions of flow recirculation in the lower liquid play a key role in destabilizing the interface. In multi-layer flows without particles, it has been suggested that the interfacial instability can be explained by the out-of-phase vorticity disturbances that appear as a result of the vorticity advection by the base flow (see Hinch Reference Hinch1984; Charru & Hinch Reference Charru and Hinch2000; Govindarajan & Sahu Reference Govindarajan and Sahu2014). For the destabilizing case discussed above ($Re=0.05$, $\mu _{1}/\mu _{2}=10$, $\eta =0.95$, $Q_{1}/Q_{2} = 3.17$, $\xi =0.90$, $L=30$, $Ca^{-1}=1$), the particle-free configuration is stable because the surface tension stabilizes short-wavelength interfacial disturbances. Also, long-wavelength disturbances remain stable (Charru & Hinch Reference Charru and Hinch2000) as the downstream convection of vorticity remains relevant only in the thicker, more viscous lower liquid layer. Interestingly, figure 12 suggests that the assumption of negligible advected vorticity is no longer valid in the thinner, less viscous liquid layer in the presence of particles. This figure shows a particle-induced vertical flow motion that extends to the lower and more viscous liquid. This effect is not observed in figure 11, when the particle stabilized the interface. An in-depth analysis of the physical mechanisms that trigger the instability is of unquestionable interest. However, that study is out of the scope of this work.

As the surface tension of the interface is reduced (lower values of $Ca^{-1}$), the maximum deformation of the interface becomes larger, achieving deformation amplitudes ${g = \varGamma (x)-\varGamma _{L/2}}$ of the order of the size of the particle. An example of that effect is shown in figure 14 for $Re=0.05$, $\mu _{1}/\mu _{2}=10$, $Q_{1}/Q_{2} = 0.35$ ($\eta = 0.60$), $\xi =0.90$, $L = 14$ and $Ca^{-1} = 1$. Keeping the other parameters unchanged, the range $L\in [5,40]$ gives a very small radius of curvature of the interface, which is what hinders the numerical calculation of the problem. This region is identified in figure 9(b) by using shading. This problem becomes especially acute when the particle is immersed in the less viscous fluid and the interface becomes unstable to long-wavelength perturbations (Yih Reference Yih1967). An example of this behaviour is included in figure 14(a), where we plot the deformation of the interface $g$ for different values of the surface tension $Ca^{-1}$, with $\varGamma _{L/2}=\varGamma (L/2)$. As shown in this figure, the maximum value of $g$ increases, creating a sharp gradient near the particle that sharpens further as $Ca^{-1}$ is reduced. In figure 14(b), we plot the evolution of the maximum value of the normal vector in the axial direction $n_{x}$ at the interface $\varGamma$. As the surface tension lowers, $n_{x}$ increases to become $n_x \simeq |\boldsymbol {n}|$ for sufficiently small values of the surface tension. This causes the interface deformation gradient ${\rm d}g/{\rm d}\kern0.7pt x$ to become infinite, anticipating the multi-evaluation of the function $g$ at a range of values of $x$. This behaviour is illustrated in the inset of figure 14(a). The information depicted in this figure suggests, therefore, that the breakdown of the code at $Ca^{-1}=1$, identified in figure 9(b) with a black star, is associated with wave steepening of the interface that would need specific numerical techniques in order to be described.

Figure 14. (a) Interface deformation $g/d$ for different values of the inverse of the capillary number $Ca^{-1}$. (b) Evolution with the inverse of the capillary number of the maximum value of the normal vector at the interface $\varGamma$ in the axial direction $n_{x}$. The rest of the parameters are $Re = 0.05$, $\mu _{1}/\mu _{2}=10$, $Q_{1}/Q_{2}=0.35$ ($\eta =0.60$), $\xi =0.90$, $L=14$.

The wave steeping of the interface illustrated in figure 14 suggests what could be expected beyond the linear regime considered in the present work. In this regard, Boeck et al. (Reference Boeck, Li, López-Pagés, Yecko and Zaleski2007) showed that for a two-layer channel flow without the train of particles, the nonlinear temporal instability associated with the interfacial mode usually results in the formation of ligaments. Valluri et al. (Reference Valluri, Naraigh, Ding and Spelt2010) described the nonlinear spatio-temporal instability in two-layer Poiseuille flow, also illustrating the formation of ligaments. The nonlinear interaction of these ligaments with the particle might trigger an unexpected behaviour that is difficult to anticipate. That question opens the door to new research that is, undoubtedly, of great interest for future study. All the same, different numerical studies in viscosity-stratified flows have concluded that the linear global analysis can anticipate the initial nonlinear performance (see e.g. Tammisola et al. Reference Tammisola, Lundell, Schlatter, Wehrfritz and Söderberg2011; Govindarajan & Sahu Reference Govindarajan and Sahu2014), with a final saturated state that, even if different, is closely related to the linear origin.