Coating flows down a vertical fibre: towards the full Navier–Stokes problem

Abstract We revisit the dynamics of a thick liquid film flowing down a vertical fibre. Instead of deriving a long-wave model, we directly solve the Navier–Stokes equations using a domain mapping technique and the exact steady travelling wave solutions are explored using a dynamical systems theory approach. Three distinct flow regimes, labelled ‘a’, ‘b’ and ‘c’, observed in previous experiments (Kliakhandler et al., J. Fluid Mech., vol. 429, 2001, pp. 381–390) are investigated. Flow regime ‘a’ refers to a steady flow state in which large droplets are separated by a long thin film. Flow regime ‘b’ is a necklace-like flow. In flow regime ‘c’, a cyclic process of droplet coalescence and breakup was observed. By matching the mean flow rates of the travelling wave solutions and experimental data, our travelling wave solutions show an excellent agreement with flow regimes ‘a’ and ‘b’. The time-periodic flow regime ‘c’ is compared with direct numerical simulation of the Navier–Stokes equations. A snapshot of the simulation shows a remarkable similarity to an experimental image and the discrepancy of mean wave speed and maximal wave height between our numerical simulation and experimental data is negligible.

can correctly predict the dynamic behaviour (droplet coalescence and breakup process), the improvement in the predicted wave speed is not significant. A plausible reason for the inaccurate prediction of the wave speed of flow regime 'c' by the single-equation model is the neglected inertial term, which can be important because of the large size of the droplet. Ruyer-Quil et al. (2008) and Novbari & Oron (2009) took inertia into consideration and proposed a weighted residual model -a two-equation model -that aims at improving the prediction of the nonlinear dynamics when the Reynolds number is moderate. However, the two-equation model also failed to predict the correct dynamics of flow regimes 'b' and 'c'. But the two-equation model was more successful for thin liquid films on an inclined plane in that it compares well with the full Navier-Stokes equations and experiments (Denner et al. 2018;Tomlin, Cimpeanu & Papageorgiou 2020). This perhaps is because the wavelengths in these problems are typically long. It should be stressed that the single-equation model or the two-equation model becomes immediately 'wrong' because the long-wave assumption breaks down when the film is thick, e.g. flow regime 'b'. In addition, in some industrial applications, fluid inertia cannot be neglected and the droplet spacing is short, e.g. desalination ; the reduced models cannot be applied and an accurate prediction of the flow dynamics is in demand.
To investigate the nonlinear dynamics when the reduced models are not valid, the Navier-Stokes equations need to be solved using a proper numerical method, which can also help to validate the working regimes of the reduced models. However, solving the Navier-Stokes equations with a moving interface is challenging. In recent years, popular numerical methods, such as the volume of fluid method, level set method and phase field method, have been widely used to handle the moving boundary problem. These methods solve a multi-phase flow problem by modelling the surface tension as a body force within a very narrow layer that separates the liquid phase and air phase. The dynamics of the air phase is coupled with the liquid phase, which cannot be neglected. For example, the volume of fluid method has been applied to simulate a liquid film flowing on a flat plane (Denner et al. 2018;Tomlin et al. 2020) and a liquid film coating a vertical fibre (Sadeghpour et al. 2017, when the fluid inertia cannot be neglected. An advantage of the volume of fluid method, or phase field method, is that it can deal with the case when the fluid interface cannot be described by a function, e.g. dripping of liquid films down an inclined plane due to the Rayleigh-Taylor instability (Tomlin et al. 2020). However, the reduced models were derived from a single-phase flow problem, which neglects the dynamics of the air phase, and it is straightforward to compare the reduced model with its prototype such that we can assess the validity of the long-wave assumption and influence of the neglected streamwise diffusion. But if the dynamics of the air phase is important (e.g. a turbulent air flow), we have to compare the model equations with the full Navier-Stokes equations of a two-phase flow system. In our present study, we focus on the Navier-Stokes equations of a single liquid phase flow and neglect the dynamics of the air phase, such that the popular numerical methods of two-phase flows are not employed. Investigating the nonlinear dynamics of a single liquid phase flow with a moving interface requires a special numerical method. Over the past decades, the finite element method (FEM) and domain mapping technique have been widely used. Bach & Villadsen (1984) investigated a two-dimensional liquid film falling down a vertical wall and solved the problem using the FEM in terms of the Lagrangian velocities. Using the FEM, Salamon, Armstrong & Brown (1994) performed a comprehensive study of the travelling waves in vertical thin films, which was compared against an asymptotic long-wave model and a boundary layer model. Ramaswamy, Chippada & Joo (1996) used the FEM with an arbitrary Lagrangian-Eulerian formulation to study the complex wave interactions in a vertically falling film. Their numerical results showed that the tendency of the film to break up can be arrested by wave merging and splitting processes. However, the grid has to be updated at every time step using the FEM, which slows down the simulation. Very recently, Richter & Bestehorn (2019) used a domain mapping method to simulate the dynamics of a vibrating liquid film on a horizontal plane. The domain mapping method is to map the time-varying domain into a time-independent rectangular domain, which has been recently extended to a liquid film flowing down a rotating cylinder by Liu & Ding (2020). The main advantage of the domain mapping method is that there is no need to update the grid at every time step.
In the present paper, we revisit the problem of coating flows down a vertical fibre by solving the Navier-Stokes equations without the assumption of long-wave dynamics. The rest of the paper is organized as follows. Mathematical formulation and numerical method are presented in § 2, which is followed by results and discussion in § 3. The travelling wave solutions are explored using a dynamical systems theory approach, which are compared with experimental studies of flow regimes 'a' and 'b'. Flow regime 'c' is a time-periodic travelling wave (or defined as relative periodic orbit by Ding & Willis (2019)) which cannot be represented by a steady travelling wave. However, it is very challenging to find an exact relative periodic orbit of the Navier-Stokes equations to compare with flow regime 'c'. Instead, we perform numerical simulation by designing an initial condition, which shows that the system evolves into a time-periodic state and is in excellent agreement with the experimental study. Conclusions are presented in § 4.

Mathematical formulation and numerical method
We consider a Newtonian fluid, of constant viscosity μ and density ρ, flowing down a vertical fibre of radius r = a under gravity g. The air phase is assumed inviscid in the present study and its dynamics is neglected. The surface tension is σ , which is constant in the present study. We directly start from the dimensionless governing equations by choosing the mean radius R and the capillary length L = σ/ρgR as the radial and axial length scales, V = ρgR 2 /μ as the velocity scale and L/V and ρgL as the time scale and pressure scale.
In the present study, we restrict ourselves to axisymmetric flow and a two-dimensional problem is considered. Then, choosing cylindrical coordinates, the dimensionless governing equations are where La = σρR/μ 2 is the Laplace number, = R/L = ρgR 2 /σ is the Bond number and u = ue r + we z is the velocity. At the surface r = S(z, t), the stress is balanced by the surface tension: The moving surface satisfies the following kinematic condition: where α ≡ a/R is the dimensionless radius of the fibre. There is no slip and penetration at the fibre wall.
To solve the Navier-Stokes equations (2.1)-(2.3) subject to a moving boundary at r = S(z, t), we employ a domain mapping technique adapted from Richter & Bestehorn (2019) (see Liu & Ding 2020): where h = S − α is the thickness of the film. Hence, the domain of fluid is mapped into a rectangular regionr ×z = [0, 1] × [0, l], where l is the axial length of the domain, which is also the wavelength of a typical travelling wave in the present study. The Chebyshev-Fourier spectral method was used to solve the problem. Ten Chebyshev modes and 100 Fourier modes were utilized to ensure numerical accuracy in the present study. Details of the numerical method can be found in the appendix of Liu & Ding (2020). It should be mentioned that the present study considers a periodic domain such that the nozzle effect and irregular bead coalescence dynamics in the convective instability regime will not be characterized. To investigate these complex situations, one has to consider a non-periodic long domain in the simulation (see Sadeghpour et al. 2017Sadeghpour et al. , 2019. The exact travelling wave solutions can be explored by tracking the recurrent solutions in the system (Cvitanović et al. 2018): where X encodes the variables (u, w, S). For travelling wave solutions, T is fixed at T = 1 in the present study. The spatial shift s defines how far the wave runs in the axial direction. Hence, the wave speed c of the travelling wave is defined as c = s/T. This dynamical systems theory approach has been widely used to find equilibria and periodic orbits in turbulent flows (Viswanath 2007;Chandler & Kerswell 2013;Budanur et al. 2017) and in liquid film flows (Ding & Willis 2019).
To improve the guess of (X 0 (t,r,z), s), one has to solve the following equations (Ding & Willis 2019): To fix the shift distance s, we require δX 0 to have no component that shifts X 0 in the axial direction: We used a Jacobian free method to approximate (∂X T /∂X 0 )δX 0 as where Δ = 10 −6 X 0 / δX 0 is a small empirical parameter. Equations (2.9) and (2.10) constitute a linear algebraic system A x = b wherein x = (δX 0 , δs) T and b = −(X T − X 0 , 0) T . The linear system can be solved using a GMRES method, which approximates the solution in the Krylov subspace K n = {b, A b, A 2 b, . . . , A n−1 b}. Typically, 20 GMRES iterations will yield a converged solution, i.e. the relative error A The Newton step is terminated when X T − X 0 / X 0 < 10 −6 . Our numerical method has been validated against the long-wave model and tested near the cut-off wavenumber predicted by linear stability, which showed that the predicted wave speed is in excellent agreement with the linear stability analysis.

Results and discussion
Note that the radial length scale R in previous studies (Kliakhandler et al. 2001;Craster & Matar 2006) was determined using the mean flow rate: The mean flow rate was measured by feeding back liquids into the upstream reservoir such that the liquid depth in the reservoir was maintained (Kliakhandler et al. 2001). This radial length R was considered as the initial mean radius of the liquid film in previous studies (Kliakhandler et al. 2001;Craster & Matar 2006). Here, we found that R does not really correspond to the initial mean radius of the liquid film, and the initial radius of the film has a significant influence on the prediction of the wave speed and wave height. Previous numerical studies have considered many other factors, such as the flow rate, the inlet conditions and the shape of the nozzle (Duprat, Ruyer-Quil & Giorgiutti-Dauphiné 2009; Sadeghpour et al. 2017), that affect the predictive capacity of reduced models. However, these studies failed to accurately predict the dynamics in different flow regimes, e.g. the wave speeds of droplets in flow regimes 'b' and 'c'. Actually, in the numerical simulations, the flow rate Q is varying with time once the interface becomes unstable because of the varying velocity profile and interface. This means that the flow rate Q at the laminar state will be different from the average flow rate of a saturated travelling wave state. Hence, there should be a special initial mean radius ( / = R) corresponding to the mean flow rate measured in experiments. Actually, the initial mean radius of the film is an unknown parameter in experiments (only the flow rate is given). Hence, to make a reasonable comparison between the travelling wave solution and experimental measurements, the flow rate of the travelling wave solution, rather than the initial laminar state, should be consistent with experimental study. However, before examining the effect of initial mean radius on the solution, we start from the situation thatS(z, 0) = 1, i.e. the initial radius of the film is equivalent to the radial length scale R, and compare our results of the Navier-Stokes equations with the previous long-wave models.
The travelling wave solutions are tracked using a branch continuation method (pseudo-arclength continuation method or parametric continuation method) starting from the Hopf bifurcation of the base state X = (0, (2 ln(r/α) − r 2 + α 2 )/4, 1) T . Figure 2 shows that the wave speed c and maximal wave amplitude increase with the wavelength l. This is in agreement with the long-wave models that larger droplets are running faster. Craster & Matar (2006) showed that there was a small capillary wave in front of the large droplet in the travelling wave solutions of flow regimes 'a' and 'c', which can be reduced by including the inertial effect using a weighted-residual model (Ruyer-Quil et al. 2008). In the present study, we found that these small capillary waves are absent in front of the large drops (see figure 2a,c), which is also because we included the inertial effect in the present study. However, when La, a or is large, the small capillary wave will appear in the travelling wave solution.
We also noticed that there are some other important issues remaining unclarified in flow regimes 'a' and 'c'. For regime 'a', most previous theoretical works seemingly agree well with experiment. However, we found that the wavelength of regime 'a' in some theoretical works (e.g. Ruyer-Quil et al. 2008;Novbari & Oron 2009) is not the same as the experimental value. In the experimental work by Kliakhandler et al. (2001), it was reported that the distance between the two large beads was 30 mm, which was used by Ruyer-Quil et al. (2008) and Novbari & Oron (2009). But as shown in images of regime 'a', it is apparent that the wavelength is approximately 22.5 mm, which is adopted in the present study. The travelling wave solutions showed that the wave speeds predicted by the Navier-Stokes equations for flow regimes 'a' and 'c' are faster than the wave speeds predicted by Craster & Matar (2006), which is due to the inclusion of inertia (Novbari & Oron 2009) (see table 1). The predicted wave speed c of regime 'c' by the Navier-Stokes equations is slightly lower than that of Novbari & Oron (2009), because the streamwise diffusion terms are retained in the present study, which can slow down the motion of droplets (Ruyer-Quil et al. 2008;Ding & Willis 2019). The predicted wave speed of regime 'b' is better than that for all the previous reduced-order models, which is because the long-wave assumption is not valid in this case. However, the predicted wave speeds of flow regimes 'b' and 'c' are still much higher than the experimental values, which will be examined later.
When the Bond number is not small and the Laplace number is large, neither the single-equation model nor the two-equation model is valid. To explore the dynamics in the regimes where these models do not apply, we solve the Navier-Stokes equations.
Experiment (Kliakhandler et al. 2001 Figure 3 indicates that the droplet size increases as La increases, and thereby a slightly faster running speed c, which agrees well with the prediction of the long-wave model of Craster & Matar (2006) and Ding & Willis (2019). However, when the fibre radius is small, e.g. a = 0.25, we observe a different phenomenon in that the wave speed can be slightly reduced for a larger droplet when La 12, which, perhaps, is due to the circulation flow in the large droplet (Liu & Ding 2020). Furthermore, we examined the dependence of the travelling wave solution on the Bond number by fixing the Laplace number at La = 10. Clearly, for > 0.3, both the streamwise diffusion effect and inertial effect are important. The destabilizing surface tension effect becomes weaker for larger . Both the wave height and the wave speed become lower as increases (see figure 4), which suggests that the strong streamwise diffusion effect and the weakened surface tension effect dominate the inertial effect in the present case.
To examine why the Navier-Stokes equations predict 'worse' dynamics of flow regimes 'a' and 'c' than the asymptotic model of Craster & Matar (2006), we calculated the mean flow rate of the travelling wave solutions. Figure 2 shows that the flow rate Q clearly varies with the wavelength l, demonstrating that the laminar flow rate is very different 1185. This clearly indicates the true initial mean value of S is less thanS = 1 and the true initial flow rate of the unperturbed flow rate is different from the flow rate of a saturated flow. Hence, to make a reasonable comparison between theoretical study and experimental measurement, we have to adjust the initial value ofS such that the mean flow rate of the travelling wave solution can be consistent with the experimental flow rate. In figure 5, we present the mean flow rate versus the the initial value ofS for flow regimes 'a' and 'b' (the result for flow regime 'c' is not presented because it is not a steady travelling wave solution). Obviously, a small change in the initial mean radiusS will have a significant impact on the flow rate of a steady travelling wave, which should be an important reason for the poor predictions of the droplet speeds in flow regimes 'a', 'b' and 'c'.
By doing a bisectional search of initial value ofS, we are able to match the mean flow rate of the travelling wave solutions and the experimental flow rates. The flow fields are compared with the experimental photos in figure 6(a,b). Furthermore, as shown in table 1, both the predicted wave speed and the wave height are improved over the above study of initial conditionS = 1. The predicted wave speed of flow regime 'a' is reduced, which is closer to the experimental value than that of Craster & Matar (2006). And the predicted wave speed of the travelling wave solution for flow regime 'b' is improved to c = 0.4437. The predicted wave speed of flow regime 'b' is still approximately 30 % higher than the experimental value. A possible reason for this 30 % discrepancy is that flow regime 'b' is not stable (Craster & Matar 2006), and hence a transient wave speed was measured in experiments. Another possible reason is that the liquid was contaminated during the experiments, and the influence of interfacial contaminants on dynamics should be investigated in further work. A third possible reason is that the dynamics of the air phase might be important in this regime, which cannot be neglected and the 'full' Navier-Stokes equations of the air phase and the liquid phase should be solved. Flow regime 'c' is a time-periodic state, which was suggested to be compared with a relative periodic orbit rather than a travelling wave solution (Ding & Willis 2019). Although Ding & Willis (2019) predicted the cyclic coalescence and breakup dynamical process in flow regime 'c', the mean flow rate of their solution does not match the experimental result and the mean wave speed of their relative periodic orbit is approximately 30 % higher than experiment. Unfortunately, we failed to obtain a converged exact relative periodic solution of the Navier-Stokes equations. To obtain the time-periodic state that can be compared with the experiment, we should start the simulation with a specific initial condition. In fact, the solutions of the Navier-Stokes equations are not unique. For the given parameters of flow regime 'c', other types of solutions do exist, such as multi-hump steady travelling waves and relative periodic solutions with different numbers of small beads (Ding & Willis 2019). Hence, in order to compare with experiment, we designed a numerical experiment using the relative periodic solution of Ding & Willis (2019) as initial condition (the flow rate is also averaged in a time period of the solution:Q = (1/T) T 0 Q dt). A time stepping study showed that the saturated solution is time-periodic (see figure 7), and a typical profile of the interface shows a structure remarkably similar to the experimental photo (see figure 6c). The cyclic process shows that the large droplet catches up with and consumes its front small beads, and small droplets are regenerated to the rear of the large droplet due to the Rayleigh-Plateau instability. In addition, since the time-marching method converges to the periodic state, it implies that the time-periodic solution is stable. Moreover, we computed the average wave speed c of the time-periodic state and recorded the maximal wave height and found that the relative error with the experimental data is approximately 3 %, showing that the result 914 A30-11 of the long-wave model has been significantly improved. Therefore, the present study suggests that, to predict the nonlinear dynamics in different flow regimes, one can keep the flow rate unchanged in the numerical simulation (this can be done by introducing a free time-dependent parameter in the momentum equation, such that the flow rate at each time step can be adjusted), and the flow rate should be fixed and controlled by a pump rather than maintaining the fluid depth in the reservoir in experiments. When the flow rate in the numerical study (either the Navier-Stokes equations or the reduced models) matches the flow rate in the experiment, the dynamics can be predicted correctly.

Conclusions
This paper has revisited the dynamics of a thick film coating a vertical fibre. The Navier-Stokes equations are solved using a domain mapping technique. A dynamical systems theory approach has been employed to find the exact travelling wave solutions, which were compared with experimental studies. We have shown that larger droplets move faster and the small capillary wave in front of the large droplet, appearing in the solution of long-wave models, can be eliminated in our solution of the Navier-Stokes equations.
Our numerical results also demonstrated that the wave speed and height can be enhanced by the inertial effect but reduced by the streamwise diffusion effect. A crucial observation by us was that the mean radius computed from the flow rate of a saturated state, which was chosen as the initial mean radius of the liquid film in previous studies, does not correspond to the true mean radius of the film at the initial state. We found that even a small change in the initial radius of the film can have a significant impact on the flow rates of the steady travelling wave solutions, which accounts for the poor predictions of the wave speeds in experiments. By doing a bisectional search of the initial thickness of the film, we obtained travelling wave solutions of the same flow rates as those in the experimental studies of Kliakhandler et al. (2001). Comparing with the experimental data obtained by Kliakhandler et al. (2001), the travelling wave solutions showed excellent agreements with experimental flow regimes 'a' and 'b'. This suggests that, when predicting the dynamics of a flow regime using a reduced model, it is important to match the flow rate rather than the fluid volume between the predicted solution and experiment.
The time-periodic state, flow regime 'c', is compared with a snapshot from a series of time-stepping simulations, which also showed an excellent agreement with the experimental study. The discrepancy of the mean wave speed and maximal height between the numerical solution and the experimental data has been improved to 3 %. Interestingly, our numerical simulation also showed that the coalescence and breakup process between a large droplet and small beads is a cyclic and stable process.