Liquid films falling down a vertical fiber: modeling, simulations and experiments

We present a control-volume approach for deriving a simplified model for the gravity-driven flow of an axisymmetric liquid film along a vertical fiber. The model accounts for gravitational, viscous, inertial and surface tension effects and results in a pair of coupled one-dimensional nonlinear partial differential equations for the film profile and average downward velocity as functions of time and axial distance along the fiber. Two versions of the model are obtained, one assuming a plug-flow velocity profile and a constant thin boundary layer thickness to model the drag force on the fluid, the other approximating the drag using the fully-developed laminar velocity profile for a locally uniform film. A linear stability analysis shows both models to be unstable to long waves or short wavenumbers, with a specific wavenumber in that range having a maximal growth rate. Numerical simulations confirm this instability and lead to nonlinear periodic traveling wave solutions which can be thought of as chains of identical droplets falling down the fiber. Physical experiments are also carried out on such a system using safflower oil as the working liquid and a taut fishing line as the fiber. A machine learning scheme is used to find the best set of parameters in the laminar flow model to match the experimental results to the simulations. Good agreement is found between the two, with parameter values that are quite close to their original estimates based on the approximate values of the physical parameters.


Introduction
Liquid film flows along vertical cylindrical fibers exhibit complex and unstable interfacial dynamics with distinct regimes. Driven by the effects of Rayleigh-Plateau instability and gravity, a wide range of dynamics can be observed experimentally. These include the formation of discontinuous bead-like droplets, periodic traveling wave-like patterns, and irregularly coalescing droplets. The study of these dynamics has widespread applications in heat and mass exchangers, desalination , and particle capturing systems (Sadeghpour et al. 2017), attracting much attention over the past two decades.
Depending on flow rate, liquid choice, fiber radius, and inlet geometry, three typical flow regimes have been observed (Kalliadasis & Chang 1994;Ji et al. 2020): (a) the convective instability regime, where bead coalescence happens repeatedly; (b) the traveling wave regime, where a steady train of beads flows down the fiber at a constant speed; and (c) the isolated droplet regime, where widely spaced large droplets are separated by small wave patterns. If other system parameters are fixed, and flow rate is varied from high to low, this can lead to flow regime transitions from (a) to (b), and eventually to (c). Further analysis of the traveling wave patterns in regime (b) is expected to provide insights into many engineering applications that utilize steady trains of beads.
For small flow rates and thin films, classical lubrication theory is typically used to model the dynamics of axisymmetric flow on a cylinder. When the fluid film thickness is significantly smaller than the cylinder radius, Frenkel (1992) proposed a weakly nonlinear thin-film equation to calculate the evolution of film thickness h (or the height of the film) and capture both stabilizing and destabilizing effects of the surface tension in the dynamics. This evolution equation was further studied by Kalliadasis & Chang (1994), Chang & Demekhin (1999), and Marzuola et al. (2019). Craster & Matar (2006) developed an asymptotic model which relaxes the thin film assumption, instead requiring that the film thickness be smaller than the capillary length. Kliakhandler et al. (2001) extended the thin film model to consider thick layers of viscous fluid by introducing fully nonlinear curvature terms. Recently Ji et al. (2019) investigated a family of full lubrication models that incorporate slip boundary conditions, fully nonlinear curvature terms, and a film stabilization mechanism. The film stabilization term, Π(h) = −A/h 3 with A > 0, is added to the pressure and is motivated by the form of disjoining pressure widely used in lubrication equations (Reisfeld & Bankoff 1992) to describe the wetting behavior of a liquid on a solid substrate, and the scaling parameter A > 0 is typically selected based on a stable liquid layer in the coating film dynamics. Numerical investigations of experimental results in  showed that compared to previous studies, the combined physical effects better describe the propagation speed and the stability transition of the moving droplets.
For higher flow rates where inertial effects are significant, coupled evolution equations of both the film thickness and local flow rate are developed (Trifonov 1992;Ruyer-Quil et al. 2008Novbari & Oron 2009). These equations incorporate inertia effects and streamwise viscous diffusion based on the integral boundary-layer approach. Recently, Ji et al. (2020) further extended a weighted-residual integral boundary-layer model to incorporate the film stabilization mechanism to address the effects of the inlet nozzle geometry on the downstream flow dynamics. Finally, Liu & Ding (2021) have solved the full Navier-Stokes equations for film flow down a fiber directly using a domain mapping technique and have been able to reproduce the various flow regimes with remarkable accuracy.
In this work we present a careful derivation of a simple new two-equation model, not starting from the Navier-Stokes equations, but based on a control volume analysis of the conservation of mass and momentum equations. While the approach is better justified if the axial velocity profile is a plug flow (i.e., uniform in the cross section as might be the case for a well-mixed turbulent flow), we can treat the drag force on the liquid film by the fiber wall in the laminar regime as well, obtaining a simple model that is suitable for viscous low Reynolds number flows. The next section provides the detailed derivation and is followed by the linear stability analysis of the system, showing that wavenumbers in a finite interval near zero are linearly unstable. Simulations of the full nonlinear equations with periodic boundary conditions in the axial direction show the emergence of finite amplitude steady traveling waves. We also carry out physical experiments on this system using a simple setup with safflower oil and fishing lines and capture images of the droplets that travel down the fiber. We show that our model can match the experimental results closely, with the best set of parameters obtained using machine learning, trained on a large set of simulation results with randomly chosen parameters near the physical range.

Model Derivation via Control Volume Analysis
In this section we derive our model for an axisymmetric liquid film flowing down an infinitely long cylindrical fiber. Our approach is based on a control volume analysis of the conservation of mass and momentum equations, in which the axial velocity is replaced by a mean velocity that is uniform in the cross section but varies with axial distance and time. Assuming such a plug-flow profile greatly simplifies the derivation. However, one of the key terms that relates the viscous drag force on the fluid by the fiber is actually treated more carefully to make it consistent with the laminar flow profile for fully-developed flow down the fiber. Even if the flow is truly closer to a plug flow-e.g., in the high-Reynolds number turbulent regime where mixing causes the profile to be more uniform-we can still account for a drag force exerted between the solid surface of the fiber and the flowing film, proportional to the flow velocity, with some constant empirical coefficient related to a thin boundary layer thickness. As such, we end up with two versions of the model, one appropriate for low Reynolds number laminar flow and the other better suited to the high Reynolds number regimes. The models will appear quite similar though the scaling and the functional relation between the mean velocity and film thickness will be different between the two. Before deriving the model, it helps to compare and contrast these two cases in more detail, in the simpler situation when the flows are fully developed.

Plug Flow
This case is simple to analyze. Consider a cylindrical fiber of radius R and a liquid film whose interface is at distance H from the fiber axis, resulting in a liquid film of thickness H − R. Suppose that the fluid is falling down the fiber under the influence of gravity at uniform speed U . At steady state (terminal draining velocity), the weight of any portion of the liquid between two axial locations is balanced by the drag force exerted by the solid surface of the fiber on the liquid. The weight of the liquid between two axial locations x 1 and x 2 , with ∆x = x 2 − x 1 , is given by ρgπ(H 2 − R 2 )∆x. If the shear stress at the fiber surface is denoted by τ rx , the drag force exerted on that portion of liquid would be 2πRτ rx ∆x. Based on a dimensional reasoning, the form of the shear stress could be assumed to be in which parameter ℓ is some quantity with units of length. It could be thought of as some measure of an extremely thin boundary layer thickness that might be separating the plug flow region with velocity U from the fiber surface on which a no-slip boundary condition would exist. Of course, we ignore the boundary layer region when assuming plug flow, but still account for the drag force that the fiber exerts on the liquid. By balancing the weight of the liquid with the drag force, we can obtain a relationship between the flow speed U and the film thickness H. The result is in which h = H/R is the ratio of liquid film radius to the fiber radius. If we assume parameter ℓ to be constant, the velocity scale can be chosen to be U o = ρgRℓ/2µ and the dimensionless draining velocity u = U/U o would be given by u = f (h) = h 2 − 1. We will compare this quadratic expression for the draining velocity as a function of h with the result for fully developed viscous flow obtained below. We will find that this function f (h) increases much more rapidly as h increases away from 1, as compared to the situation with viscous laminar flow.

Viscous Laminar Flow
For the case of fully-developed laminar flow down the fiber, the velocity profile u(r) can be obtained by integrating the axial component of the Navier-Stokes equation which reads µ r d dr (r du dr ) + ρg = 0 .
The boundary conditions are that u(R) = 0 (no slip on the fiber surface) and u ′ (H) = 0 (zero shear stress at the free surface). The resulting velocity profile is given by The mean velocity U can be calculated using the definition U = H R ru(r)dr/ H R rdr resulting in The shear stress at the fiber surface τ rx = µu ′ (R) can be expressed as before in the form but with length parameter ℓ now depending on h and given by As such, the main difference between the plug flow model and the viscous laminar flow one is that in the former, ℓ is treated as a constant, whereas in the latter, it depends on the film thickness. With the proportionality constant between the shear stress and the mean velocity being dependent on h, the functional form of the dependence of the mean draining velocity on film thickness is quite different. In particular, the dimensionless mean velocity, now scaled with velocity scale U 1 = 2ρgR 2 /µ, would be given by where function I(h) is given by Eq. (2.1); this can be compared to the result for plug flow, which was u(h) = f (h) = h 2 − 1. The mean velocity u varies much more as h increases away from 1 for the plug flow model than for the laminar flow case. Close to As such, for thin liquid films, the rate of increase of draining velocity with increasing film thickness is much stronger in the plug flow model than in the laminar flow one.

Control Volume Analysis
In order to derive the equations of motion for a falling film in which the film thickness varies with axial distance and time, i.e., H = H(x, t), we use a control volume approach as depicted in Figure 1. We assume the velocity in the film to be uniform in the crosssection (interpreted as the mean velocity in the laminar case), but allow the latter to vary with axial location and time as well: U = U (x, t). We consider a control volume consisting of the portion of the fluid between two axial locations x and x + ∆x, as shown in the figure. Denote the cross-sectional area of the fluid at any axial position and time The integral form of the conservation of mass in the region between x and x + ∆x reads equating the rate of change of mass to the rate at which mass enters the control volume at position x minus the rate at which it leaves at position x+∆x. Based on the intermediate value theorem from calculus, the left-hand side of this equation can be written as Dividing both sides of the equation by ∆x and taking the limit ∆x → 0 results in the equation for conservation of volume, as expected. Since A(x, t) = π(H 2 (x, t) − R 2 ), we can rewrite this equation as Moving on to the conservation of linear momentum in the axial direction, one can similarly equate the rate of change of total linear momentum in the control volume to the net rate at which momentum flows into the control volume plus the sum of the forces in the axial direction acting on the fluid in that volume. This equation takes the form The terms on the right-hand side of this equation have the following physical interpretations: The first two terms provide the net rate at which momentum enters the control volume across the two boundaries, the next term is the weight of the volume of fluid in the control volume, the next two capture the contribution from the pressure force acting on the two cross-sections, followed by the two terms that account for any viscous normal stress at those same cross-sections, the next term is the drag force exerted on the fluid by the solid surface of the fiber, and finally, the last two terms capture the effect of surface tension acting on the perimeter of the free surface (since surface tension is tangent to the interface, to project it onto the axial direction, we need the cosine of the angle that the tangent vector makes with the axial direction in those terms). Points ξ, ξ ′ and ξ ′′ are somewhere in the interval [x, x + ∆x]; their precise location becomes irrelevant as ∆x tends to zero. Upon dividing this equation by ∆x and taking the limit ∆x → 0, we get the differential equation Using the conservation of volume equation, the left-hand side of the last equation can be simplified to ρA(∂U/∂t + U ∂U/∂x). Also, we substitute µU/ℓ for the shear stress τ rx and 2µ∂U/∂x for the normal viscous stress τ xx . Upon dividing the entire equation by the cross-sectional area A(x, t) we thus obtain In this equation, the cross-sectional area is given by A(x, t) = π(H 2 (x, t) − R 2 ), and since tan(θ) = ∂H/∂x (the slope of the free surface), the cosine of that angle is given by cos(θ) = 1/ 1 + H 2 x in which subscript refers to a partial derivative. The pressure within the film, p(x, t), is taken to be uniform in the cross section and related by the Young-Laplace equation to the curvature of the free surface, namely p(x, t) = σκ(x, t), in which σ is the surface tension and the curvature κ is given in this geometry by x ) 3/2 , with subscripts referring to partial derivatives. Note that ordinarily the pressure in the fluid would be written as p = p o +σκ in which p o is the constant pressure in the air outside the interface. However, in calculating the force on the control volume, the contribution of the force due to p o acting all around the control volume (including on the curved free surface) integrates to zero, so that constant part of the pressure is omitted.
The pressure term in the momentum equation can be written as a sum of two terms: Interestingly, the second term on the right-hand side is exactly equal to the surface tension term on the right-hand side of the momentum equation, namely the term 2πσ A ∂(H cos θ) ∂x , so those two terms cancel each other leaving simply ∂p/∂x on the left-hand side of the momentum equation. The above cancellation is a consequence of a relationship that appears to be purely geometrical, involving the curvature κ and the rates of change of area and the perimeter multiplied by the cosine factor, namely: κ∂A/∂x = 2π∂(H cos θ)/∂x in which cos θ = (1 + H 2 x ) −1/2 . After this simplification, the momentum equation further divided by density ρ becomes , upon choosing the velocity scale U o = gRℓ/2ν and defining the dimensionless velocity u = U/U 0 , and upon scaling time with U o /g, i.e., witht = gt/U o , the first term on the left-hand side and the first two terms on the right-hand side would yield a dimensionless equation of the form .
Such an equation would hold if all x-derivatives were absent. It would suggest that for a given dimensionless film thickness h, the velocity u of the film would relax exponentially in time to its terminal velocity f (h) = h 2 − 1 with a relaxation time of order one in dimensionless timet. Carrying the scaling further by nondimensionalizing the axial distance x and curvature κ with the fiber radius R so thatx = x/R andκ = Rκ, we obtain the fully nondimensional form of the axial momentum equation which, upon dropping the hats for clarity, reads (2.5) in which subscripts represent partial derivatives. The dimensionless curvature appearing in this equation is given by (2.6) Three dimensionless parameters, called a, b and c, also appear in this equation, given respectively by: In its original form, parameter a is seen to be the square of the Froude number and parameter b is the reciprocal of the Bond or Eötvös number. Parameter c is the ratio of the characteristic boundary layer thickness to the fiber radius.
Using the same scaling, the dimensionless form of the conservation of volume equation takes the form: (2.8) 2.3. Laminar flow case For the fully-developed laminar flow model, the above derivation proceeds similarly though in an approximate sense. The dimensionless mean velocity is now given by u(h) = U/U 1 = f 1 (h) as shown in Eq. (2.2) with U 1 = 2ρgR 2 /µ. This modifies the scaling and some of the dimensionless parameters in the model. Also, the balance of the gravity force and the radial derivatives in the viscous term produces the term g[1 − u/f 1 (h)] on the right-hand side of the momentum equation where f 1 (h) is defined within Eq. (2.2). With those changes, and upon scaling time with U 1 /g and length with R, the dimensionless momentum equation takes a similar form to the one for the plug-flow model, i.e., with parameter b staying the same, while dimensionless parameters a 1 and c 1 differ from their earlier counterparts, now being given by (2.10) The key assumption underlying this derivation is that the velocity in the axial direction can be replaced within most of the terms in the momentum equation by its mean over the liquid cross section. In particular, the square of the velocity on the left-hand side of the equation is replaced by the square of the mean velocity. Also, the pressure field is still assumed to be uniform in the cross section. Despite the approximate nature of this model, we shall see that it is successful in matching the experimental results fairly well with parameter values that are quite close to their theoretical values. The dimensionless equation for conservation of mass (or volume) looks the same as Eq. (2.8) but with parameter a replaced by a 1 in this case: (2.11)

Summary
To summarize, our one-dimensional two-equation model for an axisymmetric liquid film falling down a vertical fiber consists of the equations for the conservation of mass (2.8) and axial momentum (2.5) for plug flow, or the corresponding pair (2.11) and (2.9) for laminar flow. The dependent variables are the dimensionless film radius h(x, t) and the mean axial velocity u(x, t). The dimensionless film thickness is given by h(x, t) − 1. The parameters for the plug flow case are given in Eq. (2.7) and those for laminar flow in Eq. (2.10). Dimensionless curvature is given by the expression in Eq. (2.6). We have used both of the forms given in that equation successfully in our numerical simulations. For the plug flow model, we have the function f (h) = h 2 − 1; for the laminar flow model, it is replaced by f 1 (h) = 4h 4 ln(h) − 3h 4 + 4h 2 − 1 / 16(h 2 − 1) . When performing numerical simulations, we choose a domain of dimensionless length L in space, so that x ∈ [0, L], and solve the equations up to a chosen final dimensionless time of T , so that t ∈ [0, T ]. We apply periodic boundary conditions in x. Parameters a, b, c (or a 1 , b, c 1 ) are all positive and represent the effects of inertia, surface tension and axial viscous diffusion respectively.

Linear Stability Analysis
In this section, we conduct a linear stability analysis about constant solutions of the system (2.8)-(2.5) (plug flow) or (2.11)-(2.9) (laminar flow). Note that any constant h 0 and u 0 that satisfy u 0 = f (h 0 ) (or u 0 = f 1 (h 0 )) is a solution of the system. We define the perturbed solution in the form below: where h 0 and u 0 are constants. Small parameter ǫ is introduced for bookkeeping purposes only. At O(ǫ) we derive a linearized system for the leading perturbations as for the plug flow case. For laminar flow, f (h 0 ) and f ′ (h 0 ) should be replaced by f 1 (h 0 ) and f ′ 1 (h 0 ), and a and c are replaced by a 1 and c 1 . These being linear constant-coefficient equations, one can seek exponential solutions for h 1 and u 1 in the form: whereĤ andÛ are complex amplitudes, ℜ denotes the real part, k is the wavenumber (assumed real) and α is the growth rate. Upon substitution of these exponential forms into the linearized equations, we obtain the linear system The above form is for the plug flow case; for laminar flow we replace f , a and c with f 1 , a 1 and c 1 .
To have a non-trivial solution of the system (3.1) forĤ andÛ , the determinant of the coefficient matrix needs to be zero. This provides a quadratic equation for the growth rate α, which can be solved analytically but involves lengthy expressions that are not displayed here. Given h 0 , u 0 = f (h 0 ), and parameters a, b and c, we obtain the two roots of the quadratic equation α 1 and α 2 (which are complex valued) as functions of the wavenumber k. The procedure is identical for the laminar flow model, although function f 1 (h) is involved and parameters a 1 and c 1 = 4. Figure 2 provides plots of the real parts of the two growth rates versus wavenumber k. These are for the laminar flow model with parameters: h 0 = 2.5, a 1 = 1 and b = 11. These values are chosen since they are quite close to those for the experiments described in Sections 5 and 6. Note that if the real part of either growth rate is positive, waves of those wavenumber grow and the system is linearly unstable. In the figure, we see that one of the roots does indeed have a positive real part over a range of wavenumbers k ∈ (0, k max ), with a maximal growth rate occurring for some wavenumber in that interval. We thus see that uniform solutions are unstable to perturbations of small wavenumber or long wavelength. Figure 3 is obtained for the plug-flow model with the set of parameters indicated in the caption, approximating the hypothetical high Reynolds number case discussed in Section 7. We still find that one of the growth rates has a positive real part over a range of wavenumbers near zero. However, in this case the actual rate of growth is much higher. As a result, numerical simulations of the plug-flow model are more challenging since the function f (h) varies a lot more than f 1 (h) and the growth rate of the instabilities is also quite a bit higher than for the laminar flow model. The stability result was checked against numerical simulations. We did a comparison between two simulations, one using the full nonlinear model and the other using the linearized one to see the effects of nonlinearity on the film thickness evolution. In the simulation, we took an initial profile h(x, 0) = h 0 + ǫ sin(kx) with a small ǫ and the domain length was chosen as L = 2π/k. We chose the wavenumber k as the most unstable wavenumber for h 0 . We found that the rate of increase of the maximum film height does follow the simple exponential function (the logarithm of the perturbation growing linearly in time) for the linearized model, while for the full nonlinear model the growth slows down as time increases as a result of nonlinear interactions.

Simulations
Using COMSOL Multiphysics we carried out many simulations of both the plug flow system (2.8)-(2.5) and the laminar flow one (2.11)-(2.9) for various sets of parameters. For the results reported in this section, we took the domain x ∈ [0, L] with L = 20 and assumed periodic boundary conditions in x. For the initial condition, we took h(x, 0) = h 0 + 0.1 sin(2πx/L) for the film profile and u(x, 0) = f (h(x, 0)) or u(x, 0) = f 1 (h(x, 0)) for the initial velocity. We integrated the equations to a final time ranging from several hundreds to several thousands until a steady traveling wave profile was obtained. We found that for smaller values of parameter a it takes longer to reach a steady traveling wave shape. Focusing on the shape of the traveling wave profile, we took the final steady shape and centered its peak in the middle of the interval in order to be able to make the following comparison plots. We see that the plug flow model leads to a more pronounced peak as compared to the laminar flow model. On the right panel of the same figure, we vary the parameter h 0 while keeping the other ones fixed in the laminar flow model. As h 0 increases, not only does one obtain a more prominent peak, the speed of the resulting traveling wave also increases substantially. In Section 6 we explain how to obtain the speed of the traveling wave from a single snapshot of the steady height and velocity profiles. Parameters a 1 and b also affect the shape of the traveling wave, but not as significantly as the film radius h 0 . As seen in Figure 5, varying a 1 (left panel) or b (right panel) does impact the shape of the height profile and the speed of the traveling wave, but much less so than parameter h 0 .

Experiment setup
The experimental setup is shown below in Figure 6. Our working fluid is safflower oil with parameter values: density ρ = 0.928 g/cm 3 , absolute viscosity µ = 0.0654 Pa·s, surface tension σ = 0.025 N/m; viscosity was extrapolated to our working temperature of 13 • C starting with the data from Diamante & Lan (2014); density and surface tension data were found online at https://cameochemicals.noaa.gov/chris/OSF.pdf. Safflower oil was placed in a cup with a hole drilled in the center. A nylon fishing line was passed through the hole and tied to a hanging weight at the bottom in order to maintain a vertical line. The line was threaded through a wooden rod that was placed horizontally across the top of the cup to hold the fiber in place. Our setup was modeled after those presented in other papers examining droplet flow on a vertical fiber, e.g., Kliakhandler et al. (2001) and Craster & Matar (2006). One main difference in our setup is the use of larger diameter fishing lines on which different regimes of coating flows were observed. Cups with hole sizes ranging from 1.6 to 2.4 mm were used in order to control the flow rates, and both 1 and 2 mm diameter fishing lines were used. Snapshots of the flow were photographed with a green screen background using a Canon EOS Rebel T6s at a shutter speed of 1/4000 s. Video was also captured in order to determine droplets velocity. Figure 7: a) nearly uniform coating flow over a 2-mm diameter fishing line with a 2.4-mm hole; b) and c) non-uniformly spaced droplets on a 1-mm diameter fishing line with the diameter of the hole being 1.6 mm; d) two sets of doublets on a 1-mm diameter fishing line with a hole diameter of 1.6 mm; e) uniformly spaced four-droplet train on a 1-mm diameter fishing line with a hole diameter of 2 mm; f) uniformly distributed train of four droplets on a 1-mm diameter fishing line with the hole diameter being 2.4 mm (the reference ruler on the left shows 1 mm tick marks, with 1/2 mm ones at the very top).

Observations
We lubricated our fishing lines with the oil before running the experiments since on dry fishing lines we would see a very strong instability near the advancing side of the leading drop, causing tiny droplets to be scattered in all directions while propagating down the line.
For the 2-mm diameter fishing line through the 2.4-mm size hole we observed a nearly uniform coating flow. We could see a very slight waviness of this film near the edge so we were certain that some very small amplitude waves were present; see a) in Figure 7. This almost uniform coating flow behavior was not mentioned in prior experimental publications such as Kliakhandler et al. (2001) and Craster & Matar (2006), most likely due to the thinner lines they used. For the 1-mm diameter fishing line, with hole diameters ranging from 1.6 to 2.4 mm, different distributions of droplets flowing down the line were observed. For the 1.6-mm diameter hole with the 1-mm diameter fishing line we observed non-evenly spaced trains of droplets: see b) and c) in Figure 7, and we also obtained interesting doublet configurations: e.g., d) in Figure 7. The reason the distances are so different in this regime is that droplets interact with each other as they do not have the same size and speed. For the 2-mm hole we observed uniformly spaced droplets, see e) in Figure 7, and for the 2.4-mm hole we also observed uniformly spaced droplets that were bigger than the previous case with slightly longer periods: see f) in Figure 7. The distributions of droplets we observed on the 1-mm line are qualitatively similar to ones reported in Kliakhandler et al. (2001) and Craster & Matar (2006): in Figure 7, b), c) and d) correspond to the Convective Instability Regime, and e) and f) correspond to the In all experiments we noticed that a steady "ring" of fluid always developed around the cup exit hole and new droplets always detached from this ring without changing the shape of the ring.

Application of Machine Learning to Model Validation
Using the approximate physical properties of the oil mentioned earlier (ρ = 928 kg/m 3 , µ = 0.0654 Pa·s and σ = 0.025 N/m), with g = 9.8 m/s 2 and for the 1-mm diameter fiber whose radius is R = 0.5 mm, for a nominal film thickness H − R = 0.75 mm (estimated visually), the velocity scale U 1 = 2gR 2 /ν ends up being 6.95 cm/s and the dimensionless parameters for the laminar flow model turn out to be a 1 = 0.987, b = 11.0 and h 0 = 2.5. Also, if we define the Reynolds number as Re=U 1 R/ν we get Re= 0.493, which is low enough that the laminar flow model should be reasonable. However, since the physical properties are approximate and the film thickness is not exactly measured, we apply Machine Learning to find the closest set of parameters that match the model predictions to the experimental observations.
Using COMSOL we created a labeled set of 182 data files in which a final snapshot of the steady traveling wave profile h(x) and mean velocity u(x) were saved as a function of discretized x at 961 equally spaced nodal points including the end points on the interval x ∈ [0, L] with L = 20. In those data sets, random values of parameters a 1 , b and h 0 were chosen in the neighborhood of the original estimates. Note that parameter c 1 is fixed at c 1 = 4 in the laminar flow model. We took uniformly distributed random values of h 0 ∈ [2.2, 2.8], of a 1 ∈ [0.1, 1.1], and of b ∈ [10, 13]. In the simulations, we started with the initial condition h(x, 0) = h 0 + 0.1 sin(2πx/L) and integrated to a long enough time that a steady traveling wave profile was reached. The label of each data file included values of the random parameters that generated it. While all three parameters varied at random in these simulations, the set of profiles obtained were roughly similar to those displayed in Figures 4 and 5 in which two of the parameters were held fixed while the third varied.
We then trained our learning algorithm on this data set so that we would be able to "predict" the values of the dimensionless parameters if we were given a certain discretized profile h(x). That way, we could discretize the experimentally observed profile and find out the set of parameter values that would generate that profile. For the supervised learning we used the so-called Gradient Boosting Regressor (GBR) based on an ensemble of decision trees. We applied the Holdout Method for cross validation by taking 20% of our data as a test set.
To compare our numerical simulation result (steady traveling wave h function) with an experimental shape of a droplet we zoomed in on the droplet photo and traced the edge by manually inserting marker points using Mathematica (red points in Figure 8) to discretize the profile. We then used interpolation of the captured coordinates of the marker points to represent this profile in our numerical 961-point data format and used these experimental data as an input for our trained GBR model to obtain the following predictions for the parameter values: a = 0.95 and b = 12.32. The value of h 0 could be obtained readily through a conservation of volume constraint. Namely, since using the experimental profile on the right-hand side yielded h 0 = 2.29. We used the initial perturbation sin 2πx/L to speed up the convergence of simulations to a steady traveling wave on the computational domain x ∈ [0, L] with L = 20 as the wavenumber k = L/(2π) is close to the most unstable mode in the linearized about h 0 model. We compare the experimental results and the numerical simulations by moving the maximum height for both to the same point of the grid. The bottom right panel in Figure 8 displays the two curves and the good agreement between them.
In the process of training the GBR (see Figure 9) we discovered that different tailored data sets could be used to improve the accuracy for the prediction of coefficients a 1 and b. The optimal approach for obtaining the coefficient a 1 is to use as training data the first 6 coefficients of the Fourier series of the periodic traveling wave h(x), and for the coefficient b to use as training data the first 6 coefficients of the Fourier series of the function 1/(h(x) 2 − 1). This difference was a good indication that the shape of the traveling wave h(x) is mostly defined by the value of the coefficient a 1 and the mean velocity u(x) is mostly governed by the value of the coefficient b. The derivation below explains the relation between the speed u(x) and the reciprocal of the area 1/(h(x) 2 − 1).
Using the traveling wave ansatz that h(x, t) and u(x, t) only depend upon the dimensionless traveling wave coordinate z = x − V t, with u = u(z) and h = h(z), the conservation of mass equation (2.11) results in the ordinary differential equation: −V (h 2 ) ′ + a 1 (u(h 2 − 1)) ′ = 0 whose solution after some algebraic manipulation provides Here, primes denote derivatives with respect to z. So u(z) does depend directly on the 1/(h(z) 2 − 1). This computations also suggest a simple method of finding the traveling wave velocity V from just one snapshot of the profiles h(x) and u(x) at some late time.
One can find the dimensionless traveling wave speed V by just fitting the profile of u(x) to a function of the form u(x) = C 1 + C 2 1 h(x) 2 −1 and using the best fit coefficient C 1 to define V = a 1 C 1 , see Figure 10. To obtain the dimensional traveling wave speed from the coefficient C 1 , we multiply the latter by velocity scale U 1 .

Discussion
One of the main modeling challenges related to droplets sliding down a fiber is the mismatch between the experimentally observed velocity of the drops and that obtained from numerical simulations based on a model. In the models by Kliakhandler et al. (2001) and Craster & Matar (2006) Figure 10: The plots on the left illustrate the perfect fit obtained when u(x) is fit to a function of the form C 1 + C 2 /(h(x) 2 − 1). The blue curve shows 1/(h 2 − 1) while the other curve shows u(x) and its best fit that appear superposed. The plots on the right show the height and velocity profiles at two different times and confirm that the traveling wave speed V is matched with high accuracy.
traveling wave droplet velocity and the experiments. Our laminar flow model provides closer values of the traveling wave velocity to the experiments for smaller values of h 0 . However, for the value h 0 = 2.29 obtained above, the predicted velocity is about 6 cm/s, whereas the experimentally observed droplet velocity seemed to be about half of that. The other challenging part in modeling is to obtain clear qualitative transition criteria between the various regimes observed in experiments. Another still unsolved puzzle is the influence of the size or geometry of the hole (i.e., the source) on the distribution of droplets and their dynamics while flow rate of the fluid is kept constant. While in this paper we focused most of our attention on the laminar flow model appropriate for low Reynolds numbers, our plug flow model should apply to well-mixed possibly turbulent flows with a velocity profile closer to plug flow. Just to see whether this is physically feasible, consider a hypothetical system with the following assumed parameters: suppose the fiber radius R is 2 mm and the film thickness T = H − R is about the same size as the fiber radius. Take the working fluid to be water whose viscosity is much less than the oil. We can approximate the corresponding Reynolds number for a steady state flow of this type. The mean wall shear stress τ is expressed in terms of the Darcy-Weisbach friction factor f D and average fluid velocity U as The force balance between the drag force from the wall and gravity gives us 2πRτ = ρgπ((R + T ) 2 − R 2 ) which can be solved to obtain τ = 78.4 Pa (we take the density of water to be ρ = 1000 kg/m 3 and its viscosity to be µ = 0.001 kg/(m s)). The Colebrook-White correlation for a smooth surface relates the friction factor to the Reynolds number by We substitute the expression for f D in terms of τ and U , and Re = ρUT µ and obtain a transcendental equation for U , whose solution yields the mean velocity and corresponding Reynolds number as U = 11.31m/s, Re = 4.5 × 10 4 .
The result shows that under some practical assumptions, the film flow on fiber could be in a turbulent regime. Under the assumptions above, the parameters a, b, c in our model would have values a ≈ 102, b ≈ 1.84, c ≈ 0.072 . Exploration of such turbulent regimes and their experimental investigations are left for future work.