3.1. Time-dependent simulations of the lubrication equation (2.7)

The lubrication equation (2.7) is solved numerically using a pseudospectral scheme (e.g. Trefethen Reference Trefethen2000) on a periodic domain $[-L,L]$ and the unconditionally stable Crank–Nicolson method is used for the time discretisation. The discretised nonlinear system is solved at each time level by using Newton's method. The (unknown) flow rate $Q(t)$ is calculated at each time step by imposing the zero-pressure-drop condition (A12), which is done numerically by setting the integral of the pressure gradient (A11) to zero and solving for $Q(t)$ (integrals are approximated by the trapezoidal rule). In all the simulations presented in this work, the (half) domain length is fixed to $L=28$, which is sufficiently large to capture the phenomena of interest. Typically, $512$ grid points are used to discretise the spatial domain, so the grid resolution is $2L/512=0.0194$, and a time step of size $0.05$ is found to be sufficient to produce accurate results.

The initial condition corresponds to a sinusoidal disturbance superimposed on a flat interface at uniform height $h_0$, and is given by

(3.1)\begin{equation} h(x,0) = h_0\left[ 1 + h_A \cos\left(\frac{n{\rm \pi} x}{L}\right) \right], \end{equation}

where $h_A$ is a constant and $n$ is an integer. In what follows we describe the dynamics for the particular case $h_A=0.2$ and $n=1$. Following extensive simulations with other choices of these parameters (some further comment on which is provided below), including combinations of Fourier modes and amplitudes $h_A$ selected from a random distribution, this was found to be representative of the dynamics that is observed. Unless otherwise stated, the parameters used in the numerical computations are $h_0=0.2$, $m=0.5$, $C=10^{-3}$ and $B=1$.

Interfacial profiles at different stages of the simulation starting from the initial condition (3.1) and using the aforementioned parameter set are shown in figure 1. At early times a number of protrusions, hereinafter referred to as fingers, appear on the interface. In this case six fingers emerge, as can be seen in figure 1(a–c) corresponding to $t=15,20,35$. Intuition suggests that the number of fingers will be dependent on the size of the Bond number, since a greater density disparity between the two fluids should correspond to a greater susceptibility to interfacial instability. While in the present results it takes approximately 35 time units for all six fingers to clearly form at the interface, we note that for larger values of $B$ more fingers develop and also these emerge earlier in the simulation.

Figure 1.

Interfacial profiles obtained at different times (a) $t=15$, (b) $t=20$, (c) $t=35$, (d) $t=50$, (e) $t=100$ and ( f) $t=200$ for the parameter values $h_0=0.2$, $m=0.5$, $C=10^{-3}$ and $B=1$. Starting from an initial condition of the form (3.1) with $h_A=0.2$ and $n=1$ (also shown with a broken line in panel a), six distinct fingers are seen to develop at early times in the simulation. These then get closer to the lower wall and slowly drift to the right due to background shear flow. The tallest finger travels faster and eventually merges with the finger directly in front of it. The merging dynamics continues to take place until a single wide finger remains that almost touches both channel walls. The solution seen in ( f) has been shifted to the right by a distance $L=28$ for better visualisation of the final finger. The time evolution of the interfacial profile can be seen in the supplementary movie available at https://doi.org/10.1017/jfm.2022.1070.

The number of fingers formed can be estimated using the growth rates determined from linear theory (e.g. Tseluiko & Papageorgiou Reference Tseluiko and Papageorgiou2007). These are calculated from (2.7) by writing $h=h_0+A\exp (\mathrm {i} Kx+\sigma t)+ \text {c.c.}$ (complex conjugate), and linearising on the basis that $A/h_0\ll 1$ with a view to calculating $\sigma$ for a prescribed wavenumber $K$ (more details are provided in Appendix B). In figure 2, the linear growth rate, $\mathrm {Re}(\sigma )$, is plotted against the spatial frequency $k = L K /{\rm \pi}$ for various values of the Bond number $B$. In the context of the numerical simulations of the evolution equation (2.7), the variable $k$ is taken to be an integer corresponding to the number of wavelengths that fit exactly into the computational domain of length $2L$. The fastest growing mode in the numerical simulation corresponds to the spatial frequency closest to that which maximises the growth rate, namely $k_{max} = (L/{\rm \pi} )\sqrt {B/2}$ (see Appendix B). For example, when $L=28$ we have $k_{max}\approx 6$ for $B=1$ and $k_{max}\approx 9$ for $B=2$, the former prediction being in line with the six fingers that are observed in figure 1(c). The cut-off wavenumber for linear instability is $K_c=\sqrt {B}$ so that the range of unstable wavenumbers depends solely on the Bond number (see also Oron et al. Reference Oron, Davis and Bankoff1997). Accordingly, in the simulations we must choose integer $k< LK_c/{\rm \pi}$ in order to observe instability.

Figure 2.

Linear growth rate for different values of the Bond number, plotted against $k=({L}/{{\rm \pi} })K$, where $K$ is the wavenumber and $L=28$ is the (half) domain length. Parameter values for $h_0$, $m$, $C$ are the same as in figure 1. A positive growth rate is supported for wavenumbers in the range $K\in (0,\sqrt {B})$, so interfacial instability is enhanced as the value of $B$ increases. The maximum growth rate occurs at $k_{{max}}=(L/{\rm \pi} )\sqrt {B/2}$.

Finger formation similar to that seen in figure 1 has also been observed in a related study of channel flow by Mavromoustaki, Matar & Craster (Reference Mavromoustaki, Matar and Craster2010). In the absence of the shear flow (so when the upper wall is stationary), the finger pattern in the interfacial profile retains the symmetry of the initial condition (3.1) approximately $x=0$ (see related studies on static liquid film configurations, for example those by Yiantsios & Higgins (Reference Yiantsios and Higgins1989), Bertozzi & Pugh (Reference Bertozzi and Pugh1998), Tseluiko & Papageorgiou (Reference Tseluiko and Papageorgiou2007), Wang, Mählmann & Papageorgiou (Reference Wang, Mählmann and Papageorgiou2009), Barannyk et al. (Reference Barannyk, Papageorgiou, Petropoulos and Vanden-Broeck2015) and Wang & Papageorgiou (Reference Wang and Papageorgiou2018)). Here the symmetry is broken by the superimposed shear flow and consequently the interfacial profile tends to drift in the positive $x$ direction. The individual fingers (as seen in the various panels of figure 1) travel more or less holding their shape, size and speed. Generally speaking, taller fingers have greater volume and move faster, as is confirmed in § 3.3. As a result the tallest finger eventually catches up with the finger directly in front of it, and the pair subsequently merge to form a larger finger whose volume is approximately the sum of the two individual volumes. An example can be seen in figure 1(c,d), where the interfacial profiles are displayed at $t=35,50$. After this first coalescence, the tallest finger has gained sufficient volume to approach the upper channel wall (see figure 1d). The merging process continues so that the number of fingers gradually diminishes until only a single, rather wide, finger remains (see figure 1f); this finger essentially combines the individual volumes of all fingers from the earlier stage of the simulation, with any deficit accounted for in the thin lower fluid film that comprises the remainder of the flow domain. Note that despite the apparently rather steep slope on either side of the finger, the scales are such that the gradient here is modest, so that the legitimacy of the assumed lubrication approximation ought to be maintained. With reference to alternative choices of initial condition mentioned above, broadly similar dynamics is observed, including the occurrence of merging events, the only difference being the number of final-state fingers (akin to that at $t=200$ in figure 1) that are ultimately obtained. Similar coarsening dynamics has been reported by Glasner & Witelski (Reference Glasner and Witelski2003) and Kitavtsev (Reference Kitavtsev2014) who investigated dewetting of thin films, as well as by Chang, Demekhin & Kalaidin (Reference Chang, Demekhin and Kalaidin2000) and Razis et al. (Reference Razis, Edwards, Gray and van der Weele2014) in their studies of roll wave structures propagating down an inclined plane.

The merging events just described can also be seen in the evolution of the displacement norm $\|h-h_0\|_2$ shown in figure 3(a), where each of the sudden upward steps in the curve corresponds to the merging of a pair of fingers. The final merging event occurs at around $t=150$ when the norm reaches the final plateau, and after this time the solitary finger propagates downstream essentially retaining the same form for the remainder of the simulation. Also shown in figure 3 is the trace of the interfacial minimum and maximum during the simulation. The minimum first decreases until around $t=35$ when the interface almost makes contact with the bottom wall; at the same time the maximum increases until the interface almost touches the upper wall (this occurs at around $t=50$; see also figure 1c,d). At later times the interfacial minimum and maximum settle down to a near-constant value; in fact they exhibit small-amplitude and apparently random fluctuations approximately a mean level, as can be seen from figure 3(d,e). This irregular large-time behaviour is reminiscent of the chaotic dynamics seen in the evolution of thin films (e.g. Kawahara Reference Kawahara1983; Kawahara & Toh Reference Kawahara and Toh1988; Kalogirou & Papageorgiou Reference Kalogirou and Papageorgiou2016). We explain this by noting that at large times the near-saturated state features a thin film of the lower fluid over much of the domain, and a thin film of the upper fluid confined between the top of the finger and the upper wall. The evolution of these thin films can be described by Kuramoto–Sivashinsky-type equations (see Appendix C and the related studies by Hooper & Grimshaw (Reference Hooper and Grimshaw1985), Charru & Fabre (Reference Charru and Fabre1994), Tilley et al. (Reference Tilley, Davis and Bankoff1994) and Kalogirou et al. (Reference Kalogirou, Cimpeanu, Keaveny and Papageorgiou2016)), which are well known for demonstrating chaotic behaviour (Papageorgiou, Maldarelli & Rumschitzki Reference Papageorgiou, Maldarelli and Rumschitzki1990). Both films are linearly unstable with a growth rate of the order of $10^{-3}$; its small size is attributable to the $h_0^3$ and $(1-h_0)^3$ terms in the growth rate formula (B2). In related studies where the flow evolves from a static initial state, the fluid film tends to approach the wall in infinite time (Tseluiko & Papageorgiou Reference Tseluiko and Papageorgiou2007). Such a touch-down scenario appears to be absent here, and we believe the near saturation combined with thin-gap chaotic dynamics is connected with the background shear flow. In Appendix D we provide an analytical argument to support the numerical results seen in figures 1 and 3, in particular that the solution remains bounded such that $0< h<1$ for all time.

Figure 3.

(a) Time evolution of interfacial solution norm, (b,c) interfacial minimum and maximum, including (d,e) a close-up evolution at later times of the simulation. Parameter values for $h_0$, $m$, $C$, $B$ remain the same as in figure 1.

The final interfacial wave profile in figure 1( f), also shown in figure 4(a), features a capillary ridge and a depression on the top-right and the bottom-right part of the finger, respectively (these features are more obviously seen in the sketch in figure 11; see also in figure 13b). The interfacial velocity, shown with a red thin line in figure 4(a), becomes large in magnitude at the two points $x_1$ and $x_2$ where the interface is closest to each wall (similar behaviour was observed in Papageorgiou, Petropoulos & Vanden-Broeck (Reference Papageorgiou, Petropoulos and Vanden-Broeck2005), Wang et al. (Reference Wang, Mählmann and Papageorgiou2009) and Barannyk et al. (Reference Barannyk, Papageorgiou, Petropoulos and Vanden-Broeck2015)). To understand the source of these peaks we plotted the magnitude of the various terms that define the interfacial velocity in (A8) to determine which term dominates. We found that the $\mathcal {F}$ term, defined in (2.8d), experiences (positive or negative) peaks as can be seen in figure 4(b), but these are most pronounced where the interface is closest to each wall (i.e. at the points labelled as $x_1$ and $x_2$). We note that the $\mathcal {F}$ term is related to $F_\chi$ in (A7b) and that this term dominates over the other terms in the pressure gradient expression (A11). The large negative peaks in the profile of $\mathcal {F}$ lead to large negative spikes in pressure gradient, which drives the rapid flow near the walls (note that the flat parts of the $\mathcal {F}$ curve are much smaller in magnitude, on the scale shown, but not constant).

Figure 4.

(a) Profile $h$ of the interfacial wave (blue line) and interfacial horizontal velocity $u_I=u|_{y=h}$ (red line) and (b) term $\mathcal {F}=C^{-1}(Bh_x + h_{xxx})$ which dominates in the pressure gradient distribution, plotted at time $t=200$. Parameter values for $h_0$, $m$, $C$, $B$ are the same as in figure 1. The interfacial velocity is seen to become large in magnitude at the two points $x_1$ and $x_2$ where the interface takes the form of a capillary ridge or depression and is closest to the upper/lower channel walls. This is accompanied by a drop in the pressure gradient as evident by the dips in the profile of term $\mathcal {F}$.

3.2. Physical discussion

Intuition might suggest that the action of gravity will lead to the heavier fluid displacing the lower, lighter fluid. In a static configuration for which there is no preference for flow to the right or left, one would anticipate that a localised interfacial disturbance will tend to migrate upwards as fluid is drawn in from either side by hydrostatic force. Since the shear flow removes the left/right symmetry, the possibility arises of a finite-amplitude travelling wave that is sustained by a tripartite combination of shear, gravity and surface tension. Indeed, this has been demonstrated in the weakly nonlinear regime by Babchin et al. (Reference Babchin, Frenkel, Levich and Sivashinsky1983). Examples of fully nonlinear waves in related unstably stratified problems include travelling waves on fluid films flowing down the underside of an inclined plane – see, for example, Kofman et al. (Reference Kofman, Rohlfs, Gallaire, Scheid and Ruyer-Quil2018) and Zhou & Prosperetti (Reference Zhou and Prosperetti2022).

Considering moderate- to large-amplitude waves, we now examine more closely one of the travelling-wave solutions obtained as the large time limit of the initial value problem for the lubrication model equation (2.7) solved in § 3. In figure 5 we show the final-state interfacial profile when $B=0.1$ together with the corresponding streamline pattern in a frame of reference that is travelling at the wave speed. Figure 6 shows a more highly deformed interfacial profile computed for $B=0.4$. We note the striking qualitative similarity between the interfacial profile seen in figure 5 and that computed for a viscous film flowing on the underside of an inclined plate (see e.g. Kofman et al. Reference Kofman, Rohlfs, Gallaire, Scheid and Ruyer-Quil2018; Zhou & Prosperetti Reference Zhou and Prosperetti2022).

Figure 5.

(a) Flow configuration and (b) interfacial curvature for $B=0.1$ and parameters $h_0$, $m$, $C$ fixed to the same values as in figure 1. In (a), the colour corresponds to the vertical velocity $v_1$ in the lower fluid, the thick line indicates the location of the interface and thin lines shown within the fluids are the streamlines.

Figure 6.

(a) Flow configuration and (b) interfacial curvature for $B=0.4$ and parameters $h_0$, $m$, $C$ fixed to the same values as in figure 1. The colour in (a) corresponds to the vertical velocity $v_1$ in the lower fluid. In (a,c–f), the thick line indicates the location of the interface and the thin lines shown within the fluids are the streamlines. Two clockwise-moving eddies exist in each fluid (panel a), and the streamlines are seen to cross at two stagnation points at $(x,y)=(15.6,0.5)$ and $(22.85,0.53)$, indicated by red dots in (e, f). Near the transition regions where the interface ascends/descends from one wall to the other, a total of four smaller counterclockwise-moving eddies are generated due to the draught from the surrounding fluids, as shown in (c,d).

The speed of the travelling frame is 0.18 for figure 5 and 0.465 for figure 6, both of which values are smaller than the upper wall speed, which is unity according to our non-dimensionalisation. Consequently, in both figures the flow is from right to left in the lower part of the channel and from left to right in the upper part of the channel. Focusing on the results in figure 5, located in the main part of the upward bulge in the lower fluid is an eddy with fluid circulating in the clockwise direction. Critically, the interface is not symmetric about the vertical line through the eddy centre. The fluid is drawn upward on the left-hand side of the eddy by hydrostatic force. On the right-hand side of the eddy the large curvatures that occur around $x=-2.5$ and $x=5$ suggest that the capillary force is sufficiently strong to combat the hydrostatic effect and to drive the fluid downwards. The whole structure is sustained by the background shear flow which flushes the rest of the lower fluid away to the left.

For fairly small Bond number, the effect of surface tension is relatively strong and only reasonably mild interfacial curvature is needed to provide the necessary capillary pressure to support the eddy motion. For larger Bond number the relative effect of surface tension is weaker and much larger curvatures are needed to provide the necessary capillary pressure. Thus we see the notably sharper features in the interfacial profile for $B=0.4$ in figure 6 at approximately $x=17$ and $x=21$, and the much stronger curvatures at the same locations, roughly a factor of ten larger than the corresponding values for the $B=0.1$ case. Interestingly, the flow structure is more complex for the larger Bond number. There are many stagnation points within a single flow period: eight are located at the centre of eddies (note the four very slender eddies adjacent to the interface, which can be seen in figure 6c,d), two more are located at local saddle points in the streamline pattern evident in figure 6(e,f) and several others are located on the interface itself.

3.3. Travelling-wave solutions

The results of the numerical simulation discussed in the previous section indicated that the interface evolves into a travelling wave. Although the flow was found to be chaotic in the narrow gaps at the walls, this does not significantly affect the rest of the wave profile so that overall the interfacial travelling wave propagates with approximately fixed form and constant speed in the direction of the background shear flow. Motivated by this, in this section we seek travelling-wave solutions to the evolution equation (2.7) using the continuation software AUTO-07p (Doedel & Oldman Reference Doedel and Oldman2009). Introducing the travelling-wave coordinate $\xi =x-ct$, with wave speed $c>0$, equation (2.7) becomes

(3.2)\begin{equation} - ch_\xi + q_\xi = 0. \end{equation}

Assuming steady flow in the travelling-wave frame, we take the total flow rate $Q$ to be constant.

To prepare the ground for the computations in AUTO, we define the scaled independent variable $X=(\xi +L)/(2L)$, so that $X\in [0,1]$, and rewrite (3.2) as the first-order system of ordinary differential equations $\boldsymbol {U}_X = 2L\mathcal {A}(\boldsymbol {U})$, where

(3.3a,b)\begin{align} \boldsymbol{U} \equiv (U_1, U_2, U_3, U_4) \equiv (h, h_X, h_{XX}, h_{XXX}) \quad \text{and} \quad \mathcal{A}(\boldsymbol{U}) \equiv \left(U_2, U_3, U_4, \mathcal{N}(\boldsymbol{U})\right). \end{align}

The nonlinear function $\mathcal {N}(\boldsymbol {U})$ in (3.3a,b) is found by rearranging equation (3.2) for $h_{XXXX}$; see also the expression for $q$ in (2.7b) and the definition of $\mathcal {F}$ in (2.8d). Periodic boundary conditions are imposed on $h, h_X, h_{XX}$. Three integral constraints are imposed to fix the volume of each fluid in a domain period, to break the translational invariance of the system (Sandstede Reference Sandstede2002; Champneys & Sandstede Reference Champneys and Sandstede2007) and to enforce the zero-pressure-drop condition (A12). The parameter continuation is initiated close to the bifurcation point for small-amplitude periodic waves using the approximate solution $h(X) = h_0 + 10^{-3}\cos (2{\rm \pi} X)$. The remaining parameters are determined so that the initial guess is close to neutral stability, according to predictions from linear theory (see Appendix B).

The number of free parameters in the continuation is determined by the relation $N_{cont}=N_{bc}+N_{int}-N_{dim}+1$, where $N_{bc}$ is the number of boundary conditions, $N_{int}$ is the number of integral constraints and $N_{dim}$ is the size of the system (Doedel & Oldman Reference Doedel and Oldman2009). Here we have $N_{bc}=N_{int}=3$ and $N_{dim}=4$ (cf. (3.3a,b)), and hence $N_{cont}=3$. We perform continuation in terms of the travelling-wave speed $c$, the flow rate $Q$ and one other geometrical or physical parameter, such as the (half) domain length $L$ or the undisturbed lower fluid thickness $h_0$, while holding the rest of the parameters fixed.

Interfacial wave profiles obtained for different values of $B$ are shown in figure 7 (the parameters $h_0$, $m$ and $C$ are set to the same values as in the previous section). Evidently the interface forms an increasingly pronounced finger as the Bond number increases. This results in a near segregation of the two fluids as the finger almost touches the walls at sufficiently large $B$. This is consistent with the time-dependent calculations of the previous section. In fact the profile for $B=1$ is almost identical to the final wave state in figure 1(f), the only difference being that the thin-film regions near each wall are in the present case flat and steady. The capillary ridge and the depression at the leading edge of the finger are also apparent.

Figure 7.

Interfacial travelling-wave profiles computed using AUTO for various values of $B$. The parameters $h_0$, $m$, $C$ are fixed to the same values as in figure 1.

Figure 8 shows the variation of the travelling-wave speed $c$ and the flow rate $Q$ with the Bond number $B$, for different values of the undisturbed lower fluid thickness $h_0$, capillary number $C$ and viscosity ratio $m$. Strikingly, for sufficiently large Bond number, both $c$ and $Q$ reach a plateau at the same limiting value. Furthermore, this limiting value increases with $m$ but seems to be independent of both $h_0$ and $C$. Reducing $h_0$ delays the approach to the plateau, i.e. the wave speed becomes constant at a larger value of the Bond number. These results may be used to corroborate the comments made in the previous section on finger merging. From figure 8(a) it is clear that $c$ depends strongly on $h_0$ before the plateau is attained. In a domain of fixed length $2L$, increasing $h_0$ increases the volume of the lower fluid leading to a taller finger. According to figure 8, roughly speaking $c$ increases with $h_0$ so that a taller finger travels faster. Once $h_0$ is sufficiently large for the travelling wave to (almost) touch the upper wall, the finger widens as $h_0$ increases but the wave speed remains approximately constant.

Figure 8.

Variation of wave speed $c$ and flow rate $Q$ with Bond number $B$. The various curves shown are obtained for different values of (a,b) undisturbed lower fluid thickness $h_0$, (c,d) capillary number C or (e, f) viscosity ratio m, as quoted in the respective legends (legends refer to both panels in the same row). Unless given in the legends, parameters $h_0$, $m$, $C$ take the same values as in figure 1.

3.4. Direct numerical simulations

Thus far our observations have been based on solutions to the model long-wave evolution equation (2.7). We have also considered values of the Bond number at which, strictly speaking, the model equation is outside of its range of validity (see the derivation in Appendix A). To provide some reassurance that our results are nonetheless relevant to a real flow, in this section we consider the full system (2.1), (2.3) and (A4) with no long-wave approximation. We solve this system numerically using the finite-element object-oriented multi-physics library oomph-lib (Heil & Hazel Reference Heil and Hazel2006a,Reference Heil and Hazelb). The flow set-up is otherwise the same and we again start our simulations from the initial condition (3.1). Some additional details about the numerical implementation in oomph-lib are provided in Appendix E.

Generally speaking the interfacial displacement evolves in time in a manner similar to that predicted by the lubrication model. For sufficiently large Bond number this includes the early finger formation, the ensuing merging dynamics and the subsequent near-wall touch-down/touch-up behaviour. The ultimate fully merged profile comprising a single wide finger is shown in figure 9. Superimposed onto this profile is that found from the travelling-wave solution to the lubrication model (cf. § 3.3) and the large-time profile obtained from the time-dependent simulation of the lubrication model (cf. § 3.1). Despite the moderate size of the Bond number (here $B=1$ and the lubrication model strictly requires $B\ll 1$) there is excellent agreement between all three sets of results. The only minor discrepancies occur in the upper part of the finger which is slightly wider in the DNS profile (figure 9a), and in the thin film near to the lower wall (figure 9c,d) where small interfacial ripples are observed in both sets of time-dependent calculations.

Figure 9.

Comparison of interfacial profiles obtained from time-dependent calculations of the lubrication model, travelling-wave calculations of the same model and DNS of the Navier–Stokes equations. Each solution has been shifted to the right by an appropriate distance for better visualisation. The solution over the whole computational domain in shown in (b) with magnifications thereof illustrated in (a,c,d). Parameter values for $h_0$, $m$, $C$, $B$ are the same as in figure 1.

Inspired by the good agreement seen in figure 9 between model results and DNS, we performed extensive numerical simulations for a range of Bond numbers to establish the predictive power of the model. A comparison between numerical results obtained from time-dependent simulations of the model evolution equation (2.7) and DNS can be found in figure 10. In figure 10(a), the maximum absolute error is shown as a function of the Bond number, found by calculating the difference between the two final-time interfacial profiles. An overall trend of growing error for increasing Bond number is evident, indicating that the model will eventually (for much larger $B$) not be able to capture the ‘true’ dynamics. Figure 10(b) demonstrates that the amplitude of the two solutions is almost identical for all Bond numbers considered. For $B\gtrapprox 0.4$, the interfacial structure almost touches both channel walls and the amplitude approaches the value 1 (cf. figure 8a). This allows us to conclude that the apparent increase in error at higher Bond numbers is mainly due to width disparity between the fingers, similar to that seen in figure 9(a).

Figure 10.

Comparison between numerical results obtained from time-dependent simulations of the model evolution equation (2.7) and DNS, for various Bond numbers. (a) Maximum absolute error found by calculating the difference between the two final-time interfacial profiles (model and DNS). (b) Solution amplitude found by taking the difference between the interfacial maximum and minimum at the final time. The parameters $h_0$, $m$, $C$ are fixed to the same values as in figure 1.

4.1. Transition regions ${{T_R}}$ and ${{T_L}}$

We suppose there exist two transition regions in which (asymptotically) the two walls are connected to each other by a section of the interface which meets them at zero contact angle. In both of these regions the interface has a relatively steep slope. We seek a solution in the form

(4.2a–c)\begin{equation} h = \mathscr{H}(\zeta),\quad {\rm with}\quad \xi = B^{{-}1/2}\zeta \quad \text{in }{{T_R}}; \quad \xi ={-}\mathscr{L} - B^{{-}1/2}\zeta \quad \text{in }{{T_L}}, \end{equation}

where $\mathscr {H}$ is to be found, $\zeta =\textit {O}{\left (1 \right )}$ and $\mathscr {L}$ is the width of region ${F_U}$ (see figure 11). We note that the scaling of $B^{-1/2}\ll 1$ for the $\xi$ coordinate was chosen to retain the effect of surface tension at leading order, by ensuring the balance of the two terms on the right-hand side of (2.8d). In the coordinate rescalings (4.2b,c), the origin of $\xi$ is assumed to be in region ${T_R}$ and hence (4.2c) is essentially a shift of (4.2b) by a distance $\mathscr {L}$. The value of $\mathscr {L}$ can be estimated analytically by noticing that almost all the volume of the lower fluid is included in the finger. The volume of the finger is approximately $\mathscr {L}$ (the dimensionless height of the channel is 1), and this should be equal to the lower fluid volume in the undisturbed state. Therefore at leading order, $\mathscr {L}=2Lh_0$, where $2L$ is the domain length and $h_0$ is the undisturbed lower fluid height.

We assume that $m$, $C$ and $Q$ are of order 1, and that $c=\textit {O}{\left (1 \right )}$ and $\mathscr {C}=o{\left (B^{3/2} \right )}$ (these latter scalings will be confirmed later). Then (4.1) becomes at leading order

(4.3)\begin{equation} \mathscr{H}_\zeta + \mathscr{H}_{\zeta \zeta \zeta} = 0, \end{equation}

noting that $\mathscr {D}$ and $H_2$ are positive for $0<\mathscr {H}<1$. Assuming a solution in which the interface meets each wall with zero contact angle, then

(4.4)\begin{equation} \mathscr{H} = \frac{1}{2}(1 - \sin\zeta), \quad \text{over } -\frac{\rm \pi}{2}<\zeta<\frac{\rm \pi}{2}. \end{equation}

For later reference, it is important to note that

(4.5a)$$\begin{align} \mathscr{H} &\sim \frac{1}{4}\left(\zeta - \frac{\rm \pi}{2}\right)^2, \quad \text{as } \zeta \to \frac{\rm \pi}{2}, \end{align}$$

(4.5b)$$\begin{align}\mathscr{H} &\sim 1 - \frac{1}{4}\left(\zeta + \frac{\rm \pi}{2}\right)^2, \quad \text{as } \zeta \to -\frac{\rm \pi}{2}. \end{align}$$

In both transition regions the limit $\zeta \to \pm {\rm \pi}/2$ corresponds to approaching the film near the lower/upper wall, respectively.

4.2. Flat-film regions ${F_D}$ and ${F_U}$

In the flat thin-film regions against the lower/upper walls (these regions lie outside/inside the finger), the film is of approximately uniform thickness and we take, respectively,

(4.6a,b)\begin{equation} h = B^{{-}1}H_{D} \quad \text{or} \quad h = 1 - B^{{-}1}H_{U}, \end{equation}

where $H_{D}$ and $H_{U}$ are both constant. The scaling of $B^{-1}$ is chosen for matching purposes as we will see later. In region ${F_D}$, (2.7b) implies $q\sim B^{-2}$, and so the right-hand side dominates in (4.1) and the flat-film thickness is found to be

(4.7)\begin{equation} H _{D} = \frac{\hat{\mathscr{C}}}{c}, \end{equation}

where we have defined

(4.8)\begin{equation} \mathscr{C}=B^{{-}1}\hat{\mathscr{C}}. \end{equation}

In region ${F_U}$, we expand $q$ in terms of $B^{-1}$, apply the scaling (4.8) in (4.1) and disregard terms of $\textit {O}{\left (B^{-2} \right )}$, resulting in

(4.9)\begin{equation} Q - B^{{-}1}H_{U} = c - B^{{-}1}\left( cH_{U} + \hat{\mathscr{C}} \right). \end{equation}

Hence the appropriate balance at $\textit {O}{\left (1 \right )}$ and $\textit {O}{\left (B^{-1} \right )}$ gives, respectively,

(4.10a,b)\begin{equation} Q = c \quad \text{and} \quad H_{U} = \frac{\hat{\mathscr{C}}}{1-c} = \frac{c}{(1-c)}H_{D}, \end{equation}

on using (4.7). The leading-order relationship between $c$ and $Q$ in (4.10a) has already been noted in the results of figure 8.

4.3. Connection regions ${C_{RD}}$ and ${C_{RU}}$

Here we write

(4.11a–c)\begin{equation} h = B^{{-}1}\tilde H(z) \quad \text{or} \quad h = 1-B^{{-}1}\hat H(z), \quad \text{and} \quad \xi ={\pm} \left( B^{{-}1/2}\frac{\rm \pi}{2} + B^{{-}1}z \right), \end{equation}

where (4.11a,b) are introduced in regions ${C_{RD}}/{C_{RU}}$ (near the lower/upper walls), respectively, and the plus/minus signs in (4.11c) apply in the respective regions ${C_{RD}}/{C_{RU}}$.

Taking Taylor expansions with respect to $B^{-1}$ in the various functions used to define the flux $q$ in (2.7b), and substituting in (4.1), provides the leading-order balance in each connection region. In particular, accounting for (4.8), in region ${C_{RD}}$ we find at $\textit {O}{\left (B^{-1} \right )}$

(4.12a)\begin{equation} \frac{1}{3C}\tilde H^3 \tilde H_{zzz} = c \tilde H - \hat{\mathscr{C}}, \end{equation}

and in region ${C_{RU}}$ we recover (4.10a) at $\textit {O}{\left (1 \right )}$, while at $\textit {O}{\left (B^{-1} \right )}$ we get

(4.12b)\begin{equation} \frac{1}{3mC}\hat H^3\hat H_{zzz} = (1-c)\hat H - \hat{\mathscr{C}}. \end{equation}

Making use of (4.7) and (4.10b), we rescale by writing $\{\tilde H,\hat H\} = \{H_{D},H_{U}\} H$ and $z = \beta _{\{{D},{U}\}} Z$, where $\beta _{D}=H_{D}/(3Cc)^{1/3}$ and $\beta _{U}=H_{U}/(3mC(1-c))^{1/3}$. Then (4.12) both reduce to the Landau–Levich equation (Landau & Levich Reference Landau and Levich1942, see also Bretherton Reference Bretherton1961)

(4.13)\begin{equation} H^{'''} = \frac{H-1}{H^3}, \end{equation}

where a prime denotes a derivative with respect to $Z$.

The far-field conditions follow by matching to the solution in transition region ${T_R}$ and to the flat-film solution at the walls. Reconciling the scalings (4.2b) and (4.11c), and using (4.5) and (4.6), the required matching conditions in regions ${C_{RD}}/{C_{RU}}$ are

(4.14a,b)\begin{equation} H \sim \frac{1}{2}\mathscr{A}_{\{{D},{U}\}} Z^2 \quad \text{as } Z\to-\infty \quad \text{and} \quad H \sim 1 \quad \text{as } Z\to \infty, \end{equation}

where

(4.15a,b)\begin{equation} \mathscr{A}_{D} = \frac{\beta_{D}^2}{2H_{D}} \quad \text{and} \quad \mathscr{A}_{U} = \frac{\beta_{U}^2}{2H_{U}}. \end{equation}

The first and second conditions in (4.14a,b) effect the matches with the transition region ${T_R}$ and with the flat thin-film regions ${F_D}$, ${F_U}$, respectively. The connecting region problems near the right transition region thus require the solution of the Landau–Levich equation (4.13) subject to the matching conditions (4.14a,b). This problem appears in a number of related studies including the dragging of a liquid by a moving plate (Landau & Levich Reference Landau and Levich1942), the motion of long bubbles in horizontal tubes (Bretherton Reference Bretherton1961) and the capillary draining of annular films in horizontal tubes (Hammond Reference Hammond1983; Lister et al. Reference Lister, Rallison, King, Cummings and Jensen2006) and vertical tubes (Jensen Reference Jensen2000). It has a one-parameter family of solutions, meaning that $\mathscr {A}_{\{{D},{U}\}}$ cannot be determined at this stage.

4.4. Connection regions ${C_{LD}}$ and ${C_{LU}}$

Working as in the previous section, we write

(4.16a–c)\begin{align} h = B^{{-}1}\tilde{\mathscr{H}}(z) \quad \text{or} \quad h = 1-B^{{-}1}\hat{\mathscr{H}}(z), \quad \text{and} \quad \xi ={-}\mathscr{L} \pm \left({-}B^{{-}1/2}\frac{\rm \pi}{2} + B^{{-}1}z \right), \end{align}

where $\mathscr {L}$ was introduced in (4.2c). The scalings (4.16a,b) are applied in regions ${C_{LD}}/{C_{LU}}$ (near the lower/upper walls), respectively, and the plus/minus signs in (4.16c) apply in the respective regions ${C_{LD}}/{C_{LU}}$.

Introducing the same scalings as before, $\{\mathscr {H},\tilde{\mathscr{H}}\} = \{H_{D},H_{U}\} H$ and $z = \beta _{\{{D},{U}\}} Z$, we again obtain the Landau–Levich equation (4.13). The matching conditions are

(4.17a,b)\begin{equation} H \sim \frac{1}{2}\mathscr{A}_{\{{D},{U}\}} Z^2 \quad \text{as } Z\to\infty \quad \text{and} \quad H \sim 1 \quad \text{as } Z\to-\infty. \end{equation}

The asymptotic problem given by (4.13) and (4.17a,b) admits a unique solution and therefore $\mathscr {A}_{D}=\mathscr {A}_{U}$. We find numerically that these constants are equal to $0.643$, which agrees with the value quoted elsewhere in the literature (e.g. Bretherton Reference Bretherton1961; Jensen Reference Jensen2000; Lister et al. Reference Lister, Rallison, King, Cummings and Jensen2006).

4.5. Comparison with numerical solutions

Since $\mathscr {A}_{D}=\mathscr {A}_{U}$, it follows from (4.15a,b) that an explicit formula can be obtained for the wave speed at leading order (using also the definitions of $H_{D}$ and $H_{U}$ from (4.7) and (4.10b)), given by

(4.18)\begin{equation} c = \frac{m^{2/5}}{1+m^{2/5}}. \end{equation}

It is surprising that $c$ depends only on the viscosity ratio, but this is in line with the results presented in § 3.3, in particular figure 8. The relationship (4.18) is plotted in figure 12. The wave speed varies rapidly with $m$ when $m<1$ and the lower fluid is the more viscous. Note that $0< c<1$ and so the wave speed is less than the wall speed (which has been used to non-dimensionalise velocities); in particular $c=1/2$ when $m=1/2$ so that the wave propagates at half the wall speed for equal-viscosity fluids. These remarks agree with the results in figure 8.

Figure 12.

Variation of the interfacial wave speed $c$ against the viscosity ratio $m$, according to the asymptotic result in (4.18).

Since $\mathscr {A}_{D}=\mathscr {A}_{U}=0.643$, it follows that the as yet undetermined constant $\hat{\mathscr{C}}$ is given by

(4.19)\begin{equation} \hat{\mathscr{C}} = 2.675 \, C^{2/3} c^{5/3}, \end{equation}

where $c$ is given by (4.18). This result may be combined with (4.6), (4.7) and (4.10b) to yield an estimate for the thin-film thickness in regions ${F_D}$, ${F_U}$ for a particular value of $B$.

We conclude this section by comparing the predictions of the present asymptotic analysis with numerical solutions of the lubrication model and of the full equations via DNS. To facilitate this comparison, we fix the parameters to the same values that were used in the numerical computations in § 3, namely $h_0=0.2$, $m=0.5$, $C=10^{-3}$, and we take a moderate value for the Bond number such as $B=1$. Setting $m=0.5$ in (4.18) yields the leading-order prediction for the wave speed $c=0.43$. The numerically computed speed for the lubrication model is approximately $c=0.4725$; this was obtained from a time-dependent simulation of the lubrication model (2.7) and it was confirmed by a travelling-wave calculation as described in § 3.3; see, in particular, figure 8. The corresponding wave speed obtained from DNS is approximately $c_{DNS}=0.4723$. Once again we note that the lubrication model agrees well with full DNS even for Bond numbers outside its strict range of validity.

Using these results and taking the same parameter values in (4.19), (4.7), (4.10b) and then (4.6), we find the individual film thicknesses to be $h_{FD}=0.015$ and $h_{FU}=0.012$ near the down/up walls, respectively, while the numerically computed values from the lubrication model are 0.0176 and 0.0126, respectively, and the corresponding DNS values are approximately 0.0175 and 0.0127. Clearly the results obtained from simulations of the lubrication model and DNS are almost identical, and they are in good agreement with the asymptotic results.

Finally we compare the numerical solution of the asymptotic problem in each connection region with that obtained via a time-dependent simulation of the lubrication model. Numerically solving the Landau–Levich equation (4.13) with far-field conditions (4.14a,b) or (4.17a,b) gives the solutions shown with the red curves in figure 13 (the results have been scaled back to the original variables). The solution obtained from time-dependent simulations is superimposed with black curves (this is the same profile shown in figure 1f), providing good agreement between the two sets of results.

Figure 13.

Solution to the asymptotic problem in the (a) left and (b) right connection regions, plotted for the same parameter values used in previous figures ($h_0=0.2$, $m=0.5$, $C=10^{-3}$, $B=1$), and shown with red dashed curves. The black solid curves correspond to the time-dependent results. The asymptotic results have been rescaled back to the original variables, $h=h(\xi )$, by using the appropriate values for $\beta _{\{{D},{U}\}}$ and scalings in (4.11), and have been also shifted horizontally in order to match the spatial location of the connection regions as seen in figure 1( f).