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 Reference Nusselt1923) 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. (Reference Tseluiko, Blyth and Papageorgiou2013) 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 Reference Benney1966) and on a weighted-residual method (following the approach of Ruyer-Quil & Manneville (Reference Ruyer-Quil and Manneville2000)); 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 (Reference Gjevik1970) 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. Reference Scheid, Ruyer-Quil, Thiele, Kabov, Legros and Colinet2005).

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. (Reference Thompson, Gomes, Pavliotis and Papageorgiou2015).

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$ .

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 (Reference Malamataris, Vlachogiannis and Bontozoglou2002) and Dietze, Leefken & Kneer (Reference Dietze, Leefken and Kneer2008). 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.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 Reference Benjamin1957). 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.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. Reference Scheid, Ruyer-Quil, Thiele, Kabov, Legros and Colinet2005), 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.