Instabilities and bifurcations of liquid films flowing down a rotating fibre

Abstract We consider the dynamics of a gravity-driven flow coating a vertical fibre rotating about its axis. This flow exhibits rich dynamics including the formation of bead-like structures and different types of steady or oscillatory travelling waves driven by a Rayleigh–Plateau mechanism modified by the presence of gravity and rotation. Linear stability shows that the axisymmetric mode dominates the instability when the rotation is slow, which allows us to derive a two-dimensional model equation under the long-wave assumption. The spatio-temporal dynamics and nonlinear wave solutions are then investigated by the model equation. The spatio-temporal stability analysis showed that the absolute instability is enhanced by the rotation. Steady travelling-wave states and relative periodic states are observed in the numerical simulations of the model equation, which show that the rotation tends to suppress the formation of relative periodic states. To examine this, a linear stability analysis of steady travelling waves is performed, indicating that the rotation has a stabilizing effect on the steady travelling waves. This result is adverse to the destabilizing effect of rotation on the linear stability of initially uniform films. A bifurcation analysis shows that the relative periodic state is born from the instability of steady travelling wave, which represents the coalescence and breakup process between a large droplet and a serial of much smaller droplets.


Introduction
The coating of liquid film on a fibre has given rise to considerable scientific interest because of its relevance to many industrial applications, for example, draining, coating of insulation on a wire, and the protection of tube walls (Quéré 1999). The exterior coating flow driven by gravity down a vertical fibre is a typical unstable open flow in which formation of droplets, or beads, due to the well-known Rayleigh-Plateau mechanism (Rayleigh 1892) inevitably occurs.
Two different types of instabilities emerge in the exterior coating flows. One is the classical spatio-temporal disorder induced by the Kapitza instability mode (Kapitza & Kapitza 1949) of films falling down vertical planes, hereinafter referred to as 'K' mode. The wavy dynamics on a liquid film flow down an inclined or vertical plate has given rise to considerable interest up to date. For more related works to the dynamics of a falling film on inclined plate, we refer the readers to the monograph by Kalliadasis et al. (2012). The other type is the emergence of very regular drop-like wave patterns due to the Rayleigh-Plateau instability, hereinafter referred to as 'RP' mode (Kliakhandler, Davis & Bankoff 2001;Craster & Matar 2006). Quéré (1990) was the first to perform experiments to investigate the drop formation in the process of drawing wires from liquid baths. He found that for a thick film on a slender fibre, small undulations grow into large droplets by swallowing the other ones, and quickly fall, leaving behind them a thick film which breaks into droplets. However, for a thin film on a thick fibre, the process of drop formation may be arrested by the mean flow.
The problems of the dynamics of exterior coating flows on fibres have been extensively investigated using asymptotic models. Frenkel (1992) derived a thin film model using a lubrication approximation for the coating flow wherein the fibre radius a is assumed to be much larger than the film thickness h. Kalliadasis & Chang (1994) studied the mechanism for drop formation of the problem using Frenkel's equation, and found that when the thickness of the coating film exceeds a critical thickness h c , the interface can evolve into highly localized drops due to the Rayleigh-Plateau mechanism. Kliakhandler et al. (2001) performed experiments to examine the drop formation process in the case where the film thicknesses are of the order of the fibre radii. Instead of using Frenkel's thin film model, a thick film model was proposed by Kliakhandler et al. (2001). Note that the thick film model is not asymptotically consistent, which was addressed by Craster & Matar (2006). The modelling work mentioned above neglects inertia effects. Therefore, these models are only valid for small Reynolds numbers of Re ∼ O(1) and cannot account for the 'K' hydrodynamic instability mechanism. Ruyer-Quil et al. (2008) took into account inertia and streamwise viscous diffusion and formulated two coupled evolution equations for both the film thickness and volumetric flow rate. This model is valid for moderate Reynolds numbers, both small and O(1) aspect ratios of h/a.
Instability characteristics in films can be categorized by the location where instability growth can be visually detected. Transitions between different wave regimes in coating flows can be understood by the concept of the absolute (AI) and convective instabilities (CI) which was first developed in the field of plasma physics (Briggs 1964;Bers 1975) and later extended to the problems of fluid mechanics (Huerre & Monkewitz 1990). In the laboratory frame of reference, the instability is said to be convective if the disturbance is advected and eventually dies out if it is not sustained continuously at a fixed point. On the contrary, absolutely unstable flows display an intrinsic self-sustaining feature: the perturbation is so strongly amplified that it contaminates the entire flow region. Duprat et al. (2007) have studied both experimentally and theoretically the absolute and convective stabilities for a viscous film flowing down a vertical fibre. The result showed that at large or small film thicknesses, the instability is convective. Whereas, the 'RP' mode may trigger an absolute instability in an intermediate range of film thicknesses when the fibre radius is sufficiently small.
Recently, the problems of coating flows arising in more complex situations, for example in the presence of thermocapillarity (Liu & Liu 2014), subject to electric field (Ding & Willis 2019), in contact with countercurrent gas flow (Dietze & Ruyer-Quil 2015), flowing over a slippery surface (Halpern & Wei 2017;Ji et al. 2019), have been extensively investigated due to their relevance to technological applications. Except for the attempts in the complex situations mentioned above, rotation is an effective way to control the stability and dynamics of coating flows. The coating flow outside or inside a rotating cylinder arises in many industrial applications owing to the need to control, protect and functionalize surfaces (Ruschak & Scriven 1976;Moffatt 1977;Ashmore & Hosoi 2003). Dynamics of a viscous film on the outer surface of a horizontal rotating cylinder was first studied by Moffatt (1977). The flow of a liquid layer around the inside surface of a rotating horizontal drum or cylinder, referred to as the rimming flow, has been extensively investigated experimentally and numerically. Many interesting effects, including two-dimensional steady states (Ruschak & Scriven 1976;Johnson 1988), three-dimensional instabilities and time-periodic flows (Thoroddsen & Mahadevan 1997;Hosoi & Mahadevan 1999) and the effect of surface tension (Pougatch & Frigaard 2011;Benilov & Lapin 2013;Groh & Kelmanson 2014) have been represented in the literature.
The works related to the coating flow or rimming flow over a rotating cylindrical substrate mentioned above focus on the case of horizontally placed cylinders. However, a careful look at the present literature indicates that the works on the coating flows on a vertical rotating cylinder are very limited. Dávalos-Orozco & Ruiz-Chavarría (1993) studied the linear stability of a fluid layer flowing down the inside and outside of a rotating vertical cylinder. Two additional approximations, i.e. the small wavenumber and the small Reynolds number, are made in their analysis. The result showed that for flow outside the cylinder, the centrifugal force may trigger the non-axisymmetric with the azimuthal wavenumber m = 1 as the dominant mode of instability. The work by Dávalos-Orozco & Ruiz-Chavarría (1993) work does not apply to the nonlinear regime of the problem. Chen, Chen & Yang (2004) investigated the stability of a condensate film on the outer surface of a rotating cylinder in the framework of the lubrication approximation and showed that the effect of increasing rotational speed is destabilizing. Rietz et al. (2017) experimentally investigated the evolution of a thin film on the outside of a vertical rotating cylinder of large radius. The experimental results indicated that the Rayleigh-Taylor instability induced by destabilizing body force result in the emergence of rivulets. Directly after destabilization of two-dimensional waves into rivulets, detachment of droplets from the cylinder has been observed.
For the coating flows on a cylinder, three different nonlinear flow regimes were observed in experimental studies, depending on the flow rates, when the fibre is stationary (Kliakhandler et al. 2001), which were labelled flow regime 'a', 'b' and 'c'. Flow regime 'a' consists of long-spaced big droplets, which can be represented by a steady travelling-wave solution (Craster & Matar 2006). Flow regime 'b' is a necklace-like structure -close and regularly spaced droplets, which cannot be described by the long-wave model. Flow regime 'c' is a time-dependent flow, in which a big droplet coexists with several small droplets. Very recently, Ding & Willis (2019) used a relative periodic solution to describe the flow regime 'c', which showed a very similar dynamics to the experimental observation. However, it remains unknown that how the rotation influences the nonlinear dynamics, i.e. steady travelling-wave solutions (Craster & Matar 2006) and relative periodic solutions (Ding & Willis 2019), which will be explored in the present study.
In the present work, we consider the problem of a gravity-driven coating flow outside a vertical rotating fibre. Firstly, we will perform a linear stability analysis of the Navier-Stokes equations and examine the stabilities of axisymmetric and non-axisymmetric disturbances. The linear stability analysis demonstrates that the Coriolis effect is negligible for the physical parameters considered in the present system.
We investigated the effect of rotation on the stability of the most unstable mode. At low frequency of rotation, the most unstable mode is in the form of axisymmetric disturbance. For axisymmetric cases, neglecting the Coriolis effect allows us to derive a long-wave model, which is valid for both thin and thick films. To validate the long wave model, we compare the linear stability analysis of the long-wave model with the linear stability analysis of the Navier-Stokes equations. Further, the nonlinear evolution of the long-wave model for the axisymmetric case is compared with the results of direct numerical simulation (DNS). At high frequency of rotation, the most unstable mode is achieved by streamwise uniform and non-axisymmetric disturbance. In this case, no asymptotic model can be obtained for thick films. We studied the nonlinear evolution for non-axisymmetric disturbance by DNS in the appendices.
The rest of the present paper is organized as follows. In § 2 the mathematical formulation of the physical model is presented. The linear stability analysis of the three-dimensional (3-D) problem of the Navier-Stokes equations is performed to evaluate the contribution of Coriolis forces and to examine the most unstable mode for various rotation frequencies.
Further, the evolution equation is derived in the framework of long-wave approximation for the axisymmetric problem. In § 3 we study the linear stability of the long-wave problem. The effect of rotation on the time growth rate and frequency of the disturbance is studied by a temporal stability analysis. A spatio-temporal stability analysis is performed to examine the characteristics of absolute and convective instabilities. In § 4, we investigate the travelling-wave solutions and their instabilities. The simulation of travelling waves perturbed by the small disturbance shows that the travelling wave loses stability to a 'relative' periodic solution. Using a bifurcation analysis with branch continuation, the relative periodic solutions are tracked as parameter changes. In § 6, conclusions are provided.

Mathematical formulation
As shown in figure 1, a Newtonian fluid, of dynamic viscosity μ and density ρ, flows down a vertical fibre of radius r = a under gravity g. The cylinder rotates about its axis (z-axis) at an angular speed Ω. The initial radius of the fluid ring measured from the centre of the fibre is r = R. Let (r, θ, z) be the rotating cylindrical coordinates moving with the fibre.
The dynamics of the flow are governed by the Navier-Stokes equations, in which p is the pressure and u is the velocity vector with u, v, w denoting the components in r, θ , z directions. In the rotating coordinates, two additional forces, i.e. the centrifugal ρΩ 2 re r and Coriolis force 2ρΩe z × u, have been added as external forces in the governing equations. At the fibre wall, r = a, the velocity satisfies no-penetration and no-slip conditions At the free surface, r = S(θ, z, t), the balance of the normal stress is where the surface tension, σ , is constant. The balances of the shear stress are t 1 · T · n = 0, (2.7) t 2 · T · n = 0, (2.8) where T is the stress tensor T = μ ∇u + (∇u) T , (2.9) t 1 , t 2 and n are the unit vectors tangential and normal to the interface and 2H is the principal curvature of the surface (2.14) The kinematic boundary condition on the free surface is (2.15) or in conservative form (2.16) 2.1. Scalings The non-dimensionalization is based on the viscous-gravity velocity, V ≡ ρgR 2 /μ, and the capillary length, L = σ/ρgR. We assume that the radius of the fluid ring, R, is much smaller than the characteristic length L. The scalings are as follows: where the aspect ratio parameter = R/L is a small number, and ' * ' denotes non-dimensional variables. Here, we assumed that the azimuthal velocity is also small, For convenience, we drop the * decoration and the dimensionless Navier-Stokes equations are expressed as where the Reynolds number is defined as Re = ρVL/μ, and the parameters related to rotation are Ω N = Ω 2 R 2 /gL and α = V/ΩR. This non-dimensionalization naturally results in = ρgR 2 /σ , which is also referred to as the Bond number. The choice of small implies a small Bond number and surface-tension-dominated case. The full set of dimensionless parameters are , a, Ω N , α and Re.
We choose the set of parameters used in the experiments by Kliakhandler et al. (2001), in which castor oil with density ρ = 0.961 × 10 3 kg m −3 , kinematic viscosity ν ≡ μ/ρ = 4.44 × 10 −4 m 2 s −1 , surface tension σ = 3.1 × 10 −2 kg s −2 , is used for coating. The radius of the fluid ring in initial state, R, is approximately 1 mm. For this set of physical parameters, the dimensionless parameter is less than 0.3 and the Reynolds number Re ∼ 10 −2 . We fixed the radius of the initial surface to be approximately 1 mm and change the angular frequency of rotation Ω = 2πf , in which f is the frequency measured in cycles per second or Hertz. As f changes in the range from 10 to 30 which is easily realized in experiments, Ω N = 0.1223, 0.4890, 1.1003, α = 0.3545, 0.1772, 0.1182 for f = 10, 20 and 30. When deriving the long-wave evolution equation, we assume the parameter Ω N has order O(1).

Linear stability
In the framework of normal mode analysis, the variables in the full Navier-Stokes equations are decomposed in the form of [u, v, w, p, S] in which '¯', denotes the base state, 'ˆ', denotes the Fourier amplitudes of the disturbances, the integer m and the real number k denote the azimuthal and streamwise wavenumbers and ω denotes the frequency. The basic states of u, w and p arē The governing equations for the disturbances are as follows: where D = d/dr. At the fibre wall (r = a), the boundary conditions arê At the perturbed surface, boundary conditions are projected to r = 1 and after linearising, we havep The fully linearized Navier-Stokes equations are solved by a Chebyshev collocation method. Figure 2 shows the dispersion relations of the axisymmetric mode (m = 0) and the non-axisymmetric mode (m = 1). At f = 0, as shown in figure 2(a,d) for a = 0.3 and 0.8, the axisymmetric mode is the most unstable. As f = 10 in figure 2(b,e), the non-axisymmetric mode with m = 1 becomes unstable in the long-wave range, but the maximum growth rate of m = 1 is less than that of m = 0. This indicates that the axisymmetric mode dominates the instability when the rotating speed is low. As f increases to 20, as shown in figure 2(c, f ) the asymmetric mode m = 1 becomes the most unstable mode.
In figure 3, we compare the dispersion relation of axisymmetric disturbances predicted by the fully linearized Navier-Stokes equations with that of neglecting the Coriolis forces, i.e. setting α to be zero in (2.27) and (2.28). Figure 3 shows that the influence of the Coriolis forces on the stability is negligible when the rotating speed is low such that Ω N ∼ O(1). Therefore, in the following discussions, the Coriolis effect has been neglected.
We will examine the nonlinear evolution of the most unstable mode for both low and high frequencies of rotation. At low frequency of rotation, the most unstable is axisymmetric, i.e. m = 0. In this case, we derive an asymptotic model equation to describe the dynamics of the flow. In the fast rotation case, the results of the linear stability analysis have shown that the m = 0 is not the most unstable. We should note that there exists different transitional routes in the present system, which depends on the initial conditions. When the rotation is slow, if the initial condition takes the non-axisymmetric form and has no streamwise component (the m = 0 component), it can develop into a non-axisymmetric state because it is unstable. However, when the initial condition has a streamwise component, we observe fthat the final saturated state is axisymmetric, i.e. the non-axisymmetric mode is suppressed by the axisymmetric mode at small rotation number (see appendix B). When the rotation is fast, if we increase the rotation gradually  and allow the axisymmetry mode to appear first, the non-axisymmetric mode (m = 1) will be bypassed, which has also been demonstrated by the simulation of a 3-D thin film model in appendix B. However, for a uniform film, if one suddenly rotates the cylinder at a high speed, the non-axisymmetric mode may appear first, which is confirmed by direct numerical simulation in appendices A and B. The axisymmetric state, however, might also be unstable to azimuthal disturbances, when the film is very thin. For example, Rietz et al. (2017) observed that axisymmetric ring waves lose instability to a rivulet structure. A similar transitional scenario has also been found in a liquid film flowing underneath an inclined plane (Charogiannis et al. 2018). It should be noted that the Rayleigh-Plateau mechanism is very weak and negligible in Rietz et al. (2017) due to the small thickness of film. In the present work, both the Rayleigh-Plateau mechanism and the Rayleigh-Taylor instability are important because the film is thick. It will be very challenging to investigate the rivulet instability in the present case because deriving a simple model equation is not possible. However, the present axisymmetric model can be applied to understand the nonlinear wave dynamics on the crests of rivulets formed on the fibre, e.g. the interaction between different waves leading to an oscillating state.

A long-wave model
The linear stability analysis shows that the Coriolis effect is negligible at small rotation frequency, Ω N = O(1). The axisymmetric unstable mode dominates the primary instability when the rotating frequency is low f ≤ 10 and surface tension effect is dominant (small ). Assuming 2 1, Re ∼ O(1) and slow rotation, we can neglect the inertial terms and the Coriolis force.
Under the assumption of axisymmetry, the leading-order Navier-Stokes equation is given by At the fibre wall r = a, w = 0. (2.37) The leading-order tangential and normal stress balances at the surface (r = S) are It is notable that the term, 2 S zz , which involves the highest derivative of S is included in (2.39). In a linear analysis, the inclusion of this term is vital to give the correct cutoff wavenumber. Conventionally, this term is kept in long-wave theories of the problem with a cylindrical free surface, such as jets, threads and liquid bridges (Eggers & Dupont 1994;Craster & Matar 2006). An understandable explanation for the inclusion of this term is that the streamwise curvature is sufficiently large in some regions at the interface such as the steep front edge of a solitary hump. So, the term which contains the highest derivative of S is sufficiently large and cannot be neglected even though it is multiplied by 2 . For more justification of the inclusion of this term, we refer the reader to Craster & Matar (2006). An alternative choice is to keep the full curvature in the equation of the balance of the normal stresses in other related topics, for example, free-surface flows with large slopes (Snoeijer 2006) and rim or coating flows (Lopes, Thiele & Hazel 2018). From (2.35) and (2.36), we obtain the pressure (2.41) Substituting w into the kinematic boundary condition yields an evolution equation for S(z, t) given by (2.42) Equation (2.42) is identical to Craster & Matar (2006) when the rotation is turned off.
2.4. Thin film limit For the case in which the film is thin relative to the fibre, i.e. S(z, t) = a + h(z, t), where the film thickness h a, expanding equation (2.42) for small h results in the evolution equation for the thin film For very thin fluid films, a → 1 and 1 + h/a → 1, the evolution equation can be further simplified as At Ω N = 0, (2.44) is reduced to the evolution equation of thin film by some authors (Frenkel 1992;Kalliadasis & Chang 1994). The roles of capillary forces and rotation on the stability of the film are manifested in (2.44). The terms ∂ z h, Ω N ∂ z h, 2 ∂ zzz h denote the contributions of azimuthal curvature, the rotation and the streamwise curvature, respectively. It is well known that the azimuthal curvature and the streamwise curvature are destabilizing and stabilizing, respectively (Kliakhandler et al. 2001;Craster & Matar 2006). From (2.44) for the thin film limit, it can be seen that the rotation is destabilizing the interface through the azimuthal curvature.
3. Linear stability analysis of the long-wave model 3.1. Dispersion relations Let us now consider the linear stability of the problem. A small periodic disturbance in the streamwise direction is imposed on the film such that the film thickness can be decomposed into a base state componentS, and a small disturbance of normal mode with amplitudeŜ, in which ω is the complex frequency and k is the wavenumber, andS = 1. Substituting (3.1) into (2.42) yields the dispersion relation At Ω N = 0, this dispersion relation is identical to that of Craster & Matar (2006) for a coating flow down a fibre. In the thin film limit as a → 1, the dispersion relation corresponding to (2.44) tends to We use ω r and ω i to denote the real part and imaginary part of ω. Here, ω i and ω r correspond to its growth or decay rate and the frequency of the disturbance. The phase speed of the disturbance is defined as c = k −1 ω r . From (3.2), the phase speed is positive and independent of the wavenumber k and the rotation parameter Ω N . This means all disturbances with different wavenumbers propagate downstream at a same speed. Craster & Matar (2006) have tested the validity of their model by comparing the results of the dispersion relation predicted by the long-wave model and the linear stability characteristics of the Stokes flow equations. The result showed that the model equation provides a reasonably good approximation of the full Navier-Stokes equations. In the present paper, we will test the validity of the model (2.42) with the presence of rotation. The details of the stability analysis of the full Navier-Stokes equations are presented in § 2.2. In order to know quantitatively the effect of rotation, in figure 4 we plot the curves of the dispersion relations (3.2) of the long-wave model and those via fully linear stability analysis. We compared the dispersion relations for both models for two typical cases, thick film with a = 0.5 and thin film with a = 0.9. Inspection of the dispersion curves indicates that the prediction of the long-wave model is in reasonable agreement with that predicted by the Navier-Stokes equations for various Ω N . For a fixed k, both models predict that the growth rate ω i increases with the increase of Ω N . In figure 4(a), we observed that for thick film the time growth rate predicted by the long-wave model is obviously higher than that predicted by the linear stability analysis. However, in figure 4(b), the agreement between both models is excellent for the case of thin film.

Absolute and convective instabilities
For a general hydrodynamics stability problem, the wavenumber and frequency satisfy a dispersion relation in the form of D(k, ω) = 0. (3.5) In this subsection, we carry out a spatio-temporal stability analysis and both the wavenumber k and the frequency ω are complex numbers. The stability properties of the perturbations are determined by the long-time behaviour of an impulsive response to a localized excitation. The solution of the impulsive response can be expressed in the form of where the Bromwich contour F in the ω-plane is a horizontal line lying above all the singularities to satisfy causality, and the integration path A lies inside the analyticity band around the k-axis. The nature of CI and AI can be identified by examining the long-time behaviour of the wavenumber k 0 along the ray z/t at a fixed spatial location. This complex k 0 has, by definition, a zero group velocity ∂ω ∂k (k 0 ) = 0, (3.7) and ω 0 = ω(k 0 ) is called the absolute frequency. If Im(ω 0 ) > 0/Im(ω 0 ) < 0, the flow is said to be absolutely/convectively unstable. It should be noted that the saddle point k 0 must be identified by the Briggs-Bers collision criterion (Briggs 1964;Bers 1975), i.e. the saddle point must be a pinch point produced by the collision of two distinct spatial branches of solutions coming from the upper and lower half-k-planes. In the pinching process, when ω moves from above along the vertical line down to ω 0 , the two spatial branches k 1 (ω) and k 2 (ω) originate in this movement on opposite sides of the real k-axis and collide at k = k 0 and ω = ω 0 . For more details on the topics related to AI-CI problems, we refer the reader to a good review article by Huerre & Monkewitz (1990). The most relevant work to the present problem is the study on the AI-CI characteristics of an exterior coating film down a vertical fibre (Duprat et al. 2007). Very recently, an experimental study by Sadeghpour et al. (2019) also showed that the water-vapour transfer rate in the coating flow is dependent on the AI-CI characteristics. For example, when the flow is absolutely unstable, the water-vapour transfer rate is generally promoted by increasing the liquid flow rate. However, in the convectively unstable situation, the water-vapour transfer rate does not increase with the liquid flow rate due to the occurrence of irregular beads. From the dispersion relation, (3.2), the complex frequency ω can be expressed as a polynomial of k: (3.10) Through the transformation k = αk, β = (1 + Ω N )(α ) −2 ,ω = ω(αf 2 ) −1 , where α = [(16f 2 )/(3f 1 2 )] 1/3 , the dispersion relation reduces to a standard form as A straightforward analysis based on (3.11) then shows that the instability becomes absolute when β > β c = [(9/4)(−17 + 7 √ 7)] 1/3 ≈ 1.5 (Duprat et al. 2007). Therefore, applying β = (1 + Ω N )(α ) −2 to the condition β > β c yields the following implicit condition for absolute instability: (3.12) In figure 5(a,b), we present the AI-CI boundaries for various Ω N in the -a plane and for various a in the -Ω N plane, respectively. The range of parameter a is (0, 1). The model is valid for small , however, when presenting the results we extend the range of to 0.5. As shown in figure 5(a), for each curve with the decrease of a or the flow regime can switch from convective instability to absolute instability. With the increase of Ω N , the absolutely unstable regime extends in the -a plane. In figure 5(a), as Ω N = 0 the convective regime occupies a large portion of area in the -a plane. With the increase of Ω N , some CI regimes become absolutely unstable and the boundary extends counterclockwise in the -a plane.
In figure 5(b), it is shown that for each value of a with the increase of Ω N the flow switches from the CI to the AI regime. The results of figure 5 indicate that the increase of rotation promotes the absolute instability.

Nonlinear evolutions
In this section, we will study the nonlinear evolution from the initial state seeded with small disturbances. Equation (2.42) subject to periodic boundary conditions was solved numerically using a Fourier spectral method. The solution of S(z, t) is expanded by the Fourier series in which N is the number of Fourier modes,ŝ n (t) is the time-dependent coefficient and l is the computational length. A Fourier collocation is used to provide the discretization in space for periodic problem. We employed an implicit Gear's method for stiff problems for the time stepping and the relative error is 10 −6 at each time step (Ding & Willis 2019).

AI-CI instability
In this subsection, we performed numerical simulations on the nonlinear evolution of (2.42) to examine the effect of rotation on the transition from CI to AI regimes. We choose the set of parameters of a = 0.3 and = 0.4 which lies in the CI regime for Ω N = 0 and in the AI regime for Ω N = 0.5. The initial condition is seeded with S(z, 0) = 0.1 exp[−(z − 20) 2 /2] and the length of computational domain l = 100.
In figure 6, we present the spatio-temporal diagram for various Ω N , which may be helpful for us to understand the effect of rotation on the AI-CI characteristics. For Ω N = 0, the flow is linearly convectively unstable. In figure 6(a) it is obvious that the disturbance propagates downstream and decays at z = 20. As Ω N increases to 0.5, as shown in figure 5, the flow is absolutely unstable as the point of (0.4, 0.3) in the -a plane is below the AI-CI boundary of Ω N = 0.5. In figure 6(b), as Ω N = 0.5, the flow becomes absolutely unstable as the disturbance is amplified both downstream and upstream.

Transient simulations in long domains
In order to know the properties of the solutions which are naturally selected, we performed numerical simulations in long domains subject to a pseudo-random initial condition. Two typical sets of parameters of a and are chosen. Kliakhandler et al. (2001) and Craster & Matar (2006) have performed simulations on the nonlinear evolution for three regimes in Kliakhandler et al. (2001) and all of them are in the absolute instability regime. The first set of parameters, a = 0.255 and = 0.292 chosen in this paper, which falls in the absolute instability regime, is identical to that used by Kliakhandler et al. (2001) and Craster & Matar (2006). The second set of parameters, a = 0.55 and = 0.3, however, is picked from the convectively unstable regime when Ω N = 0, which becomes absolutely unstable when Ω N > 0.2.
The system becomes saturated after long-time evolution. In figure 7, we present the results of the first set of parameters, which is in the absolutely unstable regime. In figure 7(a-d), similar-sized large droplets are moving as a bound state. The heights of these droplets are very close but the gaps between them can vary. The increase of Ω N leads to the formation of more pronounced and closely spaced beads with smaller sizes. Moreover, the gaps between different beads become flatter and thinner. Figure 8 shows the film profiles for different values of Ω N , which demonstrates that the rotation has a significant influence on the flow pattern. In these figures, the height of the beads increases with the increase of Ω N . Comparing with figure 8(a,c), in which Ω N increases from Ω N = 0 to Ω N = 0.5, the height of the beads in figure 8(c) is obviously larger than that in figure 8(a). However, the large droplet numbers are reduced from six to two. Hence, the gap between the large drops in figure 8(c) is much longer than that in figure 8(a). This long film between the two large droplets in figure 8(c) is unstable and breaks into several tiny beads. The large droplet, however, will catch up with and consume these small beads ahead and leave behind a flat film where small beads will regenerate again.

Travelling-wave solutions and their stabilities
5.1. Travelling-wave solutions Experiments have shown different types of steadily propagating solutions in the form of highly organized droplets or bead-like solutions separated by long space (Kliakhandler et al. 2001). In this subsection, we will examine how rotation influences the characteristics of these steady travelling-wave solutions.
We solve the governing equations by moving to a travelling-wave coordinate. In the moving system, it is convenient to define where c is the wave speed of the moving coordinate. The derivatives with respective to t and z can be expressed as For a travelling-wave solution, the time derivation of S with respect to τ is zero, and the profile of the interface has the form of S = S(ξ ). The nonlinear equation (2.42) becomes Integrating equation (5.5) once, we have where q is an integral constant. Equation (5.7) subject to periodic boundary conditions will be solved in the region of −l/2 < ξ ≤ l/2, where the length of the computational domain, l, is the wavelength of the travelling-wave solution. To fix c, we impose the mass conservative condition To break the translational symmetry, one can enforce the real (or imaginary) part of the first Fourier mode to zero such that q can be determined. A code for Newton iteration based on Ding & Willis (2019) was developed to solve the present problem. The solution can be tracked from the Hopf bifurcation point (ω i = 0) with branch continuation. Typical profiles of the travelling-wave solutions are shown in figure 9. As shown in figure 9, the droplet height is increased by the axial rotation. In Craster & Matar (2006) and Ding & Willis (2019), an observation is that higher droplets are moving faster than lower droplets. However, in the present study, it is found that this conclusion holds only if the axial rotation is slow (see figure 10a, in which the wave speed c increases with Ω N ). When the rotation becomes fast (Ω N > 0.5), we observe that large droplets can move slower than small droplets (see figure 10b), which is contrary to the conclusions of Craster & Matar (2006) and Ding & Willis (2019).
To account for the physical mechanism of the abnormal phenomenon, we plot the flow field in figure 11    15 model, which considered the streamwise dissipation effect, and they found that the wave speed of the travelling-wave solution in flow regime 'c' is lower than Craster & Matar (2006). As also discussed by Ding & Willis (2019) that viscous dissipation effect can slow the droplet motion, it is suggestive that the motion of large droplets is retarded by this strong reversal flow in the hump.

Stability of the travelling-wave solution
Numerical simulations of the present problem showed that not all the travelling waves can be obtained through the evolution of an initial condition. It is natural to examine whether the travelling-wave solutions with different wavelengths are stable.
We performed a linear stability analysis on the travelling-wave solutions in the moving frame at the wave speed c. The profile of the interface is the travelling-wave solutions perturbed by an infinitesimal disturbance The linearized evolution equation is written as (5.14) The curves of the leading eigenvalue of the travelling-wave solutions versus the wavelength are shown in figure 12. For Ω N = 0, there exists a critical wavelength l c ≈ 3.5 beyond which the travelling-wave solutions are unstable. As l > l c , a conjugate complex pair of eigenvalues are found. In figure 12(a), the curves of the growth rate consists of several segments, indicating that mode switching occurs as l increases. This is also reflected by the non-continuous curves of temporal frequency ω r as shown in figure 12(b). As shown in figure 12(a), with the increase of Ω N the critical wavelength  l c increases. As Ω N increases from 0 to 0.2, the growth rate decreases as l < 12 and increases as l > 12. As Ω N increases from 0.2 to 0.6, the effect of rotation is stabilizing.
In figure 12(a), for Ω N = 0.6 the unstable region is confined in a limited band of l.
When Ω N increases further, the travelling-wave solutions become completely stable. In figure 12(b), the frequency of the most unstable mode increases with the increase of Ω N . The results in figure 12 showed that the effect of rotation on the travelling-wave solution is adverse to the linear stability of the trivial base stateS = 1 in which the rotation is destabilizing.
It is interesting to examine the structures of the most unstable mode of different segments in figure 12(a). In figure 13, the profiles of the most unstable mode are plotted for two typical sets of parameters. The structures of the disturbances are similar in these two figures. The disturbance consists of a series of ripples occupying the whole spatial domain. The disturbed surface is in the form of a large bead accompanied by a series of smaller beads. We also examine the structures of the most unstable mode for various parameters of l and Ω N . The results of the linear stability analysis indicate that the structures of different segments are similar. The slight difference is that the numbers of small ripples may increase with l. It is obvious that the profiles in figure 13 provide initial states involving the interaction between a large bead and several smaller ones.
In the experiments by Kliakhandler et al. (2001), the flow regime 'c' presents dynamic behaviours different from the steady travelling waves, in which a complicated interaction between a large droplet and several small droplets causes an oscillatory flow. It is interesting to investigate how the rotation influences the dynamics of the oscillatory flow, i.e. flow regime 'c'. To answer this question, we perform numerical simulations on the evolution initiated by a travelling-wave solution superposed of the most unstable mode. Before we study the nonlinear evolution by numerical simulations, it is helpful to discuss the nature of the stability of travelling-wave solutions. We classify the stability of travelling-wave solutions as being either convectively unstable or absolutely unstable. If the travelling-wave solutions are convectively unstable, they can tolerate localized disturbances and the flow ultimately evolves into a state in which the travelling waves interact with small ripples. If the travelling-wave solutions are absolutely unstable, these solutions can be destroyed by localized disturbances and evolve into a new state after a long time of evolution. As discussed in figure 13, the eigenvalues corresponding to the most unstable modes are in the form of a pair of conjugated complex numbers. This means the most unstable mode is convectively unstable and can escape the travelling wave downstream or upstream. For related discussion on the stability of a two-dimensional pulse on electrified falling films, see Blyth et al. (2018). In figure 14, the space-time diagrams of a long-time evolution are plotted for several typical cases. In figure 14(a) for l = 4 and Ω N = 0 in which the travelling-wave solution is unstable as shown in figure 12, the oscillatory behaviour of the height of the large droplets is obvious. In this space-time diagram, it is observed that the oscillation is both temporally and spatially periodic. In figure 14(b) for l = 4 and Ω N = 0.4, there is no oscillation in the space-time diagram, indicating the system is a steady travelling wave. In figure 14(c,d), the space-time diagrams are plotted for two typical cases with l = 7 in which the travelling-wave solutions are unstable. The characteristics of the structures in figure 14(a,c,d) are qualitatively similar. In these figures, smaller ripples arrange periodically in the z-t plane, and the loci of the large droplets oscillate periodically. The flow can be considered as an oscillatory travelling wave, which conforms to the fact that the travelling-wave solutions are convectively unstable. The numerical simulations suggest that the time-periodic state is a relative periodic orbit (RPO), which has been reported by Ding & Willis (2019) in a non-rotating system. However, it is unclear how the rotation affects the mean speed and the period of the RPO. We will proceed to § 5.3 to answer this question.

Exact RPOs
To find the exact RPOs of the present problem, we search for the recurrent solutions in the system S(z, t) = S(z + s, t + T), (5.15) where s is the spatial shift along the axis of the fibre and T is the time period. Each RPO has its own particular period T, while a travelling wave has arbitrary period T. Given an It is impossible to write down the explicit Jacobian A but an alternative way to solve this problem is to work in the Krylov subspace spanned by K n = {b, Ab, Ab 2 , . . . , Ab n−1 }. To evaluate the result of multiplication by the Jacobian, we use the approximation (5.20) in which ε is an empirical small parameter and typically ε = 10 −4 S 0 / δS 0 . The generalized minimal residual method (GMRES) iterations are continued until (5.21) where tol is the tolerance chosen in the range of 10 −3 -10 −4 . Typically, it takes 10 to 20 GMRES iterations to get a converged solution to δX 0 . The Newton iteration is terminated when Failure of the Newton iteration is commonplace when the initial guess is poor. To improve the convergence domain, a hookstep method is applied to move δX 0 into the trust region. We adapted the Newton-Krylov-hookstep algorithm by Ding & Willis (2019) and implemented the pseudo-arclength continuation technique to it for the present problem. The RPO is tracked from the Hopf bifurcation of the steady travelling-wave solutions, which is continued to other parameters via the pseudo-arclength continuation technique.
In figure 15(a), for Ω N = 0 the first RPO branch bifurcates from the travelling-wave branch at l ≈ 3.5, which corresponds to the first unstable mode as shown in figure 12. The bifurcation point of the second RPO branch locates at l ≈ 4.5. As the parameter Ω N = 0.3, as shown in figure 15(a), the first RPO branch bifurcates from the travelling-wave branch at l ≈ 4.0 with c ≈ 1.5. As shown in figure 15(a), at a given wavelength l, the RPO solution propagates slower than the travelling-wave solution. In this figure, it is found that both the first and the second branches of RPO propagate faster when Ω N = 0.3. This indicates that rotation can not only promote the mean wave speed of travelling-wave solutions but also the RPO. In figure 15(b), at a given domain size l, the time period of the second RPO solution is obviously smaller than that of the first RPO solution. In figure 15(b), it is found that the value of T of Ω N = 0.3 is much smaller than the case of Ω N = 0. The decreasing of the period means that the coalescence and breakup between a large droplet and a small droplet becomes faster. This suggests that the droplet coalescence and breakup process in flow regime 'c' can be promoted by the axial rotation.
For higher rotation speed, the travelling waves become more stable due to the rotation as discussed in § 5.2. Thus as Ω N exceeds a critical value, the RPO branch will disappear because the travelling waves are stable. As shown in figure 16, the first PRO branch is plotted for Ω N = 0.5. The relation between the RPO branch and the travelling-wave branch is similar to that in figure 15(a). In this figure, the wave speed c of the RPO is closer to that of travelling-wave solution than the cases of lower rotations in figure 15(a). We have not tracked the RPOs at higher rotating speeds because, with the increase of Ω N , the oscillation becomes weak and the RPOs are very close to the travelling waves.

Summary and conclusions
In the present paper, we have investigated the dynamics of a coating flow driven by gravity over a fibre rotating about its axis. The evolution equation for the surface is derived in the framework of the long-wave theory with an assumption that the characteristic radius of fluid ring is much smaller than the capillary length scale. The parameters a, and Ω N are involved in the evolution equation.
A linear stability analysis is performed to examine the influence of rotation on the stability of axisymmetric disturbances. The result of linear stability shows that the rotation is destabilizing without changing the phase speed of the disturbance. The growth rate and the wavenumber of the most unstable mode, λ m and k m , are proportion to (1 + Ω N ) 2 and (1 + Ω N ) 1/2 , respectively. We also performed a spatio-temporal stability analysis to examine the effect of rotation on the characteristics of the absolute and convective instabilities. The results showed that the increase of Ω N promotes the absolute instability. The result of the spatio-temporal stability analysis was validated by comparison with numerical simulations on the nonlinear evolution of the long-wave model equation. The results of numerical simulations indicate that with the increase of Ω N the convectively unstable flow can become absolutely unstable. This result is in excellent agreement with the linear spatio-temporal stability analysis.
The travelling-wave solutions were solved using the Newton method with branch continuation. The effect of rotation on the structure and speed of these solutions was analysed. When the rotation is slow, Ω N < 0.5, we observed that the travelling-wave speed is promoted by the axial rotation because the height of the droplet becomes larger. However, the travelling-wave speed starts to decrease as the rotation speed increases further because of the strong reversal flow in the wave hump, which retards the motion of droplets.
Linear stability analysis of the travelling-wave solutions was investigated, which showed that the travelling-wave solutions could lose their stabilities to an oscillatory mode when the wavelength exceeds the critical value. With the increase of Ω N , the travelling-wave solutions become more stable. When the travelling wave is unstable, the disturbance evolves into a near recurrent solution, which was identified as an RPO. We demonstrated that there exists series of RPOs, which are born from the travelling-wave solutions through Hopf bifurcations. Numerical solutions also demonstrate that axial rotation accelerates the mean wave speed of the RPOs and reduce their periods, indicating that axial rotation can promote the droplet coalescence and breakup process in the system. 1 + 2 (∂zS) 2 2 h ∂rw The rectangular domain is discretized by the Fourier-Chebyshev method, grid [z 1 ,z 2 , . . . ,z M ] × [r 1 ,r 2 , . . . ,r N ] with M and N mesh points inz andr directions, respectively. The temporal evolution of the problem is obtained by two steps: For more details of the procedure of the time stepping, we refer the readers to the work by Richter & Bestehorn (2019) in which an algorithm based on a nonlinear coordinate transformation is used for direct numerical simulations of liquid films under horizontal and vertical external vibrations. In order to validate our long-wave model, we have compared the travelling-wave solutions of the long-wave model with that of the Navier-Stokes equations for several sets of parameters. The results of the Navier-Stokes equations are the saturated states of the most unstable modes after a long-time evolution. In figure 17, we present the results of these two models for Ω N = 0 and 0.5. As shown in figure 17(a,b), the profiles of the interface of the long-wave mode are similar to those of the Navier-Stokes equations. In figure 18, we compare the maximum of the position of the interface, S max , of the long-wave model and Navier-Stokes for various Ω N . As shown in this figure, the results of the long-wave model are in good agreement with those of the Navier-Stokes equations. This means that the long-wave model is valid for the parameters in the present study.
in which the principal curvature We introduced the coordinate transformation t =t, θ =θ, r = h(θ,t)r + a, (A 24) where the thickness of the film h = S(θ, t) − a. The first-order derivatives The numerical method for solving the problem is similar to the axisymmetric case.
In the numerical simulations, the initial state is given by S(θ, 0) = 1 + δ cos(θ ), (A 27) in which δ is a small number. At the initial time, the maximum and the minimum of S locates at θ = 0 and θ = π, respectively. Because of the destabilizing effect of rotation, the film thickness increases with time near the region of θ = 0 and the decreases gradually on the opposite side. After a long-time evolution, the film thickness at θ = π approaches zero. This process corresponds to the fluid detachment from the substrate. As the minimum thickness approaches zero, the stiffness of the numerical procedure increases dramatically. In our simulations, we terminate the process of time stepping when the minimum thickness is close to 0.01. In figure 19, we present the profile of the free surface of the evolution of the non-axisymmetric mode with m = 1 and k = 0 when the film is close to the detachment. As shown in figure 19(a) for a = 0.8, the profile of the free surface is close to a circle. However, as shown in figure 19(b) for a = 0.3, the profile evolves into a flat lobe due to the destabilizing body force. The simulations on the disturbance with m = 1 and k = 0 showed the process of detachment from the substrate. However, these streamwise uniform surface is unstable due to the Rayleigh-Plateau mechanism. We can expect that in the process of detachment, there exists the formation of non-axisymmetric beads. The dynamics of this complex process needs further investigation on the stability analysis of the full linearized Navier-Stokes equations.  Firstly, we study the case of a small rotation number. We choose the parameters Ω N = 0.122 and ε = 0.3 as in figure 2. For a non-axisymmetric initial condition (there is no streamwise component), h = 1 + δ sin(θ ) or h = 1 + δ sin(θ + kz), the flow is unstable and a non-axisymmetric and a streamwise uniform state can develop. But if there is a streamwise component (we set h = 1 + δ sin(θ ) + δ sin(kz)), we observe that the final state is axisymmetric which indicates that the non-axisymmetric mode is suppressed by the axisymmetric mode (see figure 20). This study indicates that the axisymmetric state will be observed in experiments when the rotation is slow, because the axisymmetric mode is the most unstable and grows fastest. We also studied the case of a large rotation number. We choose Ω N = 0.5 and ε = 0.3 as in figure 2. We tested two different initial conditions, h = 1 + δ sin(θ ) + δ sin(kz)) and h = 1 + δ sin(kz). For the axisymmetric initial condition, we observe that the state remains axisymmetric. However, as expected, a 3-D state develops when the system is perturbed by the non-axisymmetric initial condition (see figure 21). However, if we rotate the cylinder after the formation of axisymmetric beads (or perturb the saturated axisymmetric state by non-axisymmetric disturbances), we did observe that the system reverts back to the axisymmetric state in the parametric space considered in the present study (see figure 22). This demonstrates that we may observe an axisymmetric film in experiments if we gradually increase the rotation speed or rotate the cylinder after the formation of axisymmetric beads.