Next Article in Journal / Special Issue
Patterns and Their Large-Scale Distortions in Marangoni Convection with Insoluble Surfactant
Previous Article in Journal
Efficient Wildland Fire Simulation via Nonlinear Model Order Reduction
Previous Article in Special Issue
Approximate Analytical Analysis of Unsteady MHD Mixed Flow of Non-Newtonian Hybrid Nanofluid over a Stretching Surface
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Liquid Films Falling Down a Vertical Fiber: Modeling, Simulations and Experiments

1
Institute of Mathematical Sciences, Claremont Graduate University, Claremont, CA 91711, USA
2
College of Engineering, University of California, Berkeley, CA 94720, USA
*
Author to whom correspondence should be addressed.
Fluids 2021, 6(8), 281; https://doi.org/10.3390/fluids6080281
Submission received: 14 July 2021 / Revised: 5 August 2021 / Accepted: 6 August 2021 / Published: 12 August 2021
(This article belongs to the Special Issue Thin Liquid Films: From Theory to Applications)

Abstract

:
We provide a new framework for analyzing the flow of an axisymmetric liquid film flowing down a vertical fiber, applicable to fiber coating flows and those in similar geometries in heat exchangers, water treatment, and desalination processes. The problem considered is that of a viscous liquid film falling under the influence of gravity and surface tension on a solid cylindrical fiber. Our approach is different from existing ones in that we derive our mathematical model by using a control-volume approach to express the conservation of mass and axial momentum in simple and intuitively appealing forms, resulting in a pair of equations that are reminiscent of the Saint-Venant shallow-water equations. Two versions of the model are obtained, one assuming a plug-flow velocity profile with a linear drag force expression, and the other using the fully-developed laminar velocity profile for a locally uniform film to approximate the drag. These can, respectively, model high- and low-Reynolds number regimes of flow. Linear stability analyses and fully nonlinear numerical simulations are presented that show the emergence of traveling wave solutions representing chains of identical droplets falling down the fiber. Physical experiments with safflower oil on a fishing line are also undertaken and match the theoretical predictions from the laminar flow model well when machine learning methods are used to estimate the parameters.

Graphical Abstract

1. Introduction

1.1. Background

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 [1,2,3,4,5]. 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 [6,7,8], and particle capturing systems [9], 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 [10,11,12,13]: (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 [14] proposed a weakly nonlinear thin-film equation to calculate the evolution of film thickness θ ( x , t ) (in this paper we use variables H ( x , t ) and R to denote the film and fiber radii, respectively, with θ ( x , t ) = H ( x , t ) R ; in dimensionless form h = H / R and θ / R = h 1 ; x is the axial coordinate and t is time) and capture both stabilizing and destabilizing effects of the surface tension in the dynamics. This evolution equation was further studied by Kalliadasis and Chang [10], Chang and Demekhin [15], Marzuola et al. [16], and Ji et al. [17]. Craster and Matar [18] 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. [19] extended the thin film model to consider thick layers of viscous fluid by introducing fully nonlinear curvature terms. Marangoni effects resulting in different types of wave instabilities for liquid droplets sliding down a wire were examined asymptotically and numerically in Chinju et al. [20], Duprat et al. [21], Duprat et al. [22], and Ding et al. [23]. Recently Ji et al. [24] 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 [25] 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 Ji et al. [24] 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 [26,27,28,29]. These equations incorporate inertia effects and streamwise viscous diffusion based on the integral boundary-layer approach. Recently, Ji et al. [13] 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 and Ding [30] 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. Existence and stability of traveling waves in falling down liquid thin film models were studied in Alekseenko et al. [31], Bertozzi et al. [32], Bertozzi and Shearer [33], and Bertozzi et al. [34].

1.2. Problem Overview and Approach

The key problem is thus to obtain an intuitive model that can capture the physics of flow down a fiber, even at high flow rates when the flow might be of high Reynolds number and possibly turbulent, accounting for the effects of inertia, gravity, surface tension and viscous drag on the film. Having such a model should lead to a better understanding of fiber coating flows or those that occur in heat and mass exchangers. For instance, it has recently been proposed to capture water drops from humid air through condensation of drops on cotton fibers [8]. A model that can be readily extended to include condensation effects would thus be valuable. This is relatively straightforward to do via the control-volume approach which we are providing.
Existing models have had some success in predicting some of the experimental observations, though discrepancies still exist and the some terms are added to the models in an ad hoc manner. The one-equation lubrication models reviewed above have the shortcoming that they do not account for the inertial effects and are appropriate for low Reynolds number regimes. The two-equation models that have been used previously do include effects of inertia, but they are usually based on assuming a particular local velocity profile (our model makes a similar assumption) and to date they have only been applied to laminar fully-developed profiles, while we also include plug flow profiles that better model high Reynolds number flows. Furthermore, the formulations that use the volumetric flow rate rather than velocity as the dependent variable yield equations that have many terms that are hard to interpret physically. In contrast, our approach yields an axial momentum equation that is very intuitive and includes all the terms one would expect, similar to the Saint-Venant shallow-water equations.
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.

2. 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 as depicted in Figure 1. 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.

2.1. Fully-Developed Flow

2.1.1. 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 . Here, ρ is the density of the liquid and g is the gravitational acceleration. If the shear stress at the fiber surface is denoted by τ r x , the drag force exerted on that portion of liquid would be 2 π R τ r x Δ x . Based on a dimensional reasoning, the form of the shear stress could be assumed to be
τ r x = μ U ,
in which parameter is some quantity with units of length and μ is the viscosity of the liquid. 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 radius H. The result is
U = ρ g R 2 μ ( h 2 1 ) ,
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 = ρ g R / 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.

2.1.2. Viscous Laminar Flow

For the case of fully-developed laminar flow down the fiber, the velocity profile u ( r ) in a cylindrical coordinate system can be obtained by integrating the axial component of the Navier–Stokes equation which reads
μ r d d r ( r d u d r ) + ρ 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
u ( r ) = ρ g R 2 4 μ 1 ( r R ) 2 + 2 ( H R ) 2 ln ( r R ) .
The mean velocity U can be calculated using the definition U = R H r u ( r ) d r / R H r d r resulting in
U = 2 ρ g R 2 μ I ( h ) ( h 2 1 )
with h = H / R as before and
I ( h ) = 1 16 4 h 4 ln ( h ) 3 h 4 + 4 h 2 1 .
The shear stress at the fiber surface τ r x = μ u ( R ) can be expressed as before in the form
τ r x = μ U ( h )
but with length parameter now depending on h and given by
( h ) R = 4 I ( h ) ( h 2 1 ) 2 .
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 ρ g R 2 / μ , would be given by
u ( h ) = U U 1 = I ( h ) h 2 1 = f 1 ( h ) ,
where function I ( h ) is given by Equation (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 h = 1 , f ( h ) 2 ( h 1 ) whereas f 1 ( h ) ( h 1 ) 2 . 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.

2.2. 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 shown in Figure 1. We assume the velocity in the film to be uniform in the cross-section (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 ) . Here x denotes the axial coordinate along the fiber center-line and t is time. 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 x by A ( x , t ) = π ( H 2 ( x , t ) R 2 ) .
The integral form of the conservation of mass in the region between x and x + Δ x reads
d d t x x + Δ x ρ A ( x , t ) d x = ρ A U | x ρ A U | x + Δ x ,
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
x x + Δ x ρ A t ( x , t ) d x = ρ A t ( ξ , t ) Δ x ,
where ξ is somewhere in the interval [ x , x + Δ x ] . Dividing both sides of the equation by Δ x and taking the limit Δ x 0 results in the equation
A t + ( U A ) x = 0
for conservation of volume, as expected. Since A ( x , t ) = π ( H 2 ( x , t ) R 2 ) , we can rewrite this equation as
2 H H t + ( U ( H 2 R 2 ) ) x = 0 .
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
ρ Δ x t ( A U ) | ξ = ρ ( A U 2 ) | x ρ ( A U 2 ) | x + Δ x + ρ g Δ x A | ξ + ( p A ) | x ( p A ) | x + Δ x + ( A τ x x ) | x + Δ x ( A τ x x ) | x 2 π R Δ x τ r x | ξ + 2 π σ ( H cos ( θ ) ) | x + Δ x 2 π σ ( H cos ( θ ) ) | x .
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
ρ ( A U ) t + ρ ( A U 2 ) x = ρ g A ( p A ) x + ( τ x x A ) x 2 π R τ r x + 2 π σ ( H cos θ ) x .
Using the conservation of volume equation, the left-hand side of the last equation can be simplified to ρ A ( U / t + U U / x ) . Furthermore, we substitute μ U / for the shear stress τ r x and 2 μ U / x for the normal viscous stress τ x x . Upon dividing the entire equation by the cross-sectional area A ( x , t ) we thus obtain
ρ ( U t + U U x ) + 1 A ( p A ) x = ρ g + 2 μ A x ( A U x ) 2 π μ R U A + 2 π σ A ( H cos θ ) x .
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 x 2 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 κ (more properly called twice the mean curvature) is given in this geometry by
κ ( x , t ) = ( 1 + H x 2 H H x x ) H ( 1 + H x 2 ) 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. To obtain the expression for curvature, it is easiest to define the function S ( r , x , t ) = r H ( x , t ) whose zero surface S = 0 defines the interface in cylindrical coordinates. The outward unit normal to the interface is given by n ^ = S / | S | and twice the mean curvature can then be obtained by taking the divergence, κ = · n ^ , and evaluating it at the interface.
The pressure term in the momentum equation can be written as a sum of two terms:
1 A ( p A ) x = p x + σ κ 1 A A x .
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 x 2 ) 1 / 2 . Note that in our approach, the effect of surface tension has been included fully, both in the way it affects the pressure within the film as a consequence of the Young–Laplace equation, and as an extra force in the axial direction where the sides of the control volume cut across the interface. The combined effects, together with the exact expression for the curvature of the interface, give rise to an overall surface tension term ( σ κ / ρ ) / x below. After this simplification, the momentum equation further divided by density ρ becomes
U t + x 1 2 U 2 + σ κ ρ = g 2 π ν R U A + 2 ν A x A U x .
Here ν = μ / ρ is the kinematic viscosity of the fluid. Since A = π ( H 2 R 2 ) = π R 2 ( h 2 1 ) = π R 2 f ( h ) , upon choosing the velocity scale U o = g R / 2 ν and defining the dimensionless velocity u = U / U 0 , and upon scaling time with U o / g , i.e., with t ^ = g t / 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
u t ^ = 1 u f ( h ) .
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 time t ^ .
Carrying the scaling further by nondimensionalizing the axial distance x and curvature κ with the fiber radius R so that x ^ = x / R and κ ^ = R κ , we obtain the fully nondimensional form of the axial momentum equation which, upon dropping the hats for clarity, reads
u t + ( a u 2 / 2 + b κ ) x = [ 1 u / f ( h ) ] + c ( h 2 1 ) 1 [ ( h 2 1 ) u x ] x ,
in which subscripts represent partial derivatives. The dimensionless curvature appearing in this equation is given by
κ = 1 h ( 1 + h x 2 ) 1 / 2 x h x ( 1 + h x 2 ) 1 / 2 = ( 1 + h x 2 h h x x ) h ( 1 + h x 2 ) 3 / 2 .
Three dimensionless parameters, called a, b and c, also appear in this equation, given, respectively, by:
a = Fr 2 = U o 2 R g = g R 2 4 ν 2 , b = 1 Bo = σ ρ R 2 g , c = 2 ν U o R 2 g = R .
Parameter a is seen to be the square of the Froude number, Fr = U o / R g , and parameter b is the reciprocal of the Bond or Eötvös number, Bo = ρ g R 2 / σ . Parameter c is the ratio of axial viscous to gravitational forces and reduces to 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 h h t + a [ u ( h 2 1 ) ] x = 0 .

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 Equation (2) with U 1 = 2 ρ g R 2 / μ . This modifies the scaling and some of the dimensionless parameters in the model. Furthermore, 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 Equation (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.,
u t + ( a 1 u 2 / 2 + b κ ) x = [ 1 u / f 1 ( h ) ] + c 1 ( h 2 1 ) 1 [ ( h 2 1 ) u x ] x ,
with parameter b staying the same, while dimensionless parameters a 1 and c 1 differ from their earlier counterparts, now being given by
a 1 = U 1 2 R g = 4 g R 3 ν 2 , c 1 = 2 ν U 1 g R 2 = 4 .
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. Furthermore, 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 Equation (9) but with parameter a replaced by a 1 in this case:
2 h h t + a 1 [ u ( h 2 1 ) ] x = 0 .

2.4. 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 (9) and axial momentum (6) for plug flow, or the corresponding pair (12) and (10) 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 Equation (8) and those for laminar flow in Equation (11). Dimensionless curvature is given by the expression in Equation (7). 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 ) = 4 h 4 ln ( h ) 3 h 4 + 4 h 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.

2.5. Relation to Previous Two-Equation Models

To put our new two-equation model into context, we relate it to the two-equation model of Ruyer-Quil et al. [27] who, in turn, have related their model to the ones that had been derived previously. For this purpose, we start with the dimensional versions of the two conservation equations, Equations (3) and (5). The latter can be written in the form:
U t + x U 2 2 = g 1 U F ( H ) σ ρ κ x + 2 ν A x A U x .
Function F ( H ) provides the correct mean velocity expression for a given film radius H; it can also be expressed in terms of the film thickness. Ruyer-Quil et al. [27] express their two-equation model in terms of the flow rate Q ( x , t ) and film thickness θ ( x , t ) , defined by:
θ ( x , t ) = H ( x , t ) R , Q ( x , t ) = U ( x , t ) A ( x , t ) ,
in which A = π ( H 2 R 2 ) is the cross-sectional area of the fluid and U ( x , t ) is the mean axial fluid velocity, which is being used in our model. Actually, they take the flow rate per perimeter of the fiber, q ( x , t ) = Q ( x , t ) / ( 2 π R ) , which we shall ultimately use below. With these changes of variables, the conservation of mass equation Equation (3), namely A t + Q x = 0 , can be written as:
( 1 + θ / R ) θ t + q x = 0 ,
which is identical to that of Ruyer-Quil et al. [27]. Subscripts refer to partial derivatives here.
To transform the conservation of momentum Equation (13), we first replace U with Q / A on its left-hand side and use the conservation of mass equation to help simplify it and obtain: U t + ( U 2 / 2 ) x = A 1 Q t + A 2 ( Q 2 ) x Q 2 A 3 A x . Next, we expand the axial diffusion term after substiting U = Q / A to get: A 1 ( A U x ) x = A 1 [ Q x x Q x ( ln A ) x Q ( ln A ) x x ] . At this point, to make the comparison easier, we focus on the limit when the film thickness θ is small compared to the fiber radius R, i.e., θ / R 1 , which is one of the limits considered in Ruyer-Quil et al. [27]. In that limit A 2 π R θ to leading order, ( ln A ) x θ 1 θ x and ( ln A ) x x θ 1 θ x x θ 2 ( θ x ) 2 . Making these substitutions in the momentum equation, using the variable q = Q / 2 π R and multiplying the entire equation by θ results in:
q t = g [ θ q / F ( θ ) ] ( σ θ / ρ ) κ x 2 q θ 1 q x + q 2 θ 2 θ x + 2 ν ( q x x θ 1 q x θ x q θ 1 θ x x + q θ 2 θ x 2 ) .
In this form, it becomes apparent that all the terms in the momentum equation of Ruyer-Quil et al. [27] are reproduced through this transformation. The coefficients in front of these terms do have the same signs as well, though their values (after nondimensionalization) are somewhat different while being comparable in magnitude.

3. Linear Stability Analysis

In this section, we conduct a linear stability analysis about constant solutions of the system (6)–(9) (plug flow) or (12)–(10) (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:
h ( x , t ) = h 0 + ϵ h 1 ( x , t ) + O ( ϵ 2 ) u ( x , t ) = u 0 + ϵ u 1 ( x , t ) + O ( ϵ 2 )
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
u 1 t + x a u 0 u 1 b ( h 1 h 0 2 + 2 h 1 x 2 ) = f ( h 0 ) h 1 f ( h 0 ) u 1 f ( h 0 ) + c 2 u 1 x 2 h 1 t + x u 0 h 1 + 1 2 u 1 ( h 0 1 h 0 ) = 0
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:
h 1 = { H ^ e i k x + α t } u 1 = { U ^ e i k x + α t }
where H ^ and U ^ 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
i ( 2 a k u 0 ) + 2 α i k a ( h 0 1 h 0 ) i ( k b h 0 2 + b k 3 ) f ( h 0 ) f ( h 0 ) i ( a k u 0 ) + α + c k 2 + 1 f ( h 0 ) H ^ U ^ = 0 .
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 (17) for H ^ and U ^ , 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 Section 5 and Section 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 ( k x ) 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.

4. Simulations

Using COMSOL Multiphysics we carried out many simulations of both the plug flow system (6)–(9) and the laminar flow one (12)–(10) 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.
Figure 4 and Figure 5 display families of steady traveling wave profiles obtained after longtime simulations of our model equations. On the left panel of Figure 4 we show the height profile for the laminar flow model (solid lines) for two sets of parameter values, as well as for the plug flow model (dashed lines) also for two different sets of parameter values. 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 .

5. Experimental Results

5.1. 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 and Lan [35]; density and surface tension data were found online at https://cameochemicals.noaa.gov/chris/OSF.pdf, accessed on 15 February 2021. 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. [19] and Craster and Matar [18]. 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.

5.2. 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. [19] and Craster and Matar [18], 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. [19] and Craster and Matar [18]: in Figure 7, (b)–(d) correspond to the Convective Instability Regime, and (e) and (f) correspond to the Rayleigh–Plateau and Isolated Droplets Regimes. Table 1 lists the conditions of the experiments that were carried out at temperatures of about 13–15 C.
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.

6. Numerical Data Analysis

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 = 2 g R 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 . Furthermore, 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. 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 Figure 4 and Figure 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
0 L ( h 0 + 0.1 sin ( 2 π x / L ) ) 2 d x = 0 L h ( x , t ) 2 d x ,
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. For calibrating the parameters based on the film profile, the model performance was measured by the mean squared error shown in the figure. 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 (12) results in the ordinary differential equation:
V ( h 2 ) + a 1 ( u ( h 2 1 ) ) = 0
whose solution after some algebraic manipulation provides
u ( z ) = V a 1 + C h ( z ) 2 1 .
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 .

7. 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. [19] and Craster and Matar [18] this mismatch of velocities was close to 40%. In Ji et al. [24] an artificial stabilizing term was introduced to improve the match between the 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
τ = 1 8 f D ρ U 2
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
1 f D = 2 log ( 2.51 Re f D ) .
Note that other similar correlations could also be used. The goal is just to get a rough estimate of the parameters in this regime. We substitute the expression for f D in terms of τ and U, and Re = ρ U T μ and obtain a transcendental equation for U, whose solution yields the mean velocity and corresponding Reynolds number as
U = 11.31 m / s , R e = 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.
Using our approach, it may even be possible to treat transitional flows (in between laminar and fully turbulent). In the axial momentum Equation (13), the term U / F ( H ) multiplying g on the right-hand side represents the drag force by the fiber on the fluid. We have seen that the function F ( H ) has different forms for low- and high-Reynolds number flows. In particular F ( H ) varies much more rapidly with film thickness for plug flow versus laminar flow. More generally, the dependence of the drag force on U might also be nonlinear. The current model represents a linearization of that force in U. However, for a fully turbulent flow, one might expect that drag force to become quadratic in U. In fluid mechanics, often the drag coefficient or friction factor become relatively constant at high Reynolds number, in which limit the drag is quadratic in U. We can thus replace the term U / F ( H ) in Equation (13) with a more general functional form F ( U , H ) that allows different exponents on U and more general dependence on H to better model transitional or fully turbulent flows.

8. Summary

We presented 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 were 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 showed 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 were also carried out on such a system using safflower oil as the working liquid and a taut fishing line as the fiber. Numerical data analysis was used to find the best set of parameters in the laminar flow model to match the experimental results to the simulations. Good agreement was found between the two, with parameter values that are quite close to their estimates based on the approximate values of the physical parameters.

Author Contributions

Conceptualization, A.N. and M.C.; methodology, Y.R., A.N., L.D. and M.C.; software, Y.R., A.N. and M.C.; validation, Y.R., A.N. and M.C.; formal analysis, Y.R., A.N. and M.C.; investigation, Y.R., A.N., L.D. and M.C.; resources, M.C. and A.N.; data curation, M.C.; writing—original draft preparation, Y.R.; writing—review and editing, A.N. and M.C.; visualization, Y.R., L.D., A.N. and M.C.; supervision, A.N. and M.C.; project administration, A.N. and M.C.; funding acquisition, N/A. All authors have read and agreed to the published version of the manuscript.

Funding

This research received no external funding.

Informed Consent Statement

Not applicable.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Quéré, D.; di Meglio, J.M.; Brochard-Wyart, F. Making van der Waals films on fibers. EPL (Europhys. Lett.) 1989, 10, 335. [Google Scholar] [CrossRef]
  2. Quéré, D. Thin films flowing on vertical fibers. EPL (Europhys. Lett.) 1990, 13, 721. [Google Scholar] [CrossRef]
  3. Quéré, D. Fluid coating on a fiber. Annu. Rev. Fluid Mech. 1999, 31, 347–384. [Google Scholar] [CrossRef]
  4. Giorgiutti-Dauphiné, F.; Duprat, C.; Ruyer-Quil, C.; Hulin, J.P.; Trevelyan, P.M.; Kalliadasis, S. Experimental and numerical study of film flows down fibers at moderate Reynolds numbers. In Proceedings of the APS Division of Fluid Dynamics Meeting Abstracts, Tampa Bay, FL, USA, 19–21 November 2006; Volume 59, p. HD–002. [Google Scholar]
  5. Halpern, D.; Wei, H.H. Slip-enhanced drop formation in a liquid falling down a vertical fibre. J. Fluid Mech. 2017, 820, 42–60. [Google Scholar] [CrossRef]
  6. Nozaki, T.; Kaji, N.; Mori, Y.H. Heat transfer to a liquid flowing down vertical wires hanging in a hot gas stream: An experimental study of a new means of thermal energy recovery. In Proceedings of the International Heat Transfer Conference Digital Library, Kyongju, Korea, 23–28 August 1998; Begell House Inc.: Danbury, CT, USA, 1998. [Google Scholar]
  7. Zeng, Z.; Sadeghpour, A.; Warrier, G.; Ju, Y.S. Experimental study of heat transfer between thin liquid films flowing down a vertical string in the Rayleigh-Plateau instability regime and a counterflowing gas stream. Int. J. Heat Mass Transf. 2017, 108, 830–840. [Google Scholar] [CrossRef] [Green Version]
  8. Sadeghpour, A.; Zeng, Z.; Ji, H.; Ebrahimi, N.D.; Bertozzi, A.; Ju, Y. Water vapor capturing using an array of traveling liquid beads for desalination and water treatment. Sci. Adv. 2019, 5, eaav7662. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  9. Sadeghpour, A.; Zeng, Z.; Ju, Y.S. Effects of nozzle geometry on the fluid dynamics of thin liquid films flowing down vertical strings in the Rayleigh–Plateau regime. Langmuir 2017, 33, 6292–6299. [Google Scholar] [CrossRef] [Green Version]
  10. Kalliadasis, S.; Chang, H.C. Drop formation during coating of vertical fibres. J. Fluid Mech. 1994, 261, 135–168. [Google Scholar] [CrossRef]
  11. Ruyer-Quil, C.; Kalliadasis, S. Wavy regimes of film flow down a fiber. Phys. Rev. E 2012, 85, 046302. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  12. Yu, L.; Hinch, J. The velocity of large viscous drops falling on a coated vertical fibre. J. Fluid Mech. 2013, 737, 232–248. [Google Scholar] [CrossRef]
  13. Ji, H.; Sadeghpour, A.; Ju, Y.S.; Bertozzi, A.L. Modelling film flows down a fibre influenced by nozzle geometry. J. Fluid Mech. 2020, 901, R6. [Google Scholar] [CrossRef]
  14. Frenkel, A. Nonlinear theory of strongly undulating thin films flowing down vertical cylinders. EPL (Europhys. Lett.) 1992, 18, 583. [Google Scholar] [CrossRef]
  15. Chang, H.C.; Demekhin, E.A. Mechanism for drop formation on a coated vertical fibre. J. Fluid Mech. 1999, 380, 233–255. [Google Scholar] [CrossRef]
  16. Marzuola, J.L.; Swygert, S.R.; Taranets, R. Nonnegative weak solutions of thin-film equations related to viscous flows in cylindrical geometries. J. Evol. Equ. 2019, 2019, 1–23. [Google Scholar] [CrossRef] [Green Version]
  17. Ji, H.; Taranets, R.; Chugunova, M. On travelling wave solutions of a model of a liquid film flowing down a fibre. arXiv 2020, arXiv:2006.04994. [Google Scholar]
  18. Craster, R.; Matar, O. On viscous beads flowing down a vertical fibre. J. Fluid Mech. 2006, 553, 85–105. [Google Scholar] [CrossRef]
  19. Kliakhandler, I.; Davis, S.H.; Bankoff, S. Viscous beads on vertical fibre. J. Fluid Mech. 2001, 429, 381–390. [Google Scholar] [CrossRef]
  20. Chinju, H.; Uchiyama, K.; Mori, Y.H. String-of-beads flow of liquids on vertical wires for gas absorption. AIChE J. 2000, 46, 937–945. [Google Scholar] [CrossRef]
  21. Duprat, C.; Ruyer-Quil, C.; Kalliadasis, S.; Giorgiutti-Dauphiné, F. Absolute and convective instabilities of a viscous film flowing down a vertical fiber. Phys. Rev. Lett. 2007, 98, 244502. [Google Scholar] [CrossRef] [Green Version]
  22. Duprat, C.; Ruyer-Quil, C.; Giorgiutti-Dauphiné, F. Spatial evolution of a film flowing down a fiber. Phys. Fluids 2009, 21, 042109. [Google Scholar] [CrossRef] [Green Version]
  23. Ding, Z.; Liu, Z.; Liu, R.; Yang, C. Breakup of ultra-thin liquid films on vertical fiber enhanced by Marangoni effect. Chem. Eng. Sci. 2019, 199, 342–348. [Google Scholar] [CrossRef]
  24. Ji, H.; Falcon, C.; Sadeghpour, A.; Zeng, Z.; Ju, Y.; Bertozzi, A. Dynamics of thin liquid films on vertical cylindrical fibres. J. Fluid Mech. 2019, 865, 303–327. [Google Scholar] [CrossRef]
  25. Reisfeld, B.; Bankoff, S. Non-isothermal flow of a liquid film on a horizontal cylinder. J. Fluid Mech 1992, 236, 167196. [Google Scholar] [CrossRef]
  26. Trifonov, Y.Y. Steady-state traveling waves on the surface of a viscous liquid film falling down on vertical wires and tubes. AIChE J. 1992, 38, 821–834. [Google Scholar] [CrossRef]
  27. Ruyer-Quil, C.; Treveleyan, P.; Giorgiutti-Dauphiné, F.; Duprat, C.; Kalliadasis, S. Modelling film flows down a fibre. J. Fluid Mech. 2008, 603, 431–462. [Google Scholar] [CrossRef]
  28. Ruyer-Quil, C.; Trevelyan, S.; Giorgiutti-Dauphiné, F.; Duprat, C.; Kalliadasis, S. Film flows down a fiber: Modeling and influence of streamwise viscous diffusion. Eur. Phys. J. Spec. Top. 2009, 166, 89–92. [Google Scholar] [CrossRef]
  29. Novbari, E.; Oron, A. Energy integral method model for the nonlinear dynamics of an axisymmetric thin liquid film falling on a vertical cylinder. Phys. Fluids 2009, 21, 062107. [Google Scholar] [CrossRef]
  30. Liu, R.; Ding, Z. Coating flows down a vertical fibre: Towards the full Navier–Stokes problem. J. Fluid Mech. 2021, 914, A30. [Google Scholar] [CrossRef]
  31. Alekseenko, S.; Nakoryakov, V.; Pokusaev, B. Wave formation on vertical falling liquid films. Int. J. Multiph. Flow 1985, 11, 607–627. [Google Scholar] [CrossRef]
  32. Bertozzi, A.L.; Münch, A.; Shearer, M. Undercompressive shocks in thin film flows. Phys. D Nonlinear Phenom. 1999, 134, 431–464. [Google Scholar] [CrossRef]
  33. Bertozzi, A.L.; Shearer, M. Existence of undercompressive traveling waves in thin film equations. SIAM J. Math. Anal. 2000, 32, 194–213. [Google Scholar] [CrossRef] [Green Version]
  34. Bertozzi, A.L.; Münch, A.; Shearer, M.; Zumbrun, K. Stability of compressive and undercompressive thin film travelling waves. Eur. J. Appl. Math. 2001, 12, 253–291. [Google Scholar] [CrossRef] [Green Version]
  35. Diamante, L.M.; Lan, T. Absolute viscosities of vegetable oils at different temperatures and shear rate range of 64.5 to 4835 s−1. J. Food Process. 2014, 2014, 1–6. [Google Scholar] [CrossRef] [Green Version]
Figure 1. Schematic plot of a liquid film on a fiber. The axial coordinate along the fiber axis is x, and the radial distance from that axis is r.
Figure 1. Schematic plot of a liquid film on a fiber. The axial coordinate along the fiber axis is x, and the radial distance from that axis is r.
Fluids 06 00281 g001
Figure 2. Plots of ( α 1 ) and ( α 2 ) versus wavenumber k for the laminar-flow model with h 0 = 2.5 , u 0 = f 1 ( h 0 ) , a 1 = 1 , and b = 11 . While ( α 2 ) stays negative for all k, α 1 has a positive real part over a finite range of wavenumbers k near k = 0 , exhibiting a maximum growth rate at a wavenumber close to 0.3.
Figure 2. Plots of ( α 1 ) and ( α 2 ) versus wavenumber k for the laminar-flow model with h 0 = 2.5 , u 0 = f 1 ( h 0 ) , a 1 = 1 , and b = 11 . While ( α 2 ) stays negative for all k, α 1 has a positive real part over a finite range of wavenumbers k near k = 0 , exhibiting a maximum growth rate at a wavenumber close to 0.3.
Fluids 06 00281 g002
Figure 3. Plots of ( α 1 ) and ( α 2 ) versus wavenumber k for the plug-flow model with h 0 = 2 , u 0 = f ( h 0 ) , a = 100 , b = 2 and c = 0.05 . While ( α 2 ) stays negative for all k, α 1 has a positive real part over a finite range of wavenumbers k near k = 0 .
Figure 3. Plots of ( α 1 ) and ( α 2 ) versus wavenumber k for the plug-flow model with h 0 = 2 , u 0 = f ( h 0 ) , a = 100 , b = 2 and c = 0.05 . While ( α 2 ) stays negative for all k, α 1 has a positive real part over a finite range of wavenumbers k near k = 0 .
Fluids 06 00281 g003
Figure 4. The left panel illustrates the difference in the maximal height of the traveling waves for the laminar model (solid lines: blue corresponds to a 1 = 0.1 , b = 11 and yellow corresponds to a 1 = 1.5 , b = 13 ) and for the plug flow model (dashed lines: blue corresponds to a = 0.2 , b = 10 , c = 1 and yellow corresponds to a = 0.4 , b = 12 , c = 3 ) using h 0 = 2.29 in both models. The right panel shows the variation of the shape of the traveling wave for different h 0 values as it ranges from 1.2 to 3.0 for fixed values of a 1 = 1 and b = 11 in the laminar flow model. The non-dimensional velocity of the traveling wave increases from 0.17 to 2.23 as h 0 increases.
Figure 4. The left panel illustrates the difference in the maximal height of the traveling waves for the laminar model (solid lines: blue corresponds to a 1 = 0.1 , b = 11 and yellow corresponds to a 1 = 1.5 , b = 13 ) and for the plug flow model (dashed lines: blue corresponds to a = 0.2 , b = 10 , c = 1 and yellow corresponds to a = 0.4 , b = 12 , c = 3 ) using h 0 = 2.29 in both models. The right panel shows the variation of the shape of the traveling wave for different h 0 values as it ranges from 1.2 to 3.0 for fixed values of a 1 = 1 and b = 11 in the laminar flow model. The non-dimensional velocity of the traveling wave increases from 0.17 to 2.23 as h 0 increases.
Fluids 06 00281 g004
Figure 5. The left panel shows the dependence of the shape of the traveling wave on the parameter a 1 as we increase it from 0.1 (blue) to 1.25 (yellow) for fixed values h 0 = 2.5 and b = 11 in the laminar model with the non-dimensional velocity of the traveling wave decreasing from 1.338 to 1.331 . The right panel shows the dependence of the shape of the traveling wave on the parameter b as we change it from 10 to 13 for fixed values a 1 = 1 and h 0 = 2.5 for the laminar flow model; the non-dimensional velocity of the traveling wave decreases from 1.34 to 1.32 in this case.
Figure 5. The left panel shows the dependence of the shape of the traveling wave on the parameter a 1 as we increase it from 0.1 (blue) to 1.25 (yellow) for fixed values h 0 = 2.5 and b = 11 in the laminar model with the non-dimensional velocity of the traveling wave decreasing from 1.338 to 1.331 . The right panel shows the dependence of the shape of the traveling wave on the parameter b as we change it from 10 to 13 for fixed values a 1 = 1 and h 0 = 2.5 for the laminar flow model; the non-dimensional velocity of the traveling wave decreases from 1.34 to 1.32 in this case.
Fluids 06 00281 g005
Figure 6. Experimental setup scheme.
Figure 6. Experimental setup scheme.
Fluids 06 00281 g006
Figure 7. (a) Nearly uniform coating flow over a 2 mm diameter fishing line with a 2.4 mm hole; (b,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).
Figure 7. (a) Nearly uniform coating flow over a 2 mm diameter fishing line with a 2.4 mm hole; (b,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).
Fluids 06 00281 g007
Figure 8. (a) Zoomed photo of the droplets with edge detection indicated below by the series of red dots. (b) Match between the edge of the drop and a numerical simulation corresponding to the parameters found by Machine Learning: h 0 = 2.29 , b = 12.32 and a = 0.95 .
Figure 8. (a) Zoomed photo of the droplets with edge detection indicated below by the series of red dots. (b) Match between the edge of the drop and a numerical simulation corresponding to the parameters found by Machine Learning: h 0 = 2.29 , b = 12.32 and a = 0.95 .
Fluids 06 00281 g008
Figure 9. (a) Plot of actual (blue) vs. predicted (yellow) values of a 1 found by the GBR. (b) Accuracy of the parameter a 1 value on a test set vs. the number of first Fourier coefficients used (blue: h, black: 1 / ( h 2 1 ) , yellow: both). (c) Learning curve: the accuracy of the parameter a 1 vs. the size of the training dataset (yellow: 961 grid values, blue: first 6 Fourier coefficients). (d) Plot of actual (blue) vs. predicted (yellow) values of b found by the GBR. (e) Accuracy of parameter b on a test set vs. the number of first Fourier coefficients used (blue: h, black: 1 / ( h 2 1 ) , yellow: both). (f) Learning curve: the accuracy of the parameter b vs. the size of the training dataset (yellow: 961 grid values, blue: first 6 Fourier coefficients).
Figure 9. (a) Plot of actual (blue) vs. predicted (yellow) values of a 1 found by the GBR. (b) Accuracy of the parameter a 1 value on a test set vs. the number of first Fourier coefficients used (blue: h, black: 1 / ( h 2 1 ) , yellow: both). (c) Learning curve: the accuracy of the parameter a 1 vs. the size of the training dataset (yellow: 961 grid values, blue: first 6 Fourier coefficients). (d) Plot of actual (blue) vs. predicted (yellow) values of b found by the GBR. (e) Accuracy of parameter b on a test set vs. the number of first Fourier coefficients used (blue: h, black: 1 / ( h 2 1 ) , yellow: both). (f) Learning curve: the accuracy of the parameter b vs. the size of the training dataset (yellow: 961 grid values, blue: first 6 Fourier coefficients).
Fluids 06 00281 g009
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.
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.
Fluids 06 00281 g010
Table 1. Experimental conditions.
Table 1. Experimental conditions.
Hole Size (mm)Flow Rate (g/s)Droplet Size (mm)Gap Size (mm)Droplet Velocity (mm/s)
1.60.00673-5various10
20.03664315
2.40.08335730
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Ruan, Y.; Nadim, A.; Duvvoori, L.; Chugunova, M. Liquid Films Falling Down a Vertical Fiber: Modeling, Simulations and Experiments. Fluids 2021, 6, 281. https://doi.org/10.3390/fluids6080281

AMA Style

Ruan Y, Nadim A, Duvvoori L, Chugunova M. Liquid Films Falling Down a Vertical Fiber: Modeling, Simulations and Experiments. Fluids. 2021; 6(8):281. https://doi.org/10.3390/fluids6080281

Chicago/Turabian Style

Ruan, Yadong, Ali Nadim, Lekha Duvvoori, and Marina Chugunova. 2021. "Liquid Films Falling Down a Vertical Fiber: Modeling, Simulations and Experiments" Fluids 6, no. 8: 281. https://doi.org/10.3390/fluids6080281

APA Style

Ruan, Y., Nadim, A., Duvvoori, L., & Chugunova, M. (2021). Liquid Films Falling Down a Vertical Fiber: Modeling, Simulations and Experiments. Fluids, 6(8), 281. https://doi.org/10.3390/fluids6080281

Article Metrics

Back to TopTop