1 Introduction
In the fluid mechanics literature, it is well known that similarity solutions can describe the time-dependent spreading of thin viscous films, which thus gives this nonlinear model problem great utility. A similarly instructive problem from the elasticity literature concerns the bending of a beam due to external forces and moments, which is described by the Euler–Bernoulli equation and is nonlinear for large changes in local orientation of the beam. It is then natural to couple these two classical prototype problems from the mechanics literature to consider how gravitational forces from a viscous film spreading over a flexible beam can deflect the beam and so modify the shape and propagation rate of the liquid film. We study this coupled fluid–elastic dynamics problem using experiments and theory and identify several distinct limits where there are similarity solutions for the spreading rate and the beam deformation.
The general topic of elastohydrodynamics concerns problems where fluid flow is coupled to the deformation of an elastic boundary (Dowson & Ehret Reference Dowson and Ehret1999; Gohar Reference Gohar2001). Examples include the flow-induced deformation of an elastic object or boundary during collision (Davis, Serayssol & Hinch Reference Davis, Serayssol and Hinch1986), droplet generation in a soft microfluidic device (Pang et al. Reference Pang, Kim, Liu and Stone2014), and the lift force on a sedimenting object generated by sliding motions accompanied by elastic deformation (Sekimoto & Leibler Reference Sekimoto and Leibler1993; Skotheim & Mahadevan Reference Skotheim and Mahadevan2005; Salez & Mahadevan Reference Salez and Mahadevan2015). There are many natural examples related to a local flow-induced deformation, e.g. ejection of fungal spores from an ascus (Fritz et al. Reference Fritz, Seminara, Roper, Pringle and Brenner2013), biological tribology (articular cartilage) (Mow, Ratcliffe & Poole Reference Mow, Ratcliffe and Poole1992), and raindrop impact on a leaf (Gart et al. Reference Gart, Mates, Megaridis and Jung2015; Gilet & Bourouiba Reference Gilet and Bourouiba2015). On the other hand, elastohydrodynamics also describes the movement of a flexible solid object interacting with a surrounding flow, for example a microswimmer (Wiggins et al. Reference Wiggins, Riveline, Ott and Goldstein1998; Tony, Lauga & Hosoi Reference Tony, Lauga and Hosoi2006), an elastic fibre in a microchannel (Wexler et al. Reference Wexler, Trinh, Berthet, Quennouz, du Roure, Huppert, Lindner and Stone2013), or a flapping flag (Shelley & Zhang Reference Shelley and Zhang2011).
Several previous studies have analysed the flow of a rivulet along a prescribed inclined or curved substrate, for example Duffy & Moffatt (Reference Duffy and Moffatt1995, Reference Duffy and Moffatt1997), Wilson & Duffy (Reference Wilson and Duffy2005), and Leslie, Wilson & Duffy (Reference Leslie, Wilson and Duffy2013). Here our focus is a situation where the substrate geometry is unknown in advance, and indeed is strongly coupled to the flow. In our recent study (Howell, Robinson & Stone Reference Howell, Robinson and Stone2013), we developed a two-dimensional model for steady gravity-driven thin-film flow over a flexible cantilever. In this paper, we analyse the flow of a liquid rivulet along a flexible narrow beam, extending our previous study to include time dependence and variations in the shape of the rivulet cross-section. We study theoretically and experimentally the time dependence of liquid propagation and beam deformation. The flow is modelled using lubrication theory and the substrate is modelled as an Euler–Bernoulli beam. The related problem of flow of a layer of viscous fluid below an elastic plate has been analysed for example by Flitton & King (Reference Flitton and King2004), Lister, Peng & Neufeld (Reference Lister, Peng and Neufeld2013), and Hewitt, Balmforth & De Bruyn (Reference Hewitt, Balmforth and De Bruyn2015), while flow over an elastic membrane without bending stiffness was studied theoretically and experimentally by Zheng, Griffiths & Stone (Reference Zheng, Griffiths and Stone2015).
The paper is organised as follows. In § 2 we present the experimental method and a large number of results for the beam deflection and rivulet propagation distance as functions of time. The experiments vary the bending modulus and length, width, and thickness of the beam, and the flow rate of the liquid. In § 3 we describe the governing equations and boundary conditions for the beam shape and the liquid film profile, demonstrating that the problems for the liquid spreading and the beam deformation are intimately coupled. We find that the dynamics generically falls into one of two regimes, namely a ‘small-deflection’ regime and a ‘large-deflection’ regime. We obtain similarity solutions to describe the time-dependent liquid propagation and the beam deflection for the different regimes. We thus find three different power laws exhibited by the system during different time periods: (i) at early times when the liquid just begins to deform the beam; (ii) at intermediate times when the beam deflection increases rapidly in response to the weight of the liquid film; (iii) at late times when the beam has sagged close to vertical. We show that the experimental data collapse under scalings provided by the theoretical similarity solutions, and are then consistent with the theoretically predicted power laws. Finally, we discuss the results and draw conclusions in § 4.
2 Experiments
2.1 Experimental setup
We performed experiments for liquid flow over a flexible cantilever. The experimental setup is shown in figure 1. The end of a thin elastic beam was fixed at a wall and a constant flow rate was applied by a syringe pump (model: NE-1000, New Era Pump, USA). In this study, we considered the effects of varying the flow rate
$q$
, as well as the Young’s modulus
$E$
of the beam, and the beam shape (i.e. length
$L$
, width
$2b$
, and thickness
$w$
, as shown in figure 1). For the liquid, we used glycerol (VWR International), which has dynamic viscosity
${\it\mu}=1.0~\text{Pa}~\text{s}$
, density
${\it\rho}=1260~\text{kg}~\text{m}^{-3}$
, and surface tension
${\it\gamma}=62.0\pm 0.5~\text{mN}~\text{m}^{-1}$
. To clearly observe the liquid propagation during the experiment, we added a red food dye (Innovating Science) to the liquid. The physical properties of the final liquid were measured at room temperature (
$T=298$
K) with a rheometer (Anton-Paar MCR 301 with the CP 50 geometry) for the viscosity and with a conventional goniometer (Theta Lite, Biolin Scientific) for the surface tension.
Figure 1. Experimental apparatus. (a) Side view: a thin elastic beam of length
$L$
is fixed at the left-hand wall and a constant flow rate
$q$
is injected along the beam. The liquid-wetted length is
${\it\sigma}(t)$
and the deflection angle at the advancing front is
${\it\phi}(t)$
, where times
$t_{0}<t_{1}<t_{2}<t_{3}$
. Also,
${\it\psi}(s,t)$
is the local deformation angle, where
$s$
is arc length, while
$\boldsymbol{s}$
and
$\boldsymbol{n}$
are the unit tangent and normal vectors, respectively. (b) Front view: the cross-sectional shape of the liquid rivulet on the beam, where
$2b$
and
$w$
denote, respectively, the beam width and thickness;
$R$
is the radius of curvature of the liquid–air interface and
$2{\it\alpha}$
is the opening angle. The
$\hat{z}$
-axis is in the direction of the normal
$\boldsymbol{n}$
.
Polycarbonate (PC) and polyether ether ketone (PEEK) were used as the material for the beam. To vary the bending stiffness, we prepared various thicknesses (
$w=0.076{-}0.38$
mm) and widths (
$2b=3{-}8$
mm) of PC and PEEK materials (McMaster-Carr, NJ, USA). We obtained the Young’s modulus of each material by measuring the self-deflection of the beam due to its own weight (Crandall et al.
Reference Crandall, Lardner, Archer, Cook and Dahl1978). The Young’s moduli of PEEK and PC were measured as
$E\approx 2.4$
and 3.5 GPa, respectively, which are consistent with the physical property values of the materials provided by the vendor. The two materials were initially covered by a protective film; before each experiment we removed the protective film and the beam was rinsed with distilled water and dried with nitrogen gas.
The deformation of the beam by the flowing liquid was observed from the side and top views, as shown in figure 2, using two CMOS colour USB cameras (EO USB 2.0 with Nikon 1 V1 lens) with a frame rate of 1, 10, or 17 f.p.s., and a spatial resolution of
$1280\times 1024$
pixels. We measured the liquid propagation length
${\it\sigma}(t)$
and the deflection angle
${\it\phi}(t)$
at the advancing front, as defined in figure 1(a). To extract these quantities from the raw images, we performed image- and post-processing by using Matlab 2014a. We measured the evolution of
${\it\sigma}(t)$
and
${\it\phi}(t)$
up to the time when the liquid reached the end of the beam and began to drip.
Figure 2. Examples of side and top views of liquid flow over an elastic beam. (a) A small beam deflection case with
$E=2.4$
GPa,
$q=1.4\times 10^{-8}~\text{m}^{3}~\text{s}^{-1}$
,
$L=100$
mm,
$w=0.51$
mm, and
$2b=7$
mm. (b) A large beam deflection case with
$E=3.6$
GPa,
$q=2.2\times 10^{-8}~\text{m}^{3}~\text{s}^{-1}$
,
$L=50$
mm,
$w=0.076$
mm, and
$2b=4$
mm.
2.2 Experimental results
We investigated beam deformation and liquid propagation along the flexible beam while a constant flow rate is applied at the base. Two typical examples of how the beam deformation and liquid film evolve over time are displayed in figure 2 for two different values of the bending stiffness
$Ebw^{3}/6$
, namely (a)
$1.84\times 10^{-4}~\text{Pa}~\text{m}^{4}$
and (b)
$5.31\times 10^{-7}~\text{Pa}~\text{m}^{4}$
, respectively (see also supplementary movies 1 and 2 available at http://dx.doi.org/10.1017/jfm.2016.258). In case (a), the relatively stiff beam suffers only a small deflection, such that the angle
${\it\phi}(t)<{\rm\pi}/6$
up until the time when the liquid reaches the end of the beam; this is an example of what we refer to below as the ‘small-deflection’ regime. Figure 2(b) shows the evolution of a much less stiff beam, which soon sags until the deflection angle
${\it\phi}(t)$
approaches
${\rm\pi}/2$
and the liquid flow is close to vertical. Below we refer to this more dramatic behaviour as the ‘large-deflection’ regime.
For the small-deflection regime, we summarise experimental conditions and results as shown in figure 3. The flow rate
$q$
, Young’s modulus
$E$
, and the beam dimensions (
$L,2b,$
and
$w$
) are all varied, as listed in figure 3(a), while the bending stiffness in each case is sufficient to keep the deflection angle less than
${\rm\pi}/6$
throughout an experiment. Figure 3(b,c) shows the time evolution of the liquid propagation length
${\it\sigma}(t)$
and the deflection angle
${\it\phi}(t)$
. Initially,
${\it\phi}(t)$
remains close to zero, and the liquid spreads steadily, with
${\it\sigma}(t)$
apparently close to linear in
$t$
. However, the angle
${\it\phi}(t)$
then increases rapidly, which in turn causes a rapid acceleration in the front position
${\it\sigma}(t)$
.
Figure 3. Small-deflection results. (a) Experimental parameters. (b) Time evolution of the liquid propagation length
${\it\sigma}(t)$
(m). (c) Time evolution of the deflection angle
${\it\phi}(t)$
(rad) at the advancing front.
Next, we present experimental results of the large-deflection regime in figure 4. The flow rate and beam geometry are again varied, as shown in figure 4(a), and the corresponding time evolutions of
${\it\sigma}(t)$
and
${\it\phi}(t)$
are shown in figures 4(b) and 4(c), respectively. Compared with the results in figure 3, the beams used here are thinner such that the deflection angle exceeds
${\rm\pi}/6$
and, indeed, approaches
${\rm\pi}/2$
. In some cases, the beam is initially slightly deformed by its weight, and there is also an angle-measurement error of approximately
$3^{\circ }\approx 0.05$
rad. Thus, for some cases the beam deflection angle
${\it\phi}(t)$
appears to start from a non-zero value at
$t=0$
s.
Figure 4. Large-deflection results. (a) Experimental parameters. (b) Time evolution of the liquid propagation length
${\it\sigma}(t)$
(m). (c) Time evolution of the deflection angle
${\it\phi}(t)$
(rad) at the advancing front.
In the following section we present a theoretical model that describes the behaviour shown in figures 2–4 and allows the experimental results to be explained and quantitatively analysed.
3 Mathematical theory
3.1 Governing equations
We use Cartesian coordinates
$(x,z)$
as shown in figure 1, with the
$z$
-axis pointing vertically upwards and the beam clamped at
$x=0$
; the width of the beam lies in the
$y$
-direction. We parametrise the deformation of the beam in the
$(x,z)$
-plane using arc length
$s$
and time
$t$
, such that
$$\begin{eqnarray}\frac{\partial x}{\partial s}=\cos {\it\psi},\quad \frac{\partial z}{\partial s}=-\sin {\it\psi},\end{eqnarray}$$
where
${\it\psi}(s,t)$
is the local angle made by the beam with the
$x$
-axis (see the definitions in figure 1
a).
Let
$A(s,t)$
denote the cross-sectional area of a thin liquid film flowing over the top of the beam. A one-dimensional mass-conservation equation for the liquid is then
$$\begin{eqnarray}\frac{\partial A}{\partial t}+\frac{\partial Q}{\partial s}=0,\end{eqnarray}$$
where
$Q(s,t)$
is the flux of liquid along the beam. We assume that a known constant flux
$q$
is supplied at the upstream end, so that
$Q(0,t)\equiv q$
.
The tangential and normal components of the external force per unit length exerted on the beam are denoted by
$f_{s}$
and
$f_{n}$
. The Euler–Bernoulli equations governing the beam deformation are then given by
$$\begin{eqnarray}\frac{\partial T}{\partial s}+N\frac{\partial {\it\psi}}{\partial s}=-f_{s},\quad \frac{\partial N}{\partial s}-T\frac{\partial {\it\psi}}{\partial s}=-f_{n},\quad EI\frac{\partial ^{2}{\it\psi}}{\partial s^{2}}=N,\end{eqnarray}$$
where
$T$
and
$N$
are the tension and shear force in the beam, and
$$\begin{eqnarray}EI=\frac{Ebw^{3}}{6}\end{eqnarray}$$
is the bending stiffness.
To close the model, we need constitutive relations for the flux
$Q$
and the components of the force/length
$(f_{s},f_{n})$
in terms of
$A$
and
${\it\psi}$
. Our aim in this study is to find a tractable model that adequately captures the behaviour observed in experiments and is amenable to mathematical analysis. To this end, we make a number of assumptions to obtain relatively simple closed-form constitutive relations. First, we neglect the contribution of the beam’s own weight to the stress components
$f_{s}$
and
$f_{n}$
. In the experiments, the beam does sag somewhat by itself, e.g. see figure 4(c), but this self-induced deflection is small compared to the subsequent deflection once the fluid is injected, and we have found that including the weight of the beam in the theory makes very little difference to the results. We thus obtain the following expressions for the components of the force/length exerted on the beam by the fluid:
$$\begin{eqnarray}f_{s}={\it\rho}gA\sin {\it\psi}-A\frac{\partial P}{\partial s},\quad f_{n}=-{\it\rho}gA\cos {\it\psi},\end{eqnarray}$$
where
${\it\rho}$
is the density of the fluid,
$g$
is the acceleration due to gravity, and
$P(s,t)$
is the fluid pressure measured at the beam surface.
Constitutive relations relating the pressure
$P$
and flux
$Q$
to
$A(s,t)$
and
${\it\psi}(s,t)$
may be formally derived using lubrication theory in the limit where the fluid layer is relatively thin. The simplified relations
$$\begin{eqnarray}P(s,t)=\left(\frac{3{\it\gamma}}{2b^{3}}\right)A,\quad Q(s,t)=\frac{9A^{3}}{70{\it\mu}b^{2}}\left({\it\rho}g\sin {\it\psi}-\frac{3{\it\gamma}}{2b^{3}}\frac{\partial A}{\partial s}\right)\end{eqnarray}$$
are derived in appendix A in the asymptotic limit where the fluid layer is relatively thin and the Bond number,
$\mathit{B}$
, is small, i.e.
$$\begin{eqnarray}\frac{A}{b^{2}}\ll 1\quad \text{and}\quad \mathit{B}=\frac{{\it\rho}gb^{2}}{{\it\gamma}}\ll 1.\end{eqnarray}$$
It must be acknowledged that neither of these assumptions holds uniformly in the experiments. For example, based on the experimental conditions, we estimated that
$A/b^{2}\simeq 1$
and
$0.4\lesssim \mathit{B}\lesssim 3$
. Nevertheless, we believe that the approximations (3.6) are qualitatively reasonable and we will use them henceforth.
Combining (3.2), (3.3), (3.5), and (3.6), our final model equations are
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\partial A}{\partial t}+\frac{9}{70{\it\mu}b^{5}}\frac{\partial }{\partial s}\left[A^{3}\left({\it\rho}gb^{3}\sin {\it\psi}-\frac{3{\it\gamma}}{2}\frac{\partial A}{\partial s}\right)\right]=0, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\partial T}{\partial s}+N\frac{\partial {\it\psi}}{\partial s}+{\it\rho}gA\sin {\it\psi}-\frac{3{\it\gamma}A}{2b^{3}}\frac{\partial A}{\partial s}=0, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\partial N}{\partial s}-T\frac{\partial {\it\psi}}{\partial s}-{\it\rho}gA\cos {\it\psi}=0, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle EI\frac{\partial ^{2}{\it\psi}}{\partial s^{2}}=N, & \displaystyle\end{eqnarray}$$
$A$
,
${\it\psi}$
,
$T$
, and
$N$
. The corresponding boundary conditions are 
$$\begin{eqnarray}\displaystyle & \displaystyle A^{3}\frac{\partial A}{\partial s}+\frac{140{\it\mu}b^{5}q}{27{\it\gamma}}={\it\psi}=0\quad \text{ at }s=0, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle A=A^{3}\frac{\partial A}{\partial s}=N=T=\frac{\partial {\it\psi}}{\partial s}=0\quad \text{ at }s={\it\sigma}(t), & \displaystyle\end{eqnarray}$$
$s={\it\sigma}(t)$
denotes the moving front of the spreading rivulet. The conditions (3.9a
) arise from the prescribed flux
$q$
and horizontal clamping at
$s=0$
. The free-boundary conditions (3.9b
) arise from kinematic conditions for the liquid layer and from the imposition of no applied force or bending moment to the free end of the beam. The problem is closed by requiring the initial condition
${\it\sigma}(0)=0$
.
3.2 Small-deflection regime
3.2.1 Normalised problem
While the deflection angle
${\it\psi}$
is relatively small, the beam equations may be linearised and the problem (3.8) is then approximated by
$$\begin{eqnarray}\frac{\partial A}{\partial t}+\frac{9}{70{\it\mu}b^{5}}\frac{\partial }{\partial s}\left[A^{3}\left({\it\rho}gb^{3}{\it\psi}-\frac{3{\it\gamma}}{2}\frac{\partial A}{\partial s}\right)\right]=0,\quad EI\frac{\partial ^{3}{\it\psi}}{\partial s^{3}}={\it\rho}gA,\end{eqnarray}$$
where we have eliminated the force components
$T$
and
$N$
. The boundary conditions (3.9b
) in terms of
$A$
and
${\it\psi}$
are
$$\begin{eqnarray}\displaystyle & \displaystyle A^{3}\frac{\partial A}{\partial s}+\frac{140{\it\mu}b^{5}q}{27{\it\gamma}}={\it\psi}=0\quad \text{at }s=0, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle A=A^{3}\frac{\partial A}{\partial s}=\frac{\partial {\it\psi}}{\partial s}=\frac{\partial ^{2}{\it\psi}}{\partial s^{2}}=0\quad \text{at }s={\it\sigma}(t). & \displaystyle\end{eqnarray}$$
The simplified problem (3.10)–(3.11) may be normalised by defining the dimensionless variables
$$\begin{eqnarray}\displaystyle \tilde{A}=\left(\frac{729}{9800}\right)^{1/8}\left(\frac{{\it\rho}^{2}g^{2}{\it\gamma}^{3}}{{\it\mu}^{4}q^{4}Eb^{18}w^{3}}\right)^{1/16}A, & & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle \tilde{{\it\psi}}=\left(\frac{2}{11\,025}\right)^{1/8}\left(\frac{{\it\rho}^{10}g^{10}E^{3}b^{22}w^{9}}{{\it\mu}^{4}q^{4}{\it\gamma}^{9}}\right)^{1/16}{\it\psi}, & & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\tilde{{\it\sigma}}=\left(\frac{4{\it\rho}^{2}g^{2}b^{2}}{{\it\gamma}Ew^{3}}\right)^{1/4}{\it\sigma},\quad \tilde{s}=\left(\frac{4{\it\rho}^{2}g^{2}b^{2}}{{\it\gamma}Ew^{3}}\right)^{1/4}s,\end{eqnarray}$$
$$\begin{eqnarray}\tilde{t}=\left(\frac{1458}{1225}\right)^{1/8}\left(\frac{{\it\rho}^{10}g^{10}q^{12}}{{\it\mu}^{4}{\it\gamma}E^{5}b^{10}w^{15}}\right)^{1/16}t.\end{eqnarray}$$
The rescaled variables satisfy the problem (3.10)–(3.11) with all the coefficients equal to unity, i.e.
$$\begin{eqnarray}\displaystyle \frac{\partial \tilde{A}}{\partial \tilde{t}}+\frac{\partial }{\partial \tilde{s}}\left[\tilde{A}^{3}\left(\tilde{{\it\psi}}-\frac{\partial \tilde{A}}{\partial \tilde{s}}\right)\right]=0,\quad \frac{\partial ^{3}\tilde{{\it\psi}}}{\partial \tilde{s}^{3}}=\tilde{A}, & & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle \tilde{A}^{3}\frac{\partial \tilde{A}}{\partial \tilde{s}}+1=\tilde{{\it\psi}}=0\quad \text{at }\tilde{s}=0, & & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle \tilde{A}=\tilde{A}^{3}\frac{\partial \tilde{A}}{\partial \tilde{s}}=\frac{\partial \tilde{{\it\psi}}}{\partial \tilde{s}}=\frac{\partial ^{2}\tilde{{\it\psi}}}{\partial \tilde{s}^{2}}=0\quad \text{at }\tilde{s}=\tilde{{\it\sigma}}(t). & & \displaystyle\end{eqnarray}$$
In figure 5, we re-plot the small-deflection experimental results from figure 3 using the normalised variables (3.12), and demonstrate that there is indeed a reasonable collapse of the data.
Figure 5. Small-deflection experimental results from figure 3 plotted using the normalised variables defined in (3.12): (a)
$\tilde{{\it\sigma}}$
versus
$\tilde{t}$
and (b)
$\tilde{{\it\phi}}$
versus
$\tilde{t}$
. The asymptotic solutions in the small- and large-time limits are, respectively,
$\tilde{{\it\sigma}}\sim 1.00101\tilde{t}^{4/5}$
(solid),
$\tilde{{\it\sigma}}\sim 0.0434638\tilde{t}^{4}$
(dashed),
$\tilde{{\it\phi}}\sim 0.129117\tilde{t}^{13/5}$
(dash-dotted), and
$\tilde{{\it\phi}}\sim 0.000302834\tilde{t}^{9}$
(dash-double-dotted).
3.2.2 Small-time limit
As
$\tilde{t}\rightarrow 0$
, we expect
$\tilde{{\it\psi}}\ll |\partial \tilde{A}/\partial \tilde{s}|$
in (3.13a,b
). In this limit, the problem becomes mathematically equivalent to a classical gravity current on an effectively horizontal substrate (Huppert Reference Huppert1982b
). While a gravity current is driven by hydrostatic pressure proportional to film height, in the present problem an analogous role is played by the capillary pressure proportional to the cross-sectional area
$A$
. The corresponding behaviour of the solution to the problem (3.13) is described by a similarity solution of the form
$$\begin{eqnarray}\tilde{A}(\tilde{s},\tilde{t})=\tilde{t}^{1/5}f({\it\eta}),\quad \tilde{{\it\psi}}(\tilde{s},\tilde{t})=\tilde{t}^{13/5}g({\it\eta}),\quad {\it\eta}=\frac{\tilde{s}}{\tilde{t}^{4/5}},\end{eqnarray}$$
where
$f$
satisfies the ordinary differential equation (ODE)
$$\begin{eqnarray}f^{\prime \prime }+\frac{3(f^{\prime })^{2}}{f}+\frac{4{\it\eta}f^{\prime }}{5f^{3}}-\frac{1}{5f^{2}}=0\end{eqnarray}$$
and the boundary conditions
$$\begin{eqnarray}f^{3}(0)f^{\prime }(0)=-1,\quad f(c)=\lim _{{\it\eta}\rightarrow c}f^{3}({\it\eta})f^{\prime }({\it\eta})=0.\end{eqnarray}$$
The constant
$c$
is to be determined as part of the solution, and the position of the free boundary is then given by
$\tilde{{\it\sigma}}(\tilde{t})\sim c\tilde{t}^{4/5}$
as
$\tilde{t}\rightarrow 0$
. The deflection of the beam is determined a posteriori from
$$\begin{eqnarray}g({\it\eta})=\frac{1}{2}\int _{0}^{{\it\eta}}f({\it\xi}){\it\xi}^{2}\,\text{d}{\it\xi}+\frac{1}{2}\int _{{\it\eta}}^{c}f({\it\xi}){\it\eta}(2{\it\xi}-{\it\eta})\,\text{d}{\it\xi}.\end{eqnarray}$$
The numerical shooting technique used to solve this problem is outlined in § B.1, and the resulting solutions for
$f({\it\eta})$
and
$g({\it\eta})$
are plotted in figure 6. The area profile resembles a classical gravity current (Huppert Reference Huppert1982a
,Reference Huppert
b
), with a cube-root singularity at the moving touch-down location
${\it\eta}=c$
. From these solutions we read off the values
$c\approx 1.00101$
,
$f(0)\approx 1.32628$
, and
$g(c)\approx 0.129117$
. Hence, in the small-deflection regime, for small times the position of the advancing front and the maximum deflection angle at the front are given asymptotically by
$$\begin{eqnarray}\tilde{{\it\sigma}}(\tilde{t})\sim 1.00101\tilde{t}^{4/5},\quad \tilde{{\it\phi}}(\tilde{t})=\tilde{{\it\psi}}(\tilde{{\it\sigma}}(\tilde{t}),\tilde{t})\sim 0.129117\tilde{t}^{13/5}\quad \text{as }\tilde{t}\rightarrow 0.\end{eqnarray}$$
The predicted power laws (3.18) for
$\tilde{{\it\sigma}}(\tilde{t})$
and
$\tilde{{\it\phi}}(\tilde{t})$
are shown in figure 5, using solid and dash-dotted lines, respectively. There appears to be a good fit for the behaviour of
$\tilde{{\it\sigma}}$
, so long as the deflection angle remains small. The fit for
$\tilde{{\it\phi}}$
is also quite good for a range of intermediate times. The significant departures observed at very small values of
$\tilde{t}$
are due to the small initial deflection of the beam under its own weight, which is not included in our model, as well as angle-measurement errors, as explained in § 2.2.
3.2.3 Large-time limit
The limiting behaviour (3.18) describes the evolution while the beam deflection remains small enough to have a negligible influence on the spreading of the liquid. As
$\tilde{t}$
increases, the coupling between liquid flow and beam deformation becomes important. Eventually, as
$\tilde{t}\rightarrow \infty$
, the non-dimensional flux term in square brackets in (3.13) is dominated by
$\tilde{A}^{3}\tilde{{\it\psi}}$
. In this case the limiting behaviour is described by a similarity solution of the form
$$\begin{eqnarray}\tilde{A}(\tilde{s},\tilde{t})=\tilde{t}^{-3}f({\it\eta}),\quad \tilde{{\it\psi}}(\tilde{s},\tilde{t})=\tilde{t}^{9}g({\it\eta}),\quad {\it\eta}=\frac{\tilde{s}}{\tilde{t}^{4}},\end{eqnarray}$$
where
$f$
and
$g$
satisfy the ODEs
$$\begin{eqnarray}(f^{3}g-4{\it\eta}f)^{\prime }+f=0,\quad g^{\prime \prime \prime }-f=0.\end{eqnarray}$$
The corresponding boundary conditions, including the imposed flux, are
$$\begin{eqnarray}g({\it\eta})\rightarrow 0,\quad f({\it\eta})^{3}g({\it\eta})\rightarrow 1\quad \text{as }{\it\eta}\rightarrow 0,\quad g^{\prime }(c)=g^{\prime \prime }(c)=0.\end{eqnarray}$$
Again the constant
$c$
is to be determined as part of the solution, and the large-
$\tilde{t}$
behaviour of the free boundary is then given by
$\tilde{{\it\sigma}}(\tilde{t})\sim c\tilde{t}^{4}$
. To close the problem, we note that a constant liquid flux imposes the net conservation equation
$$\begin{eqnarray}\int _{0}^{c}f({\it\eta})\,\text{d}{\it\eta}=1.\end{eqnarray}$$
By integrating (3.20a
) with respect to
${\it\eta}$
, this integral condition may equivalently be stated as the boundary condition
$$\begin{eqnarray}f(c)^{2}g(c)=4c.\end{eqnarray}$$
The boundary-value problem (3.20)–(3.23) is solved using a shooting method outlined in § B.2, and the resulting solutions for
$f({\it\eta})$
and
$g({\it\eta})$
are plotted in figure 7. We note that
$f({\it\eta})$
decreases as
${\it\eta}$
increases from zero, attains a minimum value of approximately 20.3181 at
${\it\eta}\approx 0.0179634$
, and then increases again as
${\it\eta}$
approaches
$c$
. This behaviour reflects well the non-monotonic profiles for the film thickness observed in the experimental results, as shown in figure 2. However, the problem (3.20)–(3.23) predicts that
$f({\it\eta})\sim 3.64271{\it\eta}^{-1/3}$
as
${\it\eta}\rightarrow 0$
, implying that the cross-sectional area diverges towards the origin; also, we are unable to impose the condition
$f(c)=0$
corresponding to the condition
$\tilde{A}=0$
at the advancing front. Both of these apparent difficulties can be resolved by analysing asymptotic boundary layers near
$\tilde{s}=0$
and
$\tilde{s}=\tilde{{\it\sigma}}(\tilde{t})$
, as demonstrated in Howell et al. (Reference Howell, Robinson and Stone2013) for the steady version of the problem.
From the numerical solutions plotted in figure 7, we read off the values
$c\approx 0.0434638$
,
$f(c)\approx 23.9603$
,
$g^{\prime }(0)\approx 0.0206883$
, and
$g(c)\approx 0.000302834$
. Hence, in the small-deflection regime, for large times the position of the advancing front and the maximum deflection angle are given asymptotically by
$$\begin{eqnarray}\tilde{{\it\sigma}}(\tilde{t})\sim 0.0434638\tilde{t}^{4},\quad \tilde{{\it\phi}}(\tilde{t})\sim 0.000302834\tilde{t}^{9}\quad \text{as }\tilde{t}\rightarrow \infty .\end{eqnarray}$$
The power laws predicted in (3.24) are shown in figure 5, using dashed and dash-double-dotted lines, respectively. We observe that these power laws do give a reasonable fit to the dramatic increase in the deflection angle and consequent rapid movement of the rivulet along the beam.
3.3 Large-deflection regime
3.3.1 Normalised problem
The power laws (3.24) are valid in an intermediate regime where there is significant feedback between the beam deflection and the liquid flow, but the deflection angle remains relatively small. However, if the beam is sufficiently long, then the assumption that
${\it\psi}\ll 1$
must eventually fail, so that the nonlinear terms in
${\it\psi}$
that were neglected in the linearised problem (3.10) become significant. However, when
${\it\psi}=O(1)$
, the capillary terms involving spatial derivatives of
$A$
in the governing equations (3.8) become negligible compared with the gravitational terms (see Howell et al.
Reference Howell, Robinson and Stone2013), and the equations may be simplified to
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\partial A}{\partial t}+\frac{9{\it\rho}g}{70{\it\mu}b^{2}}\,\frac{\partial }{\partial s}(A^{3}\sin {\it\psi})=0, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\partial T}{\partial s}+N\frac{\partial {\it\psi}}{\partial s}+{\it\rho}gA\sin {\it\psi}=0, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\partial N}{\partial s}-T\frac{\partial {\it\psi}}{\partial s}-{\it\rho}gA\cos {\it\psi}=0, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle EI\frac{\partial ^{2}{\it\psi}}{\partial s^{2}}=N. & \displaystyle\end{eqnarray}$$
As in Howell et al. (Reference Howell, Robinson and Stone2013), a first integral of (3.25b,c ) allows us to write
$$\begin{eqnarray}T=F\sin {\it\psi},\quad N=-F\cos {\it\psi},\end{eqnarray}$$
where
$F$
is the vertical component of stress in the beam. The leading-order large-deflection equations (3.25b
)–(3.25d
) therefore reduce to
$$\begin{eqnarray}\frac{\partial F}{\partial s}=-{\it\rho}gA,\quad EI\frac{\partial ^{2}{\it\psi}}{\partial s^{2}}=-F\cos {\it\psi},\end{eqnarray}$$
which, with (3.25a
), form a closed system for
$A$
,
${\it\psi}$
, and
$F$
. The boundary conditions for
${\it\psi}$
and
$F$
are
$$\begin{eqnarray}{\it\psi}=0\quad \text{at }s=0,\quad F=\frac{\partial {\it\psi}}{\partial s}=0\quad \text{at }s={\it\sigma}(t),\end{eqnarray}$$
corresponding to horizontal clamping at
$s=0$
and zero applied force and moment at the free end of the beam.
Now that the highest spatial derivatives of
$A$
have been neglected, it is impossible to satisfy exactly the boundary conditions for
$A$
at
$s=0$
and
$s={\it\sigma}(t)$
. Instead, we impose the net flux conditions
$$\begin{eqnarray}\frac{9{\it\rho}gA^{3}\sin {\it\psi}}{70{\it\mu}b^{2}}\rightarrow q\quad \text{as }s\rightarrow 0,\quad \frac{\text{d}{\it\sigma}}{\text{d}t}=\frac{9{\it\rho}gA^{2}\sin {\it\psi}}{70{\it\mu}b^{2}}\quad \text{at }s={\it\sigma}(t).\end{eqnarray}$$
The full boundary conditions for
$A$
may be imposed by analysing asymptotic boundary layers near
$s=0$
and
$s={\it\sigma}(t)$
, in which the spatial derivatives of
$A$
regain their significance, as shown in Howell et al. (Reference Howell, Robinson and Stone2013).
Now the problem (3.25a ), (3.27)–(3.29) may be normalised by introducing the new dimensionless variables
$$\begin{eqnarray}\displaystyle \hat{A}=\left(\frac{9{\it\rho}g}{70{\it\mu}qb^{2}}\right)^{1/3}A,\quad \hat{{\it\sigma}}=\left(\frac{1680{\it\rho}^{2}g^{2}{\it\mu}q}{w^{9}bE^{3}}\right)^{1/9}{\it\sigma}, & & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle {\hat{s}}=\left(\frac{1680{\it\rho}^{2}g^{2}{\it\mu}q}{w^{9}bE^{3}}\right)^{1/9}s,\quad \hat{t}=\left(\frac{4374{\it\rho}^{5}g^{5}q^{7}}{1225{\it\mu}^{2}w^{9}b^{7}E^{3}}\right)^{1/9}t, & & \displaystyle\end{eqnarray}$$
with respect to which the governing equations read
$$\begin{eqnarray}\frac{\partial \hat{A}}{\partial \hat{t}}+\frac{\partial }{\partial {\hat{s}}}(\hat{A}^{3}\sin {\it\psi})=0,\quad \frac{\partial ^{2}{\it\psi}}{\partial {\hat{s}}^{2}}+\cos {\it\psi}\int _{{\hat{s}}}^{\hat{{\it\sigma}}}\hat{A}(s^{\prime },t)\,\text{d}s^{\prime }=0\end{eqnarray}$$
subject to
$$\begin{eqnarray}\displaystyle {\it\psi}\rightarrow 0,\quad \hat{A}^{3}\sin {\it\psi}\rightarrow 1 & & \displaystyle \text{as }{\hat{s}}\rightarrow 0,\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle \frac{\partial {\it\psi}}{\partial {\hat{s}}}=\frac{\text{d}\hat{{\it\sigma}}}{\text{d}\hat{t}}-\hat{A}^{2}\sin {\it\psi}=0 & & \displaystyle \text{at }{\hat{s}}=\hat{{\it\sigma}}(\hat{t}).\end{eqnarray}$$
${\it\psi}$
is
$O(1)$
, we expect the new scalings (3.30) to collapse the experimental data: this prediction will be confirmed below.
As
$\hat{t}\rightarrow 0$
, the solution of the problem (3.31)–(3.32) may be described by a similarity solution that is equivalent to the large-
$\tilde{t}$
solution (3.19). This result just confirms that the small- and large-deflection regimes are mutually consistent for intermediate values of
${\it\psi}$
.
3.3.2 Large-time limit
At large values of
$\hat{t}$
, assuming that the beam is sufficiently long, the weight of the fluid causes the beam to sag until it is approximately vertical. To study this limit, we write
${\it\psi}={\rm\pi}/2-{\it\chi}$
, where
$0<{\it\chi}\ll 1$
: it will transpire that
${\it\chi}$
is exponentially small. The governing equation (3.31a
) for
$\hat{A}$
thus becomes
$$\begin{eqnarray}\frac{\partial \hat{A}}{\partial \hat{t}}+\frac{\partial }{\partial {\hat{s}}}(\hat{A}^{3})=0,\end{eqnarray}$$
which is subject to
$$\begin{eqnarray}\hat{A}^{3}\rightarrow 1\quad \text{as }{\hat{s}}\rightarrow 0,\quad \frac{\text{d}\hat{{\it\sigma}}}{\text{d}\hat{t}}=\hat{A}^{2}\quad \text{at }{\hat{s}}=\hat{{\it\sigma}}(\hat{t}).\end{eqnarray}$$
The relevant large-
$\hat{t}$
limiting solution of the problem (3.33) is
$$\begin{eqnarray}\hat{A}({\hat{s}},\hat{t})=1,\quad \hat{{\it\sigma}}(\hat{t})=\hat{t}.\end{eqnarray}$$
With
$\hat{A}$
given by (3.34), the deflection equation (3.31b
) reduces to a form of the Airy equation for
${\it\chi}$
:
$$\begin{eqnarray}\frac{\partial ^{2}{\it\chi}}{\partial {\hat{s}}^{2}}=(\hat{{\it\sigma}}-{\hat{s}}){\it\chi}.\end{eqnarray}$$
Given
$\partial {\it\chi}/\partial {\hat{s}}=0$
at
${\hat{s}}=\hat{{\it\sigma}}$
, the solution of (3.35) is
$$\begin{eqnarray}{\it\chi}({\hat{s}},\hat{t})=\frac{3^{1/6}{\it\Gamma}(2/3)C}{2}[\text{Bi}(\hat{{\it\sigma}}-{\hat{s}})+\sqrt{3}\text{Ai}(\hat{{\it\sigma}}-{\hat{s}})],\end{eqnarray}$$
where
$\text{Ai}$
and
$\text{Bi}$
denote Airy functions and
$C(t)={\it\chi}(\hat{{\it\sigma}},\hat{t})$
is an arbitrary integration function, equal to the value of
${\it\chi}$
at the advancing front.
To determine
$C$
, and thus the deviation of the deflection from vertical, we have to match with an inner region near
${\hat{s}}=0$
in which
${\it\psi}$
rapidly adjusts from
$0$
to almost
${\rm\pi}/2$
. In this region, to lowest order the deflection equation (3.31b
) reduces to
$$\begin{eqnarray}\frac{\partial ^{2}{\it\psi}}{\partial {\hat{s}}^{2}}+\hat{{\it\sigma}}\cos {\it\psi}=0.\end{eqnarray}$$
The solution of (3.37) subject to the boundary and matching conditions
$$\begin{eqnarray}{\it\psi}=0\quad \text{at }{\hat{s}}=0,\quad {\it\psi}\rightarrow {\rm\pi}/2\quad \text{as }{\hat{s}}\rightarrow \infty\end{eqnarray}$$
is
$$\begin{eqnarray}{\it\psi}({\hat{s}},\hat{t})=\frac{{\rm\pi}}{2}-4\tan ^{-1}\left((\sqrt{2}-1)\text{e}^{-{\hat{s}}\sqrt{\hat{{\it\sigma}}}}\right).\end{eqnarray}$$
Finally, we get an expression for
$C$
by matching (3.39) with (3.36):
$$\begin{eqnarray}C=\frac{8\sqrt{{\rm\pi}}(\sqrt{2}-1)\hat{{\it\sigma}}^{1/4}\text{e}^{-2\hat{{\it\sigma}}^{3/2}/3}}{3^{1/6}{\it\Gamma}(2/3)}.\end{eqnarray}$$
In conclusion, when the beam sags to a nearly vertical configuration, we predict that the liquid front should grow linearly with time, i.e.
$\hat{{\it\sigma}}(\hat{t})\sim \hat{t}$
, and that the deflection angle and normalised free-boundary position should satisfy the relation
$$\begin{eqnarray}{\it\phi}\sim \frac{{\rm\pi}}{2}-\frac{8\sqrt{{\rm\pi}}(\sqrt{2}-1)\hat{{\it\sigma}}^{1/4}\text{e}^{-2\hat{{\it\sigma}}^{3/2}/3}}{3^{1/6}{\it\Gamma}(2/3)}.\end{eqnarray}$$
In figure 8(a), we re-plot the large-deflection results for
${\it\sigma}$
from figure 4(b) using the normalised variables defined in (3.30). We find that the data collapse onto a single curve, which agrees quite well with the linear behaviour predicted by (3.34), although with an
$O(1)$
disagreement in the prefactor. We discuss this disparity further in § 4. To test the predicted relation (3.41), in figure 8(b) we plot
$K-\log ({\rm\pi}/2-{\it\phi})$
versus
$\hat{{\it\sigma}}$
, where
$K$
is used as shorthand for the constant
$$\begin{eqnarray}K=\log \left(\frac{8\sqrt{{\rm\pi}}(\sqrt{2}-1)}{3^{1/6}{\it\Gamma}(2/3)}\right)\approx 1.284.\end{eqnarray}$$
Again we observe a dramatic collapse of the data in figure 8(b), as well as approximate convergence towards the asymptotic behaviour
$$\begin{eqnarray}K-\log \left(\frac{{\rm\pi}}{2}-{\it\phi}\right)\sim \frac{2}{3}\hat{{\it\sigma}}^{3/2}-\frac{1}{4}\log \hat{{\it\sigma}}\end{eqnarray}$$
corresponding to (3.41), which is indicated by a dashed curve.
Figure 8. Large-deflection results from figure 4 plotted using the normalised variables defined in (3.30). (a)
$\hat{{\it\sigma}}$
versus
$\hat{t}$
; the power law
$\hat{{\it\sigma}}=\hat{t}$
predicted by the large-time asymptotic analysis is indicated using a solid line. (b)
$K-\log ({\rm\pi}/2-{\it\phi})$
versus
$\hat{{\it\sigma}}$
, where
$K$
is defined by (3.42); the predicted behaviour (3.43) is indicated by the dashed curve.
4 Discussion and conclusions
We have studied both experimentally and theoretically the flow of a thin liquid rivulet along a flexible beam that is fixed at one end. The propagation of the liquid and the deflection of the beam are intimately coupled: the weight of the liquid causes the beam to bend, which, in turn, determines the effective body force driving the spreading of the liquid. Thus, this problem naturally combines two classical nonlinear mechanics problems in fluid mechanics and elasticity.
In analysing the problem mathematically, two distinct limits for the beam deflection were identified. In the ‘small-deflection’ limit, the contributions to the liquid flux from the slope of the beam and of the free surface are comparable, but the beam equations may be linearised. In the ‘large-deflection’ limit, the full nonlinear beam equations must be solved, but the liquid flux is dominated by the large beam slope. In either case, the mathematical model may be simplified and then made parameter free by a suitable normalisation. We demonstrated that the scalings thus predicted by the theory provide a very good collapse of a wide range of experimental data.
We found three distinct limiting solutions to the mathematical models obtained in the small- and large-deflection limits. The resulting power-law solutions for the position of the liquid front and the beam deflection are collected in table 1. The ‘small-time’ solution is valid while the beam deflection is so small as to have a negligible influence on the liquid, which therefore spreads as if on a horizontal substrate. The ‘intermediate-time’ solution occurs when the beam deflection is large enough to dominate the spreading of the liquid, but still small enough for the beam equations to be linearised. Finally, the ‘large-time’ solution emerges when the liquid has spread so far as to weigh the beam down almost to the vertical.
By comparison with the time scales used to normalise the problem in (3.12) and (3.30), we infer that the corresponding ranges for the dimensionless time
$t$
are given by
$$\begin{eqnarray}\displaystyle & \displaystyle \text{small time: }t\ll t_{\mathit{small}}=\left(\frac{{\it\mu}^{4}{\it\gamma}E^{5}b^{10}w^{15}}{{\it\rho}^{10}g^{10}q^{12}}\right)^{1/16}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \text{intermediate time: }t_{\mathit{small}}\ll t\ll t_{\mathit{large}}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \text{large time: }t\gg t_{\mathit{large}}=\left(\frac{{\it\mu}^{2}w^{9}b^{7}E^{3}}{{\it\rho}^{5}g^{5}q^{7}}\right)^{1/9}. & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\frac{t_{\mathit{small}}}{t_{\mathit{large}}}=\frac{{\it\gamma}^{1/16}{\it\mu}^{1/36}q^{1/36}}{{\it\rho}^{5/72}g^{5/72}b^{11/72}w^{1/16}E^{1/48}}.\end{eqnarray}$$
For the experimental parameter values, we find that
$t_{\mathit{small}}/t_{\mathit{large}}$
is in the range 0.6–0.85, that is, smaller than one but not very small. This perhaps helps to explain why the intermediate regime appears to persist only briefly in figure 5.
Table 1. Asymptotic solutions for liquid propagation along a flexible beam and the deformation of the beam. In this table,
${\it\sigma}$
and
${\it\phi}$
represent the length of the advancing liquid rivulet and the deflection angle of the beam, respectively, and
$t$
represents time.
Figures 5 and 8 demonstrate that the power laws listed in table 1 agree quite well with experimental results. However, there is some discrepancy in the prefactors. This is probably due to the simplified constitutive relations (3.6) for the liquid pressure and flux used in our mathematical analysis. The dramatic collapse of the experimental data and the apparent agreement with the predicted power-law exponents both support our claim that the relations (3.6) contain the relevant physics and exhibit the right qualitative behaviour. However, as pointed out in § 3.1, these relations are strictly valid only if the Bond number
$\mathit{B}$
and the ratio
$A/b^{2}$
are both small, neither of which is universally true in the experiments.
If the Bond number is not assumed to be small, then, under the lubrication approximation, the free surface of the liquid layer satisfies the Young–Laplace equation, balancing the capillary and hydrostatic pressures. Provided
$A/b^{2}$
is small, the relation between the base pressure
$P$
and the cross-sectional area
$A$
may then in principle be expressed in terms of hyperbolic functions (as in Paterson, Wilson & Duffy Reference Paterson, Wilson and Duffy2013). On the other hand, if
$A/b^{2}$
is not small, implying that the liquid layer is not thin, then in general the flux
$Q$
can only be found numerically, by solving Poisson’s equation for the liquid velocity along the beam. In principle, one can address each of these mathematical complications in a full computational solution of the problem, but it would seem to preclude any possibility of finding universal analytical predictions like those listed in table 1.
As shown in appendix A, one can relatively easily calculate the first corrections to the leading-order constitutive relations (3.6) when
$A/b^{2}$
and
$\mathit{B}$
are small but non-zero, namely
$$\begin{eqnarray}\displaystyle P(s,t) & {\sim} & \displaystyle \left(\frac{3{\it\gamma}}{2b^{3}}\right)A\left[1-\frac{27}{40}\,\frac{A^{2}}{b^{4}}+\frac{2}{5}\mathit{B}\cos {\it\psi}+\cdots \!\right],\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle Q(s,t) & {\sim} & \displaystyle \frac{9A^{3}}{70{\it\mu}b^{2}}\left({\it\rho}g\sin {\it\psi}-\frac{\partial P}{\partial s}\right)\left[1-\frac{2}{5}\,\frac{A^{2}}{b^{4}}-\frac{1}{45}\mathit{B}\cos {\it\psi}+\cdots \!\right].\end{eqnarray}$$
${\it\psi}\rightarrow 0$
, we therefore find that 
$$\begin{eqnarray}Q(s,t)\sim -\frac{27{\it\gamma}}{140{\it\mu}b^{5}}A^{3}\frac{\partial A}{\partial s}\left[1+\frac{17}{45}\mathit{B}-\frac{97}{40}\frac{A^{2}}{b^{4}}\right].\end{eqnarray}$$
Thus, inclusion of the transverse gravitational term proportional to
$\mathit{B}$
increases the spreading rate, while the geometric correction proportional to
$A^{2}/b^{4}$
decreases the spreading rate. It is conceivable that the combination of these effects could help to explain the discrepancy observed in figure 5(a), where the theory appears consistently to over-predict the spreading rate by a factor of 2–3. In the large-time regime where
${\it\psi}\rightarrow {\rm\pi}/2$
and the pressure gradient becomes negligible, we instead have
$$\begin{eqnarray}Q(s,t)\sim \frac{9{\it\rho}g}{70{\it\mu}b^{2}}A^{3}\left[1-\frac{2}{5}\frac{A^{2}}{b^{4}}\right].\end{eqnarray}$$
The leading-order term is equivalent to equation (1) of Wilson & Duffy (Reference Wilson and Duffy2005), and we observe that the geometric correction always decreases the spreading rate. This result is consistent with the observation in figure 8(a) that the simplified theory persistently over-predicts the spreading rate by a factor of around 5–10.
Finally, we note that the wettability of the substrate to the working fluid appears to give rise to a rather large advancing contact angle. In figure 2, for example, we observe a blunt free-surface profile and the formation of a noticeable bulge near the advancing front of the liquid film. Our simplified thin-film model is unlikely to capture accurately the quantitative behaviour of this localised structure. It may be that capillary effects near the advancing contact line limit the propagation of the front such that it lags behind the spreading rate of the thin film, resulting in accumulation of liquid into the observed bulge near the front.
Acknowledgement
We are very grateful to an anonymous referee, whose insightful suggestions resulted in significant improvements to this paper.
Supplementary movies
Supplementary movies are available at http://dx.doi.org/10.1017/jfm.2016.258.
Appendix A. Derivation of constitutive relations
Here we sketch the derivation of the constitutive relations (3.6) for the base pressure
$P$
and the flux
$Q$
in the rivulet. A schematic of the cross-section of the rivulet is shown in figure 1(b). The
$y$
- and
$\hat{z}$
-axes are parallel and normal respectively to the upper surface of the beam, which is at
$\hat{z}=0$
. Note the distinction between
$\hat{z}$
and the vertical coordinate
$z$
defined in figure 1(a); they are related by
$$\begin{eqnarray}z=-\int _{0}^{s}\sin {\it\psi}\,\text{d}s+\hat{z}\cos {\it\psi}.\end{eqnarray}$$
The free surface is denoted by
$\hat{z}=h(y)$
, where the parametric dependence upon time
$t$
and arc length
$s$
along the beam has been temporarily suppressed.
Under the assumptions of lubrication theory, the pressure in the rivulet is purely hydrostatic, and the free-surface profile
$h(y)$
satisfies the Young–Laplace equation
$$\begin{eqnarray}\frac{{\it\gamma}h^{\prime \prime }(y)}{[1+h^{\prime }(y)^{2}]^{3/2}}={\it\rho}gh(y)\cos {\it\psi}-P.\end{eqnarray}$$
The solution of (A 2) subject to
$h^{\prime }(0)=h(b)=0$
determines
$h(y)$
and hence
$$\begin{eqnarray}A=\int _{-b}^{b}h(y)\,\text{d}y\end{eqnarray}$$
in terms of
$P$
and
${\it\psi}$
; inversion of this relation then in principle gives
$P$
as a function of
$A$
and
${\it\psi}$
.
The velocity
$u$
in the
$s$
-direction satisfies Poisson’s equation in the form
$$\begin{eqnarray}{\it\mu}\left(\frac{\partial ^{2}u}{\partial y^{2}}+\frac{\partial ^{2}u}{\partial \hat{z}^{2}}\right)=\frac{\partial P}{\partial s}-{\it\rho}g\sin {\it\phi}.\end{eqnarray}$$
The imposition of zero slip at the base and a zero shear stress at the free surface leads to the boundary conditions
$$\begin{eqnarray}u=0\quad \text{at }\hat{z}=0,\quad \frac{\partial u}{\partial \hat{z}}-h^{\prime }(y)\frac{\partial u}{\partial y}=0\quad \text{at }\hat{z}=h(y).\end{eqnarray}$$
The solution of (A 4) subject to (A 5) in principle determines
$u$
and hence
$$\begin{eqnarray}Q=\int _{-b}^{b}\int _{0}^{h(y)}u(y,\hat{z})\,\text{d}\hat{z}\,\text{d}y\end{eqnarray}$$
in terms of
$A$
,
$P$
, and
${\it\psi}$
.
To obtain the simplified expressions (3.6), we assume that the rivulet is thin and that gravity is subdominant to surface tension, so that the cross-sectional Bond number is small. We formalise these assumptions by non-dimensionalising the above equations and boundary conditions as follows:
$$\begin{eqnarray}y=b{\tilde{y}},\quad \{\hat{z},h\}={\it\epsilon}b\{\tilde{z},\tilde{h}\},\quad P=\left(\frac{{\it\epsilon}{\it\gamma}}{b}\right)\tilde{P},\quad u=\frac{{\it\epsilon}^{2}b^{2}}{{\it\mu}}\left({\it\rho}g\sin {\it\phi}-\frac{\partial P}{\partial s}\right)\tilde{u} ,\end{eqnarray}$$
where
${\it\epsilon}\rightarrow 0$
in the limit of a thin rivulet. Henceforth the tildes will be dropped to reduce clutter. We also define
$$\begin{eqnarray}\mathit{B}\cos {\it\phi}={\it\epsilon}^{2}{\it\beta}\end{eqnarray}$$
and suppose that
${\it\beta}=O(1)$
as
${\it\epsilon}\rightarrow 0$
: this conveniently ensures that gravitational and geometric corrections enter at the same order.
The Young–Laplace equation (A 2) becomes
$$\begin{eqnarray}\frac{h^{\prime \prime }(y)}{[1+{\it\epsilon}^{2}h^{\prime }(y)^{2}]^{3/2}}={\it\epsilon}^{2}{\it\beta}h(y)-P,\end{eqnarray}$$
which is subject to
$h(\pm 1)=0$
. The cross-sectional area is then given by
$$\begin{eqnarray}\frac{A}{{\it\epsilon}b^{2}}=2\int _{0}^{1}h(y)\,\text{d}y.\end{eqnarray}$$
We then write
$h$
and
$P$
as asymptotic expansions in powers of
${\it\epsilon}^{2}$
, i.e.
$$\begin{eqnarray}h(y)\sim h_{0}(y)+{\it\epsilon}^{2}h_{1}(y)+\cdots \,,\quad P\sim P_{0}+{\it\epsilon}^{2}P_{1}+\cdots \,.\end{eqnarray}$$
Equation (A 9) may be solved successively for
$h_{0},h_{1},\ldots \,$
, and then the condition (A 10) determines
$P_{0},P_{1},\ldots \,$
. After halting this procedure at order
${\it\epsilon}^{2}$
and returning to dimensional variables, we find the approximation (4.3a
) for
$P$
. The first term corresponds to the model (3.6) used in the body of the paper. The following two terms are the first corrections arising from the nonlinear geometry and from gravity, respectively.
Next we solve for the normalised velocity
$u(y,z)$
, which satisfies the problem
$$\begin{eqnarray}\displaystyle \frac{\partial ^{2}u}{\partial z^{2}}+{\it\epsilon}^{2}\frac{\partial ^{2}u}{\partial y^{2}}=-1, & & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle u=0\quad \text{at }z=0,\quad \frac{\partial u}{\partial z}={\it\epsilon}^{2}h^{\prime }(y)\frac{\partial u}{\partial y}\quad \text{at }z=h(y). & & \displaystyle\end{eqnarray}$$
As above, we solve by writing
$u$
as an asymptotic expansion in powers of
${\it\epsilon}^{2}$
, and the normalised flux is then given by
$$\begin{eqnarray}Q=2\int _{0}^{1}\int _{0}^{h(y)}u(y,z)\,\text{d}z\,\text{d}y.\end{eqnarray}$$
We truncate the expansion at
$O({\it\epsilon}^{2})$
and return to dimensional variables to obtain the approximation (4.3b
) for
$Q$
. Again the leading term gives the model (3.6), and the subsequent terms give the first corrections in
$A/b^{2}$
and
$\mathit{B}$
.
Appendix B. Solution of numerical shooting problems
B.1 Small deflection, small
$\tilde{t}$
We have to solve the ODE (3.15) subject to the boundary conditions (3.16). We first make the problem autonomous via the transformation
$$\begin{eqnarray}{\it\eta}=c\text{e}^{-{\it\xi}},\quad f({\it\eta})={\it\eta}^{2/3}F({\it\xi}),\end{eqnarray}$$
so that
$F({\it\xi})$
satisfies the ODE
$$\begin{eqnarray}F^{\prime \prime }+\frac{3{F^{\prime }}^{2}}{F}-\frac{4F^{\prime }}{5F^{3}}-\frac{13F^{\prime }}{3}+\frac{10F}{9}+\frac{1}{3F^{2}}=0\end{eqnarray}$$
and the initial conditions
$$\begin{eqnarray}F({\it\xi})\rightarrow 0,\quad F^{\prime }({\it\xi})F({\it\xi})^{3}\rightarrow 0\quad \text{as }{\it\xi}\rightarrow 0.\end{eqnarray}$$
There is a unique solution of this initial-value problem, with the asymptotic behaviour
$$\begin{eqnarray}F({\it\xi})\sim \left(\frac{12{\it\xi}}{5}\right)^{1/3}\left\{1+\frac{47{\it\xi}}{96}+\frac{8983{\it\xi}^{2}}{64\,512}+\cdots \!\right\}\quad \text{as }{\it\xi}\rightarrow 0.\end{eqnarray}$$
We use this behaviour to integrate from a small positive value of
${\it\xi}$
. The initial condition
$f(0)^{3}f^{\prime }(0)=-1$
then allows us to determine both
$c$
and the value of
$f(0)$
from the far-field behaviour of
$F({\it\xi})$
, using
$$\begin{eqnarray}\text{e}^{-2{\it\xi}/3}F({\it\xi})\rightarrow c^{-2/3}f(0),\quad \text{e}^{-5{\it\xi}/3}F({\it\xi})^{3}\left(F^{\prime }({\it\xi})-{\textstyle \frac{2}{3}}F({\it\xi})\right)\rightarrow c^{-5/3}\quad \text{as }{\it\xi}\rightarrow \infty .\end{eqnarray}$$
We thus obtain the values
$f(0)\approx 1.32628$
and
$c\approx 1.00101$
. The normalised deflection angle
$g({\it\eta})$
is then determined by the integral (3.17), from which we find that
$g(0)\approx 0.129117$
.
The numerical solutions thus obtained for the functions
$f({\it\eta})$
and
$g({\it\eta})$
are plotted in figure 6.
B.2 Small deflection, large
$\tilde{t}$
The small-deflection large-
$\tilde{t}$
problem from § 3.2 leads to the system of ODEs (3.20) and boundary conditions (3.21), (3.23) for the similarity solution variables
$f({\it\eta})$
and
$g({\it\eta})$
. We now make the problem autonomous by defining
$$\begin{eqnarray}{\it\eta}=c\text{e}^{-{\it\xi}},\quad f({\it\eta})={\it\eta}^{-2/3}F({\it\xi}),\quad g({\it\eta})={\it\eta}^{7/3}G({\it\xi}),\end{eqnarray}$$
so that
$F$
and
$G$
satisfy the ODEs
$$\begin{eqnarray}F^{\prime }=\frac{F(1+3F^{2}G^{\prime }-F^{2}G)}{3(4-3F^{2}G)},\quad G^{\prime \prime \prime }-4G^{\prime \prime }+\frac{13}{3}G^{\prime }-\frac{28}{27}G+F=0,\end{eqnarray}$$
and boundary conditions
$$\begin{eqnarray}G^{\prime }(0)={\textstyle \frac{7}{3}}G(0),\quad G^{\prime \prime }(0)={\textstyle \frac{49}{9}}G(0),\quad F(0)=2G(0)^{-1/2}.\end{eqnarray}$$
The conditions (3.21) at
${\it\eta}=0$
transform to the far-field conditions
$$\begin{eqnarray}G({\it\xi})\sim g^{\prime }(0)c^{-4/3}\text{e}^{4{\it\xi}/3},\quad F({\it\xi})\sim g^{\prime }(0)^{-1/3}c^{1/3}\text{e}^{-{\it\xi}/3}\quad \text{as }{\it\xi}\rightarrow \infty .\end{eqnarray}$$
We therefore use
$G(0)$
as a shooting parameter to get
$$\begin{eqnarray}G^{\prime \prime }({\it\xi})-{\textstyle \frac{5}{3}}G^{\prime }({\it\xi})+{\textstyle \frac{4}{9}}G({\it\xi})\rightarrow 0\quad \text{as }{\it\xi}\rightarrow \infty\end{eqnarray}$$
(corresponding to
$G({\it\xi})\text{e}^{7{\it\xi}/3}\rightarrow 0$
), and then use (B 9) to infer the values of
$g^{\prime }(0)$
and
$c$
.
By following this procedure, we obtain the values
$$\begin{eqnarray}G(0)\approx 0.455938,\quad c\approx 0.0434638,\quad g^{\prime }(0)\approx 0.0206883.\end{eqnarray}$$
The corresponding value of the film area and the normalised angle at the advancing front are then given by
$$\begin{eqnarray}f(c)=2c^{-2/3}G(0)^{-1/2}\approx 23.9603,\quad g(c)=c^{7/3}G(0)\approx 0.000302834.\end{eqnarray}$$
The resulting numerical solutions for
$f({\it\eta})$
and
$g({\it\eta})$
are plotted in figure 7.
