The Moffatt-Pukhnachev flow: a new twist on an old problem

TL;DR

This paper investigates the dynamics of a viscous film on a rotating cylinder with oscillating angular velocity, revealing complex fractal structures, conditions for overturning, and asymptotic behaviors at frequency extremes.

Contribution

It introduces a numerical and analytical study of oscillating cylinder flows, uncovering fractal-like structures and conditions for film overturning and periodicity.

Findings

Fractal-like structures in amplitude-frequency space.

Overturning occurs at finite time for general initial conditions.

High-frequency oscillations can delay overturning and produce periodic solutions.

Abstract

The flow of a thin viscous film on the outside of a horizontal circular cylinder, whose angular velocity is time-periodic with specified frequency and amplitude, is investigated. The constant angular velocity problem was originally studied by Moffatt (1977) and Pukhnachev (1977). Surface tension is neglected. The evolution equation for the film thickness is solved numerically for a range of oscillation amplitudes and frequency. A blow-up map charted in amplitude-frequency space reveals highly intricate fractal-like structures exhibiting self-similarity. For a general initial condition the film surface reaches a slope singularity at a finite time and tends to overturn. The high-frequency and low-frequency limits are examined asymptotically using a multiple-scales approach. At high frequency the analysis suggests that an appropriate choice of initial profile can substantially delay the overturning time, and even yield a time-periodic solution. In the low-frequency limit it is possible to construct a quasi-periodic solution that does not overturn if the oscillation amplitude lies below a threshold value. Above this value the solution tends inexorably toward blow-up. It is shown how solutions exhibiting either a single-shock or a double-shock may be constructed in common with the steadily rotating cylinder problem.

Equations

$$r=a+ϵa{G}^{-1/2}\hat{y}.$$

$$u\mapsto {ϵ}^3{G}^{-1/2}a{\Omega }_{\ast }u,v\mapsto {ϵ}^{2}a{\Omega }_{\ast }\hat{v},p\mapsto \rho gap,t\mapsto {ϵ}^{-2}{\Omega }_{\ast }^{-1}t.$$

$$0=p_y,0=-p_{\theta }-\cos \theta +v_{yy},0=u_y+v_{\theta }.$$

$$h_t+v{h}_{\theta }-u=0,$$

$$p=0,v_y=0,$$

$$u=0,v=\Omega (t),$$

$$\Omega (t)=1+b\cos \sigma t$$

$$h_t+Q_{\theta }=0,$$

$$Q={\int }_0^{h}vdy=\Omega (t)h-\frac{1}{3}{h}^3\cos \theta$$

$$h(\theta ,0)={\text{h}}_0(\theta )$$

$$h_t+{\left(\Omega h-\frac{1}{3}{h}^3\cos \theta \right)}_{\theta }=0,$$

$$\Omega h-\frac{1}{3}{h}^3\cos \theta =Q,$$

$$\frac{dh}{d\tau }=-\frac{1}{3}{h}^3\sin \theta ,\frac{d\theta }{d\tau }=\Omega -h^2\cos \theta \equiv M(\theta ),\frac{dt}{d\tau }=1,$$

$$P(Q)={\int }_0^{2\pi }\frac{1}{\Omega -h^2\cos \theta }d\theta .$$

$$Q^{\ast }=\frac{2}{3}{\Omega }^{3/2}.$$

$$V={\int }_0^{2\pi }h(\theta )d\theta .$$

$$V^{\ast }=12{\Omega }^{1/2}{\int }_{{\hat{h}}_m}^1\frac{dx}{(x+2{)}^{1/2}(x^3+3x-2{)}^{1/2}},$$

$$g_t+(\Omega -h^2\cos \theta )g_{\theta }-2h{g}^2\cos \theta +2h^{2}g\sin \theta +\frac{1}{3}{h}^3\cos \theta =0.$$

$${\chi }_n^{\ast }={\int }_0^{2\pi }{\chi }_{n}d\theta$$

$$h_t+\frac{\partial }{\partial \theta }\left(\frac{\delta {\chi }_1^{\ast }}{\delta h}\right)=0,$$

$${\eta }_t+(M_s\eta {)}_{\theta }=0,$$

$$\psi (\theta )=\frac{1}{M_s}{\mathrm{e}}^{-\mathrm{i}\omega \nu (\theta )},\nu (\theta )\equiv {\int }_{\pi }^{\theta }\frac{d\xi }{M_s(\xi )}.$$

$$\frac{1}{M_s(\xi )}=\sum _{n=-\infty }^{\infty }{a}_n{\mathrm{e}}^{\text{i}n\xi },$$

$$a_0=\frac{1}{2\pi }P(Q),$$

$$\eta (\theta ,t)=\frac{1}{\Omega -h_s^2\sin \theta }U\left(t-{\int }_{{\theta }_0}^{\theta }\frac{d\xi }{M_s(\xi )}\right),$$

$$\mathbf{J}\frac{d\mathbf{u}}{dt}=\nabla H_p+\mathbf{p},\mathbf{p}=\left(\begin{array}{c}-b\cos \sigma t\\ 0\end{array}\right),\mathbf{J}=\left(\begin{array}{cc}0& -1\\ 1& 0\end{array}\right),$$

$$-H_p=h-\frac{1}{3}{h}^3\cos \theta .$$

$$\theta =-\frac{\pi }{2},h={\text{h}}^{\ast }$$

$$\mathcal{H}=\{h(\theta ):{\Omega }_{0}h-\frac{1}{3}{h}^3\cos \theta ={\mathcal{Q}}_0\text{mbox}forsome{\mathcal{Q}}_0\in [0,2/3]\},$$

$$h_t+\sigma h_T+{\left(\Omega (T)h-\frac{1}{3}{h}^3\cos \theta \right)}_{\theta }=0,$$

Full text

1 Introduction

Moffat [1] studied the dynamics of a thin viscous liquid film coating the outside of a horizontal circular cylinder that is rotating about its axis at a constant rate. This problem has important practical applications ranging from industrial processes (Ribatski & Jacobi [3]) to applications in the art world (Herczynski et al. [4]). In the absence of rotation the film will eventually drip under the action of gravity. Moffatt showed that rotation prevents dripping if the angular velocity of the cylinder exceeds a threshold value that depends on the kinematic viscosity of the liquid, the radius of the cylinder, and the acceleration due to gravity. This is in line with everyday experience: we can prevent honey from dripping off a spoon by rotating the spoon. However, it is tricky to rotate the spoon continuously at a constant rate, and in practice one tends rather to twist it back and forth, endowing the spoon with a time-dependent angular velocity.

Working on the basis of lubrication theory, Moffatt [1] derived a nonlinear model equation for the film thickness and showed that it has a steady solution if the aforementioned threshold criterion is met. Working around the same time, Pukhnachev [2] derived a more general version of the governing equation that incorporated the effect of surface tension, and also demonstrated the existence and uniqueness of a steady solution. The steady solution describes a film profile which is stationary in the laboratory frame. Assuming counterclockwise rotation, the film exhibits a bulge in thickness on the right side of the cylinder and is thinner on the left side (see figure 1). Moffatt [1] also described the results of experiments that revealed the importance of transverse instability, which manifests as a sequence of liquid lobes spaced out along the axis of the cylinder.

Numerous papers have followed examining various aspects of a problem that has proved to be very rich from a dynamics perspective. Hinch and Kelmanson [5] used asymptotic methods to show that surface perturbations decay and drift over a four time-scale cascade. Hinch, Kelmanson & Metcalfe [6] probed these results further, focusing on shock formation in the zero surface-tension case and providing an estimate of the shock formation time. Hansen & Kelmanson [7] used a boundary integral formulation to compute surface profiles under conditions of Stokes flow, allowing for films of arbitrary thickness and including surface tension. Peterson et al. [8] carried out a comprehensive linear stability analysis and revealed the parameter regimes in which steady-state solutions are stable. Duffy and Wilson [9] analysed both attached films and curtain flows (for which a film falls onto the cylinder from above, curves around, and falls off the bottom), developed analytical approximations, and identified critical flow transitions. Wray and Cimpeanu [10] used reduced-order modelling techniques to effectively capture the key features of the dynamics when inertia is included. Karabut [11] explored two distinct flow regimes depending on the angular velocity, which contributed to the classification of flow behaviours.

Evans, Schwartz and Roy [12, 13] presented both two- and three-dimensional models for coating flow on a rotating cylinder. They examined unsteady behaviours and instabilities. Gorla [14] focused on the rupture dynamics of non-Newtonian (power-law) films, addressing industrially relevant coating processes where film integrity is critical. Lopes et al. [15] introduced a new model equation which was derived using a gradient method based on Onsager’s variational principle and which includes the full expression for the surface curvature in the capillary stress term. These authors also compared the model predictions for steady flow with full computations of the Stokes equations. Kelmanson [16] extended the Moffatt–Pukhnachev model by including inertial effects, thereby increasing its relevance to moderate Reynolds number regimes.

Weak solutions exhibiting shocks were first introduced for the steady Moffatt-Pukhnachev flow by Johnson [17]; see also Badali et al. [18] and Benilov et al. [19]. Such solutions allow the cylinder to support a greater liquid volume than the smooth Moffatt solutions [1]. The stability of these weak solutions was studied by O’Brien [20] and Villegas-Díaz at al. [21], who demonstrated that stable configurations occur only when the shock is located in the fourth quadrant of the plane.

The problem of rimming flow, in which the liquid film coats the inside of the cylinder, is also of interest. As highlighted by Lopes [22] the thin-film equations for rimming flow and for the exterior flow problem are identical up to a certain order of approximation in the film thickness parameter, although key differences appear at higher order. Johnson [17] presented an analysis of steady-state coating flows inside rotating horizontal cylinders. O’Brien and Gath [23] identified the formation and position of shocks in rimming flows, addressing the occurrence of sharp transitions in film thickness. O’Brien [20] further contributed a linear stability analysis of rimming flows, describing conditions under which small disturbances may be amplified. Villegas-Diaz, Power and Riley [21, 24] examined the stability of rimming flows to two-dimensional perturbations, combining analytical and numerical techniques and exploring the impact of surface shear on flow stability.

The present paper is devoted to the study of a viscous liquid film that coats the outside of a horizontal cylinder which rotates at a constant rate onto which is superimposed oscillations of a certain amplitude and frequency. As in the related works discussed above, we herein assume that the thickness of the liquid film is everywhere much smaller than the radius of the cylinder and employ lubrication theory to derive a generalisation of the Moffatt-Pukhnachev equation, which incorporates a time-dependent modulation to the rotation rate. Although we discuss the subsequent dynamics in the context of the film coating the exterior of the cylinder, the aforementioned equivalence (up to some order in the film thickness parameter) between this and the rimming flow problem means that our observations carry over to the latter problem to some extent.

The paper is organised as follows. In section 2 we derive the thin-film equation that forms the basis of our model. In section 3, we briefly review the Moffatt-Puckhnachev flow for constant rotation, taking a dynamical system perspective. In section 4, the dynamical system formed from the characteristic equations of our model equation is studied and discussed. In section 5 we present solutions to the model equation and study the asymptotic limits of high frequency and low frequency oscillations. Finally, in section 6 we summarise our results.

2 Problem statement

We consider the flow of a liquid film of viscosity $\mu$ and density $\rho$ that coats the exterior of a circular cylinder of radius $a$, as is illustrated in figure 1. The motion in the liquid is driven by the downwards force of gravity, which acts in the negative $y$ direction, and by torsional rotations of the cylinder whose angular velocity is a prescribed periodic function of time, $t$. The flow is assumed to be two-dimensional in the $xy$-plane of the cylinder cross-section. The dynamics are described with reference to plane polar coordinates $(r,\theta)$ centred at the cylinder axis with $\theta=0$ aligned with the horizontal. Using thin-film theory Pukhnachev [2] and Moffatt [1] derived an equation for the film thickness when the cylinder rotates at a constant speed. Our first goal is to derive a modified version of this equation which accounts for a general angular velocity of the cylinder.

We define the dimensionless thin-film coordinate, $\hat{y}$, which measures radial distance from the surface of the cylinder, such that, with reference to figure 1,

[TABLE]

Here we have introduced the dimensionless group $G=\rho ga/\mu\Omega_{*}$, where $g$ is the acceleration due to gravity and $\Omega_{*}$ is a reference angular velocity, and the film thickness parameter $\epsilon=\bar{h}/a$, where $\bar{h}$ is the mean film thickness. Henceforth we shall assume that $\epsilon\ll 1$ so that the average film thickness $\bar{h}$ is small in comparison with the cylinder radius $a$.

The velocity $(u,v)$ in the $(r,\theta)$ directions, the pressure $p$, and time $t$ are non-dimensionalised by making the replacements

[TABLE]

To leading order in $\epsilon$ the dimensionless Navier-Stokes equations are

[TABLE]

The inertia terms have been neglected in these equations. This is justified provided that $\epsilon^{2}Re\ll G$, where the Reynolds number $Re=\rho\Omega_{*}\bar{h}^{2}/\mu$.

At the film surface, $y=h(\theta,t)$, we impose the kinematic condition,

[TABLE]

as well as the normal and tangential stress conditions

[TABLE]

respectively. The pressure in the air outside of the film has been taken to be zero. Moreover the contributions of both the viscous normal stress and the capillary stress in the first condition in (5) have been neglected. These assumptions are justified provided that $\epsilon^{2}G^{-1}\ll 1$ and $\epsilon G^{-3/2}C^{-1}\ll 1$, respectively, where $C=\mu\Omega_{*}a/\gamma$ is the capillary number with $\gamma$ the coefficient of surface tension. The boundary condition on the cylinder, $y=0$, is

[TABLE]

where $\Omega(t)$ is assumed to take the form

[TABLE]

for given constants $b>0$ and $\sigma>0$. If $b=0$ then the cylinder is rotating at a constant rate; this is the case originally studied by Moffatt and Pukhnachev and we shall henceforth refer to it as the MP problem.

It would appear, then, that there are two relevant time scales in the problem: first, there is the time taken for a fluid particle to complete one circuit of the cylinder under steady rotation, i.e. for the MP problem which has a steady flow solution (see [1]), and, second, there is the time scale associated with the modulational frequency $\sigma$.

Integrating the governing equations (3), applying the boundary conditions (5)-(6), and inserting the resulting expressions into the kinematic condition (4), we obtain the evolution equation

[TABLE]

where

[TABLE]

is the dimensionless flux in the film. The initial condition is

[TABLE]

for some appropriate choice of the function $\textup{h}_{0}$.

In summary we aim to solve (8) with initial condition (10) to determine the surface profile $h(\theta,t)$ for different choices of the parameters $b$ and $\sigma$. Before doing this, however, it is instructive to recall the salient details of the steady flow for the MP problem. We do this in the next section from a dynamical systems perspective which provides a novel and intuitive way of visualising the solution space.

3 The steady Moffatt-Pukhnachev flow

The steady MP problem is recovered by setting $b=0$ so that the cylinder is rotating at a constant rate. Although it is strictly equal to unity in this case, we find it convenient to retain $\Omega$ in the relevant equations to facilitate later discussion. The governing equation (8) is

[TABLE]

with initial condition given by (10). Moffatt [1] showed that, if a certain criterion is met, there exists a steady solution describing a fully attached film. Integrating the steady version of (11) with respect to $\theta$, we obtain

[TABLE]

where $Q$ coincides with its definition in (9) and is herein constant.

The constant $Q$ level curves for (12) are shown in the phase portrait in figure 2. Critical to note in this figure is the separatrix curve, shown with a solid gold line, which divides level curves that correspond to physical, fully- attached solutions (those below the solid gold curve) from those that correspond to unphysical solutions (those above it) for which $h$ blows up; in the latter case the blow-up has the local form $h\sim a|\theta_{c}-\theta|^{-1/2}$ for constant $a$, where $\theta_{c}=-\pi/2$ or $\pi/2$. Solutions with shocks, which include sections of level curves above and below the solid separatrix, are also possible and have been discussed by previous workers (e.g. Johnson[17], O’Brien & Gath [23]). We will touch upon these later for the case $b>0$.

A parametric description of the level curves is obtained by solving the characteristic equations for equation (11), namely

[TABLE]

where the independent variable $\tau$ varies continuously along a characteristic. It is straightforward to check that (12) is a first integral of (13). We define the period $P(Q)$ of a steady orbit to be the time taken for $\theta$ to change by $2\pi$ radians. By integrating the second equation in (13), we find

[TABLE]

This fixes a base frequency associated with the steady flow, given by $\omega^{\ast}(Q)=2\pi/P(Q)$. The dependence of $\omega^{\ast}$ on $Q$ is graphed in figure 3(b) where it can be seen that $\omega^{\ast}\in[0,1]$ and it is monotone decreasing in $Q$. We note the limits $\omega^{\ast}\to 1$ ($P\to 2\pi$) as $Q\to 0$, and $\omega^{\ast}\to 0$ ($P\rightarrow\infty$) as $Q\to Q^{*}$, that is as the separatrix in Figure 2 is approached.

Viewed as a two-dimensional dynamical system (13) has a saddle point at $(h,\theta)=(1,\,0)$. Its stable and unstable manifolds are the level curves of (12) with $Q=Q^{*}$, where

[TABLE]

They are shown with solid gold and dashed gold lines in figure 2. An initial condition for (13) at $\tau=0$ which corresponds to a point $(h,\theta$) lying beneath the solid separatrix will trace out a physically acceptable, fully attached solution for some $Q<Q^{*}$. The requirement that $Q\leq Q^{*}$ coincides with the criterion given by Moffatt [1].

In Figure 3(a) we show the variation with the dimensionless flux $Q$ of the dimensionless fluid volume,

[TABLE]

The plot assumes a smooth solution with no shocks. The maximum volume, attained when $Q=Q^{*}$, is given by

[TABLE]

where $\hat{h}_{m}$ is the minimum value of $h/\Omega^{1/2}$, which occurs at $\theta=\pm\pi$. In the case of $\Omega=1$, $\hat{h}_{\mathrm{min}}=p^{1/3}-1/p^{1/3}\approx 0.5961$, where $p=1+\sqrt{2}$. The integral (17) can be calculated exactly but the resulting expression in terms of elliptic integrals is unwieldy and is not included here. Numerically we calculate $V^{*}/(2\pi\Omega^{1/2})\approx 0.70708140$, which agrees with the value quoted by O’Brien & Gath [23].

If at $t=0$ the starting profile $\textup{h}_{0}(\theta)$ in the initial condition (10) coincides with one of the phase plane trajectories with $Q\leq Q^{*}$ in figure 2, $H(\theta)$, say, then $h(\theta,t)=H(\theta)$ for all $t\geq 0$. Consider instead a starting profile $\textup{h}_{0}$ that does not coincide with a steady solution such as that shown with a dashed red line in figure 2. If all or part of the profile crosses the solid gold separatrix then $h(\theta,t)$ will blow up at a finite-time singularity. Assume that $\textup{h}_{0}(\theta)$ is everywhere underneath the solid gold separatrix. Then $h(\theta,t)$ will be confined between the two osculating trajectories which are just tangent to the maximum and minimum of $\textup{h}_{0}$ (these trajectories are shown with heavier solid blue lines in figure 2). This is clear since each point on the initial profile $\textup{h}_{0}$ must traverse one of the level curves sandwiched between the two osculating curves. During the subsequent motion, the film slope $g\equiv h_{\theta}$ satisfies the equation obtained by differentiating (11)

[TABLE]

By the argument given above, $h$ is bounded, and the $g^{2}$ term is expected to drive the solution toward a finite-time slope singularity in a region where $\cos\theta>0$ (a similar assertion was made by Moffatt [1] via a slightly different argument). This infinite slope singularity heralds the onset of film overturning and shock formation, a phenomenon captured by the analysis of Hinch, Kelmanson & Metcalfe [6].

For the time-dependent problem it is worth recording that, with $\Omega$ constant, equation (8) possesses an infinite set of conserved quantities. We fix $\Omega=1$ and define, for integer $n\geq 0$, $\chi_{n}=\int Q^{n}dh$, where $Q$ is defined in (9) and where we treat $h$ and $\theta$ as being independent in the integration. Then

[TABLE]

is a conserved quantity for (8). To see this, differentiate $\chi_{n}^{*}$ with respect to $t$ and use the fact that $h_{t}=-Q_{\theta}$ and the $2\pi$-periodicity of $Q$ in $\theta$. The case $n=0$ corresponds to volume conservation, but $n\geq 1$ do not have an obvious physical interpretation. We also note that with $\Omega$ constant (8) can be put into the Hamiltonian form

[TABLE]

where $-\chi_{1}^{*}$ plays the role of the Hamiltonian.

3.1 Stability of the steady solution

The linear stability of the steady solutions described above for constant $\Omega$ has been discussed by O’Brien [20] for rimming flow, and the same analysis carries over here. We review briefly the essential details as these will prove useful in the ensuing analysis.

Writing $h=h_{s}(\theta)+\eta(\theta,t)$, where $h_{s}(\theta)$ is a steady solution of the MP problem for $Q<Q^{*}$, and $\eta(\theta,t)$ is a small perturbation, we substitute into (11). Neglecting higher order terms,

[TABLE]

where $M_{s}(\theta)=\Omega-h_{s}^{2}\cos\theta$ and, we emphasise, $\Omega$ is constant. Assuming that $Q<Q^{*}$ it is straightforward to show that $M>0$ (see Appendix A).

Since the coefficients in (21) are $2\pi$-periodic in $\theta$, we can use Floquet theory to justify writing $\eta=\textup{e}^{i\omega t}f(\theta)+\mbox{c.c.}$, where c.c. means complex conjugate, and where the $2\pi$-periodic in $\theta$ function $f(\theta)$ and the constant $\omega$ are to be found. Substituting into (21), and integrating, we find $f(\theta)=c\psi(\theta)$ for arbitrary constant $c$, where

[TABLE]

The integrand in (22) can be expressed as the Fourier series

[TABLE]

with $a_{n}=\overline{a}_{-n}$. The zeroth mode has

[TABLE]

where $P(Q)$ was defined in (14). It is clear from (22) that the required periodicity of $f$ is assured only if $\omega a_{0}=m$ for $m\in\mathbb{N}\cup\{0\}$. Since $a_{0}$ is real it follows that $\omega$ is real and the MP solution is neutrally stable (O’Brien [20]). This neutral stability was also previously noted by Villegas-Díaz et al. [21], who used the method of characteristics to obtain the solution of (21),

[TABLE]

where $\theta_{0}$ is an arbitrary constant. The function $U$ is set by the disturbance to the steady surface profile at $t=0$. Villegas-Díaz et al. [21] also discussed the stability of the steady MP solution for the special case $Q=Q^{*}$.

The fact that linear perturbations are stable, taken together with the observation from the previous section that any perturbation from a steady solution will lead to a finite-time slope singularity, makes clear that nonlinearity plays an important role in the dynamics, even for arbitrarily small perturbations. In passing it is interesting to note that an extended form of (11) that incorporates higher order terms in the lubrication approximation was derived by Benilov et al. [25] and was also shown to have neutrally stable steady solutions. Despite the neutrality of its eigenmodes, Benilov et al. [25] showed that the linearisation of the extended equation about a steady state admits a so-called ‘explosive instability’: despite the linearised problem yielding an infinite number of bounded harmonic modes (which would normally be taken to imply stability) it supports explosive disturbances that blow up in finite time.

4 The characteristic dynamical system

In general a numerical approach is required to handle the case of a time-dependent angular velocity, $\Omega=\Omega(t)$. In this section we carry out a numerical investigation of the dynamical system (13) with $\Omega(t)$ given by (7). We leave any interpretation of the results in the context of the cylinder flow problem until later sections.

The problem may be put into the form of the time-periodically perturbed Hamiltonian system

[TABLE]

where $\bm{u}=(h,\theta)^{T}$ and $\nabla=(\partial_{h},\,\partial_{\theta})^{T}$. The Hamiltonian, $H_{p}$, is such that

[TABLE]

The form of $\Omega(t)$ is given in (7). When $b=0$ the perturbation vanishes, $\bm{p}=\bm{0}$, and the system (26) is integrable. This corresponds to the steady flow MP problem discussed in section 3.

We integrate (26) numerically using the Matlab routine ode45. The initial condition is taken to be

[TABLE]

at $t=\tau=0$. We build up a behavioural map of the solution space distinguishing between solutions that remain bounded and those that blow up in finite time. The latter is detected by testing when $h$ exceeds a selected value. This is sufficient to give an accurate picture: blow-up is initiated when the solution trajectory of (26) latches onto the unstable manifold (leftmost dashed part of the gold separatrix in figure 2) and it occurs very rapidly thereafter and on a much shorter timescale than the time period of the cylinder oscillations. In our computations we deemed blow-up to have occurred when $h\geq 5.0$.

Figure 4 shows the behavioural map in the $b\sigma$-plane for $\textup{h}^{*}=1/\sqrt{3}$. The blue shading indicates the blow-up time with darker blue corresponding to later blow-up times. In the white region no blow-up was encountered in $0\leq t\leq 3500$, and we take this to mean that the solution remains regular and bounded. Since the chosen $\textup{h}^{*}$ is below $Q^{*}=2/3$, the solution at $b=0$ traverses one of the level curves below the solid gold separatrix in figure 2 and is therefore bounded. The complexity of the map is apparent and some of the features are reminiscent of the complex structures created by discrete one-dimensional maps, including highly intricate boundaries and the emergence of apparent self-similarity upon zooming in to certain parts of the picture. The latter phenomenon is seen in the various subsidiary panels in figure 4.

The red marker point at $(\sigma,b)=(0,0.317)$ indicates the threshold for blow-up predicted by the small-$\sigma$ analysis to be discussed in section 5.2. When $\sigma$ is small the solution is found to be quasiperiodic in nature. This is illustrated in figure 5 for the point $(b,\sigma)=(0.25,0.01)$ in figure 4. The integration was carried out up to $t=5000.0$. Panel (a) shows the time signal over the last $1000$ time units, and panel (b) shows a return map with $(h^{m}_{i},h^{m}_{i+1})$, where $h^{m}_{i}$ is the $i$th local maximum of the time signal, shown with dots. The appearance of an almost complete closed loop in the return map is the classic hallmark of quasiperiodic dynamics (e.g. Guckenheimer & Holmes [26]).

A particularly interesting feature of the map is that there appears to be some sort of resonance manifesting as a sequence of sharp protrusions extending to the left. With $\textup{h}^{*}=1/\sqrt{3}$ the unforced oscillator (viz (26) with $\bm{p}=\bm{0}$) has the natural frequency $\omega^{*}(1/\sqrt{3})=0.877$. The uppermost protrusion in figure 4 has its apex at a value of the forcing frequency, $\sigma$, that is close to this. Moreover, the various protrusions below occur at values of $\sigma$ that are close to rational multiples of $\omega^{*}$. The fact that these apparently resonant values of $\sigma$ do not quite coincide with $\omega^{*}$ (or a rational multiple thereof) is presumably attributable to the nonlinearity of the underlying oscillator.

5 The Moffatt-Pukhnachev flow with periodic modulation

Since the original partial differential equation, (8), for the cylinder flow problem is hyperbolic, the appearance of shocks is expected. Therefore, while indicative, the results of the previous section, which followed single-trajectory solutions of the Hamiltonian system (26), should be interpreted with some care in the context of the rotating cylinder problem. In this section we study solutions to the initial value problem (8), (10) with $\Omega(t)$ given in (7). It will be important to distinguish between initial conditions that correspond to solutions of the equivalent instantaneous steady problem and those that do not. To this end we define the class of functions

[TABLE]

where $\Omega_{0}=\Omega(0)=1+b$. If we write $h=\Omega_{0}^{1/2}H$ and $\mathcal{Q}_{0}=q_{0}\Omega_{0}^{3/2}$ then the restriction in $\mathscr{H}$ reduces to $H-\frac{1}{3}H^{3}\cos\theta=q_{0}$, contours of which correspond to solutions of the steady MP problem and are shown in figure 2.

Intuitively we might expect the $b\neq 0$ oscillations to cause the film surface to overturn and become multi-valued signifying the breakdown of (8). This is indeed the case, and for a general choice of modulation frequency $\sigma$ and amplitude $b$ the solution $h(\theta,t)$ develops an infinite slope singularity at some time $t^{*}>0$. Sample profiles are shown in figure 6 for the case $\mathcal{Q}_{0}=\sqrt{2}/3$ (black curve) taking $b=0.25$ and $\sigma=0.5$. The numerical computations were carried out by integrating the characteristic system (26) forward in time from the starting profile using the Matlab routine ode45. A parametric representation is introduced in which $(\theta(\xi_{i},t),h(\xi_{i},t))$ is tracked in time for $i=1,\ldots,N$, where $\xi_{i}$ is one of $N$ equally-spaced points in the range $[0,2\pi)$. The onset of overturning is detected by monitoring if there is a sign change in $\theta_{\xi}$, the differentiation in $\xi$ being done with spectral accuracy using a FFT. Typically we found that using a grid with $128$ equally-spaced points in $\xi$ is sufficient to get an accurate solution; the Matlab integrator ode45 uses an adaptive time step that, in the results to be presented, typically varies between $10^{-3}$ and $10^{-5}$.

In figure 7(a-c) we show how the overturning time $t^{*}$ varies with $\sigma$ for two cases chosen so that $\textup{h}_{0}\in\mathscr{H}$ in one and $\textup{h}_{0}\notin\mathscr{H}$ in the other. The modulation amplitude is set at $b=0.3$. The asymptotic results shown in this figure will be discussed in a later section. It is quite striking, particularly in figure 7(a), that the $t^{*}(\sigma)$ curve exhibits numerous jump discontinuities. The presence of these discontinuities has been confirmed by meticulous numerical computations, including a very careful convergence study using up to $N=1024$ points in $\xi$ and a time step of size $10^{-5}$.

To explain the jump discontinuities, and with reference to the system (26), we regard the Hamiltonian $H_{p}(\theta,h)$ given in (27) as the stream function for a two-dimensional incompressible fluid flow in the $\theta h$-plane. In this interpretation the flow occupies the whole plane and is not confined to the domain of the liquid for the cylinder problem. Solving (26) elicits the trajectories of individual particles advected passively within this flow, with the liquid surface corresponding to a material line. As it moves with the flow, this material line is stretched and distorted by the flow’s strain and vorticity fields. These fields, which are independent of time, are illustrated in figure 6(b): the arrows show the direction of dominant strain corresponding to the eigenvector of the rate of strain tensor with positive eigenvalue (we recall that since the rate of strain tensor is symmetric its two eigenvalues, which are real, sum to zero, and its two eigenvectors are mutually orthogonal). In general, material elements tend to align with the dominant strain eigenvector whilst being rotated one way or another according to the sign of the vorticity. The constant vorticity contours are solid where the vorticity is positive, promoting counterclockwise rotation, and dashed where the vorticity is negative, promoting clockwise rotation. (For $b=0$ a steady solution has a surface profile that stays fixed as the competitive effects of strain and rotation are in perfect balance; see the solid red line in figure 6b.)

Certain parts of the strain-vorticity field promote steepening of the material surface and others promote flattening. The region around $\theta=-\pi/2$ presents a particular danger zone for surface steepening. Here the magnitude of the vorticity is small and the dominant strain eigenvectors are almost vertical. In contrast, the region around $\theta=\pi/2$ has weak vorticity and almost horizontal strain vectors, and therefore this region strongly encourages flattening of the surface. It appears that we get a jump discontinuity in $t^{*}$ at a certain $\sigma$ because, just below this frequency, the surface overturns near to $\theta=-\pi/2$; but for a slightly larger $\sigma$ the overturning is just avoided and the surface must travel a further distance through a less dangerous region before overturning is finally induced in the danger zone some time units later.

In figure 7(d) we show the overturning time $t^{*}$ for the case $b=0.5$ and $\sigma=0.6$ and for a range of different initial conditions. The latter are given by (10) with $\textup{h}_{0}\in\mathscr{H}$ and $\mathcal{Q}_{0}=q_{0}\Omega_{0}^{3/2}$. The panel shows $t^{*}$ plotted against $q_{0}$, and we see the occurrence of jump discontinuities at certain values of $q_{0}$. Some of these are highlighted in the inset. Since $q_{0}\to 0$ corresponds to taking an initial profile that is close to the wall and almost flat (see figure 2), we expect the overturning time to diverge in this limit. The dashed red line in figure 7(d) suggests that this is indeed what occurs and, moreover, it happens such that $t^{*}\sim q_{0}^{-6}$.

To further understand the dynamics it is instructive to consider the flow in the high frequency ($\sigma\gg 1$) and low frequency ($\sigma\ll 1$) limits. These are examined in the following subsections.

5.1 High frequency limit ($\sigma\gg 1$)

Our numerical results suggest that the flow is periodic in time if $\sigma$ is sufficiently large, and this motivates an investigation of the dynamics when $\sigma\gg 1$. In this limit there are two naturally disparate time scales in the problem: an $O(1)$ time scale associated with the steady part of $\Omega$, that is the continuous rotation, and an $O(1/\sigma)$ timescale associated with the rapid superimposed oscillations. Having said this, we might intuit that for an initial condition corresponding to a steady Moffatt-Pukhnachev solution, the former timescale is essentially removed and the flow will develop on timescales of the rapid oscillation.

Keeping the preceding comments in mind, we perform a multiple-scales analysis incorporating both time scales. Let $h=h(\theta,t,T)$ with $T=\sigma t$ acting as the fast time variable and $t$ acting as the slow time variable. Assuming in the usual way that $t$ and $T$ are independent (8) becomes

[TABLE]

where $\Omega(T)=1+b\cos T$. We expand by writing

[TABLE]

Introducing this expansion into (30), at $O(\sigma)$ we find that $h_{0T}=0$, which implies $h_{0}=h_{0}(\theta,t)$. At order $O(1)$ we obtain $h_{1T}=-bh_{0\theta}\cos T-F(\theta,t)$, where

[TABLE]

Integrating with respect to $T$,

[TABLE]

where $A_{1}(\theta,t)$ is an arbitrary function of integration. In order to avoid secular terms, we demand that $F=0$. Then the leading order term $h_{0}$ satisfies the original equation (8) with constant unit forcing frequency, $\Omega=1$, that is it corresponds to a solution of the constant rotation rate Moffatt-Pukhnachev problem. To prevent $h_{0}(\theta)$ developing a slope singularity (see section 3), we take $h_{0}\equiv h_{s}(\theta)$, where the steady solution $h_{s}$ satisfies (12) with $\Omega=1$ for some flux $Q<2/3$.

Proceeding, we have

[TABLE]

At order $O(\sigma^{-1})$,

[TABLE]

Substituting (34) into the right hand side of (35), we see that secular terms in $h_{2}$ will not arise if the $T$ independent terms vanish, that is if

[TABLE]

where $\mathcal{M}=1-h_{0}^{2}\cos\theta$. This is the same as the equation that governs the linear stability of the MP problem, namely (21). Since $h_{0}$ has been chosen to be a steady solution of the MP problem, according to the results of section 3.1, $A_{1}$ is $t$-periodic with frequency $\omega=m/a_{0}$, where $a_{0}$ is given by (24) ($\Omega=1$) and $m\in\mathbb{N}\cup\{0\}$. Hence we have

[TABLE]

where $\lambda$ is an arbitrary constant. The $2\pi$-periodic function $\psi(\theta)$ was given in (22).

Integrating (35),

[TABLE]

where $A_{2}(\theta,t)$ is an arbitrary function of integration. At order $O(\sigma^{-2})$

[TABLE]

The secularity condition requiring that the $T$-independent terms on the right hand side of (39) be eliminated takes the form

[TABLE]

where

[TABLE]

We seek a solution in the form $A_{2}(\theta,t)=\alpha_{20}(\theta)+\mathrm{e}^{2\mathrm{i}\omega t}\alpha_{22}(\theta)+\mbox{c. c.}$, with $\alpha_{20}$ and $\alpha_{22}$ required to be $2\pi$-periodic in $\theta$. It is clear from (40) that $\alpha_{20}(\theta)$ has this property. Solving for $\alpha_{22}(\theta)$ using an integrating factor, we find that it is $2\pi$-periodic if

[TABLE]

This holds if $\omega=0$ (so that $\psi\propto 1/\mathcal{M}$) in which case $A_{1}=\kappa_{1}/\mathcal{M}(\theta)$ and $A_{2}=(S(\theta)+\kappa_{2})/\mathcal{M}(\theta)$ for arbitrary constants $\kappa_{1}$ and $\kappa_{2}$. Thus the slow timescale $t$-dependence drops out to the current order of approximation, as was anticipated.

In the light of the preceding discussion, and recognising that the first term on the right hand side of (34) represents a translation in $\theta$ by an amount $\sigma^{-1}\Omega_{T}$, we may now revise the expansion (31) to read

[TABLE]

This represents a time-periodic asymptotic solution valid when $\sigma\gg 1$. Suppose that we compute a numerical solution to the governing equation (8) for a particular choice of parameter values and taking $\textup{h}_{0}(\theta)$ in the initial condition (10) to be a steady solution $h_{s}(\theta)$. Since $\Omega_{T}(T=0)=0$, we can fit (43) to the initial profile at $t=T=0$ up to the second order in $\sigma$ by setting $\kappa_{1}=0$. However, at second order there is no way to choose the constant $\kappa_{2}$ such that $h_{2}(T=0)=0$. We therefore have a discrepancy between the numerical solution and the asymptotic approximation of $O(\sigma^{-2})$. In figure 8 we demonstrate agreement between the numerical computation and the asymptotic solution (43) with $\kappa_{1}=0$. Here the numerical solution was obtained by integrating the characteristic system (26) forward in time from the starting profile $h_{s}(\theta)$, where $h_{s}(\theta)$ solves (12) with $\Omega=1$ and $Q=0.5$. The vertical axis shows the norm $h_{N}$ defined to be the maximum value of $h-h_{s}(\theta+\sigma^{-1}\Omega_{T})$ over $\theta\in[0,2\pi)$ during the time integration period $t\in[0,10/\sigma]$. A best fit curve $h_{N}=0.1/\sigma^{2}$ is shown with a dashed red line.

The above discussion suggests that, for large $\sigma$ at least, for a sufficiently carefully chosen initial condition it should be possible to obtain a time-periodic solution to (8) which does not overturn. However, overturning is expected for a general initial condition and we can use our theory to estimate the overturning time, $t^{*}$. In general the starting profile $\textup{h}_{0}(\theta)$ in the initial condition (10) will not exactly meet the requirement that $F\equiv 0$, and so we expect transient growth in $T$ according to (33). If we identify $\textup{h}_{0}$ with $h_{s}$, where $h_{s}$ solves (12) with $\Omega=1$, it follows from (32) that $F=O(\sigma^{-2})$, and the asymptotic theory predicts transient growth so that the uniformity of the expansion (31) is destroyed when $t=O(\sigma^{2})$. If instead we identify $\textup{h}_{0}$ with a profile that does not solve (12) with $\Omega=1$, then $F=O(1)$ and the expansion (31) fails when $t=O(1)$. We interpret the failure of the expansion as the signature of overturning. This viewpoint is supported by our numerical solutions to the full governing equation (8). For the overturning times reported in figure 7 we see that in panel (a), for which $F=O(1)$, $t^{*}$ approaches a constant as $\sigma\to\infty$ in agreement with the preceding remarks. In panel (b), for which $F=O(\sigma^{-2})$, $t^{*}$ grows apparently like $\sigma^{2}$. Evidently the onset of overturning can be considerably delayed in the high frequency limit by a judicious choice of the initial condition.

5.2 Low frequency limit ($\sigma\ll 1$)

The numerical results shown in figure 4 suggest that when $\sigma\ll 1$ blow-up occurs when $b$ exceeds a threshold value. This motivates an analysis in the low frequency limit. As in the previous subsection, we follow a multiple scales approach, in this case with $t$ as the fast time scale and $T=\sigma t$ as the slow time scale.

It is convenient at the outset to rescale the film thickness, writing $h(\theta,t,T)=\Omega^{1/2}H(\theta,t,T)$, where $\Omega(T)=1+b\cos T$. We assume $|b|<1$ so that $\Omega>0$ for all $T$. We then posit the expansion

[TABLE]

assuming $\sigma\ll 1$. Inserting into (30) we obtain at leading order, $O(1)$,

[TABLE]

Following the discussion in section 3, for a general initial condition we expect that $H_{0}$ will reach a slope singularity at a finite time $t$. However the solution will remain regular if at $t=0$ the profile $H_{0}$ coincides with a $t$-independent solution of (45), namely one that satisfies

[TABLE]

for some $q(T)$ to be determined later. We therefore insist that $H_{0}(\theta,T)$ solves (46). We note that this is essentially the same as the cubic equation (12) for the steady MP problem, with $T$ playing the role of a parameter, and with the relationship $Q=\Omega^{3/2}q$. A question of interest, then, is whether $q(T)$ can reach the threshold value of $2/3$ identified in section 3, and so drive the leading order solution toward blow-up.

At first order, $O(\sigma)$,

[TABLE]

where $M_{0}(\theta,T)=1-H_{0}^{2}\cos\theta$. This essentially presents a forced version of the linear stability equation (21). The solution is

[TABLE]

where the function $U_{1}$ is arbitrary, $z=\Omega(T)t-\nu_{0}(\theta)$, and

[TABLE]

We require $H_{1}$ to be $2\pi$-periodic in $\theta$. To check this we first differentiate (46) with respect to $T$, and rearrange to get

[TABLE]

Integrating this with respect to $\theta$,

[TABLE]

With reference to (23) and (24) we notice that $\nu_{0}(\theta)$ is composed of $a_{0}\theta$ plus a periodic function of $\theta$. Rearranging (46), using the definition of $M_{0}$ and integrating with respect to $\theta$, we find

[TABLE]

Using (51) and (52) we may write

[TABLE]

If $q<2/3$ then $H_{0}$ satisfying (46) is bounded and $2\pi$-periodic in $\theta$. The right hand side of (53), and hence $H_{1}$, has the same properties if

[TABLE]

where

[TABLE]

The solution for $H_{1}$ is then given by

[TABLE]

It is clear from (56) that $H_{1\theta}$ is bounded. Furthermore, $H_{1T}$ is bounded for all $T$ provided that $|q_{T}|<\infty$ by virtue of (50).

The solution to the problem at $O(\sigma^{n})$ has the form

[TABLE]

with $U_{n}$ arbitrary, and $R_{n}=-(\Omega^{1/2}H_{(n-1)})_{T}+(F_{n}\cos\theta)_{\theta}/3$, where the $F_{n}(H_{0},H_{1},\ldots,H_{n-1})$ are known (for example, $F_{2}=3H_{0}H_{1}^{2}$). Let us assume that the $U_{n}$ are given for $n\geq 0$, for example by setting an appropriate initial condition at $t=T=0$. Since, if $|q_{T}|<\infty$, $H_{n}$ and $H_{n\theta}$ are both bounded for $n=0$, $1$, it is clear that the same property holds for $n\geq 2$. Given these remarks we expect the expansion (44) to remain uniform as time increases.

The nonlinear ordinary differential equation (54) determines $q(T)$. It has the invariant $\Omega^{1/2}\int_{0}^{2\pi}H_{0}d\theta$, which means that fluid volume is conserved at leading order (this can be shown by differentiating (46) with respect to $T$ and then by straightforward manipulations). Furthermore we can show that when $2/3-q\ll 1$

[TABLE]

Inserting this result into (54) we can then see that $q_{T}\to 0$ as $q\to 2/3$ and, moreover, $q_{TT}\to-\infty$ in the same limit. This behaviour is demonstrated numerically in figure 9(a,b) for $b=0.5$ and initial condition $q(0)=0.48$. Good agreement between the approximation (58), valid when $2/3-q\ll 1$, and the numerically computed $\mathcal{H}(q)$, is shown in figure 9(d).

When $q=2/3$ the leading order profile $H_{0}$ is described by the the heteroclinic connection along the separatrix that connects the saddle points at $\theta=0$, $2\pi$ (see figure 1). In this case $\nu_{0}(\theta)\to\infty$ as $\theta\to 0$. A slight rearrangement of (53) is

[TABLE]

With the help of (46) written for $q=2/3$ we can deduce that

[TABLE]

the inequality being clear given that the integral represents one half of the total fluid volume on the cylinder. Setting $q=2/3$ in (59), and taking the limit $\theta\to 0$, we observe that the terms in the large curved bracket on the right hand side of (59) must vanish. Therefore in this case (cf. (54) for $q<2/3$) we have

[TABLE]

and it follows that $q=2/3$ for all $T$. Unlike (54), this equation does not conserve the quantity $\Omega^{1/2}\int_{0}^{2\pi}H_{0}d\theta$ since in this case the integral has a fixed value and $\Omega$ varies with $T$.

Next we integrate (54) numerically starting from some $q(0)$. For $q(0)<q^{*}(b)$, where $q^{*}(b)$ is to be determined, $q$ is periodic in time $T$. For $q(0)=q^{\ast}(b)$, we find that $q(\pi)=2/3$ and $q_{T}(\pi)=0$. If $q(0)>q^{*}(b)$ then $q(T_{s})=2/3$ with $q_{T}(T_{s})>0$ at some $T_{s}=T_{s}(q_{0},b)>0$. This heralds the breakdown of the approximation since a smooth solution of (46) ceases to exist (see the discussion in section 3). Accordingly $q^{*}(b)$ descirbes a curve in the $(b,q_{0})$ plane which divides regular solutions of (54) from those which blow up. The implicit nonlinear relation for $q^{\ast}(b)$,

[TABLE]

follows on integrating (54) from $T=0$ to $T=\pi$. We solve this using Newton’s method to obtain the $q^{*}(b)$ curve shown with a black solid line in figure 10.

As an alternative to integrating (54) numerically, a simple approximation to $q(T)$ can be found by first noting that the one-term Taylor expansion about $q=0$, given by $\mathcal{H}\approx q$, holds good unless $q$ is very close to $2/3$ (this is confirmed by comparison with the numerical solution to (62)). With this approximation (54) simplifies to $q_{T}=-(\Omega_{T}/2\Omega)q$, which can be integrated exactly to yield

[TABLE]

where $q_{0}=q(0)$ and $\Omega_{0}=\Omega(0)=1+b$. This result immediately provides the approximation to the $q^{*}(b)$ curve shown with a solid red line in figure 10,

[TABLE]

In figure 11 we compare the approximation (63) with the solution of the full equation (54) for the case $b=0.5$ and for a number of different initial conditions, $q(0)$. Evidently the approximation performs very well with the greatest discrepancy occurring at the maxima for the larger values of $q(0)$.

It is instructive to compare the present small-$\sigma$ asymptotic results with numerical solutions to the full governing equation (8). First, we note that for $q_{0}<q_{*}(b)$ the dependence of the solution on both $t$ and $T$ indicates that the small-$\sigma$ solution will in general be quasiperiodic, and this is consistent with the comment made earlier in section 4 (see also figure 5). Figure 12 shows time snapshots of the surface profiles obtained by integrating the full governing equation (8) numerically for a variety of initial conditions (10) with $\textup{h}_{0}\in\mathscr{H}$. The leading order asymptotic solution, using the approximation (63) for $q(T)$, is included for comparison. The agreement between the two sets of results is generally excellent, but with some discrepancy at the surface maximum seen when $q$ gets close to $2/3$.

In figure 7(a) we showed the numerically computed overturning time $t^{*}$ for the case $b=0.5$ and $q_{0}=0.5$. These parameter values correspond to a point in the white regular region in figure 10. We can match the expansion (44) to the initial condition used in figure 7(a) by taking $H_{0}(\theta,0,0)=\textup{h}_{0}$, and $H_{n}(\theta,0,0)=0$ for $n\geq 1$ (the latter can be effected by choosing the arbitrary functions $U_{n}$ in (57) appropriately). Accordingly the asymptotics predict no overturning in the limit $\sigma\to 0$ and this is in line with the very rapid growth of $t^{*}$ observed in figure 7(a) in this limit. In contrast figure 7(b) shows the overturning time for the same value of $b$ but with the initial profile solving $\textup{h}_{0}-(1/3)\textup{h}_{0}^{3}\cos\theta=1/2$. Matching the asymptotic expansion (44) to this initial condition, at leading order we find that (46) is not satisfied at $T=0$ and, consequently, the solution of the leading order equation (45) will develop a slope singularity at a finite time, as was discussed in section 3. The predicted overturning time $t^{*}=242$ obtained by integrating (45) with $\Omega$ fixed at its initial value $\Omega(0)=1+b$, is shown with a dashed line in figure 7(b).

The blow-up that occurs at $T_{s}=T_{s}(q_{0},b)$ in the blue region in figure 10 (i.e. for $q^{*}(b)<q_{0}<2/3$) can be interpreted as the overturning of the film profile in the full problem. At $T=T_{s}$ we expect the expansion (44) to disorder, yielding the approximation to the overturning time for the full problem $t^{*}=T_{s}(q_{0},b)\sigma^{-1}$ valid for $\sigma\ll 1$. In figure 7(c) we showed the overturning time for the case $b=0.5$, $q_{0}=\sqrt{2}/3=0.471$ corresponding to a point in the blow-up region in figure 10. With this value of $q_{0}$, integrating (54) we compute $T_{s}=2.270$ and hence the prediction $t^{*}=2.270\sigma^{-1}$. This is shown with a red dashed line in figure 7(c). At sufficiently small $\sigma$ It is in excellent agreement with the numerically computed overturning time for the full governing equation.

Finally, we note that the present results are consistent with the blow-up map shown earlier in figure 4. For this figure, trajectories solving (26) were computed with the initial condition (28) and $\textup{h}^{*}=1/\sqrt{3}$. This means that, in the context of the present small $\sigma$ analysis, the relevant trajectories start on the $(\theta,H_{0})$ curve with $q_{0}=1/\sqrt{3(1+b)}$. This curves is shown with a red dashed line in figure  10, and where it intersects the black curve representing $q^{*}(b)$ gives the threshold $b$ value for small-$\sigma$ blow-up in figure 4. The intersection is found to occur at $b=0.317$. This is shown with a red marker in figure 4, where it can be seen to be in good agreement with the boundary curve between regular and blow-up behaviour at small $\sigma$.

5.2.1 Single shock solution for $q=2/3$

Under the small-$\sigma$ assumption, according to (61) if $q(0)=2/3$ then $q=2/3$ for all $T$. If we take the surface profile for $q=2/3$ to be $H_{0}$ satisfying (46) and given by the solid gold separatrix curve in figure 2, a contradiction arises since conservation of mass at leading order requires $\Omega^{1/2}\int_{0}^{2\pi}H_{0}d\theta$ and $\Omega$ is $T$-dependent. We can circumvent the mass conservation issue by following Johnson [17] to introduce a shock at a location $\theta^{\textup{s}}(T)>0$. The film profile follows the solid separatrix curve in figure 2 over $-\pi\leq\theta\leq\theta^{\textup{s}}$ (proceeding smoothly through the saddle at $\theta=0$ onto the dashed part of the separatrix curve where $H_{0}>1$) and then, at $\theta=\theta^{\textup{s}}$, it drops vertically down to the solid part of the separatrix below and continues along this up to $\theta=\pi$.

A shock solution allows the cylinder to support a larger liquid volume than a smooth solution. The liquid volume increases monotonically as $\theta_{s}$ varies from [math] for which $V^{*}/2\pi\approx 0.70708140$, up to $\theta_{s}=\pi/2$ for which value the film has infinite thickness at the shock but the finite volume

[TABLE]

where $(a_{1},b_{1})=(H_{m},1)$, $(a_{2},b_{2})=(1,\infty)$, $(a_{3},b_{3})=(H_{m},2/3)$. Here $H_{m}$ is the minimum of $H$ on the separatrix which occurs at $\theta=\pm\pi$. We calculate $V_{\text{max}}/(2\pi\Omega^{1/2})\approx 1.102317$, which agrees with the value quoted by Villegas-Diaz et al. [24].

Suppose that we start at $T=0$ with a leading order profile $H_{0}$ featuring a shock located at $\theta=\theta^{\textup{s}}(0)$ with $0<\theta^{\textup{s}}(0)\leq\pi/2$. The fluid volume $V(T)=\int_{0}^{2\pi}H_{0}d\theta$. As time increases the shock location $\theta^{s}(T)$ must adjust to ensure mass conservation, i.e. such that $\Omega^{1/2}V(T)=(1+b)^{1/2}V(0)$. The requirement that $|\theta^{\textup{s}}|<\pi/2$ imposes the constraint

[TABLE]

This places an upper limit on the modulation amplitude, $b$, for a given initial volume.

5.2.2 Double shock solution for $q<2/3$

Up to this point the discussion for $q(T)<2/3$ has assumed that the leading order term $H_{0}$ is such that $\Omega^{1/2}H_{0}\in\mathscr{H}$ at $t=T=0$ with $H_{0}$ a smooth function (i.e. no shocks). Johnson [17] demonstrated that one can construct double-shock solutions for which the film profile exhibits a bulge in fluid thickness enclosed by the two shocks. Referring to the phase portrait in figure 2, we see that a double-shock solution can be put together by following one of the blue trajectories at the bottom of the figure for a chosen $Q$, and then, somewhere in the interval $-\pi/2<\theta<\pi/2$, jumping up to and then back down from the equivalent $Q$ curve at the top of the figure, and then continuing along the original trajectory to complete the profile. It is not necessary that the jumps be located symmetrically about $\theta=0$.

A case of particular interest is that of zero flux, $Q=0$, indicated by the broken blue line in figure 2. For $Q=0$, and in the absence of shocks, the only static solution is a film of zero thickness. For the double-shock solution a solitary mass of fluid with compact support is carried around by the cylinder, with the fluid inside recirculating (Johnson [17] describes how to extend the thin-film model to include next order terms and allow for a boundary layer around each shock to facilitate this recirculation). The maximum scaled volume that can be supported in this case can be expressed by

[TABLE]

where $F$ is the incomplete elliptic integral of the first kind (e.g. Olver et al. [27]).

Assume for our small-$\sigma$ problem that a double shock is present in the leading order profile $H_{0}$. Let the jumps be located at $\theta=\theta_{L}<0$ and $\theta_{R}>0$ with $|\theta_{L,R}|\leq\pi/2$. Since $\Omega^{1/2}\int_{0}^{2\pi}H_{0}d\theta$ is automatically conserved as $q$ evolves according to (54), the shock locations $\theta_{L,R}$ do not depend on $T$ and so remain fixed. If $q(0)<q^{*}(b)$ then $q(T)$ is $T$-periodic and will oscillate in the range $(0,2/3)$ with the double shock profile adjusting accordingly. If $q(0)\geq q^{*}(b)$ then $q(T)\to 2/3$ as $T\to T_{1}$, for some finite $T_{1}$, with $q_{T}(T_{1})>0$; in this case, we jump to the single shock solution described in the previous subsection.

6 Conclusion

We have examined the dynamics of a thin viscous liquid film flowing on the outside of a horizontal cylinder which is rotating about its axis with an angular velocity that includes a steady part and a time-periodic part. Surface tension has been neglected. If the time-periodic component of the angular velocity is removed, the problem reduces to that studied by Moffatt [1] and Pukhnachev [2], and subsequent authors.

Our brief review of the constant rotation rate Moffatt-Pukhnachev problem included a novel perspective on describing steady solutions using a phase-plane analysis. This perspective makes plain visually the steady solution space, including both smooth solutions and those exhibiting shocks. It also makes clear graphically the threshold on the rotation rate for a steady solution to exist, and the onset of dripping that is expected if this threshold is exceeded. A separatrix in the relevant phase portrait delineates a boundary representing the extreme continuous steady solution which supports the maximum fluid volume (although larger volumes are possible for solutions with shocks).

The introduction of a time-periodic part into the angular velocity fundamentally alters the behaviour of the system. In this case it appears that in general it is not possible to choose an initial condition such that the slope of the film profile remains finite for all time. Rather it seems that the dynamics inevitably progress to a point where the film starts to overturn and a shock is formed. Numerical investigation of the characteristic dynamical system revealed a highly intricate and complex blow-up map showing the tendency toward film overturning or blow-up depending on the amplitude and frequency of the angular velocity modulation. Notably, sharp protrusion-like structures in these maps appear near to the base frequency $\omega^{*}$ of the Moffatt-Puchnachev system and, seemingly, at rational multiples thereof.

To develop a deeper understanding of the dynamics we also performed asymptotic analyses in the limits of high and low forcing frequency. In latter case a multiple-scales analysis suggested that the flow is time-periodic if the initial condition is chosen sufficiently carefully. For a general initial condition the film profile will overturn; but this overturning can be considerably delayed by choosing the initial condition to coincide with a Moffatt-Puchnachev steady profile. A multiple-scales approach was also used in the low-frequency limit. In this case a secularity condition was derived, which yielded a threshold modulation amplitude beyond which a singularity is reached at some time. This singularity, which heralds the breakdown of the multiple-scales expansion, was interpreted as indicating the blow-up of the solution to the full problem, and confirmation of this point was provided through comparison with numerical solutions of the full problem. Regular, bounded solutions exist below the threshold amplitude, and it was shown that these correspond to quasiperiodic solutions of the full problem. Single-shock and double- shock solutions for which the film jumps instantaneously in height, were also discussed. In both cases the abrupt change in height can be smoothed by introducing boundary layers as discussed by Johnson [17].

In summary our results have revealed a complex interplay between the amplitude and frequency of the torsional movement of the cylinder and the stability and regularity of the liquid film. As in Moffatt’s original work we have herein neglected the effect of surface tension. The influence of this force on the dynamics of the film in the presence of a modulated rotation rate is the subject of our ongoing research.

Appendix

We demonstrate that $M=\Omega-h_{s}^{2}\cos\theta$ is positive for all $\theta\in[0,2\pi)$, where $h_{s}(\theta)$ satisfies

[TABLE]

for constant $\Omega$, provided that $Q<(2/3)\Omega^{3/2}$. Writing $h_{s}=\Omega^{1/2}H$, $Q=\Omega^{3/2}q$, and manipulating (68), we have

[TABLE]

where $q<2/3$. The case $\cos\theta=0$ is clear since then $H=q$. Assuming $\cos\theta>0$ and writing $H=(\cos\theta)^{-1/2}R$, (68) becomes

[TABLE]

where $\lambda=(\cos\theta)^{1/2}q<2/3$ and $g(x)\equiv x-x^{3}/3$. The cubic equation (70) has three real roots, and it’s clear that the root of interest is that for which $0<R<1$ (refer to Figure 2; the other positive root has $R>1$ and will produce blow-up). Evidently $g(x)>2x/3$ for $0<x<1$ since over this range $g(x)-2x/3=x(1-x^{2})/3>0$. It follows that $\lambda>2R/3$ and, consequently, $H<3q/2$ so that $M>0$. A similar, but slightly simpler, argument applies in the case when $\cos\theta<0$.

Acknowledgments

This work was partially supported by the research project, PID2020-115961RB-C31 financed by MCIN/AEI/10.13039/501100011033. AJBL would like to thank the Spanish Ministry of Science, Innovation and Universities for the financial support provided by the Fellowship PRE2021-099112 that allowed his research stay in the School of Mathematics at the University of East Anglia.