3.1. Long-wave model

To derive an LW model, it is assumed that interfacial waves have a long wavelength when compared with the thickness of the film. Applying the classical LW expansion of the velocity and pressure fields, first proposed by Benney (Reference Benney1966), Tseluiko & Kalliadasis (Reference Tseluiko and Kalliadasis2011) obtain

(3.1) $$\begin{eqnarray}\displaystyle & & \displaystyle h_{t}+2h^{2}h_{x}+\left(\left(\frac{8}{15}\mathit{Re}\,h^{6}-{\displaystyle \frac{2}{3}}\cot {\it\theta}h^{3}\right)h_{x}+\frac{1}{3}\mathit{Re}\,\mathit{We}\,h^{3}h_{xxx}\right)_{x}\nonumber\\ \displaystyle & & \displaystyle \quad +\,\frac{\hat{{\it\tau}}_{w0}}{2}\left(h^{2}\left(1+{\displaystyle \frac{8\mathit{Re}}{15}}h^{3}h_{x}+{\it\nu}\hat{{\it\tau}}_{w1}[h]\right)\right)_{x}=0.\end{eqnarray}$$

By ignoring the terms multiplying $\hat{{\it\tau}}_{w0}$ in the parentheses that represent the effect of turbulent gas flow, we recover the classic Benney’s equation for the case of a falling film. Although this nonlinear partial differential equation is probably the simplest possible model that correctly predicts the linear stability threshold for long waves, its disadvantage is that it suffers finite-time blowup when $\mathit{Re}$ exceeds a limiting value.

3.2. Weighted integral-boundary-layer (WIBL) model

The shortcomings of the LW model can be overcome with the methodology developed by Ruyer-Quil & Manneville (Reference Ruyer-Quil and Manneville1998, Reference Ruyer-Quil and Manneville2000) to model falling films on inclined planes in the region of small-to-moderate Reynolds numbers. The starting point for this derivation is the first-order boundary-layer equations. Combination of these equations with a weighted-residual technique where the velocity components are expanded in terms of polynomial test functions satisfying the boundary conditions yields the residual. Application of this methodology to our problem and use of a simple Galerkin projection with just one test function yield a set of coupled equations for the evolution of film thickness ( $h$ ) and the local flow rate ( $q$ ):

(3.2) $$\begin{eqnarray}\displaystyle q_{t} & = & \displaystyle -\frac{17}{7}\frac{q\,q_{x}}{h}+\frac{9}{7}\frac{q^{2}\,h_{x}}{h^{2}}-\frac{5}{2\,\mathit{Re}}\frac{q}{h^{2}}+\frac{5}{3\,\mathit{Re}}h-\frac{5\cot {\it\theta}}{3\mathit{Re}}hh_{x}+\frac{5\,\mathit{We}}{6}hh_{xxx}\nonumber\\ \displaystyle & & \displaystyle +\,\hat{{\it\tau}}_{w0}\left(\frac{5}{4\,\mathit{Re}}-\frac{19\hat{{\it\tau}}_{w0}}{672}h^{2}h_{x}-\frac{19}{336}hq_{x}-\frac{5}{112}qh_{x}+\frac{5\,\text{sign}\left(\hat{{\it\tau}}_{w0}\right){\it\nu}}{4\mathit{Re}}\hat{{\it\tau}}_{w1}[h]\right)\!,\qquad\end{eqnarray}$$

(3.3) $$\begin{eqnarray}h_{t}+q_{x}=0.\end{eqnarray}$$

The terms multiplying $\hat{{\it\tau}}_{w0}$ correspond to the influence of the turbulent gas on the liquid film dynamics, and by setting $\hat{{\it\tau}}_{w0}$ to zero we recover the case of a falling film (Ruyer-Quil & Manneville Reference Ruyer-Quil and Manneville2000). This WIBL model predicts the correct value for the linear stability threshold, but is also applicable for moderately large values of $\mathit{Re}$ . We note that although in many applications involving gas–liquid flows the liquid Reynolds number is large, there are applications where it can take moderate values, for example in desiccant cooling systems (see, e.g. Ali, Vafai & Khaled Reference Ali, Vafai and Khaled2003; Mouza, Pantzali & Paras Reference Mouza, Pantzali and Paras2005) and in the cooling of microelectronic devices (see, e.g. Kabov et al. Reference Kabov, Lyulin, Marchuk and Zaitsev2007). In addition, there are several experiments that are preformed at moderate liquid Reynolds numbers (see, e.g. Semyonov Reference Semyonov1944; Clift, Pritchard & Nedderman Reference Clift, Pritchard and Nedderman1966; Zapke & Kröger Reference Zapke and Kröger2000b ).

4.1. Formulation of the OS problem

To obtain the base-state solution, we shall assume the flow to be steady and unidirectional with the interface being flat. The governing equations simplify to

(4.1a ) $$\begin{eqnarray}\displaystyle -\widehat{P}_{0x}+u_{0yy}+2 & = & \displaystyle 0,\end{eqnarray}$$

(4.1b ) $$\begin{eqnarray}\displaystyle \widehat{P}_{0y}+2\cot {\it\theta} & = & \displaystyle 0.\end{eqnarray}$$

(4.2a,b ) $$\begin{eqnarray}u_{0}(0)=0\quad \text{and}\quad u_{0y}(1)=\hat{{\it\tau}}_{w0}.\end{eqnarray}$$

By integrating equation (4.1a ) twice with respect to $y$ and applying the boundary conditions (4.2), we obtain

(4.3) $$\begin{eqnarray}u_{0}=\left({\displaystyle \frac{\widehat{P}_{0x}}{2}}-1\right)(y^{2}-2y)+\hat{{\it\tau}}_{w0}y.\end{eqnarray}$$

Further, by making use of the relation (2.20) and taking the channel to be much wider than the liquid film thickness, (4.3) simplifies to

(4.4) $$\begin{eqnarray}u_{0}=-(y^{2}-2y)+\hat{{\it\tau}}_{w0}y.\end{eqnarray}$$

The base-state stream function can be obtained by integrating (4.4) with respect to $y$ :

(4.5) $$\begin{eqnarray}{\it\psi}_{0}=-\left(\frac{y^{3}}{3}-y^{2}\right)+\frac{\hat{{\it\tau}}_{w0}}{2}y^{2}.\end{eqnarray}$$

We eliminate the pressure terms from the Navier–Stokes equations (2.8) and (2.9) by cross-differentiation, followed by making use of the definition of the stream function, $u={\it\Psi}_{y}$ , $v=-{\it\Psi}_{x}$ , to obtain

(4.6) $$\begin{eqnarray}\left(\frac{\partial }{\partial t}+\frac{\partial {\it\Psi}}{\partial y}\frac{\partial }{\partial x}-\frac{\partial {\it\Psi}}{\partial x}\frac{\partial }{\partial y}\right){\rm\nabla}^{2}{\it\Psi}=\frac{1}{\mathit{Re}}{\rm\nabla}^{4}{\it\Psi}.\end{eqnarray}$$

As the next step, we introduce infinitesimal perturbations to the base flow by writing ${\it\Psi}={\it\psi}_{0}+{\it\epsilon}{\it\psi}_{1}+O({\it\epsilon}^{2})$ , where ${\it\epsilon}\ll 1$ . Substitution of this in (4.6) and retention of only the linear terms yield

(4.7) $$\begin{eqnarray}\left(\frac{\partial }{\partial t}+u_{0}\frac{\partial }{\partial x}\right){\rm\nabla}^{2}{\it\psi}_{1}-u_{0}^{\prime \prime }\frac{\partial {\it\psi}_{1}}{\partial x}=\frac{1}{\mathit{Re}}{\rm\nabla}^{4}{\it\psi}_{1},\end{eqnarray}$$

where a prime denotes differentiation with respect to $y$ and $u_{0}$ is the base-state velocity given by $u_{0}={\it\psi}_{0}^{\prime }$ . Further, by writing the perturbation ${\it\psi}_{1}$ in the form of normal modes, ${\it\psi}_{1}={\it\phi}(y)\text{e}^{\text{i}k(x-ct)}$ , (4.7) results in the following classical OS equation:

(4.8) $$\begin{eqnarray}{\it\phi}^{\prime \prime \prime \prime }-2k^{2}{\it\phi}^{\prime \prime }+k^{4}{\it\phi}=\text{i}k\mathit{Re}[(u_{0}-c)({\it\phi}^{\prime \prime }-k^{2}{\it\phi})-u_{0}^{\prime \prime }{\it\phi}].\end{eqnarray}$$

The no-penetration and no-slip conditions at the bottom wall respectively imply

(4.9a,b ) $$\begin{eqnarray}{\it\phi}(0)=0,\quad {\it\phi}^{\prime }(0)=0.\end{eqnarray}$$

The perturbation to the flat interface can be represented as $y=h+{\it\epsilon}{\it\eta}(x,t)$ . Therefore, the waveform of the interface disturbance is given by ${\it\eta}(x,t)=A_{1}\text{e}^{\text{i}k(x-ct)}$ , where $A_{1}$ is a complex quantity. By substitution of this in (2.13), we obtain the linearised form of the kinematic condition,

(4.10) $$\begin{eqnarray}{\it\phi}(1)=(c-1-\hat{{\it\tau}}_{w0})A_{1}.\end{eqnarray}$$

From the tangential stress balance condition (2.14), we have

(4.11) $$\begin{eqnarray}{\it\phi}^{\prime \prime }(1)+k^{2}{\it\phi}(1)-[|\hat{{\it\tau}}_{w0}|{\it\nu}{\it\tau}_{w1}(k/{\it\nu})+2]A_{1}=0.\end{eqnarray}$$

The normal stress balance condition (2.17) becomes

(4.12) $$\begin{eqnarray}\displaystyle & & \displaystyle {\it\phi}^{\prime \prime \prime }(1)+k\left[\text{i}\mathit{Re}(c-1-\hat{{\it\tau}}_{w0})-3k\right]{\it\phi}^{\prime }(1)+\text{i}k\mathit{Re}\hat{{\it\tau}}_{w0}{\it\phi}(1)\nonumber\\ \displaystyle & & \displaystyle \quad -\,k[2\text{i}\cot {\it\theta}+\text{i}k^{2}\mathit{We}\mathit{Re}+2k\hat{{\it\tau}}_{w0}+\text{i}|\hat{{\it\tau}}_{w0}|{\it\nu}P_{w1}(k/{\it\nu})]A_{1}=0.\end{eqnarray}$$

The OS equation (4.8) together with the boundary conditions (4.9)–(4.12) constitute a generalised eigenvalue problem, with $c$ the eigenvalue.

4.2. Dispersion relations for low-dimensional models

To understand the linear stability of the flat-film solution ( $h=1$ ) in an LW model, we linearise (3.1) by substituting $h=1+f$ , where $|\,f|\ll 1$ , and retaining only the first-order terms we obtain

(4.13) $$\begin{eqnarray}f_{t}+af_{x}+bf_{xx}+cf_{xxxx}+d(\hat{{\it\tau}}_{w1}[\,f])_{x}=0,\end{eqnarray}$$

where

(4.14a−d ) $$\begin{eqnarray}a=2+\hat{{\it\tau}}_{w0},\quad b=\frac{8\mathit{Re}}{15}-\frac{2\cot {\it\theta}}{3}+\frac{4\mathit{Re}\hat{{\it\tau}}_{w0}}{15},\quad c=\frac{\mathit{Re}\mathit{We}}{3},\quad d=\frac{|\hat{{\it\tau}}_{w0}|{\it\nu}}{2}.\end{eqnarray}$$

The linearised LW equation (4.13) represents the evolution of an infinitesimal perturbation $f$ . We look for solutions in the form of normal modes, $f\propto \text{e}^{s(k)t+\text{i}kx}$ , to obtain the following dispersion relation:

(4.15) $$\begin{eqnarray}s(k)=-\text{i}ak+bk^{2}-ck^{4}-\text{i}dk{\it\tau}_{w1}\left(k/{\it\nu}\right).\end{eqnarray}$$

In a similar fashion, by substituting $h=1+f$ and $q=2/3+\hat{{\it\tau}}_{w0}/2+g$ , where $|\,f|\ll 1$ and $|g|\ll 1$ in (3.2) and (3.2), the linearised form of the WIBL model is obtained as

(4.16) $$\begin{eqnarray}f_{t}+g_{x}=0,\end{eqnarray}$$

(4.17) $$\begin{eqnarray}\displaystyle g_{t} & = & \displaystyle -\left(\frac{34}{21}+\frac{61}{48}\hat{{\it\tau}}_{w0}\right)g_{x}-\frac{5}{2\mathit{Re}}g+\frac{5\mathit{We}}{6}f_{xxx}+\left(\frac{4}{7}+\frac{139}{168}\hat{{\it\tau}}_{w0}+\frac{13}{48}\hat{{\it\tau}}_{w0}^{2}-\frac{5\cot {\it\theta}}{3\mathit{Re}}\right)f_{x}\nonumber\\ \displaystyle & & \displaystyle +\,\frac{5}{2\mathit{Re}}(2+\hat{{\it\tau}}_{w0})f+\frac{5\,\text{sign}(\hat{{\it\tau}}_{w0}){\it\nu}}{4\mathit{Re}}\hat{{\it\tau}}_{w0}\hat{{\it\tau}}_{w1}[\,f].\end{eqnarray}$$

As with the LW model, we look for solutions in the form of normal modes, $(\,f,g)\propto \text{e}^{s(k)t+\text{i}kx}$ , to obtain the following quadratic dispersion relation:

(4.18) $$\begin{eqnarray}\displaystyle & & \displaystyle s^{2}(k)+\frac{1}{336\mathit{Re}}\left[840+\text{i}k\mathit{Re}\left(427\hat{{\it\tau}}_{w0}+544\right)\right]s(k)+\frac{5\mathit{We}}{6}k^{4}\nonumber\\ \displaystyle & & \displaystyle \quad +\,\frac{1}{336\mathit{Re}}\left[-\mathit{Re}\left(192+278\hat{{\it\tau}}_{w0}+91\hat{{\it\tau}}_{w0}^{2}\right)+560\cot {\it\theta}\right]k^{2}\nonumber\\ \displaystyle & & \displaystyle \quad +\,\text{i}\frac{5}{4\mathit{Re}}\left[4+2\hat{{\it\tau}}_{w0}+\text{sign}\left(\hat{{\it\tau}}_{w0}\right)\hat{{\it\tau}}_{w0}{\it\nu}{\it\tau}_{w1}\left(k/{\it\nu}\right)\right]k=0.\end{eqnarray}$$

The above set of dispersion relations, (4.15), (4.18), along with the OS system, (4.8)–(4.12), will be utilised in the stability calculations that will be presented in the following sections.

5.1. Temporal stability

In order to analyse the temporal stability of the system, the wavenumber $k$ is taken to be real whereas $c$ is complex. This implies that the real part of $c$ gives the phase speed, whereas $kc_{i}$ is the growth rate of the perturbation, with $c_{i}$ being the imaginary part of the complex quantity $c$ .

Figure 3.

The growth rate, $s_{r}$ , versus the Reynolds number, $\mathit{Re}$ , for a falling liquid film when ${\it\Theta}=0$ , ${\it\theta}={\rm\pi}/4$ , $k=0.001$ and $\mathit{Ka}=1988.5$ .

5.1.1. Case I: no gas flow

To verify our numerical procedure, we first switch off the gas flow and examine whether we recover the linear stability results for a falling film (see, e.g. Yih Reference Yih1963). This is realised by setting $\hat{{\it\tau}}_{w0}=0$ , ${\it\tau}_{w1}(k/{\it\nu})=0$ and $P_{w1}(k/{\it\nu})=0$ . The asymptotic solution for small $k$ is given by

(5.3a ) $$\begin{eqnarray}\displaystyle {\it\phi} & = & \displaystyle (z+1)^{2}+k\Big[{\textstyle \frac{1}{240}}\text{i}\mathit{Re}z^{5}-{\textstyle \frac{1}{48}}\text{i}\mathit{Re}z^{4}+{\textstyle \frac{1}{6}}\left(-{\textstyle \frac{3}{4}}\text{i}\mathit{Re}+\text{i}\cot {\it\theta}\right)z^{3}\nonumber\\ \displaystyle & & \displaystyle +\,{\textstyle \frac{1}{2}}\left(-{\textstyle \frac{41}{120}}\text{i}\mathit{Re}+{\textstyle \frac{2}{3}}\text{i}\cot {\it\theta}+2\right)z^{2}+\left(-{\textstyle \frac{17}{240}}\text{i}\mathit{Re}+{\textstyle \frac{1}{6}}\text{i}\cot {\it\theta}+2\right)z+1\Big]+\cdots \,,\qquad\end{eqnarray}$$

(5.3b ) $$\begin{eqnarray}\displaystyle c & = & \displaystyle 2+\text{i}k\left({\textstyle \frac{8}{15}}\mathit{Re}-{\textstyle \frac{2}{3}}\cot {\it\theta}\right)+\cdots \,,\end{eqnarray}$$

while the growth rate is

(5.4) $$\begin{eqnarray}s_{r}=-\text{Re}(\text{i}kc)\equiv k^{2}\left({\textstyle \frac{8}{15}}\mathit{Re}-{\textstyle \frac{2}{3}}\cot {\it\theta}\right).\end{eqnarray}$$

The flat-film solution becomes unstable when $s_{r}>0$ , from which one can deduce the classical result that the critical Reynolds number at which the instability sets in is

(5.5) $$\begin{eqnarray}\mathit{Re}_{c}={\textstyle \frac{5}{4}}\cot {\it\theta}.\end{eqnarray}$$

Figure 3 shows the growth rate of the disturbance, $s_{r}$ , as a function of the Reynolds number, $\mathit{Re}$ , with the wavenumber $k=0.001$ when ${\it\theta}={\rm\pi}/4$ and $\mathit{Ka}=1988.5$ . Our computations by two different techniques, namely ‘fast Fourier transform (FFT)-based’ and ‘matrix-based’, agree very well with the asymptotic result. For details on the numerical scheme, the reader is referred to appendix A. Beyond the critical Reynolds number ( $\mathit{Re}_{c}\equiv 5\cot {\it\theta}/4$ ), the flat-film solution becomes unstable and the growth rate increases with increasing $\mathit{Re}$ . As we increase the wavenumber ( $k=0.01$ ) further, the growth rate increases for any given $\mathit{Re}$ , but the numerics still predicts the critical Reynolds number very well (as shown in figure 4).

Figure 4.

The growth rate, $s_{r}$ , versus the Reynolds number, $\mathit{Re}$ , for a falling liquid film when ${\it\Theta}=0$ , ${\it\theta}={\rm\pi}/4$ , $k=0.01$ and $\mathit{Ka}=1988.5$ .

Figure 5 shows the growth rates corresponding to the dominant normal mode for the case of a vertically falling liquid film ( ${\it\theta}={\rm\pi}/2$ ) for various $\mathit{Re}$ . As can be observed from the plot, the vertically falling liquid film is unstable to long-wavelength disturbances (for any $\mathit{Re}$ ), and the band of unstable wavenumbers increases as $\mathit{Re}$ is increased. As the wavenumber is increased beyond a particular value, the growth rate becomes negative, implying that the short-wave disturbances are damped.

Figure 5.

The growth rate, $s_{r}$ , versus the wavenumber, $k$ , for a falling liquid film when ${\it\Theta}=0$ , ${\it\theta}={\rm\pi}/2$ and $\mathit{Ka}=1988.5$ . The solid curves correspond to the solutions from the ‘FFT-based’ method and the open circles to the solutions from the ‘matrix-based’ technique.

Another way to assess the linear stability of the falling film system is to consider the neutral stability curves. In figure 6, we compare the neutral curves computed from our OS formulation (with the gas flow switched off) with those of Kalliadasis et al. (Reference Kalliadasis, Ruyer-Quil, Scheid and Velarde2012) for the parameter values $\mathit{Ka}=528.8$ and ${\it\theta}={\rm\pi}/2$ , and we find excellent agreement.

Figure 6.

The neutral stability curve for a vertically falling film ( ${\it\theta}={\rm\pi}/2$ ) when ${\it\Theta}=0$ and $\mathit{Ka}=528.8$ . The solid curve corresponds to our OS system, whereas the closed circles correspond to the results of Kalliadasis et al. (Reference Kalliadasis, Ruyer-Quil, Scheid and Velarde2012).

Finally, we compare the neutral curves obtained by solving the dispersion relations of the LW model (4.15) and the WIBL model (4.18) with the OS system for the vertical configuration ( ${\it\theta}={\rm\pi}/2$ ) and $\mathit{Ka}=1988.5$ . From figure 7, it can be observed that the LW and WIBL models correctly predict the linear stability threshold for $\mathit{Re}$ up to $\simeq 10$ , beyond which the agreement becomes poorer as the inertial effects become important, which are not captured by the low-dimensional models. It can also be observed that the results for the WIBL model are in much better agreement with the OS results than those of the LW model.

Figure 7.

The neutral stability curve for a vertically falling film ( ${\it\theta}={\rm\pi}/2$ ) when ${\it\Theta}=0$ and $\mathit{Ka}=1988.5$ . The solid curve corresponds to our OS system, the closed circles correspond to the LW model and the dashed curve corresponds to the WIBL model.

5.1.2. Case II: with gas flow

Following a few reasonable approximations in the gas flow (Lin et al. Reference Lin, Tseluiko and Kalliadasis2014), the asymptotic solution for the falling film in the presence of counter-current gas is given by

(5.6a ) $$\begin{eqnarray}\displaystyle {\it\phi} & = & \displaystyle (z+1)^{2}+k\Big[{\textstyle \frac{1}{480}}\text{i}\mathit{Re}(2+\hat{{\it\tau}}_{w0})z^{5}-{\textstyle \frac{1}{96}}\text{i}\mathit{Re}(2+\hat{{\it\tau}}_{w0})z^{4}-{\textstyle \frac{1}{48}}\text{i}\big(3\mathit{Re}(2+\hat{{\it\tau}}_{w0})\nonumber\\ \displaystyle & & \displaystyle +\,\cot {\it\theta}\big)z^{3}+{\textstyle \frac{1}{480}}(-\text{i}\mathit{Re}(2+\hat{{\it\tau}}_{w0})+160\text{i}\cot {\it\theta}+480)z^{2}\nonumber\\ \displaystyle & & \displaystyle +\,{\textstyle \frac{1}{480}}(-17\text{i}\mathit{Re}(2+\hat{{\it\tau}}_{w0})+80\text{i}\cot {\it\theta}+960)z+1\Big]+\cdots \,,\end{eqnarray}$$

(5.6b ) $$\begin{eqnarray}\displaystyle c & = & \displaystyle 2+\hat{{\it\tau}}_{w0}-2k|\hat{{\it\tau}}_{w0}|+\text{i}k\left({\displaystyle \frac{4\mathit{Re}}{15}}(2+\hat{{\it\tau}}_{w0})-{\displaystyle \frac{2}{3}}\cot {\it\theta}\right)+\cdots \,.\end{eqnarray}$$

Therefore, the growth rate is given by

(5.7) $$\begin{eqnarray}s_{r}=k^{2}\left({\displaystyle \frac{4\mathit{Re}}{15}}(2+\hat{{\it\tau}}_{w0})-{\displaystyle \frac{2}{3}}\cot {\it\theta}\right)\!.\end{eqnarray}$$

The details of the solution procedure are provided in appendix B. A comparison of the growth rates and phase speeds of the disturbance obtained from asymptotic and numerical approaches is displayed in figure 8, which shows a good agreement. However, the agreement deteriorates when the wavenumber of the disturbance waves is taken to be sufficiently large. Therefore, we resort to numerical treatment of the OS equations and the low-dimensional models.

Figure 8.

(a) The growth rate, $s_{r}$ , versus the gas parameter, ${\it\Theta}$ , and (b) the phase speed, $c_{r}$ , versus the gas parameter, ${\it\Theta}$ , for a falling liquid film in the presence of counter-current gas when ${\it\theta}={\rm\pi}/2$ , $k=10^{-5}$ , $N=0.4$ and $\mathit{Ka}=1988.5$ .

First, we present the neutral curves of the LW model (figure 9), the WIBL model (figure 10) and the OS problem (figure 11). We can observe that the stability boundaries obtained from the LW and WIBL models are slightly different from the OS results, especially for larger values of ${\it\Theta}$ . It is also evident that the results for the WIBL model are in better agreement with the OS results than the results for the LW model, at least for ${\it\Theta}$ up to the value 2.5. We can observe that for the LW model, the region of stability increases with increasing gas shear for larger values of $\mathit{Re}$ . The behaviour is similar for the WIBL model but starts to change for larger values of ${\it\Theta}$ . The OS results also show similar behaviour for smaller values of ${\it\Theta}$ , but this behaviour drastically changes for larger values of ${\it\Theta}$ . This qualitative difference between the results for the reduced model equations and the OS results is attributed to the fact that a consistent derivation of the reduced models (see Tseluiko & Kalliadasis Reference Tseluiko and Kalliadasis2011) is made under assumptions which imply that the normal stress imposed by the turbulent gas is of higher order and is therefore neglected, whereas for the considered helium–methanol system and for sufficiently large values of the gas parameter ${\it\Theta}$ , it turns out that the energy production due to the normal stress imposed by the turbulent gas dominates over the energy production due to the tangential stress imposed by the turbulent gas (as is evident from § 5.2).

Figure 9.

The neutral stability curves from the LW model for a vertically falling film in contact with counter-current gas flow for various values of the gas shear ${\it\Theta}$ , when ${\it\theta}={\rm\pi}/2$ , $\mathit{Ka}=1988.5$ and $N=0.4$ .

Figure 10.

The neutral stability curves from the WIBL model for a vertically falling film in contact with counter-current gas flow for various values of the gas shear ${\it\Theta}$ , when ${\it\theta}={\rm\pi}/2$ , $\mathit{Ka}=1988.5$ and $N=0.4$ .

Figure 11.

The neutral stability curves from the OS equations for a vertically falling film in contact with counter-current gas flow for various values of the gas shear ${\it\Theta}$ , when ${\it\theta}={\rm\pi}/2$ , $\mathit{Ka}=1988.5$ and $N=0.4$ .

To achieve a better comparison between the various models, we plot the stability boundary in the $(k,\,{\rm\Theta})$ -plane for two different values of the Reynolds number, $\mathit{Re}=5$ and 18. It can be seen from figure 12 that both the LW and the WIBL model show good agreement with the OS results for $\mathit{Re}=5$ , but as the Reynolds number is increased to $\mathit{Re}=18$ , the WIBL model overall shows better agreement with the OS results than the LW model.

Figure 12.

The neutral stability curves from the different formulations: the solid curves correspond to the OS equation, the open circles to the LW model and the closed circles to the WIBL model, when ${\it\theta}={\rm\pi}/2$ , $\mathit{Ka}=1988.5$ and $N=0.4$ . The black and grey correspond to $\mathit{Re}=5$ and $\mathit{Re}=18$ , respectively.

5.2. Energy analysis

To investigate the mechanism of instability in the gas–liquid system, we adopt a mechanical energy balance formulation (Smith & Davis Reference Smith and Davis1982; Hooper & Boyd Reference Hooper and Boyd1983; Kelly et al. Reference Kelly, Goussis, Lin and Hsu1989). The energy balance equation is derived by multiplying the linear equation for the perturbation evolution by $(u_{1},v_{1})$ and further by integrating across the interface thickness and over a wavelength $(L_{x}=2{\rm\pi}/k_{m})$ , where $k_{m}$ is the wavenumber corresponding to the maximum growth rate. The main ideas in the derivation are given below. The linearised equations and the boundary conditions governing the evolution of perturbations are given by

(5.8) $$\begin{eqnarray}\displaystyle & \mathit{Re}\left(u_{1t}+u_{0}u_{1x}+v_{1}u_{0y}\right)=-p_{1x}+u_{1xx}+u_{1yy}, & \displaystyle\end{eqnarray}$$

(5.9) $$\begin{eqnarray}\displaystyle & \mathit{Re}\left(v_{1t}+u_{0}v_{1x}\right)=-p_{1y}+v_{1xx}+v_{1yy}, & \displaystyle\end{eqnarray}$$

(5.10) $$\begin{eqnarray}\displaystyle & u_{1x}+v_{1y}=0, & \displaystyle\end{eqnarray}$$

(5.11) $$\begin{eqnarray}\displaystyle & u_{1}=v_{1}=0,\quad \text{at}~y=0, & \displaystyle\end{eqnarray}$$

(5.12) $$\begin{eqnarray}\displaystyle & {\it\eta}_{t}+\left(1+\hat{{\it\tau}}_{w0}\right){\it\eta}_{x}=v_{1},\quad \text{at}~y=1, & \displaystyle\end{eqnarray}$$

(5.13) $$\begin{eqnarray}\displaystyle & u_{1y}+v_{1x}=[|\hat{{\it\tau}}_{w0}|{\it\nu}{\it\tau}_{w1}(k/{\it\nu})+2]{\it\eta},\quad \text{at}~y=1, & \displaystyle\end{eqnarray}$$

(5.14) $$\begin{eqnarray}p_{1}=[|\hat{{\it\tau}}_{w0}|{\it\nu}P_{w1}(k/{\it\nu})+2\cot {\it\theta}]{\it\eta}-\mathit{Re}\mathit{We}{\it\eta}_{xx}-2\hat{{\it\tau}}_{w0}{\it\eta}_{x}-2u_{1x},\quad \text{at}~y=1,\end{eqnarray}$$

where ${\it\eta}$ corresponds to the deflection in the free surface about $h=1$ . By multiplying (5.8) by $u_{1}$ and (5.9) by $v_{1}$ , adding the results and integrating over the film thickness and the wavelength ( $L_{x}$ ) we obtain

(5.15) $$\begin{eqnarray}\displaystyle & & \displaystyle {\displaystyle \frac{\mathit{Re}}{2L_{x}}}{\displaystyle \frac{\text{d}}{\text{d}t}}\int _{0}^{1}\int _{0}^{L_{x}}\left(u_{1}^{2}+v_{1}^{2}\right)\text{d}x\text{d}y\nonumber\\ \displaystyle & & \displaystyle \quad =-{\displaystyle \frac{\mathit{Re}}{L_{x}}}\int _{0}^{1}\int _{0}^{L_{x}}u_{0y}u_{1}v_{1}\text{d}x\text{d}y-{\displaystyle \frac{1}{L_{x}}}\int _{0}^{1}\int _{0}^{L_{x}}u_{1}p_{1x}\text{d}x\text{d}y-{\displaystyle \frac{1}{L_{x}}}\int _{0}^{1}\int _{0}^{L_{x}}v_{1}p_{1y}\text{d}x\text{d}y\nonumber\\ \displaystyle & & \displaystyle \qquad +\,{\displaystyle \frac{1}{L_{x}}}\int _{0}^{1}\int _{0}^{L_{x}}(u_{1}u_{1xx}+u_{1}u_{1yy}+v_{1}v_{1xx}+v_{1}v_{1yy})\text{d}x\text{d}y.\end{eqnarray}$$

Integration of (5.15) by parts followed by substitution of $p_{1}$ from (5.14) results in

(5.16) $$\begin{eqnarray}\displaystyle & & \displaystyle {\displaystyle \frac{\mathit{Re}}{2L_{x}}}{\displaystyle \frac{\text{d}}{\text{d}t}}\int _{0}^{1}\int _{0}^{L_{x}}\left(u_{1}^{2}+v_{1}^{2}\right)\text{d}x\text{d}y\nonumber\\ \displaystyle & & \displaystyle \quad =-{\displaystyle \frac{\mathit{Re}}{L_{x}}}\int _{0}^{1}\int _{0}^{L_{x}}u_{0y}u_{1}v_{1}\text{d}x\text{d}y\nonumber\\ \displaystyle & & \displaystyle \qquad -\,{\displaystyle \frac{1}{L_{x}}}\int _{0}^{L_{x}}v_{1}\big[(|\hat{{\it\tau}}_{w0}|{\it\nu}P_{w1}(k/{\it\nu})+2\cot {\it\theta}){\it\eta}-\mathit{Re}\mathit{We}{\it\eta}_{xx}-2\hat{{\it\tau}}_{w0}{\it\eta}_{xx}-2u_{1x}\big]\big|_{y=1}\text{d}x\nonumber\\ \displaystyle & & \displaystyle \qquad +\,{\displaystyle \frac{1}{L_{x}}}\int _{0}^{1}\int _{0}^{L_{x}}\left(u_{1}u_{1xx}+u_{1}u_{1yy}+v_{1}v_{1xx}+v_{1}v_{1yy}\right)\text{d}x\text{d}y.\end{eqnarray}$$

The above equation can be further simplified by substituting the square of the continuity equation (5.10). This yields

(5.17) $$\begin{eqnarray}E_{K}=E_{ST}+E_{GP}+E_{R}+E_{SH}+E_{D}+E_{T1},\end{eqnarray}$$

where

(5.18) $$\begin{eqnarray}E_{K}={\displaystyle \frac{1}{2L_{x}}}{\displaystyle \frac{\text{d}}{\text{d}t}}\int _{0}^{1}\int _{0}^{L_{x}}\left(u_{1}^{2}+v_{1}^{2}\right)\text{d}x\text{d}y\end{eqnarray}$$

is the rate of change of ‘kinetic energy’ of the disturbance, which is proportional to the growth rate. The disturbance energy production due to surface tension forces is given by

(5.19) $$\begin{eqnarray}E_{ST}={\displaystyle \frac{\mathit{We}}{L_{x}}}\int _{0}^{L_{x}}v_{1}\big|_{1}{\it\eta}_{xx}\text{d}x.\end{eqnarray}$$

The rate of work done by the pressure gradient (due to gas flow and gravity) in deforming the interface is

(5.20) $$\begin{eqnarray}E_{GP}=-{\displaystyle \frac{1}{L_{x}\mathit{Re}}}\int _{0}^{L_{x}}v_{1}\big|_{1}[|\hat{{\it\tau}}_{w0}|{\it\nu}P_{w1}(k/{\it\nu})+2\cot {\it\theta}]{\it\eta}\text{d}x,\end{eqnarray}$$

whereas the term $E_{R}$ corresponds to the rate of energy transfer from the mean flow of the liquid to the disturbance by means of Reynolds stresses, and is given by

(5.21) $$\begin{eqnarray}E_{R}=-{\displaystyle \frac{1}{L_{x}}}\int _{0}^{1}\int _{0}^{L_{x}}u_{0y}u_{1}v_{1}\text{d}x\text{d}y.\end{eqnarray}$$

The work done by the disturbance shear force at the interface is represented by

(5.22) $$\begin{eqnarray}E_{SH}={\displaystyle \frac{1}{L_{x}\mathit{Re}}}\int _{0}^{L_{x}}u_{1}\big|_{1}\left[u_{1y}+v_{1x}\right]_{1}\text{d}x,\end{eqnarray}$$

whereas the rate of dissipation of energy by viscous forces is

(5.23) $$\begin{eqnarray}E_{D}=-{\displaystyle \frac{1}{L_{x}\mathit{Re}}}\int _{0}^{1}\int _{0}^{L_{x}}[2(u_{1x})^{2}+\left(u_{1y}+v_{1x}\right)^{2}+2\left(v_{1y}\right)^{2}]\text{d}x\text{d}y.\end{eqnarray}$$

Finally, the work done by the applied shear stress at the interface due to counter-current gas flow is given by

(5.24) $$\begin{eqnarray}E_{T1}={\displaystyle \frac{2\hat{{\it\tau}}_{w0}}{L_{x}\mathit{Re}}}\int _{0}^{L_{x}}v_{1}\big|_{1}{\it\eta}_{x}\text{d}x.\end{eqnarray}$$

The various energy quantities (5.18)–(5.24) in the kinetic energy equation (5.17) can be evaluated by making use of the solution obtained by solving the eigenvalue problem (A 7). The variables involved in the energy terms are derived once the eigenvalue $c$ and the eigenfunction ${\it\phi}$ are obtained by solving the OS problem. In order to perform the integration along the cross-streamwise direction, we employ a pseudo-spectral approach, whereas integration along the streamwise direction is carried out using the trapezoidal rule.

In this energy analysis, we are interested in the energy composition for the most dominant normal mode; therefore we consider wavenumbers corresponding to the maximum growth rate. The various energy terms in the kinetic energy equation are evaluated for increasing values of the gas flow rate ${\it\Theta}$ using the numerical methodology described above, and the results, scaled with the kinetic energy $E_{K}$ , are shown in figure 13. For the case of a vertically falling film in the absence of gas flow ( ${\it\Theta}=0$ ), the main contribution to energy production comes from the disturbance shear force at the interface $\widehat{E}_{SH}$ , whereas some of this energy is dissipated by the viscous forces $\widehat{E}_{D}$ . As the gas flow is switched on, the energy production has a contribution from the applied shear stress $\widehat{E}_{T1}$ in addition to the applied pressure gradient due to gas flow and gravity $\widehat{E}_{GP}$ . (We note that for the vertical configuration of the channel considered here, the contribution due to gravity in $\widehat{E}_{GP}$ vanishes.) It can be observed that $\widehat{E}_{T1}$ dominates over $\widehat{E}_{GP}$ for smaller values of the gas flow rate, whereas $\widehat{E}_{GP}$ is dominant for larger gas flow rates. The contribution from the Reynolds stresses $\widehat{E}_{R}$ is negligible for the moderate values of the Reynolds number considered in this study.

Figure 13.

Plots of various energy terms in (5.17) normalised with the total kinetic energy $E_{K}$ as a function of the gas shear ${\it\Theta}$ , when $\mathit{Re}_{0}=10$ , ${\it\theta}={\rm\pi}/2$ , $\mathit{Ka}=1988.5$ and $N=0.4$ .

5.3. Absolute and convective instabilities

In this section, we outline our methodology to analyse absolute and convective instabilities by considering the dispersion relation obtained from the LW model. Our analysis follows that by Fokas & Papageorgiou (Reference Fokas and Papageorgiou2005). We note first that the dispersion relation $s(k)$ is not analytic at zero and, therefore, it is not possible to study absolute and convective instabilities in a straightforward way by analysing the saddle point of the real part of the dispersion relation and applying the Briggs pinching criterion, see e.g. Huerre & Monkewitz (Reference Huerre and Monkewitz1990). To begin, let us assume that $f(x,\,0)=f_{0}(x)$ is the given initial condition. Therefore, the solution of the linearised equation can be written as

(5.25) $$\begin{eqnarray}f(x,\,t)=\frac{1}{2{\rm\pi}}\int _{-\infty }^{\infty }\text{e}^{s(k)t+\text{i}kx}\hat{f}_{0}(k)\hspace{2.22198pt}\text{d}k,\end{eqnarray}$$

where $\hat{f}_{0}(k)$ is the Fourier transform of the initial condition:

(5.26) $$\begin{eqnarray}\hat{f}_{0}(k)=\int _{-\infty }^{\infty }\text{e}^{-\text{i}kx}f_{0}(x)\hspace{2.22198pt}\text{d}x.\end{eqnarray}$$

There are two possibilities. If $\text{Re}[s(k)]\leqslant 0$ for all $k$ , the flat-film solution is linearly stable. On the other hand, if $\text{Re}[s(k)]>0$ at least for some $k$ , then the flat-film solution is unstable, and the instability could be either absolute or convective. Noting that ${\it\tau}_{w1}(-{\it\alpha})=\overline{{\it\tau}_{w1}({\it\alpha})}$ for real ${\it\alpha}$ , where the bar denotes complex conjugation (see Vellingiri et al. Reference Vellingiri, Tseluiko, Savva and Kalliadasis2013), the growth rate can be written as

(5.27) $$\begin{eqnarray}s(k)=\left\{\begin{array}{@{}ll@{}}-\text{i}ak+bk^{2}-ck^{4}-\text{i}\text{d}k\hspace{2.22198pt}\overline{{\it\tau}_{w1}(k/{\it\nu})}\equiv s_{-}(k)\quad & \text{if}~k<0,\\ -\text{i}ak+bk^{2}-ck^{4}-\text{i}\text{d}k\hspace{2.22198pt}{\it\tau}_{w1}(k/{\it\nu})\equiv s_{+}(k)\quad & \text{if}~k\geqslant 0.\end{array}\right.\end{eqnarray}$$

The growth rates are related by $s_{-}(-k)=\overline{s_{+}(k)}$ . We denote the real part of the growth rate as $s_{r}(k)\equiv \text{Re}[s_{+}(k)]=bk^{2}-ck^{4}+\text{d}k\hspace{2.22198pt}\text{Im}[{\it\tau}_{w1}(k/{\it\nu})]$ . The integral (5.25) can be rewritten as

(5.28) $$\begin{eqnarray}f(x,\,t)=\frac{1}{2{\rm\pi}}[I(x,\,t)+J(x,\,t)],\end{eqnarray}$$

where

(5.29) $$\begin{eqnarray}I(x,\,t)=\left(\int _{-\infty }^{-k_{1}}+\int _{k_{1}}^{\infty }\right)\text{e}^{s(k)t+\text{i}kx}\hat{f}_{0}(k)\hspace{2.22198pt}\text{d}k\end{eqnarray}$$

and

(5.30) $$\begin{eqnarray}J(x,\,t)=\left(\int _{-k_{1}}^{0}+\int _{0}^{k_{1}}\right)\text{e}^{s(k)t+\text{i}kx}\hat{f}_{0}(k)\hspace{2.22198pt}\text{d}k,\end{eqnarray}$$

with $k_{1}$ being the critical wavenumber below which the flat-film solution is unstable. For fixed $x/t$ , $I(x,\,t)$ goes to zero as $t$ increases. Therefore, it is sufficient to analyse the other integral, $J(x,\,t)$ :

(5.31) $$\begin{eqnarray}\displaystyle J(x,\,t) & = & \displaystyle \int _{-k_{1}}^{0}\text{e}^{s_{-}(k)t+\text{i}kx}\hat{f}_{0}(k)\hspace{2.22198pt}\text{d}k+\int _{0}^{k_{1}}\text{e}^{s_{+}(k)t+\text{i}kx}\hat{f}_{0}(k)\hspace{2.22198pt}\text{d}k\end{eqnarray}$$

(5.32) $$\begin{eqnarray}\displaystyle & = & \displaystyle \int _{0}^{k_{1}}\text{e}^{s_{-}(-k)t-\text{i}kx}\hat{f}_{0}(-k)\hspace{2.22198pt}\text{d}k+\int _{0}^{k_{1}}\text{e}^{s_{+}(k)t+\text{i}kx}\hat{f}_{0}(k)\hspace{2.22198pt}\text{d}k.\end{eqnarray}$$

By making use of the fact that $s_{-}(-k)=\overline{s_{+}(k)}$ for real $k$ , and that $\hat{f}_{0}(-k)=\overline{\hat{f}_{0}(k)}$ as $f_{0}(x)$ is a real function, we obtain

(5.33) $$\begin{eqnarray}\displaystyle J(x,\,t) & = & \displaystyle \int _{0}^{k_{1}}\left[\text{e}^{\overline{s_{+}(k)}t-\text{i}kx}\overline{\hat{f}_{0}(k)}+\text{e}^{s_{+}(k)t+\text{i}kx}\hat{f}_{0}(k)\right]\hspace{2.22198pt}\text{d}k\end{eqnarray}$$

(5.34) $$\begin{eqnarray}\displaystyle & = & \displaystyle \int _{0}^{k_{1}}\left[\,\overline{\text{e}^{s_{+}(k)t+\text{i}kx}\hat{f}_{0}(k)}+\text{e}^{s_{+}(k)t+\text{i}kx}\hat{f}_{0}(k)\right]\hspace{2.22198pt}\text{d}k\end{eqnarray}$$

(5.35) $$\begin{eqnarray}\displaystyle & = & \displaystyle 2\text{Re}\left[K(x,\,t)\right],\end{eqnarray}$$

where

(5.36) $$\begin{eqnarray}K(x,\,t)=\int _{0}^{k_{1}}\text{e}^{s_{+}(k)t+\text{i}kx}\hat{f}_{0}(k)\hspace{2.22198pt}\text{d}k.\end{eqnarray}$$

The function $s_{+}(k)$ is defined for real non-negative values of $k$ , since ${\it\tau}_{w1}({\it\alpha})$ is obtained by solving equation (3.34) in Tseluiko & Kalliadasis (Reference Tseluiko and Kalliadasis2011) with homogeneous boundary conditions (3.37) and (3.38) in the same study for positive values of ${\it\alpha}$ . Noticing that the coefficients in (3.34) in Tseluiko & Kalliadasis (Reference Tseluiko and Kalliadasis2011) are analytic functions of ${\it\alpha}$ , we can obtain an analytic continuation of ${\it\tau}_{w1}({\it\alpha})$ from the non-negative real semi-axis to the whole complex plane by solving equation (3.34) in Tseluiko & Kalliadasis (Reference Tseluiko and Kalliadasis2011) with homogeneous boundary conditions (3.37) and (3.38) in the same study when ${\it\alpha}$ is an arbitrary complex number. Then the analysis of $K(x,\,t)$ depends on the topology of the curves $s_{r}(k)=0$ in the complex $k$ -plane, where $s_{r}(k)$ denotes the real part of $s_{+}(k)$ . There are two possibilities.

Case I: there is part of the curve $s_{r}(k)=0$ that connects the origin $(0,0)$ and $(k_{1},0)$ in the $(k_{r},k_{i})$ -plane, with $k_{r}$ and $k_{i}$ denoting the real and imaginary parts of $k$ respectively. In addition, if $\hat{f}_{0}(k)$ can be analytically continued off the real axis so that $\hat{f}_{0}(k)$ is analytic in a neighbourhood of the curve $s_{r}(k)=0$ connecting the points $(0,0)$ and $(k_{1},0)$ (this is the case, for example, when $f_{0}(x)$ has a sufficiently fast exponential decay as $|x|\rightarrow \infty$ ), the integration contour for $K(x,t)$ can be deformed to follow that part of the curve $s_{r}(k)=0$ that connects $(0,0)$ to $(k_{1},0)$ . The Riemann–Lebesgue lemma then implies that for a fixed $x$ , $K(x,t)\rightarrow 0$ as $t\rightarrow \infty$ , and we have convective instability.

Case II: the other possibility is that it is impossible to connect the origin with $(k_{1},\,0)$ using part of the curve $s_{r}(k)=0$ in the $(k_{r},k_{i})$ -plane. In such a scenario, there must exist a saddle point (also called the pinch point) of $s_{r}(k)$ , at which $s_{r}(k)$ is positive, so that the integration contour can be deformed into a steepest descent path passing through this saddle point. Here, $K(x,\,t)$ is seen to grow for a fixed $x$ as $t\rightarrow \infty$ (Hinch Reference Hinch1991), and the type of instability in this case is absolute. We also note here that the saddle point corresponds to the mode with zero group velocity.

The details of the numerical scheme to compute the curves $s_{r}(k)=0$ in the $(k_{r},k_{i})$ -plane for the OS problem are provided in appendix C. We begin by presenting in figure 14 the curves $s_{r}(k)=0$ in the $(k_{r},k_{i})$ -plane for various values of the gas shear, ${\it\Theta}$ , for $\mathit{Re}=5$ , corresponding to the dispersion relation (4.15) obtained from the LW model. It can be observed from figure 14 that for sufficiently small values of ${\it\Theta}$ , the origin is connected to $(k_{r},0)$ . Therefore, we have convective instability. This is to be expected, as any localised disturbance introduced to the system is convected downstream, in the positive $x$ direction. The situation is similar for large values of ${\it\Theta}$ , but the disturbance is convected in the negative $x$ direction. However, there exists a band of gas flow rates ( $1.235\leqslant {\it\Theta}\leqslant 2.1$ ) for which the origin is not connected to $(k_{r},0)$ by $s_{r}(k)=0$ , and the instability is absolute.

Figure 14.

The curves of $s_{r}(k)=0$ in the complex wavenumber plane $(k_{r},k_{i})$ from the LW model for various values of the gas shear, ${\it\Theta}$ , when $\mathit{Re}_{0}=10$ , ${\it\theta}={\rm\pi}/2$ , $\mathit{Ka}=1988.5$ and $N=0.4$ .

Figure 15 shows the curves $s_{r}(k)=0$ corresponding to the quadratic dispersion relation (4.18) obtained from the WIBL model. We note that for the WIBL model there exist two branches, $s_{r1}(k)$ and $s_{r2}(k)$ . However, for the parameter space under consideration, only for one of these branches does there exist a band of unstable wavenumbers. We denote this branch by $s_{r}(k)$ . The results obtained using the WIBL model are similar to those obtained from the LW model with only minor differences.

Figure 15.

The curves of $s_{r}(k)=0$ in the complex wavenumber plane $(k_{r},k_{i})$ from the WIBL model for various values of the gas shear, ${\it\Theta}$ , when $\mathit{Re}_{0}=10$ , ${\it\theta}={\rm\pi}/2$ , $\mathit{Ka}=1988.5$ and $N=0.4$ .

Finally, we present the curves $s_{r}(k)=0$ in the $(k_{r},k_{i})$ -plane corresponding to the dominant mode in the OS problem (figure 16), which also turns out to be the only unstable mode for the parameter values considered. The results are again similar to those presented earlier for the LW and WIBL models, with the region of absolute instability predicted to be $1.21\leqslant {\it\Theta}\leqslant 2.1$ . This is not very surprising as the liquid flow rate that we have considered is sufficiently small for the low-dimensional models to have good agreement with the OS results (see also figure 12, from the temporal analysis). Moreover, we find that the WIBL model has a better quantitative agreement with the OS equations than the LW model in the flow regime considered.

Figure 16.

The curves of $s_{r}(k)=0$ in the complex wavenumber plane $(k_{r},k_{i})$ by solving the OS problem for various values of the gas shear, ${\it\Theta}$ , when $\mathit{Re}_{0}=10$ , ${\it\theta}={\rm\pi}/2$ , $\mathit{Ka}=1988.5$ and $N=0.4$ .

We note, in addition, that the nonlinear computations of the WIBL model of Tseluiko & Kalliadasis (Reference Tseluiko and Kalliadasis2011) reveal that the flooding point (defined as the point at which there appear large-amplitude nearly stationary waves) is close to ${\it\Theta}\sim 2.5$ . Interestingly, this flooding point is also close to the upper limit of the absolute instability regime ( ${\it\Theta}\sim 2.1$ ) predicted by our analysis.

5.3.1. Influence of the inclination angle

We analyse the influence of the inclination angle on the onset of absolute and convective instabilities by performing the analysis for the OS problem for various values of the inclination angle ${\it\theta}$ , namely ${\rm\pi}/3$ , ${\rm\pi}/4$ and ${\rm\pi}/6$ , with the modified Reynolds number being fixed at $\mathit{Re}_{0}=10$ . Table 1 shows the regimes of absolute and convective instabilities for the various angles (the results for ${\rm\pi}/2$ are shown for the various models, whereas for the other angles the results are obtained by performing the analysis for the OS problem only). It should be noted that as the inclination angle of the channel is decreased, the liquid flow rate decreases as well and the stabilisation of the interface is enhanced by the normal component of gravity resulting from the decrease of the inclination angle. It can be observed that the increase in the inclination angle results in a delay of the onset of absolute instability.

Table 1.

Absolute and convective instability regimes computed using various models when $\mathit{Re}_{0}=10$ , $\mathit{Ka}=1988.5$ and $N=0.4$ .

Table 2.

Comparison with experiments of Zapke & Kröger (Reference Zapke and Kröger2000a ) with air and methanol as working fluids (see their figure 8). Here, $V_{\ell }$ denotes the average liquid velocity, $V_{g,ZK}$ corresponds to the average flooding gas velocity in Zapke & Kröger (Reference Zapke and Kröger2000a ) and $V_{g,low}$ and $V_{g,up}$ correspond to the lower and upper limits of absolute instability respectively, computed using our analysis.

5.3.2. Comparison with experiments

In order to verify the predictions of our model, we consider air and methanol as working fluids in a vertical channel; these fluids were also considered in the experiments by Zapke & Kröger (Reference Zapke and Kröger2000a ) (see their figure 8). As the parameters fall in the regime where our WIBL model was shown to have a good agreement with the OS problem, we consider the WIBL model to compute the limits of absolute instability, which we believe has a close connection to flooding. The comparison shown in table 2 indicates a fair agreement. The experimental flooding point turns out to be close to the region of absolute instability (moreover, it turns out that the flooding point lies between the absolute instability boundaries, except for the last line in the table), but it has to be noted that the trend is slightly different, i.e. the gas velocities for the absolute instability boundaries monotonically increase as the liquid velocity increases, whereas in the experiments of Zapke & Kröger (Reference Zapke and Kröger2000a ) the flooding gas velocity is found to decrease with the liquid velocity (which appears to be counterintuitive). However, confidence limit error bars were not provided for the experiments of Zapke & Kröger (Reference Zapke and Kröger2000a ), and, therefore, it is not really clear how accurate the experimental findings are. Another explanation for the differences between our theoretical predictions and the experiments may be that our analysis assumes two-dimensional flow, whereas two-dimensionality may not be valid in the experiments. Moreover, the observed differences may be related to the way in which Zapke & Kröger (Reference Zapke and Kröger2000a ) define flooding and their interpretation of flooding, namely they state that flooding occurs when a vortex (large-amplitude wave) formed at the inlet of the channel has grown to such an extent that it obstructs the flow of the gas (i.e. the height of the wave is comparable with the height of the channel), resulting in the liquid travelling up past the injection point. (We note that Zapke & Kröger (Reference Zapke and Kröger2000a ) used a channel with a rectangular cross-section of width 10 mm and height 50 mm.) However, the width of the channel considered in our model is assumed to be very large compared with the liquid film thickness, so that the amplitude of the standing wave formed at the point of wave reversal is not comparable with the channel width. Notwithstanding the above mentioned disparities, our linear theory appears to be useful in predicting the nonlinear phenomenon of flooding, but clearly further detailed experiments are needed to test our theoretical predictions carefully.

Figure 17.

Time-dependent computation of the linearised WIBL model subject to a localised disturbance, showing convective instability for (a) ${\it\Theta}=1.0$ and (b) ${\it\Theta}=2.2$ , with $\mathit{Re}_{0}=10$ , ${\it\theta}={\rm\pi}/2$ , $\mathit{Ka}=1988.5$ and $N=0.4$ . The three-dimensional plot is inclined here for better visualisation.

5.4. Time-dependent computations

We further explore the absolute and convective instability regimes predicted by our linear analysis in § 5.3, by solving the linearised WIBL model (4.16) and (4.17) as an initial-value problem subject to a localised initial condition. This is carried out numerically using a scheme that is based on a Fourier pseudo-spectral representation of the derivatives along the spatial direction and an adaptive stepping along time, e.g. Tseluiko et al. (Reference Tseluiko, Saprykin, Duprat, Giorgiutti-Dauphiné and Kalliadasis2010). We choose the initial condition of the interface disturbance, $f(x,t)$ , to be localised about $x=0$ , namely, we choose the Gaussian function $f(x,0)=0.1\exp (-0.5x^{2})$ . In addition, we choose $g(x,0)=0$ . The liquid Reynolds number is taken to be $\mathit{Re}=5$ , consistent with that in the previous section. The equations are integrated on a long periodic domain, $[-750,750]$ , for a sufficiently long time, say $t\approx 500$ , to determine the evolution of the wavepackets.

We begin with a relatively small value of the gas shear, ${\it\Theta}=1.0$ , and the corresponding evolution of the disturbance is depicted as a space–time plot in figure 17(a). It can be observed that the wavepacket starts propagating downstream along the positive $x$ direction, therefore the instability is convective. On the other side of the gas flow spectrum, e.g. for ${\it\Theta}=2.2$ , the instability is again convective, with the wavepacket propagating upstream along the negative $x$ direction (figure 17 b). This can be further confirmed by examining the wave profile at a sufficiently long time instant and the time evolution of the disturbance, $f(x,t)$ , at $x=0$ , see the top and bottom panels of figure 18 respectively. From figure 18, it can be ascertained that the wave profiles for ${\it\Theta}=1.0$ and ${\it\Theta}=2.2$ are convected in the positive and negative $x$ directions (see the top panels), and also the disturbance, $f(0,t)$ eventually decays to zero (see the bottom panels). This is consistent with the analytical results discussed above. For intermediate values of the gas flow rate, e.g. for ${\it\Theta}=1.5$ and ${\it\Theta}=1.65$ shown in figure 19, the wavepackets invade both the negative and the positive $x$ directions, and for longer times they are expected to contaminate the whole domain. This can also be established from figure 20, from which it can be concluded that the disturbance $f(0,t)$ does not decay to zero. These results agree with the band of absolutely unstable gas flow rates, $1.2\leqslant {\it\Theta}\leqslant 2.1$ , predicted by our analysis of various models in § 5.3.

Figure 18.

Time-dependent computations of the linearised WIBL model depicting wave profiles at $t=500$ and the time evolution of the disturbance at $x=0$ for (a) ${\it\Theta}=1.0$ and (b)  ${\it\Theta}=2.2$ , with $\mathit{Re}_{0}=10$ , ${\it\theta}={\rm\pi}/2$ , $\mathit{Ka}=1988.5$ and $N=0.4$ .

Figure 19.

Time-dependent computation of the linearised WIBL model subject to a localised disturbance, showing convective instability for (a) ${\it\Theta}=1.5$ and (b) ${\it\Theta}=1.65$ , with $\mathit{Re}_{0}=10$ , ${\it\theta}={\rm\pi}/2$ , $\mathit{Ka}=1988.5$ and $N=0.4$ . The three-dimensional plot is inclined here for better visualisation.

Figure 20.

Time-dependent computations of the linearised WIBL model depicting wave profiles at $t=500$ and the time evolution of the disturbance at $x=0$ for (a) ${\it\Theta}=1.5$ and (b)  ${\it\Theta}=1.65$ , with $\mathit{Re}_{0}=10$ , ${\it\theta}={\rm\pi}/2$ , $\mathit{Ka}=1988.5$ and $N=0.4$ .