Falling liquid films with blowing and suction

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 with 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 nonuniform steady states is investigated for perturbations of arbitrary wavelengths, 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 suction amplitude increases.


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 1916) 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 = ρ 2 gh 3 s sin θ/2µ 2 (see (2.4) also) where ρ is the fluid density, g is the gravitational acceleration, θ is the angle of inclination to the horizontal (see Figure 1), and µ is the viscosity of the fluid; the flow becomes linearly 2 A. B. Thompson, D. Tseluiko and D. T. Papageorgiou unstable to long waves when R > R c = (5/4) cot θ 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.
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), in the sense that far away from critical in the unstable regime solutions can become unbounded in finite time (Pumir et al. 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. 2004). 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 spatiotemporal 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. 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 underpredicts the critical Reynolds number. Using a weighted-residual method, Ruyer-Quil & Manneville (2000) recently developed a new two-degree 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.
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 nonuniform 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 et al. 2009;Ogden et al. 2011;Samanta et al. 2011Samanta et al. , 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 et al. 1991;Renardy & Sun 1994), whereas an electric field applied normal to the interface can destabilize the flow and in fact drive it to chaotic spatiotemporal dynamics even below critical R < R c (Tseluiko & Papageorgiou 2006. 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 critical Reynolds numbers, Falling liquid films with blowing and suction 3 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 et al. 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 nonuniform 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 nonuniform 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 Blyth & Bassom 2012;Ogden et al. 2011). 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 et al. (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 break up streamlines. 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 et al. (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. The magnitude of the blowing/suction is assumed to vary spatially along the planar substrate and hence modifies the no penetration boundary condition there, but we assume that there is no slip along the substrate. We envision that such a model would be appropriate to experimental setups where tiny slits on the substrate would provide the conduit for fluid to enter and leave the wall. Such  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 wall inclined at angle θ 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).
mechanisms 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 increase the total mass and hence the mean thickness at any given time. An increase in thickness would consequently increase the local (in time) 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.

Nondimensionalisation and scaling
We wish to determine the evolution of a falling liquid film, with mean thickness h s , on a slope inclined at angle θ to the horizontal. We model the liquid as a Newtonian fluid of constant dynamic viscosity µ and density ρ, and the air as a hydrodynamically passive region of constant pressure p a . The coefficient of surface tension across the air/liquid Falling liquid films with blowing and suction 5 interface is γ. 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 which are coupled to the incompressibility condition 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 = ρgh 2 s sin θ/2µ (Nusselt 1916) as the velocity scale, h s /U s as the time scale, and µ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 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.5) At the interface, y = h(x, t), the tangential and normal components of the dynamic stress balance condition yield respectively. The kinematic boundary condition on the interface can be written in integral form as where q(x, t) is the stream-wise flow rate, defined as (2.9)

Long-wave evolution equations
We now seek solutions with wavelength L 1, and we introduce the long-wave parameter δ = 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 weightedresidual method (following the approach of Ruyer-Quil & Manneville 2000); derivations of both models are presented in appendix A. We assume that cot θ = O(1). To retain both inertia and surface-tension effects, we additionally assume that R = O(1) and C = O(δ 2 ). We choose the canonical scaling F = O(δ) 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.8), and so we have (2.10) 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.11) The first-order weighted-residual approach (see § A.2) instead yields an evolution equation for q: The Benney and weighted-residual equations are identical when R = 0. For R = 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, despite the additional complexity, the weighted residual equations in fact support bounded solutions across a greater range of parameter space (Scheid et al. 2004).

Choice of blowing and suction function
In time-dependent evolution, 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 case, 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 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, as the mean layer height does not change with time.
Most of the results we present are for fixed mean layer height, but we will also present some results for fixed flow rate, and we note that there is a mapping between results for these two conditions.
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.13) This function F (x) has the symmetry that the transformation A → −A is recovered by Falling liquid films with blowing and suction 7 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.
The long wave equations (2.10), (2.11) and (2.12) also apply if F is unsteady, so long as F varies on dimensionless timescales no shorter than O(1/δ), or the dimensional timescale h s /(δU s ). Time derivatives of the vertical velocity, and hence F , would feature in the equations at the next order in δ. 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.

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 & Oldman 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 pseudospectral 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 F(u,u) = 0 (3.1) so that steady solutions u 0 satisfy F(u 0 , 0) = 0. In order to determine the linear stability of a steady solution, we suppose that We now expand (3.1) for small , to obtain which is a generalised eigenvalue problem for λ and v, where J is the Jacobian matrix and M 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 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 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 smallamplitude states.

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 equations, 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.

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π/L, and seek a steady solution for h and q when |A| 1 of the form The mass conservation equation (2.10) immediately suppliesQ = exp(imx)/(im). However, the complex constantĤ differs between the two long-wave models. The Benney result, obtained by linearising (2.11), iŝ while the weighted residual equation (2.12) giveŝ The resulting solutions forĤ are illustrated in figure 2 for L = 10 and L = 40. The two expressions (4.2) and (4.3) are equal only if R = 0; otherwise both the magnitude and phase ofĤ 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 = δM and C = δ 2 C, and expand for small δ; we find that both models yield which is the expected order of agreement given that terms beyond the second order in δ were neglected in the derivation of each model. We note from (4.4) that the magnitude ofĤ 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 (4.4) also reveals a phase shift between h and F , which tends to π/2 as δ → 0. 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.10) to yield where we have used the fact that the spatially-averaged flux is even in A. We then expand the interface height in A as We solve the flux equation, either (2.11) or (2.12), at O(A) and O(A 2 ) to determine the constants p 1 , p 2 , p 3 , p 4 and Q 2 . The resulting coefficients are somewhat lengthy and so we do not list them here. In the long wave limit m = δM , C = δ 2 C, with δ 1, both the Benney and weighted residual equations yield so that small amplitude, long wave blowing and suction always increases the mean downslope flux.  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πx/L), L = 10, C = 0.05, θ = π/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 ( ).

Steady solutions as A increases
As A increases, the film profile deviates from the uniform state, and the nonlinearities in equations (2.11) and (2.12) 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, θ and C subject to constrained mean layer height 1.
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, 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 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. An alternative solution measure is the mean fluxq, 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 fixedq = 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 fixedh = 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.

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 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 fixedq = 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 A. B. Thompson, D. Tseluiko and D. T 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. 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 subcritical nature of 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).

Streamfunction and flow reversal
In addition to altering the interface height, the imposed suction also fundamentally affects the arrangement of streamlines within the fluid film. As the flow is incompressible, we can write the velocity in terms of a streamfunction ψ(x, y) such that (u, v) = (ψ y , −ψ x ). Under the weighted-residual formulation, we have for steady solutions, we can write the solution for ψ, up to the addition of an arbitrary constant, as The same result applies to O(δ) in the Benney equations.

13
According to (4.9), there are no stagnation points in 0 < y ≤ h if q > 0, and points along the wall are stagnation points if and only if F (x) = 0. In the absence of 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 Substituting (4.10) into (4.9) yields If F (x 0 )q(x 0 ) > 0, the stationary point at x 0 is a saddle point of ψ, and the streamlines are locally straight lines, such that The streamlines become increasingly vertical as q(x 0 ) → 0. In contrast, if F (x 0 )q(x 0 ) < 0, the stationary point at x 0 is a local extremum of ψ, and the streamlines are locally elliptical, with For any smooth, periodic, non-zero F (x) with mean zero, there must be some stagnation points x 0 with F (x 0 ) > 0 and others with F (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 * 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 (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 (x) > 0 are saddle points of ψ, while those with F (x) < 0 are extrema of ψ. 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 (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 injection layer thickens.
At A = A * , the stagnation point at x = x 0 , y = 0 switches from a saddle point to an extremum of ψ. The corresponding change in the flow field is illustrated in the transition from figure 6(b) to 6(c). For A > A * , q(x) is negative for some x, and hence by (4.8) 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 * . 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). (e) There are no steady states for A = 0.8; this is a final snapshot before drying. Figure 6: Solutions for the interface shape, velocity field and stagnation points for θ = π/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 (x) < 0 and with q(x)F (x) > 0 are indicated by and respectively. In (a-b), q > 0 everywhere, so a streamline emanating from the stagnation point with q(x)F (x) > 0 divides fluid into particles which never meet the wall (yellow), and particles which enters and leaves through the wall (blue). For A > A * = 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.

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(ikx + λt). The Benney model yields the eigenvalue λ directly, as For the weighted residual model, the linear stability analysis yields a quadratic equation for λ: In both models, the stability threshold for fixed k is 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 λ = λ 0 = −2ik at R = R H , and positive real part for R > R H . We note that (5.2) also has a second root for λ, but this root always has negative real part when R > 0.

Perturbations of arbitrary wavelength about a periodic base state
When the base state is spatially uniform, the eigenmodes can be written as (exp(ikx + λt)), and by calculating λ for all real k, we can determine linear stability to 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 forced equations, and consider the evolution of arbitrary small perturbations, so that The weighted residual equivalent of (5.6) additionally involves a term proportional toq 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ĥ where the Floquet wavenumber k is real, andĥ k andq 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 (λ) < 0 for all eigenfunctions when k = 0. The base solution is linearly stable to perturbations of all wavelengths if (λ) < 0 for all eigenfunctions for each real k ∈ [0, π/L].

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 λ for R close to R H , and for small A: We are concerned with eigenvalue equal to λ 0 = −2ik at R = R H , which is a root of both the Benney and weighted residual characteristic equations. We can evaluate the term ∂λ/∂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 → −A can be recovered by translation in x by a distance L/2, and the original eigenmode exp(ikx + λt) has no preferred position; thus the leading order contribution for small A is O(A 2 ). The effect of the imposed suction is encapsulated in the coefficient Z, which depends on the suction wavenumber m, the perturbation wavenumber k, and C and θ. In order to calculate W and 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 and the base solution H, Q 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 functionŝ h k andq k are periodic with wavenumber m and are constant when A = 0, and can be expanded for small A aŝ 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 h k (0) = 1. (5.13) To determine the unknown constants, we set R = R H , 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 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, n, C and θ; 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 {λ} = 0, so that (5.14) Figure 7 shows the effect of small A on the stability boundary, as quantified by the weighted-residual results for W . For θ = π/4, we find that forcing via blowing and suction at wavelength L ≡ 2π/m = 10 destabilises short wavelength perturbations, but has little effect on long wavelength perturbations. However if the forcing wavelength L ≥ 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. For θ = π/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 → 0, and so the imposed forcing in fact slightly destabilises long wave perturbations. For θ = 3π/4, the behaviour for small perturbation wavenumber is reversed from the θ = π/4 case; long wave forcing is always destabilising for k = 0, and becomes increasingly destabilising as L is increased.
In figure 8, we survey the sign of W for k = 0 and k = m, with m and C varying, for three values of θ, for both the Benney and weighted-residual models. The two models concur that for each inclination angle θ, there are some regions of (C, m) parameter space where the imposition of blowing and suction stabilises the uniform state, and other regions where such forcing destabilises the flow. For θ = π/4, both models predict that long wave forcing, with m → 0, is stabilising to perturbations with Floquet wavenumber k = 0 and with wavenumber k = m. For θ = π/2, forcing is destabilising in the limit k → 0 for all C and m, but has a stabilising effect on perturbations with k = m if m is not too large. For θ = 3π/4, the long wave forcing is destabilising to perturbations both with k = 0 and with k = 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 long-wavelength blowing and suction and long perturbation wavelength, setting m = δM , k = δK and C = δ 2 C, and expanding the full expression for Z with δ 1. Both the 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 θ = π/4, π/2, 3π/4 are shown in rows (a-c) respectively.

Benney and weighted-residual results yield
(5.15) Substituting (5.15) into (5.14) yields the leading order long-wave expansion (5.16) The prediction (5.16) is plotted in figure 7 for comparison to the non-long-wave results obtained by direct evaluation of the Maple expressions for Z and ∂λ/∂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 R * (m) = min k∈R 5 4 cot θ + 5k 2 8C + A 2 W (m, k, C, θ) . (5.17) Using the long wave expansion for W given by (5.16), we obtain The term in square brackets is positive, and so the minimum of the expression occurs

Falling liquid films with blowing and suction
19 at k = 0. As expected, this implies that long-wave perturbations are the first to become unstable as R is increased, both for A = 0 and also for small A in the long wave limit. Proceeding with the long wave approximation, we evaluate the minimum of (5.18) as while if perturbations restricted are restricted to those with wavenumber k = m, we find The stability boundaries for perturbations with all 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 θ = π/4, R * is positive in the absence of suction, and long-wave blowing and suction can either decrease or increase R * , depending on the values of θ, m and C, and here it is reasonable to consider an optimisation strategy. To maximise the stabilising effect at fixed A for θ < π/2, we should choose m as small as possible, and in fact (5.19) predicts that R * → +∞ as m → 0 with A fixed. Regarding constraints on the magnitude of film height variations, we note that these are of order A/m for A 1 (see § 4.1) , and so we can approximately constrain height variations by keeping α = A/m fixed. We find that the largest R * is again obtained as m → 0, but in contrast to the fixed A case, R * is bounded as m → 0.
The governing equations are valid for 0 < θ < π, with θ = π/2 corresponding to flow down a vertical wall, and π/2 < θ < π corresponding to flow along the underside of an inclined plane. In the absence of blowing and suction, R * = 0 for θ = π/2, and R * < 0 for π/2 < θ < π. Negative Reynolds numbers are not physically relevant, but consideration of (5.19) may still give some indication of the effect of blowing and suction when θ > π/2. According to (5.19), when θ ≥ π/2, imposing blowing and suction at long wavelengths always reduces R * , and so is destabilising. Our other results for θ = 3π/4, shown in figures 7 and 8, also suggest that suction destabilises long wave perturbations. Figure 9 shows stability regions as R and A vary, at fixed values of C and θ, for imposed suction with wavelength L = 10 and L = 40. For L = 10, introducing finite amplitude blowing and suction decreases the critical R for linear stability 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.

Finite A stability results
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ĥ = e ikx with k = 2π/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 0 < k 2 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 2 c are unstable, and there is a finite k 2 with maximum growth rate (λ). However, it is possible that for larger A, imposing blowing and suction can 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. For the cases shown in figure 9, we find that at small A, the boundary for stability to perturbations of arbitrary wavelengths 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. Results from time dependent simulations in and around the island are shown in figure 14. 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.

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 question, requiring finite amplitude results for the steady solutions and for their linear stability operator. We can address this task numerically, by choosing a suction wavenumber m, increasing R slightly from R 0 = 1.25 cot θ, and testing the stability of steady solutions beginning from A = 0 until A is too large for steady states. Figure 10 shows numerical results from this search for the case θ = π/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.
Firstly, at large wavelengths with L ≈ 200, imposing blowing and suction allows the critical Reynolds number for stability to all perturbations to increase from the unforced value of 1.25 to around 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 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, but there is no such stable island in the results for L = 40.

Travelling waves in the absence of suction
In the absence of suction, the system can support travelling waves, propagating at a constant speed U without changing form. We can write the solution as where U is the unknown constant wave speed. H and Q must satisfy where a prime indicates a derivative with respect to ζ. If H(ζ) 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 θ = π/4, C = 0.05, R = 3. 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 :

4)
H 2 (ζ) = r 1 cos 2kζ + r 2 sin 2kζ, q(x, t) = U H(ζ) + S, (6.5) We expand the equations for 1, and must solve the equations at up to O( 3 ) in order to determineR. 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 ifR > 0 and subcritical ifR < 0. The Benney equations yield which may be positive or negative. However, for the weighted-residual equations, we find R = 4410C 2 + 6670k 2 C 2 cot 2 θ + 12235k 4 C cot θ + 4450k 6 2352Ck 4 . (6.8) When θ < π/2, (6.8) yieldsR > 0 and so the Hopf bifurcation in the weighted residual model is supercritical. In the long-wave limit, with k = δK and C = δ 2 C, both (6.7) and (6.8) yieldR which is positive. The quantityR 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. 2004), and so the value ofR does not necessarily prescribe the stability of finite-amplitude travelling waves.

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 final part 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 amplitude of the imposed blowing and suction is increased.
We now 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 where the behaviour ofĥ andq is in some way related to the periodicity of the original travelling wave with respect to ζ 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), h t − f (x) +q x = 0 (6.11) and q + w 6 (ζ)q t = w 0 (ζ)ĥ + w 1 (ζ)ĥ x + w 2 (ζ)ĥ xxx + w 3 (ζ)q + w 4 (ζ)q x + w 5 (ζ)f (x), (6.12) 14) For non-zero f , we can integrate forward in time from a given initial conditionĥ(x, 0) =ĥ 0 (x),q(x, 0) =q 0 (x), to determineĥ(x, t) andq(x, t).
The system (6.11) and (6.12) features some terms that are known functions of ζ, and others that are known functions of x, but the system is autonomous with respect to t. We therefore choose to rewrite the equations (6.11) and (6.12) in terms of x and ζ: (6.16) We obtain a system in two variables, x and ζ, with equations − Uĥ ζ − f (x) +q x +q ζ = 0 (6.17) and q − U w 6 (ζ)q ζ = w 0 (ζ)ĥ + w 1 (ζ)(ĥ x +ĥ ζ ) + w 2 (ζ)(ĥ xxx + 3ĥ xxζ + 3ĥ xζζ +ĥ ζζζ ) +w 3 (ζ)q + w 4 (ζ)(q x +q ζ ) + w 5 (ζ)f (x). (6.18) We now regard ζ and x as independent variables. Under this transformation, x remains as a purely spatial variable, but ζ has a dual identity, incorporating both spatial and timelike components. The statement of initial conditions becomes more complicated in the new variables:ĥ so that the initial conditions are spread across the whole range of ζ.
If the Fourier expansion of f (x) is (6.20) we can write the general solution of (6.17) and (6.18) aŝ where the functions J m (ζ), K m (ζ), M m (ζ) and N m (ζ) are periodic in ζ. The remaining terms, h * and q * , 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 underlying travelling wave. We can calculate the periodic functions J m , K m , ..., 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 * , q * to the homogeneous problem eventually decay to zero and soĥ andq are limit cycles, periodic in time.
We now calculateĥ 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 − U J m + mN m + M m = 1, (6.23) +w 3 M m + w 4 (M m + mN m ) + w 5 (6.25) and (6.26) 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 ζ with period L. Thus we obtain the solution Regardless of the relation between the travelling wave period L and the blowing/suction wavenumber m, we see that if J m (ζ) and K m (ζ) are periodic with period L, then the functionĥ(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ĥ(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(ζ), J m (ζ) and K m (ζ), and also the reconstructed fieldĥ, are shown in figure 11. For the travelling wave solution shown in figure 11, we find that the functions J m (ζ) and K m (ζ) display little relative variation with ζ. This means that the fieldĥ(x, t) is essentially steady. 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 wavelength of the blowing and suction.
The composite solution shown in figure 11 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 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 can cause the underlying travelling wave to double in spatial period, and likely change shape and speed. Eventually we would expect to reach the state where the dominant wave has period 60, and the linear perturbation fieldĥ 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 initial condition 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 with wavelength equal to the initial forcing, 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 state regardless of initial conditions.

Initial value problems
Even in the absence of blowing and suction, a thin liquid film falling along an inclined plane can display rich behaviour, including pattern formation, transition to chaos, travelling waves, and pulse-like structures. Many of the instabilities that dominate the dynamics occur at long wavelength, and so the system is frequently studied using longwave models. However, a well-known feature of the Benney model is that the interface can exhibit finite-time flow up at Reynolds numbers larger than critical (Pumir et al. 1983;Ruyer-Quil & Manneville 2000), which is not replicated in full Navier-Stokes simulations. The weighted-residual equations were developed partly to avoid this blow-up behaviour, and are indeed better behaved that the Benney models. The two models agree when the system parameters are close to neutral stability.
The introduction of forcing in the form of suction boundary conditions introduces another potential mechanism for finite-time blow-up, whereby the film thins to zero in finite time, as illustrated in figure 4. This thinning occurs even at zero Reynolds number, and so arises in both Benney and weighted-residual models. The blow-up need not occur at the same wavelength as the forcing function; for example it can be triggered by subharmonic-perturbations to a periodic steady state. 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.
Blow-up is generally avoided if R is not too large, and blowing and suction is applied with sufficiently small amplitude, but the film dynamics can still interact with the imposed suction. 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 time-periodic behaviour manifests as nonlinear limit cycles which depend 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 non-zero perturbation wave speed even when the interface shape is steady.
The applied suction can also lead to states which are neither steady nor limit cycles. 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 case shown in figure 14 , and varying initial conditions: we set h(x, 0) = 1 + 0.1 cos(2πnx/60) and q(x, 0) = 2/3 + 0.2 cos(2π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. 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 in a regular manner, as shown in figure 14(c).

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 function F (x, t). We derived and studied two long-wave models, based on the first-order Benney and weightedresidual formulations, for this system. 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 F must themselves be non-uniform. We calculated the interface shape for small A, and showed that this differs between the two models when R = 0, but that the results agree in the long-wave limit. There is a phase difference between the shape of the blowing and suction and interface height, which originates from the preferred down-slope direction of the base flow.
As A increases, the interface becomes increasingly non-uniform, and the nonlinear terms in the governing equations become increasingly important. Two notable transitions occur as A increases: firstly all steady states feature regions with negative volume-flux q, and secondly steady solutions disappear altogether for sufficiently large A. Time-dependent evolution beyond the existence of steady states shows that the film can dry in finite time, as the fixed rate of fluid removal via F is not matched by the supply of fluid.
We analysed the behaviour of the streamfunction in order to interpret steady solutions with negative q. We found that imposing non-zero, non-uniform F leads to a series of isolated stagnation points along the wall, wherever F (x) = 0. These stagnation points are either saddle points or local extrema in the streamfunction, depending on the sign of F (x) and q(x, t) at these points. Connections between the stagnation points allow us to demarcate the boundary between fluid which enters and leaves through the wall, and 'propagating' fluid which never reaches the wall. At A = 0, all fluid propagates down the wall, but the height of the recirculating layer near the wall increases with A. For large enough A, all steady solutions have negative 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.
Steady solutions need not have the same spatial periods as the imposed suction, and subharmonic steady solutions can emerge via pitchfork bifurcations from the harmonic solution branches. In our calculations, we have only found subcritical pitchfork bifurcations, leading to unstable subharmonic steady states. 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.
In the absence of suction, the primary instability of the uniform state is to eigenfunctions with complex eigenvalues, leading to Hopf bifurcations which do not appear in the bifurcation structure for steady solutions. If F is spatially periodic, then the steady states are also periodic, as are the coefficients of the linear equation governing the evolution of small perturbations. This means that perturbations must propagate through a heterogeneous environment, and also that the eigenmodes are no longer simply exponentials of the form exp(ikx + λt). Instead, according to the Floquet-Bloch form, the eigenmodes can be written as g(x) exp(ikx + λt), where g(x) is a complex-valued function, spatially periodic with the same period as the base state. We calculated linear stability 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 introducing spatially-periodic suction with amplitude A can increase or decrease the critical Reynolds number for linear stability to perturbations of both classes. Due to the symmetries of the system, the correction to the eigenvalue, and hence to the critical Reynolds number, is even in A.
For small A, we calculated the critical Reynolds number to O(A 2 ) for perturbations of arbitrary wavenumber. Imposing blowing and suction at very long wavelength increases the critical Reynolds number for both classes of perturbations when θ < π/2. Due to the quadratic dependence of the critical R on perturbation wavenumber, the first eigenmodes to become unstable as R increases are those with infinite wavelength, both when A = 0 and for small A, In contrast, for the case of flow underneath an inclined plane, where θ > π/2, small-amplitude long-wave blowing and suction always reduces the critical Reynolds number, and so has a destabilising effect on the flow. However, as the small A expansion is valid only close to the critical Reynolds number, which is negative for θ > π/2, the predicted dependence of stability properties on A may not hold for physically realisable flows.
Determination of the largest R that can be stabilised in an infinite domain by imposing steady suction is difficult as it is not clear that long-wavelength waves are always the most unstable, and both the steady solutions and their corresponding linear stability operators Falling liquid films with blowing and suction 31 must be computed at finite A. We find that the weighted-residual model predicts large increases in the critical Reynolds number for stability to perturbations of all wavelengths, from 1.25 to 4 or more, if blowing and suction is introduced with wavelengths of around 200. However, the imposed suction appears to become less effective at stabilising the flow if applied with very long wavelength.
In the absence of suction, the system can support travelling waves which propagate at a constant speed U without changing form. Periodic travelling waves are periodic with respect to both space and time, but their dependence on x and t can be expressed in terms of a single variable ζ = x − U t. Imposing blowing and suction as F (x) introduces heterogeneity to the system, and so we expect travelling waves to become limit cycles that are periodic in both x and t, with explicit dependence on both variables. We calculated the O(A) corrections to the travelling wave, and showed that the equations are closely related to those governing linear stability of the travelling wave. We demonstrated that limit cycles at O(A) can be decomposed by rewriting the equations in terms of the independent variables x and ζ. The resulting solution for h(x, t) and q(x, t) is periodic with respect to 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. We showed that the predictions of the O(A) perturbation analysis are in good agreement with the results of fully nonlinear simulations. We also conducted initial-value calculations to investigate competition between the period of the travelling wave and of the imposed suction. We found that the same state is eventually reached, but this competition can persist over a large number of cycles.
We observed in § 2.3 that the derivation of our equations is unchanged if the blowing and suction profile additionally depends on time with sufficiently long timescale. Allowing time-dependence of the blowing and suction leads to the possibility of using F to explicitly control the flow, and thus to stabilise otherwise unstable states by acting in response to the development of perturbations. For example, with perfect implicit knowledge of the system, we could choose F (x, t) in order to drive the system towards a particular steady state. However, the additional independent degree of freedom that emerges in the weighted-residual model may be important when exploring control of the system, as we cannot choose F (x, t) to simultaneously specify h(x, t) and q(x, t). Under the assumption of small deviations from a uniform state, with small F , both the weightedresidual equations and the Benney equations reduce via a weakly nonlinear analysis to a forced version of the Kuramoto-Sivanshinsky equations, for which optimal controls have been successfully applied . Work is in progress (Thompson et al. 2015) to explore the effectiveness of feedback control strategies based on the nonlinear long-wave models derived in this paper.