Direct Simulation on the Dynamics of Liquid Films Flowing Down a Fiber

Direct simulations on the dynamics of a liquid film coating on the outer surface of a vertical fiber are performed in the present paper. A domain mapping technique has been used on solving the Navier-Stokes problem in stream function form and predicts the dynamics of the various flow regimes with remarkable accuracy. We investigated the morphologies that develop in the coating flow down fiber and analyze the effects of Bond number and fiber diameter on the maximum liquid film height and drop speed. The results showed that the Bond number and fiber diameter have a significant effect on the instability. The transition between the Plateau-Rayleigh and convective regimes occurs at the conditions when the flow transitions from absolutely to convectively unstable.

The study of the dynamical properties and scales of jet and film flow fragmentation continues to reveal many novel phenomena that continue to attract the research interest of academics (Craster and Matar (2009); Liu and Liu (2014); Xie et al. (2021)). The thin film flow on a cylindrical surface is an interesting class of mechanical problems related to liquid film flow. Rayleigh first discovered the Rayleigh mechanism, where surface tension causes instability at the interface (Rayleigh (2010)), leading to the spontaneous formation of droplets or water droplets from the liquid film. Kapitza (1949) were the first to study the dynamics of a liquid film flowing on a stationary vertical cylindrical solid. Two types of waves, so-called "periodic waves" and "isolated waves", were experimentally observed by applying perturbations to the incoming velocity of the liquid film surface. Kliakhandler et al. (2001) assuming that the film thickness is much smaller than the characteristic capillary length rather than the fiber radius, the proposed model reveals the rich kinetic behavior of the liquid film flow under low inertial forces and describes the droplet flow formed by three different modes of interfacial instability when the flow rate is changed. Craster and Matar (2006) conducted a theoretical study of the vertical fall of a liquid film along a fiber, focusing on the assumption that the Navier-Stokes equation is simplified when the film thickness is much smaller than the fiber radius and obtained a mathematical model based on the simplified hydrodynamic equation of lubrication approximation.
The flow of liquid films along vertical cylindrical fibers shows complex and unstable interfacial dynamics with different regimes. As illustrated in Fig. 1 (Ji et al. (2019)), accoding to the nozzle size, fiber radius, flow rate and liquid properties, Kliakhandler et al. (2001) categorized fully developed droplet patterns into three main regimes: (A) the convective instability regime, in which faster-moving droplets collide into slower-moving ones, and bead coalescence occurs repeatedly; (B) the Rayleigh-Plateau regime, where sliding droplets at a constant speed and bead spacing; and (C) the isolated droplet regime, in which where large droplets are widely spaced traveling beads separated by smaller droplets depending on the flow rate, nozzle size, fiber radius, liquid choice, and inlet geometry. When other system parameters are fixed, and the flow rate is varied from high to low, this may result in flow regime transitions from (A) to (B) and eventually to (C).
Previous numerical studies have considered many other factors, such as the inlet conditions, the flow rate, and the shape of the nozzle Ji et al. (2020)). However, these studies failed to accurately predict the dynamics in different flow regimes, such as the wave speeds of droplets in flow regimes 'B' and 'C'. Actually, simulating the liquid film flow phenomenon in 3D space requires a very large amount of computational resources, which is reflected in several places: (1) Using a discrete 3D computational domain with a body grid means that very many unknown physical quantities, including fluid volume fraction, fluid velocity, fluid temperature, and other derived physical quantities need to be stored and calculated; (2) The huge number of unknown physical quantities leads to a large matrix size, which means that the solution is inefficient and requires a lot of computation time; (3) The complexity of the model also causes the problem to become enormously difficult to solve. Liu and Ding solved the full Navier-Stokes equations for film flow down a fiber directly and have reproduced the various flow regimes with remarkable accuracy (Liu and Ding (2020)). We revisit the problem of film flows down a vertical fiber by direct numerical simulation of the dynamics process using a domain mapping technique and without the assumption of long-wave dynamics.
The present paper is organized as follows. In section 2, the mathematical formulation of the physical model is presented. In section 3, the stream function and numerical method are presented. In Section 4, we present the numerical results and compare them with experimental studies of flow regimes 'A', 'B', and regime 'C', which shows that the results are in excellent agreement with the experimental study. In Section 5, we summarize the results and present the conclusions.

Governing Equations
As shown in Fig. 2, the problem is that a Newtonian fluid of constant viscosity and density flows down a vertical fiber of radius r = a under gravity g. The air phase is assumed inviscid, and its dynamics are neglected in the present study.
The Navier-Stokes equations govern the dynamics of the flow,

Boundary Conditions
Periodic boundary conditions are assumed in the vertical direction. The no-slip boundary conditions are used for the fiber wall r = a , and the velocity satisfies no penetration and no-slip conditions At the free surface r = S(z, t) , the balance of the shear stress is the balance of the normal stress is where t and n are the unit vectors tangential and normal to the interface 2H is the principle curvature of the surface in which The kinematic boundary condition on the free surface is in the form of a conservation equation

Scaling and Reduction
The non-dimensionalization is based on the viscous-gravity velocity V ≡ gR 2 ∕ and the character length L ≡ ∕ gR . The scalings are as follows, in which * denotes non-dimensional variables.
Dropping the * decoration for convenience, the dimensionless Navier-Stokes equations are expressed as where Reynolds number by Craster and Matar (2006) is very small, Re = R∕ 2 the inertial terms are neglected. The surface tension is . The parameter = Bo = R∕L = gR 2 ∕ is the Bond number and Bo∼10 −1 in the present study.
The boundary conditions at r = a are The boundary conditions at the surface:

Stream Function
In the axisymmetric case, the velocity components can be expressed in terms of a stream function , We can eliminate the pressure from the two equations of motion and substitute the above expressions for u and w, and give the linearized differential equation for , in which The boundary conditions at r = a are where and We can eliminate p using by performing ∇ s in which

Mapping
We introduced the coordinate transformation where k × z is the computational space, the thickness of the film h = S − a. This set of mapping transforms an arbitrarily shaped flow region to a rectangular region. The first-order derivatives and the second-order derivatives At the surface are

Algorithm
The idea of solving partial differential equations by spectral methods is as follows: firstly, discretize the spatial variables in the solution domain; secondly, use the Fourier spectrum or Chebyshev spectrum to derive the partial differentiation of the matrix treatment function to the space; after which the partial differential equations are transformed into a system of ordinary differential equations; finally, use the time stepping method to obtain the numerical solution at the specified inscription by numerical iteration with a given boundary edge strip. The computational steps for solving the differential equations using the Chebyshev spectral method are as follows: Step 1: Discrete the direction of Z-axis into N equally spaced locations within the computational interval. To improve computational accuracy and efficiency, Chebyshev non-uniform point discretization is used in the radial direction, and then select an initial state data.
Step 2: Given the boundary conditions and solving for the liquid film thickness of the lattice points in the computational interval.
Step 3: Calculate the pressure at each grid point by solving a sparse linear system.
Step 4: End the loop and save the results if the relative error of the velocity is less than the set value, or repeat if the velocity does not converge.
Step 5: Return to step 2 and solve the loop 200000 iterations using the time stepping method.

Grid-independent Validation
In the process of numerical calculation, the size of the grid number will have a certain influence on the calculation results. In order to reduce the influence of the grid number (41) on the calculation results of numerical simulation, it is necessary to verify that the calculation here is independent of the grid number. In this paper, structural grids are used in the computational domain, and three numbers of grids, 32x10, 64x20, and 128x30, are used to test the simulations under the same conditions. The calculation results of grid-independent verification are shown in Table 1. The results show that the calculation results of the three grids tend to be the same as shown in Fig. 3, and the maximum relative error is only 0.26% . The computer used for the calculation is a Lenovo brand machine with Intel Core i7 6700 series CPU 3.4GHZ, DDR4 memory 8GB, Windows 10 operating system, and Matlab2016 as the calculation software. The time spent on calculating the three grids differs greatly, and the time spent on the calculation increases from 15 to 1292 minutes when the number of grids doubles. And it should be noted that the calculation results are unstable when the number of networks in the vertical direction is less than 20.

Compared with Experimental Results
In the numerical calculation of simulating the liquid film flow of the nozzle containing fiber core, the surface wave thickness and velocity of the liquid film are two important parameters. In order to verify the validity of the model and calculation method constructed in this paper, the calculation results are compared and analyzed with the experimental results.
In this section, simulations are conducted to assess the validity of the model for simulating the dynamics process for liquid films flowing down a fiber. Figure 4 shows the profiles of surface travelling wave compared with numerical solution computed using a domain mapping technique. The results show that under the same working conditions, the calculated result data is more consistent with the trend of experimental data Craster and Matar (2006), and the difference of liquid film thickness at different locations is smaller.
To test the new approach for direct numerical simulation of liquid films, we have carried out numerical simulations and made a reasonable comparison between theoretical study and experimental measurement. We compute the corresponding wave speed c and droplet profile. The  Table 2  Flow regime 'C' occurs at a very low flow rate, and a large droplet is observed to coexist with several tiny beads. The travelling wave solutions showed that the wave speeds predicted by the model for flow regimes 'A' and 'C' are faster than the wave speeds predicted by Craster and Matar (2006), which may be due to the inclusion of inertia.

Analysis of the Surface Waves
The liquid film flow state is a non-stationary process in which the thickness, velocity, and behavior of the free interface of the liquid film change with time, and these parameters have an important influence on the flow characteristics of the flow film. Therefore, in this section, we investigate the influence of instability on predicting of the dynamics of the thin film flowing. Data processing is performed on the calculated results from simulations so that the thickness, velocity, and flow state of the liquid film can be analyzed and studied. Here we see that the maximum liquid film height significantly decreases with increased fiber diameter. Figure 5(b) plots the variation pattern of the drop speed c against the fiber diameter, showing a rapid decrease in velocity with increased fiber diameter. The dependence of the drop speed on fiber diameter agrees well with the theoretical prediction. Figure 6 plots the variation pattern of the profiles of interface wave, showing that increased fiber diameter leads to decreased maximum liquid film height, as expected. Where increases to the fiber diameter, especially for = 0.3, will result in transitions from the isolated to Plateau-Rayleigh to convective regime. The transition between the Plateau-Rayleigh and convective regimes occurs at the conditions when the flow transitions from absolutely to convectively unstable. Figure 7 shows the variation pattern of the maximum liquid film height and speed c with Bo number. It is obvious that there is a tendency for the maximum surface height to decrease as the Bo number increases. The speed c decrease as Bo increases; the relative importance of surface tension decreases, and the amplitude of the droplets also decrease. At the same time, we also discovered increasing Bo results in a large quantitative change in the profile, especially when Bo = 0.3, the isolated droplet starting transition to the convective instability regime. We have also investigated the effect of Bo on speed c. Figure 7(b) shows solutions for regimes A, B, and C. When increasing Bo number, the speed c reduces; thus increasing Bo can be interpreted as a decrease in the relative significance of surface tension and gravitational forces.
We provide numerical predictions for parameter values that correspond to different Bo values. The interfacial profiles were obtained as shown in Fig. 8. Comparison of Fig. 8(a)-(c), respectively, reveals that there is a tendency that the surface thickness decreases with the increase of Bo number, and the large Bo number makes the liquid film surface tend to be stable. It is also found that the effect is more obvious in flow regime 'A' and regime 'B'. It is important to note that when Bo = 0.3, the stabilization is completed. It is also notable that when plotted as S versus Z, the droplets have a non-dimensional length of order unity, indicating that the choice of length scale in the nondimensionalization is the appropriate one.

Coating Flow with Slip Boundary
The behaviors of wall slip at the microscope will significantly affect the behaviors of film flows. In this subsection, we will use the numerical scheme to investigate the effect of wall slip on the behaviors of the coating flow, wall slip boundary conditions to be considered in the model. The slip boundary condition at the interface r = a is specialized to: where is the slip coefficient, n is the normal direction to the surface. Figure 9 plots the variation of liquid film velocity with slip coefficient Liquid film surface wave velocity and maximum liquid film height both increase with increasing slip coefficient. It shows that the slip coefficient has a significant effect on the instability of the liquid film flow, with an increase in the slip coefficient causing the surface of the liquid film to move more towards an unstable pattern.
To further investigate the evolution of the slip coefficient on the surface liquid film flow pattern, keeping the other parameters constant, by constantly changing the values of . Figure 10 plots the variation of liquid film surface waves with in two unstable modes, A and B. The results show that the slip coefficient has a greater effect on the A model. An increase in the surface wave height of the liquid film from 2.7 to 3.3 as the slip coefficient increases from 0.1 to 0.4, and for mode B, the liquid film surface wave height increases from 1.9 to 2.2.
In summary, the wall slip has a large impact on the liquid film morphology, making the liquid film more unstable. The calculation results show that the designed model is able to

Summary and Conclusions
This paper presented a stream function approach for deriving a simplified model for the gravity-driven flow of an axisymmetric liquid film along a vertical fiber. The dynamics process of thin film flows down a vertical fiber by direct numerical simulation using a domain mapping technique. An algorithm based on a nonlinear coordinate transformation is presented that allows for direct numerical solutions of the fully nonlinear Navier-Stokes equations and appropriate boundary conditions. No surface tracking is necessary. By comparing the different regimes profiles and speed values with the experimental results, it is indicated that the presented method is highly consistent with the actual value.
In order to further analyze the morphologies that develop in the flow of thin liquid films down fiber, the effects of Bond number and fiber diameter on the maximum liquid film height and drop speed are analyzed. The results showed that the Bond number and fiber diameter significantly affect the problem and can change the velocity and the maximum liquid film height of the bead. In between these parameters, the Bond number has the most effect on the instability of the fluid flow. An increase in Bond number and fiber diameter will result in transitions from the isolated to Plateau-Rayleigh to the convective regime. The transition between the Plateau-Rayleigh and convective regimes occurs at the conditions when the flow transitions from absolutely to convectively unstable. Numerical experiments have shown that the new method has better stability in the investigation into bead dynamics. Thus this work provides an initial investigation into a topic with many fruitful areas yet to be explored, which can be used as a design tool for 3D printing, heat, and mass transfer applications.