We consider the genesis and dynamics of interfacial instability in vertical gas-liquid flows, using as a model the two-dimensional channel flow of a thin falling film sheared by counter-current gas. The methodology is linear stability theory (Orr-Sommerfeld analysis) together with direct numerical simulation of the two-phase flow in the case of nonlinear disturbances. We investigate the influence of two main flow parameters on the interfacial dynamics, namely the film thickness and pressure drop applied to drive the gas stream. To make contact with existing studies in the literature, the effect of various density contrasts is also examined. Energy budget analyses based on the Orr-Sommerfeld theory reveal various coexisting unstable modes (interfacial, shear, internal) in the case of high density contrasts, which results in mode coalescence and mode competition, but only one dynamically relevant unstable interfacial mode for low density contrast. A study of absolute and convective instability for low density contrast shows that the system is absolutely unstable for all but two narrow regions of the investigated parameter space. Direct numerical simulations of the same system (low density contrast) show that linear theory holds up remarkably well upon the onset of large-amplitude waves as well as the existence of weakly nonlinear waves. For high density contrasts, corresponding more closely to an air-water-type system, linear stability theory is also successful at determining the most-dominant features in the interfacial wave dynamics at early-to-intermediate times. Nevertheless, the short waves selected by the linear theory undergo secondary instability and the wave train is no longer regular but rather exhibits chaotic motion. The same linear stability theory predicts when the direction of travel of the waves changes — from downwards to upwards. We outline the practical implications of this change in terms of loading and flooding. The change in direction of the wave propagation is represented graphically in terms of a flow map based on the liquid and gas flow rates and the prediction carries over to the nonlinear regime with only a small deviation.

## I. INTRODUCTION

Vertical two-phase flows of thin liquid films sheared by a counter-current gas are prototypical for many technical applications, such as absorption and distillation (using structured packing), evaporation and condensation. In these applications, not only mass and heat transfer but also operational limits are closely linked to the prevailing hydrodynamics. The flow in the two phases, in turn, is largely determined by the interactions between gas and liquid at its interface. Although gas-sheared liquid films have been part of active research on fundamental and practical level for several decades, the rich interfacial dynamics are still not fully understood. These rich dynamics emerge as the liquid interface becomes unstable, leading to the development of waves. Depending on the flow rates, these waves can form highly complex structures, which may give rise to the breakup of the regular wave train, ligament formation, or droplet entrainment.

One of the first to investigate a channel flow with two superimposed fluid layers from a theoretical point was Yih,^{1} who used asymptotic expansion to solve the Orr-Sommerfeld (OS) eigenvalue problem associated with the temporal linear stability in the long wavelength limit for equal densities and layer thicknesses. He found that viscosity stratification alone can cause interfacial instability at arbitrarily small Reynolds numbers, which is also referred to as Yih mechanism. Yiantsios and Higgins^{2} extended the linear stability analysis by accounting for short waves as well as effects due to surface tension and gravity. They observed that, depending on the choice of parameters, the flow is receptive to a short-wave instability at low Reynolds numbers and, moreover, to a shear-mode instability (Tollmien-Schlichting mechanism) for sufficiently large Reynolds numbers. The same approach is extended in the present work to vertical counter-current flow configurations.

To classify the various types of instabilities arising in parallel two-phase flow, Boomkamp and Miesen^{3} analysed by which mechanism energy is transferred from the primary flow to growing disturbances, thereby verifying that both the Yih and the shear mode are important routes to interfacial instability. In fact, a combination of these two mechanisms, also referred to as internal mode, represents a further possible route. We use this approach in the present work for vertical counter-current flows and apply it in cases where a number of modes are active — and where the classification-type of each mode is changing — as the many parameters in the problem are varied throughout the full parameter space. Further classification of parallel flow instability is made by way of the absolute/connective dichotomy,^{4} which we also pursue in the present context of counter-current vertical flows.

In this way, linear stability analysis is shown to be an effective technique to understand the genesis of interfacial instability. However, its validity is limited to disturbances with infinitesimal amplitude. As these disturbances grow, nonlinear effects become important and have to be taken into account. Due to the complexity of nonlinear stability only few general theories exist.^{5} Nonetheless, a variety of modelling techniques have been proposed to describe the development of the interfacial waves up to a finite amplitude. Some of these techniques impose *a priori* assumptions about the wave dynamics under investigation, like long-wave or lubrication approximation, resulting in model equations such as the Kuramoto-Sivashinsky equation or depth-averaged integral equations. A large number of studies on the nonlinear dynamics of interfacial flows are based on these kind of models.^{6–9} However, given that their range of applicability is generally not known in advance, they may at times produce incorrect results (e.g., the erroneous prediction of absolute instability in a falling film, as highlighted by Brevdo *et al.*^{10}). This is especially the case for flow regimes involving large pressure fluctuations and potentially large-amplitude waves, for which there is a major necessity to gain fundamental understanding.

By contrast, weakly nonlinear theories based on either the Stuart-Landau or the Ginzburg-Landau equations dispense with assumptions that cannot be confirmed *a priori* and are capable of matching experimental observations.^{11–13} Such “pure” weakly nonlinear theories are therefore appropriate to complement direct numerical simulations (DNS) of the full Navier-Stokes equations, which, in turn, are guided by basic OS analysis, to study interfacial instability in a rigorous manner. Note that Ginzburg-Landau theory can be applied far from criticality (by criticality we mean that a single mode is barely linearly unstable, and all other modes are linearly stable). Proximity to criticality is sufficient for the theory to apply, but is not necessary. Instead, the necessary condition for the theory to apply is separation of time scales — the growth rate in linear theory of the fundamental mode should be well separated from the time scales of the other modes in the problem.^{14} This is the case in the examples considered throughout this work.

The aforementioned techniques have helped to shed light on the genesis and development of interfacial instability in shear flows and a good understanding of the mechanisms at play has been gained. However, a substantial amount of the available literature is dedicated to horizontal flows or flows down an inclined plane. Flow dynamics specific to a vertical configuration have not received as much attention even though the same methods are applicable. Phenomena related to (partial) liquid flow reversal due to the counter-current gas flow, i.e., (the onset of) flooding, influence the design, optimization as well as operation of technical processes to a great extent and have so far mainly been investigated experimentally.^{15–19} Amongst the theoretical studies of such flows, the works by Trifonov,^{20,21} Tseluiko and Kalliadasis^{22} as well as Dietze and Ruyer-Quil^{9} are worth pointing out. The first three studies consider a liquid film sheared by a turbulent gas stream and analyse the evolution of the interface by solving the governing equations of gas and liquid phase independently (under appropriate assumptions), therefore following the “quasi-laminar” approach of Miles^{23} and Benjamin.^{24} In contrast, Dietze and Ruyer-Quil focussed on a laminar gas phase and included the interfacial coupling of the two phases in a complete manner using the weighted-residual integral boundary-layer (WRIBL) method. Even though the authors also used DNS to validate their low dimensional model, studies on interfacial instability in parallel shear flows employing direct numerical simulation are relatively scarce. To further elucidate the rich dynamics of vertical counter-current gas-liquid flows, we therefore use simulations of the full Navier-Stokes equations together with the semi-analytical methods described above. It is worth noting that applying a low to moderate density contrast to this type of flow is popular in the simulation literature for a number of reasons, such as less complex system dynamics, which helps to pinpoint dynamically relevant mechanisms, but also numerical convenience, see, e.g., the works by Scardovelli and Zaleski,^{25} Boeck *et al.*^{26} and Fuster *et al.*^{27} In this respect, we will make contact to the existing literature by considering low and high density contrasts in our analysis and report on its influence on the system behaviour.

This work is organized as follows. We give a description of the investigated problem and outline the employed methods in Sec. II. Results regarding temporal linear stability of the system are discussed in Sec. III. Spatio-temporal behaviour with respect to absolute/convective instability of the linearized dynamics is presented in Sec. IV. Section V discusses nonlinear wave dynamics. Finally, we give concluding remarks in Sec. VI.

## II. PROBLEM DESCRIPTION AND COMPUTATIONAL METHODS

In this analysis, we consider the dynamics of a gas-liquid flow in a vertical channel, described schematically in Fig. 1. The two continuous phases are separated by an, initially, flat interface. A pressure gradient Δ*p*/Δ*x* > 0 in vertical direction counteracts gravity. We investigate cases in which the balance between gravity and pressure gradient gives rise to counter-current flow, with gas flowing in the direction of decreasing pressure and liquid flowing in the direction of gravity. Figure 1 also shows the development of a linear small-amplitude wave at the interface. Typically, the evolution of interfacial waves depends sensitively on the details of the mean flow.

Both fluid layers exhibit steady, spatially uniform, laminar, and incompressible flow along the vertical rectangular channel. To describe this two-dimensional flow, we use a Cartesian coordinate system, $ x , z $, in which the flat interface is located at *z* = 0 and the confining channel walls are located at *z* = − *d _{L}* and at

*z*=

*d*, respectively. Within these boundaries, the fully developed liquid and gas layer occupy the regions −

_{G}*d*≤

_{L}*z*≤ 0 and 0 ≤

*z*≤

*d*, respectively.

_{G}### A. Base flow and linear stability analysis

With the above mentioned conditions, the Navier-Stokes equations describing the fluid flow in both phases reduce to standard balances between pressure, viscous, and gravitational forces, which are subject to no-slip condition at the channel walls, *z* = − *d _{L}* and

*z*=

*d*, as well as continuity of tangential stress at the interface,

_{G}*z*= 0. To write the governing equations in nondimensional form, we introduce the following dimensionless variables (without tildes) and scalings:

where $x= x , z $ and $u= u , w $ are the coordinate and velocity vector, *H* is the channel height, *V _{p}* is an inertial pressure scale, Δ

*p*/Δ

*x*is a positive pressure gradient,

*τ*is the interfacial shear stress,

_{int}*V*

_{*}is the gas-side interfacial friction velocity, and

*t*denotes time. Further, the following dimensionless parameters arise:

Here, *μ _{j}* is the dynamic viscosity and

*ρ*the density of the respective phase $ j = L , G $, whereas

_{j}*δ*is the relative thickness of the respective fluid layer. The Reynolds numbers

_{j}*Re*,

_{p}*Re*, and

_{g}*Re*

_{τ}, in turn, relate to the applied pressure drop, to gravity and to interfacial shear. A Weber number

*We*accounts for surface tension. With this rescaling, the velocity profile for the undisturbed base flow in nondimensional form reads

Subsection 1 of the Appendix gives the dimensional counterpart to Eq. (3). The ratio of pressure and gravity Reynolds number results in a Froude number, which, in combination with the pressure scale in Eq. (1), represents a measure for the effect of applied pressure drop relative to gravity acting on the gas layer,

To gain insight into the development of small disturbances *η*(** x**,

*t*) centred around the flat interface

*z*= 0, we examine the linear stability of the interface in terms of a standard Orr-Sommerfeld-type analysis (a comprehensive formulation of this analysis can be found in Subsection 2 of the Appendix). In this approach, the governing equations are linearized around the base flow, Eq. (3), and infinitesimally small perturbations in the associated streamfunction $\psi z $ are assumed to have wave-like solutions of the form Ψ(

*x*,

*z*,

*t*) = e

^{i(αx−ωt)}

*ψ*(

*z*), where

*α*=

*α*+ i

_{r}*α*is the complex wavenumber and

_{i}*ω*=

*ω*+ i

_{r}*ω*is the complex angular frequency (both dimensionless). The imaginary parts of these quantities,

_{i}*α*and

_{i}*ω*, denote the spatial and temporal growth rates, respectively. If

_{i}*ω*> 0, the interface is considered temporally unstable and developing (initially) sinusoidal waves propagate with the phase velocity

_{i}*v*=

_{p}*ω*/

_{r}*α*.

_{r}### B. Nonlinear direct numerical simulation

To capture and analyse the development of the flow beyond the linear regime, we use the two-phase flow solver presented by Ó Náraigh *et al.*^{29} for direct numerical simulation of the full Navier-Stokes equations. This in-house solver is level set method^{30} based and uses the continuum surface force (CSF) formulation to model surface tension effects.^{31} In the level set formalism, the governing equations read as

Here, $\varphi x , t $ is the level set function indicating the phase in which a point ** x** lies (liquid phase for

*ϕ*< 0, gas phase for

*ϕ*> 0). Hence, the zero level set, $\varphi x , t =0$, represents the interface $\eta x , t $. The level set function also determines the unit vector $ n \u02c6 $ normal to the interface and the interface curvature

*κ*in Eq. (5d). As the level set function differentiates between the two phases, it is also used to identify the respective density and viscosity through the expressions $\rho = H \u03f5 \varphi +r 1 \u2212 H \u03f5 \varphi $ and $\mu = H \u03f5 \varphi +m 1 \u2212 H \u03f5 \varphi $. The function $ H \u03f5 \varphi $ is a regularised Heaviside function, which is smooth in an interval $ \u2212 \u03f5 , \u03f5 $ around the interface with

*ϵ*set equal to 1.5 times the grid spacing. This interval also supports the regularised delta function $ \delta \u03f5 \varphi =d H \u03f5 \varphi /d\varphi $.

To discretize Eq. (5), we use an isotropic marker-and-cell (MAC) grid with a spacing that resolves the height of the undisturbed liquid film with at least 30 grid points. Additionally, all simulations have been checked for convergence. On the implemented grid, vector quantities are defined at cell faces and scalar quantities are defined at the respective cell centres. A third-order Adams-Bashforth scheme is used to treat the convective derivative, while the momentum fluxes are treated in a flux-conservative fashion employing a (semi-implicit) combination of the Crank-Nicholson and third-order Adams-Bashforth methods.^{32} Pressure and the associated incompressibility of the flow are treated using standard projection method.^{33} We use a combination of Jacobi’s method and successive over-relaxation on a red-black scheme to evaluate the predictor and corrector step. Moreover, a third-order (fifth-order accurate) weighted essentially non-oscillatory (WENO) scheme^{34} is used to advect the level set function *ϕ*, which is subsequently reinitialised applying a Hamilton-Jacobi equation and the algorithm formulated by Russo and Smereka.^{35}

At the domain boundaries, we apply standard no-slip and no-penetration conditions at the confining walls, *z* = 0 and *z* = *H* (note that the coordinate system underlying Eq. (5) is shifted by −*δ _{L}* compared to the one shown in Fig. 1), as well as periodicity in

*x*-direction. The pressure is decomposed as $p= p \u0303 + \Delta p / \Delta x x$, where $ p \u0303 $ satisfies the periodic boundary conditions in

*x*-direction and Δ

*p*/Δ

*x*in the present context is a positive, constant, dimensionless pressure gradient. Solving the standard force balance in both phases (Subsection 1 of the Appendix) numerically gives the initial velocity field (base flow).

The solver is implemented in Fortran 90 using MPI (Message Passing Interface) to decompose the computational domain (Fig. 1) in *x*-direction. This parallelization scheme makes efficient use of the architecture of the UK’s “Advanced Research Computing High End Resource” (ARCHER, http://www.archer.ac.uk) on which we run our high-resolution simulations.

## III. TEMPORAL STABILITY ANALYSIS

In this section, we restrict ourselves to the study of the temporal evolution of an (initially) infinitesimally small perturbation of the liquid interface. This temporal framework provides deep insight into the onset of interfacial waves and, thus, the complex dynamics of vertical films sheared by counter-current gas flows. It further enables us to map flow regimes typical for such systems.

We analyse the temporal stability for the following two distinct cases:

High density contrast: We assume a liquid density of

*ρ*= 1000 kg m_{L}^{−3}, corresponding to a gas-liquid flow typified by an air and water combination. We demonstrate below that the large density contrast leads to a complicated characterization of the instability, consisting of competing and coalescing linearly unstable modes.Low density contrast: To make contact with the existing simulation literature,

^{25–27}we study a low-density-contrast case, with*ρ*= 10 kg m_{L}^{−3}. Studying this test case is of further benefit, since it corresponds to a system without mode competition. This provides a “clean” database of linear stability results, which can be used as an unambiguous benchmark for direct numerical simulations.

For both considered cases, we further assume a liquid dynamic viscosity of *μ _{L}* = 500 ⋅ 10

^{−6}Pa s and a surface tension of

*γ*= 1 ⋅ 10

^{−3}N m

^{−1}, while the dynamic viscosity and density on the gas side are

*μ*= 10 ⋅ 10

_{G}^{−6}Pa s and

*ρ*= 1 kg m

_{G}^{−3}, respectively. The flow is confined by a channel of the height

*H*= 0.01 m. These values result in density ratios of

*r*= 1000 and 10, a viscosity ratio of

*m*= 50 as well as a gravity Reynolds number of

*Re*= 313 and leaves the relative film thickness

_{g}*δ*and the Froude number

_{L}*Fr*, which can be related to liquid and gas flow rates, respectively, as the remaining parameters to determine the two-phase flow. The chosen value of the surface tension is not reflective of true liquid/gas systems (e.g., an air-water system at 20 °C has a surface tension of about 73 ⋅ 10

^{−3}N m

^{−1}; a methanol-air system of about 23 ⋅ 10

^{−3}N m

^{−1}). Instead, the value is chosen so as to throw into sharp relief some particular features of the dispersion relation in the linear theory, which also appear in high-surface-tension cases, but in a less clear-cut fashion. Where appropriate, results at higher surface-tension-values (which are qualitatively similar) are discussed and compared with the presented results for

*γ*= 1 ⋅ 10

^{−3}N m

^{−1}.

Throughout the linear stability analysis, i.e., both density-ratio cases, we consider *δ _{L}* ∈ [0.02, 0.14], whereas the Froude number is varied within the interval $Fr\u2208 1 . 05 , 13 $ for the case of high density contrast. This corresponds to an absolute film thickness

*d*ranging from 0.2 ⋅ 10

_{L}^{−3}m to 1.4 ⋅ 10

^{−3}m and an applied pressure drop Δ

*p*/Δ

*x*in the range of 10.8 Pa m

^{−1}–1657.9 Pa m

^{−1}. Within this parameter space, we determine the temporal dispersion relation $ \omega temp =\omega \alpha r , \alpha i = 0 $ numerically on a grid with Δ

*δ*= 0.005 and Δ

_{L}*Fr*= 0.25. The associated eigenvalue problem, Eq. (A24), is solved for $ \alpha r \u2208 0 . 05 , \alpha c $ with Δ

*α*= 0.05, where

_{r}*α*is the cutoff wavenumber beyond which $ \omega i temp <0$. In the low-density-contrast case, we apply the interval $Fr\u2208 1 . 05 , 1 . 55 $, corresponding to $\Delta p/\Delta x\u2208 10 . 8 , 23 . 6 Pa m \u2212 1 $, with Δ

_{c}*Fr*= 0.025 and determine the dispersion relation for $ \alpha r \u2208 0 . 05 , 15 $ with Δ

*α*= 0.01. For each parameter set $ \delta L , Fr $, we further identify the pair $ \alpha m temp , \omega i temp $ that maximizes

_{r}*ω*as the linearly most unstable mode.

_{i}For both density-ratio cases, our analysis shows that the temporal growth rate of the most unstable mode always attains positive values (Fig. 2), which means that the interface is inherently unstable. However, the same figure demonstrates that the behaviour of the growth rate of the most unstable mode throughout the parameter space is markedly different for the two density-ratio cases. Additional contrasts can be seen from consideration of the wavelength of the most unstable wave, $ \lambda m =2\pi / \alpha m temp $, and its extent relative to the respective film thickness, *λ _{m}*/

*δ*. Figure 3(a) shows the abrupt changes in

_{L}*λ*/

_{m}*δ*in various regions of the parameter space for

_{L}*r*= 1000. Note that the dashed line in Fig. 3 indicates the locus of zero phase velocity (i.e., standing waves) as will be discussed in more detail in Sec. III B.

The rich dynamics observed for the high-density-ratio case can be explained by looking at the dispersion curve of representative points in the parameter space (Fig. 4). It becomes apparent that multiple unstable modes are active, which are coalescing and competing with each other. It is this mode interaction that causes the complex structures presented above. In contrast, the case with low density ratio shows mostly one unstable mode throughout the considered parameter space. Even though further unstable modes are observed, these modes are very weak and appear only occasionally. Hence, for this case, the system behaviour is less intricate (Figs. 2(b) and 3(b)).

### A. Energy budget

To identify the character of these unstable modes and hence the driving force of the instability in these scenarios, we characterize the energy transfer from the base flow by decomposing the disturbance kinetic energy into production and dissipation terms. However, as both presented parameter spaces are quite extensive, we restrict ourselves to the detailed analysis of the representative scenarios listed in Table I. The mode characteristics for all other parameter sets are deduced by inspection.

r . | Scenario . | δ
. _{L} | Fr
. | We
. | Re
. _{L} | Re
. _{G} | $ \alpha m temp $ . | v _{p,OS}
. | $ \omega i , OS temp $ . | $ \omega i , DNS temp $ . |
---|---|---|---|---|---|---|---|---|---|---|

1000 | T1d | 0.13 | 3.00 | 8.829 | 27 261 | 36 323 | 40.05 | 16.47 | 1.6853 | n/a |

T2d | 0.08 | 7.50 | 55.181 | 3 125 | 350 800 | 161.90 | 0.82 | 8.8447 | n/a | |

T3d | 0.08 | 3.00 | 8.829 | 6 166 | 48 332 | 40.85 | 5.97 | 1.7662 | 1.7783 | |

10 | dT4d | 0.08 | 1.157 | 1.313 | 0.384 | 2 142 | 3.99 | 0.05 | 0.3662 | 0.3668 |

dT5s | 1.179 | 1.363 | 0.351 | 2 471 | 4.29 | 0.00 | 0.4669 | 0.4634 | ||

dT6u | 1.201 | 1.415 | 0.317 | 2 806 | 4.59 | −0.05 | 0.5829 | 0.5790 | ||

dT7u | 1.319 | 1.706 | 0.126 | 4 710 | 6.19 | −0.30 | 1.4347 | 1.3806 |

r . | Scenario . | δ
. _{L} | Fr
. | We
. | Re
. _{L} | Re
. _{G} | $ \alpha m temp $ . | v _{p,OS}
. | $ \omega i , OS temp $ . | $ \omega i , DNS temp $ . |
---|---|---|---|---|---|---|---|---|---|---|

1000 | T1d | 0.13 | 3.00 | 8.829 | 27 261 | 36 323 | 40.05 | 16.47 | 1.6853 | n/a |

T2d | 0.08 | 7.50 | 55.181 | 3 125 | 350 800 | 161.90 | 0.82 | 8.8447 | n/a | |

T3d | 0.08 | 3.00 | 8.829 | 6 166 | 48 332 | 40.85 | 5.97 | 1.7662 | 1.7783 | |

10 | dT4d | 0.08 | 1.157 | 1.313 | 0.384 | 2 142 | 3.99 | 0.05 | 0.3662 | 0.3668 |

dT5s | 1.179 | 1.363 | 0.351 | 2 471 | 4.29 | 0.00 | 0.4669 | 0.4634 | ||

dT6u | 1.201 | 1.415 | 0.317 | 2 806 | 4.59 | −0.05 | 0.5829 | 0.5790 | ||

dT7u | 1.319 | 1.706 | 0.126 | 4 710 | 6.19 | −0.30 | 1.4347 | 1.3806 |

Thus, following the approach of Boomkamp and Miesen,^{3} the rate of change of kinetic energy (per unit width in the spanwise direction) in both phases *KIN*_{L,G} can be decomposed into

The terms *DISS*_{L,G} represent energy losses due to viscous dissipation, whereas *REY*_{L,G} denotes wave-induced Reynolds stresses transferring energy between the base flow and the perturbation in the bulk of the two phases. Lastly, *NOR* and *TAN* describe the work done per unit time (also per unit width) by normal and tangential stresses at the interface (for a precise explanation of all of the terms in Eq. (6), see Subsection 2 of Appendix). The energy budgets for the investigated scenarios are given in Table II. The aim here is to use the different positive terms in Eq. (6) together with the shape of the associated streamfunction as a way of classifying the different instability mechanisms at work.

r . | Scenario . | α
. | KIN
. _{L} | KIN
. _{G} | DISS
. _{L} | DISS
. _{G} | REY
. _{L} | REY
. _{G} | NOR
. | TAN
. |
---|---|---|---|---|---|---|---|---|---|---|

1000 | T1d | 1.30 | 0.02 | 0.98 | −0.01 | −1.39 | −0.04 | 2.42 | 0.00 | 0.02 |

T1d | 6.20 | 1.00 | 0.00 | −0.33 | −0.01 | 1.31 | 0.02 | 0.00 | 0.01 | |

T1d | 40.05 | 0.51 | 0.49 | −0.46 | −31.37 | −0.05 | −4.82 | −0.30 | 37.98 | |

T2d | 69.50 | 0.76 | 0.24 | −0.31 | −18.65 | 1.06 | −3.67 | −0.02 | 22.57 | |

10 | dT4d | 3.99 | 0.06 | 0.94 | −1.09 | −13.93 | 0.00 | −1.77 | −0.46 | 18.25 |

dT5s | 4.29 | 0.09 | 0.91 | −1.21 | −12.24 | 0.00 | −1.19 | −0.44 | 16.08 | |

dT6u | 4.59 | 0.12 | 0.88 | −1.36 | −10.67 | 0.00 | −0.72 | −0.42 | 14.17 | |

dT7u | 6.19 | 0.32 | 0.68 | −1.48 | −5.44 | 0.00 | 0.44 | −0.32 | 7.79 |

r . | Scenario . | α
. | KIN
. _{L} | KIN
. _{G} | DISS
. _{L} | DISS
. _{G} | REY
. _{L} | REY
. _{G} | NOR
. | TAN
. |
---|---|---|---|---|---|---|---|---|---|---|

1000 | T1d | 1.30 | 0.02 | 0.98 | −0.01 | −1.39 | −0.04 | 2.42 | 0.00 | 0.02 |

T1d | 6.20 | 1.00 | 0.00 | −0.33 | −0.01 | 1.31 | 0.02 | 0.00 | 0.01 | |

T1d | 40.05 | 0.51 | 0.49 | −0.46 | −31.37 | −0.05 | −4.82 | −0.30 | 37.98 | |

T2d | 69.50 | 0.76 | 0.24 | −0.31 | −18.65 | 1.06 | −3.67 | −0.02 | 22.57 | |

10 | dT4d | 3.99 | 0.06 | 0.94 | −1.09 | −13.93 | 0.00 | −1.77 | −0.46 | 18.25 |

dT5s | 4.29 | 0.09 | 0.91 | −1.21 | −12.24 | 0.00 | −1.19 | −0.44 | 16.08 | |

dT6u | 4.59 | 0.12 | 0.88 | −1.36 | −10.67 | 0.00 | −0.72 | −0.42 | 14.17 | |

dT7u | 6.19 | 0.32 | 0.68 | −1.48 | −5.44 | 0.00 | 0.44 | −0.32 | 7.79 |

In general, positivity of the *TAN* term corresponds either to the viscosity-contrast mechanism of Yih^{1} or to the density contrast. Both viscosity and density contrasts are relevant here. The latter causes a jump in the curvature of the base-state velocity profile at the interface, which in turn induces disturbance shear stresses 𝕋_{xz,j} to satisfy continuity of total shear stress.^{3} The combination of these two mechanisms produce a viscosity-gravity-induced instability which is characterized by a streamfunction compactly supported around the location of the undisturbed interface (hence, this mode is also referred to as an “interfacial mode”). Positivity of *REY _{G}* or

*REY*corresponds to an instability whose streamfunction extends into the gas or liquid layer, respectively, and whose shape is similar to that observed for a Tollmien-Schlichting mode in a single-phase flow. These are referred to herein as shear modes. Occasionally, several source terms in the energy budget are relevant, in which case the shape of the streamfunction is used to classify the instability conclusively.

_{L}Decomposition of the energy budget in scenario T1d (Fig. 4(f)) uncovers three unstable (active) modes. Consideration of the budget in Table II (first three line therein) reveals that the most-dangerous mode is an interfacial one and that the other two are shear modes associated with either fluid layer. This conclusion is supported by the shape of the streamfunction in each case: the streamfunction of the interfacial mode is localized near the undisturbed interface location (Fig. 5(c)), while the streamfunction of the two shear modes extends throughout one fluid layer or the other (Figs. 5(a) and 5(b)). The interfacial mode is active across the entire parameter space of the high-density-contrast case (Fig. 4). The shear mode (first line in Table II) in the gas layer is also active for all values of *δ _{L}* above a minimum Froude number ranging from 1.46 to 2.34 for increasing film thickness. Above this threshold, the mode becomes stronger with increasing gas flow, though thicker films exhibit slightly lower growth rates. The opposite trend is observed for the wavenumber at maximum growth rate, which shifts towards lower values with increasing

*Fr*and, in general, higher values for thicker films. Contrasting the behaviour of the shear mode in the gas layer, the liquid layer shear mode (second line in Table II) is only active for thick enough films $ \delta L > 0 . 067 $. While this threshold value increases with increasing Froude number, the maximum growth rate of this mode drops, thereby rendering it less unstable for higher gas flow rates. Further complexity is added to the system dynamics by the fact that the energy contribution to this particular mode due to tangential stresses at the interface becomes larger for increasing

*Fr*and reaches significant values $ TAN \u2248 1 $. Eventually,

*TAN*becomes the dominant term but with

*REY*still important enough to overcome the restoring effects in both phase. This combination of energy sources driving the instability is also known as internal mode.

_{L}^{3}

A further internal mode emerges at Froude numbers ranging from 5.21 to 5.55 for increasing *δ _{L}*, e.g., Fig. 4(g) or scenario T2d (fourth line in Table II, its streamfunction shown in Fig. 5(d)). As this mode grows stronger with higher gas flow rates, the wavenumber $ \alpha m temp $ at its maximum growth rate shifts towards lower values (compare Figs. 4(g) and 4(j)). Furthermore, this mode starts to coalesce with the dominant interfacial mode, leading to a second “hump” in the dispersion curve of the interfacial mode (Fig. 4(h)). During the formation of this second hump, the dispersion curve broadens with the position of maximum growth rapidly shifting towards higher wavenumbers (shorter wavelengths) as

*Fr*increases (Fig. 3(a)). Yet, at the same time, the growth rate increases only slightly, resulting in the plateau-like region at moderate Froude numbers observed in Fig. 2(a). For even higher

*Fr*, the coalescing internal mode undergoes a change in identity itself and forms a second active interfacial mode, thereby effectively “splitting” the initial dispersion curve into two separate ones (Figs. 4(k) and 4(l)). Although mode coalescence occurs throughout the entire range of investigated film thicknesses, splitting of the dispersion curve is not observed for low values of

*δ*(Fig. 4(j)). It is further worth mentioning that additional modes can become unstable, which are, however, mostly temporarily active and thus of minor significance.

_{L}The multitude of active and coexisting modes leads not only to mode coalescence but also to a certain amount of mode competition. In the high-density-contrast case, this mode competition mainly occurs at low Froude numbers, where the growth rate associated with the interfacial mode is still comparatively low. The shear modes at low and high values of *δ _{L}*, on the other hand, exhibit stronger growth and therefore constitute the dominant mode (Figs. 4(a) and 4(c)). As

*Fr*increases, the interfacial mode picks up strength and supersedes the shear modes as the dominating mode (Figs. 4(d) and 4(f)), which results in a jump of the wavenumber $ \alpha m temp $ associated with the maximum growth rate. This jump corresponds exactly to the sharp jumps in the contour plot in Fig. 3(a) at both low and high values of

*δ*.

_{L}As mentioned above, the low-density-ratio case exhibits mainly one unstable mode, which is due to the density and viscosity contrasts of the two fluids (fourth line in Table II). Although both differences account for energy transferred towards the disturbed flow, the contribution related to the viscosity contrast dominates. Hence, this mechanism is, in general, consistent with the so-called Yih mode.^{1} It is further apparent that the relative fraction of kinetic energy associated with the liquid phase increases with increasing Froude number. This rise, together with an enhanced energy dissipation, can be linked to more agitation in the liquid film as we will show in Section V. The amount of energy dissipated in the gas phase, on the other hand, seems to drop, which is counter-intuitive for an increased *Fr*. Yet, dissipation in the gas does increase in absolute terms but at a slower rate as the total kinetic energy. That, in turn, leads to the seemingly decreasing rate of dissipation in the gas phase. The same effect can be seen for *NOR* and *TAN*.

Another result worth noting is the change in sign of *REY _{G}* as

*Fr*increases, turning wave-induced Reynolds stresses from an additional dissipative energy “sink” into an energy “source” (scenario dT7u). This, in turn, can be explained by an unstable Tollmien-Schlichting mode appearing in the gas stream that delivers energy to the interfacial instability, thereby suggesting a transition to turbulence in the bulk of the gas phase. Therefore, this particular scenario can conceptually not be regarded as a “pure” Yih-type instability any more but tends towards an internal mode. This positive contribution of wave-induced Reynolds stresses to the instability above a certain Froude number occurs throughout the entire parameter space but the threshold decreases for thicker liquid films. However, it has to be emphasized that the tangential stresses doing work at the interface are the dominant driving force of the instability in all presented scenarios for the low-density-ratio case (lower half of Table II). Across the entire parameter space, the maximum growth rate $ \omega i temp $ of this interfacial mode increases with increasing

*δ*and

_{L}*Fr*(Fig. 2(b)), whereas the associated wavenumber $ \alpha m temp $ increases with increasing Froude number but decreases for thicker films, which is in agreement with results presented by Dietze and Ruyer-Quil.

^{9}Overall, the system is predominantly receptive to long-wave instability under low-density-ratio conditions (Fig. 3(b)).

### B. Flow regimes

Figure 6 shows the phase velocity *v _{p}* of the fastest growing wave developing on the interface for both density-ratio cases. It becomes apparent that the parameter space is divided into two regimes: one in which developing waves exhibit a positive phase velocity and another in which the phase velocity is negative. With respect to the chosen coordinate system (Fig. 1), these regions correspond to waves propagating downwards and upwards, respectively. The vanishing of

*v*at the demarcation between these regimes (dashed lines in Fig. 6) relates therefore to a standing wave and is herein referred to as the

_{p}*loading curve*. This demarcation between downward- and upward-travelling waves is important from a practical point of view because it is related to the onset of

*flooding*, which is understood as the partial upward flow of the liquid phase

^{15}(flooding itself is regarded as the zero net flow of the liquid

^{36}). Because the demarcation refers to the direction of travel of waves only, it does not by itself imply flooding. Yet, it implies the possibility of flooding, since upward-travelling nonlinear waves may form complicated nonlinear structures leading to ligaments and droplet entrainment, of which the latter will promote transport of the liquid in the upwards direction.

The flow map for the case with high density ratio (Fig. 6(a)) reflects the complex interfacial dynamics described above. In regions with mode competition and a changing dominant mode, the phase velocity changes drastically due to the jump of the wavenumber $ \alpha m temp $ associated with the linearly most unstable mode. For low Froude numbers, this leads to an “island” of upward-travelling disturbances amidst the sea of downward-travelling waves at the thin-film end of the parameter space (Fig. 6(b)), while a region of comparatively slow waves occurs at the high-film-thickness end. Once *v _{p}* < 0 for the most unstable mode, we can expect that exponentially growing wave to travel upwards and overwhelm all other signals. In contrast, the phase velocity changes smoothly in the low-density-ratio case due to the absence of mode competition (Fig. 6(c)).

To illustrate the different factors that determine the shape of the loading curve, Fig. 6 also shows the curve of zero interfacial velocity (dotted line) obtained from the undisturbed base flow. In the high-density-ratio case, the two curves coincide (with exception of the island region) and information from the base state alone is sufficient to determine the loading curve (Fig. 6(a)). This can be explained by looking at the expression for the phase velocity, which can be decomposed as $ v p = u 0 , int + u 1 \alpha $, where *u*_{0,int} denotes the interfacial velocity of the undisturbed base flow. At the linearly most unstable mode, the difference |*v _{p}* −

*u*

_{0,int}| decreases with increasing

*r*(Fig. 7), except at very large density ratios

*r*> 1000, whereupon the difference saturates at a small residual value. Thus, the phase speed of the (linearly) most unstable interfacial mode can be well approximated by

*u*

_{0,int}for high density ratios. Conversely, at low density ratios, determination of the phase speed requires information not only from the base state but also from the full eigenvalue problem, as the loading curve is substantially modified by the latter (Fig. 6(c)).

We have carried out further investigations at higher surface-tension values to understand qualitatively the behaviour of the linearized system at conditions more indicative of an air-water-type system. The two broad behaviour types obtained in the present study carry over to higher values of surface tension. Not surprisingly, the shear modes are virtually unaffected by an order-of-magnitude change in surface tension. On the other hand, the growth rate of the interfacial mode as well as the corresponding wavenumber are reduced by the same increase in surface tension. For the *r* = 10 case, this leads to a simple “shift” in the dispersion relation of the interfacial mode to longer wavelengths. In contrast, higher surface tension promotes further mode competition between the interfacial and shear modes for the high-density-ratio case, leading to a more complicated flow map (not shown) than the one presented in Fig. 6(a); with larger “islands” of negative phase velocity and relatively slow waves (bottom corners of the flow pattern map) due to mode competition and the increasing prominence of the shear modes.

## IV. ABSOLUTE AND CONVECTIVE INSTABILITY OF THE LINEAR DYNAMICS

Complementing the temporal analysis of the interfacial instability, we study the response of the system to a small localized impulsive disturbance imposed on the liquid interface. If the disturbance grows in-situ, the system is absolutely unstable. Conversely, if the disturbance grows but only as it is convected away from the source, the system is called convectively unstable. In this section, we restrict ourselves to the study of the low-density-ratio case. The reason for this is the complex system behaviour associated with the high-density-ratio case caused by mode competition and mode coalescence, which complicates the use of the methods outlined below. This is discussed further at the end of the present section.

To determine the spatio-temporal nature of the instability, the complex dispersion relation $D \alpha , \omega =0$, which is obtained by solving the eigenvalue problem of Eq. (A18) numerically for a range of complex wavenumbers *α* = *α _{r}* + i

*α*, has to be evaluated against conditions essential for absolute instability as outlined by Huerre and Monkewitz

_{i}^{4}: (i) a positive imaginary part of the angular frequency $ \omega i , 0 \u2254 \omega i \alpha S >0$ at a saddle point

*α*in the complex

_{S}*α*-plane, where

*α*solves d

_{S}*ω*/d

*α*= 0, forms the necessary condition; (ii) to satisfy the sufficient condition, spatial branches $ \alpha \xb1 \omega $ that originate from opposite halves of the

*α*-plane have to coalesce at the saddle point

*α*, forming a genuine pinch point (this coalescence corresponds to the formation of a branch point/cusp at

_{S}*ω*

_{i,0}in the complex

*ω*-plane

^{37}). Meeting both conditions will result in growth of the disturbance at its source with the

*absolute growth rate*

*ω*

_{i,0}.

Although the described procedure seems straightforward, inspection of the results of a spatio-temporal Orr-Sommerfeld (ST-OS) analysis (Fig. 8) requires great care in order to avoid misinterpretation due to the complexity of the dispersion relation that can arise from the multivalued nature of the eigenvalue problem as well as specifics of the investigated problem, such as applied boundary conditions or multiple unstable temporal modes. Hence, to conclude the character of the instability correctly, it is necessary to examine the global topography of *ω _{i}* in the complex

*α*-plane for each relevant set of parameters. As we are interested in identifying convective/absolute instability (C/A) transition, the described method is not practical given the large parameter space considered herein. Instead, we use an approximation technique (Quadratic Approximation - QA) based on an analytical connection between temporal and spatio-temporal frequencies presented by Ó Náraigh and Spelt.

^{38}

The approximation technique is based on the following identity for the analytic continuation of the growth rate *ω*_{i} into the complex plane:

where *c*_{g} = d*ω*_{r}/d*α*_{r} is the group velocity in a standard temporal analysis and $ \omega i temp $ is the temporal growth rate in the same. Equation (7) is a consequence of the Cauchy-Riemann conditions on *ω*(*α*) viewed as an analytic function on an appropriate open subset in the complex plane.^{38} We truncate this series at quadratic order in *α*_{i} to yield

We apply the conditions for a saddle point d*ω*/d*α* = 0 to Eq. (8). By Cauchy-Riemann, this implies that ∂*ω*_{i}/∂*α*_{r} = ∂*ω*_{i}/∂*α*_{i} = 0 for a saddle point. Hence, in the quadratic approximation, for a saddle point, the following simultaneous equations are satisfied:

Solution of these equations for *ω*_{i} and *ω*_{r} yields the quadratic approximation for the location of the saddle point *α _{S}* and hence, an estimate for the absolute growth rate

*ω*

_{i0}=

*ω*

_{i}(

*α*) for a given set of flow parameters. Note that all of these estimates are based on results from a temporal linear stability analysis only, meaning that it is straightforward to make these estimates using standard temporal linear stability theory, and one apparently circumvents the pitfalls associated with the full spatio-temporal linear stability analysis outlined above. However, one must be cautious in applying this methodology, as the quadratic approximation is, strictly speaking, only valid inside a disc of convergence with radius

_{S}*R*, where

*R*is the minimum distance from the centre of the complex Taylor series $ \alpha r = \alpha m temp , 0 $ to the nearest singularity of $\omega \alpha $.

Like in the purely temporal framework, we study several scenarios in more detail (Table III). Using spatio-temporal Orr-Sommerfeld analysis and quadratic approximation as outlined above, we determine the saddle point *α _{S}* as well as the corresponding absolute growth rate

*ω*

_{i0}for all listed scenarios and further perform direct numerical simulations for the scenarios ST1 and ST4. In the following, we want to compare and discuss the obtained results for scenario ST4 in more detail.

Scenario . | δ
. _{L} | Fr
. | Re
. _{p} | We
. | Method . | α _{r,S}
. | α _{i,S}
. | ω _{i0}
. |
---|---|---|---|---|---|---|---|---|

ST1 | 0.08 | 1.075 | 337 | 1.134 | QA | 3.34 | −2.24 | −0.0877 |

DNS | n/a | n/a | <0 | |||||

ST2 | 0.08 | 1.157 | 362 | 1.313 | QA | 4.01 | 0.24 | 0.3632 |

ST-OS | 4.02 | 0.24 | 0.3632 | |||||

ST3 | 0.08 | 1.179 | 369 | 1.363 | QA | 4.35 | 0.66 | 0.4422 |

ST-OS | 4.43 | 0.65 | 0.4423 | |||||

ST4 | 0.08 | 1.201 | 376 | 1.415 | QA | 4.72 | 1.02 | 0.5205 |

ST-OS | 4.87 | 0.97 | 0.5221 | |||||

DNS | n/a | n/a | 0.5269 |

Scenario . | δ
. _{L} | Fr
. | Re
. _{p} | We
. | Method . | α _{r,S}
. | α _{i,S}
. | ω _{i0}
. |
---|---|---|---|---|---|---|---|---|

ST1 | 0.08 | 1.075 | 337 | 1.134 | QA | 3.34 | −2.24 | −0.0877 |

DNS | n/a | n/a | <0 | |||||

ST2 | 0.08 | 1.157 | 362 | 1.313 | QA | 4.01 | 0.24 | 0.3632 |

ST-OS | 4.02 | 0.24 | 0.3632 | |||||

ST3 | 0.08 | 1.179 | 369 | 1.363 | QA | 4.35 | 0.66 | 0.4422 |

ST-OS | 4.43 | 0.65 | 0.4423 | |||||

ST4 | 0.08 | 1.201 | 376 | 1.415 | QA | 4.72 | 1.02 | 0.5205 |

ST-OS | 4.87 | 0.97 | 0.5221 | |||||

DNS | n/a | n/a | 0.5269 |

As shown earlier already, the global topography of *ω _{i}* in the complex

*α*-plane (Fig. 8) is rather complex. First, the dispersion relation contains two saddle points, of which both may be dynamically relevant. The confinement of the flow by the channel walls has, furthermore, created a discrete pole (not shown) on the imaginary axis,

^{38}$ \alpha r , \alpha i = 0 , 3 . 34 $, which has implications on the character of the saddle point closer to that particular singularity. Lastly, the multivalued nature of the dispersion relation becomes apparent by the branch cut just below the real axis. Although these features make the final characterisation of the instability more difficult, the saddle point at $ \alpha r , \alpha i = 4 . 87 , 0 . 97 $ clearly appears as a result of the coalescence of spatial branches emanating from opposite half-planes. It is therefore also a pinch point and contributes to spatio-temporal growth at a rate of

*ω*

_{i0}= 0.5221. The positive value of

*α*

_{i}indicates that spatial growth of the disturbance happens for

*x*< 0 (upwards), which is confirmed by DNS (Fig. 9(a)). On the other hand, the saddle point at $ \alpha r , \alpha i = 0 . 62 , 0 . 65 $ is not a pinch point. In fact, the spatial branch $ \alpha + \omega $ (above the saddle point) is a closed curve, whereas the $ \alpha \u2212 \omega $ branch does not (entirely) originate from the negative half-plane. Hence, this saddle point is dynamically irrelevant regarding absolute instability.

As mentioned before, the singularity closest to the position of the temporally most unstable mode $ \alpha r = \alpha m temp , 0 $ determines the disc of convergence in which the quadratic approximation can be applied with confidence. In scenario ST4, the confinement pole at $ \alpha r , \alpha i = 0 , 3 . 34 $ limits the outermost radius *R* of this disc to about 5.68. Equation (8) is thus convergent across the relevant section of the complex *α*-plane depicted in Fig. 8. The approximated position of the pinching saddle point as well as the corresponding growth rate agree well with spatio-temporal OS analysis, see Table III.

We further compare the impulse response for the parameters of scenario ST4 captured by DNS against linear theory (ST-OS). Use of the DNS with periodic boundary conditions is justified, provided a long channel is used, such that a comparison between the theory and the numerics can be made before the onset of finite-size effects. We carefully check in what follows that our findings are robust to such effects. In contrast to simulations within the purely temporal framework, a Gaussian pulse, centred around *x*_{0} = 16, with a height of 1 ⋅ 10^{−3} and a standard deviation of 1 ⋅ 10^{−1} is applied to the otherwise flat interface position as initial condition to study the spatio-temporal nature of the flow. Growth of this perturbation at its source then constitutes absolute instability. However, the developing disturbances will inevitably contaminate the pulse source due to the implemented streamwise periodic boundary conditions. To delay this contamination sufficiently, the channel length is set to *L _{x}* = 20. Using the pulse norm

where *w* denotes the wall-normal velocity, the space-time plot in Fig. 9(a) shows the temporal evolution of the perturbations along the streamwise coordinate. It can be seen that, starting from its initial position at *x*_{0} = 16, the pulse is convected upwards, which is in accordance with the results from linear theory mentioned above. This upward motion triggers a pressure disturbance moving ahead of the pulse, which is visible on the left-hand side of Fig. 9(a). However, this “shock” decays sufficiently before it re-enters the computational domain, thus avoiding early contamination of the instability source. Further information from this plot is extracted in Fig. 9(b), where the growth of the instability at its source is given. After an initial transient period, the growth rate of the instability follows the value predicted by ST-OS and QA closely, therefore validating the DNS results. Thereafter, a second, more violent, pressure shock develops as a result of two merging waves (Fig. 10) at *t* ≈ 11.5 and travels upwards to eventually contaminate the pulse source at around *t* = 14.0 (disturbances on the right-hand side of Fig. 9(a)).

The excellent agreement that has been established between linear theory (ST-OS) and quadratic approximation in ST4 can also be observed for the scenarios ST2 and ST3 (Table III). For scenario ST1, we were not able to determine the absolute growth rate of the system by means of OS analysis and saddle point method due to the multivalued nature of the associated eigenvalue problem. In fact, the dynamically relevant saddle point appears in the lower half-plane below the branch cut arising near the real axis. To analyse this saddle point, a laborious reconstruction of the corresponding part of the Riemann surface from discrete eigenvalues would have to be carried out. However, we avoid this procedure by using the quadratic approximation, which indicates scenario ST1 to be convectively unstable. A further direct numerical simulation confirms this result (not shown).

Given these results, we apply the QA to identify the spatio-temporal nature of the system throughout the wide parameter space considered in the low-density-ratio case. The calculated growth rates *ω*_{i,0} are shown in Fig. 11, where the dashed lines demarcate the C/A boundaries. It becomes apparent that the system is absolutely unstable in almost the entire domain with exception of two narrow bands at the low-Froude-number and low-film-thickness end of the parameter range, respectively. In view of these results and those of Sec. III B, it appears that C/A transition and the onset of upward-travelling waves are not closely related for the case of density ratio *r* = 10.

Regarding the accuracy of these results, it has to be mentioned that the relative error $ \omega QA \u2212 \omega ST-OS / \omega ST-OS $ between QA and linear theory increases with decreasing film thickness and increasing Froude number. At the same time, the quadratic approximation generally underestimates the “true” value of *ω*_{i,0} as obtained by ST-OS, which leads to a larger regime of absolute instability (towards thin films) than displayed in the above figure.

As mentioned at the beginning of this section, the mode coalescence and mode competition encountered in the case with high density ratio (Fig. 4) makes it difficult to identify the nature of the instability in a spatio-temporal framework with both the present semi-analytical methods. To study this manifestation of the instability, alternative techniques, such as linearized DNS^{39} (wherein the linearized equations of motion are solved as a Cauchy problem, but still within the framework of the Orr-Sommerfeld linear operators) or construction of series solutions of the underlying stability problem,^{40–42} might be more suitable. A detailed investigation of this particular regime, i.e., high density contrast, will be left to future work.

## V. NONLINEAR WAVE DYNAMICS

As temporal linear stability analysis permits the analysis of infinitesimally small interface perturbations only, we carry out direct numerical simulations to study the evolution of these disturbances up to finite amplitudes for both low and high density contrast. To that end, we use the in-house solver described in Sec. II B with streamwise-periodic boundary conditions and initial interface location

where *A*_{0} is some small initial amplitude (herein *A*_{0} = 1 ⋅ 10^{−3}), *N* is the number of linearly unstable modes initialized, $ \alpha n =n 2 \pi / L x $ is a wavenumber in streamwise direction and *φ _{n}* is a random phase. Even though periodic boundary conditions do not reflect the behaviour of a real system, it is appropriate to consider this setup as it allows for a rigorous comparison of DNS results with linear theory as well as unambiguous results of the Fourier transform taken of the interfacial wave at finite times.

^{29}

### A. Low density contrast

To allow for rigorous comparison with linear theory, we perform direct numerical simulations for each flow regime of the low-density-contrast case (lower part of Table I). For that purpose, we set the wavenumber *α _{n}* in Eq. (11) equal to the linearly most unstable mode $ \alpha m temp $ of the respective scenario. Figure 12(a) shows the

*L*

^{2}-norm of the wall-normal velocity perturbation

*w*over time for scenario dT4d and it can be seen that the disturbance grows with the rate $ \omega i temp $ predicted by Orr-Sommerfeld theory. Comparing the streamfunction $ \psi j z $ at the crest of the developing wave also yields excellent agreement (Fig. 12(b)). Also the other scenarios follow the theoretical predictions in terms of growth rate in an equally convincing fashion (Table I). It can further be seen in Fig. 12(a) that the wave initially grows exponentially before nonlinear effects gain importance at about

*t*= 8.0. Beyond that point, wave growth slows down and eventually saturates, leading to a steady nonlinear wave of constant amplitude travelling along the interface. The nonlinear saturation in Fig. 12 occurs almost immediately after the period of exponential growth ends. This is because of the large separation in time scales between the linearly unstable fundamental and the linearly stable harmonics. Because the harmonics are damped in the linear theory, they become rapidly slaved to the fundamental, which in turn leads to prompt saturation of the same.

Figure 13(a) depicts the early stage evolution of the interface up to saturation for the scenario with downward-travelling wave, dT4d. It can be appreciated that the wave developing on the interface moves indeed downwards as is predicted by linear theory. More insight into the flow features of this scenario is given by Fig. 14(a). Plotted in a fixed frame of reference, a large anticlockwise rotating recirculation zone, positioned on the “downwind” side of the wave, is revealed in the gas phase. This vortex is sandwiched between the main flow of the gas and a thin layer of gas along the entire length of the interface, which is dragged downwards by the liquid due to interfacial shear stress. The liquid phase itself does not exhibit any major disturbances.

According to linear theory, the loading point of the system, with *δ _{L}* = 0.08, is reached by increasing the Froude number to

*Fr*= 1.179 (scenario dT5s, see Fig. 6(c)). At this point, the phase velocity of the wave vanishes as gravitational and lift forces on the wave hump balance each other, resulting in a standing wave. The evolution of this wave in its early stage is depicted in Fig. 13(b). It can be apprehended that the initial upward motion of the wave slows down after it is saturated and eventually comes to a complete halt around the position last indicated in the plot. As the system evolves further, the saturated wave begins to move downwards with the bulk of the liquid (not shown here). This is attributed to nonlinear effects which lead to an anticipated deviation from the system behaviour predicted by linear theory. However, it is worth mentioning that the deviation is relatively small (

*v*= 0.0103) and, hence, linear theory gives a good indication of the development and behaviour of the liquid interface for the loading scenario (dT5s). Moreover, the increased Froude number leads to a decrease of the mean liquid streamwise velocity, which, in turn, reduces the extent of the thin layer of downward gas flow along the liquid interface, causing the gas-side vortex to move closer to the interface. On the liquid side, contrary to the scenario with downward-travelling wave, a small anticlockwise rotating vortex occurs in the wave body, which indicates a negative (upwards) streamwise velocity of the liquid at the interface in the vicinity of the wave crest and is in agreement with the results of Trifonov.

_{p}^{20}However, the bulk of the liquid, which refers to the region of the film between channel wall and wave trough, remains largely unaffected.

By increasing the Froude number beyond the loading point, as in scenario dT6u, shear on the wave body due to the enhanced gas flow overcomes the gravitational forces, causing interfacial waves to travel upwards (Fig. 13(c)). Similar to the previous scenarios, the developing wave saturates and forms a coherent structure. However, with higher Froude numbers, the general trend of increasing growth rate and decreasing amplitude becomes apparent. Furthermore, the change of direction in which the interfacial wave propagates leads to a more agitated liquid film, especially in the wave body (Fig. 14(c)). Compared to scenario dT5s, an extended region of liquid near the wave crest experiences upward movement with the wave body and becomes recirculated in an enlarged eddy. The bulk of the liquid, on the other hand, essentially remains unaltered also in this scenario.

After having discussed the influence of the Froude number, i.e., gas flow rate, on the gas-liquid flow in a qualitative manner, we now want to turn our attention to a qualitative description of the nonlinear wave dynamics. To understand the genesis of these, we compute the amplitude of the Fourier modes of the interface height for each time step. The amplitudes of the first three harmonics are shown in Fig. 15(a). Matching the rate given by Orr-Sommerfeld theory, the first harmonic ($ \alpha 1 = \alpha m temp =3.99$) grows exponentially fast at the beginning. At the same time, the higher harmonics (*α*_{2} = 7.98, *α*_{3} = 11.97) are linearly stable, as predicted by the same theory. However, these harmonics eventually also undergo exponential growth, whereat the *n ^{th}* harmonic grows at a rate $n \omega i temp \alpha 1 $. As time goes by, the growth rate of all harmonics decreases and the amplitudes saturate simultaneously. It is thus apparent that the dynamics of the higher harmonics are “slaved” to the first harmonic. This temporal development of the amplitudes is perfectly consistent with weakly nonlinear theory.

^{11}Using this theory, one can also model the temporal development of the first harmonic using the Stuart-Landau (SL) equation.

^{43,44}Including up to quintic-order terms, the absolute value of the finite amplitude follows

where *μ*_{1} is twice the temporal growth rate $2 \omega i temp $ and *μ*_{2,3} is twice the real part of the Landau coefficients.^{5} These coefficients have been fitted numerically for the different scenarios discussed above and are listed in Table IV together with the root mean square deviation (RMSD) and R^{2}-value of the best fit, illustrating the excellent agreement between theory and direct numerical simulations. The negative value of *μ*_{2} confirms that the instability is saturating in all scenarios. For scenario dT4d, the relatively large value of *μ*_{3} points towards a longer transition phase, with decreasing growth of the wave amplitude, leading from the initial stage of exponential growth to saturation in this scenario. In contrast, the wave in case dT5s, dT6u, and dT7u reaches its saturated state faster, which is reflected by the vanishing of *μ*_{3} in Eq. (12). The comparatively low R^{2}-value in scenario T7u can be explained by an overshoot of the amplitude in the transition phase due to transient effects before the wave settles at a stable equilibrium amplitude. In principle, one could compute the coefficients of Eq. (12)*a priori*. This has been done for some fairly simple single-phase flows.^{5} For two-phase flows, the calculations rapidly become very complex and analytical progress is difficult. Therefore, we use Eq. (12) not for *a priori* prediction, but rather as an explanation and a theoretical description of the nonlinear saturation. In conclusion, Fig. 15 and the surrounding discussion establish that the waves created via the linear-stability mechanisms saturate at a finite amplitude. In this context, the saturated travelling waves can be regarded as the end-point of the dynamics and there is no ligament formation or droplet entrainment.

Scenario . | μ _{1}
. | μ _{2}
. | μ _{3}
. | RMSD . | R^{2}
. |
---|---|---|---|---|---|

dT4d | 0.7344 | −867.1 | 252 355 | 1.1306 ⋅ 10^{−3} | 0.9932 |

dT5s | 0.9290 | −913.9 | 0 | 8.1751 ⋅ 10^{−4} | 0.9933 |

dT6u | 1.1582 | −1410.1 | 0 | 5.7185 ⋅ 10^{−4} | 0.9957 |

dT7u | 2.7613 | −5466.6 | 0 | 9.0486 ⋅ 10^{−4} | 0.9749 |

Scenario . | μ _{1}
. | μ _{2}
. | μ _{3}
. | RMSD . | R^{2}
. |
---|---|---|---|---|---|

dT4d | 0.7344 | −867.1 | 252 355 | 1.1306 ⋅ 10^{−3} | 0.9932 |

dT5s | 0.9290 | −913.9 | 0 | 8.1751 ⋅ 10^{−4} | 0.9933 |

dT6u | 1.1582 | −1410.1 | 0 | 5.7185 ⋅ 10^{−4} | 0.9957 |

dT7u | 2.7613 | −5466.6 | 0 | 9.0486 ⋅ 10^{−4} | 0.9749 |

### B. High density contrast

To gain insight in the nonlinear wave dynamics under high-density-contrast conditions, we perform direct numerical simulations with the system parameters corresponding to Fig. 4(e), scenario T3d (Table I). This system exhibits two distinct linearly unstable modes, one long-wave shear mode and a short-wave interfacial mode, which are accounted for by *α*_{REYG} = 1.549 and *α _{TAN}* = 40.27 in Eq. (11), respectively. It is due to the short-wave nature of the interfacial mode combined with the high inertia of the liquid phase that the flow characteristics near the interface (Fig. 16) differ significantly from those observed in the low-density-contrast case. Unlike the scenarios presented in Sec. V A, no prominent features, such as recirculation zones, emerge within the wave train. However, a vortex layer with the familiar “cat’s eye structure” develops at the demarcation between the bulk of the gas and a thin gas layer dragged downwards by the liquid due to interfacial shear stress. Thereby, the vortices are pinched between two consecutive high-pressure regions, which are forming on the “upwind” side of the short-wave crests. Further snapshots indicate that this vortex layer is unstable to secondary instability. Hence, the wave form in Fig. 16 should not be regarded as a quasi-steady state.

Similar to the low-density-ratio case, the nonlinear dynamics of the interfacial mode appear to be consistent with Stuart-Landau theory, albeit these dynamics develop faster due to the higher growth rate of the first harmonic (*α _{TAN}* = 40.27). Evidence is provided in Fig. 17. Consistency with the Stuart-Landau theory is especially clear in the time-frequency domain. In particular in Fig. 17(b), where a power-spectral density

is shown (*T* corresponds to the duration of the simulation). A well-defined global maximum is observed at *ω _{r}* = 239.8, corresponding to the frequency of the most-unstable mode at

*α*= 40.27 in the spatial domain. Successive maxima at

*ω*=

_{r}*n*(239.8) with

*n*= 2, 3, … indicate that the harmonics of the most-unstable mode are slaved to the fundamental.

The maxima in the power-spectral density function, although well defined, are by no means sharp. The broadening of the maxima here is a sign that the wave train is not strictly periodic, and that a very large number of frequencies is present in the dynamics. Crucially, the broadness of the lines is a function of the simulation time: for longer simulations, more and more frequencies come into play (albeit that the same maxima remain dominant throughout). This indicates the onset of chaos. Thus, the Stuart-Landau theory manifests itself as the leading-order approximation to a chaotic dynamics, albeit that the wave train in Fig. 16 is subject to secondary instability.

Other second-order effects are key to understanding the secondary instability. The vortex layer in Fig. 16(a) is not steady but breaks down at later times (Fig. 16(b)). In particular, a long-wave perturbation to the vortex layer (wavelength on the domain scale) is seen to coincide with the significant growth of the long-wave shear mode. Therefore, it would appear that a complicated secondary instability sets in, involving a destabilization of the vortex layer by perturbations that are fed by the long-wavelength linearly unstable subdominant mode. (This mode is subdominant in the sense that it is linearly unstable but its growth rate is not as large as that of the most-unstable mode.) At this stage, individual waves steepen even further, up to the point where wave overturning is imminent (Fig. 16(b)). Summarizing, it is clear that secondary instability may inhibit the operation of the system in a quasi-steady laminar state at high density ratios.

## VI. CONCLUSIONS

We have presented a comprehensive study on two-dimensional laminar flow of a vertical film sheared by laminar counter-current gas in a confined channel. This study tries to further elucidate the nature of interfacial instability in such two-phase flows using several complementary methods, namely Orr-Sommerfeld analysis, energy budget analysis as well as high resolution direct numerical simulations. Two sample systems have been selected for investigation: one with a high density contrast (*ρ _{L}*/

*ρ*= 1000) and a second with a low density contrast (

_{G}*ρ*/

_{L}*ρ*= 10). In both cases, the same viscosity contrast (

_{G}*μ*/

_{L}*μ*= 50) and comparatively low surface tension (

_{G}*γ*= 1 ⋅ 10

^{−3}N m

^{−1}) were used. In our study, we focussed on analysing the influence of liquid film thickness and applied pressure drop on the development of interfacial waves.

Temporal linear stability analysis reveals that the liquid interface is inherently unstable for both cases. In the system with high density ratio, short-wave instability is predominant, whereas the low-density-contrast case tends to favour long-wave instability. Furthermore, the instability is governed by a multitude of coexisting linearly unstable modes (interfacial mode, shear mode in both phases, internal mode) under high-density-ratio conditions, where the latter two modes indicate the onset of turbulence in the bulk of one of the phases. In contrast, the low-density-ratio system exhibits only one linearly unstable mode, which is consistent with the Yih mode^{1} and a manifestation of the viscosity contrast.

Additionally, we use the phase velocity of the linearly most unstable mode to identify two distinct flow regimes, which we characterize by the direction of propagation of travelling waves. A standing wave marks the boundary of these regimes and is associated with the loading point. A flow map has been constructed to illustrate these regimes in a precise manner; for the high density contrast, the map reveals an island of upward-travelling disturbances amidst a sea of downward-travelling waves, which is the result of mode competition (shear mode in the gas layer vs. interfacial mode).

We have further determined the nature of the instability in the spatio-temporal framework under low-density-ratio conditions. Using this information, we established a second flow map indicating the transition from convective to absolute instability (C/A). Besides standard Orr-Sommerfeld analysis, we also adopted the analytic connection between temporal and spatio-temporal growth rates in the linear regime as presented by Ó Náraigh and Spelt.^{38} This approach, which is based on analytical continuation, circumvents possible difficulties in identifying the absolute growth rate that are associated with the multivalued nature of the eigenvalue problem and specifics of the problem at hand. Compared with OS analysis and DNS, this method shows good agreement and is therefore appropriate to accurately estimate the absolute growth rate of the instability. We find that the system is absolutely unstable in most parts of the parameter space considered herein with the exception of two narrow bands at low applied pressure drop and low film thickness, respectively. In the high-density-ratio case, mode competition and mode coalescence between the multiple linearly unstable modes hindered the mapping of convective and absolute instability in the respective parameter space. For such systems, linearized DNS^{39} or series solutions of the underlying stability problem^{40–42} may allow for more conclusive results.

To assess the development of interfacial waves up to finite amplitudes, we perform direct numerical simulations for parameters within the established flow regimes of the low-density-ratio case, using a level set method based solver that has been developed in-house. These simulations show excellent agreement with linear theory during the stage of exponential wave growth and confirm the determined flow patterns. The simulations also show saturation of the waves once nonlinear effects become important. A Fourier analysis reveals that the growth of the higher harmonics of the interfacial waves is coupled to that of the fundamental in a fashion which is consistent with weakly nonlinear theory. The growth of the first harmonic also agrees well with the Stuart-Landau model, thus underpinning the weakly nonlinear nature of the instability and, moreover, suggesting the existence of a supercritical bifurcation. Regarding high-density-contrast conditions, the dynamics of the (short-wave) interfacial mode appear to be similar to those observed under low density ratio. However, direct numerical simulations suggest that the emergence of an additional (long-wave) shear mode triggers a secondary instability, which leads to a chaotic wave train showing signs of imminent wave overturning.

Although we have focussed on laminar-laminar cases only in this work, the ideas and results contained herein can be extended to turbulent gas streams. In this context, the most accurate and appropriate methodology is full-scale DNS, which is the target of future work. However, short of full-scale DNS, a quasi-laminar model may be assumed for the linear stability of the two-phase flow^{39,45} or, alternatively, a weighted-residual integral boundary-layer model,^{22,46,47} which also models nonlinear interfacial waves. These approaches both have their own advantages and shortcomings, and the aim of future work will be to confirm these reduced-dimensional modelling approaches with evidence from accurate DNS, towards which the present work is an initial contribution.

In summary, the combination of generic complementary (semi-)analytical and numerical methods presented herein yields a comprehensive and convincing characterization of interfacial instability in vertical counter-current gas-liquid flows. We are therefore confident that this rigorous methodology can be employed to further elucidate the dynamics of parallel shear flows in a wide range of technically relevant parameter regimes, such as flows with high viscosity contrast, high surface tension or non-steady gas flow. Beyond that, the outlined approach may be used as a starting point for the future study of heat and mass transfer phenomena in vertical gas-liquid flows.

## Acknowledgments

The authors gratefully acknowledge the financial support of Sulzer Chemtech Ltd as well as the Scottish Funding Council through the Scottish Energy Technology Partnership (Project No. ETP24). This work made use of the facilities of ARCHER, the latest UK National Supercomputing Service (http://www.archer.ac.uk). ARCHER provides a capability resource to allow researchers to run simulations and calculations that require large numbers of processing cores working in a tightly coupled, parallel fashion. The ARCHER Service is provided by the EPSRC, NERC, EPCC, Cray Inc., and the University of Edinburgh. Further access to the facilities of ARCHER was provided through the ARCHER Resource Allocation Panel, Project No. e174. The optimization of the in-house solver was funded under the ARCHER Embedded Computational Science and Engineering (eCSE) Service.

### APPENDIX: FULL FORMULATION OF THE LINEAR STABILITY ANALYSIS AND NUMERICAL SOLUTION

#### 1. Base flow velocity profile

Here, we give a detailed formulation of the governing equations underlying the linear stability analysis undertaken in Secs. III and IV as well as a description of the numerical method used to solve the corresponding generalized complex eigenvalue problem.

Under the assumption of a steady, spatially uniform, laminar, and incompressible flow in both phases, the Navier-Stokes equations governing the undisturbed base flow reduce to standard balances between pressure as well as viscous and gravitational forces. For the liquid film, this balance can be written as

where $ u \u0303 0 $ is the dimensional base flow velocity in the respective phase (tildes denote dimensional quantities). Equation (A1) is subject to no-slip at the liquid side channel wall, $ z \u0303 =\u2212 d L $, and continuity of tangential stress at the interface, $ z \u0303 =0$,

Thus, integration of Eq. (A1) yields

as the velocity profile for the liquid film. The interfacial velocity, which constitutes one of the boundary conditions of the gas layer, reduces to

The velocity profile for the laminar gas layer is derived analogous to the liquid layer. We therefore write the force balance as

which is subject to continuity of velocity and shear stress at the interface,

Applying these interfacial conditions in Eq. (A5) yields

as the gas side velocity profile. To determine the interfacial shear *τ _{int}*, we apply the no-slip condition at the gas side channel wall, $ u \u0303 0 z \u0303 = d G =0$,

In summary, the velocity profile of the undisturbed base flow in dimensional form reads

#### 2. Perturbation equations and energy budget

As mentioned in Sec. II A, we introduce a small disturbance that shifts the flat interface from *z* = 0 to *z* = *η*, where $ \eta \u226a1$. This (dimensionless) wave elevation gives rise to perturbations in the flow field of the following form:

where the subscript zero denotes base flow quantities and the *δ* quantities are infinitesimally small perturbations. We use these variables of the flow field and obtain the linearized Navier-Stokes equations in both phases $ j = L , G $,

where $ r L , r G = r , 1 $, and $ m L , m G = m , 1 $. The pressure is eliminated from Eq. (A11) by taking the curl of both sides, which yields the linearized equation for the component of vorticity *δω _{y}* out of the plane generated by the wall-normal and streamwise directions,

where

It is further convenient to use the streamfunction representation (*δu*, *δw*) = (∂Ψ/∂*z*, − ∂Ψ/∂*x*) for the two-dimensional disturbance velocity field, hence

Assuming a wave-like solution of the form Ψ(*x*, *z*, *t*) = e^{i(αx−ωt)}*ψ*(*z*) for the streamfunction, Equations (A12) and (A14) lead to the Orr-Sommerfeld equations governing the stability of the liquid interface,

where *α* = *α _{r}* + i

*α*and

_{i}*ω*=

*ω*+ i

_{r}*ω*are the complex wavenumber and angular frequency, respectively. This equation is subject to no-penetration and no-slip conditions at both channel walls,

_{i} Furthermore, conditions for continuity of velocity and tangential stress as well as a jump in the normal stress due to surface tension are applied to match the streamfunction across the interface at *z* = 0 (we use the notation *c* = *ω*/*α*),

Using operator notation to rewrite Eqs. (A15)-(A17) yields

which highlights the generalized eigenvalue problem associated with the stability problem.

To understand the physical mechanism that causes the instability, we perform an energy-budget analysis. To find the energy budget of the system, we multiply Equation (A11) by *δ*** u** and integrate over a single wavelength in the streamwise

*x*-direction and over the entire wall-normal

*z*-domain. The right-hand side of Equation (A11) is first rewritten as the divergence of the stress tensor ∇ ⋅ 𝕋, where

In deriving the Orr-Sommerfeld equation, we find that the pressure terms have a representation in terms of the streamfunction $ \psi j z $,

Thus, we multiply Equation (A11) by *δ*** u** and integrate over the unit cell, $ 0 , \lambda \xd7 \u2212 \delta L , 0 $ for the liquid and $ 0 , \lambda \xd7 0 , \delta G $ for the gas layer. Here

*λ*= 2

*π*/

*α*is the wavelength of the disturbance (working with a sinusoidal disturbance gives periodic boundary conditions in the lateral direction). Using Gauss’ theorem on the stress term, we obtain the energy relation for each phase,

where the positive sign corresponds to the liquid phase and the negative sign to the gas phase (the integration limits are obvious and are not stated explicitly here). We sum up over *j* in Equation (A21) to obtain, in a standard fashion, the energy-budget relation

where

Lastly, the term “*INT*” is related to the following interfacial conditions:

which is decomposed into normal and tangential contributions,

where

and

Note that the tangential contribution can be further decomposed to highlight the effect of the viscosity contrast,

where $\eta = ( \delta u L \u2212 \delta u G ) / [ u 0 \u2032 ( 0 + ) \u2212 u 0 \u2032 ( 0 \u2212 ) ] $ is the height of the perturbed interface (*cf*. Equation (A17b)) and where we have used the continuity of tangential stress at the interface to write 𝕋_{zx,L} = 𝕋_{zx,G} ≔ 𝕋_{zx} at *z* = 0. Thus, provided the absolute value of the phase difference between *η* and the tangential stress does not exceed *π*/2, a viscosity contrast *m* > 1 implies that the *TAN* term is a source of instability.

#### 3. Numerical method

We solve the eigenvalue equation (A18) numerically by employing a standard Chebyshev collocation method^{28} with an approximation for the streamfunction of the form

where $ T n \u22c5 $ is the *n*^{th} Chebyshev polynomial of the first kind. We substitute this trial solution into Eq. (A18) and evaluate the result at *N _{L}* +

*N*− 6 Gauss-Lobatto collocation points. Together with the eight boundary and interfacial conditions, this yields

_{G}*N*+

_{L}*N*+ 2 linear equations in as many unknowns. In matrix terms, we solve

_{G} where $v= a 0 , \u2026 , a N L , b 0 , \u2026 , b N G T $. We use a standard linear algebra package (MATLAB^{®}) to solve this equation, thereby adjusting the number of collocation points $ N L + 1 , N G + 1 $ until convergence is reached.