Instability of a liquid sheet with viscosity contrast in inertial microfluidics

Abstract Flows at moderate Reynolds numbers in inertial microfluidics enable high throughput and inertial focusing of particles and cells with relevance in biomedical applications. In the present work, we consider a viscosity-stratified three-layer flow in the inertial regime. We investigate the interfacial instability of a liquid sheet surrounded by a density-matched but more viscous fluid in a channel flow. We use linear stability analysis based on the Orr–Sommerfeld equation and direct numerical simulations with the lattice Boltzmann method (LBM) to perform an extensive parameter study. Our aim is to contribute to a controlled droplet production in inertial microfluidics. In the first part, on the linear stability analysis we show that the growth rate of the fastest growing mode \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\xi ^{*}$$\end{document}ξ∗ increases with the Reynolds number \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\text {Re}$$\end{document}Re and that its wavelength \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\lambda ^{*}$$\end{document}λ∗ is always smaller than the channel width w for sufficiently small interfacial tension \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Gamma $$\end{document}Γ. For thin sheets we find the scaling relation \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\xi ^{*} \propto mt^{2.5}_{s}$$\end{document}ξ∗∝mts2.5, where m is viscosity ratio and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$t_{s}$$\end{document}ts the sheet thickness. In contrast, for thicker sheets \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\xi ^{*}$$\end{document}ξ∗ decreases with increasing \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$t_s$$\end{document}ts or m due to the nearby channel walls. Examining the eigenvalue spectra, we identify Yih modes at the interface. In the second part on the LBM simulations, the thin liquid sheet develops two distinct dynamic states: waves traveling along the interface and breakup into droplets with bullet shape. For smaller flow rates and larger sheet thicknesses, we also observe ligament formation and the sheet eventually evolves irregularly. Our work gives some indication how droplet formation can be controlled with a suitable parameter set \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\{\lambda ,t_{s},m,\Gamma ,\text {Re}\}$$\end{document}{λ,ts,m,Γ,Re}. Graphical Abstract Supplementary Information The online version supplementary material available at 10.1140/epje/s10189-021-00147-1.


Introduction
Miniaturized flow devices in the form of a lab-on-achip [1] are often employed for processing fluid flows on the micron scale [2]. Lab-on-a-chip microfluidic applications are used in cell biology [3], chemical synthesis [4], and for manipulating multi-component flows [5], to name but a few. Standard microfluidic devices operate in the Stokes flow regime, while only recently inertial microfluidic platforms have emerged [6]. Their flows at moderate Reynolds numbers enable high throughput and inertial focusing [7,8] in order to develop manipulation techniques for biomedical applications. Motivated by this, a plethora of research has been carried out on inertial microfluidics in the last decade [9][10][11][12][13] including our own studies on the manipulation of soft capsules and solid particles using the inertial lift force [14][15][16].
Recently, instabilities of single-phase flow in different geometries have also been investigated in the inertial regime with the aim to enhance fluid mixing [17,18]. In this article, we use linear stability analysis and lattice Boltzmann simulations to investigate the viscositydriven instability of a multi-component microfluidic flow at finite Reynolds numbers. We let a liquid sheet a e-mail: kuntal.h.patel@tu-berlin.de (corresponding author) stream at the center of a microchannel surrounded by a flowing liquid of larger viscosity and same density and monitor its instability towards modulated interfaces and droplet breakup. Figure 1a shows how the instability develops along the flow direction in a sufficiently long channel. In contrast, in our theoretical investigation we will assume periodic boundary conditions. Such three-layer configurations with two interfaces are commonly encountered in two-phase microfluidic flows [19].
Yih [20] first showed that the fluid-fluid interface in two-layer Couette and Poiseuille flows with viscosity contrast is unstable irrespective of the value of the Reynolds number. Later studies concentrated on interface perturbations with small wavelengths [21,22]. In general, instabilities in viscosity-stratified flows can occur either due to the direct presence of the fluid interface but also due to bounding walls. Boomkamp and Miesen presented an energy budget analysis for the unstable Yih or interface mode, which is triggered by the discontinuity of the shear rate at the interface [23]. Already single-phase flows become unstable at sufficiently large Reynolds numbers due to the presence of bounding walls which cause destabilizing Reynolds stresses. The resulting shear or Tollmien-Schlichting modes also exist in viscosity-stratified flows. Different energy contributions in the energy budget analysis of Boomkamp and Miesen [23] were quantified for twolayer channel flows by Valluri et al. [24] using linear stability analysis. Various nonlinear mechanisms governing the instability of viscosity-stratified flows were reported byÓ Náraigh et al. [25] using three-dimensional direct numerical simulations. Recently, Kalogirou et al. presented the interface dynamics of a thin viscous film adjacent to a wall in a two-layer channel flow with small viscosity contrast [26].
In addition to planar configurations, also core-annular flows in cylindrical channels have been investigated [27][28][29][30][31]. A recent linear stability analysis of core annular flows by Sahu [32] showed the existence of an unstable mode different from Yih and Tollmien-Schlichting modes, which Mohammadi and Smits [33] had also reported earlier in their linear stability analysis of twolayer Couette flows. Redapangu et al. [34] considered a two-phase flow in an inclined channel with the fluidfluid interface of two immiscible fluids normal to the channel walls. In their numerical simulations they then studied how one fluid intrudes the other so that a very irregular three-layer flow arose. For more details on the instability of viscosity-stratified flows, we refer the reader to the comprehensive review article by Govindarajan and Sahu [35].
Viscosity-stratified flows naturally occur in microfluidics when droplets are generated. We review some relevant work. Kurdzinski et al. [36] working in the inertial regime reported different behavior of the central stream in a three-layer configuration of miscible fluids. With increasing Reynolds number they observed a disturbed, a broken, an oscillating, and a stable central stream. In their experiments at low to moderate Reynolds numbers, Hu and Cubaud [37] studied two-layer flows of miscible and immiscible fluids. They observed a linear relation between the wave frequency and the interface velocity. They could describe the dispersion relationship of the interfacial wave using capillary theory in the long-wave regime. Sengupta et al. [38] performed linear stability analysis of miscible fluids with viscosity contrast in rotationally actuated microfluid platforms in order to study fluid mixing. Already in 2007, Guillot et al. [39] investigated the formation of jets and droplets in core-annular microfluidic flows of immiscible fluids at low Reynolds numbers. Microfluidic experiments of Hu and Cubaud [40] addressed the formation of droplets via dripping and jetting in a quadratic microchannel with a less viscous center fluid. Moreover, at relatively high flow rates they only observed waves along the interface. Most recently, Dinh and Cubaud [41] also studied the effect of surface tension continuing work in Ref. [40] but with coaxial microchannels. A detailed review of various active and passive drop generation techniques in microfluidics can be found in Ref. [42]. In particular, external forces are used for active drop generation [43].
In the present work, we investigate the viscositydriven instability of the fluid-fluid interface in a threelayer Poiseuille flow. We perform our work in the framework of inertial microfluidics also with the idea to contribute to controlled droplet generation by tuning fluid properties and geometry appropriately. We perform a thorough parameter study by varying sheet thickness, viscosity ratio, interfacial tension, and flow rate. In particular, we are also interested how the instability develops for different wavelengths of the interfacial perturbation as control parameter. In an experimental setting as sketched in Fig. 1b, an interfacial perturbation with a defined wavelength can be imposed by applying an electric field [44] or an oscillating pressure drop [45] across the inlet of the center fluid. Note, these mechanisms are required only for initial perturbation. The remaining process is solely governed by hydrodynamics. We also note that in this article, we study a twodimensional flow configuration, which can be realized in microchannels with a large aspect ratio.
In the first part of our numerical study we perform a linear stability analysis based on the Orr-Sommerfeld equation and obtain dispersion relations of the unstable interface mode for various combinations of sheet thickness, viscosity ratio, interfacial tension, and Reynolds number. Furthermore, we discuss the structure of the eigenvalue spectra obtained by solving the Orr-Sommerfeld equation in the form of a generalized eigenvalue problem. In the second part, we present direct numerical simulations based on the lattice Boltzmann method (LBM). Starting from an interfacial perturbation of wavelength λ, for thin central sheets we observe two distinct dynamic states of the unstable fluid-fluid interface: traveling interfacial waves and generation of droplets via breakup of the fluid interface. Instead, for thicker sheets we observe interfacial waves but also ligament formation, where the sheet ultimately develops a very irregular shape instead of breaking up into droplets. Thus our study demonstrates the advantage of thin central sheets for controlled droplet genera- Computational setup for studying the instability of a liquid sheet. A microfluidic channel of width w and length λ is considered, where λ is also the wavelength of the interfacial perturbation. The liquid sheet of thickness ts, density ρ, and viscosity μ2 is placed at the center of the microfluidic channel. The rest of the channel is filled with the fluid of the same density and viscosity μ1, where μ1 > μ2. The surface tension of the fluid interface is σ tion. We also perform direct numerical simulations with multi-mode perturbations and show that they result in traveling waves with irregular shape.
The outline of the present article is as follows. In Sect. 2, we present the computational setup and the theory of the linear stability analysis along with the modeling of multi-component flows using the lattice Boltzmann method. We also benchmark our solvers for the linear stability analysis and the lattice Boltzmann simulations. Section 3 presents the numerical results of the linear stability analysis and the direct numerical simulations. We conclude with a summary and final remarks in Sect. 4.

Problem setup and numerical methodology 2.1 Computational domain
The microfluidic setup for the present problem is shown in Fig. 2. We consider a microfluidic channel with noslip boundary condition at the bounding top and bottom walls of the channel and periodic boundary condition in the x-direction. The width and length of the channel are set to w and λ, respectively. λ is also the wavelength of the interfacial perturbation to be studied in the following. To simplify the numerical implementation of the linear stability analysis, we will only consider perturbations symmetric about the channel center with an appropriate symmetry boundary condition at the channel center. However, direct numerical simulations are performed in the entire domain with noslip boundary condition at the top and bottom channel walls. This allows us to capture a possible symmetry breaking of the liquid sheet driven by the instability.
We place the liquid sheet of thickness t s , density ρ, and viscosity μ 2 at the center of the microfluidic channel and fill the remaining part with another more viscous fluid of viscosity μ 1 but matched by density. Solving the Navier-Stokes equation with constant magnitude of the pressure gradient γ and interfacial boundary conditions so that flow velocity U and shear stress are continuous at the two interfaces, we arrive at the flow profile symmetric about the channel center at y = 0: with Finally, to characterize the strength of the flow profile and the surface tension σ of the interface, when introducing dimensionless quantities in our equations, we define the Reynolds number and the inverse capillary number as where U c is a characteristic flow velocity. We choose U c = γw 2 /8μ 1 , which according to Eqs. (2) and (3) is the maximum flow velocity for t s → 0. We also take μ 2 as the viscosity scale and for later use introduce the viscosity ratio m = μ 1 /μ 2 . From now on, we rescale all lengths by the channel width w and, in particular, use the dimensionless parameter t s /w → t s to quantify the sheet thickness. In the following we consider two types of perturbations of the fluid-fluid interface. For the single-mode perturbation the interface is perturbed symmetrically about the channel center with an initial amplitude of 10 −3 . For the multi-mode perturbation we superimpose multiple cosine waves with random phase and amplitude 10 −3 .

Computational linear stability analysis
To set up the linear stability analysis, we introduce a small perturbation of the base solution introduced in Eq. (2) and linearize the mass continuity equation and the Navier-Stokes equations in this perturbation. We name the non-dimensional perturbation for pressure p and we use u, v for the x, y components of the perturbations in the velocity flow field. Using the fact that the base solution satisfies the Navier-Stokes equations, we arrive at the linearized equations written in nondimensional form: where m 1 = m = μ 1 /μ 2 and m 2 = 1 refer to fluids 1 and 2, respectively, and U n are the dimensionless flow profiles introduced in Eqs. (2) and (3). We now formulate the Orr-Sommerfeld equation for both fluid phases, which is the basis for the linear stability analysis. For the velocity components in fluid 1 and 2 we introduce the stream functions ψ n such that the perturbations are written as For stream function and pressure perturbation in both fluids, we make a plane wave ansatz along the channel axis and varying amplitude along the y direction: [ψ n , p n ] = [ϕ n (y), g n (y)]e ik(x−ct) , where k and c represent the real wave number and complex wave speed, respectively. The base flow is linearly unstable if growth rate ξ = Im(kc) is positive. Upon substituting the ansatz functions in the linearized Eqs. (6) and (7) and eliminating the pressure gradient terms, we arrive at the well-known Orr-Sommerfeld equations for fluids 1 and 2: (a n U n − c)(ϕ n − k 2 ϕ n ) − a n U n ϕ n where a 1 = 1 and a 2 = m = μ 1 /μ 2 . The linear stability analysis of the present problem requires the solution of Eq. (10) along with the following boundary conditions: • No-slip boundary condition at the walls at y = ±1/2: • Symmetry boundary condition at y = 0 for symmetric perturbations about the channel center: • The continuity of u at both interfaces: • The continuity of v at both interfaces: • The continuity of shear stress at both interfaces: • The continuity of normal stress at both interfaces: In the present work, we employ the Chebyshev collocation method [46,47] to obtain a numerical solution of the Orr-Sommerfeld equation. In the Chebyshev collocation method ϕ n (y) is approximated as where T n i (η n ) = cos(i cos −1 (η n )) is the i th Chebyshev polynomial of the first kind. The application of the Chebyshev collocation method requires a linear transformation of the global y-coordinate (see Fig. 2) into local coordinates η 1 , η 2 of fluids 1 and 2 such that η 1 , η 2 ∈ [−1, 1]: Derivatives of ϕ n (y) needed in Eqs. (10)-(16) are obtained by differentiating Eq. (17) and using the linear transformation. Once the number of collocation points N 1 and N 2 are fixed, the local coordinates of N n − 3 interior collocation points are Now, we evaluate the Orr-Sommerfeld equation at the interior collocation points in the liquid sheet and the surrounding fluid using the trial solution of Eq. (17) and obtain N 1 +N 2 −6 equations with N 1 +N 2 +2 unknowns. The system of equations is closed by introducing the eight boundary conditions of Eqs. (11)- (16). Finally, the entire problem can be presented in the form of a generalized eigenvalue problem: where [A] and [B] are complex square matrices of dimension N 1 + N 2 + 2 and Φ is a vector of the same dimension. In the present work, we use the eigenvalue solver provided by Matlab to obtain the solution of Eq. (20).

Lattice Boltzmann method for multi-component flows
Direct numerical simulations of multi-component flows [48] require the solution of the mass and momentum conservation laws in each fluid component. Moreover, one also needs to trace the temporal evolution of the interface separating different fluid components. We start with the latter and then address lattice Boltzmann simulations of the fluid flow.

Interface dynamics in lattice Boltzmann simulations
Several techniques exist to access the time evolution of interfaces. On the one hand, they include tracking methods based on Lagrangian coordinates such as front tracking [49] and Lagrangian-Eulerian methods [50]. On the other hand, there also exist capturing methods based on Eulerian coordinates such as the volume-of-fluid [51], the diffuse [52] or sharp [53] interface level-set, and the phase-field [54] method. Interface capturing methods employ a scalar field to distinguish between two fluid components and have the inherent capability to realize the breakup of the fluid-fluid interface. In the present work, we implement the phase-field model based on the conservative Allen-Cahn equation (CACE): where φ ∈ [0, 1] is the phase-field variable and φ = 0 and 1 characterize the two pure fluid components. Furthermore,û is the fluid velocity field, M a mobility, and W quantifies the thickness of the fluid-fluid interface. The Allen-Cahn equation can also be derived from a generalized advection equation for interfaces [55], which was then reformulated as the above continuity equation [56].
In steady state, the phase-field equation (21) is solved by the hyperbolic tangent profile where ζ is the coordinate normal to the fluid-fluid interface [55]. For more details on the implications, which the hyperbolic tangent profile has in descriptions of the interface different from the phase-field method, please refer to Ref. [57]. The Allen-Cahn equation can be discretized using techniques from traditional computational fluid dynamics (CFD) such as finite-volume and finite-element methods. However, for the present work the mesoscale approach based on the lattice-Boltzmann method (LBM) is more appropriate [58][59][60]. In the last two decades the lattice Boltzmann method has emerged as a computationally efficient alternative to traditional CFD techniques, mainly due to the capability of handling multi-physics problems involving soft-matter and fluid flow and the relative ease of parallel implementation [58,59,61].
We implement the phase-field lattice Boltzmann equation proposed by Fakhari et al. [62] on a uniform grid. The method uses the distribution function h i for the phase-field variable φ, where the index i refers to a specific lattice-Boltzmann velocity vector c i . The distribution evolves through successive steps of collision, and advection, where the equilibrium phase-field distribution function obeys τ φ in Eq. (23) represents the relaxation time, which can be related to the mobility as M = τ φ /3. Note that while deriving the Allen-Cahn equation (21) [55,56], the curvature-driven motion is discarded using the counter-term approach proposed by Folch et al. [63]. Hence, M is a purely numerical parameter in our calculations. Note, the present LBM simulations are performed using the D2Q9 lattice with 9 discrete velocity vectors c i = {(0, 0), cyc(±1, 0), (±1, ±1)} and corresponding weights w i = {4/9, 1/9, 1/36}. Now, we can calculate the phase-field variable φ since it is the zeroth moment of the distribution function h i : Finally, we compute the gradient of the phase-field variable ∇φ in second-order accuracy using the isotropic

Lattice Boltzmann model for fluid flow
The lattice Boltzmann equation governing the transport of the distribution function f i for the discrete velocity vectors c i and in the presence of body forces F can be written as [66] ∂f i ∂t where Λ is a collision operator to be introduced below and f eq i is the discrete Maxwell distribution: In general, the density ρ in f eq i is calculated as an average over the two fluid densities ρ 1 and ρ 2 using the phase field, ρ = φ(ρ 1 − ρ 2 ) + ρ 2 , where ρ 1 > ρ 2 . However, in our case we consider density-matched fluids (see Fig. 2). Finally, the last term in Eq. (28) is a forcing scheme proposed by He et al. [67].
The body force F can be decomposed into three parts [66]: wherep is the pressure in lattice units and F b is an external force acting on the fluid. The first term ∇( ρ 3 − p) arises because close to the interface the LB fluid is no longer ideal. The force F s is due to the surface tension of the fluid interface and in the present work following Ref. [68] is modeled as with the chemical potential [69,70] An alternative definition of the surface tension force uses F s = −φ∇μ φ [71,72]. However, Eq. (31) is the most suitable way for implementing F s in a computer code since it avoids the necessity to calculate thirdorder derivatives of φ. Following Fakhari et al. [65], we discretize ∇ 2 φ in Eq. (32) using which is accurate up to second-order.
To improve the numerical stability of our lattice Boltzmann scheme [62,66,73], we introduce a new distribution function Substituting it in Eq. (28) and using Eq. (30) for the body force, we obtain the pressure-evolution lattice Boltzmann equation, which looks the same as Eq. (28) but with a modified body force F g : where in our case ρ 1 = ρ 2 and the respective term in F g vanishes.
In the present work, we employ the multiple-relaxation-time (MRT) collision operator whereŜ is the diagonal relaxation matrix of the following form [62]: Following the Gram-Schmidt procedure [59], it is transformed using the transformation matrix for the D2Q9 lattice: In the relaxation matrixŜ, τ = 3ν + 0.5 is the relaxation time and ν is the kinematic viscosity. To obtain the relaxation time within the fluid-fluid interface, we average the inverse time accordingly using the phase field φ, where τ 1 and τ 2 are the respective relaxation times of the two fluid phases.
As usual, one solves the pressure-evolution lattice Boltzmann equation (36) in two steps [62,74]: where we introduced the pressure distribution functionḡ and the expression forḡ eq i (x, t) is obtained by replacing g i andḡ i by their respective equilibrium distributions g eq i andḡ eq i in the previous formula. Note, the transformation of the pressure distribution function (42) in the collision and the following advection steps is performed to make the scheme explicit in time. [62,65].

Advection step
The directional derivative of the phase-field variable, c i · ∇φ, which appears in the collision step when calculating (c i −û) · F g , is discretized using the mixed finite difference approximation [62,75], The isotropic central difference scheme with the same accuracy is used for computing the modified equilibrium distribution functionḡ eq i (x, t): Finally, the macroscopic quantities such as momentum density ρû and pressurep are calculated from the pressure distribution functionḡ i : In our lattice Boltzmann simulations, we set the interface width W to 4 lattice units and use a bounce-back rule to implement a no-slip boundary condition at the channel walls [59,62].

Simulation parameters
The instability of the liquid sheet is governed by its thickness, the wavelength of the interfacial perturbation, viscosity ratio, Reynolds number, and inverse capillary number. In this subsection, we discuss the values of these parameters, which we take for the linear stability analysis and direct numerical simulations. Before proceeding, we introduce the mean flow velocity U b in addition to the previously defined velocity scale U c = γw 2 /8μ 1 . The velocity scale U c is convenient in theory for making the governing equations nondimensional. However, to connect the present numerical simulations to microfluidic experiments, the mean flow velocity is more appropriate since it directly determines the volumetric flow rate Q = U b w. Thus, for the following discussions we use the new Reynolds number Re and inverse capillary number Γ based on the velocity scale U b . They are readily linked to Re c and Γ c introduced in Eq. (4) using the flow profiles from Eqs. (2) and (3) to calculate U b from the flow rate Q: We will perform an extensive parameter study to obtain a complete overview on the behavior of unstable interfacial perturbations using realistic parameter values. In the linear stability analysis we vary the Reynolds number Re from 10 to 500 in steps of 5 and the reduced wavelength λ of the interfacial perturbation from 0.01 to 5 in steps of 1.25 × 10 −3 , quasi continuously. This will enable us to present the color coded growth rate of an unstable perturbation in the Re-λ plane for two values of the non-dimensional sheet thickness, t s = 0.1 and 0.13, two viscosity ratios, m = 30 and 100 and four values of the inverse capillary number, Γ = {10 −2 , 10 −1 , 1, 10}. We will also investigate larger sheet thicknesses t s up to 0.5 at Re = 100 and Γ = 0.01. In direct numerical simulations on the time evolution of interfacial single-mode perturbations, we study all the combinations arising from λ = {0.25, 0.5, 0.75, 1} for Re = 500, with the same parameter values of Γ, t s , and m, as mentioned before. The specific choices of the wavelength λ and Reynolds number Re is based on the outcome of the linear stability analysis and will become apparent in the next section. Exemplary simulations are also carried out at Re=100. Moreover, we use the same values of t s , m, and Re for direct numerical simulations of multi-mode perturbations using Γ = 10 −2 .
In typical inertial microfluidic applications, where particles are manipulated using the inertial lift force, the channel Reynolds number varies in the range 1 − 100. [76]. We are well in this range, when we define a proper Reynolds number. For example, when we concentrate on the thin liquid sheet, we need to replace in Eq. (47) the channel width w by the sheet thickness t s . Or to describe the outer fluid using the higher viscosity μ 1 is more appropriate. An alternative would be to replace the viscosity μ 2 by averaging the viscosity over the whole channel cross section with the two fluid components, which then gives the Reynolds number Re[t s (m − 1) + 1]m −1 . In the present setup for all the possible combinations of Re, t s , and m this multicomponent Reynolds number is below 100.
In experiments one of the possible ways to realize the setup shown in Fig. 2 is by appropriately merging the outlets of three different microchannels containing the top (fluid 1), middle (fluid 2), and bottom (fluid 1) layers, as sketched in Fig. 1b. In this case the size of the microchannel outlet delivering fluid 2 will determine the sheet thickness t s . The small thickness-towidth ratio t s chosen in our work is achievable in experiments since microchannels with widths ranging from 16 to 500 μm were employed in experiments [77,78]. Furthermore, the viscosity ratio m can roughly range from 1 to 1000, which justifies our choice of m = 30 and 100. For example, pure water (μ = 1 cP) and different silicone oils (μ = 9.3, 97, and 971 cP) cover a range of the viscosity ratio from ca. 10 to 1000 [79]. This makes our choice of high viscosity ratios realistic. Finally, in microfluidic experiments at Re 1 the inverse capillary number Γ can assume values of order 1 or larger, e.g., by choosing ethanol and silicone oil with m = 45 and σ = 0.65 mN/m [40]. In an inertial microfluidic setup, the bulk velocity U b is larger and inverse capillary numbers, Γ ∝ U −1 b , below one are feasible. Hence, also the range of Γ in the present work is justified.

Validation and grid convergence
To validate our solver for the linear stability analysis, we consider a two-layer flow with just one interface in a channel of height 1, where the bottom layer has viscosity μ 1 and thickness h. The viscosity of the less viscous top layer is μ 2 . Note, this problem can be transformed into our original problem (see Fig. 2) by replacing the no-slip boundary condition with the symmetry boundary condition at the top domain boundary. We address two cases with {Re, h, m, Γ} = {500, 0.15, e = 2.718, 0} and {50, 0.3, 30, 0.01}. In both cases, we use 51 collocation points in each fluid layer. The dispersion curves, growth rate ξ versus wave number k, obtained from our solver for both cases are plotted in Fig. 3 along with the dispersion curve reported by Sahu and Govindarajan [80] for the first case and the fastest growing mode from Valluri et al. [81] for the second case. It is evident from Fig. 3 that our results are in very good agreement with published works. In particular, our dispersion curve and the one from Ref. [80] lie exactly on top of each other.
We now focus our attention on the validation and grid convergence of the present lattice-Boltzmann solver for the setup shown in Fig. 2  As shown in Fig. 4, the linear stability analysis and lattice Boltzmann simulations with N = 512 and 1024 produce nearly identical results for the early growth of the interface disturbance. In the first case [see Fig. 4a], lattice Boltzmann simulations suggest that the perturbation grows initially and then reaches a plateau. In contrast, in the second case [see Fig. 4b] we observe the breakup of the liquid sheet. Finally, we conclude that the grid resolution of N = 512 is good enough to capture the interfacial instability in the present problem using lattice Boltzmann simulations.

Results and discussion
We begin with the results obtained from the linear stability analysis for the perturbed liquid sheet. Based on our findings, we then present results from lattice Boltzmann simulations for specific cases with single and multi-mode perturbations.

Parameter study of the interfacial mode
We first look at how our material parameters influence the instability of the liquid sheet when symmetric perturbations with wavelength λ are applied at both interfaces. For this we performed an extensive parameter study and plot in Fig. 5 the color-coded growth rate ξ in the Re-λ plane for specific parameter combinations of the sheet thickness t s , viscosity ratio m, and inverse capillary number Γ. Stable perturbations are indicated by white regions. We find that for all Re-λ combinations, perturbations are unstable for flows with relatively weak interfacial tension [see while stable modes at small λ appear in the columns 2-4, where Γ is larger. Interestingly, in all cases except when larger white regions with stable modes exist, the wavelength associated with the fastest growing mode is always below the channel width, which is one in our reduced units. Moreover, the growth rate of the fastest growing mode always increases with the Reynolds number. The growth rate strongly depends on the wavelength for values smaller than the channel width. For wavelengths larger than the channel width, it reaches a plateau for constant Re.
Increasing the viscosity ratio m or the sheet thickness t s , makes the fluid-fluid interface clearly more unstable, which can be qualitatively described as ξ ∝ m a t b s (a, b > 0). This shows that the present interfacial instability is strongly governed by the viscosity contrast. Most notably, increasing the sheet thickness from 0.1 to 0.13 seems to be a small variation, but it has a significant effect on the growth rate, as rows 3 and 4 in Fig. 5 show. The observed dependence is made more quantitative upon rescaling the growth rate of the fastest growing mode, ξ * → ξ * /mt 2.5 s . Values of exponents a and b in the scaling law are obtained based on the best fit. The rescaled values plotted versus Re nicely fall on a master curve for the same Γ = 0.01 and Re up to ca. 200 (see Fig. 6). The master curve is represented by a cubic polynomial fit and shown as a dashed line. While the curve for t s = 0.13 and m = 100 strongly deviates from the master curve beyond Re ≈ 200, the curves for the other parameters still roughly follow it. In particular, the rescaling with t 2.5 s demonstrates the sensitivity of the most unstable mode on the sheet thickness.
We now focus our attention on investigating the instability over an extensive range of sheet thickness and viscosity ratio keeping Reynolds and inverses capillary numbers fixed. In order to study sheet thicknesses up to values of 0.5 relevant in microfluidic devices, we choose Re = 100 which keeps the relevant multicomponent Reynolds number (see discussion in Sect. 2.4) below 100. In Fig. 7, we plot the color-coded growth rate of the fastest growing mode ξ * in the m-t s plane. Interestingly, in contrast to what we observed so far that with increasing t s the instability of the interface becomes stronger, we now identify a critical sheet thickness t * s beyond which the growth rate decreases. We speculate that this is due to interactions with the channel walls. The critical sheet thickness for different viscosity ratios is marked using filled circles in Fig. 7 and described by the fit curve t * s = 0.7m −0.2675 (solid black line). It quantifies the slow decrease of t s with increasing m.
To explore the large range of sheet thicknesses further, we plot in Fig. 8 the wavelength of the fastest growing mode λ * versus t s for viscosity ratios m = 50 and 100. While for small t s below the critical thickness, the wavelength increases noticeably and roughly linearly, its variation above t * s is much weaker. Similar trends are observed for the other viscosity ratios. Secondly, in Fig. 6 we showed that the growth rate ξ * of the fastest growing mode scales with mt 2.5 s up to Re ≈ 200. In Fig. 9, we plot ξ * /m versus t s for the full range from 0.1 to 0.5 for several m values. On the left of the maximum of ξ * /m, i.e., for small t s , we confirm the scaling, while beyond the maximum it is no longer valid. In Fig. 9, we see that the scaling law starts to deviate as the sheet thickness approaches the critical value. For a fixed viscosity ratio, the value of the critical sheet thickness decreases with the increase in the Reynolds number (not shown here), which results in the early deviation of the scaling law at higher values of Re. This can be the possible cause for the departure of t s = 0.13 m = 100 curve beyond Re=200 in Fig. 6.

Exemplary eigenvalue spectra
So far we have concentrated on the unstable interface mode. For sufficiently small reduced surface tension Γ it always occurs for a nonzero viscosity contrast, i.e., for m = 1 [20,23]. Now, we present the full eigenvalue spectrum of modes, which we obtain for a given wavelength λ, when solving the Orr-Sommerfeld equations  Fig. 5), we choose two wavelengths: λ = 0.25, which is close to the fastest growing mode, and λ = 2.5, which is less unstable. Figure 10 represents the complex wave speed of the perturbation modes introduced in Eq. (9), where Im(c) = ξ/k with k = 2π/λ is proportional to the growth rate and Re(c) is the wave velocity.
In both cases, we observe only one unstable mode [red dot with positive Im(c)], which clearly belongs to the interface becoming unstable. All the other modes are stable [blue dots with negative Im(c)]. For λ = 0.25 and small |Im(c)|, we identify A and P branches, which are reminiscent of the eigenvalue spectrum observed in the single-phase channel flow of a Newtonian fluid [47]. They form the typical upper left and upper right part of a Y-shaped structure in the eigenvalue spectrum called Tollmien-Schlichting modes and center modes, respectively. The modes travel with the respective speeds Re(c)/U 0 → 0 and Re(c)/U 0 → 1. As expected, for stable Poiseuille flow they are also stable.
The unstable interface mode is called Yih mode, named after the scientist who first described these modes quantitatively also in a Poiseuille flow with viscosity contrast [20,23]. For λ = 0.25, the wave speed is close to the flow velocity of the interface, thus in the co-moving frame of the interface, the perturbation wave is static and only the amplitude grows. This changes, when we look at the wavelength λ = 2.5, larger than the channel width. First of all, in the stable modes we do not observe A and P branches. Instead, we recognize two vertical branches: one with wave speed very close to the interface velocity and another branch with a larger wave speed. Interestingly, the unstable interface mode has a wave speed close to the centerline velocity U 0 . Hence, in a co-moving frame of the fluid-fluid interface, the wave is traveling with growing amplitude.

Energy-budget analysis
The unstable interface mode reveals itself also when looking at the energy budget of the disturbance flow [23,24,82]. Here, one considers the kinetic energy of the disturbance flow, which in our case amounts to The base flow solution of the Navier-Stokes equations is stable if the disturbance kinetic energy E(t) > 0 converges to zero with time, meaning the base flow is restored. To decide on this, one calculates the time derivative dE/dt from the Navier-Stokes equations as outlined in "Appendix A." For a flow geometry with interfaces such as the three-layer flow, one obtains with four energy rates on the RHS associated with different mechanisms. The first two terms already exist in a single-phase flow. The rate q dis < 0 denotes the typical energy dissipation through viscous friction, while q Rey results from the convective derivative in the Navier-Stokes equations. It couples Reynolds shear stresses to the base flow and when positive signals the presence of an unstable Tollmien-Schlichting mode. The third and the fourth terms occur when applying the continuity of normal and tangential shear stresses at the fluid interface, respectively. In the energy-budget analysis, the destabilizing disturbance mode can be assigned to one of the four energy rates driving the instability.
To illustrate the energy budget analysis, we apply it to the unstable mode identified in Fig. 10a. The stream function ϕ is illustrated in the inset of Fig. 11. Since the x component of the velocity derives from the derivative of the stream function, u n ∝ dϕ n /dy, one realizes that the disturbance flow is concentrated at the interface. When calculating the different energy rates, we skipped the exponential time factor. Clearly, dissipation due to viscous friction (q dis ) is overcome by energy q tan fed into the system by the tangential stresses at the interface due to the viscosity contrast. Thus the interface mode is unstable. Negligible contributions arise from q Rey > 0 and q nor < 0. From a mathematical view, the destabilizing energy rate q tan is connected to the discontinuity in u n ∝ dϕ n /dy (see "Appendix A"), which has its origin in the viscosity contrast.

Fluid-fluid interface close to the channel wall
As a side remark, we consider the case where the thin liquid sheet is put in contact with one of the channel walls. Thus the current three-layer system transforms into a two-layer configuration with only one interface. We can also treat this geometry with our formalism by changing the relevant boundary condition. Figure 12 Fig. 12 Comparison of the dispersion curves, growth rate ξ versus wavelength λ, for two-and three-layer configurations. Other parameters are indicated in the plot compares the growth rates as a function of λ for the exemplary case of Re = 500, t s = 0.1, m = 100, and Γ = 0.01. While the three-layer configuration shows the behavior as discussed before, we observe that the presence of the wall close to the interface makes the interface mode stable for wavelengths smaller than the channel width, but even at larger λ the growth rate is considerably smaller. Interested readers may refer to Ref. [21] for further discussion on the neutral stability of interface modes in two-layer configurations.

Direct numerical simulations
In this subsection, we present results obtained from direct numerical simulations (DNS) using the lattice Boltzmann method. We carry out DNS for two different values of the Reynolds number: Re = 100 and 500. We begin with results from DNS of flows with Re = 500 since, compared to Re = 100, it requires less computational time to obtain conclusive results for the small sheet thicknesses t s = 0.1 and 0.13. In a second step, Fig. 13 The liquid sheet is either stable against an imposed perturbation of wavelength λ, develops traveling waves along the two interfaces, or breaks up into droplets depending on inverse capillary number Γ, sheet thickness ts, and viscosity ratio m. The Reynolds number is Re = 500 based on the results for Re = 500, we simulate flows at Re = 100 for specific values of the perturbation wavelength and interfacial tension.

Flows with Re = 500
First, we discuss the case of single-mode perturbations and then proceed to multi-mode perturbations. The results of the linear stability analysis in Fig. 5 indicated strong variations of the growth rate ξ for λ ≤ 1. In the following, we thus consider the wavelengths λ = {0.25, 0.5, 0.75, 1} and use the same combinations of t s , m, and Γ as in Fig. 5.
In Fig. 13, we summarize the dynamic states into which the perturbed interface evolves in the λ-Γ plane. Only for the smaller viscosity ratio and the highest reduced surface tension Γ, the interface is stable in agreement with the findings in Fig. 5. Unstable perturbations evolve into traveling waves in the range of smaller wavelengths, which extends to larger values if Γ is increased. In contrast, at larger wavelengths the interface breaks up and droplets form. We now proceed to a detailed discussion of the two dynamic states of propagating interface waves and droplet formation. Figure 14 shows examples of steady-state wave patterns, which we observe for different λ and Γ. For flows with weak interfacial tension (Γ = 0.01) and small wavelength a sawtooth-shaped fluid interface develops. For Γ = 1 and λ = 0.5, we obtain an interface shaped like an hourglass, reminiscent of experimental results reported by d'Olce et al. [83] for core annular flows of two miscible fluids with a viscosity contrast. At our maximal value of the interfacial tension (Γ = 10) and perturbation wavelength equal to the channel width (λ = 1), we see two liquid blobs connected by a thin liquid thread [see Fig. 14c and video V1]. For all these cases the interface patterns preserve the wavelength of the initial perturbation. However, we observe an exception for parameters t s = 0.1, m = 100, and Γ = 10. For this particular setting, the initial perturbation ultimately develops into a secondary disturbance with half the initial wavelength λ = 1 (see video V2). At intermediate times we observe two weak peaks as shown in the top inset of Fig. 15. They oscillate relative to each other in the stream-wise and cross-stream directions. Upon plotting their relative coordinates in both directions, we obtain an inward spiral, as shown in the main plot. In steady state (blue dot in Fig. 15), both peaks attain the same lateral position and a distance of 0.5 along the flow. Thus, the initial perturbation has developed into a an interfacial wave with half the initial wavelength (see the bottom inset in Fig. 15). Finally, in Fig. 16 we show the disturbance velocity field in the y-direction for the wave pattern in Fig. 14c. Since the blob moves to the right, it has to shuffle fluid from the front to the back. The cross-stream motion of the viscous fluid in Fig. 16 can potentially be utilized to enhance heat transfer in liquid-liquid microfluidic flows [84].
In the third dynamic state, the unstable interface breaks up and develops into bullet-shaped droplets of nearly the same aspect ratio around 2. We show this in Fig. 17 for specific parameters but observe similar behavior for other parameter values. However, the transient dynamics towards the final droplet can differ, as we illustrate in Fig. 18 where we show different stages of droplet generation for three representative cases. For λ = 0.5 the breakup of the interface occurs due to the merger of the approaching wave crests. In contrast, for λ = 1 and m = 30, a thin liquid thread of less viscous fluid breaks into multiple satellite droplets, which then coalesce with the primary droplet (see video V3). At the higher viscosity ratio of m = 100, we initially see an additional modulation with half the prescribed wavelength similar to what we described in Fig. 15. Instead of observing satellite droplets, we now obtain a droplet with a slender tail [see Fig. 18(d 3 )]. The tail eventually breaks and merges with the droplet behind it. During this process, the entrapment of a viscous droplet in the less viscous primary droplet occurs, as shown in Fig. 18(e 3 ). Eventually, the entrapped droplet escapes from the tip of the primary droplet. Interestingly, in all three cases the bullet-shaped form of the droplets always develops at approximately the same time t ≈ 20 (see Fig. 18) irrespective of their timing of sheet breakup. In Fig. 19, we summarize the breakup times t b of the liquid sheet for different parameter sets with relatively weak interfacial tension (Γ = 0.01). At constant λ we observe that an increase in either sheet thickness t s or viscosity ratio m leads to an earlier breakup, which coincides with the increase in the growth rate of the interfacial mode. We do not see any substantial change of this behavior at higher values of Γ. The viscosity-driven breakup of the liquid sheet in the present work can be utilized for a controlled droplet production in inertial microfluidics. At constant Re, droplet size and breakup time can be adjusted by selecting an appropriate parameter set {λ, t s , m, Γ}.
In the end we briefly discuss a generic case where we disturb the liquid sheet with a multi-mode pertur- bation. We use the same values of viscosity ratio m and sheet thickness t s from earlier cases of single-mode perturbation and set the inverse capillary number to Γ = 0.01. We choose a channel with length λ = 3 and width w = 1 and disturb the two interfaces symmetrically with a superposition of 15 planar waves: where θ n is a randomly generated phase. For a representative case the time evolution of the maximum displacement A(t) of one interface determined in the lattice-Boltzmann simulations is plotted in Fig. 20 along with the superposition of the 15 time-evolving modes from linear stability analysis. The dashed line indicates the growth rate of the fastest growing mode. The LSA and DNS curves in Fig. 20 indicate three different stages of time evolution. In the beginning, both curves roughly lie on top of each other. The instability grows more slowly in comparison to the fastest growing mode (the dashed black line in Fig. 20). But then it speeds up since small wavelengths become effective and enhance the growth rate; from Fig. 5 we know that waves of small wavelengths grow faster than those with long wavelengths. In particular, in the LSA curve one recognizes that between t ≈ 0.5 and 0.7 the fastest growing mode dominates before all the other modes further speed up the growth. The time evolution enters the second stage at t ≈ 0.7, where both curves start to depart from each other obviously due to nonlinearities, which only the DNS curve takes into account. Interestingly, after t ≈ 1 both curves grow with the same rate. In the third and final stage starting at t ≈ 2 the DNS curve reaches a plateau as indicated in the bottom inset of Fig. 20, while the LSA curve grows further. The snapshot of the fluid-fluid interface from this stage shows a sawtooth-like interface profile as observed earlier in single-mode simulations but now the sequence is irregular (see video V4). We observe similar dynamics for other combinations of the sheet thickness t s and viscosity ratio m.

Flows with Re = 100
We further explore here the possibility to generate droplets, now at a smaller Reynolds number Re = 100. We concentrate on a sheet thickness set to the critical value t * s , where the growth rate is largest, and consider our two viscosity ratios m = 30 and 100. Furthermore, we choose the perturbation wavelength λ = 1 and Γ = 0.01.
For m = 100 and t * s = 0.2, only traveling interfacial waves form and the fluid interface does not break. However, in the second case, for m = 30 and t * s = 0.28, the traveling wave develops ligaments close to the wave crest, as shown in Fig. 21a. Eventually, the ligaments connect back to the liquid sheet and thereby small droplets of the more viscous fluid are enclosed in the sheet [see Fig. 21b]. The formation and reconnection of ligaments continues while the small droplets get advected downstream and eventually merge with the interface, as shown in Fig. 21c (see video V5). This leads to a small local disturbance at the thinner part of the sheet [Fig. 21d]. Eventually, the mirror symmetry of the liquid sheet breaks and the fluid-fluid interface exhibits an irregular topology [ Fig. 21e, f]. A similar phenomenon of ligament formation and irregular evolution has recently been observed in microfluidic experiments of Hu and Cubaud [40].
Interestingly, for thinner sheets with t s = 0.1 and {λ, m} = {1, 30} we do observe the formation of a droplet although the instability is less pronounced compared to the discussed example. This indicates that a strong growth rate does not guarantee droplet formation. Instead, it is important that the sheet thickness is small enough to realize the breakup of the sheet such that droplets are able to form.

Concluding remarks
Lab-on-a-chip devices based on inertial microfluidics operate in the regime of moderate Reynolds numbers where fluid flow is still regular. While conventional inertial microfluidic applications manipulate soft capsules and solid particles using the inertial lift force (see, for example, Ref. [16]), in this article we focused on multicomponent flows in the inertial regime. Concretely, we studied the motion of a liquid sheet at the center of a microchannel surrounded by a flowing liquid of larger viscosity and monitored its instability towards traveling interfacial waves and droplet breakup. Such a configuration is ubiquitous in microfluidic applications involving multi-component flows.
In the first part, we presented dispersion relations for the interfacial mode of the three-layer configuration. The computational linear stability analysis was based on the Orr-Sommerfeld equation. We carried out an extensive parameter study to quantify the effect of sheet thickness t s , viscosity ratio m, interfacial tension Γ, and Reynolds number Re on the growth rate. In particular, we observed that the growth rate of the fastest growing mode ξ * increases with the Reynolds number Re and that its wavelength λ * is always smaller than the channel width w for sufficiently small inverse capillary number Γ. Furthermore, we could quantify the dependence of ξ * on m and t s by the scaling law ξ * ∝ mt 2.5 s in the case of thin sheets, moderate Reynolds numbers, and at weak interfacial tension. In contrast, for thick sheets ξ * decreases with increasing sheet thickness t s or viscosity ratio m due to the presence of the channel walls. Upon examining eigenvalue spectra obtained by solving the generalized eigenvalue problem of Eq. (20), we concluded that Yih modes drive the present interfacial instability [20].
In the second part, we narrowed down the range of material parameters based on the linear stability analysis and presented numerical solutions of the full Navier-Stokes equations using lattice Boltzmann (LBM) simulations. At Re = 500 the thin liquid sheet was stable only for the smaller viscosity ratio and the highest reduced surface tension Γ. Unstable interfaces either evolve into traveling waves or for wavelengths λ ≥ 0.5w they break up and ultimately form droplets when interfacial tension is sufficiently small. The droplets eventually assume a bullet-shaped form. We found that for different perturbation wavelengths λ and viscosity ratios m the temporal evolution towards the droplet differed in the observed interface dynamics and the way the sheet ultimately broke up. However, the droplet formation always ended roughly at the same time. In the introduction, we briefly explained how the wavelength of the interfacial perturbation can be controlled in an experiment.
In practical applications droplet size and breakup time can be tuned by selecting a suitable parameter set {λ, t s , m, Γ} at constant Re. In slower flows at Re = 100 and for thicker liquid sheets, we observed the formation of viscous ligaments at the wave crest, which eventually break up into small droplets and ultimately the shape of the sheet becomes irregular. Thus, we conclude that controlled droplet formation is better achieved with thin sheets. Finally, in direct numerical simulations of multi-mode perturbations we demonstrated how the interface develops a non-regular shape.
In the present work, we investigated the viscositydriven instability of the fluid-fluid interface in the context of inertial microfluidics and demonstrated that it can be exploited for controlled droplet production. Such droplets are utilized in food and pharmaceutical industries [85]. Moreover, microfluidic droplets are also employed to encapsulate biological cells [86] and as chemical reactors in lab-on-a-chip devices [87].

Supplementary information
The online version contains supplementary material available at https://doi.org/ 10.1140/epje/s10189-021-00147-1. and S± indicate the interfaces at y = ±ts/2. We also used the continuity of the y component of the velocity, u2 = u 1 2 = u 2 2 , while the x components are different due to the viscosity contrast.