3.1. Liquid-side WRIBL model

The liquid film (subscript $L$), with density $\rho _L$, dynamic viscosity $\mu _L$ and surface tension $\sigma$, is governed by the dimensionless continuity and Navier–Stokes equations written in Cartesian coordinates $x$ and $y$ (figure 4):

(3.1a)$$\begin{gather} \partial_xu_L + \partial_yv_L = 0, \end{gather}$$

(3.1b)$$\begin{gather}\epsilon(\partial_tu_L + u_L\,\partial_xu_L + v_L\,\partial_yu_L)={-} \epsilon\, \partial_xp + \frac{1}{Re_L}\,(\epsilon^2\,\partial_{xx}u_L + \partial_{yy}u_L) + \frac{\sin(\phi)}{Fr ^2}, \end{gather}$$

(3.1c)$$\begin{gather}\epsilon^3(\partial_tv_L + u_L\,\partial_xv + v_L\,\partial_yv_L) ={-} \epsilon\,\partial_yp_L + \frac{1}{Re_L}\,(\epsilon^4\,\partial_{xx}v_L + \epsilon^2\,\partial_{yy}v_L) - \epsilon\,\frac{\cos(\phi)}{Fr ^2}, \end{gather}$$

where $Re_L=\rho _L\mathcal {U}_L\mathcal {L}/\mu _L$ and $Fr=\mathcal {U}_L/\sqrt {\mathcal {L}g}$ denote the liquid Reynolds number and Froude number, and where we have applied the scaling

(3.1d)\begin{equation} u_L=\frac{u^\star_L}{\mathcal{U}_L},\quad v_L=\frac{v_L^\star}{\epsilon\mathcal{U}_L},\quad x=\epsilon\,\frac{x^\star}{\mathcal{L}},\quad y=\frac{y^\star}{\mathcal{L}},\quad t=\epsilon t^\star\,\frac{\mathcal{U}_L}{\mathcal{L}},\quad p_L=\frac{p_L^\star}{\rho_L\mathcal{U}_L^2}. \end{equation}

Here, we have introduced the long-wave parameter $\epsilon =\mathcal {L}/\varLambda ^\star$, which relates the cross-stream length scale $\mathcal {L}$ to the streamwise length scale given by the wavelength $\varLambda ^\star$. For the purpose of the current derivation, it suffices to say that the scales $\mathcal {L}$ and $\mathcal {U}_L$ are representative of the film thickness $h^\star$ and streamwise liquid velocity $u_L^\star$. In § 3.3, we will rescale our problem and make the final choice for $\mathcal {L}$ and $\mathcal {U}_L$.

The system is closed with the boundary conditions at $y=0$

(3.1e)\begin{equation} u_L=v_L=0, \end{equation}

the kinematic condition

(3.1f)\begin{equation} \left.{v_L}\right|_{y=h}=\left.{u_L}\right|_{y=h} \partial_{x} h + \partial_{t} h, \end{equation}

and the inter-phase stress coupling conditions at $y=h$

(3.1g)\begin{align}& -\partial_yu_L\,\frac{1}{1+\epsilon^2\,\partial_xh^2}\, (-\epsilon^4\,\partial_xh^2\,\partial_xv_L -4 \epsilon^2\, \partial_xh\,\partial_xu_L + \epsilon^2\,\partial_xv_L)= \frac{\varPi_\mu \varPi_u}{\varPi_L}\,T_G, \end{align}

(3.1h)\begin{align} & \epsilon P_L + \frac{2}{1+\epsilon^2\,\partial_xh^2}\,\frac{1}{Re_L}\, (\epsilon^4\,\partial_xh^2\,\partial_xu_L - \epsilon^4\,\partial_xh\, \partial_xv_L -\epsilon^2\,\partial_xu_L - \epsilon^2\,\partial_xh\,\partial_yu_L) \nonumber\\ &\quad -\epsilon^3\,We\,\partial_{xx}h = \frac{1}{Re_G}\,\frac{\varPi_\rho \varPi_u^2}{\varPi_L}\, P_G, \end{align}

where $We=\sigma/(\rho_L\mathcal{L}\mathcal{U}_L^2)$ denotes the Weber number. The liquid–gas coupling enters through $T_G$ and $P_G$, which are scaled as follows:

(3.2a,b)\begin{equation} T_G=\frac{\mathcal{L}_G}{\mu_G \mathcal{U}_G}\,T_G^\star,\quad P_G=\underline{\epsilon}\,\frac{\mathcal{L}_G}{\mu_G \mathcal{U}_G}\,P_G^\star, \end{equation}

where $\mathcal {L}_G$, $\mathcal {U}_G$ and $\underline {\epsilon }=\mathcal {L}_G/\varLambda ^\star =\epsilon \varPi _L$ denote the gas-side cross-stream length scale, velocity scale and long-wave parameter, which will be defined in § 3.2. As a result, the gas Reynolds number $Re_G=\rho _G\mathcal {U}_G\mathcal {L}_G/\mu _G$, the velocity scale ratio $\varPi _u=\mathcal {U}_G/\mathcal {U}_L$, the length scale ratio $\varPi _L=\mathcal {L}_G/\mathcal {L}$, and the viscosity and density ratios $\varPi _\mu =\mu _G/\mu _L$ and $\varPi _\rho =\rho _G/\rho _L$ enter (3.1g) and (3.1h).

Next, we apply the WRIBL approach to derive two evolution equations involving the local instantaneous liquid flow rate $q(x,t)$ and the film thickness $h(x,t)$. In principle, we follow the same steps as Samanta (Reference Samanta2014), except that we account for the gas pressure $P_G$, which plays an important role in our current configuration, allow $P_G$ and $T_G$ to vary in space and time, and account for turbulence in the gas.

First, the governing equations (4.1a,b) are truncated at ${O}(\epsilon ^2)$, except for inertial terms, which are truncated at ${O}(Re_L\,\epsilon )$. Next, we eliminate $p$ from (3.1b) via an integrated form of (3.1c) using (3.1h). Then we substitute for the streamwise velocity $u$ ($v$ is obtained from (3.1a)) the decomposition

(3.3)\begin{equation} u_L=\hat{u}_L+\epsilon u_L^{(1)}, \end{equation}

where the base profile $\hat {u}_L$ is governed by

(3.4)\begin{equation} \left.\begin{gathered} \partial_{yy}\hat{u}_L=\text{const.},\quad \left.{\partial_y\hat{u}_L}\right|_{y=h}= \frac{\varPi_\mu\varPi_u}{\varPi_L}\,T_G,\quad \left.{\hat{u}_L}\right|_{y=0}=0,\\ \int_0^{h(x,t)}\hat{u}_L\,\mathrm{d}y=q_L(x,t). \end{gathered}\right\} \end{equation}

Finally, the unknown ${O}(\epsilon )$ velocity correction $\epsilon u^{(1)}_L$ is eliminated from the problem by multiplying the truncated form of (3.1b) with a weight function $w(y)$, integrating the result across the film thickness $h(x,t)$, and applying the tangential inter-phase coupling condition (3.1g). The weight function $w$ satisfies

(3.5a–c)\begin{equation} \partial_{yy}w=\text{const.},\quad \left.{w}\right|_{y=0}=0,\quad \partial_yw|_{y=h}=0. \end{equation}

As a final result, we obtain the integral momentum equation for the liquid film,

(3.6a) \begin{align} & \partial_t q_L+\frac{17}{7}\,\frac{q_L}{h}\,\partial_xq_L-\frac{9}{7}\, \frac{q^2_L}{h^2}\,\partial_xh=\frac{5}{6}\,We\,h\,\partial_{xxx}h+ \frac{5}{6}\,Fr^{{-}2}\,h\left\{\sin\left(\phi\right)-\cos(\phi)\,\partial_xh\right\} \nonumber\\ &\quad +Re_L^{{-}1}\left\{-\frac{5}{2}\,\frac{q_L}{h^2}+4\,\frac{q_L}{h^2}\, \partial_xh^2-\frac{9}{2h}\,\partial_xq_L\,\partial_xh-6\,\frac{q_L}{h}\,\partial_{xx}h+ \frac{9}{2}\,\partial_{xx}q_L\right\} \nonumber\\ &\quad +\frac{\varPi_\mu\varPi_u}{\varPi_L}\,T_G\left\{Re_L^{{-}1} \left[\frac{5}{4}+\frac{h}{6}\,\partial_{xx}h+\frac{1}{2}\,\partial_{x}h^2\right]- \frac{5}{112}\,q_L\,\partial_xh-\frac{19}{336}\,\partial_xq_L h\right\} \nonumber\\ &\quad -\frac{19}{672}\,\frac{\varPi_\mu^2\varPi_u^2}{\varPi_L^2}\,h^2\,\partial_xh T_G^2- \frac{5}{6}\,Re_G^{{-}1}\,\frac{\varPi_\rho\varPi_u^2}{\varPi_L}\,h\,\partial_xP_G \nonumber\\ &\quad +\frac{\varPi_\mu\varPi_u}{\varPi_L}\left\{\partial_xT_G\left[Re_L^{{-}1}\, \frac{3}{4}\,h\,\partial_{x}h-\frac{15}{224}\,hq_L\right]-\frac{25}{1344}\, \frac{\varPi_\mu\varPi_u}{\varPi_L}\,h^3T_G\,\partial_xT_G-\frac{h^2}{48}\, \partial_tT_G\right\}, \end{align}

to which is added an integral continuity equation obtained by integrating (3.1a) across the liquid film and applying (3.1f):

(3.6b)\begin{equation} \partial_th+\partial_xq_L=0. \end{equation}

In the limit $T_G=\partial _xP_G=0$, (3.6a) reduces to (41) from Ruyer-Quil & Manneville (Reference Ruyer-Quil and Manneville2000). In the limit $\partial _xT_G=\partial _tT_G=\partial _xP_G=0$, it collapses with (3.9) from Samanta (Reference Samanta2014), except for a typo in the $T_G\,\partial _xh^2$ term of that reference. Here, we will neglect the terms involving $\partial _xT_G$ and $\partial _tT_G$, but we will account for the space and time variation of $T_G(x,t)$ and $P_G(x,t)$ in the remaining terms. This amounts to a quasi-developed approach. See Appendix C for a justification of this approximation.

Versus the model of Tseluiko & Kalliadasis (Reference Tseluiko and Kalliadasis2011), which is based on a linear representation of the gas response, our model accounts for the gas pressure $P_G$, which plays a role for the confinement considered here. It also accounts for streamwise viscous diffusion in the liquid, which is known to affect the dynamics of precursory capillary ripples (Ruyer-Quil & Manneville Reference Ruyer-Quil and Manneville2002).

3.2. Gas-side asymptotic long-wave model

We represent the turbulent flow of the gas (subscript $G$), with density $\rho _G$ and dynamic viscosity $\mu _G$, in two dimensions via the (dimensionless) Reynolds-averaged continuity and steady Navier–Stokes (RANS) equations, written here in the Cartesian gas-side coordinates $x$ and $\underline {y}$ (see figure 4) as

(3.7)\begin{equation} \partial_x u_G+\partial_{\underline{y}}v_G=0, \end{equation}

(3.8a)\begin{align} {\underline{\epsilon}}( u_G\,\partial_x u_G + v_G\,\partial_{\underline{y}} u_G) &={-}\frac{1}{{Re_G}}\,\partial_x p_G + \frac{\varPi_L}{\varPi_u^2}\,\frac{\sin(\phi)}{Fr^2} + \frac{1}{{Re_G}} \left\{\partial_{\underline{y} \underline{y}} u_G+{\underline{\epsilon}}^2\,\partial_{xx} u_G\right\} \nonumber\\ &\quad + \frac{1}{{Re_G}} \left\{\partial_{\underline{y}} \left(\frac{\mu_{t}}{\mu_G}\,\partial_{{\underline{y}}} u_G \right)+ {\underline{\epsilon}}^2\,\partial_x \left(\frac{\mu_{t}}{\mu_G}\, \partial_{x} u_G \right)\right\}, \end{align}

(3.8b)\begin{align} {\underline{\epsilon}}^3( u_G\,\partial_x v_G + v_G\, \partial_{{\underline{y}}} v_G) &={-}\frac{1}{{Re_G}}\, \partial_{{\underline{y}}} p_G+{\underline{\epsilon}} \frac{\varPi_L}{\varPi_u^2}\,\frac{\cos(\phi)}{Fr ^2} + \frac{1}{{Re_G}} \left\{{\underline{\epsilon}}^2\,\partial_{{\underline{y}} {\underline{y}}} v_G+{\underline{\epsilon}}^4\,\partial_{xx} v_G\right\} \nonumber\\ &\quad + \frac{1}{{Re_G}} \left\{{\underline{\epsilon}}^2\,\partial_{\underline{y}} \left(\frac{\mu_{t}}{\mu_G}\,\partial_{{\underline{y}}} v_G \right)+{\underline{\epsilon}}^4\,\partial_x \left(\frac{\mu_{t}}{\mu_G}\,\partial_{x} v_G \right)\right\}, \end{align}

where $\mu _{t}$ denotes the turbulent viscosity, $Re_G=\rho _G\mathcal {U}_G\mathcal {L}_G/\mu _G$ is the gas Reynolds number, and we have applied the scaling

(3.9a–e)\begin{equation} u_G=\frac{u^\star_G}{\mathcal{U}_G},\quad v_G=\frac{v^\star_G}{\underline{\epsilon}\mathcal{U}_G},\quad x=\underline{\epsilon}\,\frac{x^\star}{\mathcal{L}_G},\quad \underline{y}=\frac{\underline{y}^\star}{\mathcal{L}_G},\quad p_G=p^\star_G\,\frac{\underline{\epsilon}\mathcal{L}_G}{\mu_G\mathcal{U}_G}, \end{equation}

introducing the gas-side long-wave parameter $\underline {\epsilon }=\mathcal {L}_G/\varLambda ^\star$. For the gas-side reference scales, we choose once and for all

(3.10a,b)\begin{equation} \mathcal{L}_G=H^\star,\quad \mathcal{U}_G=\frac{q_{G 0}^\star}{H^\star}, \end{equation}

where $q_{G 0}^\star$ is the nominal gas flow rate per unit width of the primary flow (subscript $0$), thus $\mathcal {U}_G$ corresponds to the superficial gas velocity. We have scaled pressure with a measure for the viscous shear stress, in contrast to (3.1d), where the dynamic pressure was used.

The turbulent viscosity $\mu _t$ is formulated via the mixing-length approach (Prandtl Reference Prandtl1925):

(3.11)\begin{equation} \frac{\mu_t}{\mu_G}=Re_G\,l_t^2\left|\partial_{\underline{y}}u_G\right|, \end{equation}

where $l_t=l_t^\star /\mathcal {L}_G$ denotes the dimensionless mixing length. At this point, a remark about choosing a turbulent viscosity model, such as (3.11), is in order. Luchini & Charru (Reference Luchini and Charru2019) have shown that such models cannot fully reproduce the momentum redistribution induced by wall perturbations to a parallel turbulent flow. Nonetheless, comparisons with different experiments (Zilker, Cook & Hanratty Reference Zilker, Cook and Hanratty1977; Frederick & Hanratty Reference Frederick and Hanratty1988) have shown that turbulent viscosity models based on the van Driest equation, which will be introduced in (3.22), capture satisfactorily the linear (Russo & Luchini Reference Russo and Luchini2016) and nonlinear (Tseluiko & Kalliadasis Reference Tseluiko and Kalliadasis2011; Camassa et al. Reference Camassa, Ogrosky and Olander2017) responses of the wall shear stress. Thus such models allow us to account adequately for the inter-phase coupling in our current configuration.

We assume a large gas/liquid velocity contrast $\varPi _u\gg 1$, which warrants two simplifications. First, we have neglected time derivatives in (3.8), as

(3.12)\begin{equation} {O}\left\{\frac{\partial_{t^\star} u^\star_G}{u^\star_G\, \partial_{x^\star} u^\star_G}\right\}=\frac{1}{\varPi_u}\ll 1, \end{equation}

assuming that the time scale is dictated by the waviness of the liquid film, i.e. $\mathcal {T}=\varLambda ^\star /\mathcal {U}_L$. Second, we set zero-velocity conditions at the film surface $\underline {y}=d$:

(3.13a)\begin{equation} u_G=v_G=0. \end{equation}

Thus, from the point of view of the gas, the film surface is represented as a frozen wavy wall (Tseluiko & Kalliadasis Reference Tseluiko and Kalliadasis2011). Our system is closed via a symmetry condition at $\underline {y}=0$:

(3.13b)\begin{equation} \partial_{\underline{y}}u_G=v_G=0. \end{equation}

The ultimate aim of the gas-side model, to be derived next, is to obtain the inter-phase coupling quantities in (3.6a), which are evaluated at $\underline {y}=d$, implying $l_{t}=0$:

(3.14a)\begin{gather} P_G=p_G, \end{gather}

(3.14b)\begin{align} T_G=\frac{T_G^\star}{\mu_G\mathcal{U}_G/\mathcal{L}_G} &={-}\partial_{\underline{y}} u_G-\left\{\epsilon^2\,\partial_xd^2-1\right\}^{{-}1} \left\{2\epsilon^2\,\partial_xd\,\partial_{\underline{y}} v_G\right. \nonumber\\ &\quad \left.{}+2\epsilon^2\,\partial_xd\,\partial_xu_G+\epsilon^2\,\partial_xv_G-\epsilon^4\, \partial_xd^2\,\partial_xv_G\right\}. \end{align}

Following Camassa et al. (Reference Camassa, Ogrosky and Olander2017), we introduce the curvilinear coordinates $\eta$ and $\xi$ (see figure 4), which will facilitate the account of turbulence:

(3.15a,b)\begin{equation} \eta=\underline{y}\,\frac{\bar{d}}{d},\quad \xi=x+\underline{\epsilon}\,F(\xi,\eta), \end{equation}

where $\bar {d}$ denotes the spatial average of $d$, and orthogonality implies

(3.16)\begin{equation} \partial_\eta F=\underline{\epsilon}\,\frac{d}{\bar{d}^2}\,\partial_\xi d\, \frac{1}{\underline{\epsilon}\,\partial_\xi F-1}. \end{equation}

Red dashed lines in figure 4 represent curves of constant $\eta$ and $\xi$, where

(3.17)\begin{equation} \left.{\partial_x\underline{y}}\right|_\eta={-}\left.{\partial_{\underline{y}}x}\right|_\xi= \frac{\underline{y}}{\bar{d}}\,\partial_xd. \end{equation}

Next, we recast the governing equations (3.8) and (3.13) in the curvilinear coordinate system (tilde symbol), using the projection rules

(3.18a,b)$$\begin{gather} u_G=\tilde{u}_G+{O}(\underline{\epsilon}^2),\quad v_G=\tilde{v}_G+\frac{\partial_\xi d}{\bar{d}}\,\eta \tilde{u}_G+{O}(\underline{\epsilon}), \end{gather}$$

(3.19a,b)$$\begin{gather}\partial_x=\partial_{\xi}-\eta\,\frac{\partial_\xi d}{d}\,\partial_\eta+ {O}(\underline{\epsilon}), \quad \partial_{\underline{y}}= \frac{\bar{d}}{d}\,\partial_\eta+{O}(\underline{\epsilon}^2), \end{gather}$$

and truncate the result at ${O}(\underline {\epsilon }^1)$. Upon eliminating the pressure variable $p$ in (3.8a) via an appropriate integration of (3.8b), we obtain

(3.20a)\begin{align} &\frac{\partial_\xi d}{d}\,\partial_\xi \tilde{u}_G+\partial_\xi \tilde{u}_G+ \frac{\bar{d}}{d}\,\partial_\eta \tilde{v}_G= 0, \end{align}

(3.20b)\begin{align} \underline{\epsilon}\tilde{u}_G\,\partial_\xi \tilde{u}_G+\underline{\epsilon}\, \frac{\bar{d}}{d}\,\tilde{v}_G\,\partial_\eta \tilde{u}_G &={-}\frac{1}{Re_G}\,\partial_\xi P_G+\frac{\varPi_L}{\varPi_u^2}\, \frac{1}{Fr^2}\left(\sin(\phi)+\underline{\epsilon}\cos(\phi)\,\partial_\xi d\right) \nonumber\\ &\quad +\frac{1}{Re_G}\,\frac{\bar{d}^2}{d^2}\left\{\partial_{\eta\eta} \tilde{u}_G+ \partial_\eta\left[\frac{\tilde{\mu}_t}{\mu_G}\,\partial_\eta\tilde{u}_G\right]\right\}, \end{align}

where $P_G=p_G|_{\eta =\bar {d}}$, and $\tilde {\mu }_t$ satisfies

(3.21)\begin{equation} \frac{\tilde{\mu}_t}{\mu_G}=\frac{d}{\bar{d}}\,Re_G\,\tilde{l}_t^2 \left|\partial_\eta \tilde{u}_G\right|, \end{equation}

with $\tilde {l}_t=l_t\bar {d}$/$d$.

In this curvilinear formulation, the variation of the mixing length $\tilde {l}_t$ is expressed in terms of $\eta$, i.e. normal to the film surface, thus correlations for parallel flows can be used. Following Tseluiko & Kalliadasis (Reference Tseluiko and Kalliadasis2011), we employ the van Driest equation (Van Driest Reference Van Driest1956)

(3.22)\begin{equation} \tilde{l}_t=\kappa\left(\bar{d}-\eta\right)\left\{1-\exp\left[\sqrt{\left|T_{G 0}\right|Re_G}\, \frac{\eta-\bar{d}}{A}\right]\right\}, \end{equation}

where $A=26$, $\kappa =0.41$ is the von Kármán constant, and $T_{G 0}$ denotes the primary-flow tangential stress, obtained by evaluating (3.30) in the limit $\underline {\epsilon }=0$, which intervenes in the traditional scaling based on the friction velocity $\mathcal {U}^+$:

(3.23a,b)\begin{equation} \mathcal{U}^+{=}\left\{\rho_G^{{-}1}\left|T^\star_{G 0}\right|\right\}^{1/2},\quad \mathcal{L}^+{=}\frac{\mu_G}{\rho_G \mathcal{U}^+}. \end{equation}

Finally, the boundary conditions (3.13) become

(3.24a,b)\begin{equation} \left.{\tilde{u}_G}\right|_{\eta=\bar{d}}=\left.{\tilde{v}_G}\right|_{\eta=\bar{d}}=0,\quad {\partial_\eta \tilde{u}_G}|_{\eta=0}=\left.{\tilde{v}_G}\right|_{\eta=0}=0. \end{equation}

The boundary-value problem (BVP) given by (3.20) and (3.24a,b) is solved order by order based on a regular expansion in $\underline {\epsilon }$ around $\underline {\epsilon }=0$ (Camassa et al. Reference Camassa, Ogrosky and Olander2017):

(3.25a)\begin{align}\tilde{u}_G=\tilde{u}_G^{(0)}+\underline{\epsilon}\tilde{u}_G^{(1)}+{O}(\underline{\epsilon}^2),\end{align}

(3.25b)\begin{align}\tilde{v}_G=\tilde{v}_G^{(0)}+ \underline{\epsilon}\tilde{v}_G^{(1)}+{O}(\underline{\epsilon}^2),\end{align}

(3.25c)\begin{align}P_G=P_G^{(0)}+\underline{\epsilon}P_G^{(1)}+{O}(\underline{\epsilon}^2).\end{align}

The zeroth-order problem is obtained by inserting (3.25a) into (3.20) and (3.24a,b), and then truncating at ${O}(\underline {\epsilon }^0)$. We anticipate a solution in the form of the product ansatz

(3.26)\begin{equation} \tilde{u}_G^{(0)}=g_0(\xi)\,U_0(\eta)=\frac{\bar{d}}{d}\,U_0(\eta), \end{equation}

which leads to the variable-separated zeroth-order momentum equation

(3.27a)\begin{equation} \frac{d^3}{\bar{d}^3}\left\{\frac{1}{Re_G}\,\partial_\xi P_G^{(0)}- \frac{\varPi_L}{\varPi_u^2}\,\frac{1}{Fr^2}\sin(\phi)\right\}=\frac{1}{Re_G}\, \partial_{\eta\eta}U_0+\partial_\eta\left\{\tilde{l}_t^2\,\text{sgn} \left(\partial_\eta U_0\right)\left(\partial_\eta U_0\right)^2\right\}=C_0, \end{equation}

subject to the boundary conditions

(3.27b)\begin{equation} \left.{U_0}\right|_{\eta=\bar{d}}={\partial_\eta U_0}|_{\eta=0}=0, \end{equation}

where we have employed the signum function $\textrm {sgn}$ to substitute $|\partial _\eta U_0|=\textrm {sgn}(\partial _\eta U_0)\,\partial _\eta U_0$, and the separation constant $C_0$ is obtained from the gauge condition

(3.27c)\begin{equation} \int_0^d\tilde{u}_G^{(0)}\,{\rm d}\tilde{y}=\int_0^{\bar{d}}U_0\,{\rm d}\eta=\frac{q_{G 0}}{2}=\frac{1}{2}. \end{equation}

At the next order, i.e. ${O}(\underline {\epsilon }^1)$, we obtain in a similar way

(3.28a)\begin{align} & \frac{d^2}{\bar{d}^2}\left\{\frac{1}{Re_G}\,\partial_\xi P_G^{(1)}\, \frac{d}{\partial_\xi d}-\frac{\varPi_L}{\varPi_u^2}\,\frac{1}{Fr^2}\cos(\phi)\,d\right\} = U_0^2+\frac{1}{Re_G}\,\partial_{\eta\eta}U_1 \nonumber\\ &\quad +\partial_\eta\left\{\tilde{l}_t^2\,\text{sgn}\left(\partial_\eta U_0\right) \partial_\eta U_0\,\partial_\eta U_1\right\}=C_1, \end{align}

(3.28b)$$\begin{align}& \left.{U_1}\right|_{\eta=\bar{d}}={\partial_\eta U_1}|_{\eta=0}=0, \end{align}$$

(3.28c)$$\begin{align}&\int_0^d\tilde{u}_G^{(1)}\,{\rm d}\tilde{y}=\int_0^{\bar{d}}U_1\,{\rm d}\eta=0, \end{align}$$

where we have employed the product ansatz

(3.29)\begin{equation} \tilde{u}_G^{(1)}=g_1(\xi)vU_1(\eta)=\frac{\partial_\xi d}{d}\,U_1(\eta), \end{equation}

and the separation constant $C_1$ is obtained from (3.28c).

The two BVPs (3.27) and (3.28) are solved numerically for $U_0$, $U_1$, $C_0$ and $C_1$ via the continuation software Auto07P (Doedel Reference Doedel2008). The solution is obtained for a given $\bar {d}$ on a fixed domain spanning $0\le \eta \le \bar {d}$. Based on this, the coupling quantities $T_G$ and $\partial _x P_G$, which appear in the liquid-side model (3.6a), are constructed readily at ${O}(\underline {\epsilon }^1)$:

(3.30a)\begin{align} T_G&={-}\frac{\bar{d}}{d}\left.\partial_\eta \tilde{u}_G\right|_{\eta=\bar{d}}+ {O}(\underline{\epsilon}^2)={-}\frac{\bar{d}}{d}\left\{\left.\partial_\eta \tilde{u}_G^{(0)}\right|_{\eta=\bar{d}}+\underline{\epsilon}\left.\partial_\eta \tilde{u}_G^{(1)}\right|_{\eta=\bar{d}} \right\}+{O}(\underline{\epsilon}^2), \nonumber\\ &={-}\frac{\bar{d}^2}{d^2}\left\{\left.\partial_\eta U_0\right|_{\eta=\bar{d}}+ \frac{\partial_{x^\star} d^\star}{\bar{d}}\left.\partial_\eta U_1\right|_{\eta=\bar{d}} \right\}+ {O}(\underline{\epsilon}^2), \end{align}

(3.30b)\begin{align} \partial_x P_G&=\partial_\xi P_G^{(0)}+\underline{\epsilon}\,\partial_\xi P_G^{(1)}+ {O}(\underline{\epsilon}^2) \nonumber\\ &=Re_G\left\{\frac{\bar{d}^3}{d^3} \left(C_0+C_1\,\frac{\partial_{x^\star}d^\star}{\bar{d}}\right)+\frac{\varPi_L}{\varPi_u^2} \frac{1}{Fr^2}\left(\sin(\phi)+\cos(\phi)\,\partial_{x^\star}d^\star\right)\right\}+ {O}(\underline{\epsilon}^2), \end{align}

where we have used the velocity expansion (3.25a):

(3.31)\begin{equation} \tilde{u}_G=\frac{\bar{d}}{d}\,U_0+\frac{\partial_{x^\star}d^\star}{d}\,U_1+ {O}(\underline{\epsilon}^2). \end{equation}

Importantly, at fixed $\bar {d}$, $T_G$ and $\partial _x P_G$ in (3.30) depend only on $d=D-h/\varPi _L$, which varies with $x$ and $t$. By contrast, Samanta (Reference Samanta2014) assumed $T_G=\textrm {const}$ and $\partial _x P_G=0$.

In contrast to the gas-side description of Demekhin (Reference Demekhin1981) and Tseluiko & Kalliadasis (Reference Tseluiko and Kalliadasis2011), (3.30) is obtained from a long-wave not small-wave amplitude expansion. Thus it works better when the liquid holdup is larger, whereas the cited models work better when the liquid holdup is small, i.e. $h^\star /\bar {d}^\star \to 0$.

As a result of our frozen-interface assumption ($\varPi _u{\gg }1$) expressed via (3.13), one would obtain exactly the same relations for the functions $U_0$ and $U_1$ appearing in (3.30) should one apply no-slip and no-penetration conditions at $y^\star =H^\star$ instead of a symmetry condition at $y^\star =D^\star$. This is because the BVPs for $U_0$ (3.27) and $U_1$ (3.28) would remain symmetrical in that case. Thus, up to the order of expansion of our WRIBL-LW model, our symmetry condition (3.13b) is valid without loss of generality.

3.3. Rescaling

For the remainder of the paper, we rescale streamwise lengths by setting $\epsilon =\underline {\epsilon }=1$, and we choose

(3.32a–c)\begin{equation} \mathcal{L}=\mathcal{L}_G=H^\star,\quad \mathcal{U}_L=\frac{q_{L 0}^\star}{H^\star},\quad \mathcal{U}_G=\frac{q_{G 0}^\star}{H^\star}. \end{equation}

This implies $\varPi _L=1$, i.e. all lengths are now scaled with the channel height $H^\star$. We recall that $q_{L 0}^\star$ and $q_{G 0}^\star$ are the primary-flow liquid and gas flow rates per unit width, thus $\mathcal {U}_L$ and $\mathcal {U}_G$ are the superficial velocities. The corresponding Reynolds numbers are

(3.33a,b)\begin{equation} Re_L=\frac{q_{L 0}^\star}{\nu_L},\quad Re_G=\frac{q_{G 0}^\star}{\nu_G}, \end{equation}

where $\nu _L=\mu _L/\rho _L$ and $\nu _G=\mu _G/\rho _G$.

At some places, we will rescale quantities with the natural scales

(3.34a–c)\begin{equation} \mathcal{L}_\nu=\nu_L^{2/3}g^{{-}1/3},\quad \mathcal{U}_\nu=(\nu_L g)^{1/3},\quad \mathcal{T}_\nu=\frac{\mathcal{L}_\nu}{\mathcal{U}_\nu}=\nu_L^{1/3}g^{{-}2/3}. \end{equation}

3.4. Model computations

We perform three types of numerical computations based on our WRIBL-LW model (3.6), (3.30): linear stability calculations, nonlinear computations of TWS, and nonlinear computations of spatially evolving falling liquid films.

To obtain the linear stability formulation, we perturb the dependent variables $q_L$ and $h$ around their primary flow values $q_{L 0}$ and $h_0$:

(3.35a)$$\begin{gather} q_L=q_{L 0}+\check{q}_L(x,t)=q_{0}+\hat{q}\exp\{{\rm i}(kx-\omega t)\}, \end{gather}$$

(3.35b)$$\begin{gather}h=h_0+\check{h}(x,t)=h_0+\hat{h}\exp\{{\rm i}(kx-\omega t)\}, \end{gather}$$

where the check mark denotes infinitesimal perturbations, $\omega$ denotes the angular frequency, and $\hat {q}_L=\hat {h}\omega /k$ follows from (3.6b). Surface waves resulting from the Kapitza instability grow spatially, but a counter-current gas flow can cause the onset of AI. Both phenomena can be captured via a spatial stability formulation (Vellingiri et al. Reference Vellingiri, Tseluiko and Kalliadasis2015). Thus we will usually (but not exclusively) assume $k \in \mathbb {C}$ and $\omega \in \mathbb {R}$, with

(3.36)\begin{equation} k = k_r + {\rm i} k_i, \end{equation}

where $k_r=2{\rm \pi} /\varLambda$ is the physical wavenumber, and $-k_i$ is the spatial growth rate.

The film surface perturbation (3.35b) translates to the gas-side problem via

(3.37)\begin{equation} d=d_0+\check{d}=d_0+\hat{d}\exp\{{\rm i}(kx-\omega t)\},\quad\hat{d}={-}\frac{\hat{h}}{\varPi_L}. \end{equation}

Inserting this in (3.30) and then linearizing yields the linear responses of the inter-phase coupling quantities:

(3.38a)$$\begin{gather} T_G=T_{G 0}+\check{T}_G=T_{G 0}+\hat{T}_G\exp\{{\rm i}(kx-\omega t)\}, \end{gather}$$

(3.38b)$$\begin{gather}P_G=P_{G 0} + \check{P}_G=P_{G 0}+\hat{P}_G\exp\{{\rm i}(kx-\omega t)\}, \end{gather}$$

with

(3.39a,b)$$\begin{gather} T_{G 0} ={-}\partial_{\eta} U_0 |_{d_0},\quad \hat{T}_G=\frac{\hat{d}}{d_0}\left\{2\,\partial_{\eta} U_0 |_{d_0}+ \varPi_L {\rm i}k\,\partial_{\eta} U_1 |_{d_0}\right\}, \end{gather}$$

(3.40a,b)$$\begin{gather}\partial_xP_{G 0}=Re_{G}\left\{C_0+\frac{\varPi_L}{\varPi_u^2}\, \frac{\sin(\phi)}{Fr^2}\right\},\quad \hat{P}_G={-}Re_{G}\,\frac{\hat{d}}{{\rm i}k} \left\{3\,\frac{C_0}{d_0}-\varPi_L\,\frac{C_1}{d_0}-\frac{\varPi_L^2}{\varPi_u^2}\, \frac{\cos\left(\phi\right)}{Fr^2}\right\}. \end{gather}$$

Introducing (3.35) and (3.38) into (3.6), and linearizing once again, yields the dispersion relation for the spatial stability problem:

(3.41) \begin{align} {DR}&={-}{\rm i}\omega^2+{\rm i}k\omega\,\frac{17}{7}\, \frac{q_{L 0}}{h_0}-{\rm i}k^2\,\frac{9}{7}\,\frac{q_{L 0}^2}{h_0^2}+ \frac{5}{6}\,Fr^{{-}2}\left\{{\rm i}k^2\cos(\phi)\,h_0-k\sin(\phi)\right\} \nonumber\\ &\quad -{\rm i}^3k^4\,\frac{5}{6}\,We\,h_0+\frac{5}{2}\,\frac{1}{Re_L}\,\frac{1}{h_0^2} \left\{\omega-2k\,\frac{q_{L 0}}{h_0}\right\}+{\rm i}^2k^3\,\frac{6}{Re_L}\, \frac{q_{L 0}}{h_0}-{\rm i}^2k^2\omega\,\frac{9}{2}\,\frac{1}{Re_L} \nonumber\\ &\quad +\frac{\varPi_\mu \varPi_u}{\varPi_L} \left\{T_{G 0} \left[{\rm i}k\omega\,\frac{19}{336}\,h_0-{\rm i}^2k^3\,\frac{1}{6}\,\frac{1}{Re_L}\,h_0+{\rm i} k^2\, \frac{5}{112}\,q_{L 0}\right]+k\,\frac{5}{4}\,\frac{1}{Re_L}\,\frac{1}{\varPi_L}\, \frac{\hat{T}_G}{\hat{d}}\right\} \nonumber\\ &\quad +{\rm i}k^2\,\frac{19}{672}\,\frac{\varPi_\mu^2 \varPi_u^2}{\varPi_L^2}\,h_0^2T_{G 0}^2+k\, \frac{5}{6}\,\frac{\varPi_\rho \varPi_u^2}{\varPi_L}\,\frac{1}{Re_{G}} \left\{\partial_xP_{G 0}-\frac{h_0}{\varPi_L}\,{\rm i}k\, \frac{\hat{P}_G}{\hat{d}}\right\}=0, \end{align}

where $\hat {d}$ will cancel, due to $\hat {T}_G\propto \hat {d}$ and $\hat {P}_G\propto \hat {d}$, according to (3.39a,b) and (3.40a,b).

To compute nonlinear TWS, we recast (3.6a) into an ordinary differential equation in terms of the wave coordinate $\gamma =x-ct$:

(3.42a)$$\begin{gather} h'''={NL}(h,h',h'';\bar{h},c,q^{MF}_L), \end{gather}$$

(3.42b)$$\begin{gather}q^{MF}_L=q_L-hc=\bar{q}_L-\bar{h}c, \end{gather}$$

where primes denote differentiation with respect to $\gamma$, bars signify averaging over the wavelength $\varLambda$ in terms of $\gamma$, $c$ denotes the nonlinear wave speed, and the superscript ${MF}$ refers to the moving reference frame. Further, (3.42b) is the integral form of (3.6b), which we have used to eliminate $q$ from (3.42a). The system is closed through periodicity boundary conditions

(3.42c)$$\begin{gather} \left.{h^{(j)}}\right|_{\gamma=0}=\left.{h^{(j)}}\right|_{\gamma=\varLambda},\quad j=0,1,2, \end{gather}$$

and it is solved for a fixed value of $\bar {q}_L$, enforced through the integral condition

(3.43)$$\begin{gather} \varLambda^{{-}1}\int_0^\varLambda q_L\,{\rm d}\gamma=\bar{q}_L. \end{gather}$$

We do this numerically via the continuation software Auto07P, after recasting (3.42a) into a dynamical system. First, we continue the fixed-point solutions ($h'=h''=h''''=0$) of (3.42a) at $q_L=q_{L 0}$ and $h=h_0$ in terms of $c$, until reaching the Hopf bifurcation of the Kapitza instability. Then, starting from this point, periodic solutions are continued in terms of a selected control parameter, e.g. the liquid Reynolds number $Re_L$. The BVPs associated with the turbulent gas flow, (3.27) and (3.28), are solved simultaneously. In addition, we solve the linear dispersion relation (3.41) for the spatially most-amplified angular frequency $\omega _{max}$:

(3.44a,b)\begin{equation} {DR}(\omega_{max},k)=0,\quad \left.{\partial_\omega k_i}\right|_{\omega=\omega_{max}}=0. \end{equation}

By imposing $f=f_{max}=\omega _{max}/2/{\rm \pi}$, TWS most likely to emerge in an experiment can be tracked.

To compute the spatial evolution of nonlinear Kapitza waves, we solve (3.6a) and (3.6b) numerically on an open domain with inlet/outlet conditions. Details of the numerical scheme are given in Appendix F3 of Kalliadasis et al. (Reference Kalliadasis, Ruyer-Quil, Scheid and Velarde2012). In particular, we apply a second-order central-differences spatial discretization and a quasi-linearized Crank–Nicolson time integration. At the liquid outlet, we impose the soft boundary conditions of Richard, Ruyer-Quil & Vila (Reference Richard, Ruyer-Quil and Vila2016). At the liquid inlet, we prescribe explicitly $h$ and $q$ at the first two grid points ($i_x=1,2$), based on the primary flow:

(3.45a)$$\begin{gather} \left.h\right|_{i_x=1}=\left.h\right|_{i_x=2}=h_0, \end{gather}$$

(3.45b)$$\begin{gather}\left.q_L\right|_{i_x=1}=\left.q_L\right|_{i_x=2}=q_{L 0}\left[1+F(t)\right], \end{gather}$$

where the function $F(t)$ allows us to apply a tailored inlet forcing,

(3.46)\begin{equation} F(t)=\varepsilon_1\sin(2{\rm \pi} ft)+\varepsilon_2 \sum_{k=1}^{N} \sin(2{\rm \pi} k\,{\rm \Delta} f\,t+\varphi_{rand}),\quad {\rm \Delta} f=2f_{c}/N. \end{equation}

The first RHS term in (3.46) constitutes a harmonic perturbation of frequency $f$, and the second RHS term mimics white noise through a series of $N=1000$ Fourier modes that are shifted by a random phase shift $\varphi _{rand}=\varphi _{rand}(k)\in [0,2{\rm \pi} ]$ and span a frequency range of twice the linear cut-off frequency $f_{c}$ (Chang, Demekhin & Kalaidin Reference Chang, Demekhin and Kalaidin1996a). When $\varepsilon _1=0$, the inlet perturbation consists of only white noise. This setting will be used to simulate the natural, noise-driven evolution of a wavy film as it would occur in a real system. In other computations, we will apply additional coherent inlet forcing by setting $\varepsilon _1> 0$.

4.1. The WRIBL-OS approach

In our WRIBL-OS approach, the linear response of the liquid film is governed by the dispersion relation (3.41), but the perturbation amplitudes $\hat {T}_G$ and $\hat {P}_G$ are now obtained from the full (steady) RANS equations (3.8). For this, we recast (3.8) in terms of the curvilinear coordinates (3.15a,b) and introduce the gas stream function $\varPsi$,

(4.1a,b)\begin{equation} \tilde{u}=\frac{\bar{d}}{d}\,\partial_{\eta}\varPsi,\quad \tilde{v}={-}\partial_{\xi}\varPsi, \end{equation}

which we perturb, along with $p_G$ and $d$, around the primary flow (subscript $0$):

(4.2a)$$\begin{gather} \varPsi=\varPsi_0+\check{\varPsi}=\varPsi_0+\psi(\eta)\exp\left\{{\rm i}(k\xi-\omega t)\right\}, \end{gather}$$

(4.2b)$$\begin{gather}p_G=p_{G 0}+\check{p}_G=p_{G 0}+\hat{p}_G(\eta)\exp\left\{{\rm i}(k\xi-\omega t)\right\}, \end{gather}$$

(4.2c)$$\begin{gather}d=d_0+\check{d}=d_0+\hat{d}\exp\left\{{\rm i}(k\xi-\omega t)\right\}, \end{gather}$$

where $k=k_r\in \mathbb {R}$, and the time dependence is included formally to account for the unsteadiness of the liquid film. Upon linearization and subtraction of the primary flow, we obtain the linearized curvilinear RANS equations ${OS}_\xi$ in the $\xi$ direction,

(4.3a) \begin{align} & Re_{G}\,{\rm i}k\left\{\psi'\varPsi_0'-\frac{\hat{d}}{d_0}\, \varPsi_0'^2-\psi\varPsi_0''\right\}+Re_{G}\,\tilde{l}_{t}\left|\varPsi_0''\right| \left\{\tilde{l}_{t}'\left[k^2\left({-}2\psi+3\,\frac{\hat{d}}{d_0}\,\eta\varPsi_0'\right) \right.\right. \nonumber\\ &\qquad \left. -\,4\psi''+6\,\frac{\hat{d}}{d_0}\,\varPsi_0''\right]+\tilde{l}_{t} \left[{-}2\psi'''+6\,\frac{\hat{d}}{d_0}\,\varPsi_0'''-2\,\frac{\varPsi_0'''}{\varPsi_0''}\,\psi'' \right.\nonumber\\ &\qquad \left.\left. +\,k^2\left(-\psi'+\frac{3}{2}\,\frac{\hat{d}}{d_0}\,\varPsi_0'+\frac{3}{2}\, \frac{\hat{d}}{d_0}\,\eta\varPsi_0''-\frac{\varPsi_0'''}{\varPsi_0''}\,\psi+\frac{3}{2}\, \frac{\hat{d}}{d_0}\,\eta\,\frac{\varPsi_0'''\varPsi_0'}{\varPsi_0''}\right)\right]\right\} \nonumber\\ &\quad ={-}{\rm i\ }k\hat{p}_G+\psi'''-3\,\frac{\hat{d}}{d_0}\,\varPsi_0'''-k^2\left\{\psi'-\frac{\hat{d}}{d_0}\, \varPsi_0'-\frac{\hat{d}}{d_0}\,\eta\varPsi_0''\right\}, \end{align}

and ${OS}_\eta$ in the $\eta$ direction,

(4.3b) \begin{align} & Re_{G}\,k^2\left\{\psi\varPsi_0'-\frac{\hat{d}}{d_0}\,\eta\varPsi_0'^2\right\}+ Re_{G}\,\tilde{l}_{t}\left|\varPsi_0''\right|{\rm i}k\left\{2\tilde{l}_{t}'\,\frac{\hat{d}}{d_0}\, \eta\varPsi_0''+\tilde{l}_{t}\left[{-}2\psi''+2\,\frac{\hat{d}}{d_0}\,\varPsi_0''\right.\right. \nonumber\\ &\quad \left.\left.{}+2\eta\,\frac{\hat{d}}{d_0}\,\varPsi_0'''+k^2\left(-\psi+\frac{3}{2}\, \frac{\hat{d}}{d_0}\,\eta\varPsi_0'\right)\right]\right\}={-}\hat{p}'_G-{\rm i}k^3 \left\{-\psi+\frac{\hat{d}}{d_0}\,\eta\varPsi_0'\right\} \nonumber\\ &\quad\ -{\rm i}k\left\{\psi''-2\,\frac{\hat{d}}{d_0}\,\varPsi_0''-\frac{\hat{d}}{d_0}\,\eta\varPsi_0'''\right\}, \end{align}

where primes denote differentiation with respect to $\eta$. The pressure perturbation amplitude $\hat {p}_G$ can be removed from (4.3a) and (4.3b) via the final gas-side Orr–Sommerfeld equation

(4.3c)\begin{equation} \partial_{\eta}{OS}_\xi-{\rm i}k\,{OS}_\eta, \end{equation}

involving only $\psi$ and its derivatives. The problem is closed with the boundary conditions (3.24a,b):

(4.3d) \begin{equation} {\psi''}|_{\eta=0}=0,\quad \left.{\psi}\right|_{\eta=0}=0,\quad {\psi'}|_{\eta=d_0}=0,\quad \left.{\psi}\right|_{\eta=d_0}=0. \end{equation}

We solve (4.3) numerically for $\psi$ with the continuation software Auto07P, starting from the analytically tractable laminar long-wave limit ($\tilde {l}_{t}=k=0$). The amplitudes of the linear perturbations of the inter-phase coupling quantities,

(4.4a,b)\begin{equation} \check{T}_G=\hat{T}_G\exp\left\{{\rm i}(k\xi-\omega t)\right\},\quad \check{P}_G=\hat{P}_G\exp\left\{{\rm i}(k\xi-\omega t)\right\}, \end{equation}

can be obtained readily by recasting (3.14) in curvilinear coordinates, inserting (4.2), and linearizing:

(4.5a)$$\begin{gather} \hat{T}_G={-}\left.{\psi''}\right|_{\eta=d_0}+2\,\frac{\hat{d}}{d_0}\left.{\varPsi_0''}\right|_{\eta=d_0}, \end{gather}$$

(4.5b)$$\begin{gather}{\rm i}k\hat{P}_G=\left\{\left.{\psi'''}\right|_{\eta=d_0}-3\,\frac{\hat{d}}{d_0} \left.{\varPsi_0'''}\right|_{\eta=d_0}+k^2\hat{d}\left.{\varPsi_0''}\right|_{\eta=d_0}\right\}. \end{gather}$$

We point out that $\psi \propto \hat {d}$, thus $\hat {d}$ once again cancels from (3.41), as it should. Also, the spatial variations prescribed in (3.37) and (4.2) are equivalent in the linear limit $\hat {d}\to 0$, where the curvilinear coordinates collapse with the Cartesian ones. Thus $\partial _x\check {P}_G=\partial _\xi \check {P}_G$.

Figures 5(a) and 5(b) represent spatial linear stability predictions obtained with our WRIBL-OS approach, based on (3.41) and (4.3c), for parameters according to the experiments of Kofman et al. (Reference Kofman, Mergui and Ruyer-Quil2017) in an $H^\star =19$ mm channel. According to figure 5(a), the maximum of the growth rate dispersion curve increases with increasing counter-current gas flow rate, until forming a pinch point at $Re_G=-8490$, where the AI limit is reached (curve with crosses).

Figure 5. Spatial linear stability calculations with the WRIBL-OS approach. Parameters based on experiments of Kofman et al. (Reference Kofman, Mergui and Ruyer-Quil2017): $Ka=3174$ (water and air I in table 1), $\phi =5^\circ$, $Re_L=45$. Circles: passive-gas limit, $\varPi _\rho =\varPi _\mu =0$ in (3.41). Values of $|Re_G|$ increase in the order diamonds, squares, crosses, pentagons. (a,c) Full model, (b,d) $\varPi _\rho =0$ in (3.41), with (a,b) $H^\star =19$ mm, (c,d) $H^\star =40$ mm. Here: (a) $Re_G=[-6234, -8145, -8490\ (\textrm {AI})]$, dashed line $Re_G=-8500$; (b) $Re_G=[-6234, -8145, -8490, -15\,000]$; (c) $Re_G=[-8145, -15\,000, -20\,430 \ (\textrm {AI\ limit})]$, dashed line $Re_G=-20\,440$; (d) $Re_G=[-8145, -15\,000, -20\,430, -35\,000]$. We have rescaled $\omega$ with $\mathcal {T}_\nu =2.207\times 10^{-3}$ s (see (3.34a–c)).

Table 1. Fluid combinations used in our computations. The Kapitza number is defined as $Ka=\sigma/(\rho_L g^{1/3} \nu_L^{4/3})$, where $\sigma$, $\rho _L$ and $\nu _L$ denote the surface tension, density and kinematic viscosity of the liquid, and $g$ designates the gravitational acceleration.

This destabilization of the liquid film is caused by the inter-phase pressure coupling, as can be deduced by confronting figure 5(a) with figure 5(b), where we have represented corresponding growth rate curves in the limit $\varPi _\rho =0$. In that case, the gas effect enters only via $T_G$, and we observe a stabilization of the liquid film at large $|Re_G|$ (compare crosses and pentagons). Models that do not account for the gas pressure $P_G$ – e.g. the weak-confinement first-order WRIBL model of Tseluiko & Kalliadasis (Reference Tseluiko and Kalliadasis2011) – may thus give qualitatively incorrect linear stability predictions for the current configuration. The same observation also holds at weaker confinement, as shown by confronting figures 5(c) and 5(d), where we have chosen $H^\star =40$ mm.

As a result of our frozen-interface assumption ($\varPi _u{\gg }1$) expressed through the last two equations in (4.3d), one would obtain exactly the same linear stability problem (4.3) should one apply no-slip and no-penetration conditions at $y^\star =H^\star$ instead of a symmetry condition at $y^\star =D^\star$. This is because the primary gas flow would remain symmetrical about the centreline of the gas layer. Thus for all linear stability calculations based on the gas-side OS BVP (4.3), our symmetry condition (4.3d) is valid analytically.

4.2. The OS-OS approach

Linear stability calculations based on our WRIBL-LW and WRIBL-OS approaches may be limited to long-wave instability modes. To capture short-wave instability modes (§ 6.2), we introduce a stability formulation based on the full Navier–Stokes equations (4.1a,b) in the liquid and the full RANS equations (4.3) in the gas.

The gas-side linear response is governed by the same equations as in the WRIBL-OS approach, i.e. (4.3) and (4.5), and we focus here on deriving the equations governing the liquid-side linear response. For this, we perturb the film thickness as

(4.6)\begin{equation} h=h_0+\check{h}=h_0+\hat{h}\exp\left\{{\rm i}k(x-ct)\right\}, \end{equation}

assuming a temporal stability formulation this time, i.e. $k\in \mathbb {R}$ and $c=c_r+\textrm {i}c_i\in \mathbb {C}$, where $c_r$ denotes the wave speed, and $k c_i$ is the temporal growth rate.

We start with the full governing equations (4.1a,b). Considering these in the limit of fully developed flow with $h=h_0$ yields the liquid primary flow

(4.7a)$$\begin{gather} u_{L 0}=\frac{1}{2}\,Re_L\left\{\varPi_\rho\,\partial_xp_{G 0}- \frac{\sin(\phi)}{Fr^2}\right\}(y^2-2yh_0)+ \frac{\varPi_\mu\varPi_u}{\varPi_L}\,T_{G 0}y, \end{gather}$$

(4.7b)$$\begin{gather}\partial_y p_{L 0}={-}\frac{\cos(\phi)}{Fr^2}. \end{gather}$$

Next, we introduce the liquid stream function $\varPhi$,

(4.8a,b)\begin{equation} u_L=\partial_{y}\varPhi,\quad v_L={-}\partial_{x}\varPhi, \end{equation}

which we perturb around the primary flow:

(4.9)\begin{equation} \varPhi=\varPhi_0+\check{\varPhi}=\varPhi_0+\phi(y)\exp\left\{{\rm i}k(x-ct)\right\}. \end{equation}

Substituting (4.8a,b) and (4.9) into (4.1a,b), linearizing with respect to $\check {\varPhi }$, subtracting the primary flow, and applying standard manipulations, we obtain the liquid-side Orr–Sommerfeld equation

(4.10a)\begin{equation} \phi^{{\rm i}v}-2k^2\phi''+k^4\phi={\rm i}k\,Re_L\left\{\left(c-u_{L 0}\right) \left(k^2\phi-\phi''\right)-\phi u''_{L 0}\right\}, \end{equation}

the boundary conditions at $y=0$

(4.10b)\begin{equation} \phi''=\phi=0, \end{equation}

and the inter-phase coupling conditions at $y=h_0$,

(4.10c)\begin{align} &\phi u_{L 0}''+\tilde{c}\left\{\phi''+k^2\phi\right\}=\varPi_\mu\varPi_u\tilde{c}\hat{\underline{T}}_G, \end{align}

(4.10d)\begin{align} & -\frac{1}{Re_L}\left\{2k^2\phi u_{L 0}'+\tilde{c} \left[3 k^2\phi'-\phi'''\right]\right\}-{\rm i} k\tilde{c}\left\{-\tilde{c}\phi'-\phi u_{L 0}\right\} \nonumber\\ &\quad +{\rm i}k\phi p_{L 0}'={\rm i}k\,\frac{\varPi_\rho\varPi_u^2}{Re_{G}}\, \tilde{c}\hat{\underline{P}}_G+{\rm i} k^3\,We\,\phi, \end{align}

where primes denote differentiation with respect to $y$, and we have introduced $\tilde {c}=c-{u_{L 0}}|_{y=h_0}$. The nonlinearity involving $\tilde {c}$ in (4.10d) can be eliminated via (4.10c). Further, $\hat {\underline {T}}_G$ and $\hat {\underline {P}}_G$ are rescaled versions of the amplitudes in (4.5):

(4.11a,b)\begin{equation} \hat{\underline{T}}_G={-}\hat{h}\,\frac{\hat{T}_G}{\hat{d}},\quad \hat{\underline{P}}_G=\hat{h}\,\frac{\hat{P}_G}{\hat{d}}, \end{equation}

where $\hat {d}$ is an arbitrary deflection amplitude used in the solution of the gas-side problem (4.3), and $\hat {h}$ is linked directly to $\phi$ via the kinematic condition (3.1f),

(4.12)\begin{equation} \hat{h}=\frac{\left.{\phi}\right|_{y=h_0}}{\tilde{c}}. \end{equation}

The rescaling in (4.11a,b) allows us to solve the gas- and liquid-side problems sequentially.

We solve the two-phase BVP comprising (4.3) and (4.10) by expanding the stream function amplitudes $\phi$ and $\psi$ in terms of Chebyshev polynomials (Boomkamp et al. Reference Boomkamp, Boersma, Miesen and Beijnon1997; Barmak et al. Reference Barmak, Gelfgat, Vitoshkin, Ullman and Brauner2016b):

(4.13a,b)\begin{equation} \phi(\zeta)=c_{L 0}+\sum_{j=1}^{Np}c_{{l}j}\,T_j(\zeta),\quad \psi(\zeta)=c_{G 0}+\sum_{j=1}^{Np}c_{{g}j}\,T_j(\zeta), \end{equation}

where $T_j$ are $j$th-degree Chebyshev polynomials of the first kind, defined on the interval $\zeta \in [-1,1]$, with

(4.14)$$\begin{gather} \zeta=2\,\frac{y}{h_0}-1\quad \text{for}\ 0\le y\le h_0, \end{gather}$$

(4.15)$$\begin{gather}\zeta=1-2\,\frac{{\eta}}{d_0}\quad \text{for}\ 0\le {\eta}\le d_0. \end{gather}$$

Thus there are $2(N_p+1)$ unknown coefficients $c_{kj}$, which are fixed by the eight conditions in (4.10b), (4.10c), (4.10d) and (4.3d), and $2(N_p+1)-8$ additional constraints obtained by evaluating the ordinary differential equations (4.10a) and (4.3c) at the inner collocation points $\zeta _{2},\ldots,\zeta _{N_p-2}$, defined according to

(4.16)\begin{equation} \zeta_{i}=\cos\left[\frac{{\rm i}{\rm \pi}}{N_p}\right] \quad\text{for all}\ {\rm i}\in[0,N_p]. \end{equation}

Instead of solving for the coefficients $c_{kj}$, we solve directly for the $2(N_p+2)$ unknowns $\phi (\zeta _i)$ and $\psi (\zeta _i)$, arranged into the solution vectors

(4.17a,b)\begin{equation} \boldsymbol{\phi}=\left[\phi(\zeta_0),\ldots,\phi(\zeta_{N_p})\right]^{{\rm T}},\quad \boldsymbol{\psi}=\left[\psi(\zeta_0),\ldots,\psi(\zeta_{N_p})\right]^{{\rm T}}. \end{equation}

Then, by making use of the Chebyshev differentiation matrix $\boldsymbol{\mathsf{D}}$ (Trefethen Reference Trefethen2000),

(4.18)$$\begin{gather} \left[\phi^{(i)}(\zeta_0),\ldots,\phi^{(i)}(\zeta_{N_p})\right]^{{\rm T}}=\boldsymbol{\mathsf{D}}^{i} \left[\phi(\zeta_0),\ldots,\phi(\zeta_{N_p})\right]^{{\rm T}}, \end{gather}$$

(4.19)$$\begin{gather}\left[\psi^{(i)}(\zeta_0),\ldots,\psi^{(i)}(\zeta_{N_p})\right]^{{\rm T}}=\boldsymbol{\mathsf{D}}^{i} \left[\psi(\zeta_0),\ldots,\psi(\zeta_{N_p})\right]^{{\rm T}}, \end{gather}$$

where $i=1, 2, 3, 4$, and $(i)$ indicates the order of differentiation with respect to $\zeta$, (4.10) is cast into a generalized eigenvalue problem in matrix form,

(4.20) \begin{equation} \underline{\boldsymbol{\mathsf{A}}}\boldsymbol{\phi}=\tilde{c}\underline{\boldsymbol{\mathsf{B}}}\boldsymbol{\phi}, \end{equation}

and (4.3) is cast into a linear system of equations,

(4.21)\begin{equation} \underline{\boldsymbol{\mathsf{C}}}\boldsymbol{\psi}=\boldsymbol{b}, \end{equation}

introducing the coefficient matrices $\underline {\boldsymbol{\mathsf{A}}}$, $\underline {\boldsymbol{\mathsf{B}}}$ and $\underline {\boldsymbol{\mathsf{C}}}$, and the inhomogeneity $\boldsymbol {b}$. With the help of MATLAB (2015), we first solve (4.21) for $\boldsymbol {\psi }$ by numerical inversion via the / operator and then (4.20) for the eigenvalues $\tilde {c}$ and eigenvectors $\boldsymbol {\phi }$ via the eig function.

Using this approach, the full set of eigenmodes is computed at once. Thus short-wave instability modes, i.e. modes with $c_i\ne 0$ at $k=0$, can be obtained readily. Once a mode has been identified at a given wavenumber $k$, it can be tracked by advancing $k$, using the function eigs, which searches for eigenvalues in the vicinity of a previous solution.

In Appendix A, we validate our OS-OS approach, (4.20) and (4.21), versus Vellingiri et al. (Reference Vellingiri, Tseluiko and Kalliadasis2015) and Schmidt et al. (Reference Schmidt, Náraigh, Lucquiaud and Valluri2016). Figure 6 confronts temporal linear stability predictions from this approach (solid lines) with predictions from our WRIBL-OS approach (symbols), for parameters similar those in figure 5(a). Agreement is good up to $|Re_G|\sim 8000$. Thus our liquid-side WRIBL description suffices to predict the gas effect on the long-wave Kapitza instability.

Figure 6. Temporal stability predictions from the OS-OS (solid curves) and WRIBL-OS (open symbols) approaches. Parameters similar to those in figure 5(a): $Ka=3174$, $H^\star =19$ mm, $\phi =5^\circ$, $Re_L=32.7$. Circles indicate passive-gas limit ($\varPi _\rho =\varPi _\mu =0$ in (3.41)); pentagons indicate $Re_G=-4123$; squares indicate ${Re_G=-6173}$; diamonds indicate $Re_G=-8220$. (a) Growth rate; (b) wave speed.