3.1 Problem set-up

Consider the turbulent flow of air through a wavy channel with mean channel half-width $\bar{L}_{0}$ and local half-width $\bar{L}(\bar{z})=\bar{L}_{0}+\bar{\unicode[STIX]{x1D6FC}}(\bar{z})$ . The time-averaged flow is governed by the incompressible Reynolds-averaged Navier–Stokes (RANS) equations

(3.1a ) $$\begin{eqnarray}\displaystyle \bar{\unicode[STIX]{x1D70C}}^{(g)}(\bar{U}\bar{U}_{\bar{x}}+\bar{W}\bar{U}_{\bar{z}}) & = & \displaystyle -\bar{P}_{\bar{x}}+\bar{\unicode[STIX]{x1D707}}^{(g)}\bar{U}_{\bar{x}\bar{x}}+\bar{\unicode[STIX]{x1D707}}^{(g)}\bar{U}_{\bar{z}\bar{z}}+\bar{\unicode[STIX]{x1D70F}}_{11\bar{x}}+\bar{\unicode[STIX]{x1D70F}}_{12\bar{z}},\end{eqnarray}$$

(3.1b ) $$\begin{eqnarray}\displaystyle \bar{\unicode[STIX]{x1D70C}}^{(g)}(\bar{U}\bar{W}_{\bar{x}}+\bar{W}\bar{W}_{\bar{z}}) & = & \displaystyle -\bar{P}_{\bar{z}}+\bar{\unicode[STIX]{x1D707}}^{(g)}\bar{W}_{\bar{x}\bar{x}}+\bar{\unicode[STIX]{x1D707}}^{(g)}\bar{W}_{\bar{z}\bar{z}}+\bar{\unicode[STIX]{x1D70F}}_{12\bar{x}}+\bar{\unicode[STIX]{x1D70F}}_{22\bar{z}},\end{eqnarray}$$

(3.1c ) $$\begin{eqnarray}\displaystyle \bar{U}_{\bar{x}}+\bar{W}_{\bar{z}} & = & \displaystyle 0,\end{eqnarray}$$

$(\bar{U},\bar{W})$

$(\bar{x},\bar{z})$

$\bar{P}$

$\bar{\unicode[STIX]{x1D70F}}_{ij}$

$\bar{\unicode[STIX]{x1D70C}}^{(g)}$

$\bar{\unicode[STIX]{x1D707}}^{(g)}$

$\bar{z}$

$\bar{x}$

$\bar{x}=\bar{L}(\bar{z})$

$\bar{x}=0$

$\bar{x}=\bar{L}(\bar{z})$

(3.2) $$\begin{eqnarray}\displaystyle \bar{U}=\bar{W}=0, & & \displaystyle\end{eqnarray}$$

and at the centre of the channel $\bar{x}=0$ we require that

(3.3) $$\begin{eqnarray}\displaystyle \bar{W}_{\bar{x}}=0. & & \displaystyle\end{eqnarray}$$

Figure 2.

Variables definition sketch for wavy channel flow.

Next, the problem can be non-dimensionalized by the following reference scales:

(3.4) $$\begin{eqnarray}\displaystyle \left.\begin{array}{@{}c@{}}\displaystyle x=\bar{x}/\bar{L}_{C},\quad z=\unicode[STIX]{x1D716}_{C}\bar{z}/\bar{L}_{C},\quad U=\bar{U}/\unicode[STIX]{x1D716}_{C}\bar{W}_{C},\\ \quad W=\bar{W}/\bar{W}_{C},\quad P=\unicode[STIX]{x1D716}_{C}\bar{P}/|\bar{T}_{C}|,\quad \unicode[STIX]{x1D70F}=\bar{\unicode[STIX]{x1D70F}}/|\bar{T}_{C}|,\end{array}\right\} & & \displaystyle\end{eqnarray}$$

where $\unicode[STIX]{x1D716}_{C}=\bar{L}_{0}/\bar{\unicode[STIX]{x1D706}}$ is an aspect ratio, the velocity scale

(3.5) $$\begin{eqnarray}\displaystyle \bar{W}_{C}=\frac{\bar{Q}_{C}}{2\bar{L}_{0}}\sqrt{\frac{F_{C}}{2}}, & & \displaystyle\end{eqnarray}$$

is the friction velocity estimated using the empirical skin friction coefficient

(3.6a,b ) $$\begin{eqnarray}\displaystyle F_{C}=0.073(Re_{C}^{(g)})^{-0.25},\quad Re_{C}^{(g)}=\frac{\bar{\unicode[STIX]{x1D70C}}^{(g)}\bar{Q}_{C}}{2\bar{\unicode[STIX]{x1D707}}^{(g)}}, & & \displaystyle\end{eqnarray}$$

and where $Re_{C}^{(g)}$ is the bulk Reynolds number and $\bar{Q}_{C}$ is the prescribed area flux of air. Subscripts of $C$ are used to denote the case of channel flow in order to avoid confusion with the analogous quantities defined for pipe flow in § 4. Note that the bulk Reynolds number and skin friction coefficient do not depend on the average channel half-width $\bar{L}_{0}$ , but only on the fluid properties and flux $\bar{Q}_{C}$ . The length and stress/pressure scales are

(3.7a,b ) $$\begin{eqnarray}\displaystyle \bar{L}_{C}=\frac{\bar{\unicode[STIX]{x1D707}}^{(g)}}{\bar{\unicode[STIX]{x1D70C}}^{(g)}\bar{W}_{C}},\quad |\bar{T}_{C}|=\bar{\unicode[STIX]{x1D70C}}^{(g)}(\bar{W}_{C})^{2}=\frac{\bar{\unicode[STIX]{x1D70C}}^{(g)}F_{C}\bar{Q}_{C}^{2}}{8\bar{L}_{0}^{2}}. & & \displaystyle\end{eqnarray}$$

Here $|\bar{T}_{C}|$ is an empirical estimate of the wall stress in flow through a flat channel of half-width $\bar{L}_{0}$ . Thus, given the ‘hardware’ parameters $\bar{L}_{0}$ , $\bar{\unicode[STIX]{x1D706}}$ , $\bar{Q}_{C}$ , $\bar{\unicode[STIX]{x1D707}}^{(g)}$ and $\bar{\unicode[STIX]{x1D70C}}^{(g)}$ , i.e. parameters that are a fixed property of the fluid or geometry (as opposed to parameters and scalings selected by initial conditions), all scales of the problem are determined. We note that recent studies have shown that liquid-infused porous surfaces are able to reduce the skin friction; drag reduction has been demonstrated for both turbulent and laminar flow (see, e.g. Solomon, Khalil & Varanasi (Reference Solomon, Khalil and Varanasi2014) and Rosenberg et al. (Reference Rosenberg, Van Buren, Fu and Smits2016), respectively).

Using this scaling, the dimensionless wall is given by $x=L(z)=L_{0}+\unicode[STIX]{x1D6FC}(z)$ , where $L_{0}=\bar{L}_{0}/\bar{L}_{C}$ is the dimensionless channel width and $\unicode[STIX]{x1D6FC}=\bar{\unicode[STIX]{x1D6FC}}/\bar{L}_{C}$ measures the local constriction of the channel; $\unicode[STIX]{x1D6FC}(z)<0$ indicates a narrow region of the channel while $\unicode[STIX]{x1D6FC}(z)>0$ indicates a relatively wide region of the channel. Substituting (3.4) into (3.1) results in

(3.8a ) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D716}_{C}^{3}(UU_{x}+WU_{z}) & = & \displaystyle -P_{x}+\unicode[STIX]{x1D716}_{C}^{2}U_{xx}+\unicode[STIX]{x1D716}_{C}^{4}U_{zz}+\unicode[STIX]{x1D716}_{C}\unicode[STIX]{x1D70F}_{11x}+\unicode[STIX]{x1D716}_{C}^{2}\unicode[STIX]{x1D70F}_{12z},\end{eqnarray}$$

(3.8b ) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D716}_{C}(UW_{x}+WW_{z}) & = & \displaystyle -P_{z}+W_{xx}+\unicode[STIX]{x1D716}_{C}^{2}W_{zz}+\unicode[STIX]{x1D70F}_{12x}+\unicode[STIX]{x1D716}_{C}\unicode[STIX]{x1D70F}_{22z},\end{eqnarray}$$

(3.8c ) $$\begin{eqnarray}\displaystyle U_{x}+W_{z} & = & \displaystyle 0.\end{eqnarray}$$

3.2 Coordinate transformation

Analysis in the usual Cartesian coordinates is difficult and, as previous studies have demonstrated, significant simplifications and improved agreement with experiments can arise if a curvilinear coordinate system is used instead. We move to a coordinate system which ‘flattens’ the wavy channel; the base flow in this flattened channel is estimated by using an effective viscosity approximated by a simple mixing-length model. This approach is similar to that used by, e.g. Benjamin (Reference Benjamin1959) and Tseluiko & Kalliadasis (Reference Tseluiko and Kalliadasis2011), but is not equivalent: here, a ‘long-wave’ approximation will be made from the outset, i.e.

(3.9) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D716}_{C}\ll 1. & & \displaystyle\end{eqnarray}$$

This results in a base flow slowly modulated in the $z$ direction to take into account the channel’s local half-width. This slowly modulated base flow is analogous to the base flow used in the locally Poiseuille modelling approach of Camassa et al. (Reference Camassa, Forest, Lee, Ogrosky and Olander2012) and is used here to estimate first-order corrections to the wall stress.

Figure 3.

Gridlines in the curvilinear coordinate system $(\tilde{x},\tilde{z})$ for flow through a wavy channel. Local basis vectors are depicted for Cartesian coordinates (black) and curvilinear coordinates (red) at one point. Points $r_{1}$ and $r_{2}$ are discussed in the text.

We undertake a coordinate transformation $(x,z)\rightarrow (\tilde{x},\tilde{z})$ , where we define

(3.10) $$\begin{eqnarray}\displaystyle \tilde{x}=\frac{L_{0}}{L(z)}x, & & \displaystyle\end{eqnarray}$$

so that $\tilde{x}=L_{0}$ is coincident with the wavy wall $x=L(z)$ , $\tilde{x}=0$ is coincident with the centre of the tube $x=0$ , and curves of constant $\tilde{x}$ are relaxed linearly to the Cartesian coordinate system as one moves away from the wall to the centre of the channel. The curves of constant $\tilde{z}$ are then defined so that they are orthogonal to curves of constant $\tilde{x}$ ; see figure 3. The linear relaxation of the curves of constant $\tilde{x}$ lies somewhat between the exponential relaxation of these curves used by Benjamin (Reference Benjamin1959) and Tseluiko & Kalliadasis (Reference Tseluiko and Kalliadasis2011), and the lack of any relaxation used by Thorsness et al. (Reference Thorsness, Morrisroe and Hanratty1978). Tildes will denote the new coordinate system throughout this section; subscripts of $1$ or $2$ denote components in the $\tilde{x}$ or $\tilde{z}$ direction, respectively (or in the $x$ or $z$ direction, respectively, depending on the context).

Given the long-wave approximation $\unicode[STIX]{x1D716}_{C}\ll 1$ , the partial derivatives with respect to $x$ and $z$ and the Laplacian $\unicode[STIX]{x0394}\equiv \unicode[STIX]{x1D6FB}^{2}$ are expressed in the new coordinate system by (A 1) and (A 2) in appendix A, respectively; likewise, the velocities $(U,W)$ and Reynolds stresses $\unicode[STIX]{x1D70F}_{ij}$ can be transformed into $(\tilde{U} ,\tilde{W})$ and $\tilde{\unicode[STIX]{x1D70F}}_{ij}$ , respectively, according to (A 3). Each of these flow variables in the new coordinate system can then be expanded in the small parameter $\unicode[STIX]{x1D716}_{C}$ ,

(3.11a ) $$\begin{eqnarray}\displaystyle \tilde{U} (x,z) & = & \displaystyle \tilde{U} _{0}(\tilde{x},\tilde{z})+\unicode[STIX]{x1D716}_{C}\tilde{U} _{1}(\tilde{x},\tilde{z})+O(\unicode[STIX]{x1D716}_{C}^{2}),\end{eqnarray}$$

(3.11b ) $$\begin{eqnarray}\displaystyle \tilde{W}(x,z) & = & \displaystyle \tilde{W}_{0}(\tilde{x},\tilde{z})+\unicode[STIX]{x1D716}_{C}\tilde{W}_{1}(\tilde{x},\tilde{z})+O(\unicode[STIX]{x1D716}_{C}^{2}),\end{eqnarray}$$

(3.11c ) $$\begin{eqnarray}\displaystyle \tilde{P}(x,z) & = & \displaystyle \tilde{P}_{0}(\tilde{x},\tilde{z})+\unicode[STIX]{x1D716}_{C}\tilde{P}_{1}(\tilde{x},\tilde{z})+O(\unicode[STIX]{x1D716}_{C}^{2}),\end{eqnarray}$$

(3.11d ) $$\begin{eqnarray}\displaystyle \tilde{\unicode[STIX]{x1D70F}}_{ij}(x,z) & = & \displaystyle \tilde{\unicode[STIX]{x1D70F}}_{ij_{0}}(\tilde{x},\tilde{z})+\unicode[STIX]{x1D716}_{C}\tilde{\unicode[STIX]{x1D70F}}_{ij_{1}}(\tilde{x},\tilde{z})+O(\unicode[STIX]{x1D716}_{C}^{2}).\end{eqnarray}$$

The essence of the locally Poiseuille approach in Camassa et al. (Reference Camassa, Forest, Lee, Ogrosky and Olander2012) is to seek velocity solutions whose leading-order components satisfy

(3.12a,b ) $$\begin{eqnarray}\displaystyle \tilde{U} _{0}(\tilde{x},\tilde{z})=0,\quad \tilde{W}_{0}(\tilde{x},\tilde{z}) & = & \displaystyle \tilde{w}_{0}(\tilde{x})g_{0}(\tilde{z}),\end{eqnarray}$$

so that the leading-order velocity profile in the new coordinate system is a parallel flow whose speed is modulated by the wall waviness present in the Cartesian coordinate system. We now make one further assumption, that the normal Reynolds stresses are $O(\unicode[STIX]{x1D716}_{C})$ , i.e. $\tilde{\unicode[STIX]{x1D70F}}_{ii_{0}}=0$ for $i=1,2$ , and adopt a turbulent viscosity model for the Reynolds shear stress $\tilde{\unicode[STIX]{x1D70F}}_{12}$ , i.e.

(3.13) $$\begin{eqnarray}\displaystyle \tilde{\unicode[STIX]{x1D70F}}_{12_{0}}(\tilde{x},\tilde{z})=\tilde{\unicode[STIX]{x1D707}}^{\ast }(\tilde{x},\tilde{z})\tilde{w}_{0}^{\prime }\tilde{g}_{0}\frac{L_{0}}{L}, & & \displaystyle\end{eqnarray}$$

where primes denote $\text{d}/\text{d}\tilde{x}$ . It is assumed that the turbulent viscosity in this new coordinate system can be approximated by the mixing-length model,

(3.14) $$\begin{eqnarray}\displaystyle \tilde{\unicode[STIX]{x1D707}}^{\ast }(\tilde{x},\tilde{z}) & = & \displaystyle \frac{L}{L_{0}}[l_{C}(\tilde{x})]^{2}|\tilde{w}_{0}^{\prime }|\tilde{g}_{0},\end{eqnarray}$$

where $\tilde{g}_{0}$ is the $\tilde{z}$ -component of the leading-order streamwise velocity. The mixing length $l_{C}$ can be modelled by the van Driest equation,

(3.15) $$\begin{eqnarray}\displaystyle l_{C}(\tilde{x})=\unicode[STIX]{x1D705}(L_{0}-\tilde{x})[1-\exp ((\tilde{x}-L_{0})/A)], & & \displaystyle\end{eqnarray}$$

where $\unicode[STIX]{x1D705}\approx 0.41$ is von Kármán’s constant, and $A=26$ is a constant determined empirically (van Driest Reference van Driest1956; Bailey et al. Reference Bailey, Vallikivi, Hultmark and Smits2014; Schlichting & Gersten Reference Schlichting and Gersten2017). While no mixing length model can adequately capture all flow features near the wall, and while various modifications have been proposed to this simple model, we elect to use the original formulation of van Driest due to its simplicity and ability to produce the main features of experimental flow near the wall over a wide range of Reynolds numbers (van Driest Reference van Driest1956).

3.3 Leading-order solution

Substitution of (A 1)–(A 3) and (3.11)–(3.15) into (3.8) and collection of the leading-order terms in $\unicode[STIX]{x1D716}_{C}$ results in

(3.16a ) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x2202}_{\tilde{x}}\tilde{P}_{0} & = & \displaystyle 0,\end{eqnarray}$$

(3.16b ) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x2202}_{\tilde{z}}\tilde{P}_{0} & = & \displaystyle \frac{\tilde{g}_{0}L_{0}^{2}}{L^{2}}\tilde{w}_{0}^{\prime \prime }+\frac{\tilde{g}_{0}^{2}L_{0}}{L}\frac{\text{d}}{\text{d}\tilde{x}}\left[l_{C}^{2}(\tilde{w}_{0}^{\prime })^{2}\right],\end{eqnarray}$$

(3.16c ) $$\begin{eqnarray}\displaystyle \left(\frac{\tilde{g}_{0}}{L}\frac{\text{d}L}{\text{d}\tilde{z}}+\frac{\text{d}\tilde{g}_{0}}{\text{d}\tilde{z}}\right)\tilde{w}_{0} & = & \displaystyle 0,\end{eqnarray}$$

$\unicode[STIX]{x2202}_{\tilde{x}}\tilde{P}_{0}$

(3.17a ) $$\begin{eqnarray}\displaystyle \tilde{w}_{0}=0 & & \displaystyle \text{at }\tilde{x}=L_{0},\end{eqnarray}$$

(3.17b ) $$\begin{eqnarray}\displaystyle \tilde{w}_{0}^{\prime }=0 & & \displaystyle \text{at }\tilde{x}=0.\end{eqnarray}$$

Equation (3.16) is satisfied by

(3.18a,b ) $$\begin{eqnarray}\displaystyle \tilde{g}_{0}=\frac{L_{0}}{L},\quad \unicode[STIX]{x2202}_{\tilde{z}}\tilde{P}_{0}=-\frac{L_{0}^{2}}{L^{3}}. & & \displaystyle\end{eqnarray}$$

The $\tilde{x}$ component of the leading-order streamwise velocity, i.e. $\tilde{w}_{0}$ , can then be found by solving

(3.19) $$\begin{eqnarray}\displaystyle -\frac{1}{L_{0}} & = & \displaystyle \tilde{w}_{0}^{\prime \prime }+\frac{\text{d}}{\text{d}\tilde{x}}[l_{C}^{2}(\tilde{w}_{0}^{\prime })^{2}],\end{eqnarray}$$

which, after integrating once, is a straightforward quadratic equation in $\tilde{w}_{0}^{\prime }$ which can be solved and integrated numerically using boundary conditions (3.17) to give the leading-order velocity profile. Figure 4 shows the velocity profile for two different Reynolds numbers. We remark again that in flows through a wavy channel the Reynolds number, and hence the skin friction coefficient, does not vary as the local channel half-width varies. This velocity profile thus forms a set of self-similar solutions that create wall stresses identical to those of the laminar solutions used with an effective viscosity in Camassa et al. (Reference Camassa, Forest, Lee, Ogrosky and Olander2012). This leading-order wall stress is given by

(3.20) $$\begin{eqnarray}\displaystyle \bar{\unicode[STIX]{x1D70F}}_{C_{0}}^{(g)} & = & \displaystyle -|\bar{T}_{C}|\frac{L_{0}^{2}}{L^{2}},\end{eqnarray}$$

and will be referred to as the locally Poiseuille (LP) model for the wall stress.

Figure 4.

Leading-order streamwise velocity profile $\tilde{w}_{0}$ in (3.19) for channel airflow with $Re_{C}^{(g)}=7850$ and $Re_{C}^{(g)}=31\,000$ . $W_{b}=1/\bar{L}_{0}\int _{0}^{\bar{L}_{0}}\tilde{w}_{0}(x)\,\text{d}x$ denotes the bulk velocity.

3.4 First-order corrections

We next extend the locally Poiseuille approach to include first-order corrections. Due to the aspect ratio assumptions made, only the streamwise velocity corrections $\tilde{W}_{1}(\tilde{x},\tilde{z})$ are necessary for estimating the wall shear stress to $O(\unicode[STIX]{x1D716}_{C})$ , and so we include only those steps needed for estimating this term.

Similar to the leading-order component, we seek a streamwise velocity correction of the form

(3.21) $$\begin{eqnarray}\displaystyle \tilde{W}_{1}(\tilde{x},\tilde{z})=\tilde{w}_{1}(\tilde{x})\tilde{g}_{1}(\tilde{z}). & & \displaystyle\end{eqnarray}$$

Substitution of (3.21) into the momentum equations and collecting terms of $O(\unicode[STIX]{x1D716}_{C})$ results in

(3.22a ) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x2202}_{\tilde{x}}\tilde{P}_{1} & = & \displaystyle 0,\end{eqnarray}$$

(3.22b ) $$\begin{eqnarray}\displaystyle -\!\frac{L_{0}^{2}}{L^{3}}\frac{\text{d}L}{\text{d}\tilde{z}}\tilde{w}_{0}^{2} & = & \displaystyle -\unicode[STIX]{x2202}_{\tilde{z}}\tilde{P}_{1}+\frac{\tilde{g}_{1}L_{0}^{2}}{L^{2}}\left[\tilde{w}_{1}^{\prime \prime }+2\frac{\text{d}}{\text{d}\tilde{x}}(l^{2}\tilde{w}_{1}^{\prime }\tilde{w}_{0}^{\prime })\right].\end{eqnarray}$$

$\tilde{U} _{1}(\tilde{x},\tilde{z})$

$x$

$\unicode[STIX]{x2202}_{z}\tilde{W}_{1}$

The $\tilde{z}$ -component $\tilde{g}_{1}$ of the streamwise velocity corrections can be found by first integrating (3.22b ) across the channel half-width resulting in

(3.23) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x2202}_{\tilde{z}}\tilde{P}_{1} & = & \displaystyle \frac{L_{0}A_{C}}{L^{3}}\frac{\text{d}L}{\text{d}\tilde{z}}+\frac{\tilde{g}_{1}L_{0}B_{C}}{L^{2}},\end{eqnarray}$$

with

(3.24a,b ) $$\begin{eqnarray}\displaystyle A_{C}=\int _{0}^{L_{0}}\tilde{w}_{0}^{2}\,\text{d}\tilde{x},\quad B_{C}=\unicode[STIX]{x2202}_{\tilde{x}}\tilde{w}_{1}|_{\tilde{x}=L_{0}}. & & \displaystyle\end{eqnarray}$$

Thus $A_{C}$ is a known constant obtained by numerically integrating the inertial terms across the channel, and $B_{C}$ is the wall stress correction term for which we seek an approximation. Substituting (3.23) into (3.22b ) gives

(3.25) $$\begin{eqnarray}\displaystyle -\!\frac{L_{0}^{2}}{L^{3}}\frac{\text{d}L}{\text{d}\tilde{z}}\left(\tilde{w}_{0}^{2}+\frac{A_{C}}{L_{0}}\right) & = & \displaystyle \frac{\tilde{g}_{1}L_{0}^{2}}{L^{2}}\left(\tilde{w}_{1}^{\prime \prime }+2\frac{\text{d}}{\text{d}\tilde{x}}[l_{C}^{2}(\tilde{w}_{1}^{\prime })(\tilde{w}_{0}^{\prime })]-\frac{B_{C}}{L_{0}}\right),\end{eqnarray}$$

and so we set

(3.26) $$\begin{eqnarray}\displaystyle \tilde{g}_{1}(\tilde{z})=\frac{1}{L}\frac{\text{d}L}{\text{d}\tilde{z}}. & & \displaystyle\end{eqnarray}$$

The $\tilde{x}$ -component $\tilde{w}_{1}$ of the streamwise velocity gradient corrections can be found in terms of the unknown pressure gradient correction $\unicode[STIX]{x2202}_{\tilde{z}}\tilde{P}_{1}$ by integrating equation (3.25) once and applying the zero shear condition at the centre of the tube; $\tilde{w}_{1}$ can then be found in terms of $\unicode[STIX]{x2202}_{\tilde{z}}\tilde{P}_{1}$ by integrating the equation a second time while applying the no-slip condition at the wall. The resulting velocity profile $\tilde{w}_{1}$ may be integrated across the channel from $\tilde{x}=0$ to $\tilde{x}=L_{0}$ ; this integral must be set to zero by selection of the pressure gradient correction $\unicode[STIX]{x2202}_{\tilde{z}}\tilde{P}_{1}$ (or equivalently selection of $B_{C}$ ) since the leading-order profile accounts entirely for the prescribed constant volume flux. The resulting dimensional first-order wall stress corrections are

(3.27) $$\begin{eqnarray}\displaystyle \bar{\unicode[STIX]{x1D70F}}_{C_{1}}^{(g)} & = & \displaystyle |\bar{T}_{C}|\frac{\bar{L}_{0}^{2}}{\bar{L}^{2}}\frac{\text{d}\bar{L}}{\text{d}\bar{z}}B_{C}^{\ast },\end{eqnarray}$$

where

(3.28) $$\begin{eqnarray}\displaystyle B_{C}^{\ast }=\frac{\bar{L}_{C}}{\bar{L}_{0}}B_{C}=\sqrt{\frac{2}{0.079}}(Re_{C}^{(g)})^{7/8}B_{C}. & & \displaystyle\end{eqnarray}$$

Having found the first-order stress corrections, the total shear stress up to $O(\unicode[STIX]{x1D716}_{C})$ may be expressed in dimensional form as $\bar{\unicode[STIX]{x1D70F}}_{C}^{(g)}=\bar{\unicode[STIX]{x1D70F}}_{C_{0}}^{(g)}+\bar{\unicode[STIX]{x1D70F}}_{C_{1}}^{(g)}$ ,

(3.29) $$\begin{eqnarray}\displaystyle \bar{\unicode[STIX]{x1D70F}}_{C}^{(g)} & = & \displaystyle |\bar{T}_{C}|\left(\frac{\bar{L}_{0}}{\bar{L}}\right)^{2}\left[-1+B_{C}^{\ast }\frac{\text{d}\bar{L}}{\text{d}\bar{z}}\right].\end{eqnarray}$$

The first term in (3.29), i.e. the leading-order stress, represents the modelled stress found using the LP model in Camassa et al. (Reference Camassa, Forest, Lee, Ogrosky and Olander2012) (adapted here for channel flow). The second term represents the first-order corrections, and we will refer to the total stress (3.29) as the first-order or long-wave model. The locally Poiseuille stress is in phase with the wall while the first-order stress (3.29) depends both on the wall and its gradient. The dimensional normal stresses at the wall, up to $O(\unicode[STIX]{x1D716}_{C})$ , are given by

(3.30) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x2202}_{\bar{z}}\bar{P}_{C}^{(g)} & = & \displaystyle \frac{|\bar{T}_{C}|}{\bar{L}_{0}}\left(\frac{\bar{L}_{0}}{\bar{L}}\right)^{3}\left[-1+(B_{C}+A_{C})\frac{\bar{L}_{C}}{\bar{L}_{0}}\frac{\text{d}\bar{L}}{\text{d}\bar{z}}\right].\end{eqnarray}$$

The coefficients $A_{C}$ and $B_{C}^{\ast }$ are numerically determined functions of the gas Reynolds number $Re_{C}^{(g)}$ defined in (3.6).

Figure 5.

Streamwise velocity profile for first-order (solid lines) and leading-order (dotted lines) models along one wavelength of a wavy channel (figure not to scale).

A sketch of the streamwise velocity profile with first-order corrections for a small-amplitude wavy wall is shown in figure 5. At locations where the wall gradient vanishes, i.e. $\text{d}\bar{L}/\text{d}\bar{z}=0$ , the velocity profile is the base profile modulated by the channel width; e.g. the velocity profile at points A and C are identical up to a factor of $\bar{L}(\bar{z}_{A})/\bar{L}(\bar{z}_{C})$ . In regions where the channel is constricting, e.g. point B where $\text{d}\bar{L}/\text{d}\bar{z}<0$ , the air flows faster near the walls and slower at the centre of the channel than it would due to the base profile alone; where the channel widens, e.g. point D where $\text{d}\bar{L}/\text{d}\bar{z}>0$ , the opposite holds. Thus constriction results in larger wall stress while widening results in smaller wall stress than that due to the base flow.

In the case of larger-amplitude wavy deformations of flat walls, and in the regions where the airflow is being constricted, i.e. $\text{d}\bar{L}/\text{d}\bar{z}<0$ , the maximum streamwise velocity occurs at a point close to the wall (not shown in figure 5). This behaviour is likely unphysical and was not seen in pipe and channel flow experiments, see, e.g. Hsu & Kennedy (Reference Hsu and Kennedy1971). We note that the same behaviour occurs in the model derived by Tseluiko & Kalliadasis (Reference Tseluiko and Kalliadasis2011) which exploits an assumed small amplitude of wall deformations – we thus refer hereafter to this model as a small-amplitude model – though the minimum amplitude required to produce this behaviour is much larger than in the first-order model. Nevertheless, our focus here is on whether the first-order model accurately captures the near-wall behaviour, and in particular the wall stress; this behaviour is considered next.

3.5 Comparison with experiments and other models

We now compare the first-order wall stress estimate (3.29) with both the leading-order locally Poiseuille estimate and the small-amplitude model of Tseluiko & Kalliadasis (Reference Tseluiko and Kalliadasis2011), as well as experiments reported in Zilker et al. (Reference Zilker, Cook and Hanratty1977). Figure 6 shows results for increasing Reynolds number with fixed wall wave amplitude of 0.079 cm (crest-to-trough amplitude 0.157 cm) and wavelength 5.08 cm in a channel with half-width 2.54 cm. Both the small-amplitude model and first-order model (3.29) match the experimental trend of smaller corrections relative to the magnitude of the leading-order stress for increasing Reynolds number. For the higher Reynolds number, $Re=31\,000$ , the first-order estimate (3.29) appears to capture both the phase shift and amplitude of the waviness-induced wall stress at least as well as the small-amplitude model and much better than the locally Poiseuille model. For the lower Reynolds number, $Re=7850$ , the small-amplitude model appears to slightly underestimate the phase shift; the first-order model slightly overestimates the phase shift. Both underestimate the amplitude of the waviness-induced stress; both models give much better agreement than the locally Poiseuille model.

Figure 6.

Wall stress along one wavelength of a wavy channel predicted by the long-wave model (solid lines), locally Poiseuille model (dotted lines), small-amplitude model of Tseluiko & Kalliadasis (Reference Tseluiko and Kalliadasis2011) (dashed lines) and found in experiments of Zilker et al. (Reference Zilker, Cook and Hanratty1977) (circles) through a channel with half-width $\bar{L}_{0}=2.54~\text{cm}$ and a sinusoidally wavy wall with crest-to-trough amplitude 0.157 cm and wavelength 5.08 cm; the ratio of amplitude to wavelength has been enlarged in the wavy wall shown in the bottom panels. (a) $Re_{C}^{(g)}=7850$ , (b) $Re_{C}^{(g)}=31\,000$ ; compare to figures 8 and 9 in Tseluiko & Kalliadasis (Reference Tseluiko and Kalliadasis2011) and figure 6 in Zilker et al. (Reference Zilker, Cook and Hanratty1977).

We note that the long-wave aspect ratio assumption used in deriving the stress estimates may not be strictly valid in the experiments of Zilker et al. (Reference Zilker, Cook and Hanratty1977), for which $\unicode[STIX]{x1D716}_{C}=L_{0}/\unicode[STIX]{x1D706}\approx 0.5$ . The reasonable agreement with experiments suggests, however, that the combination of terms neglected in the modelling may have only a small net contribution to the wall shear stress.

We next focus on the long-wave model for higher-amplitude wall perturbations; figure 7 shows results for increasing Reynolds number with fixed wall wave amplitude of 0.127 cm (crest-to-trough amplitude 0.254 cm) and wavelength 5.08 cm. The first-order model again exhibits the experimental trend of higher corrections relative to the magnitude of the leading-order stress for increased wave amplitude (a trend also seen in the small-amplitude model). In addition, the first-order model exhibits a small phase shift as the Reynolds number increases (a trend not present in the small-amplitude model). For the smaller Reynolds number, the first-order model estimates the amplitude of the corrections well; for the higher Reynolds number, the first-order model appears to overestimate the amplitude while also not accurately capturing the phase shift.

The experiments in figure 7 also show that if the wall waviness has large enough amplitude, the wall stress is nonlinear even for a sinusoidal wall profile. For example, in figure 7(b), the minimum reported experimental wall stress occurs at $x/\unicode[STIX]{x1D706}\approx 0.15$ and the maximum occurs at $x/\unicode[STIX]{x1D706}\approx 0.85$ ; in contrast, a linear wall stress would have extrema located at 0.5 wavelengths apart. This effect is of course not captured in the small-amplitude model which neglects nonlinear interactions in velocity perturbations. The first-order model, through the dependence of $\unicode[STIX]{x1D70F}^{(g)}$ on $1/\bar{L}^{2}$ , displays this nonlinear behaviour, albeit very slightly, in qualitative agreement with experiments, i.e. the distance between the location of the model wall stress minimum and maximum in figure 7(b) is greater than 0.5 wavelengths. This nonlinear behaviour can be more readily seen in the results for wavy pipe flow in § 4.

It is also notable that in the experiments seen in figure 7(a) a region of negative wall stress was observed; the first-order model does not capture this behaviour in these figures, but does display this behaviour for smaller Reynolds number (not shown). It is important to note that the first-order model was derived under the assumption that the ratio of the characteristic velocity scales in the transverse and streamwise directions was small, $\bar{u}_{c}=\unicode[STIX]{x1D716}_{C}\bar{w}_{c}$ . In these regions of negative stress this assumption is clearly violated, as the negative shear stress indicates regions of air recirculation along the wavy wall troughs. However, as is frequently the case with asymptotics, the asymptotic formula (3.29) appears to have validity outside of the range of strict asymptotic derivation for small $\unicode[STIX]{x1D716}_{C}$ values.

Figure 7.

Same as figure 6 but with wall crest-to-trough amplitude 0.254 cm and wavelength 5.08 cm; the ratio of amplitude to wavelength has been enlarged in the wavy wall shown in the bottom panel (small-amplitude model not shown). (a) $Re_{C}^{(g)}=19\,260$ ,(b) $Re_{C}^{(g)}=31\,850$ . Compare to figure 7 in Zilker et al. (Reference Zilker, Cook and Hanratty1977).

4.1 Governing equations

The governing equations are the RANS in cylindrical coordinates,

(4.1a ) $$\begin{eqnarray}\displaystyle \bar{\unicode[STIX]{x1D70C}}^{(g)}(\bar{U}\bar{U}_{\bar{r}}+\bar{W}\bar{U}_{\bar{z}}) & = & \displaystyle -\bar{P}_{\bar{r}}+\bar{\unicode[STIX]{x1D707}}^{(g)}\left(\frac{1}{\bar{r}}\unicode[STIX]{x2202}_{\bar{r}}(\bar{r}\unicode[STIX]{x2202}_{\bar{r}}\bar{U})+\bar{U}_{\bar{z}\bar{z}}-\frac{\bar{U}}{\bar{r}^{2}}\right)+\frac{1}{\bar{r}}\unicode[STIX]{x2202}_{\bar{r}}(\bar{r}\bar{\unicode[STIX]{x1D70F}}_{11})+\unicode[STIX]{x2202}_{\bar{z}}\bar{\unicode[STIX]{x1D70F}}_{12},\qquad\end{eqnarray}$$

(4.1b ) $$\begin{eqnarray}\displaystyle \bar{\unicode[STIX]{x1D70C}}^{(g)}(\bar{U}\bar{W}_{\bar{r}}+\bar{W}\bar{W}_{\bar{z}}) & = & \displaystyle -\bar{P}_{\bar{z}}+\bar{\unicode[STIX]{x1D707}}^{(g)}\left(\frac{1}{\bar{r}}\unicode[STIX]{x2202}_{\bar{r}}(\bar{r}\unicode[STIX]{x2202}_{\bar{r}}\bar{W})+\unicode[STIX]{x2202}_{\bar{z}\bar{z}}\bar{W}\right)+\frac{1}{\bar{r}}\unicode[STIX]{x2202}_{\bar{r}}(\bar{r}\bar{\unicode[STIX]{x1D70F}}_{12})+\unicode[STIX]{x2202}_{\bar{z}}\bar{\unicode[STIX]{x1D70F}}_{22},\end{eqnarray}$$

(4.1c ) $$\begin{eqnarray}\displaystyle \frac{1}{\bar{r}}\unicode[STIX]{x2202}_{\bar{r}}(\bar{r}\bar{U})+\unicode[STIX]{x2202}_{\bar{z}}\bar{W} & = & \displaystyle 0,\end{eqnarray}$$

$(\bar{U},\bar{W})$

$(\bar{r},\bar{z})$

$\bar{z}$

$\bar{r}$

$\bar{r}=\bar{R}(\bar{z})$

$1$

$2$

$\tilde{r}$

$\tilde{z}$

$r$

$z$

$\bar{r}=\bar{R}(\bar{z})$

(4.2) $$\begin{eqnarray}\displaystyle \bar{U}=\bar{W}=0, & & \displaystyle\end{eqnarray}$$

and at the centre of the tube $\bar{r}=0$ we require that the radial derivative of velocity vanish,

(4.3) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x2202}_{\bar{r}}\bar{W}=0. & & \displaystyle\end{eqnarray}$$

The problem can be non-dimensionalized by the following reference scales:

(4.4) $$\begin{eqnarray}\displaystyle \left.\begin{array}{@{}c@{}}\displaystyle r=\bar{r}/\bar{R}_{P},\quad z=\unicode[STIX]{x1D716}_{P}\bar{z}/\bar{R}_{P},\quad U=\bar{U}/\unicode[STIX]{x1D716}_{P}\bar{W}_{P},\\ \displaystyle W=\bar{W}/\bar{W}_{P},\quad P=\unicode[STIX]{x1D716}_{P}\bar{P}/|\bar{T}_{P}|,\quad \unicode[STIX]{x1D70F}=\bar{\unicode[STIX]{x1D70F}}/|\bar{T}_{P}|,\end{array}\right\} & & \displaystyle\end{eqnarray}$$

where $\unicode[STIX]{x1D716}_{P}=R_{0}/\unicode[STIX]{x1D706}$ is the aspect ratio, and the velocity scale

(4.5) $$\begin{eqnarray}\displaystyle \bar{W}_{P}=\frac{\bar{Q}_{P}}{\unicode[STIX]{x03C0}\bar{R}_{0}^{2}}\sqrt{\frac{F_{P}}{2}}, & & \displaystyle\end{eqnarray}$$

is the friction velocity estimated using the empirical skin friction coefficient

(4.6a,b ) $$\begin{eqnarray}\displaystyle F_{P}=0.079(Re_{P}^{(g)})^{-0.25},\quad Re_{P}^{(g)}=\frac{2\bar{\unicode[STIX]{x1D70C}}^{(g)}\bar{Q}_{P}}{\unicode[STIX]{x03C0}\bar{\unicode[STIX]{x1D707}}^{(g)}\bar{R}_{0}}; & & \displaystyle\end{eqnarray}$$

$Re_{P}$ is the bulk diameter Reynolds number and $\bar{Q}_{P}$ is the prescribed volume flux of air. Subscripts of $P$ are used to denote the case of pipe flow. The slight difference in skin friction coefficient $F_{P}$ from that of the channel flow case, i.e. $F_{C}$ , is due to the cylindrical geometry of the pipe flow considered in this section. The length and stress/pressure scales are

(4.7a,b ) $$\begin{eqnarray}\displaystyle \bar{R}_{P}=\frac{\bar{\unicode[STIX]{x1D707}}^{(g)}}{\bar{\unicode[STIX]{x1D70C}}^{(g)}\bar{W}_{P}},\quad |\bar{T}_{P}|=\bar{\unicode[STIX]{x1D70C}}^{(g)}(\bar{W}_{P})^{2}=\frac{\bar{\unicode[STIX]{x1D70C}}^{(g)}F_{P}\bar{Q}_{P}^{2}}{2\unicode[STIX]{x03C0}^{2}\bar{R}_{0}^{4}}. & & \displaystyle\end{eqnarray}$$

Here $|\bar{T}_{P}|$ is an empirical estimate of the wall stress in flow through a flat tube of radius $\bar{R}_{0}$ .

With the above scaling, the wall is at $r=R(z)=R_{0}+\unicode[STIX]{x1D6FC}(z)$ , where $R_{0}=\bar{R}_{0}/\bar{R}_{P}$ is the dimensionless tube radius and $\unicode[STIX]{x1D6FC}=\bar{\unicode[STIX]{x1D6FC}}/\bar{R}_{P}$ measures the local constriction of the tube; thus, $\unicode[STIX]{x1D6FC}(z)<0$ indicates a narrow region of the tube while $\unicode[STIX]{x1D6FC}(z)>0$ indicates a relatively wide location in the tube. Note that unlike the case of channel flow, the bulk Reynolds number here does depend on the average tube radius $\bar{R}_{0}$ . We will not take this dependence into account in the pipe flow model and will assume that $\unicode[STIX]{x1D6FC}(z)/R_{0}\ll 1$ for all $z$ .

Substituting (4.4) into (4.1) results in

(4.8a ) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D716}_{P}^{3}(UU_{r}+WU_{z}) & = & \displaystyle -P_{r}+\frac{\unicode[STIX]{x1D716}_{P}^{2}}{r}\unicode[STIX]{x2202}_{r}(r\unicode[STIX]{x2202}_{r}U)+\unicode[STIX]{x1D716}_{P}^{4}U_{zz}+\frac{\unicode[STIX]{x1D716}_{P}}{r}\unicode[STIX]{x2202}_{r}(r\unicode[STIX]{x1D70F}_{11})+\unicode[STIX]{x1D716}_{P}^{2}\unicode[STIX]{x2202}_{z}\unicode[STIX]{x1D70F}_{12},\end{eqnarray}$$

(4.8b ) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D716}_{P}(UW_{r}+WW_{z}) & = & \displaystyle -P_{z}+\frac{1}{r}\unicode[STIX]{x2202}_{r}(r\unicode[STIX]{x2202}_{r}W)+\unicode[STIX]{x1D716}_{P}^{2}W_{zz}+\frac{1}{r}\unicode[STIX]{x2202}_{r}(r\unicode[STIX]{x1D70F}_{12})+\unicode[STIX]{x1D716}_{P}\unicode[STIX]{x2202}_{z}\unicode[STIX]{x1D70F}_{22},\end{eqnarray}$$

(4.8c ) $$\begin{eqnarray}\displaystyle \frac{1}{r}\unicode[STIX]{x2202}_{r}(rU)+W_{z} & = & \displaystyle 0.\end{eqnarray}$$

4.2 Coordinate transformation

We now undertake a coordinate transformation $(r,z)\rightarrow (\tilde{r},\tilde{z})$ analogous to the one made for flow through a wavy channel in § 3, i.e.

(4.9) $$\begin{eqnarray}\displaystyle \tilde{r} & = & \displaystyle \frac{R_{0}}{R(z)}r,\end{eqnarray}$$

so that $\tilde{r}=0$ coincides with the centre of the tube $r=0$ and $\tilde{r}=R_{0}$ coincides with the wavy tube wall $r=R$ .

The details of how the differential operators and flow variables of (4.8) are expressed in the new coordinate system are very similar to the case of wavy channel flow and are omitted. As before, we make a long-wave assumption, i.e. $\unicode[STIX]{x1D716}_{P}\ll 1$ , and expand the flow variables in $\unicode[STIX]{x1D716}_{P}$ as in (B 1) and seek velocity components of the form

(4.10a,b ) $$\begin{eqnarray}\displaystyle & \displaystyle \tilde{U} _{0}(\tilde{r},\tilde{z})=0,\quad \tilde{W}_{0}(\tilde{r},\tilde{z})=\tilde{w}_{0}(\tilde{r})g_{0}(\tilde{z}), & \displaystyle\end{eqnarray}$$

(4.10c ) $$\begin{eqnarray}\displaystyle & \displaystyle \tilde{W}_{1}(\tilde{r},\tilde{z})=\tilde{w}_{1}(\tilde{r})g_{1}(\tilde{z}). & \displaystyle\end{eqnarray}$$

As with the channel flow case, we consider the case where the normal Reynolds stresses are $O(\unicode[STIX]{x1D716}_{P})$ , i.e. $\tilde{\unicode[STIX]{x1D70F}}_{ii_{0}}=0$ for $i=1,2$ , and adopt a turbulent viscosity model for the Reynolds shear stress $\tilde{\unicode[STIX]{x1D70F}}_{12}$ given by (B 2); it is assumed that the turbulent viscosity in this new coordinate system can be modelled by the mixing-length model in (B 3)–(B 4).

4.3 Leading-order solution and first-order corrections

Substitution of (B 1)–(B 4) and (4.10) into the dimensionless governing system (4.8) results in the leading-order equations

(4.11a ) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x2202}_{\tilde{r}}\tilde{P}_{0} & = & \displaystyle 0,\end{eqnarray}$$

(4.11b ) $$\begin{eqnarray}\displaystyle \tilde{r}\unicode[STIX]{x2202}_{\tilde{z}}\tilde{P}_{0} & = & \displaystyle \frac{\tilde{g}_{0}R^{2}}{R_{0}^{2}}\frac{\text{d}}{\text{d}\tilde{r}}\left[\tilde{r}\tilde{w}_{0}^{\prime }\right]+\frac{\tilde{g}_{0}^{2}R}{R_{0}}\frac{\text{d}}{\text{d}\tilde{r}}\left[\tilde{r}\,\tilde{l}(\tilde{r})^{2}(\tilde{w}_{0}^{\prime })^{2}\right],\end{eqnarray}$$

(4.11c ) $$\begin{eqnarray}\displaystyle \left(\frac{2\tilde{g}_{0}}{R}\frac{\text{d}R}{\text{d}\tilde{z}}+\frac{\text{d}\tilde{g}_{0}}{\text{d}\tilde{z}}\right)\tilde{w}_{0} & = & \displaystyle 0,\end{eqnarray}$$

$\unicode[STIX]{x2202}_{\tilde{r}}\tilde{P}_{0}$

$\text{d}/\text{d}\tilde{r}$

(4.12a ) $$\begin{eqnarray}\displaystyle \tilde{w}_{0}=0\ & & \displaystyle \text{at }\tilde{r}=0,\end{eqnarray}$$

(4.12b ) $$\begin{eqnarray}\displaystyle \tilde{w}_{0}^{\prime }=0\ & & \displaystyle \text{at }\tilde{r}=R_{0}.\end{eqnarray}$$

The base flow is modulated by the slowly varying $\tilde{z}$ -dependent terms

(4.13a,b ) $$\begin{eqnarray}\displaystyle \tilde{g}_{0}(\tilde{z})=\frac{R_{0}^{2}}{R^{2}},\quad \unicode[STIX]{x2202}_{\tilde{z}}\tilde{P}_{0}=\frac{2R_{0}^{3}}{R^{4}}\left[\tilde{w}_{0}^{\prime }(R_{0})\right]. & & \displaystyle\end{eqnarray}$$

The $\tilde{r}$ -component $\tilde{w}_{0}$ of the leading-order velocity is identical to that found for channel flow, but note that the $z$ -dependence (4.13) of the base flow through a wavy pipe (with scale dictated by ${\sim}1/R^{2}$ ) is different than that of the channel flow case (3.18) (with scale ${\sim}1/L$ ) due to the cylindrical geometry of the problem. This leading-order flow results in a dimensional leading-order wall stress

(4.14) $$\begin{eqnarray}\displaystyle \bar{\unicode[STIX]{x1D70F}}_{P_{0}}^{(g)} & = & \displaystyle -\frac{|\bar{T}_{P}|\bar{R}_{0}^{3}}{\bar{R}^{3}};\end{eqnarray}$$

this will again be referred to as the LP model for the wall stress.

The first-order stress corrections may be found in a manner analogous to the channel flow; the details are given in appendix B. The combined stress $\bar{\unicode[STIX]{x1D70F}}_{P}^{(g)}=\bar{\unicode[STIX]{x1D70F}}_{P_{0}}^{(g)}+\bar{\unicode[STIX]{x1D70F}}_{P_{1}}^{(g)}$ is

(4.15) $$\begin{eqnarray}\displaystyle \bar{\unicode[STIX]{x1D70F}}_{P}^{(g)} & = & \displaystyle |\bar{T}_{P}|\left(\frac{\bar{R}_{0}}{\bar{R}}\right)^{3}\left[-1+B_{P}^{\ast }\left(\frac{\bar{R}_{0}}{\bar{R}}\right)\bar{R}_{\bar{z}}\right].\end{eqnarray}$$

The first term in (4.15), i.e. the leading-order stress $\bar{\unicode[STIX]{x1D70F}}_{P_{0}}$ , is the LP modelled stress; the total stress (4.15) will be referred to as the first-order or long-wave model for the wall stress. Similar to the channel flow case, the LP stress is in phase with the wall while the first-order model (4.15) depends both on the wall and its gradient, although the exact form of this dependence is different than in the channel flow case. Lastly, the dimensional pressure gradient at the wall through first order, i.e. $\unicode[STIX]{x2202}_{\bar{z}}\bar{P}_{P}^{(g)}=\unicode[STIX]{x2202}_{\bar{z}}\bar{P}_{0}+\unicode[STIX]{x2202}_{\bar{z}}\bar{P}_{1}$ , is

(4.16) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x2202}_{\bar{z}}\bar{P}_{P}^{(g)} & = & \displaystyle \frac{2|\bar{T}_{P}|}{\bar{R}_{0}}\left(\frac{\bar{R}_{0}}{\bar{R}}\right)^{4}\left[-1+(B_{P}+A_{P})\left(\frac{\bar{R}_{P}}{\bar{R}}\right)\frac{\text{d}\bar{R}}{\text{d}\bar{z}}\right].\end{eqnarray}$$

4.4 Comparison with experiments and other models

Once again we compare the model predictions for the wall stress with experiments. Two sinusoidally varying pipes were used in Hsu & Kennedy (Reference Hsu and Kennedy1971) with different amplitudes, and air was forced to flow at two different Reynolds numbers. Figure 8 shows a comparison of the model with two of the parameter settings from Hsu & Kennedy (Reference Hsu and Kennedy1971) with mean tube radius $0.204$ ft, wall wavelength $1.667$ ft and wall wave amplitude $0.0184$ ft (crest-to-trough amplitude $0.0368$ ft). The amplitude of the stress corrections seen in experiments is matched well by the first-order model, though the model overstates the phase shift by a considerable amount. The reasons why the model is able to more accurately capture the behaviour of wavy channel flow experiments in § 3 than those in pipes such as the ones here are not immediately clear, but there are several possible explanations. For example, both the Reynolds number and wave amplitude of the wavy pipe flow experiments are substantially higher than those of the channel flow experiments. In addition, the first-order model was derived while neglecting variations in the local Reynolds number due to the wall waviness. It is also possible that use of a different empirical formula other than (B 3)–(B 4) for the mixing length ratio may result in improved agreement between model and experiment. Nevertheless, the model appears to be an improvement over the LP model in all cases, and it will be incorporated into the liquid model next.

Figure 8.

Comparison of locally Poiseuille (dashed line) and first-order (solid line) model predictions of wall shear stress for flow through a wavy pipe with experiments in Hsu & Kennedy (Reference Hsu and Kennedy1971) (circles). $Re_{P}^{(g)}=90\,000$ , $\bar{a}=0.204$ ft, wall wavelength is $1.667$ ft and wall wave amplitude is $0.0184$ ft (crest-to-trough amplitude $0.0368$ ft); the ratio of amplitude to wavelength has been enlarged in the wavy wall shown in the bottom panel.

In closing this section, we remark that the nonlinear response in the stress is visible in the model estimates in figure 8. The dependence on the axial coordinate of the wall shear stress is not given by a simple sinusoid; the distance when moving upstream from a crest to a trough in the wall shear stress is smaller than the distance moving upstream from a trough to a crest, while the wavy wall is sinusoidally varying.

5.1 Governing equations

The governing equations for the flow of the liquid film inside a vertical tube are the axisymmetric Navier–Stokes equations in cylindrical coordinates (2.1), made dimensionless here using the scalings

(5.1) $$\begin{eqnarray}\displaystyle \left.\begin{array}{@{}rcl@{}}r\ & =\ & \bar{r}/\bar{h}_{0},\quad z=\bar{z}/\bar{\unicode[STIX]{x1D706}},\quad \,u=\bar{u}/\bar{U}_{0},\\ w\ & =\ & \bar{w}/\bar{W}_{0},\quad t=\bar{t}\bar{W}_{0}/\bar{\unicode[STIX]{x1D706}},\quad p=\bar{p}\bar{h}_{0}/(\unicode[STIX]{x1D707}^{(l)}\bar{W}_{0}),\end{array}\right\} & & \displaystyle\end{eqnarray}$$

where $\bar{h}_{0}$ is the film’s mean thickness, $\bar{W}_{0}=g\unicode[STIX]{x1D70C}^{(l)}\bar{h}_{0}^{2}/\unicode[STIX]{x1D707}^{(l)}$ is the axial velocity scale, $\bar{U}_{0}=\unicode[STIX]{x1D716}\bar{W}_{0}$ is the radial velocity scale as determined by the continuity equation and $\unicode[STIX]{x1D716}=\bar{h}_{0}/\bar{\unicode[STIX]{x1D706}}\ll 1$ is the traditional long-wave aspect ratio small parameter. Substitution of (5.1) into (2.1) for the liquid flow yields

(5.2a ) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D716}^{2}Re^{(l)}(u_{t}+uu_{r}+wu_{z}) & = & \displaystyle -p_{r}+\unicode[STIX]{x1D716}\left(\frac{1}{r}\unicode[STIX]{x2202}_{r}(ru_{r})+\unicode[STIX]{x1D716}^{2}u_{zz}-\frac{u}{r^{2}}\right),\end{eqnarray}$$

(5.2b ) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D716}Re^{(l)}(w_{t}+uw_{r}+ww_{z}) & = & \displaystyle -\unicode[STIX]{x1D716}p_{z}+\frac{1}{r}\unicode[STIX]{x2202}_{r}(rw_{r})+\unicode[STIX]{x1D716}^{2}w_{zz}-1,\end{eqnarray}$$

(5.2c ) $$\begin{eqnarray}\displaystyle \frac{1}{r}\unicode[STIX]{x2202}_{r}(ru)+w_{z} & = & \displaystyle 0,\end{eqnarray}$$

$\mathit{Re}^{(l)}\equiv \bar{\unicode[STIX]{x1D70C}}^{(l)}\bar{W}_{0}\bar{h}_{0}/\bar{\unicode[STIX]{x1D707}}^{(l)}$

$\bar{z}$

Boundary conditions at the wall, $r=a$ , are no slip,

(5.3) $$\begin{eqnarray}\displaystyle u=w=0, & & \displaystyle\end{eqnarray}$$

and at the free surface, $r=R(z,t)$ , are

(5.4a ) $$\begin{eqnarray}\displaystyle w_{r}+O(\unicode[STIX]{x1D716}^{2}) & = & \displaystyle \unicode[STIX]{x1D70F}^{(g)},\end{eqnarray}$$

(5.4b ) $$\begin{eqnarray}\displaystyle -\!p+O(\unicode[STIX]{x1D716}) & = & \displaystyle -p^{(g)}+\frac{1}{Bo}\left(\frac{1}{R}-\unicode[STIX]{x1D716}^{2}R_{zz}\right),\end{eqnarray}$$

(5.4c ) $$\begin{eqnarray}\displaystyle w & = & \displaystyle R_{t}+uR_{z},\end{eqnarray}$$

(5.5a,b ) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D70F}^{(g)}=\bar{\unicode[STIX]{x1D70F}}_{P}^{(g)}\bar{h}_{0}/(\unicode[STIX]{x1D707}^{(l)}\bar{W}_{0}),\quad p^{(g)}=\bar{P}_{P}^{(g)}\bar{h}_{0}/(\unicode[STIX]{x1D707}^{(l)}\bar{W}_{0}), & & \displaystyle\end{eqnarray}$$

and

(5.6) $$\begin{eqnarray}\displaystyle Bo & = & \displaystyle \bar{W}_{0}\unicode[STIX]{x1D707}^{(l)}/\unicode[STIX]{x1D6FE},\nonumber\\ \displaystyle & = & \displaystyle \frac{\unicode[STIX]{x1D70C}^{(l)}g\bar{a}^{2}}{\unicode[STIX]{x1D6FE}}\frac{1}{a^{2}},\end{eqnarray}$$

is the Bond number, where $\unicode[STIX]{x1D6FE}$ is surface tension. We note that a term of $O(\unicode[STIX]{x1D716}^{2})$ has been retained in (5.4b ) despite being formally of higher order; this term is traditionally kept in long-wave modelling to account for uniform matched-asymptotics validity (see, e.g. Oron, Davis & Bankoff (Reference Oron, Davis and Bankoff1997), and Craster & Matar (Reference Craster and Matar2009), for discussion of this term in models for related problems). Integrating the continuity equation across the annular-sectional area of the liquid and substituting (5.4c ) yields the layer-mean evolution equation

(5.7) $$\begin{eqnarray}\displaystyle R_{t}-\frac{1}{R}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}z}\int _{R}^{a}wr\,\text{d}r=0, & & \displaystyle\end{eqnarray}$$

for the location $R(z,t)$ of the free surface. All that is needed to close the model equation (5.7) is an estimate of the axial velocity $w$ .

5.2 Leading-order solution

The flow variables $u,w,p$ and stresses $p^{(g)}$ and $\unicode[STIX]{x1D70F}^{(g)}$ may be expanded in terms of the small parameter $\unicode[STIX]{x1D716}$ as

(5.8a ) $$\begin{eqnarray}\displaystyle u & = & \displaystyle u_{0}+\unicode[STIX]{x1D716}u_{1}+O(\unicode[STIX]{x1D716}^{2}),\end{eqnarray}$$

(5.8b ) $$\begin{eqnarray}\displaystyle w & = & \displaystyle w_{0}+\unicode[STIX]{x1D716}w_{1}+O(\unicode[STIX]{x1D716}^{2}),\end{eqnarray}$$

(5.8c ) $$\begin{eqnarray}\displaystyle p & = & \displaystyle p_{0}+\unicode[STIX]{x1D716}p_{1}+O(\unicode[STIX]{x1D716}^{2}),\end{eqnarray}$$

(5.8d ) $$\begin{eqnarray}\displaystyle p^{(g)} & = & \displaystyle p_{0}^{(g)}+\unicode[STIX]{x1D716}p_{1}^{(g)}+O(\unicode[STIX]{x1D716}^{2}),\end{eqnarray}$$

(5.8e ) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D70F}^{(g)} & = & \displaystyle \unicode[STIX]{x1D70F}_{0}^{(g)}+\unicode[STIX]{x1D716}\unicode[STIX]{x1D70F}_{1}^{(g)}+O(\unicode[STIX]{x1D716}^{2}).\end{eqnarray}$$

$O(1)$

(5.9a-c ) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x2202}_{r}p_{0}=0,\quad \frac{1}{r}\unicode[STIX]{x2202}_{r}(r\unicode[STIX]{x2202}_{r}w_{0})=1,\quad \unicode[STIX]{x2202}_{r}(ru_{0})=-r\unicode[STIX]{x2202}_{z}w_{0}, & & \displaystyle\end{eqnarray}$$

with leading-order boundary conditions at the wall $r=a$ given by

(5.10) $$\begin{eqnarray}\displaystyle u_{0}=w_{0}=0, & & \displaystyle\end{eqnarray}$$

and at the free surface $r=R(z,t)$ by

(5.11a,b ) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x2202}_{r}w_{0} & = & \displaystyle \unicode[STIX]{x1D70F}_{0}^{(g)},\quad p_{0}=p_{0}^{(g)}.\end{eqnarray}$$

The leading-order stress $\unicode[STIX]{x1D70F}_{0}^{(g)}$ imparted by the gas flow at the interface can be approximated using $\bar{\unicode[STIX]{x1D70F}}_{P_{0}}^{(g)}$ found in (4.14), here non-dimensionalized by the liquid scalings,

(5.12) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D70F}_{0}^{(g)} & = & \displaystyle -\frac{\tilde{T}_{w0}a^{3}(a-1)}{R^{3}},\end{eqnarray}$$

where

(5.13) $$\begin{eqnarray}\displaystyle \tilde{T}_{w0} & = & \displaystyle \frac{|\bar{T}_{P}|(a-1)^{3}}{\unicode[STIX]{x1D70C}^{(l)}g\bar{h}_{0}a^{3}}=\frac{0.079(\unicode[STIX]{x1D707}^{(g)})^{2}}{8g\unicode[STIX]{x1D70C}^{(g)}\unicode[STIX]{x1D70C}^{(l)}\bar{a}^{3}}(Re_{P}^{(g)})^{7/4},\end{eqnarray}$$

while the leading-order pressure gradient can be approximated by

(5.14) $$\begin{eqnarray}\displaystyle p_{0z}^{(g)}=\frac{2\unicode[STIX]{x1D70F}_{0}^{(g)}}{R} & = & \displaystyle -\frac{2\tilde{T}_{w0}a^{3}(a-1)}{R^{4}}.\end{eqnarray}$$

Substituting (5.12) into (5.11) and solving (5.9)–(5.11) results in the leading-order velocity profile

(5.15a ) $$\begin{eqnarray}\displaystyle w_{0} & = & \displaystyle \frac{1}{4}\left[r^{2}-a^{2}-2\left(R^{2}+\frac{2\tilde{T}_{w0}a^{3}(a-1)}{R^{2}}\right)\log \frac{r}{a}\right],\end{eqnarray}$$

(5.15b ) $$\begin{eqnarray}\displaystyle u_{0} & = & \displaystyle \frac{1}{4r}\left(a^{2}-r^{2}+2r^{2}\log \frac{r}{a}\right)\left[R-\frac{2\tilde{T}_{w0}a^{3}(a-1)}{R^{3}}\right]R_{z}.\end{eqnarray}$$

5.3 First-order corrections

Collecting all terms of $O(\unicode[STIX]{x1D716})$ in (5.2b ) results in

(5.16a ) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x2202}_{r}p_{1} & = & \displaystyle \frac{1}{r}\unicode[STIX]{x2202}_{r}(r\unicode[STIX]{x2202}_{r}u_{0})-\frac{u_{0}}{r^{2}},\end{eqnarray}$$

(5.16b ) $$\begin{eqnarray}\displaystyle \frac{1}{r}\unicode[STIX]{x2202}_{r}(r\unicode[STIX]{x2202}_{r}w_{1}) & = & \displaystyle \mathit{Re}^{(l)}(\unicode[STIX]{x2202}_{t}w_{0}+u_{0}\unicode[STIX]{x2202}_{r}w_{0}+w_{0}\unicode[STIX]{x2202}_{z}w_{0})+\unicode[STIX]{x2202}_{z}p_{0},\end{eqnarray}$$

(5.16c ) $$\begin{eqnarray}\displaystyle \frac{1}{r}\unicode[STIX]{x2202}_{r}(ru_{1})+\unicode[STIX]{x2202}_{z}w_{1} & = & \displaystyle 0.\end{eqnarray}$$

$Re^{(l)}=O(\unicode[STIX]{x1D716})$

$Re^{(l)}=O(10^{-2})$

$w_{1}$

(5.17) $$\begin{eqnarray}\displaystyle \frac{1}{r}\unicode[STIX]{x2202}_{r}(r\unicode[STIX]{x2202}_{r}w_{1})=\unicode[STIX]{x2202}_{z}p_{0}. & & \displaystyle\end{eqnarray}$$

At $O(\unicode[STIX]{x1D716})$ , the boundary conditions at the wall ( $r=a$ ) are

(5.18) $$\begin{eqnarray}\displaystyle u_{1}=w_{1}=0, & & \displaystyle\end{eqnarray}$$

and at the free surface ( $r=R(z,t)$ ) are

(5.19a,b ) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x2202}_{r}w_{1} & = & \displaystyle \unicode[STIX]{x1D70F}_{1}^{(g)},\quad p_{1}=p_{1}^{(g)}.\end{eqnarray}$$

The first-order shear stress corrections $\unicode[STIX]{x1D70F}_{1}^{(g)}$ can be approximated by $\bar{\unicode[STIX]{x1D70F}}_{P_{1}}^{(g)}$ in (B 7),

(5.20) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D70F}_{1}^{(g)} & = & \displaystyle \frac{\tilde{T}_{w1}a^{3}(a-1)^{2}}{R^{4}}R_{z},\end{eqnarray}$$

where

(5.21) $$\begin{eqnarray}\displaystyle \tilde{T}_{w1}=B_{P}^{\ast }\tilde{T}_{w0}, & & \displaystyle\end{eqnarray}$$

with $B_{P}^{\ast }$ defined in (B 8). Substituting (5.14) and (5.20) into (5.17) and (5.19), respectively, and solving (5.17)–(5.19) for $w_{1}$ results in

(5.22) $$\begin{eqnarray}\displaystyle w_{1} & = & \displaystyle \frac{1}{4}\left[-\frac{2\tilde{T}_{w0}a^{3}(a-1)}{R^{4}}R_{z}+\frac{1}{Bo}\left(\frac{R_{z}}{R^{2}}+\unicode[STIX]{x1D716}^{2}R_{zzz}\right)\right](r^{2}-a^{2})\nonumber\\ \displaystyle & & \displaystyle +\left[\frac{\tilde{T}_{w1}a^{3}(a-1)^{2}}{R^{3}}R_{z}+\frac{\tilde{T}_{w0}a^{3}(a-1)}{R^{2}}R_{z}-\frac{1}{Bo}(R_{z}+\unicode[STIX]{x1D716}^{2}R^{2}R_{zzz})\right]\log \frac{r}{a}.\qquad\end{eqnarray}$$

The first-order long-wave model then results from substitution of $w=w_{0}+w_{1}$ into (5.7). After adopting the convention of rescaling $z$ and $t$ by $\unicode[STIX]{x1D716}$ to return to the original aspect ratio, this model equation is

(5.23) $$\begin{eqnarray}\displaystyle & & \displaystyle R_{t}+\frac{1}{2}\left[\tilde{T}_{w0}a^{3}(a-1)f_{1}(R;a)-f_{2}(R;a)\right]R_{z}\nonumber\\ \displaystyle & & \displaystyle \quad +\,\frac{1}{16R}\left[16Sf_{3}(R;a)(R_{z}+R^{2}R_{zzz})+4\tilde{T}_{w1}a^{3}(a-1)^{2}f_{2}(R;a)R^{-3}R_{z}\right]_{z}=0,\qquad\end{eqnarray}$$

where

(5.24a ) $$\begin{eqnarray}\displaystyle f_{1}(R;a) & = & \displaystyle \frac{a^{2}}{R^{4}}\left(\frac{a^{2}}{R^{2}}-1\right),\end{eqnarray}$$

(5.24b ) $$\begin{eqnarray}\displaystyle f_{2}(R;a) & = & \displaystyle a^{2}-R^{2}+2R^{2}\log (R/a),\end{eqnarray}$$

(5.24c ) $$\begin{eqnarray}\displaystyle f_{3}(R;a) & = & \displaystyle \frac{a^{4}}{R^{2}}-4a^{2}+3R^{2}-4R^{2}\log (R/a),\end{eqnarray}$$

$S=1/16Bo$

The LP model studied in Camassa et al. (Reference Camassa, Forest, Lee, Ogrosky and Olander2012) is recovered if we set $\tilde{T}_{w1}=0$ , given here for later reference,

(5.25) $$\begin{eqnarray}R_{t}+\frac{1}{2}\left[\tilde{T}_{w0}a^{3}(a-1)f_{1}(R;a)-f_{2}(R;a)\right]R_{z}+\frac{1}{16R}\left[16Sf_{3}(R;a)(R_{z}+R^{2}R_{zzz})\right]_{z}=0.\end{eqnarray}$$

In the case of no imposed airflow, i.e. $\tilde{T}_{w0}=\tilde{T}_{w1}=0$ , the model of Part 1 (Camassa et al. Reference Camassa, Ogrosky and Olander2014) studied by Camassa et al. (Reference Camassa, Marzuola, Ogrosky and Vaughn2016) is recovered, and is similar in form to models for the flow of a film down the exterior of a tube (Lin & Liu Reference Lin and Liu1975; Kliakhandler, Davis & Bankoff Reference Kliakhandler, Davis and Bankoff2001; Craster & Matar Reference Craster and Matar2006).

The coefficients $\tilde{T}_{w0}$ and $\tilde{T}_{w1}$ are functions of the gas Reynolds number $Re_{P}^{(g)}$ (which itself depends on the film thickness parameter $a$ ). Thus the model is seen to depend on three parameters: a film thickness parameter $a$ , the gas Reynolds number $Re_{P}^{(g)}$ and Bond number $Bo$ . In experiments conducted in Camassa et al. (Reference Camassa, Forest, Lee, Ogrosky and Olander2012) the Bond number was identical in each realization of the experiment, and we will likewise focus on the model’s behaviour for various values of $a$ and $Re_{P}^{(g)}$ , keeping $Bo$ fixed in § 6.

We remark that the model retains the cylindrical geometry of the problem, a feature which has been shown to be important for certain types of model behaviour in related models (Camassa & Ogrosky Reference Camassa and Ogrosky2015). Note also that the free-surface evolution according to the long-wave model depends only on the local tube radius and its gradients.

5.4 Thin-film limit

The long-wave model (5.23) was derived without assuming the film’s thickness was small relative to the tube. If the film is taken to be thin, i.e. let $\unicode[STIX]{x1D6FD}=1/a\ll 1$ , then (5.23) can be simplified by expanding each function $f_{i}$ in (5.24) in terms of $\unicode[STIX]{x1D6FD}$ ,

(5.26a ) $$\begin{eqnarray}\displaystyle f_{1} & = & \displaystyle 2h\unicode[STIX]{x1D6FD}^{3}+O(\unicode[STIX]{x1D6FD}^{4}),\end{eqnarray}$$

(5.26b ) $$\begin{eqnarray}\displaystyle f_{2} & = & \displaystyle 2h^{2}+O(\unicode[STIX]{x1D6FD}),\end{eqnarray}$$

(5.26c ) $$\begin{eqnarray}\displaystyle f_{3} & = & \displaystyle \frac{16h^{3}}{3}\unicode[STIX]{x1D6FD}+O(\unicode[STIX]{x1D6FD}^{2}),\end{eqnarray}$$

(5.27) $$\begin{eqnarray}\displaystyle h=\frac{1}{\unicode[STIX]{x1D6FD}}-R & & \displaystyle\end{eqnarray}$$

represents the thickness of the annular film. Substitution of (5.27) and (5.26) into (5.23) results in

(5.28) $$\begin{eqnarray}\displaystyle h_{t}+(\tilde{T}_{w0}h-\unicode[STIX]{x1D6FD}h^{2})h_{z}+\left[\frac{16\unicode[STIX]{x1D6FD}\tilde{S}}{3}h^{3}\left(h_{z}+\frac{h_{zzz}}{\unicode[STIX]{x1D6FD}^{2}}\right)+\frac{\tilde{T}_{w1}}{2}h^{2}h_{z}\right]_{z}=0, & & \displaystyle\end{eqnarray}$$

where $\tilde{S}=\unicode[STIX]{x1D6FD}^{2}S$ .

This thin-film limit ‘flattens’ the geometry by neglecting the local azimuthal curvature of the surface over scales of the film thickness. It is perhaps not surprising then that the $\tilde{T}_{w1}$ term in (5.28) is of a somewhat similar form to the corresponding term in the flat-film long-wave model in Tseluiko & Kalliadasis (Reference Tseluiko and Kalliadasis2011) (their equation (4.62)), with the important distinction that here the correction is specified by $h$ and its axial derivatives rather than global information through the Fourier modes of $h$ . These models are of course also distinct from one another in the destabilizing role played by capillary forces due to the azimuthal curvature in the cylindrical set-up considered here.

We note that several related thin-film models can be recovered from (5.28) under appropriate limits. If only the leading-order effects of the airflow are included (i.e. set $\tilde{T}_{w1}$ to zero) and the effects of gravity are considered negligible (i.e. the term $\unicode[STIX]{x1D6FD}h^{2}h_{z}$ in (5.28) is neglected), the model derived and studied by Kerchman (Reference Kerchman1995) is recovered. If, instead, there is no imposed airflow (i.e. $\tilde{T}_{w0}=\tilde{T}_{w1}=0$ ), the model derived by Frenkel (Reference Frenkel1992) and studied by Kalliadasis & Chang (Reference Kalliadasis and Chang1994) and Kerchman & Frenkel (Reference Kerchman and Frenkel1994) is recovered. Finally, if both gravity and the effects of airflow are considered negligible, the model derived by Hammond (Reference Hammond1983) and studied by Lister et al. (Reference Lister2006) is recovered.