2.1. Governing equations

The governing equations for the flow of each fluid are the cylindrical axisymmetric Navier–Stokes equations,

where $i=(a)$ corresponds to the annular film and $i=(c)$ to the core fluid. Here, $(\bar{u},\bar{w})$ are the velocity components in the $(\bar{r},\bar{z})$ directions respectively, $\bar{p}$ is pressure, ${\it\rho}$ is density, ${\it\mu}$ is molecular viscosity and $g$ is gravity. Subscripts denote partial derivatives. The boundary conditions at the wall, $\bar{r}=\bar{a}$ , are no-slip,

while at the free surface, $\bar{r}=\bar{R}(\bar{z},\bar{t})$ , the boundary conditions

are enforced, corresponding to the continuity of velocity, continuity of tangential stress and jump in normal stress according to the Young–Laplace law respectively. Here, ${\it\gamma}$ denotes the surface tension. In addition, at the free surface the kinematic boundary condition,

must also be satisfied. Of course, physical significance requires that at the tube centreline the velocity field be finite, which provides an additional constraint for the core flow solutions.

The flow governed by the system (2.1)–(2.8) always admits an exact solution corresponding to flow with a flat interface, $\bar{R}=\bar{R}_{0}$ ,

where the only non-zero velocity component is the axial one, $\bar{w}$ , and the only non-zero gradient is the radial one, $\partial _{\bar{r}}$ . Here, the normal stress boundary condition (2.7) reduces to $\bar{p}_{\bar{z}}^{(a)}=\bar{p}_{\bar{z}}^{(c)}\equiv \bar{p}_{\bar{z}}$ . This parallel-flow flat solution is a balance of viscous and gravitational forces and the imposed pressure gradient driving the core flow; an example of velocity profiles given by (2.9) and (2.10) is shown in figure 1(b). It should be noted that if $\bar{R}_{0}=0$ , corresponding to the tube being entirely filled with the annular fluid, (2.9) reduces to the profile of Poiseuille pipe flow; likewise, if one sets $\bar{R}_{0}=\bar{a}$ , (2.10) reduces to that of Poiseuille flow of the core fluid.

The stability of the exact solution given by (2.9) and (2.10) was studied by Hickox (Reference Hickox1971) for various values of the parameters of the problem (see also Yih Reference Yih1967 for an analysis with a similar set-up without surface tension). In order to study the evolution of these interfacial disturbances and the flow of the annular film in the simplest possible setting, we turn to a model of the flow developed using long-wave asymptotics.

2.2. Fully coupled long-wave equations

The dimensional parameters ${\it\mu}^{i}$ , ${\it\rho}^{i}$ , $g$ , ${\it\gamma}$ , $\bar{Q}^{(c)}$ , $\bar{a}$ , $\bar{R}_{0}$ and $\bar{{\it\lambda}}$ , where $\bar{{\it\lambda}}$ denotes the length of a typical interfacial disturbance, combine to form seven dimensionless parameters: (i) Reynolds number Re, (ii) Froude number Fr, (iii) capillary number $C$ , (iv) viscosity ratio $m$ , (v) density difference ratio  $l$ , (vi) film thickness parameter $a$ and (vii) aspect ratio ${\it\epsilon}$ , defined by

The dynamics of the film flow is uniquely determined by the values of these seven dimensionless parameters. We note that other parameter choices are also valid. For instance, (2.11) corresponds to a viewpoint where the core volume flux $\bar{Q}^{(c)}$ and mean film thickness $\bar{R}_{0}$ are prescribed, and the annular volume flux $\bar{Q}^{(a)}$ and pressure gradient respond accordingly. One could replace the aspect ratio $a$ with a flux-ratio parameter $Q^{(a)}=\bar{Q}^{(a)}/\bar{Q}^{(c)}$ , corresponding to a viewpoint where both fluid fluxes are prescribed and the mean film thickness and pressure gradients respond accordingly. Moreover, in other studies of core–annular flows the Weber number $\mathit{We}=C\times \mathit{Re}$ is sometimes used in place of the capillary number. Finally, we note that with the density ratio defined as in (2.11), the case where equal density fluids occupy the core and annular regions is defined by $l=0$ .

System (2.1)–(2.8) can be made dimensionless by typical scales

where $\overline{W}_{0}=2\bar{Q}^{(c)}/{\rm\pi}\bar{R}_{0}^{2}$ is the core flow centreline velocity and $\overline{U}_{0}\equiv {\it\epsilon}\overline{W}_{0}$ . After substitution of (2.12) into (2.1)–(2.8), we consider first the leading order in an asymptotic expansion for the limit of the long-wave parameter ${\it\epsilon}\rightarrow 0$ . In addition to small ${\it\epsilon}$ , we also restrict the values of $\mathit{Re}$ by requiring that in the long-wave limit the limit ${\it\epsilon}\mathit{Re}\rightarrow 0$ as well. Under these assumptions, the motion equations (2.1)–(2.3) and boundary conditions (2.4)–(2.8) made dimensionless by (2.12) yield a leading-order velocity profile $w_{LW}^{(a)}$ and $w_{LW}^{(c)}$ given in appendix A by (A 1)–(A 2); this solution is a slowly varying $z$ -dependent modulation of the exact parallel flow (2.9)–(2.10). The slowly varying annular pressure is now related to that of the core fluid by the asymptotic normal stress boundary condition (2.7)

at the long-wave interface $r=R(z,t)$ . As is customary, while the second term  $R_{zz}$ in the capillary tension force balance expressed by (2.13) is formally of higher order, it is retained to account for uniform matched-asymptotics validity of the long-wave expansion (see, e.g., Oron et al. Reference Oron, Davis and Bankoff1997).

Next, we seek approximate expressions for the pressure gradients $p_{z}^{i}$ . These are found through the dimensionless volume fluxes

which are given in expressions (A 3)–(A 4) in appendix A. Incompressibility of both core and annular fluids implies that the total volume flux $q^{(a)}+q^{(c)}$ through the tube is independent of the spatial location $z$ . In the case where the volume flux $\bar{Q}^{i}$ of each fluid is prescribed at one end of the tube, the total flux anywhere in the tube is thus

where the $Q^{i}$ are the prescribed fluxes made dimensionless by $\bar{W}_{0}\bar{R}_{0}^{2}$ . Then, with the free surface boundary condition (2.13), equations (A 3)–(A 4) determine the pressure terms $p_{z}^{(a)}$ and $p_{z}^{(c)}$ as functions of $R$ and its spatial derivatives; $p_{z}^{(a)}$ is given in (A 5) in appendix A and $p_{z}^{(c)}$ is recovered through (2.13).

An evolution equation for $R(z,t)$ follows from volume conservation of either the core or the annular fluid and the kinematic boundary condition (2.8). Focusing for definiteness on the annular fluid, this equation can be expressed in terms of the annular flux as

2.3. Decoupled long-wave models

Various model equations can be derived via different asymptotic balances between the parameters $m$ , $l$ , $a$ , $C$ and $\mathit{Re}/\mathit{Fr}^{2}$ , and the free surface location  $R$ . The model equation derived in Camassa et al. (Reference Camassa, Forest, Lee, Ogrosky and Olander2012) is recovered when the additional asymptotic orderings $m\gg (a/\tilde{R})^{4}$ and $m\gg l(a/\tilde{R})^{4}$ are satisfied, where $\tilde{R}$ is a characteristic value for the core radius, and when $C=O(1)$ and $\mathit{Fr}^{2}=O(\mathit{Re})$ . These asymptotic relations are realized by a high annular-to-core viscosity ratio and a core region that is not exceedingly small in volume (see appendix for a brief discussion of a model appropriate under a different asymptotic balance corresponding to ‘very thick-film’ problems studied elsewhere), and imply that $q^{(a)}\ll q^{(c)}$ so that the core fluid is moving much faster than the annular fluid. This difference in time scales for the two fluids allows the annular film flow to be decoupled from the core flow, in the sense that the core flow experiences the free surface essentially as a rigid wavy wall.

Under these assumptions, the pressure gradient $p_{z}^{(a)}$ simplifies to (A 8). Substitution of this approximate $p_{z}^{(a)}$ into (2.16) yields the model equation in Camassa et al. (Reference Camassa, Forest, Lee, Ogrosky and Olander2012) after a rescaling in time $t^{\prime }=t/m$ (and immediately dropping the prime for ease of notation of the time derivative),

where the functions $f_{n}$ , $n=1,2,3$ , are given by

and the parameters $F$ and $S$ by

Here, $R(z,t)=\bar{R}/\bar{R}_{0}$ is the dimensionless radial location of the interface as a function of axial location and time, where $\bar{R}_{0}$ is the mean radius of the core region, and $a=\bar{a}/\bar{R}_{0}$ is the dimensionless tube radius. The parameter $F$ represents the ratio of the free surface stress created by the core flow to gravity and $S$ is the inverse of a capillary number. We have also adopted the convention of returning to the original aspect ratio by rescaling $z$ and $t$ by ${\it\epsilon}$ as in Camassa et al. (Reference Camassa, Forest, Lee, Ogrosky and Olander2012). This results in ${\it\epsilon}$ being scaled out of (2.17), but the validity of the model equation still depends on the condition ${\it\epsilon}\ll 1$ being satisfied.

Equation (2.17) is strongly nonlinear, i.e. no smallness assumption on the free surface displacement from a reference flat configuration has been made. Moreover, the model retains the cylindrical geometry of the problem, as evidenced by the inverse powers and logarithms of  $R$ .

The dynamics of the governing equations was determined by seven dimensionless parameters, while the dynamics of (2.17) is determined by three parameters: $F$ , which represents the competition between viscous forces at the interface and gravity forces; $S$ , which keeps track of the competition between viscous forces at the interface and stresses from capillary forces; and $a$ , the film thickness parameter. Both $S$ and $F$ are dependent on the thickness of the film. All three parameters play a role in determining whether nonlinear disturbances to the free surface travel upwards or downwards. As we will see, the weakly nonlinear regime is somewhat simpler in that the direction of disturbance propagation is determined only by $F$ and $a$ . Specifically, if $F<f_{1}(1;a)/f_{2}(1;a)$ disturbances travel up the tube, i.e., the free surface stress caused by the core flow dominates the effects of gravity, otherwise disturbances travel downward due to gravity overpowering the effects of the core flow.

It should be noticed that $F$ and $S$ depend on the film thickness parameter $a$ through the Reynolds, Froude and capillary numbers as a result of the choice of length scales used to define these parameters. An alternative choice for parameters which would remove the $a$ dependence from $F$ and $S$ and incorporate it into the functions $f_{n}$ may be preferable, i.e. define new coefficients $\tilde{F}=a^{4}F$ and $\tilde{S}=a^{2}S$ which would correspond to $\tilde{f}_{2}=a^{-4}f_{2}$ and $\tilde{f}_{3}=a^{-2}f_{3}$ . This choice makes the coefficients $F$ and $S$ independent of film thickness initial data and thus constant across different realizations of an experiment with fixed ‘hardware’ set-up ${\it\mu}^{i}$ , ${\it\rho}^{i}$ , $\bar{Q}^{(c)}$ , $\bar{a}$ and ${\it\gamma}$ . However, to be consistent with previous literature, we follow the definitions given in (2.19) for the rest of this work.

In the case of equal densities of the core and annular fluid (or negligible gravity), i.e. $F=0$ , (2.17) reduces to the model equation derived by Camassa & Lee (Reference Camassa, Lee, Chwang, Teng and Valentine2006) up to a choice of scalings:

For ease of reference, it is worth labelling these two models with easily identifiable tags. Thus, in the following we will refer to the long-wave models (2.17) and (2.20) as (LW-PG) and (LW-P) respectively.

2.4. Thin-film models

If the liquid film is thin, further simplifications of (LW-PG)/(2.17) and (LW-P)/(2.20) can be expected. We recall the non-dimensional scaling $a=\bar{a}/\bar{R}_{0}$ and $R=\bar{R}/\bar{R}_{0}$ : the limit of thin films corresponds to $\bar{R}_{0}\rightarrow \bar{a}$ , so that $a$ limits to 1 from above and the spatial mean of the non-dimensional interface position $R$ also limits to unity. The thin-film limit asymptotics can then be described by a thinness small parameter $0<{\it\beta}\ll 1$ , where $a=1+{\it\beta}$ , and by shifting to a new dependent variable $h(z,t)$ , with spatial mean 1, defined by

so that ${\it\beta}h$ represents the non-dimensional height of the interface measured from the tube wall (i.e. the thickness of the annular film). Then, (LW-PG)/(2.17) can be rewritten in terms of $h$ , and each of the functions $f_{n}$ , $n=1,2,3$ , in (LW-PG) can be expanded about ${\it\beta}=0$ :

Notice that the order of the gravity term is $f_{2}=O({\it\beta}^{2})$ , while the order of the core flow term is $f_{1}=O({\it\beta})$ . Thus, asymptotically the effects of tangential stress at the interface dominate the effects of gravity in the thin-film limit. Notice also that the effects of surface tension are of order $O({\it\beta}^{3})$ ; these effects are again typically retained despite their higher order for reasons discussed in, e.g., Camassa et al. (Reference Camassa, Ogrosky and Olander2014). Substituting the expanded $f_{n}$ into the long-wave model (LW-PG) and retaining up to $O({\it\beta}^{2})$ hyperbolic terms (while retaining the $O({\it\beta}^{3})$ term of $f_{3}$ to keep the effects of surface tension) gives

where $\tilde{F}={\it\beta}(11-2F)/(2-6{\it\beta})$ and $\tilde{S}^{\ast }=8{\it\beta}^{2}S/[3(1-3{\it\beta})]$ , and time has been rescaled by $t^{\prime }=(2{\it\beta}-6{\it\beta}^{2})t$ (again having dropped primes).

If the densities of the two fluids are equal or the effects of gravity are negligible, i.e.  $F=0$ , then substitution of (2.22) into (2.17) and retention of only the $O({\it\beta})$ hyperbolic terms results in the standard form of the thin-film model for pressure-driven flow (after a rescaling in time by $t^{\prime }=2{\it\beta}t$ , dropping primes hereafter),

where $\tilde{S}=8{\it\beta}^{2}S/3$ . Equation (2.24) is the thin-film counterpart to the long-wave model (LW-P)/(2.20) and was derived directly from the governing equations and studied numerically by Kerchman (Reference Kerchman1995). For models (2.23) and (2.24), we will use, respectively, the tags (TF-PG) and (TF-P).

While a study of the general dynamical features of all the above models has to rely in large part on their numerical simulations, much insight can be gained by analysing special solutions such as travelling waves, or by studying transient growth of perturbations from exact simple solutions such as those corresponding to a flat interface, thereby determining their linear stability. We will use all these tools, starting with linearized dynamics, next.

3.1. Pressure-driven flow

Introduction of the perturbation (3.1) into the pressure-driven model equations (LW-P) and (TF-P) yields the dispersion relations

(It should be noted that (3.3) corresponds to (TF-P) prior to a rescaling of time by  $2{\it\beta}$ .)

The maximum growth rate of each model occurs for disturbances with wavenumber $k_{m}=1/\sqrt{2}$ , corresponding to a (dimensional) wavelength of $\bar{{\it\lambda}}_{m}=2{\rm\pi}\sqrt{2}\bar{R}_{0}$ , indicating that the typical wavelength of disturbances is set by the radius of the interface, in analogy to the classical result for jets due to Rayleigh. The maximum growth rate for each model is given by

Figure 2(a) shows the growth rates of each model for two values of  $a$ . The models are in close agreement with one another, with model (TF-P) producing slightly smaller growth rates. The real part of each dispersion relation describes the propagation speed of infinitesimal interfacial disturbances. Figure 2(b) shows the phase speeds for each model as a function of the film thickness parameter  $a$ ; disturbances in the thin-film model travel slower than those in the long-wave model.

Figure 2. (a) Growth rates of (LW-P) (solid lines) and (TF-P) (dashed lines) for two values of  $a$ . (b) Speed of linear perturbations as a function of  $a$ .

We remark that both dispersion relations (3.2) and (3.3) are, to within a choice of scaling, identical to the well-known Kuramoto–Sivashinsky (KS) equation (see, e.g., Kuramoto & Tsuzuki Reference Kuramoto and Tsuzuki1980, Sivashinsky & Michelson Reference Sivashinsky and Michelson1980 or Kerchman Reference Kerchman1995). Variants of the KS equation have been used to model weakly nonlinear thin films, and solutions to the KS equation are known to exhibit spatio-temporal chaos; we return briefly to this in § 4.

Since solutions to both models are unstable to a range of small wavenumbers $k$ , we now examine whether solutions to each model are convectively or absolutely unstable. This classification of instabilities is borrowed from the literature concerning jets, and has been studied in the case of gravity-driven film flows which coat tubes, e.g. Duprat et al. (Reference Duprat, Ruyer-Quil, Kalliadasis and Giorgiutti-Dauphine2007) for exterior gravity-driven liquid coatings and Camassa et al. (Reference Camassa, Ogrosky and Olander2014) for interior gravity-driven liquid coatings; see also Vellingiri, Tseluiko & Kalliadasis (Reference Vellingiri, Tseluiko and Kalliadasis2015) for a study of absolute and convective instabilities in the related problem of gas–liquid film flows in a channel geometry. It was shown by Duprat et al. (Reference Duprat, Ruyer-Quil, Kalliadasis and Giorgiutti-Dauphine2007) that there is a simple criterion for determining whether solutions to (5.3) are absolutely or convectively unstable using the method developed by Briggs (Reference Briggs1964) and Bers (Reference Bers, Rosenbluth and Sagdeev1983), and expanded by Huerre & Monkewitz (Reference Huerre and Monkewitz1990).

The criterion for absolute instability in the Frenkel model (5.3) (Frenkel Reference Frenkel1992; also see § 5 for further discussion of its relation to the core–annular models) was found by Duprat et al. (Reference Duprat, Ruyer-Quil, Kalliadasis and Giorgiutti-Dauphine2007) to be

Since the dispersion relation for the Frenkel equation,

can be transformed into (3.3) by the transformation ${\it\omega}_{F}=-2(a-1){\it\omega}_{TF\text{-}P}$ and $S^{\prime }=-(a-1)S$ , we find that interfacial perturbations to model (TF-P) are absolutely unstable if the following condition is satisfied (after a sign change to account for the different direction of disturbance propagation):

Alternatively, this condition can be expressed as a condition on the dimensional film thickness:

If condition (3.9) is not satisfied, the free surface is convectively unstable.

The thin-film dispersion relation (3.3) can be transformed into (3.2) by the transformation

Thus, the condition for interfacial disturbances to be absolutely unstable in the long-wave model is (immediately dropping asterisks)

Figure 3 shows the regions of $S$ – $a$ parameter space for which interfacial disturbances are absolutely unstable for each model. While both models are convectively unstable for small $S$ and absolutely unstable for large $S$ , there is an intermediate range of $S$ values for which disturbances in the thin-film model are absolutely unstable while disturbances in the long-wave model are only convectively unstable. This is perhaps to be expected from the results in figure 2. The thin-film model’s growth rate is slower than its long-wave counterpart (recall that $a>1$ for fluid coating the interior of a tube), which intuitively would decrease the region of parameter space in which absolute instabilities are found. However, figure 2(b) shows that the phase speed for the thin-film case is smaller than its long-wave counterpart, and this could be expected to have the opposite effect of increasing the region of absolute instability for model (TF-P). The results in figure 3 indicate that the lower phase speed is the dominant factor in determining the nature of instabilities; this also agrees with the relatively large difference in phase speeds and small difference in growth rates between the models, as seen in figure 2, for values of the radius $a$ sufficiently larger than unity.

Figure 3. Shaded regions of $S$ – $a$ parameter space indicating interfacial disturbances that are absolutely unstable for (TF-P) only (light grey) or both (LW-P) and (TF-P) (dark grey).

3.2. Full equations

We now briefly repeat the linear stability analysis for the full model equations. Introduction of perturbations (3.1) into the long-wave and thin-film models (LW-PG) and (TF-PG) respectively yields the dispersion relations

The maximum growth rates of each model are identical to (3.4) and (3.5); the wavelengths $k_{m}$ of maximum growth rate are also unchanged by the addition of the gravity term. The influence of gravity is seen solely in the speed of disturbances through the real part of each dispersion relation. Figure 4(a) shows the speed of disturbances in each model as a function of $a$ for two values of  $F$ . For each model there is a critical value $F_{c}(a)$ for which the effects of gravity and surface stress balance one another so that disturbances remain in a fixed location as long as their growth remains in the linear regime. For $F>F_{c}$ (gravity stronger than surface stress), linear disturbances move down the tube, and for $F<F_{c}$ (gravity weaker), disturbances move up the tube. For the long-wave and thin-film models these critical values are, respectively, $F_{c_{LW}}=(a^{4}-a^{2})/(a^{2}-1-2\log a)$ and $F_{c_{TF}}=(5a-3)/(2a-2)$ (these values are shown in figure 4 b). It should be noted that the thin-film model once again predicts lower velocities than the long-wave model, which is consistent with both the pressure-driven subcase of speed in model (TF-P) being lower than its long-wave counterpart and the gravity-driven subcase where the appropriate model (5.3) predicts higher (downward) speed than its long-wave counterpart as in Camassa et al. (Reference Camassa, Ogrosky and Olander2014).

Figure 4. (a) Speed of an infinitesimal disturbance to (LW-PG) (solid line) and (TF-PG) (dashed line) for two values of $F$ ; (b)  $F=F_{c_{LW}}$ (solid line) and $F=F_{c_{TF}}$ (dashed line) for (LW-PG) and (TF-PG) respectively, distinguishing between upwards and downwards flow of linear disturbances.

The dispersion relation (3.3) can be transformed into (3.13) by a similar transformation to those in § 3.1. The condition for interfacial disturbances to be absolutely unstable in the thin-film model (TF-PG) is then

The absolute value bars in (3.14) are necessary since, for $F>2.5$ , there is a range of large enough $a$ for which the speed of instabilities will be negative, due to the effects of gravity being stronger than the interfacial stress from the core flow. For given $F$ in the range $F>2.5$ there is exactly one value of $a$ , specifically $a=(2F-3)/(2F-5)$ , for which the right-hand side of (3.14) vanishes, indicating that disturbances are absolutely unstable for any $S$ ; this is due to the effects of gravity and interfacial stress exactly cancelling one another.

The long-wave model’s dispersion relation (3.12) can be recovered from (3.13) by yet another transformation; this leads to the following condition for absolute instability in the long-wave model (LW-PG):

Figure 5 shows the regions in $S$ – $a$ parameter space for which disturbances are absolutely unstable for four values of $F$ . For small values of $F$ , the long-wave model predicts convective instability for more parameter combinations than the thin-film model. For larger values of $F$ , there is a range of $a$ in which the long-wave model predicts absolute instability for more parameter combinations than the thin-film model. Unlike the thin-film model, for every $F$ in the range $F>6.391$ there are exactly two values of $a$ for which any value of $S$ results in absolute instability.

Figure 5. Shaded regions of $S$ – $a$ parameter space indicating interfacial disturbances that are absolutely unstable for only one of (light grey) or both (dark grey) (LW-PG) and (TF-PG) for four different values of  $F$ (1.0 (a), 2.5 (b), 5.0 (c), 7.5 (d)).

4.1. Transient solutions

We first seek solutions to the model equations (TF-P) and (LW-P) with periodic boundary conditions. Numerical integration of each equation is achieved by a pseudospectral method, whereby spatial derivatives and nonlinearities are calculated in Fourier and physical space respectively. A simple second-order predictor–corrector scheme for time integration is used. Due to the nonlinearities in all the models, the Fourier coefficients are appropriately dealiased. The conserved quantities $R^{2}$ or $h$ are monitored throughout the simulations. For initial conditions, the flat free surface is perturbed with either a single mode or multiple modes, i.e.  $h(z,0)=1+\sum _{k=1}^{n}{\it\delta}_{k}\cos (2{\rm\pi}kz+{\it\alpha}_{k})$ , where ${\it\alpha}_{k}$ is a random phase shift for each mode. See Camassa et al. (Reference Camassa, Forest, Lee, Ogrosky and Olander2012), Ogrosky (Reference Ogrosky2013) or Camassa et al. (Reference Camassa, Ogrosky and Olander2014) for more details on the numerical simulations.

Figure 6(a–c) contains series of snapshots showing the evolution of a slightly perturbed interface for long-wave model (LW-P) with various values of $S$ and $a$ ; figure 6(a) shows a relatively thick film ( $a=1.25$ ) and figure 6(c) shows a much thinner film ( $a=1.05$ ). The values of $S$ were chosen in such a way that all three correspond to the thin-film case where $\tilde{S}=0.3$ in order to make a direct comparison with the dynamics of solutions to equation (TF-P). In addition, the initial conditions were identical, i.e. each ${\it\delta}_{k}$ and ${\it\alpha}_{k}$ used in one simulation was identical to those used in every other simulation; the scaling for the initial condition $h(z,0)$ is shown to the right of figure 6(d). Instabilities propagate up the tube, shown here from left to right; it should be noted that successive snapshots are displayed here in a frame of reference moving with the speed of the undisturbed interface. Unlike the case of gravity-driven flow studied in Camassa et al. (Reference Camassa, Ogrosky and Olander2014), in no simulations of (LW-P) did we observe unchecked amplitude growth in finite time (corresponding to ‘choking’ of the core fluid). This highlights the coupling between the core flow and the annular flow through the stress at the free surface; the effect of a growing instability on the hyperbolic term (through the increase in free surface stress) prevents blowup. This is in line with behaviour seen in experiments conducted in Camassa et al. (Reference Camassa, Forest, Lee, Ogrosky and Olander2012): at high core volume fluxes $\bar{Q}^{(c)}$ the rapid increase in the speed of the core fluid in order to pass through a constricted opening prevents plug formation. We also remark that the slight depression in front of each wave, and the corresponding fore–aft asymmetry, is noticeable in all the wave forms in our simulations. This feature is shared among wave solutions to long-wave, thin-film, flat-film and integral models, and it is even observed in simulations of propagating liquid plugs; see, e.g., Campana et al. (Reference Campana, Ubal, Giavedoni and Saita2007) and Ubal et al. (Reference Ubal2008).

Figure 6. (a–c) Time snapshots showing the evolution of the free surface for (LW-P) in a periodic domain with (a)  $S=1.8$ and $a=1.25$ , (b)  $S=11.25$ and $a=1.1$ , and (c)  $S=45$ and $a=1.05$ . Free surface profiles are shown successively shifted at equal time intervals. Profiles are shown in the frame of reference moving with an undisturbed interface. (d) Solution to (TF-P) is shown for $\tilde{S}=0.3$ , which corresponds to each combination of $S$ and $a$ used in (a)–(c). Shown below each snapshot is the initial condition $h(z,0)$ ; its scale is shown to the right in (d). Note the different time and amplitude scales of (d) with respect to the other panels.

In figure 6(a–c), the free surface undergoes a period of instability growth governed by the linear analysis discussed in § 3. In figure 6(a), the free surface then exhibits spatio-temporal chaotic behaviour, reminiscent of solutions to the Kuramoto–Sivashinsky equation. Indeed, all the models considered here reduce to some form of the KS equation in the limit of weakly nonlinear disturbances. In figures 6(b) and 6(c), instabilities grow into the strongly nonlinear regime where a nonlinear balance of forces allows the formation of more coherent wave motion. The waves propagate up the tube while undergoing nonlinear interactions with one another, eventually settling into a steady wavetrain. The effect of dispersive terms in KS-like equations has been shown to be regularizing, leading to the formation of travelling waves, see, e.g., Chang, Demekhin & Kopelevich (Reference Chang, Demekhin and Kopelevich1993), Ogawa (Reference Ogawa1993) and Johnson et al. (Reference Johnson2012). Because of this, it is perhaps to be expected that solutions to the fully nonlinear models (LW-P) and (TF-P) saturate to a coherent travelling wave regime.

Figure 6(d) shows the evolution of a perturbed interface for thin-film model (TF-P), for the parameter choice $\tilde{S}=0.3$ . There is a period of wave collisions and mergers following the initial instability growth, but unlike the evolution depicted for the long-wave model, the largest wave consumes all others it overtakes. Kerchman (Reference Kerchman1995) found this to be the case in all his simulations with $0.15<\tilde{S}<1.0$ (see figures 4 and 5 in Kerchman Reference Kerchman1995 for results with different values of $\tilde{S}$ in this range). Thus, even for relatively thin films, i.e.  $a=1.05$ , the behaviour of the thin-film model solutions is qualitatively different from its counterpart in the corresponding long-wave model. This difference in solution behaviour can be attributed to the thin-film model’s lack of memory of the tube’s centreline. As the annular fluid waves grow in the long-wave model, impinging on the core region, the core fluid moves faster to maintain a constant volume flux. The additional stress caused by the increased velocity of the core fluid contributes to the nonlinear balance of travelling wave solutions through the hyperbolic term $f_{1}(R;a)R_{z}$ . It should be noted that for $R\ll 1$ this term behaves like $R^{-6}R_{z}$ . In the thin-film model there is no memory of the centreline of the tube, and instabilities can grow unchecked until they exhaust the supply of fluid in the rest of the domain. Only then, due to conservation of $h$ , do the instabilities saturate.

In figure 7, the initial period of instability growth shown in figure 6(b) is compared with the linear theory of (3.2). A single wavelength with wavenumber closest to $k_{m}=1/\sqrt{2}$ while still fitting in the periodic domain is shown at equally spaced time intervals in figure 7(a,b). For times $t<25$ the growth is in excellent agreement with the linear theory. Figure 7(c) shows the amplitude (trough to crest) of the wave during this initial period of growth.

Figure 7. Snapshots of a single wave growing according to the linear theory of (3.2) (dotted line) and the nonlinear model (LW-P) (solid line) at (a)  $t=0$ , 8.3, 16.7 and (b)  $t=25$ , 33.3, 41.7 with $S=11.25$ and $a=1.1$ . (c) The amplitude (distance between wave crest and trough) of the linear theory and strongly nonlinear model. Circles correspond to the model solutions shown in (a) and (b).

4.2. Travelling wave solutions

The numerical simulations show evidence of the dominant role played by travelling waves in most time evolutions governed by some of the models. As they are easily computable numerically with any desired accuracy, travelling wave solutions can thus be used to further expose the different mathematical and physical features exhibited by the models described in § 4.1, and to explore the kinematics of the underlying fluid, which is relevant for studying transport phenomena in these flows.

Figure 8. (a) Maximum value of $h$ attained in travelling wave solutions of (LW-P) and (TF-P) for various values of $\tilde{S}$ and  $a$ . (b) Solutions to the long-wave model (LW-P) for $\tilde{S}=0.5$ are plotted (solid lines) for various values of  $a$ . The solution to the thin-film model (TF-P) for $\tilde{S}=0.5$ is plotted (dashed line). The crosses in (a) correspond to profiles shown in (b).

We once again focus here on the case of pressure-driven flow without gravity. Moving to a frame of reference travelling with the wave, i.e. letting

and substituting (4.1) into the long-wave (LW-P) and thin-film (TF-P) models yields fourth-order ordinary differential equations (ODEs)

respectively. As neither (4.2) or (4.3) is solvable analytically, we seek solutions numerically as in Camassa et al. (Reference Camassa, Forest, Lee, Ogrosky and Olander2012). We do so by providing an initial guess, via a long-time simulation of the model PDEs, to an iteration scheme to solve a boundary value problem. Specifically, we let the time-dependent simulations settle to a steady series of travelling waves, from which a single wave form can be ‘clipped’. This wave is then used as an initial approximation, within an appropriate domain, to be refined with a collocation method for a two-point boundary value problem solver.

Figure 8(a) shows the amplitude (maximum value of  $h$ ) of travelling wave solutions to the long-wave model (LW-P) and thin-film model (TF-P) for two values of $\tilde{S}$ and various values of  $a$ . The long-wave model predicts wave amplitudes smaller than those of the thin-film model; even for moderate values of $a$ this difference is significant. Figure 8(b) shows several of the travelling wave solutions for $\tilde{S}=0.5$ . As expected, in the thin-film limit $a\rightarrow 1$ the long-wave model solutions approach the waveform of their thin-film model counterparts.

Next, we investigate the effect of period cell length on solutions for each model class. Figure 9 depicts travelling waves for five different period cell lengths for (a) the long-wave model (LW-P) and (b) the thin-film model (TF-P). As the period cell length increases in the long-wave model, the geometry of the tube restricts growth of the wave amplitude. In model (TF-P) the wave amplitude grows roughly proportionally to the cell length, as shown by Kerchman (Reference Kerchman1995) (see figure 10 in Kerchman Reference Kerchman1995 for similar growth with a different value of  $\tilde{S}$ ). This is consistent with the model solution behaviour described in § 4.1: in the thin-film model (TF-P), instability growth is bounded by mass conservation, while growth in model (LW-P) is only bounded primarily by mass conservation for very thin films, with thicker film growth being bounded primarily by the core fluid’s increased velocity. This difference has implications for the thickness of the substrate layer outside the wave crest as well; in the long-wave model, as the cell length increases and the geometry restricts the growth of the wave, the substrate layer increases in thickness to accommodate the additional liquid that the increasing cell length introduces. This effect is missed by the thin-film model where the waves are allowed to grow unbounded in amplitude; in fact, the substrate layer slightly decreases in thickness with increasing cell length.

Figure 9. Various wavelength travelling wave solutions to (a) the long-wave model (LW-P) with $S=4.17$ and $a=1.25$ , and (b) the thin-film model (TF-P) with $\tilde{S}=1.0$ . The mean film thickness has been normalized to $h_{0}=1$ .

It is instructive to view these differences between travelling wave solutions of the two model classes by monitoring the nonlinear balance of the terms in (4.2) and (4.3), reported here for convenience after taking the overall derivatives as

respectively, where

and

In each of these ODEs, the first two terms give the hyperbolic contributions, the third and fourth terms give the weakly nonlinear contributions due to capillary forces and the final two terms give the strongly nonlinear capillary contributions. Figure 10(a) shows the profile of a solution to both equations (4.4) for $\tilde{S}=0.5$ . The contribution of each of the six terms $\tilde{g}_{i}(Z)$ , $1\leqslant i\leqslant 6$ , to the solution of the thin-film equation in (4.4) is shown in figure 10(b) and the contribution of each of the six terms $g_{i}(Z)$ to the solution of the long-wave equation in (4.4) is shown in figure 10(c).

Figure 10. (a) Travelling wave solution to the long-wave model (LW-P) and the thin-film model (TF-P) for $\tilde{S}=0.5$ . (b) Contributions of the six terms $g_{i}(Z)$ for the thin-film model solution shown in (a). The inset shows the percentage contribution $e_{i}$ defined in (4.7). (c) Contributions of the six terms $g_{i}(Z)$ for the long-wave model solution shown in (a).

The relative contribution of each term $\tilde{g}_{i}(Z)$ in the thin-film model (4.4) can be assessed by evaluating their mean contributions

where $L$ is the solution domain length; this gives a percentage measure of the overall contribution of each term to the solution across the domain. These percentages are shown in the inset of figure 10(b). The analogous percentages for the corresponding long-wave terms in (4.4) are shown in the inset of figure 10(c). We remark that the term $g_{1}$ which represents the free surface stress due to the core flow is a much larger contributor to the nonlinear balance in the long-wave model than its counterpart $\tilde{g}_{1}$ in the thin-film model. This suggests that, as in the case of periodic solutions to the evolution equation found in the previous section, the free surface stress due to the core flow is underestimated by the thin-film model; this agrees with the larger amplitudes that this model’s travelling waves can attain, as free surface stress can act as a limiting factor to amplitude growth.

Figure 11. (a) Travelling wave solution to the long-wave model (LW-P) for $\tilde{S}=0.5$ and the thin-film model (TF-P) for $\tilde{S}=0.2$ . (b) Contributions of the six terms $g_{i}(Z)$ for the thin-film model solution shown in (a). The inset shows the percentage contribution $e_{i}$ defined in (4.7). (c) Contributions of the six terms $g_{i}(Z)$ for the long-wave model solution shown in (a).

Finally, we note that by fitting the parameter $\tilde{S}$ in the thin-film model its solutions can be brought to provide better approximations to their long-wave model counterparts. Figure 11 compares the same long-wave solution for $\tilde{S}=0.5$ as shown in figure 10 with the thin-film model solution for $\tilde{S}=0.2$ . While still significantly larger, the amplitude of the thin-film solution is now closer to that of the long-wave model. Moreover, the relative contributions of the various thin-film terms shown in figure 10(b) closely mirror the contributions of the long-wave terms shown in (c). This suggests that an alternative definition of the function $\tilde{S}(S)$ could improve the thin-film model’s ability to represent the dynamical effects due to cylindrical geometry for a wider range of film thicknesses.

4.3. Streamline patterns of travelling waves and mass transport

To examine the flow pattern of the fluid comprising the annular film in the wave frame, it useful to look at the streamfunction defined by the annular fluid’s radial and axial velocity components $u$ and $w$ respectively. The streamfunction is

for the long-wave and thin-film models respectively; here, $y=(a-r)/(a-1)$ is a stretched radial coordinate for the thin-film models. A straightforward calculation shows that

respectively.

Streamline patterns for three travelling wave solutions to model (LW-P) are depicted in figure 12, where the fluid thickness parameter $a=1.25$ for all three solutions. For small values of $S$ as in figure 12(a), the streamlines are fanned and constricted by the interfacial disturbance; such waves will be referred to as ‘shear-enhancing’ waves. For larger values of $S$ as in figure 12(b,c), the additional feature of a region of fluid recirculation, or trapped core, near the wave crest is seen by the presence of closed streamlines; we will refer to these as ‘mass-transport’ waves. The volume of the recirculation region for travelling wave solutions is shown as a function of $S$ in figure 12(d). (We note that the absence of a recirculation region in figure 12(a) is not an artefact from lack of resolution, a point considered in more detail in Camassa et al. Reference Camassa, Forest, Lee, Ogrosky and Olander2012.) The wave solutions thus undergo a transition from shear-enhancing waves to mass-transport waves as $S$ increases. For some critical value of $S$ , defined as the largest value for which no recirculation region exists (in figure 12 this is $S_{c}\approx 0.6$ ), there is exactly one stationary point (in the travelling reference frame) at the wave crest. As $S$ increases past $S_{c}$ , this one fixed point produces three fixed points: two hyperbolic points at the free surface, which are endpoints of a separatrix dividing the closed streamlines from the open ones, and one elliptic fixed point at the centre of the closed streamlines. This change of streamline pattern as the parameters vary can be viewed as a bifurcation of the characteristic equations (4.8) and (4.9) governing the fluid parcel motion.

Figure 12. (a–c) Streamlines plotted in a travelling reference frame for solutions to the long-wave model (LW-P) for varying values of $S$ (0.4 (a), 0.7 (b), 1.1 (c)) and $a=1.25$ . (d) Volume of the recirculation region as a function of  $S$ . For $a=1.25$ , $S<0.6$ , no recirculation region exists. Each cross corresponds to one of the three solutions shown.

Streamline patterns for $S=1$ and varying film thickness parameter $a$ are shown in figure 13(a–e). As the film thins, the waves again transition from shear-enhancing waves to mass-transport waves. Figure 13(f) shows the dependence of the wave volume on the film thickness. It is interesting to note that there is some film thickness for which the wave volume is maximized (in figure 13 this is $a\approx 1.1$ ). The volume of a recirculation region in films thinner than this critical thickness is limited by the amount of liquid in the domain. For thicker films, wave amplitudes are suppressed by the core flow, as discussed in § 4.1; increasing film thickness thus restricts, and eventually eliminates altogether, the size of the recirculation region.

Figure 13. (a–e) The same as figure 12(a–c) but for fixed $S=1$ and varying $a$ (1.35 (a), 1.32 (b), 1.25 (c), 1.14 (d), 1.04 (e)). (f) The volume of the recirculation region as a function of  $a$ . For thick films ( $a>1.32$ when $S=1$ ) no recirculation region exists.

Figure 14 shows the streamlines for various values of the surface tension parameter $\tilde{S}$ for the thin-film model (TF-P). Unlike solutions to (LW-P), there is no transition from shear-enhancing waves to mass-transport waves; all waves are found to be mass-transport waves. Three fixed points exist in all wave solutions: two on the free surface and one in the recirculation zone, as described above. The volume of the recirculation region vanishes as $\tilde{S}\rightarrow 0$ . This difference in streamline topology is summarized in figure 15, where the regions of mass-transport versus shear-enhancing waves for (LW-P) are shown in $S$ – $a$ parameter space. The vertical boundary on the left represents the thin-film limit $a\rightarrow 1$ .

Figure 14. (a–c) The same as figure 12(a–c) but for (TF-P) ( $S=0.4$ (a), 0.7 (b), 1.1 (c)). (d) The volume of the recirculation region as a function of  $\tilde{S}$ ; every travelling wave has a non-zero volume of fluid inside a recirculation region.

Figure 15. Regions of $S$ – $a$ parameter space where mass-transport waves or shear-enhancing waves exist for the long-wave model (LW-P).

The degree to which the recirculation region contributes to the overall annular flux in travelling wave solutions to (LW-P) can be quantified by separating the total annular flux into the portion due to the recirculation region and the portion due to the underlying substrate, i.e.

respectively, where $R^{\ast }(Z)$ denotes the radial location of the separatrix at locations $Z$ that lie within the domain of the recirculation region, and $R^{\ast }(Z)=R(Z)$ for all other  $Z$ . Each quantity in (4.12) can be averaged over one wavelength, i.e.  $\langle q_{r}^{(a)}\rangle$ and $\langle q_{s}^{(a)}\rangle$ , where $\langle \cdot \rangle =L^{-1}\int _{0}^{L}\cdot \,\text{d}Z$ . Figure 16(a) shows each of these quantities, as well as the average total flux $\langle q^{(a)}\rangle$ , for solutions with $a=1.25$ and varying $S$ ; it should be noted that in (a) the dimensionless annular volume flux has been rescaled by the viscosity ratio $m={\it\mu}^{(a)}/{\it\mu}^{(c)}$ as liquid scales are being emphasized here. The percentage that each portion of the flux contributes is shown in figure 16(b). Consistent with figure 12, no recirculation region exists for $S\leqslant S_{c}\approx 0.6$ . For larger $S$ values, as $S$ increases, the recirculation region grows in size and contributes an increasingly large fraction of the total flux. We note that the percentage of flux attributable to the recirculation region appears to converge to 50 % in the limit of large $S$ in (b), though this limiting value cannot be expected to hold for all values of the film thickness parameter  $a$ .

Figure 16. (a) The total dimensionless flux $m\langle q^{(a)}\rangle$ for travelling wave solutions to (LW-P) with $a=1.25$ . The volume fluxes due to the recirculation region $m\langle q_{r}^{(a)}\rangle$ and the substrate $m\langle q_{s}^{(a)}\rangle$ are shown as well. (b) Percentage of the total flux due to the recirculation region and substrate in (a).

The source of the difference in behaviour between the two models is worth discussing. We remark that a recirculation region can only exist when the speed of a fluid parcel somewhere on the free surface is identical to the wave speed; such a parcel would be located at one of the two fixed stagnation points lying on the free surface in the travelling reference frame. The velocity of a fluid parcel lying on an unperturbed interface in the long-wave (LW-P) and thin-film (TF-P) models is

respectively. The phase speed of a linear disturbance is the real part of ${\it\omega}/k$ in (3.2) and (3.3),

Thus, in model (TF-P) the fluid can match the wave speed, $\tilde{w}=\tilde{c}$ , and it can be expected that there exist fluid parcels moving at the speed of the wave regardless of how small the free surface disturbance is. The long-wave model (LW-P), however, has a larger disturbance speed than interface speed, i.e.  $w<c$ , guaranteeing that no recirculation region can exist for small-amplitude waves.

We stress that simply retaining one additional power of ${\it\beta}$ in the hyperbolic terms is sufficient to prevent this phenomenon of persistence of closed streamlines at all travelling wave amplitudes from occurring in the thin-film model. If $F=0$ in (TF-PG), then the phase speed of a linear disturbance in (TF-PG) is now given by $\tilde{c}_{2}=2{\it\beta}+5{\it\beta}^{2}=c+O({\it\beta}^{3})$ , while the speed of the undisturbed interface is $\tilde{w}_{2}=w$ ; thus the higher-order thin-film model (TF-PG) contains the same streamline bifurcation as seen in its long-wave counterpart.