2.1 Nusselt flow

We consider a liquid thin film on an inclined plane. The angle of the plane with respect to the horizontal is $\unicode[STIX]{x1D703}$ . The flow is supposed to be two-dimensional. The directions $Ox$ and $Oz$ are respectively the directions parallel to the plane and normal to the plane. The components of the fluid velocity field are denoted by $u$ and $w$ in the directions $Ox$ and $Oz$ respectively. The bottom wall is plane and the fluid depth is denoted by $h$ . The notations are specified in figure 1.

Figure 1.

Notations used in the text.

The liquid film is thin. If $h_{0}$ is the order of magnitude of the fluid depth and $L$ the order of magnitude of the typical distance of variations of the variables of the flow in the $Ox$ direction, then we can define a small parameter $\unicode[STIX]{x1D700}$ by

The liquid is an incompressible Newtonian fluid of density $\unicode[STIX]{x1D70C}$ , dynamic viscosity $\unicode[STIX]{x1D707}$ and kinematic viscosity $\unicode[STIX]{x1D708}=\unicode[STIX]{x1D707}/\unicode[STIX]{x1D70C}$ . Denoting the pressure by $p$ , the Cauchy stress tensor $\unicode[STIX]{x1D748}$ can be written as $\unicode[STIX]{x1D748}=-p\unicode[STIX]{x1D644}+\unicode[STIX]{x1D749}$ where $\unicode[STIX]{x1D644}$ is the identity tensor and $\unicode[STIX]{x1D749}$ is the viscous stress tensor which is related to the strain rate tensor $\unicode[STIX]{x1D63F}$ by $\unicode[STIX]{x1D749}=2\unicode[STIX]{x1D707}\unicode[STIX]{x1D63F}$ .

The boundary conditions are the no-penetration condition and the no-slip condition on the bottom and, at the free surface, the kinematic boundary condition and the dynamic boundary condition

where $\unicode[STIX]{x1D6FE}$ is the surface-tension coefficient, $R$ the radius of curvature of the free surface and $\boldsymbol{n}$ the unit normal vector at the free surface that points towards the atmosphere. The atmospheric pressure is constant and can be taken equal to zero and does not appear in the equations. With the boundary conditions, the Navier–Stokes equations admits a stationary solution with a constant depth $h_{0}$ , a parabolic velocity profile

and a hydrostatic pressure distribution

This solution is the Nusselt flow. The velocity averaged over the fluid depth is

The following relation is satisfied for the Nusselt flow:

In order to study a change of Galilean reference frame, we can consider the more general case where the bottom wall has a velocity $v_{b}$ in the $Ox$ direction. The no-slip condition is then

The pressure distribution (2.4) is unchanged and the parabolic velocity profile becomes

The Galilean invariance is studied in § 3.5. Otherwise, we suppose $v_{b}=0$ .

In the following, we consider Nusselt flow as the reference flow which is perturbed when this flow becomes unstable. We perform an asymptotic analysis of this flow in the long-wave regime.

2.2 Dimensionless form of the equations

The equations are put into dimensionless form using the depth $h_{0}$ and velocity $\langle u_{0}\rangle$ of the Nusselt flow as the characteristic length in the $Oz$ direction and characteristic speed in the $Ox$ direction. The characteristic length in the $Ox$ direction is $L=H/\unicode[STIX]{x1D700}$ and the characteristic speed in the $Oz$ direction is $\unicode[STIX]{x1D700}\langle u_{0}\rangle$ . The following dimensionless numbers can be defined, the Reynolds number

the Froude number

and the Weber number

It follows from equation (2.5) that the dimensionless number defined by $\unicode[STIX]{x1D706}=(Re\sin \unicode[STIX]{x1D703})/F^{2}$ satisfies the relation

Since $\unicode[STIX]{x1D706}$ could be different from 3 if a different characteristic speed was used instead of (2.5) to put the equations into dimensionless form, we keep $\unicode[STIX]{x1D706}$ in the following equations choosing not to replace it by 3. The following dimensionless quantities, denoted by a tilde, are defined with the scaling of the Nusselt flow and the small parameter $\unicode[STIX]{x1D700}$ defined in (2.1) which gives the ratio of transverse to longitudinal spatial scales:

For convenience, we define the following dimensionless number

This number is related to the Eötvös number (also called the Bond number) defined as $Eo=Bo=\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}gL^{2}/\unicode[STIX]{x1D6FE}$ (with $\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}=\unicode[STIX]{x1D70C}_{liq}-\unicode[STIX]{x1D70C}_{gas}\simeq \unicode[STIX]{x1D70C}$ ) by $\unicode[STIX]{x1D705}=\unicode[STIX]{x1D700}^{2}/Eo$ . The equations are then put in dimensionless form

Equation (2.15) is the dimensionless continuity equation. Equations (2.16)–(2.17) are the Navier–Stokes momentum equations. Equation (2.18) gives the relation between the components of the viscous stress tensor and those of the strain rate tensor and equations (2.19)–(2.22) are the dimensionless boundary conditions.

2.3 Successive approximations

3.1 Mass conservation

The mass equation (2.15) is averaged over the depth. With the kinematic boundary condition (2.20), we get the exact equation

If the asymptotic solution of order 0 for $\tilde{U}$ is used, i.e. if $\tilde{U} =\langle \tilde{u} _{0}\rangle$ , then equation (3.2) together with the expression (2.31) of $\langle \tilde{u} _{0}\rangle$ gives the kinematic wave equation

which shows that a depth perturbation is propagated at the celerity $\tilde{c}_{0}=\unicode[STIX]{x1D706}\tilde{h}^{2}$ which is equal to $3\tilde{U}$ . This is the celerity of the kinematic waves.

3.2 Momentum balance

The momentum balance equation (2.16) is averaged over the depth. Using the no-penetration condition (2.19), the kinematic boundary condition (2.20) and the dynamic boundary conditions (2.21) and (2.22) written

the following equation is found

Equation (2.17) is used to express the integral of the pressure of (3.6). The left-hand side of (2.17) is of $O(\unicode[STIX]{x1D700}^{2})$ and hence it does not appear in the second-order model. The integration of (2.17) with the boundary conditions leads to

With this expression of the pressure, and since (2.18) implies that $\tilde{\unicode[STIX]{x1D70F}}_{zz}=-\tilde{\unicode[STIX]{x1D70F}}_{xx}$ , equation (3.6) becomes

The two terms with integrals can be calculated with the expressions of the components of the viscous stress tensor given by (2.18) with $\tilde{\unicode[STIX]{x1D70F}}_{xz}=\unicode[STIX]{x2202}\tilde{u} /\unicode[STIX]{x2202}\tilde{z}+O(\unicode[STIX]{x1D700}^{2})$ . The result depends on the value of $\tilde{u}$ at $\tilde{z}=\tilde{h}$ , denoted by $\tilde{u} (\tilde{h})$ :

Since the terms involving $\tilde{u} (\tilde{h})$ are all of second order, $\tilde{u} (\tilde{h})$ can be evaluated at order 0. As $\tilde{u} _{0}(\tilde{h})=\unicode[STIX]{x1D706}\tilde{h}^{2}/2$ and $\langle \tilde{u} _{0}\rangle =\unicode[STIX]{x1D706}\tilde{h}^{2}/3$ (see (2.28) and (2.31)),

Equation (3.9) becomes

The two terms of $O(\unicode[STIX]{x1D700}/Re)$ can be estimated at order 0:

Written this way, this term has the form of a viscous diffusive term. Equation (3.11) can therefore be written

To write this equation at the first order of approximation, the order of magnitude of $\unicode[STIX]{x1D705}$ has to be increased to $O(1)$ (i.e. $\unicode[STIX]{x1D700}\unicode[STIX]{x1D705}_{1}$ is replaced by $\unicode[STIX]{x1D705}$ ). The averaged momentum equation reduces to

The velocity $\tilde{u}$ is decomposed as the sum of its average value $\tilde{U}$ and a fluctuation $\tilde{u} ^{\prime }$ which depends on $\tilde{x}$ , $\tilde{z}$ and $\tilde{t}$ :

By definition, $\langle \tilde{u} ^{\prime }\rangle =0$ . The quantity $\langle \tilde{u} ^{2}\rangle$ can be written

At the first order of approximation, $\langle \tilde{u} ^{\prime 2}\rangle$ could be estimated at order 0 either as a function of $\tilde{h}$

or as a function of $\tilde{U}$

In the latter case, the first-order equation (3.14) becomes

If now $\tilde{\unicode[STIX]{x1D70F}}_{xz}(0)$ is supposed to be equal to $3\tilde{U} /\tilde{h}$ , which is true at order 0 but wrong at order 1 (see (2.35)), equation (3.19) gives the Shkadov equation. The first-order correction to $\tilde{\unicode[STIX]{x1D70F}}_{xz}(0)$ cannot be neglected, even at the first order of approximation. Using the asymptotic expansion $\tilde{\unicode[STIX]{x1D70F}}_{xz}=\tilde{\unicode[STIX]{x1D70F}}_{0}+\unicode[STIX]{x1D700}\tilde{\unicode[STIX]{x1D70F}}_{1}+\unicode[STIX]{x1D700}^{2}\tilde{\unicode[STIX]{x1D70F}}_{2}+\cdots$ and the expression $\tilde{\unicode[STIX]{x1D70F}}_{0}(0)=\unicode[STIX]{x1D706}\tilde{h}$ coming from (2.27), a consistent form of equation (3.19) is

At the first order,

Since $\tilde{\unicode[STIX]{x1D70F}}_{1}(0)\neq 3\langle \tilde{u} _{1}\rangle /\tilde{h}$ ,

Therefore, the Shkadov equation is not consistent. However it is possible to derive a consistent two-equation model. For instance, the asymptotic expansions (2.33) and (2.34) allow us to write

Using this expression (3.23) and the expression (3.17), equation (3.14) becomes

which, together with equation (3.2), gives the model of Noble & Vila (Reference Noble and Vila2013) in the particular case of Newtonian fluids. Another consistent way to write this is

Using this expression and (3.17), equation (3.14) becomes

It is now possible to derive a consistent two-equation model including diffusion by using (3.13). We can write

The expression (2.39) can be written

Equation (3.13) becomes

Keeping only viscous diffusive and capillary dispersive terms, the expression (2.40) can be written

With (3.25), (3.27) and (3.30), equation (3.29) becomes

Equation (3.31) and (3.2) together constitute another two-equation model.

Another approach, first proposed by Teshukov (Reference Teshukov2007) and expanded by Richard & Gavrilyuk (Reference Richard and Gavrilyuk2012), is to treat $\langle \tilde{u} ^{\prime 2}\rangle$ as an independent variable rather than enslaving it to $\tilde{h}$ as in (3.17) or to $\tilde{U}$ as in (3.18). The introduction of this new variable guarantees the compatibility of the work–energy theorem with the mass and momentum equations (see § 3.3). Since $\langle \tilde{u} ^{\prime 2}\rangle$ is the variance of $\tilde{u}$ , the three variables of the flow are its depth, its average velocity and the variance of its velocity. In fact, it is more convenient to use

as the third variable of the model. In the particular case of a constant vorticity $\unicode[STIX]{x1D714}$ , $\unicode[STIX]{x1D711}$ is proportional to the square of the vorticity, more precisely $\unicode[STIX]{x1D711}=\unicode[STIX]{x1D714}^{2}/12$ (Teshukov Reference Teshukov2007). For this reason, $\unicode[STIX]{x1D711}$ was called enstrophy (Richard & Gavrilyuk Reference Richard and Gavrilyuk2012) because every definitions of an enstrophy variously relate it to the square of the vorticity. The dimensionless form of $\unicode[STIX]{x1D711}$ is given by $\tilde{\unicode[STIX]{x1D711}}=\unicode[STIX]{x1D711}h_{0}^{2}/\langle u_{0}\rangle ^{2}$ . With the introduction of the enstrophy, equation (3.13) can be written

The expansion of the enstrophy

leads to $\tilde{\unicode[STIX]{x1D711}}_{0}=\langle \tilde{u} _{0}^{\prime 2}\rangle /\tilde{h}^{2}$ , $\tilde{\unicode[STIX]{x1D711}}_{1}=2\langle \tilde{u} _{0}^{\prime }\tilde{u} _{1}^{\prime }\rangle /\tilde{h}^{2}$ and $\tilde{\unicode[STIX]{x1D711}}_{2}=(\langle \tilde{u} _{1}^{\prime 2}\rangle +2\langle \tilde{u} _{0}^{\prime }{\tilde{u} ^{\prime }}_{2}\rangle )/\tilde{h}^{2}$ . The calculation at order 0 gives

For the expansion up to first order, with $\unicode[STIX]{x1D705}=O(1)$ ,

For the expansion up to second order, with $\unicode[STIX]{x1D705}=\unicode[STIX]{x1D700}Re\unicode[STIX]{x1D705}_{1}$ ,

Keeping only viscous diffusive and capillary dispersive terms, the second-order term of the enstrophy expansion reduces to

The full expression of the second-order term of the enstrophy is given in appendix B.

With the introduction of a third variable, a third equation is needed. This equation comes from an energy equation which completes the mass conservation equation and the momentum balance equation.

3.3 Kinetic energy equation

The basis of the energy equation is the work–energy theorem (also called theorem of the kinetic energy). The differential form of this theorem is

where $\boldsymbol{V}$ is the velocity field. The colon symbol denotes the double contraction. Using the explicit expressions of the components of the vectors and tensors involved in (3.39), we get the equation

Its dimensionless form is

We calculate the average of this equation over the depth. The details of this long calculation are given in appendix C. Taking into account the boundary conditions and dropping all second-order terms except viscous diffusive and dispersive capillary terms, the averaged work–energy theorem can finally be written

At the first order of approximation, increasing $\unicode[STIX]{x1D705}$ to $O(1)$ , the averaged work–energy theorem reduces to

It can be seen that there is no first-order corrective term to the right-hand side of this equation. The right-hand side can be interpreted as the sum of the power of the component of the weight on the $Ox$ axis and of the power of a friction force $-3\tilde{U} /\tilde{h}$ . This is exactly the same result as that found at order 0. We recover the result that was pointed out by Luchini & Charru (Reference Luchini and Charru2010): the first-order correction to the dissipation function is identically zero. Because of this, the averaged work–energy theorem is consistent to the first order in its straightforward formulation (3.43) whereas the momentum balance equation (3.14) needs a first-order correction to the shear stress at the bottom $\tilde{\unicode[STIX]{x1D70F}}_{xz}(0)$ because of (2.35).

Using the expansion (3.15), we get $\langle \tilde{u} ^{2}\rangle =\tilde{U} ^{2}+\langle \tilde{u} ^{\prime 2}\rangle$ and

In the left-hand side of equation (3.43), at the first order of approximation, $\langle \tilde{u} ^{\prime 2}\rangle$ and $\langle \tilde{u} ^{\prime 3}\rangle$ can be estimated at order 0, either as a function of $\tilde{h}$ , (3.17) for $\langle \tilde{u} ^{\prime 2}\rangle$ and

for $\langle \tilde{u} ^{\prime 3}\rangle$ , or as a function of $\tilde{U}$ , (3.18) and

For example, enslaving these quantities to $\tilde{U}$ , equation (3.43) becomes (Luchini & Charru Reference Luchini and Charru2010):

Even if this equation is written differently, it is the same as the first-order model derived by Ruyer-Quil & Manneville (Reference Ruyer-Quil and Manneville2000) (but formulated in terms of $h$ and $q=hU$ ) and also as the energy equation derived by Usha & Uma (Reference Usha and Uma2004) by the energy integral method if the second-order terms are discarded (equation (31) of Usha & Uma (Reference Usha and Uma2004) restricted to the first order of approximation). The terms of $O(\unicode[STIX]{x1D700}/Re)$ of equation (C 23) are consistent to the corresponding terms of the second-order two-equation RQM model. The model of Usha & Uma (Reference Usha and Uma2004) is identical to the two-equation RQM model except that it adds third-order terms. This is also the case of the model of Abderrahmane & Vatistas (Reference Abderrahmane and Vatistas2007) who proposed an improvement of the third-order terms of Usha & Uma (Reference Usha and Uma2004). The model of Mudunuri & Balakotaiah (Reference Mudunuri and Balakotaiah2006) is also identical to the two-equation RQM model at the first order of approximation and is consistent to it at the second order of approximation (terms of $O(\unicode[STIX]{x1D700}/Re)$ ). In all these models, the second-order inertial terms are also neglected.

3.4 Structure of the Euler equations of compressible fluids

Using the enstrophy as a third variable, we can write $\langle \tilde{u} ^{3}\rangle =\tilde{U} ^{3}+3\tilde{h}^{2}\tilde{U} \tilde{\unicode[STIX]{x1D711}}+\langle \tilde{u} ^{\prime 3}\rangle$ . The approach in Richard & Gavrilyuk (Reference Richard and Gavrilyuk2012) was similar except that the flow was supposed to be weakly sheared. This hypothesis allowed to neglect the term with $\langle \tilde{u} ^{\prime 3}\rangle$ in the energy equation. This was sufficient to close the problem. Here, because of the no-slip condition, the flow is not weakly sheared. It is not possible to treat $\langle \tilde{u} ^{\prime 3}\rangle$ as a fourth variable because this would lead to an infinite cascade of equations and the problem would not be closed. On the other hand, we have a reference flow, which is the Nusselt flow. To close the problem, we can calculate $\langle \tilde{u} ^{\prime 3}\rangle$ as a function of $\tilde{h}$ with the asymptotic expansion. Equation (3.42) becomes

At the first order of approximation, $\langle \tilde{u} ^{\prime 3}\rangle$ can be estimated at order 0.

Under the form of (3.2), (3.33) and (3.48), the equations have not yet a satisfactory structure. Denoting

and

the equations (3.2), (3.33) and (3.48) of the system can be written

3.5 Galilean invariance

The general case of a bottom wall moving in the $Ox$ direction at a constant velocity $v_{b}$ allows us to study the Galilean invariance of the model equations. When $U$ is the fluid average velocity in the Galilean reference frame $\mathscr{R}$ for which the bottom velocity is $v_{b}$ , $U-v_{b}$ is the fluid average velocity in the Galilean reference frame $\mathscr{R}_{0}$ of the bottom wall. In $\mathscr{R}$ , the average velocity of the Nusselt flow (2.5) takes the more general form

We can define  $\tilde{v}_{b}=v_{b}/\langle u_{0}^{\mathscr{R}0}\rangle$ where $\langle u_{0}^{\mathscr{R}0}\rangle$ is the average velocity of the Nusselt flow in $\mathscr{R}_{0}$ . The velocity $\langle u_{0}^{\mathscr{R}0}\rangle$ is used to define the dimensionless numbers (2.9)–(2.11) and to put the equations into dimensionless form. At order 0, $\langle \tilde{u} _{0}\rangle =\unicode[STIX]{x1D706}\tilde{h}^{2}/3+\tilde{v}_{b}$ . The mass conservation equation takes the same form in every Galilean reference frame. The momentum balance equation in $\mathscr{R}$ can be written

where $\tilde{f}_{2}=315(\tilde{U} -\tilde{v}_{b})/(2\unicode[STIX]{x1D706}^{2}\tilde{h}^{3})+O(\unicode[STIX]{x1D700})$ . This equation is Galilean invariant. The case of the energy equation is more complex since the kinetic energy is not Galilean invariant. However, supplemented by the momentum equation, the work–energy theorem as a whole is Galilean invariant (see also the discussion in Ostapenko (Reference Ostapenko2014) in the case of the one-dimensional equations of Teshukov (Reference Teshukov2007)) and the energy equation can be put in the invariant form

where $\tilde{e}_{i}=\tilde{h}^{2}\tilde{\unicode[STIX]{x1D711}}/2+\tilde{h}\cos \unicode[STIX]{x1D703}/(2F^{2})$ , $\tilde{P}$ , $f_{1}$ , $\tilde{f}_{2}$ and $\tilde{\unicode[STIX]{x1D711}}$ are Galilean invariant. It can be noted that, for this system, neither the momentum flux $\tilde{P}$ nor the potential energy $\tilde{e}_{i}$ (which is analogous to the internal energy of compressible fluids) depends on any velocity. For two-equation models that use $\langle u^{2}\rangle =6U^{2}/5$ in $\mathscr{R}_{0}$ , both the momentum flux and the potential energy depend on the difference between the average velocity of the fluid and the velocity of the bottom wall since, in this case, in $\mathscr{R}$ , $\langle u^{2}\rangle =U^{2}+(U-v_{b})^{2}/5$ . The first term is not Galilean invariant whereas the second term is. From the discussion in Gavrilyuk & Perepechko (Reference Gavrilyuk and Perepechko1998), this means that the first term gives the kinetic energy $U^{2}/2$ and the second term is included in the potential energy as $(U-v_{b})^{2}/10$ .

4.1 Long-wave instability threshold

In order to study the neutral stability, the equations (3.2), (3.65) and (3.69) are linearized around the parameters of the Nusselt flow which are $\tilde{h}_{0}=1$ , $\tilde{U} _{0}=1$ and $\tilde{\unicode[STIX]{x1D711}}_{0}=1/5$ with $\unicode[STIX]{x1D706}=3$ . We write $\tilde{h}=\tilde{h}_{0}+h^{\prime }$ , $\tilde{U} =\tilde{U} _{0}+U^{\prime }$ and $\tilde{\unicode[STIX]{x1D711}}=\tilde{\unicode[STIX]{x1D711}}_{0}+\unicode[STIX]{x1D711}^{\prime }$ where $h^{\prime }$ , $U^{\prime }$ and $\unicode[STIX]{x1D711}^{\prime }$ are infinitesimal perturbations of the Nusselt flow. Taking perturbations of the form $[h^{\prime }\quad U^{\prime }\quad \unicode[STIX]{x1D711}^{\prime }]^{\text{T}}=[A_{1}\quad A_{2}\quad A_{3}]^{\text{T}}\exp [\text{i}k(x-ct)]$ , the three linearized equations can be written

Since the equations are linearized around the parameters of the Nusselt flow, the relaxation terms are of order 1 with respect to the perturbations (they are automatically equal to 0 at order zero). The functions $f_{1}$ and $\tilde{f}_{2}$ appear only in the relaxation terms. Consequently the linearized equations do not depend on the particular form chosen for $f_{1}$ and $\tilde{f}_{2}$ since they differ from each other by a term which is of order 1 with respect to the perturbations. This means that the corrective term is nonlinear. All consistent expressions of $f_{1}$ and $\tilde{f}_{2}$ give $f_{1}(\tilde{h}=1,\tilde{U} =1,\tilde{\unicode[STIX]{x1D711}}=1/5)=7/3$ and $\tilde{f}_{2}(\tilde{h}=1,\tilde{U} =1,\tilde{\unicode[STIX]{x1D711}}=1/5)=35/2$ . To obtain solutions different from the trivial solution, the matrix of this system must be singular. The dispersive relation is obtained by writing that the determinant of this system is equal to zero. The solution can be written

The Nusselt flow is stable if $\text{Im}(c)<0$ . At the long-wave limit ( $k\rightarrow 0$ ), the stability occurs if $F<F_{c}$ where $F_{c}$ is the critical Froude number

This long-wave stability threshold can be written with the Reynolds number. Since $Re=3F^{2}/\sin \unicode[STIX]{x1D703}$ (see (2.12)), the Nusselt flow is stable at the long-wave limit if

and it is unstable if $Re>(5\cot \unicode[STIX]{x1D703})/6$ . This is exactly the instability threshold found by Benjamin (Reference Benjamin1957) and Yih (Reference Yih1963) with the Orr–Sommerfeld equation.

Figure 2.

Neutral stability curves (wavenumber $k$ as a function of the Reynolds number $Re$ ) calculated from the present three-equation model (thick solid curve), the Orr–Sommerfeld equation (thin solid curve), the two-equation model of Ruyer-Quil & Manneville (Reference Ruyer-Quil and Manneville2000) (dotted curve) and the two-equation model of Usha & Uma (Reference Usha and Uma2004) based on energy integral method (dashed curve). (a) $\unicode[STIX]{x1D703}=4.6^{\circ }$ , $Ka=332$ ; (b) $\unicode[STIX]{x1D703}=6.4^{\circ }$ , $Ka=238$ . (c) $\unicode[STIX]{x1D703}=30.0^{\circ }$ , $Ka=238$ ; (d) $\unicode[STIX]{x1D703}=90.0^{\circ }$ , $Ka=238$ .

4.2 Neutral stability curves

The neutral stability curves are then plotted for various inclination angles $\unicode[STIX]{x1D703}$ and for some values of the Kapitza number defined by

These curves give the marginally stable states in the $(k,Re)$ plane, dividing a stable region above the curve and an unstable region below the curve. The marginally stable curve corresponds to $\text{Im}(c)=0$ .

Due to the first-order consistency of the equations, the neutral stability curve is nearly identical to the result of the Orr–Sommerfeld equation in the validity domain of the asymptotic expansion which is $kRe\ll 1$ since $k\sim \unicode[STIX]{x1D700}$ . Outside this domain, the neutral stability curve given by the model is slightly below the Orr–Sommerfeld curve but is comparable to other models such as the two-equation model of Ruyer-Quil & Manneville (Reference Ruyer-Quil and Manneville2000) (see figure 2).

4.3 Growth rate

The spatial stability problem is studied with the linearized system taking infinitesimal perturbations of the form $[h^{\prime }\quad U^{\prime }\quad \unicode[STIX]{x1D711}^{\prime }]^{\text{T}}=[A_{1}~A_{2}~A_{3}]^{\text{T}}\exp [\text{i}(kx-\unicode[STIX]{x1D714}t)]$ with $\unicode[STIX]{x1D714}$ real and $k$ complex. Writing $k=k_{r}+\text{i}k_{i}$ , we define a Reynolds number $R_{S}=3Re/2$ from the surface velocity instead of the mean velocity and a pulsation $\unicode[STIX]{x1D714}^{\ast }=2\unicode[STIX]{x1D714}/3$ to compare with Brevdo et al. (Reference Brevdo, Laure, Dias and Bridges1999) and Ruyer-Quil & Manneville (Reference Ruyer-Quil and Manneville2002). The spatial growth rate is $-k_{i}$ and is represented in figure 3 as a function of $\unicode[STIX]{x1D714}^{\ast }$ for various Reynolds numbers in the conditions of the experiments of Liu, Paul & Gollub (Reference Liu, Paul and Gollub1993): inclination angle $\unicode[STIX]{x1D703}=4.6^{\circ }$ ; kinematic viscosity $\unicode[STIX]{x1D708}=5.02\times 10^{-6}~\text{m}^{2}~\text{s}^{-1}$ ; surface tension $\unicode[STIX]{x1D6FE}=69.0\times 10^{-3}\;\text{N}~\text{m}^{-1}$ ; density $\unicode[STIX]{x1D70C}=1130\;\text{kg}~\text{m}^{-3}$ . The predicted growth rates are in agreement to those of Brevdo et al. (Reference Brevdo, Laure, Dias and Bridges1999) except at large Reynolds numbers where they are too small.

Figure 3.

Growth rate $-k_{i}$ as functions of pulsation $\unicode[STIX]{x1D714}^{\ast }$ for $R_{S}$ equal to the critical Reynolds number $R_{Sc}=15.5$ (curve 1), $R_{S}=10$ (dotted curve), $R_{S}=20$ (curve 2), $R_{S}=40$ (curve 3), $R_{S}=60$ (curve 4), $R_{S}=100$ (curve 5) and $R_{S}=200$ (curve 6), in the conditions of the experiments of Liu et al. (Reference Liu, Paul and Gollub1993) (see text). (a) Brevdo et al. (Reference Brevdo, Laure, Dias and Bridges1999); (b) calculated with the model.

5.1 Numerical scheme

The results of the previous section show that the linear properties of the model are very similar to those of the two-equation RQM model and other models based on the energy integral method (Usha & Uma Reference Usha and Uma2004; Abderrahmane & Vatistas Reference Abderrahmane and Vatistas2007). Now the simpler structure of our three-equation model allows an easier and more reliable numerical resolution than those two-equation models.

In particular, since these two-equation models feature a mass equation and an energy equation but are written and solved in terms of depth $h$ and average velocity $U$ (or discharge $q$ ) with a non-conservative momentum equation, the structure is not simple and, in particular, it is not obvious whether they admit a mathematical entropy (in classical two-equation models with equations of mass and momentum in conservative form, the mathematical entropy is the energy). On top of that, the third-order terms of Usha & Uma (Reference Usha and Uma2004) and Abderrahmane & Vatistas (Reference Abderrahmane and Vatistas2007) are very complex and are certainly extremely difficult to solve numerically.

On the other hand, the three-equation model has the same mathematical structure as the Euler equations of compressible fluids. The mass, momentum and energy equation are all in conservative form and the enstrophy plays the role of a mathematical entropy (Richard & Gavrilyuk Reference Richard and Gavrilyuk2012). This allows to use classical and simple numerical schemes which are reliable and robust. The capillary dispersive terms need a special treatment however. Without this treatment, upwind schemes are linearly unstable, whether they are explicit or implicit, and centred schemes are linearly stable if they are implicit (explicit schemes are linearly unstable) but nothing can be proved for the nonlinear stability (Noble & Vila Reference Noble and Vila2014). We used the numerical scheme of Noble & Vila (Reference Noble and Vila2014). With this treatment of the capillary term, it can be proved, in the case of a Rusanov Riemann solver, that the scheme is nonlinearly stable. An explicit scheme is nonlinearly stable under a Courant–Friedrichs–Levy (CFL) condition (see Noble & Vila (Reference Noble and Vila2014)) and a semi-implicit scheme is unconditionally nonlinearly stable (see Bresch et al. Reference Bresch, Couderc, Noble and Vila2016). The particular mathematical structure of the equations of our model gives the guarantee of the nonlinear stability. For the structure of most other models, no proof of nonlinear stability can be found. The equations are now written in dimensional form to make a comparison with experimental results. Using the variable $W$ (see (3.70) and (3.72)), the equations can be written, with ${\hat{g}}=g\sin \unicode[STIX]{x1D703}$

This is a system of four equations for the four variables $h$ , $U$ , $\unicode[STIX]{x1D711}$ and $W$ . The fourth equation (5.4) is deduced from the first equation (5.1). It is a transport equation with a second-order derivative term. The energy of this system takes the special form

where there is a new term due to capillary effects. The new variable $W$ allows a reduction of the order of the system. There are no more third derivative terms. All capillary terms take the form of second-order terms which can be handled numerically much more easily. The numerical scheme we used and the stability condition can be found in Noble & Vila (Reference Noble and Vila2014). Details are provided in appendix D.

The boundary conditions were chosen to simulate the experiments of Liu & Gollub (Reference Liu and Gollub1994). The film is 2.0 m long. A small sinusoidal disturbance of various frequency was applied at the entrance of the system. Either the depth was perturbed with a constant flow rate or the flow rate was perturbed. The final waves did not depend on this choice. At the entrance, the enstrophy was supposed to be equal to its equilibrium value and $W$ was taken equal to 0.

At the end of the channel, we used a generalization of the boundary conditions used in Richard & Gavrilyuk (Reference Richard and Gavrilyuk2013). Denoting by $N$ the total number of cells, there was a special cell $N+1$ . The flow rate at the cell $N+1$ was calculated from the variables at the cell $N$ by a formula obtained by a dimensional analysis. Two conditions were used. The first one supposed that the outflow discharge depends only on $h$ . Consequently, the boundary condition is $q_{N+1}=C\sqrt{gh_{N}^{3}}$ with $C=q_{0}/\sqrt{gh_{0}^{3}}$ where $q_{0}$ and $h_{0}$ are the values of $q$ and $h$ of the corresponding Nusselt flow. The second outflow boundary condition supposed that the outflow discharge depends both on the depth and on the enstrophy. This condition is written $q_{N+1}=C^{\prime }\,h_{N}^{2}\sqrt{\unicode[STIX]{x1D711}_{N}}$ with $C^{\prime }=q_{0}/(h_{0}^{2}\sqrt{\unicode[STIX]{x1D711}_{0}})$ where $q_{0}$ , $h_{0}$ and $\unicode[STIX]{x1D711}_{0}$ are the values of $q$ , $h$ and $\unicode[STIX]{x1D711}$ of the corresponding Nusselt flow. This second boundary condition proved to work even better than the first one and it was our final choice. No disturbance was noticed at the end of the system on the fluid depth $h$ and a very minor one was sometimes seen on the enstrophy and even then only in the very last cells. We supposed that the enstrophy and $W$ were simply advected at the end of the system, i.e. that $\unicode[STIX]{x1D711}_{N+1}=\unicode[STIX]{x1D711}_{N}$ and $W_{N+1}=W_{N}$ . We found that the simple condition $h_{N+1}=h_{N}$ was sufficient and that there is no need to write a more complex boundary condition for $h$ .

The initial condition was the Nusselt flow over the whole length of the system. However, this initial condition is not physical since the Nusselt flow is unstable in the conditions where waves appear. This non-physical initial condition leads to a non-physical first wave which can grow to huge and aberrant size. Physically, a Nusselt flow cannot exist indefinitely over an arbitrary long distance since spontaneous waves appear. This is the cause of the aberrant growth of the first wave. To prevent this problem, the initial condition was taken for a flow slightly above the instability threshold. In this case, the instability is small enough that a Nusselt flow over the whole system is an acceptable initial condition. When the first wave had gone out of the system, the discharge was progressively increased to reach the desired value.

Mesh convergence was reached. For the simulation of the experiments of Liu & Gollub (Reference Liu and Gollub1994), the mesh convergence was obtained with 8000 cells for a total length of 2.0 m. The numerical solution converges exactly towards the periodic solution calculated in the stationary case by the resolution of the corresponding system of ordinary differential equations.

5.2 Simulation of the Liu and Gollub experiments

The experiments of Liu & Gollub (Reference Liu and Gollub1994) were numerically simulated for an inclination angle of $\unicode[STIX]{x1D703}=6.4^{\circ }$ , $Re=19.33$ , a kinematic viscosity $\unicode[STIX]{x1D708}=6.28\times 10^{-6}~\text{m}^{2}~\text{s}^{-1}$ , a surface tension $\unicode[STIX]{x1D6FE}=67.0\times 10^{-3}~\text{N}~\text{m}^{-1}$ and a density $\unicode[STIX]{x1D70C}=1070~\text{kg}~\text{m}^{-3}$ .

Various forcing frequencies were used. The whole depth profile over a 2.0 m long film for a frequency of $1.5~\text{Hz}$ is shown in figure 4. After a transient part where the instability grows, a periodic solution is reached. The length of the transient part depends on how the inflow perturbation is done and in particular it depends on the amplitude of the sinusoidal perturbation. The final periodic solution depends only on the frequency of the inflow perturbation and not of its amplitude.

Figure 4.

Dimensionless depth profile of waves for a $2.0$ m long film forced at a frequency of $1.5\;\text{Hz}$ calculated with the model in the conditions of the Liu & Gollub (Reference Liu and Gollub1994) experiments (see text).

The final periodic solution obtained with a frequency of $1.5\;\text{Hz}$ is compared with the measurements of Liu & Gollub (Reference Liu and Gollub1994) in figure 5. There is a good agreement between the numerical solution and the experiments. The calculated wavelength is slightly longer than the measured one. As both frequencies were the same, this means that the numerical wave velocity is slightly greater than the measured velocity. The calculated amplitude of the wave is almost identical to the measured amplitude. The numerical curve is slightly shifted towards higher depths. The amplitudes of the capillary ripples are very close to the measured ones.

Figure 5.

Dotted curve: experimental dimensionless depth profile of waves measured by Liu & Gollub (Reference Liu and Gollub1994) with a forcing frequency of $1.5\;\text{Hz}$ . Solid curve: depth profile of waves calculated with the model in the same conditions (the $x$ -axis origin is arbitrary).

The calculated depth profiles for frequencies of $3.0~\text{Hz}$ , $4.5~\text{Hz}$ and $7.0~\text{Hz}$ are compared in figure 6 to the measured depth profiles of Liu & Gollub (Reference Liu and Gollub1994) of the same forcing frequencies. On the whole, the agreement is good. Small differences can be noticed but the overall profiles are the same. It should be noted that, for these frequencies, neither the experimental waves nor the numerically calculated waves are truly periodic in space. These waves are no longer solitary and interact significantly.

Figure 6.

Dimensionless depth profile of waves forced at various frequencies. In all cases, the abscissa is $x$ in cm. (a,c,e) Measurements of Liu & Gollub (Reference Liu and Gollub1994). Reprinted from Liu & Gollub (Reference Liu and Gollub1994), Solitary wave dynamics of film flows, Phys. Fluids, volume 6, issue 5, with the permission of AIP Publishing. (b,d,f) Numerical calculations from the three-equation model. Forcing frequency: $3.0\;\text{Hz}$ (a,b); $4.5\;\text{Hz}$ (c,d); $7.0\;\text{Hz}$ (e,f).

At a frequency of $7\;\text{Hz}$ , the obtained saturated waves are very similar to the measured ones. At a frequency of $4.5\;\text{Hz}$ , the first part of the wavy film is shown in figure 7. The formation and propagation of a small depression (figure 7) and then of a subsidiary peak behind the primary peak (figure 6) is very similar to the experimental observations with a wave amplitude slightly smaller.

Figure 7.

Dimensionless depth profile of waves forced at a frequency of $4.5\;\text{Hz}$ in the upstream part of the film. The abscissa is $x$ in cm. (a) Measurements of Liu & Gollub (Reference Liu and Gollub1994). Reprinted from Liu & Gollub (Reference Liu and Gollub1994), Solitary wave dynamics of film flows, Phys. Fluids, volume 6, issue 5, with the permission of AIP Publishing. (b) Numerical calculations from the three-equation model (abscissas are somewhat arbitrary since the growing of the waves depends on the amplitude of the forcing).

5.3 Physical relevance of the enstrophy variations

It remains to be seen if the enstrophy variations predicted by the model have any physical relevance. If the equilibrium parabolic velocity profile was reached everywhere, the expression of the enstrophy would be $g^{2}h^{2}\sin ^{2}\unicode[STIX]{x1D703}/(45\unicode[STIX]{x1D708}^{2})$ . It is interesting to study the deviation of the enstrophy from this value.

Nakoryakov et al. (Reference Nakoryakov, Pokusaev, Alekseenko and Orlov1977) and Alekseenko et al. (Reference Alekseenko, Nakoryakov and Pokusaev1985) measured the velocity profiles in a wavy liquid film on a vertical wall for $Re=12.4$ and $\unicode[STIX]{x1D708}=7.2\times 10^{-6}\;\text{m}^{2}~\text{s}^{-1}$ , forced at a particular frequency. For a wavy film with an averaged depth of $0.545\;\text{mm}$ and a wavelength of $36\;\text{mm}$ , the authors distinguished four regions in the wave. In region 1, the ‘velocity profile is described by the self-similar parabolic law’. In region 2, ‘the velocity profile is less filled as compared to the parabolic one’ while in region 3, it is more filled. In region 4, no velocity profile was found by the authors due to a great scatter of the experimental points. From the front wave to the back of the wave, there is first a region 4 then a region 3 near the maximum depth followed by a region 1, then by a region 2 and again a region 1.

In figure 8 the enstrophy profile, calculated from the model in the conditions of the experiments of Nakoryakov et al. (Reference Nakoryakov, Pokusaev, Alekseenko and Orlov1977) and Alekseenko et al. (Reference Alekseenko, Nakoryakov and Pokusaev1985), is represented along with the equilibrium enstrophy at the same depth $g^{2}h^{2}\sin ^{2}\unicode[STIX]{x1D703}/(45\unicode[STIX]{x1D708}^{2})$ in dimensionless form. From the back of the wave towards its front, it is possible to distinguish (i) a region where the enstrophy $\unicode[STIX]{x1D711}$ is almost equal to its equilibrium value at the same depth $\unicode[STIX]{x1D711}_{eq}$ ; (ii) a region where $\unicode[STIX]{x1D711}$ is greater; an intersection of the curves which corresponds to a similar situation as in (i); (iii) a region where $\unicode[STIX]{x1D711}$ is smaller than $\unicode[STIX]{x1D711}_{eq}$ ; (iv) a region where $\unicode[STIX]{x1D711}$ is again greater. These regions coincide with the regions of Nakoryakov et al. (Reference Nakoryakov, Pokusaev, Alekseenko and Orlov1977) and they are denoted in figure 8 by the numbers 1, 2, 3 and 4. Note that, at this point, region (iv) cannot be clearly distinguished from region (ii).

Figure 8.

Numerical calculations with the three-equation model in the conditions of Nakoryakov et al. (Reference Nakoryakov, Pokusaev, Alekseenko and Orlov1977) ( $Re=12.4$ , $\unicode[STIX]{x1D703}=90^{\circ }$ , $Ka=193$ ); (a): solid curve: dimensionless enstrophy; dashed curve: dimensionless equilibrium enstrophy at the same depth; (b): solid curve dimensionless profile. Regions denoted by 1, 2, 3 or 4: see text.

5.4 Noise-driven waves

Instead of a periodic sinusoidal perturbation, a random noise can be applied at the inlet. The same procedure is used as in § 5.2 with the same numerical values ( $\unicode[STIX]{x1D703}$ , $Re$ , $\unicode[STIX]{x1D708}$ , $\unicode[STIX]{x1D6FE}$ and $\unicode[STIX]{x1D70C}$ ). The perturbation is chosen in a similar way as in Chang, Demekhin & Kalaidin (Reference Chang, Demekhin and Kalaidin1996): The depth at inlet is written $h(x=0,t)=h_{N}[1+f(t)]$ with

where the phase $\unicode[STIX]{x1D711}_{n}$ is a random number taken in the interval $[0,2\unicode[STIX]{x03C0}[$ and $\unicode[STIX]{x1D714}_{c}$ is a cutoff angular frequency. The amplitude was taken as $A=0.001$ and the cutoff frequency was $100\;\text{Hz}$ . As in Chang et al. (Reference Chang, Demekhin and Kalaidin1996), the number of terms was $N=1000$ . Therefore, the perturbation had 1000 sinusoidal terms with a random phase and every frequencies between $0.1~\text{Hz}$ and $100~\text{Hz}$ with a $0.1~\text{Hz}$ increment.

Figure 9.

Depth profile (a) and enstrophy profile (b) of a 4 m long thin film perturbed by a random noise. The dashed curve in (b) is the equilibrium enstrophy at the same depth. $\unicode[STIX]{x1D703}=6.4^{\circ }$ , $Re=19.33$ ; $\unicode[STIX]{x1D708}=6.28\times 10^{-6}~\text{m}^{2}~\text{s}^{-1}$ ; $\unicode[STIX]{x1D6FE}=67.0\times 10^{-3}\;\text{N}~\text{m}^{-1}$ ; $\unicode[STIX]{x1D70C}=1070~\text{kg}~\text{m}^{-3}$ .

Numerical calculations were made on a $4~\text{m}$ long channel. The results for the depth and the enstrophy are presented in figure 9 at a time $t=50\;\text{s}$ (the initial conditions and the procedure are the same as in § 5.2). The depth variations are qualitatively similar to previous studies such as Chang et al. (Reference Chang, Demekhin and Kalaidin1996). Distances between waves tend to increase as larger waves, which are slower to develop but which propagate faster than smaller waves, tend to catch smaller waves up and to swallow them, producing coalescence.

Except when a coalescence event is about to occur, the depth profile of a wave is similar to profiles found with a periodic perturbation. This is also the case with the enstrophy profile. Details of the enstrophy profile is given in figure 10: the beginning of the channel in (a), a part in the middle of the channel in (b) and near the end in (d). Figure 10(c) shows a coalescence event which is close to its conclusion as a big wave is about to swallow a smaller one (these two waves can also be seen in the depth profile of figure 9(a) around $x=3.2~\text{m}$ ). The equilibrium enstrophy at the same depth is also represented in dashed curves in figure 9(b) and in figure 10. The difference between the actual enstrophy and its equilibrium value is large in the wavefront region and rather small elsewhere.

Figure 10.

Enstrophy profiles for a random noise perturbation at various places of a 4 m long film (see text). The dashed curve is the value of the equilibrium enstrophy at the same depth. $\unicode[STIX]{x1D703}=6.4^{\circ }$ , $Re=19.33$ ; $\unicode[STIX]{x1D708}=6.28\times 10^{-6}~\text{m}^{2}~\text{s}^{-1}$ ; $\unicode[STIX]{x1D6FE}=67.0\times 10^{-3}\;\text{N}~\text{m}^{-1}$ ; $\unicode[STIX]{x1D70C}=1070~\text{kg}~\text{m}^{-3}$ .

5.5 Wall shear stress below waves

This model allows us to calculate very easily the wall shear stress i.e. $\unicode[STIX]{x1D70F}_{xz}(z=0)$ , that appears naturally in the model. The momentum and energy equations can be written

To the first order of approximation, we can write $\tilde{\unicode[STIX]{x1D70F}}_{xz}(0)=\tilde{\unicode[STIX]{x1D70F}}_{0}(0)+\unicode[STIX]{x1D700}\tilde{\unicode[STIX]{x1D70F}}_{1}(0)$ . At order 0, $\tilde{\unicode[STIX]{x1D70F}}_{0}(0)=\unicode[STIX]{x1D706}\tilde{h}$ . The expression of $\tilde{\unicode[STIX]{x1D70F}}_{1}(0)$ is given by equations (3.55) and (3.56) using (3.63) and (3.64):

Reverting to dimensional form, the wall shear stress can be written, with the chosen expressions of $f_{1}$ and $f_{2}$ ,

Figure 11.

Thin curve: normalized wall shear stress under a wave in the conditions of Tihon et al. (Reference Tihon, Serifi, Argyriadi and Bontozoglou2006). Thick curve: normalized depth $h/h_{0}$ (normalized order-zero wall shear stress).

The wall shear stress was calculated in the conditions of the experiments of Tihon et al. (2006): $\unicode[STIX]{x1D70C}=1007~\text{kg}~\text{m}^{-3}$ ; $\unicode[STIX]{x1D708}=3.9\times 10^{-6}~\text{m}^{2}~\text{s}^{-1}$ ; $\unicode[STIX]{x1D6FE}=0.06~\text{N}~\text{m}^{-1}$ ; $\unicode[STIX]{x1D703}=5^{\circ }$ ; $Re=35$ . In figure 11, the wall shear stress normalized by its value for the Nusselt flow $gh_{0}\sin \unicode[STIX]{x1D703}/\unicode[STIX]{x1D708}$ is represented as a function of the abscissa $x$ , along with $h/h_{0}$ which is the normalized zero-order wall shear stress.

This result can be compared to the experimental measures and DNS of Tihon et al. (Reference Tihon, Serifi, Argyriadi and Bontozoglou2006). The overall curve is similar. The maximum in the front side of the main hump and the second one are notably smaller in our model. The change of sign of the wall shear stress in a slim region around the first depression ahead of the main hump is in very good agreement with both the experiments and the DNS of Tihon et al. (Reference Tihon, Serifi, Argyriadi and Bontozoglou2006). This indicates a local reversal of the flow which was also predicted by Malamataris et al. (Reference Malamataris, Vlachogiannis and Bontozoglou2002) with a DNS. As noted by Tihon et al. (Reference Tihon, Serifi, Argyriadi and Bontozoglou2006), capillary ripples induce strong stress oscillations which are disproportional to the ripples’ size, there is a stress deficit under the wave crest (deeper in our model) and then a flat profile along the tail of the wave.