Falling liquid films with blowing and suction

Abstract

Flow of a thin viscous film down a flat inclined plane becomes unstable to long-wave interfacial fluctuations when the Reynolds number based on the mean film thickness becomes larger than a critical value (this value decreases as the angle of inclination to the horizontal increases, and in particular becomes zero when the plate is vertical). Control of these interfacial instabilities is relevant to a wide range of industrial applications including coating processes and heat or mass transfer systems. This study considers the effect of blowing and suction through the substrate in order to construct from first principles physically realistic models that can be used for detailed passive and active control studies of direct relevance to possible experiments. Two different long-wave, thin-film equations are derived to describe this system; these include the imposed blowing/suction as well as inertia, surface tension, gravity and viscosity. The case of spatially periodic blowing and suction is considered in detail and the bifurcation structure of forced steady states is explored numerically to predict that steady states cease to exist for sufficiently large suction speeds since the film locally thins to zero thickness, giving way to dry patches on the substrate. The linear stability of the resulting non-uniform steady states is investigated for perturbations of arbitrary wavelength, and any instabilities are followed into the fully nonlinear regime using time-dependent computations. The case of small amplitude blowing/suction is studied analytically both for steady states and their stability. Finally, the transition between travelling waves and non-uniform steady states is explored as the amplitude of blowing and suction is increased.

1. Introduction

The flow of a viscous liquid film down an inclined plane under the action of gravity, inertia and surface tension, is a fundamental problem in fluid mechanics that has received considerable attention both theoretically and experimentally due to the richness of its dynamics and its wide technological applications, e.g. in coating processes and heat or mass transfer enhancement. For sufficiently thin films (with other parameters such as inclination angle and viscosity fixed), the uniform Nusselt solution (Nusselt 1923) is stable, but for thicker layers, the flow is susceptible to interfacial instabilities in the form of two-dimensional travelling waves which propagate down the slope, followed by more complicated time-dependent and three-dimensional behaviour. The linear stability of a uniform film was first considered by Benjamin (1957) and Yih (1963), who used an Orr–Sommerfeld analysis to show that instability first appears at wavelengths that are large compared to the undisturbed film thickness, $h_{s}$ . Using the Nusselt velocity at the free surface, we define a Reynolds number $R={\it\rho}^{2}gh_{s}^{3}\sin {\it\theta}/2{\it\mu}^{2}$ (see (2.4) also) where ${\it\rho}$ is the fluid density, $g$ is the gravitational acceleration, ${\it\theta}$ is the angle of inclination to the horizontal (see figure 1) and ${\it\mu}$ is the viscosity of the fluid; the flow becomes linearly unstable to long waves when $R>R_{c}=(5/4)\cot {\it\theta}$ and we can see that the critical Reynolds number tends to zero as the plate becomes vertical. This result also shows that, for a given angle, the flow can be made unstable by increasing the density and/or the film thickness, or decreasing the viscosity.

Figure 1. Sketch of flow domain showing coordinate system. We consider a fluid layer with mean height $h_{s}$ , bounded along $y=0$ by a rigid planar wall inclined at angle ${\it\theta}$ to the horizontal, and at $y=h(x,t)$ by a free surface. Fluid is injected through the wall, and so the normal velocity at the wall is given by the prescribed function $v=F(x,t)$ .

In order to go beyond linear theory without recourse to direct numerical simulations of the Navier–Stokes equations, a hierarchy of long-wave reduced-dimension models have been developed to analyse in detail the stability and nonlinear development of the flow (see reviews by Craster & Matar 2009; Kalliadasis et al. 2012). The simplest fully nonlinear long-wave model was developed by Benney (1966), who used an expansion in a small slenderness parameter to derive a single evolution equation for the interface height. The Benney equation is valid for Reynolds numbers near the critical value $R_{c}$ (in fact it captures exactly the linear stability threshold). However, for $R$ sufficiently above this critical value, solutions can become unbounded in finite time (Pumir, Manneville & Pomeau 1983), a phenomenon that invalidates the long-wave approximation and is not observed in numerical simulations of the Navier–Stokes equations (Oron & Gottlieb 2002; Scheid et al. 2005). The Benney equation forms a rational basis for the development of asymptotically correct weakly nonlinear models that lead to the Kuramoto–Sivashinsky (KS) equation (see Tseluiko & Papageorgiou 2006, and references therein, for example). The KS equation displays very rich dynamical behaviour, including spatio-temporal chaos. In other canonical asymptotic weakly nonlinear regimes one can derive the generalised (i.e. dispersively modified) KS equation along with electric-field induced instabilities (Tseluiko & Papageorgiou 2006, 2010). In order to overcome the near-critical restrictions of the Benney equation, Shkadov (1969) developed a coupled fully nonlinear long-wave system for the free-surface height and the local mass flux. The Shkadov model avoids finite-time singularities, but under predicts the critical Reynolds number. Using a weighted-residual method, Ruyer-Quil & Manneville (2000) recently developed a new two-degrees-of-freedom long-wave model, which has the correct stability threshold and is well behaved in the nonlinear regime.

In this paper, we consider two such reduced-dimensional models in the presence of blowing and suction, based on the Benney and the weighted-residual methodologies. These two models are the simplest models which correctly predict the critical Reynolds number in the absence of blowing and suction. The Benney equation is easier to use as it is scalar, but it is also less accurate at moderate Reynolds number and can display unphysical finite-time blow up. The weighted residual model has two coupled degrees of freedom, but is more accurate and appears to avoid blow up in numerical computations. Many of the phenomena described in this paper occur for both of these models, or indeed at zero Reynolds number, where the models are identical, and so we can have reasonable confidence that the results do not depend qualitatively on the model used.

In applications it is useful to be able to control the film dynamics. For instance, in coating applications, a stable uniform film is required, whereas in heat or mass transfer applications, efficiency is improved if the flow is non-uniform and attains increased surface area and recirculation regions. Such diverse requirements motivate the introduction of extrinsic modifications to the system in the interest of controlling the dynamics. An example of such a modification is the utilization of a heated substrate, which affects the interfacial dynamics via a combination of Marangoni effects and evaporation (Kalliadasis et al. 2012). Substrate heating thus introduces new modes of instability relating to convection, and lowers the critical Reynolds number. The substrate behaviour can also be altered by allowing chemical coatings, elastic deformations or interactions with flow through porous media (Thiele, Goyeau & Velarde 2009; Ogden, D’Alessio & Pascal 2011; Samanta, Ruyer-Quil & Goyeau 2011; Samanta, Goyeau & Ruyer-Quil 2013), which is often modelled by an effective slip condition. External fields can also be used to stabilise or destabilise the interface. Depending on the fluids used, an applied magnetic field can stabilise the flow (see Hsieh 1965; Shen, Sun & Meyer 1991; Renardy & Sun 1994), whereas an electric field applied normal to the interface can destabilize the flow and in fact drive it to chaotic spatio-temporal dynamics even below critical $R<R_{c}$ (Tseluiko & Papageorgiou 2006, 2010).

One way to modify the dynamics of falling film flows is to topographically structure the substrate. There have been several theoretical and experimental studies of film flows down wavy inclined planes (typically with sinusoidal and step topographies), aiming to explore how topography affects stability and stability criteria such as the critical Reynolds number, how substrate heterogeneity interacts with nonlinear coherent structures and from a practical perspective, how topography induces flows that can be useful in heat or mass transfer by creating regions of recirculating fluid (see for example Tseluiko, Blyth & Papageorgiou (2013), and numerous references therein). The problem is quite complex with several parameters and an overall conclusion of these studies is that topography can either decrease or increase the critical Reynolds number in different regimes.

Inhomogeneous heating of the substrate can generate non-uniform film profiles, even in the absence of inclination (Saprykin et al. 2007). Scheid et al. (2002) used an extension of the Benney model to analyse flow over a substrate with a sinusoidally varying heat distribution, where temperature is coupled to flow via Marangoni effects. They solved the equations to obtain travelling waves in the case of uniform heating, and steady non-uniform solutions in the case of non-uniform heating, and used initial value calculations to demonstrate that imposing heating is able to halt the progress of a travelling wave, leading to stable, steady interface shapes. The combination of localized heating and topography has also been considered and it has been demonstrated that features that form in the isothermal case (e.g. ridge formation ahead of step-down topography) can be removed by suitable heating (see Ogden et al. 2011; Blyth & Bassom 2012). The removal of such a ridge has also been shown to be possible by the imposition of vertical electric fields rather than heating, providing another physical mechanism for interface control (see Tseluiko et al. 2008a , b ).

Suction and injection blowing through an otherwise rigid substrate has well-known applications in stabilizing flows and changing global structures that can negatively affect performance, such as boundary-layer separation, for example. In the types of interfacial flows of interest here, there has been much more limited exploration; Momoniat, Ravindran & Roy (2010) studied the effect of imposing either suction or injection on a spreading drop and found that injection enhances ridge formation, while suction leads to cavities forming on the free surface. As the total mass is not conserved, a steady state is impossible, but they found that both injection and suction are able to disrupt the predominantly streamwise spreading that would occur in the absence of such effects. The total fluid mass is also not conserved for the flows of films and drops over porous substrates (Davis & Hocking 2000); in fact both drops and films are drawn entirely into the substrate in finite time as would be expected. In the analysis of drop evaporation, a number of studies, e.g. Anderson & Davis (1995) and more recently Todorova, Thiele & Pismen (2012) among others, have considered the steady state obtained by imposing injection through the substrate that exactly balances the mass lost to evaporation. In the latter study, the injection profile was imposed according to a Gaussian distribution, and the drop shape is largely independent of this distribution as long as the injection is not too large near the drop contact line (note that a precursor film model was used by Todorova et al. 2012). For continuous falling liquid films over porous substrates there have been linear stability studies invoking a Darcy law in the porous medium and a Beavers–Joseph boundary condition at the liquid substrate boundary (Sadiq & Usha 2008; Usha et al. 2011). It is found that an effective slippage takes place that enhances the instability in the sense that it reduces the critical Reynolds number. Slippage models were investigated further by Samanta et al. (2011) and an alternative porous medium model is proposed and analysed by Samanta et al. (2013).

In this paper we impose blowing or suction through the wall, and perpendicular to it, of fluid that is identical to that of the liquid film. We thus replace the typical no penetration condition at the boundary of the porous wall with blowing and suction applied in the form of a specified, spatially varying, normal velocity. We will assume that there is no slip along the substrate. Beavers & Joseph (1967) argued that tangential slip occurs at the boundary between a fluid-filled porous medium and fluid layer, and that the slip length scales with typical pore size. The no-slip boundary condition is therefore appropriate to those experimental set-ups for which the pores or slits on the substrate which provide the conduit for fluid to enter and leave the wall are much smaller than the typical film thickness. When this condition does not hold, it would be necessary to modify the model to include a non-zero slip length along the substrate, but we do not pursue this here. The blowing and suction of fluid through the substrate affect the total mass in the film and on physical grounds we can anticipate that a net suction would dry the substrate in finite time, whereas a net blowing would increase the total mass and hence the mean thickness at any given time. An increase in thickness would consequently increase the instantaneous Reynolds number since it is proportional to the mean thickness, hence the flow is expected to become more unstable. In this study we will consider the case of blowing/suction that conserves the total film mass (e.g. spatially periodic blowing/suction of zero mean), which is possibly the most interesting case since it sits on the boundary of the net-mass decrease or increase, and hence both stabilising and destabilising phenomena can occur depending on the parameters, as will be seen later.

The rest of the paper is organised as follows. In § 2, we discuss the governing equations and dimensionless parameters, the scaling and statement of the two long wave models and the choice of the blowing and suction function. The numerical methods used to solve these models are described in § 3. In § 4, we explore the steady states and bifurcations obtained for non-zero imposed suction, discovering a non-trivial bifurcation structure even at zero Reynolds number. We also discuss the distinctive effect of the suction on flow streamlines. Linear stability of steady solutions is discussed in § 5, with a focus on stability to perturbations of arbitrary wavelength, and thus the effect of suction on the critical Reynolds number. In § 6, we investigate the effect of imposing suction on the travelling waves which occur above the critical Reynolds number. In § 7, we review the various initial value calculations performed in this paper and provide further results. Finally, we present our conclusions in § 8.

2. Problem formulation

2.1. Non-dimensionalisation and scaling

We wish to determine the evolution of a falling liquid film, with mean thickness $h_{s}$ , on a slope inclined at angle ${\it\theta}$ to the horizontal. We model the liquid as a Newtonian fluid of constant dynamic viscosity ${\it\mu}$ and density ${\it\rho}$ , and the air as a hydrodynamically passive region of constant pressure $p_{a}$ . The coefficient of surface tension across the air–liquid interface is ${\it\gamma}$ . We assume that the film flow is two-dimensional, with no variations in the cross-stream direction. We use coordinates as illustrated in figure 1; $x$ is the down-slope coordinate, and $y$ is the coordinate in the direction perpendicular to the slope, so that the wall is located at $y=0$ and the interface is defined as $y=h(x,t)$ . We denote the velocity components in the $x$ and $y$ directions as $u$ and $v$ , respectively.

The dimensionless film flow is governed by the Navier–Stokes equations

(2.1) $$\begin{eqnarray}\displaystyle & \displaystyle R(u_{t}+uu_{x}+vu_{y})=-p_{x}+2+u_{xx}+u_{yy}, & \displaystyle\end{eqnarray}$$

(2.2) $$\begin{eqnarray}\displaystyle & \displaystyle R(v_{t}+uv_{x}+vv_{y})=-p_{y}-2\cot {\it\theta}+v_{xx}+v_{yy}, & \displaystyle\end{eqnarray}$$

which are coupled to the incompressibility condition

(2.3) $$\begin{eqnarray}u_{x}+v_{y}=0.\end{eqnarray}$$

Here we have non-dimensionalised the equations using the mean film thickness $h_{s}$ as the length scale, the Nusselt surface speed of a flat film $U_{s}={\it\rho}gh_{s}^{2}\sin {\it\theta}/2{\it\mu}$ (Nusselt 1923) as the velocity scale, $h_{s}/U_{s}$ as the time scale and ${\it\mu}U_{s}/h_{s}$ as the pressure scale. We note that this scaling, based on surface speed, is the same scaling as used by Tseluiko et al. (2013) to study the influence of wall topography on the stability of flow down an inclined plane. We define the Reynolds number $R$ and the capillary number $C$ based on the surface speed $U_{s}$ , so that

(2.4a,b ) $$\begin{eqnarray}R=\frac{{\it\rho}h_{s}U_{s}}{{\it\mu}},\quad C=\frac{{\it\mu}U_{s}}{{\it\gamma}}.\end{eqnarray}$$

We will suppose that the imposition of suction boundary conditions at the wall does not alter the no-slip condition, but does affect the no-penetration condition, so that the complete boundary conditions at the wall, $y=0$ , are

(2.5a,b ) $$\begin{eqnarray}u=0,\quad v=F(x,t).\end{eqnarray}$$

At the interface, $y=h(x,t)$ , the tangential and normal components of the dynamic stress balance condition yield

(2.6) $$\begin{eqnarray}\displaystyle & \displaystyle (v_{x}+u_{y})(1-h_{x}^{2})+2h_{x}(v_{y}-u_{x})=0, & \displaystyle\end{eqnarray}$$

(2.7) $$\begin{eqnarray}\displaystyle & \displaystyle p-p_{a}-\frac{2}{1+h_{x}^{2}}(v_{y}+u_{x}h_{x}^{2}-h_{x}(v_{x}+u_{y}))=-\frac{1}{C}\frac{h_{xx}}{(1+h_{x}^{2})^{3/2}}, & \displaystyle\end{eqnarray}$$

respectively. The kinematic boundary condition on the interface can be written as

(2.8) $$\begin{eqnarray}h_{t}=v-uh_{x}.\end{eqnarray}$$

We can use (2.3) together with (2.5) to rewrite (2.8) as

(2.9) $$\begin{eqnarray}h_{t}-F(x,t)+q_{x}=0,\end{eqnarray}$$

where $q(x,t)$ is the streamwise flow rate, defined as

(2.10) $$\begin{eqnarray}q(x,t)=\int _{0}^{h}u(x,y,t)\,\text{d}y.\end{eqnarray}$$

2.2. Long-wave evolution equations

We now seek solutions with wavelength $L$ relative to the thickness of the undisturbed fluid layer, with $L\gg 1$ , and so we introduce the long-wave parameter ${\it\delta}=1/L$ . We derive two first-order long-wave models, based on an asymptotic expansion in the long-wave parameter (a Benney-type model, see Benney 1966) and on a weighted-residual method (following the approach of Ruyer-Quil & Manneville (2000)); derivations of both models are presented in appendix A. We assume that $\cot {\it\theta}=O(1)$ . To retain inertial effects, we additionally assume that $R=O(1)$ , and following Gjevik (1970) we choose $C=O({\it\delta}^{2})$ to maintain the effect of surface tension. We choose the canonical scaling $F=O({\it\delta})$ so that the imposed suction can enter and compete with the perturbed flow and hence facilitate possible instability enhancement or reduction.

The essential task of the derivation is to estimate the flow rate $q(x,t)$ for a non-uniform film. In both models, the mass conservation equation is unchanged from (2.9), and so we have

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

However, the two models differ in their treatment of nonlinearities in the momentum equation, and thus yield different equations for  $q$ .

In the Benney equation (see § A.1), $q$ is slaved to the interface height $h$ , and is given by

(2.12) $$\begin{eqnarray}q(x,t)=\frac{2h^{3}}{3}-\frac{h^{3}}{3}\left(2h\cot {\it\theta}-\frac{h_{xx}}{C}\right)_{x}+R\left(\frac{8h^{6}h_{x}}{15}-\frac{2h^{4}F}{3}\right).\end{eqnarray}$$

The first-order weighted-residual approach (see § A.2) instead yields an evolution equation for  $q$ :

(2.13) $$\begin{eqnarray}\frac{2}{5}Rh^{2}q_{t}+q=\frac{2h^{3}}{3}-\frac{h^{3}}{3}\left(2h\cot {\it\theta}-\frac{h_{xx}}{C}\right)_{x}+R\left(\frac{18q^{2}h_{x}}{35}-\frac{34hqq_{x}}{35}+\frac{hqF}{5}\right).\end{eqnarray}$$

The Benney and weighted-residual equations are identical when $R=0$ . For $R\neq 0$ , the Benney equation has a single degree of freedom $h(x,t)$ , while the weighted-residual model has two degrees of freedom, $h(x,t)$ and $q(x,t)$ . As a result, the weighted-residual equations can in principle exhibit richer behaviour. However, it seems that because of this additional complexity, the weighted-residual equations in fact support bounded solutions across a greater range of parameter space (Scheid et al. 2005).

2.3. Choice of blowing and suction function

In time-dependent evolution with periodic boundary conditions, the mean layer height is conserved only if the imposed flux function $F$ has zero mean, and hence steady states can only exist if $F$ has zero mean. Conserved mean layer height is the natural state for numerical calculations in a periodic domain, and is sometimes called ‘closed’ conditions as there is no net flux out of the domain. However, in experiments, closed conditions are not easy to implement, and so experimental realisations more typically impose the fluid flux at the domain inlet. In this second case, known as ‘open’ conditions, there is no direct control over the mean layer height.

In the absence of blowing and suction, both the open and closed systems support uniform flow via the Nusselt solution, and thus we obtain the same scaling whether based on the mean layer height or mean flux. In order to investigate the influence of suction, we can consider steady states where the mean layer height remains fixed in time for closed conditions, or where there is no change to the mean flux for open conditions. However, the bifurcation diagrams for fixed mean layer height correspond most naturally to statements about time evolution in closed conditions, as the layer height is effectively constrained in both cases. Most of the results we present are for fixed mean layer height, but we will also present some results for fixed flow rate. In principle, by varying either $\bar{h}$ or $\bar{q}$ , we should be able to map results between the two conditions, but the bifurcations that occur as $A$ is varied, for example, may depend on whether layer height or flux is constrained.

Examination of the derivation presented in appendix A reveals that the long-wave equations (2.11)–(2.13) also apply if $F$ is unsteady, so long as $F$ varies on dimensionless time scales no shorter than $O({\it\delta})$ , or the dimensional time scale ${\it\delta}h_{s}/U_{s}$ . While all results presented in this paper are for the case of steady $F$ , we note that time derivatives of the vertical velocity, and hence $F$ , would not feature in the equations until the next order in ${\it\delta}$ . We are therefore free, at this order, to impose a time dependence on $F$ , or even choose $F$ in response to the film evolution. The latter formulation would be particularly useful in feedback control studies, which we explore in Thompson et al. (2015).

In the rest of this manuscript, we will consider the simplest functional form for $F$ with zero mean, i.e. a single harmonic mode:

(2.14) $$\begin{eqnarray}F(x)=A\cos (mx)=A\cos \left(\frac{2{\rm\pi}x}{L}\right).\end{eqnarray}$$

This function $F(x)$ has the symmetry that the transformation $A\rightarrow -A$ is recovered by translation in $x$ by a distance $L/2$ , and so translationally-invariant solution measures, such as the critical Reynolds number for onset of instability, must be even in  $A$ .

3. Numerical methods

The numerical calculations that we perform are of three types: computation of steady periodic solutions and their bifurcation structure, linear stability calculations of such steady states to perturbations of arbitrary wavelength and nonlinear time-evolution via initial value problems. We conduct these calculations using the continuation software package Auto-07p (Doedel et al. 2009) and Matlab.

The first task, of computing steady solutions and exploring their bifurcation structure, was conducted in Auto-07p by formulating the problem as a boundary-value problem with periodic boundary conditions. The Auto-07p code is spatially adaptive, and unlikely to return spurious solutions. We are particularly interested in limit point, pitchfork and Hopf bifurcations. The first two of these correspond to bifurcations of steady states, and so can be detected and tracked using the same formulation as for standard steady states. With regard to Hopf bifurcations, here we are concerned with Hopf bifurcations with respect to time, whereby a steady state becomes oscillatory in time when subject to a perturbation of fixed spatial wavelength. If the perturbation satisfies periodic spatial boundary conditions, the instability will in fact be oscillatory in both space and time. We used Auto-07p to track individual Hopf bifurcations by manually augmenting the steady system with a boundary-value problem for the spatially periodic but non-constant eigenfunction, with the temporal eigenvalue determined as part of the solution.

The linear stability calculations were performed in Matlab, and we used a pseudo-spectral method for the spatial discretisation. After spatial discretisation, the governing partial differential equation becomes a large system of coupled first-order ordinary differential equations, which we can write in the general form

(3.1) $$\begin{eqnarray}\boldsymbol{F}(\boldsymbol{u},\dot{\boldsymbol{u}})=0\end{eqnarray}$$

so that steady solutions $\boldsymbol{u}_{0}$ satisfy $\boldsymbol{F}(\boldsymbol{u}_{\mathbf{0}},\mathbf{0})=0$ . In order to determine the linear stability of a steady solution, we suppose that

(3.2) $$\begin{eqnarray}\boldsymbol{u}(t)=\boldsymbol{u}_{0}+{\it\epsilon}\boldsymbol{v}\exp ({\it\lambda}t),\quad {\it\epsilon}\ll 1.\end{eqnarray}$$

We now expand (3.1) for small ${\it\epsilon}$ , to obtain

(3.3) $$\begin{eqnarray}\unicode[STIX]{x1D645}\boldsymbol{v}+{\it\lambda}\unicode[STIX]{x1D648}\boldsymbol{v}=0,\quad \unicode[STIX]{x1D645}=\left.\frac{\partial \boldsymbol{F}}{\partial \boldsymbol{u}}\right|_{\boldsymbol{u}_{\mathbf{0}},\mathbf{0}},\quad \unicode[STIX]{x1D648}=\left.\frac{\partial \boldsymbol{F}}{\partial \dot{\boldsymbol{u}}}\right|_{\boldsymbol{u}_{\mathbf{0}},\mathbf{0}},\end{eqnarray}$$

which is a generalised eigenvalue problem for ${\it\lambda}$ and $\boldsymbol{v}$ , where $\unicode[STIX]{x1D645}$ is the Jacobian matrix and $\unicode[STIX]{x1D648}$ is the mass matrix. As very few points were needed for the spatial discretisation (we typically used 99 equally spaced points), we solved the eigenvalue problem (3.3) directly in Matlab. We used Floquet multipliers to determine linear stability to perturbations of arbitrary wavelength, and so modified the Jacobian matrix to account for these when necessary.

We note that in the Benney equations, the only time derivatives are those of interface height $h$ , while the weighted-residual equations also feature time derivatives of flux $q$ . This means that, for the same spatial discretisation, there are twice as many eigenmodes for the weighted-residual equations as for the Benney calculations. We found that the weighted-residual calculations were prone to spurious eigenmodes, which we removed by careful comparison of the eigenvalue spectrum at different spatial resolutions. We also neglected the neutrally-stable eigenmode corresponding to increasing the total volume of fluid in the domain, which arises in both sets of equations.

Time evolution calculations were always performed in a fixed spatially-periodic domain. The spatial problem was discretised via a pseudo-spectral method, while time derivatives were handled via a second-order backward finite-difference scheme (BDF2). The resulting code is fully implicit, and solved via direct Newton iteration.

The spatial discretization in the Matlab code was verified by comparing the steady solutions and bifurcation structure obtained in Matlab to those obtained in Auto-07p. Further validation was obtained by comparison of numerical results to analytical solutions for the shape of small-amplitude steady states and to analytical results for the linear stability of uniform and small-amplitude states.

4. Bifurcation structure for steady states

In the absence of blowing or suction, the only spatially-periodic steady state is a uniform film. Introducing periodic suction naturally imposes a spatial structure on the solutions, and means that any steady states must be non-uniform. When $R>0$ , the solutions and bifurcations differ between the two long-wave models. The weighted-residual model avoids the blow-up behaviour sometimes exhibited by the Benney equation, and more accurately represents the behaviour of the Navier–Stokes equations at moderate Reynolds number. We will generally present results for the weighted-residual model when the two models differ, but we note that a non-trivial bifurcation structure emerges even at zero Reynolds number, where the models are identical.

4.1. Steady solutions at small $A$

We begin by considering the effect of small-amplitude forcing, in the form of blowing and suction, on the uniform steady state $h=1$ . We choose $F=A\cos mx$ where $m=2{\rm\pi}/L$ , and seek a steady solution for $h$ and $q$ when $|A|\ll 1$ of the form

(4.1a,b ) $$\begin{eqnarray}h=1+A\,\text{Re}({\hat{H}})+O(A^{2}),\quad q={\textstyle \frac{2}{3}}+A\,\text{Re}(\hat{Q})+O(A^{2}),\end{eqnarray}$$

where $\text{Re}()$ indicates the real part. The mass conservation equation (2.11) immediately supplies $\hat{Q}=\exp (\text{i}mx)/(\text{i}m)$ . However, the complex quantity ${\hat{H}}$ differs between the two long-wave models. The Benney result, obtained by linearising (2.12), is

(4.2) $$\begin{eqnarray}{\hat{H}}_{\mathit{Benney}}=\frac{1+{\displaystyle \frac{2}{3}}\text{i}mR}{2\text{i}m+{\displaystyle \frac{2}{3}}m^{2}\cot {\it\theta}-{\displaystyle \frac{8}{15}}Rm^{2}+{\displaystyle \frac{m^{4}}{3C}}}\exp (\text{i}mx),\end{eqnarray}$$

while the weighted-residual equation (2.13) gives

(4.3) $$\begin{eqnarray}{\hat{H}}_{\mathit{WR}}=\frac{1+{\displaystyle \frac{18}{35}}\text{i}mR}{2\text{i}m+{\displaystyle \frac{2}{3}}m^{2}\cot {\it\theta}-{\displaystyle \frac{8}{35}}Rm^{2}+{\displaystyle \frac{m^{4}}{3C}}}\exp (\text{i}mx).\end{eqnarray}$$

The resulting solutions for ${\hat{H}}$ are illustrated in figure 2 for $L=10$ and $L=40$ .

Figure 2. The steady-state deflection of a uniform film for $F=A\cos (2{\rm\pi}x/L)$ for $A\ll 1$ , with $L=10$  (a) and $L=40$  (b),  $C=0.05$ , ${\it\theta}={\rm\pi}/4$ and $R=2$ . The $O(A)$ correction $\text{Re}({\hat{H}})$ is shown for the Benney equations (solid line) and the weighted-residual equations (dashed line), with ${\hat{H}}$ defined in (4.2) and (4.3), respectively. The dotted line indicates the scaled suction profile $F(x)/A$ .

The two expressions (4.2) and (4.3) are equal only if $R=0$ ; otherwise both the magnitude and phase of ${\hat{H}}$ differ between the two equations. At $L=10$ , the predictions obtained via the two models are in reasonable agreement, but they are indistinguishable when $L=40$ . Returning to the variables of the long-wave derivation, we set $m={\it\delta}M$ and $C={\it\delta}^{2}\widehat{C}$ , and expand for small ${\it\delta}$ ; we find that both models yield

(4.4) $$\begin{eqnarray}{\hat{H}}=\frac{1}{2{\it\delta}\text{i}M}+\left(\frac{\cot {\it\theta}}{6}+\frac{M^{2}}{12\widehat{C}}+\frac{R}{5}\right)+O({\it\delta}),\end{eqnarray}$$

which is the expected order of agreement given that terms beyond the first order in ${\it\delta}$ were neglected in the derivation of each model. We note from (4.4) that the magnitude of ${\hat{H}}$ is inversely proportional to $M$ at leading order, so that for fixed $A$ , the maximum perturbation to the interface grows linearly with the wavelength $L$ of the imposed blowing and suction. The long-wave expression for ${\hat{H}}$ (4.4) is in general complex, so there is a phase shift between $h$ and  $F$ . As ${\it\delta}\rightarrow 0$ , the dominant term is purely imaginary, and so the phase shift tends to ${\rm\pi}/2$ in the long-wave limit.

In order to calculate the small $A$ correction to the mean flux analytically, we must expand the steady solution $h=H(x)$ , $q=Q(x)$ to $O(A^{2})$ . For steady states with $F(x)=A\cos mx$ , we integrate the mass conservation equation (2.11) to yield

(4.5) $$\begin{eqnarray}Q(x)=\frac{2}{3}+\frac{A\sin mx}{m}+Q_{2}A^{2}+O(A^{4}),\end{eqnarray}$$

where we have used the fact that the spatially-averaged flux is even in $A$ . We then expand the interface height in $A$ as

(4.6) $$\begin{eqnarray}H(x)=1+A(p_{1}\cos mx+p_{2}\sin mx)+A^{2}(p_{3}\cos 2mx+p_{4}\sin 2mx)+O(A^{3}).\end{eqnarray}$$

We solve the flux equation, either (2.12) or (2.13), at $O(A)$ and $O(A^{2})$ to determine the constants $p_{1}$ , $p_{2}$ , $p_{3}$ , $p_{4}$ and $Q_{2}$ . The full result for $Q_{2}$ in the Benney equation is

(4.7) $$\begin{eqnarray}Q_{2}=\frac{5(45/m^{2}-28R^{2}-20R(2\cot {\it\theta}+m^{2}/C))}{Rm^{2}(64R-80(2\cot {\it\theta}+m^{2}/C))+25(4m^{2}\cot {\it\theta}(\cot {\it\theta}+m^{2}/C)+m^{6}/C^{2}+36)}\end{eqnarray}$$

and in the weighted-residual equation is

(4.8) $$\begin{eqnarray}Q_{2}=\frac{9(1225/m^{2}-264R^{2}-245R(2\cot {\it\theta}+m^{2}/C))}{Rm^{2}(576R-1680(2\cot {\it\theta}+m^{2}/C))+1225(4m^{2}\cot {\it\theta}(\cot {\it\theta}+m^{2}/C)+m^{6}/C^{2}+36)}.\end{eqnarray}$$

It is easy to check that the two results are identical if $R=0$ .

In the long-wave limit $m={\it\delta}M$ , $C={\it\delta}^{2}\widehat{C}$ , with ${\it\delta}\ll 1$ , both the Benney and weighted-residual equations yield

(4.9) $$\begin{eqnarray}Q_{2}=\frac{1}{4{\it\delta}^{2}M^{2}}+O(1).\end{eqnarray}$$

Thus we obtain agreement at the first two orders in ${\it\delta}$ . Furthermore, if blowing and suction are applied at sufficiently long wavelength, and at small amplitude, the mean down-slope flux is always increased relative to its value for a flat film. In fact, this flux increase can be explained from the leading-order flux equation: $q=2h^{3}/3+O({\it\delta})$ , with steady states satisfying $q_{x}=F(x)$ while subject to the constraint that the mean value of $h$ is 1. With non-zero $F$ , any steady states are non-uniform. In order to compare the effect of the net flux in a non-uniform flow $h(x)$ to that of a uniform film of height 1, we consider

(4.10) $$\begin{eqnarray}\int _{-\infty }^{\infty }[q(h)-q(1)]\,\text{d}x=\frac{2}{3}\int _{-\infty }^{\infty }[h^{3}-1]\,\text{d}x=\frac{2}{3}\int _{-\infty }^{\infty }[3(h-1)+(2+h)(h-1)^{2}]\,\text{d}x.\end{eqnarray}$$

The first term on the right-hand side of (4.10) vanishes because $h$ has mean value 1. The remaining term is non-negative, and so choosing any non-zero height profile, while respecting the volume constraint, leads to an increased mean flux at this order. However, the full result at $O(A^{2})$ , given by (4.7) or (4.8), can be negative if $m$ is not sufficiently small, and hence blowing and suction can increase or decrease the mean downslope flux.

4.2. Steady solutions as $A$ increases

As $A$ increases, the film profile deviates more significantly from the uniform state, and the nonlinearities in (2.12) and (2.13) can lead to bifurcations between steady solutions. Figure 3 shows the behaviour of steady solutions to the weighted-residual equations as $A$ is varied for a selection of $R$ at fixed $L$ , ${\it\theta}$ and $C$ subject to constrained mean layer height 1.

Figure 3. Bifurcation structure for steady solutions to the weighted-residual equations as $A$ increases, subject to fixed mean layer height of 1, and calculated in a periodic domain of length $2L$ . Here $F=A\cos (2{\rm\pi}x/L)$ , $L=10$ , $C=0.05$ , ${\it\theta}={\rm\pi}/4$ , with $R=0$ , 3, 6, 9 in (a–d), respectively. The shaded inset solutions lie on the dashed branch of subharmonic steady solutions, which is unstable; all other solutions are harmonic, with period $L$ . Solutions filled in black are stable, while white solutions are unstable. We also indicate pitchfork bifurcations at $A=A_{P}$ (♦), limit points at $A=A_{LP}$ (●), Hopf bifurcations within the domain of length $L$ (▪) and states with $\min q=0$ ( $\star$ ).

For each $R$ , the only steady solution at $A=0$ is the uniform state with $h=1$ . All steady solutions for $A>0$ are non-uniform, and the extent of this non-uniformity increases with $A$ when $A$ is moderate. However, for sufficiently large $A$ , there are no steady solutions that describe continuous liquid films due to a limit point in the bifurcation diagram. We can follow the steady solution branches around the limit point, and find that each branch eventually terminates when the minimum film height tends to zero, so that the layer dries up at some position. The film height profiles for these drying states are shown as insets in figure 3.

When $A$ is sufficiently large that there are no steady solutions, we explore the system behaviour by conducting time-dependent simulations in a periodic domain of length $L$ , with a typical result shown in figure 4(a). We find that the film thins rapidly in some localised regions, and is able to dry in finite time due to the fixed speed of fluid removal. We have also conducted unsteady simulations starting from slightly-perturbed states on the solution branch with the lower minimum height. We find that these steady states are unstable, and the instability usually manifests as thinning and drying behaviour, rather than tending towards the linearly stable steady states with a larger minimum film thickness.

The solutions shown in figure 3 are all computed subject to a fixed mean layer height of 1, and are plotted according to the minimum value of $h$ in a single period. As an alternative solution measure, we consider the mean flux $\bar{q}$ , averaged over one spatial period. Figure 5(a,c) shows the bifurcation structure plotted in terms of the mean flux; we find that imposing blowing and suction at finite wavelength can either increase or decrease the mean flux. Figure 5(b,d) shows bifurcation diagrams for two values of $R$ with fixed $\bar{q}=2/3$ and mean height free to vary (this condition is appropriate to experimental investigations performed under open conditions, corresponding to fixed mean down-slope flux). We find that, as was the case for the fixed $\bar{h}=1$ calculations, there is a limit point corresponding to a maximum value of $A$ with steady solutions. However, time-dependent calculations for open conditions require non-periodic boundary conditions, which we have not pursued here.

4.3. Subharmonic steady states

The steady solutions discussed in § 4.2 are limited to cases where both the steady solution and the imposed blowing and suction are periodic in terms of $q$ and $h$ with the same wavelength  $L$ . However, subharmonic steady solutions may also exist for which the steady solution is spatially periodic with period  $nL$ , where $n$ is an integer equal to 2 or greater; we calculate subharmonic steady states by continuation in a domain of length $nL$ . For each case in figure 3, and also for the fixed $\bar{q}=2/3$ calculations shown in figure 5, we have detected subharmonic solution branches with $n=2$ . In each of these cases, the subharmonic solutions emerge via a subcritical pitchfork bifurcation at $A=A_{P}$ in the steady-solution structure, and the subharmonic steady states are all unstable. The unstable eigenmode of these states results in a drying process similar to that occurring for unstable harmonic steady states.

Figure 4. Time evolution subject to periodic-boundary conditions, for $R=0$ , $L=10$ , $C=0.05$ , $F=A\cos (2{\rm\pi}x/L)$ , with (a) initial condition $h=1$ in a domain of length $L$ for $A_{LP}<A$ ( $A=0.8$ ) so that there are no steady states, and (b) initial condition $h=1+0.001\cos {\rm\pi}x/10$ in a domain of length $2L$ with $A_{P}<A<A_{LP}$ ( $A=0.72$ ) so that steady states exist but are unstable. In both cases, the film height vanishes in finite time at a single point in the domain.

Figure 5. The bifurcation structure for steady solutions to the weighted residual equations as $A$ increases, subject to fixed mean height (a,c) and fixed mean flux (b,d). Here $L=10$ , $C=0.05$ and ${\it\theta}={\rm\pi}/4$ . The symbols mark the same bifurcations as in figure 3, but we do not show Hopf bifurcations here. Solution branches terminate at a point where the minimum layer height vanishes. (a)  $\bar{h}=1$ , $R=0$ . (b)  $\bar{q}=2/3$ , $R=0$ . (c)  $\bar{h}=1$ , $R=2$ . (d)  $\bar{q}=2/3$ , $R=2$ .

As the subharmonic steady states shown in figure 3 are unstable, the only possible stable steady states are harmonic, with period $L$ . However, for $A>A_{P}$ , the period- $L$ steady states are also unstable to perturbations of wavelength $2L$ ; this is a consequence of the exchange of stability that occurs at the pitchfork bifurcation. The subharmonic instability eventually results in a film thinning and drying event, but at only one of the local minima within the domain of length $2L$ . One such drying sequence is illustrated in figure 4(b).

4.4. Streamfunction and flow reversal

In addition to altering the interface height, the imposed blowing and suction also fundamentally affect the arrangement of stagnation points and streamlines within the fluid film. Under the weighted-residual formulation, the velocity field $(u,v)=({\it\psi}_{y},-{\it\psi}_{x})$ , where ${\it\psi}$ is the streamfunction, is given by

(4.11a,b ) $$\begin{eqnarray}u={\it\psi}_{y}=\frac{3q}{h}\left(\frac{y}{h}-\frac{y^{2}}{2h^{2}}\right)+O({\it\delta}),\quad v=-{\it\psi}_{x},\end{eqnarray}$$

with boundary condition $v=F(x)$ on $y=0$ . Recalling that $F(x)=q^{\prime }(x)$ for steady solutions, we can write the solution for ${\it\psi}$ , up to the addition of an arbitrary constant, as

(4.12) $$\begin{eqnarray}{\it\psi}=-q(x)+\frac{3q(x)}{h(x)}\left(\frac{y^{2}}{2h}-\frac{y^{3}}{6h^{2}}\right)+O({\it\delta}).\end{eqnarray}$$

The same result applies to $O({\it\delta})$ in the Benney equations. According to (4.12), there are no stagnation points in $0<y\leqslant h$ if $q\neq 0$ , and points along the wall are stagnation points if and only if $F(x)=0$ . This means that, in the absence of blowing and suction, every point on the wall is a stagnation point. With non-zero blowing and suction, there are isolated stagnation points on the wall at the zeros of $F(x)$ .

We now examine the steady flow near such an isolated stagnation point, located at $x=x_{0}$ , $y=0$ , with $F(x_{0})=0$ . As $q^{\prime }(x)=F(x)$ , we can write

(4.13) $$\begin{eqnarray}q\sim q(x_{0})+F^{\prime }(x_{0})\frac{(x-x_{0})^{2}}{2}+O(x-x_{0})^{3}.\end{eqnarray}$$

Substituting (4.13) into (4.12) yields

(4.14) $$\begin{eqnarray}{\it\psi}\sim -q(x_{0})-\frac{F^{\prime }(x_{0})x^{2}}{2}+\frac{3q(x_{0})y^{2}}{2h(x_{0})^{2}}+O((x-x_{0})^{2}+(y-y_{0})^{2}).\end{eqnarray}$$

If $F^{\prime }(x_{0})/q(x_{0})>0$ , the stationary point at $x_{0}$ is a saddle point of ${\it\psi}$ , and the streamlines are locally straight lines, such that

(4.15) $$\begin{eqnarray}\frac{y^{2}}{x^{2}}=\frac{F^{\prime }(x_{0})h(x_{0})^{2}}{3q(x_{0})}.\end{eqnarray}$$

The streamlines become increasingly normal to the wall as $q(x_{0})\rightarrow 0$ . In contrast, if $F^{\prime }(x_{0})/q(x_{0})<0$ , the stationary point at $x_{0}$ is a local extremum of ${\it\psi}$ , and the streamlines are locally elliptical, with

(4.16) $$\begin{eqnarray}-F^{\prime }(x_{0})x^{2}+\frac{3q(x_{0})y^{2}}{h(x_{0})^{2}}=\text{const}.\end{eqnarray}$$

For any smooth, periodic, non-zero $F(x)$ with mean zero, there must be some stagnation points $x_{0}$ with $F^{\prime }(x_{0})>0$ and others with $F^{\prime }(x_{0})<0$ . There are no steady solutions in figure 3 with negative $h$ , but the same is not true for $q$ ; for each $R$ with solutions shown in figure 3, there is an amplitude $A^{\ast }$ above which all steady solution have $q(x)<0$ over a finite interval in $x$ . As $q_{x}=F$ , the minimum of $q$ occurs at a point $x_{0}$ with $F(x_{0})=0$ and $F^{\prime }(x_{0})>0$ .

There are no stagnation points inside the fluid domain, and so streamlines never cross. For small $A$ , $q(x)>0$ for all $x$ , and so stagnation points with $F^{\prime }(x)>0$ are saddle points of ${\it\psi}$ , while those with $F^{\prime }(x)<0$ are extrema of ${\it\psi}$ . Two examples of the flow field can be seen in figure 6(a,b). Streamlines emanate into the fluid domain from each stagnation point with $F^{\prime }(x)>0$ , and these stagnation points are connected by a streamline which separates the fluid into a layer in which fluid particles propagate down the plane, and a layer in which fluid particles must enter and leave the flow domain via injection through the walls. At $A=0$ , all particles propagate. As $A$ increases, the propagating layer thins and the recirculating layer thickens.

At $A=A^{\ast }$ , the stagnation point at $x=x_{0}$ , $y=0$ switches from a saddle point to an extremum of ${\it\psi}$ . The corresponding change in the flow field is illustrated in the transition from figure 6(b) to (c). For $A>A^{\ast }$ , $q(x)$ is negative for some $x$ , and hence by (4.11) the horizontal velocity is directed up slope for all $y$ . As a result, no particles can propagate more than one period down the plane, and so the propagating layer vanishes for $A>A^{\ast }$ . The width of the region with up-slope flow increases with $A$ , until eventually there are no further steady solutions. Nonlinear time evolution calculations at large $A$ show that the film dries at a point just upstream of the negative- $q$ stagnation point. An instantaneous snapshot of the system at the moment of drying is shown in figure 6(e).

If the mean downstream flow rate were held fixed at $2/3$ , rather than the mean film thickness constrained to 1, a largely similar sequence of events would occur to that shown in figure 6. If $F=A\cos mx$ , for steady solutions, we immediately have $q=2/3+A\sin mx/m$ , and so flow reversal would occur for any steady solutions with $A>2m/3$ .

Figure 6. Solutions for the interface shape, velocity field and stagnation points for ${\it\theta}={\rm\pi}/4$ , $C=0.05$ , $L=10$ , $R=0$ . Periodicity is enforced with $L=10$ , but solutions are plotted over two periods and are shown with aspect ratio 1. Stagnation points with $q(x)F^{\prime }(x)<0$ and with $q(x)F^{\prime }(x)>0$ are indicated by ○ and ▫, respectively. (a) Steady state for $A=0.2$ , (b)  $A=0.4$ , (c)  $A=0.5$ , (d)  $A=0.6$ , (e) there are no steady states for $A=0.8$ ; this is a final snapshot before drying. In (a,b), $q>0$ everywhere, so a streamline emanating from the stagnation point with $q(x)F^{\prime }(x)>0$ divides fluid into particles which never meet the wall (yellow), and particles which enter and leave through the wall (blue). For $A>A^{\ast }=0.46$ , all steady solutions have regions of negative $q$ , corresponding to a region of upstream flow near the stagnation point (between the vertical lines), and all fluid particles must reach the wall.

It is possible that our predictions of regions of reversed flow for steady states with large amplitude blowing and suction are related to the phenomenon of capillary backflow, as described by Malamataris, Vlachogiannis & Bontozoglou (2002) and Dietze, Leefken & Kneer (2008). Backflow is observed in large-amplitude travelling waves at moderately large Reynolds number, and occurs in the form of eddies near the wall. In contrast, in our system, the reversed flow persists through the thickness of the film, can occur at zero Reynolds number and would disappear if the blowing and suction was removed.

5. Linear stability

5.1. Stability of uniform state

In the absence of imposed suction, the only steady state with mean layer height unity is the Nusselt film solution, with $h=1$ and $q=2/3$ . We linearise about this state, and seek eigenmodes proportional to $\exp (\text{i}kx+{\it\lambda}t)$ . The Benney model yields the eigenvalue ${\it\lambda}$ directly, as

(5.1) $$\begin{eqnarray}{\it\lambda}=-2\text{i}k+k^{2}\left(\frac{8R}{15}-\frac{2}{3}\cot {\it\theta}\right)-\frac{k^{4}}{3C}.\end{eqnarray}$$

For the weighted-residual model, the linear stability analysis yields a quadratic equation for  ${\it\lambda}$ :

(5.2) $$\begin{eqnarray}\frac{2R}{5}{\it\lambda}^{2}+{\it\lambda}\left(1-\text{i}k\frac{68R}{105}\right)+2\text{i}k+k^{2}\left(\frac{2\cot {\it\theta}}{3}+\frac{k^{2}}{3C}-\frac{8R}{35}\right)=0.\end{eqnarray}$$

In both models, the stability threshold for fixed $k$ is

(5.3) $$\begin{eqnarray}R<R_{H}\equiv \frac{5}{4}\cot {\it\theta}+\frac{5}{8C}k^{2}.\end{eqnarray}$$

The linear stability threshold (5.3) is also valid for perturbations with long wavelengths in full solutions to the Navier–Stokes equations (Benjamin 1957). At $R=R_{H}$ , both (5.1) and (5.2) have a root with negative real part for $R<R_{H}$ , the value ${\it\lambda}={\it\lambda}_{0}=-2\text{i}k$ at $R=R_{H}$ and positive real part for $R>R_{H}$ . We note that (5.2) also has a second root for ${\it\lambda}$ , but this root always has negative real part when $R>0$ .

5.2. Perturbations of arbitrary wavelength about a periodic base state

When the base state is spatially uniform, the eigenmodes can be written as $\text{Re}(\exp (\text{i}kx+{\it\lambda}t))$ , and by calculating ${\it\lambda}$ for all real $k$ , we can determine linear stability for perturbations of all wavelengths. However, when non-zero suction and blowing is imposed, the base state is not uniform, and so the eigenmodes are no longer simple exponential functions. Fortunately, Floquet–Bloch theory allows us to compactly describe eigenmodes of arbitrary wavelength.

We suppose that we are given a steady solution $h=H(x)$ , $q=Q(x)$ to the equations with fixed, non-zero $F(x)$ , and consider the evolution of arbitrary small perturbations, so that

(5.4a,b ) $$\begin{eqnarray}h(x,t)=H(x)+{\it\epsilon}\,\text{Re}\{{\hat{h}}(x,t)\}+O({\it\epsilon}^{2}),\quad q(x,t)=Q(x)+{\it\epsilon}\,\text{Re}\{\hat{q}(x,t)\}+O({\it\epsilon}^{2}),\end{eqnarray}$$

with ${\it\epsilon}\ll 1$ . The mass conservation equation (2.11) becomes

(5.5) $$\begin{eqnarray}{\hat{h}}_{t}+\hat{q}_{x}=0.\end{eqnarray}$$

At $O({\it\epsilon})$ , the Benney flux equation (2.12) yields

(5.6) $$\begin{eqnarray}\hat{q}=b_{0}(x){\hat{h}}+b_{1}(x){\hat{h}}_{x}+b_{2}(x){\hat{h}}_{xxx},\end{eqnarray}$$

where

(5.7) $$\begin{eqnarray}\displaystyle & \displaystyle b_{0}=H^{2}\left(2-2H_{x}\cot {\it\theta}+\frac{H_{xxx}}{C}\right)+\frac{16R}{5}H^{5}H_{x}-\frac{8R}{3}H^{4}F, & \displaystyle\end{eqnarray}$$

(5.8a,b ) $$\begin{eqnarray}\displaystyle & \displaystyle b_{1}=-\frac{2}{3}H^{3}\cot {\it\theta}+\frac{8}{15}RH^{6},\quad b_{2}=\frac{H^{3}}{3C}. & \displaystyle\end{eqnarray}$$

The weighted-residual equivalent of (5.6) additionally involves a term proportional to $\hat{q}_{t}$ . Given a base solution $H(x)$ , $Q(x)$ and $F(x)$ , the coefficients $b_{0}$ , $b_{1}$ and $b_{2}$ are all known periodic functions of $x$ , with the same period as the base solution.

We now invoke the Floquet–Bloch form, observing that as (5.5) and (5.6) are linear equations with coefficients that are periodic in $x$ with period $L$ , the eigenfunctions can be written as

(5.9a,b ) $$\begin{eqnarray}{\hat{h}}(x,t)=\text{e}^{{\it\lambda}t+\text{i}kx}{\hat{h}}_{k}(x),\quad \hat{q}(x,t)=\text{e}^{{\it\lambda}t+\text{i}kx}\hat{q}_{k}(x),\end{eqnarray}$$

where the Floquet wavenumber $k\neq 2{\rm\pi}/L$ is real, and ${\hat{h}}_{k}$ and $\hat{q}_{k}$ are periodic functions of $x$ with period $L$ . Setting $k=0$ recovers perturbations of wavelength $L$ , while very small but non-zero $k$ corresponds to very long-wave perturbations. The solution is stable to perturbations with period $L$ if $\text{Re}({\it\lambda})<0$ for all eigenfunctions when $k=0$ . The base solution is linearly stable to perturbations of all wavelengths if $\text{Re}({\it\lambda})<0$ for all eigenfunctions for each real $k\in [0,{\rm\pi}/L]$ .

5.3. Effect of small-amplitude blowing and suction on stability

We now analyse the effect of small amplitude blowing and suction in the form $F=A\cos mx$ on the stability of eigenmodes with underlying Floquet wavenumber $k$ , which will allow us to determine how such forcing affects the critical Reynolds number. We do so by expanding the eigenvalue ${\it\lambda}$ for $R$ close to $R_{H}$ , and for small  $A$ :

(5.10) $$\begin{eqnarray}{\it\lambda}={\it\lambda}_{0}+(R-R_{H})\frac{\partial {\it\lambda}}{\partial R}+ZA^{2}+O((R-R_{H})^{2})+O((R-R_{H})A^{2})+O(A^{4}).\end{eqnarray}$$

We are concerned with the eigenvalue equal to ${\it\lambda}_{0}=-2\text{i}k$ at $R=R_{H}$ , which is a root of both the Benney and weighted-residual characteristic equations (5.1) and (5.2). We can evaluate the term $\partial {\it\lambda}/\partial R$ which appears in (5.10) by differentiating (5.1) or (5.2), as appropriate. The eigenvalue is even in $A$ because the transformation $A\rightarrow -A$ can be recovered by translation in $x$ by a distance $L/2$ , and the original eigenmode $\exp (\text{i}kx+{\it\lambda}t)$ has no preferred position; thus the leading-order contribution for small $A$ is $O(A^{2})$ .

The effect of the imposed suction on the linear stability is encapsulated in the coefficient $Z$ appearing in (5.10), which depends on the suction wavenumber $m$ and the perturbation wavenumber $k$ as well as $C$ and ${\it\theta}$ . In order to calculate $Z$ , we need to calculate both the base solution and the eigenfunction to $O(A^{2})$ . The calculation of the base solution to $O(A^{2})$ was discussed earlier in the context of steady states, and the required expansion is given by (4.5) and (4.6).

For small-amplitude steady solutions, the suction function $F(x)$ and the base solution $H(x)$ , $Q(x)$ are all periodic with wavenumber $m$ , and so the coefficients in the linearised equations (5.5) and (5.6) are also periodic with wavenumber $m$ . As a result, all eigenfunctions of (5.5) and (5.6) can be written in the Floquet form given by (5.9). The unknown functions ${\hat{h}}_{k}$ and $\hat{q}_{k}$ are periodic with wavenumber $m$ , are constant when $A=0$ and can be expanded for small $A$ as

(5.11) $$\begin{eqnarray}\displaystyle {\hat{h}}_{k}(x)=1+A(C_{1}\text{e}^{\text{i}mx}+C_{2}+C_{3}\text{e}^{-\text{i}mx})+A^{2}(D_{1}\text{e}^{2\text{i}mx}+D_{2}\text{e}^{\text{i}mx}+D_{3}+D_{4}\text{e}^{-\text{i}mx}+D_{5}\text{e}^{-2\text{i}mx}) & & \displaystyle \nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$

and

(5.12) $$\begin{eqnarray}\displaystyle \hat{q}_{k}(x) & = & \displaystyle \frac{\text{i}{\it\lambda}_{0}}{k}+A(E_{1}\text{e}^{\text{i}mx}+E_{2}+E_{3}\text{e}^{-\text{i}mx})\nonumber\\ \displaystyle & & \displaystyle +\,A^{2}(G_{1}\text{e}^{2\text{i}mx}+G_{2}\text{e}^{\text{i}mx}+G_{3}+G_{4}\text{e}^{-\text{i}mx}+G_{5}\text{e}^{-2\text{i}mx}),\end{eqnarray}$$

where the complex constants $C_{1}$ , $C_{2}$ , $C_{3}$ , $D_{1}$ , $D_{2}$ , $D_{3}$ , $D_{4}$ , $D_{5}$ , $E_{1}$ , $E_{2}$ , $E_{3}$ , $G_{1}$ , $G_{2}$ , $G_{3}$ , $G_{4}$ and $G_{5}$ are to be found along with  $Z$ . We must also choose a normalisation condition on the amplitude and phase of the eigenvector, and for this we use the condition

(5.13) $$\begin{eqnarray}{\hat{h}}_{k}(0)=1.\end{eqnarray}$$

To determine the unknown constants, we set $R=R_{H}$ , where $R_{H}$ is defined in (5.3) and is the critical Reynolds number for the stability of a flat film to perturbations of wavenumber $k$ , solve for the steady state given by (4.5) and (4.6) and then substitute the eigenvalue expansion from (5.10) and the eigenfunction from (5.9), (5.11), (5.12) into the two linearised equations (5.5) and (5.6), and also use the normalisation condition (5.13). We solve the linearised equations first at $O(A)$ and then at $O(A^{2})$ , at each order obtaining linear systems for the unknown coefficients, including $Z$ . We perform this calculation in Maple, and obtain a lengthy expression for $Z$ as a function of $m$ , $k$ , $C$ and ${\it\theta}$ ; the full result depends on whether we have used the weighted residual or Benney equations. Once $Z$ is known, we can obtain the neutral stability curve for fixed $k$ , defined by the condition $\text{Re}\{{\it\lambda}\}=0$ , so that

(5.14a,b ) $$\begin{eqnarray}R\sim R_{H}+A^{2}W,\quad W\equiv -\text{Re}\{Z\}\left(\text{Re}\left\{{\displaystyle \frac{\partial {\it\lambda}}{\partial R}}\right\}\right)^{-1}.\end{eqnarray}$$

Figure 7. Effect of small $A$ blowing and suction with wavelength $L$ on the critical Reynolds number for perturbations of underlying wavenumber $k$ as measured through the weighted-residual results for $W$ , with $C=0.05$ . The filled circles on each curve mark the point where $k=2{\rm\pi}/L$ . The dashed lines show the long-wave asymptotic estimate (5.16). (a)  ${\it\theta}={\rm\pi}/4$ , (b)  ${\it\theta}={\rm\pi}/2$ , (c)  ${\it\theta}=3{\rm\pi}/4$ .

Figure 7 shows the effect of small $A$ on the stability boundary as quantified by the weighted-residual results for $W$ . For ${\it\theta}={\rm\pi}/4$ , we find that forcing via blowing and suction at wavelength $L\equiv 2{\rm\pi}/m=10$ destabilises short wavelength perturbations, but has little effect on long wavelength perturbations. However, if the forcing wavelength $L\geqslant 20$ , $W$ is positive for all $k$ and so forcing has a stabilising effect on all perturbations. Furthermore, for long wavelength forcing, the magnitude of $W$ increases with $L$ , and the minimum value of $W$ occurs at $k=0$ ; this value is positive and increases rapidly with $L$ , so that the stabilising effect is increased for longer wavelength blowing and suction. For ${\it\theta}={\rm\pi}/2$ , forcing with wavelength $L$ has a destabilising effect on short waves for $L=10$ , and a stabilising effect for $L=20,$ 40, 80. However, $W$ tends to approximately $-6$ as $k\rightarrow 0$ , and so the imposed forcing in fact slightly destabilises long-wave perturbations. For ${\it\theta}=3{\rm\pi}/4$ , the behaviour for small perturbation wavenumber is reversed from the ${\it\theta}={\rm\pi}/4$ case; long-wave forcing is always destabilising for $k=0$ , and becomes increasingly destabilising as $L$ is increased.

Figure 8. The effect of small amplitude blowing and suction with wavenumber $m$ on perturbations with $k=0$ and with $k=m$ in the Benney and weighted-residual (WR) models, as described by the quantity $W$ defined in (5.14); $W>0$ in the shaded regions, so small-amplitude blowing and suction increases the critical Reynolds number, and has a stabilising effect on the flow. Results for ${\it\theta}={\rm\pi}/4$ , ${\rm\pi}/2$ , $3{\rm\pi}/4$ are shown in rows (a–d), (e–h) and (i–l) respectively. (a,e,i) Benney, $k=0$ ; (b, f, j) WR, $k=0$ ; (c,g,k) Benney, $k=m$ ; (d,h,l) WR, $k=m$ .

In order to gauge the behaviour of $W$ over a wider range of parameters, in figure 8, we survey the sign of $W$ for $k=0$ and $k=m$ , with $m$ and $C$ varying, for three values of ${\it\theta}$ , for both the Benney and weighted-residual models. The two models concur that for each inclination angle ${\it\theta}$ , there are some regions of $(C,m)$ parameter space where the imposition of blowing and suction has a stabilising effect on the uniform state, and other regions where such forcing destabilises the flow. For ${\it\theta}={\rm\pi}/4$ , both models predict that imposing long-wave blowing and suction (i.e.  $m\rightarrow 0$ ) is stabilising to perturbations with wavenumber $k\ll 1$ and also with wavenumber  $m$ . For ${\it\theta}={\rm\pi}/2$ , blowing and suction is destabilising in the limit $k\rightarrow 0$ for all $C$ and $m$ , but has a stabilising effect on perturbations with wavenumber $m$ if $m$ is not too large. For ${\it\theta}=3{\rm\pi}/4$ , blowing and suction with small wavenumber $m$ is destabilising to both perturbations with wavenumber 0 and with wavenumber  $m$ .

The critical Reynolds number increases quadratically with perturbation wavenumber $k$ in the absence of blowing and suction, and so for very small $A$ , the lowest critical Reynolds number must also be obtained at small $k$ . We now consider $W$ in the limit of both long-wavelength blowing and suction (small  $m$ ) and long perturbation wavelength (small  $k$ ), setting $m={\it\delta}M$ , $k={\it\delta}K$ and $C={\it\delta}^{2}\widehat{C}$ and expanding the full expression for $Z$ with ${\it\delta}\ll 1$ . Both the Benney and weighted-residual results yield

(5.15a,b ) $$\begin{eqnarray}Z\sim \frac{3\text{i}K}{4{\it\delta}M^{2}}+\left[-\frac{K^{2}\cot {\it\theta}}{2M^{2}}+\frac{K^{2}}{6\widehat{C}}-\frac{11K^{4}}{12\widehat{C}M^{2}}\right]+O({\it\delta}^{2}),\quad \text{Re}\left\{\frac{\partial {\it\lambda}}{\partial R}\right\}\sim \frac{8K^{2}{\it\delta}^{2}}{15}+O({\it\delta}^{4}).\end{eqnarray}$$

Substituting (5.15) into (5.14) yields the leading-order long-wave expansion

(5.16) $$\begin{eqnarray}W=-\frac{15}{8{\it\delta}^{2}}\left[-\frac{\cot {\it\theta}}{2M^{2}}+\frac{1}{6\widehat{C}}-\frac{11K^{2}}{12\widehat{C}M^{2}}\right]+O(1)=-\frac{15}{8}\left[-\frac{\cot {\it\theta}}{2m^{2}}+\frac{1}{6C}-\frac{11k^{2}}{12Cm^{2}}\right]+O(1).\end{eqnarray}$$

The prediction (5.16) is plotted in figure 7 (dashed lines) for comparison to the non-long-wave results obtained by direct evaluation of the Maple expressions for $Z$ and $\partial {\it\lambda}/\partial R$ .

For small-amplitude blowing and suction at fixed wavenumber $m$ , the critical Reynolds number for instability to perturbations of all real wavenumbers $k$ is

(5.17) $$\begin{eqnarray}R^{\ast }(m)=\min _{k\in \mathbb{R}}\left[\frac{5}{4}\cot {\it\theta}+\frac{5k^{2}}{8C}+A^{2}W(m,k,C,{\it\theta})\right].\end{eqnarray}$$

Using the long-wave expansion for $W$ given by (5.16), we obtain

(5.18) $$\begin{eqnarray}R^{\ast }(m)=\frac{5}{4}\cot {\it\theta}+\frac{15A^{2}}{16m^{2}}\cot {\it\theta}-\frac{5A^{2}}{16C}+\min _{k\in \mathbb{R}}\left[k^{2}\left(\frac{5}{8C}+\frac{55A^{2}}{32m^{2}C}\right)\right].\end{eqnarray}$$

The term in square brackets is positive, and so the minimum of the expression occurs at $k=0$ . This behaviour is driven by the behaviour of a flat film, where long-wave perturbations are the first to become unstable, and so have the smallest critical Reynolds number. Introducing blowing and suction will modify the eigenvalue spectrum, but the minimum critical Reynolds number must still occur for long-wave perturbations if $A$ is sufficiently small.

Proceeding with the long-wave approximation, we evaluate the minimum of (5.18) as

(5.19) $$\begin{eqnarray}R^{\ast }\sim \frac{5}{4}\cot {\it\theta}+\frac{15}{16}A^{2}\left(\frac{\cot {\it\theta}}{m^{2}}-\frac{1}{3C}\right),\end{eqnarray}$$

while if perturbations are restricted to those with wavenumber $k=m$ , we find

(5.20) $$\begin{eqnarray}R_{m=k}\sim \frac{5}{4}\cot {\it\theta}+\frac{5m^{2}}{8C}+\frac{15}{16}A^{2}\left(\frac{\cot {\it\theta}}{m^{2}}-\frac{3}{2C}\right).\end{eqnarray}$$

For comparison, we also calculate the linear stability properties of the full nonlinear steady states by numerical solution of the eigenvalue problem described in § 5.2, using Floquet multipliers in the form (5.9) to account for general perturbation wavenumber  $k$ . It is also possible to formulate an augmented boundary-value problem for the nonlinear base solution, linear eigenmode and the eigenvalue, and thus to track the neutral stability boundaries directly using Auto-07p.

The stability boundaries for perturbations with arbitrary $k$ and with $k=m$ are plotted in figure 9 for nonlinear solutions to the weighted-residual equations, small $A$ predictions obtained using the full version of $W$ and the long-wave, small $A$ predictions (5.19) and (5.20). We see that the first two boundaries are in good agreement with each other when $A$ is sufficiently small, for both $L=40$ and $L=10$ . The long-wave small-amplitude predictions are in good agreement with the other two versions of boundaries when $L=40$ , but are poor agreement when $L=10$ .

In the case ${\it\theta}={\rm\pi}/4$ , $R^{\ast }$ is positive in the absence of suction, and long-wave blowing and suction can either decrease or increase $R^{\ast }$ , depending on the values of $A$ , ${\it\theta}$ , $m$ and $C$ , and here it is reasonable to seek an optimal strategy to increase the critical Reynolds number. To maximise the stabilising effect at fixed $A$ for ${\it\theta}<{\rm\pi}/2$ , we should choose $m$ as small as possible, and in fact (5.19) predicts that $R^{\ast }\rightarrow +\infty$ as $m\rightarrow 0$ with $A$ fixed. Regarding constraints on the magnitude of film height variations, we note that these are of order $A/m$ for $A\ll 1$ (see § 4.1) , and so we can approximately constrain height variations by keeping ${\it\alpha}=A/m$ fixed. We find that the largest $R^{\ast }$ is again obtained as $m\rightarrow 0$ , but in contrast to the fixed $A$ case, $R^{\ast }$ is bounded as $m\rightarrow 0$ .

The governing equations are valid for $0<{\it\theta}<{\rm\pi}$ , with ${\it\theta}={\rm\pi}/2$ corresponding to flow down a vertical wall, and ${\rm\pi}/2<{\it\theta}<{\rm\pi}$ corresponding to flow along the underside of an inclined plane. In the absence of blowing and suction, $R^{\ast }=0$ for ${\it\theta}={\rm\pi}/2$ , and $R^{\ast }<0$ for ${\rm\pi}/2<{\it\theta}<{\rm\pi}$ . The definition of the Reynolds number $R$ in (2.4) means that we cannot construct a physical system with $R<0$ , but consideration of (5.19) may still give some indication of the effect of blowing and suction when ${\it\theta}>{\rm\pi}/2$ . According to (5.19), when ${\it\theta}\geqslant {\rm\pi}/2$ , imposing blowing and suction at long wavelengths always reduces $R^{\ast }$ , and so is destabilising. Our other results for ${\it\theta}=3{\rm\pi}/4$ , shown in figures 7 and 8, also suggest that blowing and suction have a destabilising effect on long-wave perturbations.

Figure 9. Stability properties for steady solutions as $R$ and $A$ vary at fixed ${\it\theta}={\rm\pi}/4$ and $C=0.05$ , for the weighted-residual model. Steady solutions are divided into three linear-stability categories as indicated by the legend, with finite- $A$ stability boundaries shown by solid black lines. Two small- $A$ estimates for these boundaries are shown: the full result from (5.14) correct to $O(A^{2})$ (solid white line), and the long-wave, small $A$ estimates (5.19) and (5.20) (dash-dotted line). We also indicate solutions with $\min (q)=0$ (dashed line), and the pitchfork bifurcation $A_{P}$ (dotted line). (a)  $L=10$ , (b)  $L=40$ .

5.4. Finite $A$ stability results

Figure 9 shows stability regions as $R$ and $A$ vary, at fixed values of $C$ and ${\it\theta}$ , for imposed blowing and suction with wavelength $L=10$ and $L=40$ . For $L=10$ , introducing small amplitude blowing and suction decreases the critical $R$ for linear stability both to perturbations of wavelength $L$ and to perturbations of arbitrary wavelength. In contrast, when the imposed suction has wavelength $L=40$ , increasing $A$ from zero initially increases both of these critical Reynolds numbers, and so has a stabilising effect on the base flow.

The boundary of the region stable to perturbations of fixed wavelength $L$ is defined by a Hopf bifurcation. When $A=0$ , this is the Hopf bifurcation at $R=R_{H}$ , for which the eigenmode can be written as ${\hat{h}}=\exp (\text{i}kx+{\it\lambda}t)$ with $k=2{\rm\pi}/L$ . We can use Auto-07p to track the Hopf bifurcation as $A$ varies. As $A$ is increased from zero, the eigenmode is modulated and is no longer a single Fourier mode, but remains periodic with period  $L$ .

Calculation of the stability boundary in the presence of patterned blowing and suction for perturbations of all wavelengths is not as simple as the case of tracking a single bifurcation, as we do not know a priori which wavelengths are most unstable. However, in the absence of blowing and suction, we know that long waves, with $k^{2}\ll 1$ , are the first to become unstable as $R$ is increased. For $R>R_{H}$ , there is a wavenumber cutoff, so that perturbations with $k^{2}<k_{c}^{2}$ are unstable, and there is a finite $k^{2}$ with maximum growth rate $\text{Re}({\it\lambda})$ . However, it is possible that for larger $A$ , imposing blowing and suction might destabilise at some wavelengths while stabilising at others, and so the system may first become unstable at a non-zero wavenumber. For finite $A$ , we conduct a brute force computation of stability properties over a large number of perturbation wavelengths to determine stability with respect to perturbations of all wavelengths; this involves solving (2.11) and either (2.12) or (2.13) to obtain the nonlinear steady state, and then solving the eigenvalue problem for linear stability described in § 5.2 with the perturbation wavenumber incorporated in the form of a Floquet multiplier in (5.9).

For the cases shown in figure 9, we find that at small $A$ , the boundary for stability to perturbations of arbitrary wavelength agrees well with the asymptotic results derived in the previous subsection. For both $L=10$ and $L=40$ , this boundary connects smoothly to a finite suction amplitude $A$ at $R=0$ . There is a turning point with respect to $R$ in the $L=40$ results, and so the largest $R$ at which the system is stable occurs for a non-zero $A$ . However, for $L=10$ , there is no turning point, but there is a second ‘island’ region where the steady state is stable to perturbations of all wavelengths. This region requires moderately large $A$ and $R$ . We will discuss the results of time-dependent simulations in and around the island in § 7. As the island region occurs only for $L=10$ , it is possible that it is simply an artefact of applying long-wave models at relatively short wavelengths, and we have checked that no such island region arises in similar calculations based on the Benney equation.

5.5. Optimal wavelength for blowing and suction to obtain a stable steady state

Determination of the largest $R$ that can be stabilised in an infinite domain by imposing a suitable suction profile is a complicated process, requiring finite amplitude results for the steady solutions and for their linear-stability operators. We can address this task numerically by choosing a suction wavenumber $m$ , increasing $R$ in small increments from $R_{0}=(5/4)\cot {\it\theta}$ and testing the stability of steady solutions beginning from $A=0$ until $A$ is too large for steady states to exist. Perturbations with wavenumber $k$ not equal to the wavenumber of the forcing $m$ were taken into account in the linear-stability calculations using Floquet multipliers. The smallest perturbation wavenumber $k$ tested was 0.005. Figure 10 shows numerical results from this search for the case ${\it\theta}={\rm\pi}/4$ , $C=0.05$ , using the weighted-residual model. We find two ranges of $m$ in which imposing suction can increase the critical Reynolds number.

Figure 10. Finite $A$ stability results across a range of blowing and suction wavenumber $m$ , with $F=A\cos mx$ , $C=0.05$ and ${\it\theta}={\rm\pi}/4$ . Each black dot indicates that for the given $R$ and $m$ , there is some $A$ for which there is a steady state stable to linear perturbations of all wavelengths. The dashed line is obtained by explicit tracking of the maximum stable $R$ for fixed Floquet number $k=0.01$ ; this curve follows the stability boundary relatively well until it diverges at $m\approx 0.05$ . In the small-amplitude long-wave limit, fixed ${\it\alpha}=A/m$ corresponds to fixing the amplitude of interface deflection as $m$ varies. We find that along the dashed line, ${\it\alpha}$ is between 0.5 and 1. To provide an indication of the effectiveness of the small- $A$ , long-wave stability prediction (5.19), the solid blue line shows the stability boundary that would be obtained if ${\it\alpha}=1$ , which is an overestimate of $A$ , but in fact the blue line under predicts the range of stable $R$ for small  $m$ .

Firstly, at large wavelengths with $L\approx 200$ , imposing blowing and suction allows the critical Reynolds number for stability to all perturbations to increase from the unforced value of $5/4$ to approximately 8. However, figure 10 seems to indicate a downturn in the maximum stable $R$ at very small  $m$ . From the small  $A$ , long-wave result (5.19), we on the contrary expect forcing via blowing and suction at small $m$ to have some stabilising effect, but this result does not predict the magnitude of the stabilising effect on its own, as we do not know  $A$ . Furthermore, estimates of the maximum stable $R$ typically require large- $A$ calculations, and so predictions based on the long-wave, small-amplitude expression (5.19) are not particularly useful in this case.

The second region in which steady solutions are stable to perturbations of all wavelengths in figure 10 appears when the wavenumber $m$ for blowing and suction is relatively large, and only at moderately large  $R$ . This region corresponds to the island region in the $L=10$ results in figure 9. There is no such stable island in the results for $L=40$ , or for the equivalent Benney results for $L=10$ , and so it is not clear that this island is of physical relevance.

6. Travelling waves

6.1. Travelling waves in the absence of suction

In the absence of blowing and suction, the system can support travelling waves, propagating at a constant speed $U$ without changing form. We can write the solution as

(6.1a−c ) $$\begin{eqnarray}h(x,t)=H({\it\zeta}),\quad q(x,t)=Q({\it\zeta}),\quad {\it\zeta}=x-Ut,\end{eqnarray}$$

where $U$ is the unknown constant wave speed. In the weighted-residual equations, $H$ and $Q$ must satisfy

(6.2) $$\begin{eqnarray}-UH^{\prime }+Q^{\prime }=0\end{eqnarray}$$

and

(6.3) $$\begin{eqnarray}-\frac{2}{5}RUHQ^{\prime }+Q=\frac{H^{3}}{3}\left(2-2H^{\prime }\cot {\it\theta}+\frac{H^{\prime \prime \prime }}{C}\right)+R\left(\frac{18}{35}Q^{2}H^{\prime }-\frac{34}{35}QQ^{\prime }H\right),\end{eqnarray}$$

where a prime indicates a derivative with respect to ${\it\zeta}$ . If $H({\it\zeta})$ is periodic, with period $L$ , then the travelling-wave solution is spatially periodic with period $L$ , and temporally periodic with period $T=L/U$ . We can compute large-amplitude travelling waves numerically; figure 11 shows one such travelling-wave solution for ${\it\theta}={\rm\pi}/4$ , $C=0.05$ , $R=3$ .

Figure 11. Travelling wave $H({\it\zeta})$ and perturbation functions $J({\it\zeta})$ and $K({\it\zeta})$ as defined in (6.27) for ${\it\theta}={\rm\pi}/4$ , $C=0.05$ , $R=3$ using the weighted-residual equations, and corresponding superpositions of the travelling wave $H({\it\zeta})$ and $O(A)$ perturbation ${\hat{h}}(x,t)$ as defined in (6.10) for various $A$ . $H({\it\zeta})$ has wavelength 30, but blowing and suction are applied with wavelength 20; we show solutions in a domain of length 60. The black contours indicate $h=1$ , and the same colour map is used in the first three figures. Both travelling wave and perturbation are periodic in time and space. Nonlinear time-dependent results for $A=0.02$ are shown in figure 13 for various initial conditions.

Individual travelling waves may be stable or unstable, and branches of travelling waves can undergo a range of bifurcations. However, small amplitude travelling waves connect to the uniform film state via a Hopf bifurcation at $R=R_{H}$ . This Hopf bifurcation can be supercritical or subcritical, with stable, small-amplitude travelling waves observed near the bifurcation only in the supercritical case. We can determine the criticality of the Hopf-bifurcation by solving for small-amplitude limit cycles near to the critical Reynolds number via the following expansion:

(6.4a,b ) $$\begin{eqnarray}{\it\zeta}=x-Ut,\quad H({\it\zeta})=1+{\it\epsilon}\cos k{\it\zeta}+{\it\epsilon}^{2}H_{2}({\it\zeta})+O({\it\epsilon}^{3}),\end{eqnarray}$$

(6.5a,b ) $$\begin{eqnarray}H_{2}({\it\zeta})=r_{1}\cos 2k{\it\zeta}+r_{2}\sin 2k{\it\zeta},\quad q(x,t)=UH({\it\zeta})+S,\end{eqnarray}$$

(6.6a−c ) $$\begin{eqnarray}U=U_{0}+{\it\epsilon}^{2}U_{2}+O({\it\epsilon}^{4}),\quad S=S_{0}+{\it\epsilon}^{2}S_{2}+O({\it\epsilon}^{4}),\quad R=R_{H}+{\it\epsilon}^{2}\bar{R}+O({\it\epsilon}^{4}).\end{eqnarray}$$

We expand the equations for ${\it\epsilon}\ll 1$ , and must solve the equations at up to $O({\it\epsilon}^{3})$ in order to determine $\bar{R}$ . We find that small-amplitude travelling waves always travel downstream, with speed $U=U_{0}=2$ , which is twice the velocity of particles on the surface. The bifurcation is supercritical if $\bar{R}>0$ and subcritical if $\bar{R}<0$ . The Benney equations yield

(6.7) $$\begin{eqnarray}\bar{R}=\frac{120C^{2}-60C^{2}k^{2}\cot ^{2}{\it\theta}-120Ck^{4}\cot {\it\theta}-45k^{6}}{64Ck^{4}},\end{eqnarray}$$

which may be positive or negative. However, for the weighted-residual equations, we find

(6.8) $$\begin{eqnarray}\bar{R}=\frac{4410C^{2}+6670k^{2}C^{2}\cot ^{2}{\it\theta}+12\,235k^{4}C\cot {\it\theta}+4450k^{6}}{2352Ck^{4}},\end{eqnarray}$$

which is always positive when ${\it\theta}<{\rm\pi}/2$ , corresponding to a supercritical Hopf bifurcation. In the long-wave limit, with $k={\it\delta}K$ and $C={\it\delta}^{2}\widehat{C}$ , both (6.7) and (6.8) yield

(6.9) $$\begin{eqnarray}\bar{R}\sim \frac{15\widehat{C}}{8{\it\delta}^{2}K^{4}}+O({\it\delta}^{2})=\frac{15C}{8k^{4}}+O({\it\delta}^{2}),\end{eqnarray}$$

which is positive.

The quantity $\bar{R}$ determines the criticality of the Hopf bifurcation, and so also governs the stability of small-amplitude travelling waves in the neighbourhood of the bifurcation. However, the branch of travelling waves may undergo further bifurcations (Scheid et al. 2005), and so the value of $\bar{R}$ does not necessarily prescribe the stability of finite-amplitude travelling waves.

6.2. Influence of heterogeneous blowing and suction on travelling waves

Introducing spatially-periodic suction means that the system is no longer translationally invariant, and so we cannot obtain true travelling waves for non-zero  $A$ . For the next stage of our analysis, we calculate the effect of small-amplitude suction on travelling waves, and consider in particular how travelling waves can transition to steady, stable, non-uniform states as the amplitudes of the imposed blowing and suction are increased.

We perturb the travelling wave by introducing a small-amplitude suction $F=Af(x)$ with $f(x)$ periodic with period $L_{F}$ , and $A$ small. We expect to find

(6.10a,b ) $$\begin{eqnarray}h(x,t)=H({\it\zeta})+A{\hat{h}}(x,t)+O(A^{2}),\quad q(x,t)=Q({\it\zeta})+A\hat{q}(x,t)+O(A^{2}),\end{eqnarray}$$

where the behaviour of ${\hat{h}}$ and $\hat{q}$ is in some way related to the periodicity of the original travelling wave with respect to ${\it\zeta}$ with period $L$ , and of the blowing and suction function with respect to $x$ with period  $L_{F}$ .

The weighted-residual equations for $h$ and $q$ yield, at $O(A)$ ,

(6.11) $$\begin{eqnarray}{\hat{h}}_{t}-f(x)+\hat{q}_{x}=0\end{eqnarray}$$

and

(6.12) $$\begin{eqnarray}\hat{q}+w_{6}({\it\zeta})\hat{q}_{t}=w_{0}({\it\zeta}){\hat{h}}+w_{1}({\it\zeta}){\hat{h}}_{x}+w_{2}({\it\zeta}){\hat{h}}_{xxx}+w_{3}({\it\zeta})\hat{q}+w_{4}({\it\zeta})\hat{q}_{x}+w_{5}({\it\zeta})f(x),\end{eqnarray}$$

where

(6.13) $$\begin{eqnarray}w_{0}=H^{2}\left(2-2H^{\prime }\cot {\it\theta}+\frac{H^{\prime \prime \prime }}{C}\right)-\frac{34RQQ^{\prime }}{35},\end{eqnarray}$$

(6.14a−c ) $$\begin{eqnarray}w_{1}=-\frac{2}{3}H^{3}\cot {\it\theta}+\frac{18RQ^{2}}{35},\quad w_{2}=\frac{H^{3}}{3C},\quad w_{3}=R\left(\frac{36QH^{\prime }}{35}-\frac{34HQ^{\prime }}{35}\right),\end{eqnarray}$$

(6.15a−c ) $$\begin{eqnarray}w_{4}=-\frac{34RHQ}{35},\quad w_{5}=\frac{RHQ}{5},\quad w_{6}=\frac{2RH^{2}}{5}.\end{eqnarray}$$

The coefficients $w_{i}$ , $i=0,\ldots ,6$ are functions of ${\it\zeta}$ . If we set $f(x)=0$ in these equations, we recover the equations governing the evolution of linearised perturbations to the travelling wave in the absence of suction, i.e. the equations of linear stability. For non-zero $f$ , we can integrate forward in time from a given initial condition ${\hat{h}}(x,0)={\hat{h}}_{0}(x)$ , $\hat{q}(x,0)=\hat{q}_{0}(x)$ , to determine ${\hat{h}}(x,t)$ and $\hat{q}(x,t)$ .

The system (6.11) and (6.12) features some terms that are known functions of ${\it\zeta}$ , and others that are known functions of  $x$ , but the system is autonomous with respect to  $t$ . We therefore choose to rewrite (6.11) and (6.12) in terms of $x$ and  ${\it\zeta}$ :

(6.16a−d ) $$\begin{eqnarray}{\hat{h}}_{t}\rightarrow -U{\hat{h}}_{{\it\zeta}},\quad {\hat{h}}_{x}\rightarrow {\hat{h}}_{x}+{\hat{h}}_{{\it\zeta}},\quad \hat{q}_{t}\rightarrow -U\hat{q}_{{\it\zeta}},\quad \hat{q}_{x}\rightarrow \hat{q}_{x}+\hat{q}_{{\it\zeta}}.\end{eqnarray}$$

We obtain a system in two variables, $x$ and ${\it\zeta}$ , with equations

(6.17) $$\begin{eqnarray}-U{\hat{h}}_{{\it\zeta}}-f(x)+\hat{q}_{x}+\hat{q}_{{\it\zeta}}=0\end{eqnarray}$$

and

(6.18) $$\begin{eqnarray}\displaystyle \hat{q}-Uw_{6}({\it\zeta})\hat{q}_{{\it\zeta}} & = & \displaystyle w_{0}({\it\zeta}){\hat{h}}+w_{1}({\it\zeta})({\hat{h}}_{x}+{\hat{h}}_{{\it\zeta}})+w_{2}({\it\zeta})({\hat{h}}_{xxx}+3{\hat{h}}_{xx{\it\zeta}}+3{\hat{h}}_{x{\it\zeta}{\it\zeta}}+{\hat{h}}_{{\it\zeta}{\it\zeta}{\it\zeta}})\nonumber\\ \displaystyle & & \displaystyle +\,w_{3}({\it\zeta})\hat{q}+w_{4}({\it\zeta})(\hat{q}_{x}+\hat{q}_{{\it\zeta}})+w_{5}({\it\zeta})f(x).\end{eqnarray}$$

We now regard ${\it\zeta}$ and $x$ as independent variables. Under this transformation, $x$ remains as a purely spatial variable, but ${\it\zeta}$ has a dual identity, incorporating both spatial and time-like components. The statement of initial conditions becomes more complicated in the new variables:

(6.19a,b ) $$\begin{eqnarray}{\hat{h}}({\it\zeta}=x)=h_{0}(x),\quad \hat{q}({\it\zeta}=x)=q_{0}(x)\end{eqnarray}$$

so that the initial conditions are spread across the whole range of  ${\it\zeta}$ .

If the Fourier expansion of $f(x)$ is

(6.20) $$\begin{eqnarray}f(x)=\mathop{\sum }_{m}f_{m}\cos mx+g_{m}\sin mx,\end{eqnarray}$$

we can write the general solution of (6.17) and (6.18) as

(6.21) $$\begin{eqnarray}\displaystyle & \displaystyle {\hat{h}}=\mathop{\sum }_{m}J_{m}({\it\zeta})\cos mx+K_{m}({\it\zeta})\sin mx+h^{\ast }({\it\zeta},x), & \displaystyle\end{eqnarray}$$

(6.22) $$\begin{eqnarray}\displaystyle & \displaystyle \hat{q}=\mathop{\sum }_{m}M_{m}({\it\zeta})\cos mx+N_{m}({\it\zeta})\sin mx+q^{\ast }({\it\zeta},x), & \displaystyle\end{eqnarray}$$

where the functions $J_{m}({\it\zeta})$ , $K_{m}({\it\zeta})$ , $M_{m}({\it\zeta})$ and $N_{m}({\it\zeta})$ are periodic in ${\it\zeta}$ . The remaining terms, $h^{\ast }$ and $q^{\ast }$ , satisfy the homogeneous versions of (6.17) and (6.18) obtained by setting $f(x)=0$ , which are exactly the equations for linear stability of the original travelling wave. We can calculate the periodic functions $J_{m}$ , $K_{m},\ldots ,$ corresponding to a limit cycle, for any periodic travelling wave. However, we would only expect to observe the limit-cycle behaviour in initial-value calculations if the limit cycle is stable. If the travelling wave is stable, all solutions $h^{\ast }$ , $q^{\ast }$ to the homogeneous perturbation problem eventually decay to zero and so ${\hat{h}}$ and $\hat{q}$ are limit cycles, periodic in time.

We now calculate ${\hat{h}}$ for a small-amplitude limit cycle driven by $f(x)=\cos mx$ , so that the forcing Fourier series has only a single component. The sums in (6.21) are then over a single value of $m$ , for which we obtain

(6.23) $$\begin{eqnarray}\displaystyle & \displaystyle -UJ_{m}^{\prime }+mN_{m}+M_{m}^{\prime }=1, & \displaystyle\end{eqnarray}$$

(6.24) $$\begin{eqnarray}\displaystyle & \displaystyle -UK_{m}^{\prime }-mM_{m}+N_{m}^{\prime }=0, & \displaystyle\end{eqnarray}$$

(6.25) $$\begin{eqnarray}\displaystyle M_{m}-w_{6}UM_{m}^{\prime } & = & \displaystyle w_{0}J_{m}+w_{1}(J_{m}^{\prime }+mK_{m})+w_{2}(J_{m}^{\prime \prime \prime }+3mK_{m}^{\prime \prime }-3m^{2}J_{m}^{\prime }-m^{3}K_{m})\nonumber\\ \displaystyle & & \displaystyle +\,w_{3}M_{m}+w_{4}(M_{m}^{\prime }+mN_{m})+w_{5}\end{eqnarray}$$

and

(6.26) $$\begin{eqnarray}\displaystyle N_{m}-w_{6}UN_{m}^{\prime } & = & \displaystyle w_{0}K_{m}+w_{1}(K_{m}^{\prime }-mJ_{m})+w_{2}(J_{m}^{\prime \prime \prime }+3mK_{m}^{\prime \prime }-3m^{2}J_{m}^{\prime }-m^{3}K_{m})\nonumber\\ \displaystyle & & \displaystyle +\,w_{3}N_{m}+w_{4}(N_{m}^{\prime }-mM_{m}).\end{eqnarray}$$

We solve the equations for $M=M_{m}$ , $N=N_{m}$ , $J=J_{m}$ and $K=K_{m}$ in Auto-07p, coupled to a system to find the travelling wave itself, seeking both travelling wave and perturbation periodic in ${\it\zeta}$ with period  $L$ . Thus we obtain the solution

(6.27) $$\begin{eqnarray}{\hat{h}}=J_{m}({\it\zeta})\cos mx+K_{m}({\it\zeta})\sin mx=J_{m}(x-Ut)\cos mx+K_{m}(x-Ut)\sin mx.\end{eqnarray}$$

Regardless of the relation between the travelling-wave period $L$ and the blowing/suction wavenumber $m$ , we see that if $J_{m}({\it\zeta})$ and $K_{m}({\it\zeta})$ are periodic with period  $L$ , then the function ${\hat{h}}(x,t)$ is temporally periodic with period $T=U/L$ , which is the same temporal period as the base travelling wave. However, the solution for ${\hat{h}}(x,t)$ is spatially periodic only if $L/L_{F}$ is rational, with period given by the least common multiple of $L$ and $L_{F}$ . A typical set of solutions for $H({\it\zeta})$ , $J_{m}({\it\zeta})$ and $K_{m}({\it\zeta})$ , and also the reconstructed field ${\hat{h}}$ , are shown in figure 11. For the travelling-wave solution shown in figure 11, we find that the functions $J_{m}({\it\zeta})$ and $K_{m}({\it\zeta})$ display little relative variation with ${\it\zeta}$ . This means that the field ${\hat{h}}(x,t)$ is essentially steady.

Figure 12. Contour plots of $h(x,t)$ from nonlinear time-dependent calculations (a–e) and using the small- $A$ asymptotic solutions for perturbed travelling waves given by (6.10) and (6.27) (f–j) for $R=1.7$ , $C=0.05$ , ${\it\theta}={\rm\pi}/4$ , in a domain of length 60 with suction wavelength $L_{F}=20$ . The colour map is scaled to the maximum and minimum value in each column. Here the travelling wave is stable in the absence of suction. Time-dependent results for a larger  $R$ , for which the uniform film solution is unstable to multiple perturbations, are shown in figure 13. (a, f)  $A=0$ ; (b,g)  $A=0.02$ ; (c,h)  $A=0.04$ ; (d,i)  $A=0.10$ ; (e, j)  $A=0.18$ .

Figure 12 shows a comparison between fully nonlinear time-dependent calculations and the perturbed travelling-wave solution at a value of $R$ small enough that the uniform state is unstable to only a single perturbation, and the travelling wave is stable. We obtain good agreement between the time-dependent calculations and the asymptotic predictions based on (6.10) and (6.27). Both sets of results show a smooth transition as $A$ increases from travelling waves with wavelength 60 at $A=0$ to almost-steady states at $A=0.18$ with wavelength 20, which is the wavelength of the imposed blowing and suction.

The composite solution shown in figure 11(b) has $L_{F}=20$ and $L=30$ ; the resulting perturbed field has spatial period 60. We have only considered corrections up to $O(A)$ , and in this expansion, the imposed blowing and suction cannot affect the periodicity of the underlying travelling wave. However, in the fully nonlinear time-dependent simulation, even if we begin with initial conditions that are perfectly periodic with period $L=30$ , but force at a different wavelength, such as $L_{F}=20$ , nonlinear effects will lead to a perturbation at wavelength 60 that will cause the underlying travelling wave to double in spatial period, and likely change shape and speed. Eventually, we expect to reach the state where the dominant wave has period 60, and the linear perturbation field ${\hat{h}}$ is a solution to equations where the coefficients $w_{i}$ are those for the travelling wave with wavelength 60. Figure 13 shows time-dependent simulations for forcing at wavelength 20 with $A=0.02$ in a domain of length 60, at a value of $R$ large enough that travelling waves with wavelengths 60, 30 and 20 are unstable. When the initial conditions involve only modes of wavelength 60, a periodic equilibrium state is quickly reached. For initial conditions of wavelength 30, these modes compete with the wavelength 20 forcing, but eventually a wavelength 60 state is indeed achieved. For initial conditions containing only perturbations with the same wavelength as the imposed blowing and suction, we rapidly reach an initial periodic state with three equal waves in the domain. However, after a very long time, noise in the system causes a switch to a single wave with wavelength 60, thus yielding the same eventual behaviour regardless of initial conditions.

Figure 13. Time-dependent simulations in a domain of length 60, with $F=0.02\cos (2{\rm\pi}x/20)$ , with $R=3$ , $C=0.05$ and ${\it\theta}={\rm\pi}/4$ , and varying initial conditions: we set $h(x,0)=1+0.1\cos (2{\rm\pi}nx/60)$ and $q(x,0)=2/3+0.2\cos (2{\rm\pi}nx/60)$ , with $n=1,2,3$ in (a–c) respectively. These initial conditions correspond to neither travelling waves nor steady states, so nonlinear evolution is required if the system is to reach a periodic state. These plots shows the single contour $h(x,t)=1$ , with $h>1$ in the shaded region. The same periodic state is reached eventually, regardless of the initial conditions.

7. Initial value problems

As discussed in § 4, steady states only exist if the amplitude of blowing and suction is below a critical value. In time-dependent simulations above the critical amplitude, as illustrated in figure 4, the film exhibits a new form of blow up by locally thinning to zero in finite time. This thinning occurs even at zero Reynolds number, and so arises in both Benney and weighted-residual models. Similar blow up can occur even below the critical amplitude if the interface is not initially close to a stable steady state. The blow up need not occur at the same wavelength as the forcing function; for example it can be triggered by subharmonic perturbations to an unstable, periodic steady state.

The film dynamics still interacts with the imposed blowing and suction even if blow up is avoided. In the absence of suction, the system exhibits periodic behaviour in the form of travelling waves, which propagate at a constant speed without changing form. If non-uniform suction is imposed via a non-constant function $F(x)$ , which remains fixed with respect to the frame of the wall, waves must change in form as they propagate, and so propagating waves are unsteady in any frame, with a complicated, but possibly periodic, dependence on both $x$ and  $t$ . Figure 12 shows initial value calculations in which the unforced system exhibits stable travelling waves. The imposed suction need not have the same wavelength as the travelling wave, and in figure 12, three periods of suction fit within the discretised domain. The initial value calculations show a smooth transition as $A$ increases from travelling waves at $A=0$ to an essentially steady state at $A=0.18$ . At small  $A$ , the blowing and suction causes slight perturbations to the travelling wave, while at larger  $A$ , small-amplitude disturbances propagate over a large-amplitude non-uniform steady state. The underlying flow field retains its dominant down-slope direction, and so there is still a well-defined wave speed, even when the interface shape is steady.

Figure 14. Nonlinear time-dependent calculations subject to blowing and suction $F=A\cos 2{\rm\pi}x/L$ , with $L=10$ , $C=0.05$ , ${\it\theta}={\rm\pi}/4$ . The simulations are conducted in a domain of length $6L$ , and the initial conditions for $h$ and $q$ includes a small perturbation proportional to $\sin (2{\rm\pi}x/(6L))$ . The stable, steady solution shown in (b) can be accessed by either increasing the Reynolds number from (c), or the amplitude of blowing and suction from (a). Time-dependent waves persist at long times in both (a) and (c); in (a), the waves propagate more rapidly and have more than one temporal period; in (c) the waves propagate more slowly, and there is just one propagating wave. (a)  $A=0.25$ , $R=4$ ; (b)  $A=0.45$ , $R=4$ ; (c)  $A=0.45$ , $R=2$ .

The applied suction can also lead to states which are neither steady nor time periodic. Figure 14 shows the result of initial value calculations conducted at parameters near to the stable ‘island’ shown in figure 9 for $L=10$ . This is a region of parameter space with steady solutions that are stable to perturbations of all wavelengths, but is unusual in that solutions lose stability if either $A$ or $R$ is decreased, so that sufficiently large inertia is required to maintain stability. In the case shown in figure 14(a), with $A=0.25$ and $R=4$ , the height field displays aperiodic behaviour. At each instant, six peaks corresponding to the six periods of the suction function are visible, but waves propagate over these peaks without a clear structure. We can increase the suction amplitude to $A=0.45$ to obtain steady, stable solution, shown in figure 14(b). Decreasing the Reynolds number to $R=2$ means that the steady state again loses stability, but in this case to a wavelength 60 travelling wave which propagates over the steady solution with a single time period, as shown in figure 14(c).

8. Conclusion

In this paper, we considered thin-film flow down an inclined plane modified so that fluid is injected and withdrawn through the wall according to an arbitrary periodic function $F(x,t)$ . We derived and studied two long-wave models for this system, based on the first-order Benney and weighted-residual formulations. If the average layer height is conserved, then $F$ must have zero mean in space. We then specialised to the case where $F$ is a single steady Fourier mode, $F=A\cos mx$ , and investigated the form, bifurcations and linear stability of steady states, as well as a range of fully-nonlinear time-dependent behaviours.

Any steady states subject to non-uniform blowing and suction must themselves be non-uniform. When the amplitude of the blowing and suction is small, we can calculate the deformation of the film from a uniform state analytically. The interface has the same period as the applied blowing and suction, with a phase difference due to the down-slope direction of the base flow which tends to ${\rm\pi}/2$ as the wavelength of blowing and suction is increased. The predicted interface deformation differs between the two models, but agrees at long wavelength and at small Reynolds number. We find that the spatially-averaged down-slope flux can be either increased or decreased by the blowing and suction, and the effect is $O(A^{2})$ . However, in the long-wave limit, where the local flux is determined only by the local film height, the mean flux is always increased. This is because steady states subject to blowing and suction must be non-uniform, and the increased flux in regions of thicker film outweighs the decrease in flux in the thinner regions.

For larger amplitude blowing and suction, we calculate periodic steady states numerically by formulating a boundary-value problem for the interface shape, firstly assuming that the interface has the same spatial period as the blowing and suction. As the blowing and suction amplitude $A$ increases, the interface becomes increasingly deformed and the nonlinearity of the governing equations becomes increasingly important, leading to a non-trivial bifurcation structure including multiple steady states. There is a maximum amplitude of blowing and suction for steady states, which corresponds to a limit point in the bifurcation structure. Time-dependent evolution beyond the existence of steady states shows that the film can dry in finite time, as the imposed rate of fluid removal in regions of suction is locally not matched by the supply of fluid.

We find that there is a region of moderately large $A$ for which the steady solutions have regions with negative down-slope volume flux  $q$ , and we analysed the flow field within the fluid layer to interpret these observations. The imposition of non-zero, non-uniform $F$ leads to a series of isolated stagnation points along the wall, wherever $F(x)=0$ . Connections between the stagnation points allow us to demarcate the boundary between ‘recirculating’ fluid which enters and leaves through the wall, and ‘propagating’ fluid which never reaches the wall. The height of the recirculating layer near the wall increases with the amplitude of the blowing and suction. For large enough amplitude, all steady solutions have negative flux $q$ for some  $x$ . For those $x$ with $q(x)<0$ , all flow is directed up the slope, against gravity; this flow reversal means that the propagating layer disappears. The disappearance of the propagating layer could be of practical relevance if blowing and suction is used to enhance mixing or heat transfer in the fluid; blowing and suction could be used to halt the down-slope progression and increase recirculation, followed by reducing the blowing and suction in order to enable down-slope transport.

Steady solutions need not have the same spatial period as the imposed blowing and suction, and as a test case we seek steady solutions that are periodic in a domain containing two periods of the blowing and suction. This allows us to calculate subharmonic solution branches in addition to the harmonic bifurcation structure discussed above. We find that the two types of solution branches are connected by symmetry-breaking pitchfork bifurcations. In our calculations, we have not found any subharmonic steady states that are stable when subjected to time-dependent evolution. Unstable subharmonic steady states also occur in thin-film flow with periodic topography (Tseluiko et al. 2013), but it is possible that subharmonic steady states may be stable for other parameter values. Initial-value calculations starting from near the periodic steady states sometimes show subharmonic drying, where the film dries at any one of the local minima of  $h$ .

For falling liquid films, disturbances to steady states will tend to develop into waves which propagate in the down-slope direction. Blowing and suction will alter the propagation of these disturbances, not least because the basic state of the system becomes non-uniform in the direction of propagation. We calculated the linear stability of the non-uniform steady states numerically for two classes of perturbations: firstly, those that are periodic with the same period as the blowing and suction, and secondly perturbations of any wavelength. We found that blowing and suction can increase or decrease the critical Reynolds number for linear stability to both classes of perturbation, depending on the system parameters $R$ , $C$ and ${\it\theta}$ , as well as the amplitude and wavelength of the applied blowing and suction.

When the amplitude $A$ of blowing and suction is small, we can calculate the $O(A^{2})$ correction to the critical Reynolds number analytically for both classes of perturbations. For flow above an inclined plane ( ${\it\theta}<{\rm\pi}/2$ ), we find that imposing blowing and suction at very long wavelength increases the critical Reynolds number for both classes of perturbations. Due to the quadratic dependence of the critical $R$ on perturbation wavenumber, the first perturbations to become unstable as $R$ increases are those with infinite wavelength, both when $A=0$ and for small  $A$ . The converse case of flow underneath an inclined plane ( ${\it\theta}>{\rm\pi}/2$ ) is always unstable in the absence of blowing and suction, and in fact our small-amplitude expansions suggest that the flow is made more unstable by applying long-wave blowing and suction.

Even when small-amplitude blowing and suction has a stabilising effect relative to the uniform flow, large-amplitude suction will eventually lead to symmetry-breaking bifurcations or localised drying of the film. There is then, for fixed blowing and suction wavelength, an intermediate amplitude which has the largest stabilising effect on the film, and a corresponding maximum stable Reynolds number. Calculations using the weighted-residual model predict large increases in the critical Reynolds number, from $5/4$ to 4 or more, for stability to perturbations of arbitrary wavelength if blowing and suction is applied with wavelengths of around 200 times the mean thickness of the film.

In the absence of suction, and beyond the critical Reynolds number, the film can exhibit a range of behaviour including the emergence of travelling waves which propagate in the down-slope direction at constant speed $U$ without changing form. Imposing blowing and suction as $F(x)$ introduces heterogeneity to the system, and so we expect periodic travelling waves $h=H({\it\zeta})$ where ${\it\zeta}=x-Ut$ to become limit cycles that are periodic in both $x$ and $t$ , with explicit dependence on both variables. We calculated the $O(A)$ correction to the travelling wave due to small amplitude blowing and suction, showing that the equations are closely related to those governing linear stability of the travelling wave and can be solved simply if transformed to $(x,{\it\zeta})$ coordinates. The resulting solution for $h(x,t)$ and $q(x,t)$ is periodic in time with the same time period as the underlying travelling wave, but is spatially periodic only if the ratio between the wavelength of the blowing and suction and the wavelength of the travelling wave is rational. In many cases, the $O(A)$ correction to the travelling wave is almost steady in a fixed frame, and so the small $A$ expansion allows exploration of the transition between finite-amplitude travelling waves propagating through a slowly-varying environment when the imposed blowing and suction is small, and small waves propagating over a steady, non-uniform state for larger-amplitude blowing and suction. We also conducted initial-value calculations for different initial conditions to investigate competition between the period of the travelling wave and of the imposed suction. The same state is eventually reached regardless of initial conditions, but competition between different wavelengths can persist over a long time.

The Benney and weighted-residual models that we have derived here are also valid for the case of a time-dependent blowing and suction profile, so long as the timescale for variation is not too short. This leads to the possibility of delivering feedback control by basing the blowing and suction profile on observations of the interface height. Under the assumption of small deviations from a uniform state, both the weighted-residual equations and the Benney equations reduce via a weakly nonlinear analysis to a forced version of the Kuramoto–Sivashinsky equations (KS), for which optimal feedback controls have been successfully applied (Gomes et al. 2015). However, away from this weakly nonlinear limit, the strong nonlinearities present in the Benney and weighted-residual models present a significant challenge for feedback control. In Thompson et al. (2015), we show that applying feedback in the form of blowing and suction based on observations of the interface height can be successfully used to control the film behaviour represented by the nonlinear long-wave models derived in this paper. We focus especially on the robustness of control strategies of varying complexity applied in a hierarchical manner starting with application to simple models (e.g. KS) and then used in more complicated models (such as the Benney, weighted-residual and Navier–Stokes equations).