Fingering Instability of a Gravity-Driven Thin Film Flowing Down a Vertical Tube with Wall Slippage

: Inspired by the antiwetting property of pitcher plants, specialists have designed different functional material with slippery surfaces, and a directional slippery surface has been fabricated. This paper considers a gravity-driven liquid ﬁlm coating the interior surface of a vertical tube, and different slippery lengths in the azimuthal direction and the axial direction are taken into account. The evolution equation of coating ﬂow is derived using the thin ﬁlm model, and time responses for two dimensional ﬂow are calculated. Linear stability analysis (LSA) is given based on the traveling wave solutions, demonstrating that the axial slippery effect suppresses the ﬂow instability and causes a larger traveling wave speed. Simultaneously, the azimuthal slippery effect plays a destabilizing role for perturbations with small wavenumbers and it is stabilizing for large wavenumbers. Direct simulations of the ﬁngering ﬂow patterns agree well with the linear stability analysis. Our results offer insight into the inﬂuence of wall slippage on the ﬂow stability of liquid in petroleum engineering.


Introduction
Thin film flowing and droplet spreading are commonly encountered in chemical engineering, nuclear industry and materials sciences. The coating examples show that the free interface of a liquid film easily breaks into fingers resulting from the capillary effect [1][2][3]. The understanding of the dynamics and flow stability of liquid films is essential to control heat and mass transport processes.
To model the thin film flows, lubrication approximation, a thin film model and weighted residual model have been introduced in previous studies [4][5][6]. In the modeling of the liquid film, assumption of a small averaged thickness is crucial [7,8]. These modeling methods were also applied in dynamical contact line problems [9][10][11][12]. Kondic [13] and Lin [14,15] provided a systematic analysis to deal with the contact line instability problem, and the method was easy to implement [16,17]. Based on the governing equations, the relationship between different parameters and the surface instability was analyzed by linear stability analysis (LSA) [18,19].
In recent years, the spreading and falling of the liquid films on different substrates have been received much attention. It was found the curvature and property of the substrate played important roles in the film flows [20][21][22][23][24]. The liquid film flowing down a cylinder/fiber was more unstable than that flowing on an inclined plane [25]. As a common situation in interstifial flow, thin liquid film flowing on porous medium exhibits interesting dynamical phenomena. Ding and Liu [26,27] discovered that the film flowing on a porous vertical cylinder was more unstable and the increasing permeability of the porous medium significantly enhanced the destabilizing effect. The slippery effect often has a connection with wettability and porosity of the materials [28][29][30]. Recently, inspired by the good antiwetting property of nepenthes pitcher plants, specialists designed various slippery liquid-infused porous substrates similar to the microstructures and functions of pitcher plants [31][32][33]. Moreover, liquid-infused surface with directional slippery lengths was fabricated making use of immobilized lubricant menisci [34]. It is necessary to understand the thin film flow on these biomimetic surfaces, as these structures are usually designed for water collection and transportation.
In the analysis, the influence of the porous medium was usually reduced as a slip condition [35,36]. The effect of a slip on the primary instability was investigated by Samanta et al. [37,38]. It was found the slip played a destabilizing role close to the instability onset, while the flow was stabilizing at larger values of the Reynolds number. Ding [39] investigated the coupling action of the thermocapillary effect and slippery effect, concluding that the instability was enhanced as the slippery length increased. Chao [40] studied the dynamics of falling films down a cylinder in the presence of wall slippage, and also found that wall slippage enhanced the capillary instability. All the aforementioned investigations were concentrated on a substrate with slippery length in only one direction. However, the substrates used in biomedical applications usually exhibit anisotropic behavior. Therefore, this paper mainly focuses on the influence of the anisotropic slippery length on the surface instability of thin films flowing down a vertical tube.
The rest of this paper is organized as follows: In Section 2, the mathematical model is conducted based on the thin film model and two dimensional (2D) evolution is displayed. In Section 3, linear stability analysis is carried out based on traveling wave solutions. To verify the linear stability analysis, three dimensional (3D) flow simulations with random perturbations are displayed in Section 4, and some conclusions are briefly summarized in Section 5. Figure 1 shows a gravity-driven liquid film flowing down a vertical tube, where the interior surface of the tube is assumed to be a directional slippery liquid-infused surface. A cylindrical coordinate (r, θ, z) is applied here. The radius of the tube is R and the different slippery lengths in azimuthal and axial directions are denoted by β θ and β z , respectively. Under the driven force of gravity g, the liquid film flows and breaks into fingers. The thickness of the film is defined by h(θ, z, t), where t is time. To study thin film flows with contact lines, a given slip boundary condition or a precursor model are both effective [13]. Considering that the substrate is slippery, we use the precursor model to drive the governing equation, and the thickness of the precursor film is defined by b. The liquid is viscous and incompressible, and the constant density and viscosity of the fluid are ρ and µ, respectively.

Governing Equations and Boundary Conditions
For the liquid, the continuity equation and Navier-Stokes equations are satisfied: where D/Dt represents the material derivative and ∇ 2 is the Laplace operator in cylindrical coordinate, given as follows: where u = (u, v, w) are velocities related to coordinates (r, θ, z) and p is the pressure. In the following text, the pressure of the ambient gas is denoted by p g . For the directional slippery wall, slip boundary conditions are applied in θ direction and z direction while the non-slip boundary condition is applied in the r direction: Here, β 1 and β 2 are the slippery length in θ direction and z direction, respectively, and v r = ∂v ∂r and w r = ∂w ∂r are partial derivative of velocity with respect to the coordinate r. Unless stated otherwise, the subscripts in this paper denote the partial derivative.
At the free surface r = R − h, the balances of normal and tangential stresses should be satisfied: where P = µ[∇u + ∇u T ] and P ∞ = (p − p g )I are the stress tensor and pressure tensor, respectively. ∇ s is the surface gradient and n is the unit vector normal to surface. σ is the surface tension and −∇ s · n means the averaged curvature of the free surface. The kinematic condition is given at the free interface:

Scalings and Non-Dimensionalization
In order to nondimensionalize the system, we introduce the following scales: The characteristic velocity scale W and time scale T are given as: The primes are dropped for simplicity and the nondimensionlized continuity equation and the N-S equations become: 1 r where Re = gWh 2 0 µ denotes the Reynold number. Dimensionless boundary conditions at the fluid-solid interface are expressed as: where β θ = β 1 /R and β z = β 2 /R are the dimensionless slippery lengths. The dimensionless kinematic condition (9) at the liquid-air interface in integral form remains the same: With the scalings, dimensionless boundary conditions for balances of normal and tangential stresses become: where K(h) is the mean curvature, composed by the principal curvature of the interface and the curvature of the substrate. C = µW/σ 0 is the capillary number and τ i (i = 1, 2) denotes the unit vector tangential to surface. Considering that the mean thickness of the liquid film is much smaller than the size in the other two directions, a new coordinate r = R − x is applied in the radial direction. = h 0 /R << 1 is a small parameter, where h 0 is the characteristic thickness of the thin film. The variables can be expanded as a perturbation series of [40]: Using the above solutions and neglecting the terms with higher orders of , the continuity equation and N-S equations at leading order are derived: Applying the rescaled β θ = − −1 β θ and β z = − −1 β z [40], and droping the primes, the boundary conditions at x = 0 become: The boundary conditions at the liquid-air interface become: v 0 Normal stress is balanced by the pressure and surface tension, giving: where the first and the second terms in the right hand are obtained using the Taylor expansion.
The difference between a tube and a cylinder is that the second term in Equation (28) has opposite signs. From Equations (23)-(28), the functions (v, w) of the zeroth order are obtained: Substituting u 0 and v 0 into the kinematic boundary condition of integral form: The evolution equation for the thin film flow with dynamical contact line are finally obtained: where Ca = C/ is the modified capillary number [16,40]. In this study, we mainly focus on the slippery effect, and the effect of capillary number Ca has little effect on the results. Therefore, Ca = 1 is to be set in the following discussions.

Two Dimensional Simulations
The most interesting problem we want to know is whether the traveling contact line is stable. Previous studies has discovered that the steady state shape of 2D flow matters a lot to the surface instability [13,15]. In other words, fingering instability in θ direction has essential relationships with the steady base flow independent of coordinate θ. Besides, the base flow plays an important role in traveling wave solutions and linear stability analysis. Therefore, the two dimensional flow is firstly investigated here, assuming that the evolution equation is independent of the coordinate θ. Then the evolution equation is simplifies to: An initial condition profile h 0 is applied to start the thin film flow. The initial condition is a smooth curve in the form of hyperbolic tangent function [13]: where z 0 is the position of the contact line. From the profile, it is seen that h is constant far away from the contact line. A constant flux boundary condition is approximated well by setting h(0, t) = 1 and h(z max , t) = b, where b is the thickness of the precursor layer. Equation (34) is solved using the finite difference method by a MATLAB code. Firstly, the influence of the slippery length is discussed. The time evolution of the fluid profile with different slippery length is visualized in Figure 2, giving the screenshots for t = 70 and t = 140. In the simulation, the slippery length is set to β z = 0, β z = 0.025, β z = 0.05, β z = 0.075 and β z = 0.1, respectively. Theoretically, a more significant influence will be induced to the responses if we give a larger slippery length β z . However, in fact, slip usually only occurs very close to the contact line, and this condition is a means of alleviating the stress singularity. Hence the value of the slippery length is given in a reasonable range [0, 0.1] in this study. From Figure 2, we note the flow evolves into a traveling wave state after a period, and the traveling wave solutions are to be used in the linear stability analysis (LSA). Resulting from the capillary effect, the steady traveling wave profile has a large ridge at the fluid front, and a higher capillary ridge implies a stronger contact line instability. It has been stated that the profile will tend to the traveling wave shape regardless of the initial condition as long as the boundary conditions were satisfied [16]. Comparing the two graphs, it is seen that the flow travels faster with a bigger value of β z , and the bump becomes lower while the axial slippery effect is larger. Thin films are inclined to break when the bump is high, hence the fluid profile predicts that the thin film flows more stably when the axial slippery length is longer.  Figure 3 shows the situation for different substrate curvature. Tubes of different radius are used to test the influence of the substrate curvature. In the calculation, R = 0.9, R = 1, R = 2, R = 10 and R = ∞ are employed and R = ∞ also corresponds to a vertical plane. It is seen that the capillary ridge appears at the same position after a specified time, indicating that radius of the tubes has no relationship with the traveling wave speed. Also, it can be found that a larger R corresponds to a larger capillary bump, contrary to the phenomenon in the thin film flow coating on a cylinder surface [16,25], because the substrate curvature contributes contrarily in the two cases. For a very large R > 10, the influence of the substrate curvature is insignificant. From the graphs, we can deduced that increasing substrate curvature (smaller R) enhances the stability of the thin film flow, while a higher capillary ridge is induced by the larger substrate curvature, making the flow more unstable. The contribution of the thickness of the precursor layer is investigated and displayed in Figure 4. Two phenomenons deserves attention: the flow travels slower with a thinner precursor layer and a thinner precursor layer induces a larger capillary ridge. The main reason for this is that the thinner precursor film causes more resistance than thicker films; hence, for a thinner precursor film, a lower fluid bump appears and travels more slowly.

Linear Stability Analysis
From the time responses, it is found that traveling wave solutions exist. Therefore, using the moving frame ξ = z − ct, we can define a new thickness function H(ξ), where c is the traveling wave speed. In the following analysis, H denotes the thickness function related to variable ξ while h represents the thickness function related to variables z and t. The governing evolution equation (34) becomes an ordinary differential equation related to only one variable: Integrating Equation (36) once, a third-order differential equation is obtained: where d denotes an integral constant. Equation (37) is solved effectively by a finite difference code, the traveling wave solutions can be obtained. Here, the influence of the slippery length, substrate curvature and thickness of precursor layer are concerned. The traveling wave solutions for different parameters are shown in Figure 5. The influence of the slippery effect on the nonlinear behavior is shown in Figure 5a, and it can be found that a larger slippery length is related to a faster moving interfacial wave following a lower capillary ridge. The lower bump usually corresponds to a more stable flow. For different precursor thickness b, the traveling wave profiles are displayed in Figure 5b, the wave profile for different b presents a similar structure for different β z . Figure 5c gives traveling wave profiles for different R. The capillary ridge becomes larger with a bigger R, indicating that substrate curvature enhances the stability of the thin liquid film. Comparing For the three dimensional problems, the solution can be simply assumed in the form of H(t, ξ, θ) =H(ξ) + H (t, ξ, θ) [17,23], whereH(ξ) is the traveling wave solution (also called base solution) of Equation (37) and H (t, ξ, θ) is the perturbation. Substituting H(t, ξ, θ) into the governing Equation (34) and applying Equation (36), the governing Equation (32) becomes: where Θ z (H) =H 3 3 + β zH 2 and Θ θ (H) =H 3 3 + β θH 2 .The perturbation H (t, ξ, θ) can be expressed as a superposition of a series of Fourier transforms [17,26], expressed as H (t, ξ, y) =Ĥe iqθ+σ(q)t , where q is the wavenumber and σ(q) relates to the growth rate. For every given wavenumber q, substituting H (t, ξ, y) =Ĥe iqθ+σ(q)t into Equation (38), Equation (38) transforms into a standard eigenvalue equation: Here only the largest real σ is to be used representing the largest eigenvalue of L and L is a linear operator with the matrix form: where κ = − 1 R 2 ; I represents the identity matrix, and D i (i = 1, 2, 3, 4) are finite difference matrix of i−th order. Using the difference matrix, we haveH ξ = D 1H ,H ξξ = D 2H ,H ξξξ = D 3H andH ξξξξ = D 4H . Equation (40) is solved numerically to obtain the values of σ(q) for a given wavenumber q, following the work by Kondic [13]. Given a wavenumber, the eigenvalues of the matrix L are obtained by a MATLAB code, and the growth rate is related to the largest real value.
The growth rate σ(q) associated with different cylinder radius is illustrated in Figure 6. Each curve represents the largest eigenvalue of −L versus the wavenumber q, and we can determine the range of unstable wavenumbers as (0, q c ) with σ(q c ) = 0. For small wavenumber q < q c , the flow is always unstable. The dispersion curves illustrate that a smaller cylinder radius results in a larger growth rate and a bigger peak wavenumber (the value corresponding to the maximum growth rate). A smaller cut-off wavenumber q c is also found when the cylinder radius is smaller. For R = ∞, the cylinder can be considered as a vertical plane, demonstrating that the thin film flow down a vertical tube is much more stable when given a certain perturbation. From the plots, it is seen that the substrate curvature greatly promotes the fingering instability, especially for the case of R < 1.
The dispersion curves for different thickness of precursor layer b are visualized in Figure 7. We note that a small thicknesses of the precursor layer is related to a much larger maximum growth rate, showing good consistency with the time responses. Much larger cutoff wavenumbers and peak wavenumbers are found for smaller b, demonstrating that both the strength and the wavenumber range of the surface instability are enhanced by a thinner precursor layer. Figure 8 shows the growth rate for a tube with the isotropic slippery lengths, i.e., β θ = β z . As shown in Figure 8, both the growth rate and the peak wavenumber change small with the increase of slippery length, showing that the slippery effect impedes the surface instability. From the curves, we can find that with a larger slippery length, a smaller unstable range is determined as the q c is smaller. :DYHQXPEHUT *URZWKUDWHσT β z = β θ = 0 β z = β θ = 0.025 β z = β θ = 0.05 β z = β θ = 0.075 β z = β θ = 0.1 Figure 8. Effect of slippery length on the growth rate of the perturbation (β z = β θ ).
As introduced above, the slippery lengths are often different in the azimuthal and axial directions, as the properties of the substrate material are anisotropic, such as fibers and composite materials. Figures 9 and 10 display the dispersion relations for different slippery lengths while β z = β θ is taken into consideration. In Figure 9, the slippery effect in only one direction is considered while in Figure 10, the slippery effect in two directions is considered. From the two figures, the change of the growth rate displays similarly with a difference in the values. It is notable that the two slippery lengths β θ and β z play different roles in the flow stability. From Figures 9a and 10a, a larger axial slippery length corresponds to a smaller growth rate of perturbation. Therefore, the surface instability is suppressed by the axial slippery length. In addition, the change in the growth rate is much more interesting when β θ increases. It can be seen that the growth rate changes slightly in the range [0, q t ], where q t is a turning point. If q < q t , the growth rate increases with the increase of azimuthal slippery length, while q > q t , the growth rate decreases with the increase of azimuthal slippery length. In other words, if the wavenumber of the perturbation is less than q t , the azimuthal slippery effect slightly enhances the surface instability, while the wavenumber of the perturbation exceeds q t , the azimuthal slippery effect suppresses the surface instability. Besides, for the two cases, both the peak wavenumber and the range of σ > 0 decreases with the increase of slippery length.

Numerical Simulation
In this section, the finite element method is applied to solve the governing equation in this section with the help of the open-source software Freefem++ [25,41]. FreeFem++ is a powerful programming language focused on solving partial differential equations using the finite element method. The growth rate σ(q) for imposed perturbations has been obtained by the linear stability analysis, relating to a wavenumber q. Using the wavenumber, we could calculate the most unstable mode of wavelength λ = 2π/q in the azimuthal direction. Here, a similar initial condition used in 2D simulation is chosen, i.e., a smooth surface consisting of two flat regions. The two regions are connected by a tanh function at the initial front position z 0 . Far behind the contact line, the thickness is set to h(θ, 0, t) = 1 and the thickness is set to h(θ, z max , t) = b far front of the contact line. A single mode perturbation is applied at the contact line [17], impelling the contact line to convex slightly into the streamwise direction, given as: where A 0 is the amplitude of the small perturbation. Figure 11 shows the single-mode perturbed thin film flows for different β z . The time profiles for t = 200 are shown in Figure 11a and the time profiles for t = 300 are shown in Figure 11b, and the thickness of the falling films is expressed using the color bar. To compare the influence of the two slippery lengths, a wavenumber of q = 0.2 is chosen, and related wavenumber length is 10π, meaning that the radius of the tube is R = 5. For β z = 0, β z = 0.05 and β z = 0.1, the growth rates σ(q) are 0.01058, 0.00678 and 0.00499, receptively. It has been concluded that the axial slippery length causes a larger traveling wave speed and suppresses the flow instability. From the 3D simulations, the flow travels much faster with a larger β z as expected, and the fingers for β z = 0 are longer than that for β z = 0.05 and β z = 0.1, showing good consistency with the prediction of the linear stability analysis. From the 3D plots, only one finger appears as a single-mode perturbation is applied, and the circumference is a precise length of one single wavelength.
Time evolution of single-mode perturbed fluid films flowing down vertical tubes for different β θ are exhibited in Figures 12 and 13. From the LSA, it is seen that the slippery effect plays opposite roles while different wavenumbers are applied. To verify the roles, the cases for q = 0.2 and q = 0.4 are simulated. In Figure 12, the wavenumber of 0.2 is also applied, and the radius is R = 5. From the two graphs, it is seen the traveling wave speed has no relationship with the circumferential slippery effect. When q = 0.2, the growth rate is 0.00490, 0.00495 and 0.00499, receptively for β θ = 0, β θ = 0.05 and β θ = 0.1. The contour plots in Figure 12 demonstrates that the fingers for β θ = 0 are a little longer than that for a larger β θ . However for q = 0.4, the circumference is two times of the a single wavelength, and two fingers appears in each graph, as shown in Figure 13. For q = 0.4, the growth rates are 0.00651, 0.00469 and 0.00275, corresponding to β θ = 0, β θ = 0.05 and β θ = 0.1. From the measurement, the maximum height is 1.5545, 1.5465 and 1.5333, respectively, also showing that the instability is impeded by the β θ . Figure 14 displays the time evolution of single-mode perturbed fluid films down vertical tubes of different radius R. Here R = 10, R = 15 and R = 20 is applied, and a single-mode perturbation of q = 0.4 is still applied as an initial condition. The contour plot shows the film profile at t = 200 and t = 300. The largest influence of the radius is changing the number of the fingers. It is found that with a larger radius of the cylinder, the number of fingers in the contact line is also larger. In Figure 13, the number of fingers is 2, and here the number is 4, 6 and 8, respectively. The number of the fingers is determined by the wavelength, hence the number of the fingers is the divisor of the circumference and wavelength related to the perturbation.

Conclusions
In the present paper, we investigated the fingering instability of a coating flow down a vertical tube with dynamic contact line, and the influence of the directional slippery substrate on the flow instability was analyzed. Based on the thin film theory, the governing equation for the thin film flow was derived, then the time responses for the axisymmetry model for different parameters were simulated, including the precursor thicknesses, the axial slippery length and the radius of the tube. Linear stability analysis was carried out based on the traveling wave solutions. Finally, three dimensional numerical simulations were given by the finite element method with the help of the open-source software Freefem++. Numerical simulations showed good agreement with linear stability analysis.
From the numerical simulations and the LSA, some conclusions were obtained: The thickness of the precursor layer induced a smaller capillary bump and causes the contact line to travel faster, while the radius of the tube promoted the surface instability. The slippery effect in the axial direction and the azimuthal direction acted differently thin film flows with dynamic contact lines. The axial slippery effect always played a stabilizing role while the azimuthal effect played two roles. When the wavenumber was smaller than a transfer wavenumber q t , the azimuthal effect slightly enhanced the flow instability, and when the wavenumber exceeded the critical value q t , the azimuthal effect impeded the flow instability.