Numerical Bifurcation Analysis of a Film Flowing over a Patterned Surface through Enhanced Lubrication Theory

: The bifurcation analysis of a ﬁlm falling down an hybrid surface is conducted via the numerical solution of the governing lubrication equation. Instability phenomena, that lead to ﬁlm breakage and growth of ﬁngers, are induced by multiple contamination spots. Contact angles up to 75 ◦ are investigated due to the full implementation of the free surface curvature, which replaces the small slope approximation, accurate for ﬁlm slope lower than 30 ◦ . The dynamic contact angle is ﬁrst veriﬁed with the Hoffman–Voinov–Tanner law in case of a stable ﬁlm down an inclined plate with uniform surface wettability. Then, contamination spots, characterized by an increased value of the static contact angle, are considered in order to induce ﬁlm instability and several parametric computations are run, with different ﬁlm patterns observed. The effects of the ﬂow characteristics and of the hybrid pattern geometry are investigated and the corresponding bifurcation diagram with the number of observed rivulets is built. The long term evolution of induced ﬁlm instabilities shows a complex behavior: different ﬂow regimes can be observed at the same ﬂow characteristics under slightly different hybrid conﬁgurations. This suggest the possibility of controlling the rivulet/ﬁlm transition via a proper design of the surfaces, thus opening the way for relevant practical application.


Introduction
The evolution of a thin liquid film is involved in several engineering applications. For example, in in-flight icing phenomenon, the ice accretion on the airplane wings during the take-off and the landing stages is driven by the evolution of a liquid layer fed by supercooled droplets from the clouds, and may severely affect flight safety [1]. The prediction of the liquid behavior, which may evolve as a droplets population, an ensemble of rivulets or a continuous film is crucial to estimate both the induced ice surface roughness and the extent of the runback water flow. Liquid film coating is driven by the evolution of a thin liquid layer [2], which is required to cover the solid surface as a continuous film in order to form a uniform layer, while being kept as thin as possible to ensure a proper coating efficiency. Even in chemical engineering, liquid film evolution is involved in absorption and distillation processes. In CO 2 absorption through structured packing, a liquid solvent, that falls down a collection of corrugated sheets, captures the exhaust CO 2 , which flows up through the same packed layers via chemical reaction. Since the absorption process is enhanced at maximum interfacial area between liquid solvent and gas solute, the continuous film regime is required. However, the liquid layer must be as thin as possible in order to avoid flooding condition occurrence. Empirical models, which correlates the interfacial area to the liquid hold-up inside structured packing, are available in literature [3,4]. However, such models assume that a continuous film flows through the packed layer and, thus, the effective liquid behavior, which may also arrange as a collection of rivulets, seriously affecting the efficiency, is not considered. Furthermore, thin liquid layers are also involved in: fluid dynamics inside a lubricated bearing [5][6][7], where the fluid is confined between two solid moving walls and different challenges arise, such as complex regimes map function of velocity and surface topology and roughness, which do not appear in free surface film problems; non-Newtonian fluid motion inside micro-systems [8].
On the other hand, film instability phenomena are of great interest from a mathematical point of view and have been largely investigated in literature in a number of numerical [9][10][11][12][13], analytical [14,15], and experimental [16,17] studies. The rupture of 1D film driven by capillary forces and viscous dissipation over a heterogeneous surface was numerically studied in [10], assuming lubrication approximation and modeling the surface wettability via disjoining pressure, and the possible configurations were mapped as a function of amplitude (multiplying disjoining pressure) and periodicity (pattern length) of the imposed pattern defining the heterogeneous surface. The occurrence of 2D finger instability over a heterogeneous surface was numerically investigated by Zhao and Marshall [12] assuming lubrication approximation. Dry patch generation and stability were experimentally studied in [16,17], introducing a local perturbance in order to induce the rupture of a continuous film pattern flowing down an inclined plate. Here, a film model based on enhanced lubrication theory, capable of simulate relatively high contact angles due to a full surface curvature formulation, is validated and used to numerically analyze the stability of the front of a thin film flowing over an inclined plate, characterized by an heterogeneous surface (i.e., non-uniform surface wettability). According to literature [12,13,[18][19][20][21], the surface wettability is implemented assuming the existence of a precursor film thickness and using the disjoining pressure model. The implementation of the exact interface curvature to compute the capillary pressure allows to investigate contact angles up to θ s = 75 • [20]. A in-house code, previously developed in FORTRAN language and validated by the Authors [19,20], was updated and parallelized for a shared memory machine using OpenMP library to speed up computations. First, the stability of a 1D film down a vertical plate is investigated for different surface wettabilities and the dynamic contact angle is verified with the well known Tanner-Hoffman-Voinov law. Then, parametric computations are run to conduct the bifurcation analysis of a 2D film flowing down an inclined plate, initially dry, characterized by an heterogeneous surface with non-uniform surface wettability, that induce film instability. Two different solid surface morphology are considered: a single, horizontal array of contamination spots with decreased surface wettability, θ s = 75 • , hit by a falling film, may lead to finger instability with generation of stationary or shedding rivulets, depending on the contamination spot spacing and on the flow characteristics; surface configurations occurring in practical problems are mimicked via the imposition of a random wettability distribution, with θ s up to 60 • , and the effect of the solid surface morphology on the liquid layer evolution is investigated. The first test case is aimed at testing the possibility of controlling the rivulet/film transition via a proper design of the heterogeneous surface, while the second test case is aimed at verifying the effect of solid surface unwanted heterogeneity due to roughness or fouling. Both the onset of finger instabilities and their long term evolution are considered. The former confirms literature linear stability analysis. The latter shows the actual practical implication of such instabilities: controlling the solid surface wettability properties allows to achieve different liquid layer regimes at the same flow characteristics.

Mathematical Model
Consider a thin liquid film flowing over an inclined plate, driven by gravity, surface tension, and viscous dissipation forces. Let α be the plate slope, h 0 and u 0 the undisturbed film thickness and the undisturbed film velocity, that can be expressed via Nusselt film theory [22], Lubrication theory allows to integrate continuity equation across the film thickness. Further assuming low film Reynolds number, liquid inertia can be neglected and the momentum equation gives a parabolic velocity profile, which can be averaged across the film thickness and substituted in the continuity equation, giving the governing lubrication equation [18], ∂h ∂t with h and p being the unknown film thickness and the local film pressure. The pressure field is given by the hydrostatic, capillary, and disjoining contributions [18], with y and 2 κ being the plate downhill direction and the liquid-gas free surface curvature. The disjoining pressure, which models the intermolecular forces between solid and liquid surfaces, was integrated with lubrication theory by Schwartz and Eley [18], with d h 0 and θ s being the precursor film thickness and the static contact angle. Since the disjoining coefficient B was derived by Schwartz and Eley [18] under the assumption of infinitesimal value of the precursor film thickness, the correction coefficient f , with δ = d/h 0 , was introduced by Zhao and Marshall [12] in order to extend disjoining pressure to non-infinitesimal values of δ, as it happens in numerical simulations. The small slope approximation, which is usually adopted to estimate the free surface curvature, is not accurate when the film slope is higher than 30 • [23]. However, it was proved in [20] that contact angles up to θ s = 60 • can be investigated via full implementation of the free surface curvature [24], Defining the film Bond number as the ratio between gravity and surface tension forces, and introducing the following non-dimensional quantities, with L 0 being the characteristic length scale [12], the governing lubrication equations, Equations (3), (4), and (9), are recast as Equations (13)- (15) are numerically solved on a orthogonal, structured grid of n x ×n y elements via Finite Volume Method. A in-house FORTRAN source code, previously developed and validated by the Authors [19,20], was updated. In particular, the first order upwind scheme was replaced by the second order centered scheme suggested by Diez and Kondic [9] for the discretization of film fluxes. The Alternating Direction Implicit (ADI) approximate factorization presented by Witelski and Bowen [25] was implemented for time marching. Thus, the film volumetric flux, is decomposed into two components, regrouping the higher derivatives in F, Linearizing the higher order derivatives terms in F, and applying the approximate factorization [18,25], the sparse, non-linear algebraic system can be splitted into two pentadiagonal, linear systems to be solved at each time step.
According to [18,25], the higher order cross derivatives arising in F are treated explicitly. Since first order accuracy is only ensured from two-step schemes when approximate ADI factorization is considered [25], the implicit Euler scheme was implemented in order to promote convergence. The integration time step is dynamically adjusted limiting the allowed increment in the film thickness at successive integration steps. The source code was parallelized for shared memory machines using OpenMP in order to speed up computations.
In particular, the two pentadiagonal systems from ADI factorization are decomposed into n x and n y independent sub-systems, which are assigned to different threads. The visualization of the numerical results was performed using the open source graphing utility gnuplot.

Numerical Setup
Three different configurations were considered for numerical computations: • First, the stability and dynamics of a 1D falling film down an inclined plate was investigated. A uniform substrate wettability was considered. Symmetry condition, Q ·n = 0 and ∇H ·n = 0, was applied through the lateral boundaries, X = 0 and X = L X , meaning that an infinitely wide plate is simulated. Q ·n = 1 and H = 1 were imposed through the inlet section, located at Y = 0, while fully developed flow, ∇H ·n = 0 and ∇P ·n = 0, was assumed through Y = L Y . L X = 1 was set far all the 1D computations. • Second, a 2D film flowing down an inclined plate, characterized by a single array of equally spaced contamination spots with decreased surface wettability, was simulated, Figure 1a. Defining L X as the distance between adjacent contamination spots, only a portion of width L X /2 was considered, due to the problem symmetry. Thus, symmetry condition was imposed through X = 0 and X = L X /2, while inlet condition and fully developed flow were applied at Y = 0 (top-most section) and Y = L Y (bottomest section) respectively. The low-wettable square patch of dimension L 0 defining the contamination spot was located at Y = 4, far enough from the inlet section (Y = 0) to avoid boundary effects. • Finally, the hybrid surface shown in Figure 1b, characterized by the following periodic, isotropic smooth random wettability distribution from [26], was investigated, with a jk being random coefficients falling inside a given range, ψ j,k ∈ [−π, +π] being the random phase of cosine function, m 0 and n 0 being the number of harmonics over x and y. In order to ensure physical consistency, periodic conditions were applied through X = 0 and X = L X . Inlet condition and fully developed flow were again applied through The solution was always initialized imposing a dry domain, i.e., assuming a uniform film thickness, equal to the characteristic precursor film over the whole computational domain at the initial stage of the computation, H(X, Y)| T 0 = δ.

Dynamic Contact Angle
The stability of a 1D film falling down an inclined plate was investigated, running parametric computations at decreasing film Bond number, until film rupture was observed. Different surface wettabilities were investigated, with the imposed contact angle varying in the wider range allowed by the full free surface curvature formulation, θ s ∈ [30 • , 60 • ]. The plate slope was set to α = 60 • for all the simulations. The dynamic (advancing) contact angle θ, computed as the maximum free surface slope at the apparent contact line, was also traced, in case of a stable film, for each computation, considering the self-similar free surface profile of the advancing film. According to literature [12,13,18,20], the precursor film thickness and the disjoining exponents were set to δ = 5 × 10 −2 and n = 3, m = 2, which represents a good compromise between computational costs and accuracy of the numerical solution. A preliminar grid independency analysis was conducted in the most critical configuration. Thus, setting θ s = 60 • , Bo = 5 × 10 −2 , which represents one of the most critical investigated configurations, and progressively decreasing the spatial discretization step ∆X, it was found that ∆X = 1.5 × 10 −2 ensures grid independency, since further decreasing ∆X to 7.5 × 10 −3 leads to a less than 1% deviation on the computed dynamic contact angle. Thus, a spatial discretization step of ∆X ≤ 1.5 × 10 −2 was imposed for all the computations.
In Figure 2a, both the stable film solution and the self-similar droplet deriving from film rupture are presented for a given surface wettability, θ s = 60 • , and slightly different Bond numbers. A bifurcation of the solution is, thus, observed as the stable film configuration is replaced, below the critical Bond number, by the droplet regime, induced by film rupture occurrence. Figure 2b shows the stable film solution at Bo = 5.01, but different surface wettabilities, qualitatively demonstrating that higher free surface slope, are obtained by imposing higher static contact angle in the disjoining pressure. The dynamic contact angle, traced as a function of both the surface wettability and the Bond number in case of a stable advancing 1D film, is presented in Figure 3. It is worth pointing out that the critical Bond number, which defines the transition between film instability occurrence and stable film (leftmost markers in Figure 3), depends on the imposed surface wettability, with a higher contact angle leading to higher critical Bond number (thus, the film is likely to be subject to instability phenomena at lower surface wettability). As widely reported in literature [12,27,28], the cube of the dynamic contact angle depends on the film velocity according to Tanner-Hoffman-Voinov law, with u 0 being the undisturbed film velocity and C being a constant parameter. Expressing u 0 via Nusselt film theory, Equation (1), and retaining the usually neglected logarithmic term [28], Equation (23) can be recast as Following [19,20], the computed dynamic contact angle is compared to Equation (24), with C 0,1,2 being the fitting parameters, in order to verify the effectiveness of disjoining pressure model in the wider range of contact angles allowed by the full curvature formulation. Almost a perfect agreement for θ s up to 60 • can be observed in Figure 3.

Single Array of Contamination Spots
A single array of lower wettability contamination spots, characterized by an increased value of the imposed static contact angle, was investigated, as pointed out in Section 3.1. A somewhat similar setup was investigated by Zhao and Marshall [12], who imposed a sequence of vertical strips, with different contact angle inducing fingering instability. Here, the eventual film instability is induced by localized spots. Furthermore, higher contact angles are here investigated, as allowed by the full modelization of the free surface curvature. The influence of both the dimensionless distance L X between successive contamination spots and the film Bond number, Equation (10), on the resulting liquid layer distribution were analysed, in terms of number of observed rivulets. Parametric computations were run, in the range L X ∈ [5,30] and Bo ∈ [0.11, 0.16]. The static contact angle was set to θ s = 60 • over the whole computational domain, except inside the contamination spot, where the contact angle was set to an increased value of 75 • , while the plate slope was set to α = 60 • . According to literature [12,13,18,20], the precursor film thickness was set to δ = 5 × 10 −2 , while the disjoining exponents were set to n = 3, m = 2. Spatial discretization steps of ∆X, ∆Y ≤ 1.5 × 10 −2 were imposed in order to ensure grid independency, while L Y = 20 was chosen as the plate length along the downhill direction, in order to let the liquid film develop and, thus, avoid boundary effects.
All the observed configurations, defined by the number of rivulets flowing within a periodic length L X , are summarized in the bifurcation diagram, Figure 4, as a function of the distance between neighbor contamination spots and of the Bond number. Figure 4a shows the number of rivulets at T = 10 (thus, right after the undisturbed film gets over the contamination spots). Since the liquid bulk has not reached the outlet section at T = 10, the number of rivulets was computed as the number of streamwise peaks of the evolving contact line, which is defined by the position of the local maximum of the free surface slope. It is important to point out that the number of observed rivulets mainly depends on the dimensionless distance between neighbor contamination spots, which in turn depends on the imposed film Bond number, In fact, critical lengths defining the transition between different configurations can be identified: L 01 ∼6 corresponds to the appearance of a single rivulet; L 12 ∼16 defines the transition from 1 to 2 rivulets; L 23 ∼25 defines the transition between 2 to 3 rivulets. Additionally, defining ∆L 1 = L 12 − L 01 and ∆L 2 = L 23 − L 12 , it comes that new rivulets forms at ∆L i close to the characteristic perturbation wavelength, λ = 9.4, deriving from linear stability analysis of the governing lubrication equation [12].  The above defined instabilities may evolve into stable dry patches or may be shed away, leaving a continuous film upstream of an advancing unstable front. These regimes are mapped in Figure 4b, which depicts the steady state solution. The figure offers a complex behavior, since it includes long term evolution and effective growth of the local perturbance, and the effect of the imposed restrictions (such as periodicity) which are not accounted for in the linear stability analysis. Figure 4b shows that five different configurations can be observed: fully wetted domain (N = 0), usually observed at low distance between neighbor low-wettable spots and at high Bond number; a single rivulet (N = 1), which flows between neighbor contamination spots; 2 rivulets (N = 2), one flowing between neighbor contamination spots and one wetting the contamination spot, as shown by     For a better understanding of the contact line dynamics, different screenshots of the moving contact line are taken from two different computations at Bo = 0.11 and slightly different L X . In Figure 8, which refers to the same test case as Figure 7, L X = 30, the contact line position is traced at successive instants: 3 rivulets have already formed at T = 10, Figure 8a; after the rivulets develop and become narrow, there is enough space for one more rivulet front to form below to the contamination spot and drain, Figure 8b,c; as long as the generated dry patches are stable, all the rivulets just settle down and remain as stationary rivulets, Figure 8d.
The moving contact line was traced, at the same instants, during a computation at slightly different distance between contamination spots, L X = 27: similarly to the previous case (L X = 30), 3 rivulets forms after the liquid front meets the contamination spot, Figure 9a; again, after the rivulets settle down, a new rivulet front forms below the contamination spot, Figure 9b; however, the dry patches, induced by the formation of the new, central rivulet, slowly moves downstream, generating a bottle-neck shape that evolves in a single rivulet, Figure 9c; thus, only 2 stationary rivulets are observed, Figure 9d. Furthermore, one of the stable rivulets is larger then the other one, since it collects the liquid flux of the three original rivulets that moved downstream, out of the computational domain. A similar behavior was observed, at Bo = 0.11, when L X = 21, with dry patch closure phenomenon that leads to only 2 stable rivulets, while 3 rivulets can be observed, when L X = 22, in Figure 6b. When fully wetted condition is detected, as it happens for Bo = 0.11 and L X = 20, the liquid front simply overtakes the contamination spot. Fully wetted condition is progressively incentivated at higher Bond number, see Figure 4b, as the gravitational force dominates the surface tension, ensuring stability of the liquid film. However, it is quite interesting for practical applications, which often requires the existence of stable and thin films at dominating surface tension forces, that the fully wetted condition can be obtained even at the lower Bond numbers, under restricted geometrical characteristics of the solid surface.
In order to test the consistency of the applied boundary conditions (i.e., half of the periodic length investigated, contamination spot located at X = 0 and symmetry conditions applied to X = 0 and X = L X ), a larger domain of width 2 L X (thus, including 2 contamination spots, located at X = 14.3, 34.3) with periodic conditions, applied through X = 0 and X = 2 L X , was also simulated. In fact, the latter test case allows the film to evolve in a larger domain (4 times the characteristic perturbation length λ cr from linear theory), mitigating the artificial constraints deriving from forcing the film to follow the geometrical symmetry. A configuration characterized by low Bond number, Bo = 0.10, giving a film subject to instability phenomena even when weak perturbations are introduced, was considered. As demonstrated by Figure 10, which shows the liquid layer distribution resulting from the two different computations at the same instant T = 125, the same number of rivulets per unit length is predicted, meaning that the results proposed in the bifurcation diagram, Figure 4b, are statistically consistent, although the solution is less regular and may also have some oscillations in time.

Randomly Generated Heterogeneous Surface
A general heterogeneous surface, characterized by a random, periodic distribution of the static contact angle, implemented through Equation (21), was also investigated. Such a test case is aimed to mimic the typical surfaces occurring in practical application. A large computational domain, characterized by L X = 40 and L Y = 50, was considered in order to let the induced perturbance grow without any numerical constraint. The plate slope and the Bond number were set to α = 60 • and Bo = 0.1, while the static contact angle was ranged in θ s ∈ [45 • , 60 • ] over the heterogeneous surface. The characteristics of the heterogeneous surface are imposed through the number of harmonics (m 0 , n 0 ) considered in Equation (21), which defines the wavelength parameters, Λ X = L X /m 0 , Λ Y = L Y /n 0 : in order to ensure isotropy, Λ = Λ X = Λ Y was always imposed. The precursor film thickness and the disjoining exponents were again set to δ = 5 × 10 −2 and n = 3, m = 2. A spatial discretization step of ∆X, ∆Y ≤ 2.5 × 10 −2 was imposed in order to ensure grid independency. Parametric computations were run at different values of the characteristic length Λ, defining the random surface heterogeneity. The number of rivulets, generated due to finger instability induced by the random contact angle distribution, was then traced at T = 25, in order to statistically investigate the effect of the heterogeneous surface characteristics on the liquid film evolution.
The hybrid surface wettability distributions at different magnitudes of Λ are shown in Figure 11a,b, referred to Λ = 2 × 10 −2 and Λ = 2, respectively. The corresponding film thickness distributions at T = 25 are reported in Figure 12a,b. Fully wetted condition was observed at the end of both the computations. However, different behaviors were observed, with 4 rivulets obtained when Λ = 2 × 10 −2 and 3 rivulets obtained when Λ = 2. Table 1, which reports the number of observed rivulets at T = 25, reveals that a transition from 4 rivulets configuration to 3 rivulets occurs at Λ ≥ 5 × 10 −2 . Fully wetted domain (N riv = 0 at T → ∞ in Table 1) was obtained for all the investigated values of Λ. Thus, the dry patches induced by finger instability are shed away, leaving a continuous film upstream from the unstable advancing liquid front.   80  100  3  0  1  40  50  3  0  2  20  25  3  0  5  8  10  3  0 Thus, the linear theory, which gives a characteristic perturbation wavelength of λ cr = 9.4 (i.e., 4 rivulets in the investigated test cases), as expected does not provide the effective evolution of film instability induced by local perturbance (in turns given by surface heterogeneity) when the magnitude of the surface parameter Λ is comparable to λ cr , i.e., when the actual geometrical details of the perturbation affect contact line movement.

Conclusions
The enhanced lubrication equations were numerically solved. Contact angles up to 60 • (locally increased at 75 • inside localized, lower-wettable contamination spots) were investigated thanks to the full curvature formulation, which substitutes the usually adopted small slope approximation, leading to accurate results for θ s < 30 • [23]. Film rupture occurrence was analyzed for a 1D falling film. The computed, dynamic contact angle was successfully verified with the Hoffman-Voinov-Tanner law for θ s up to 60 • . Then, parametric computations were run for a 2D film over a heterogeneous surface. It was verified that the onset of finger instability, induced by an horizontal array of localized, lower wettable contamination spots with θ s = 75 • , follows the literature linear stability analysis. The long term evolution of the induced instability, affected by the imposed restrictions, was traced under different flow conditions and contamination spot spacing. A complex behavior, including multiple configurations (continuous film and 1 to 4 rivulet regimes), was observed due to the bifurcation occurrence. However, artificially controlling the solid surface morphology (i.e., the contamination spot spacing) allows to achieve the continuous film regime at given flow characteristics (i.e., required film thickness). This is crucial in practical applications, since the whole coverage of a solid surface can be ideally obtained at the lower film thickness (usually more subject to instability phenomena): in absorption/distillation through structured packing, a proper design of the packing geometry may lead to high liquid-gas interface area even at low film thickness, enhancing mass transfer and avoiding flooding condition inside structured packing. The flow down an inclined plate with random surface distribution, which mimics practical configurations, was also investigated: the liquid film behavior can no longer be predicted via linear stability analysis when the magnitude of the characteristic surface length is comparable to the perturbation length from linear theory. Finally, it is worth pointing out that the parallelization of the solver code allowed to speed up computations, while the implemented ADI factorization allowed to reach time integration step up to 10 8 times the step size from Von Neumann stability analysis on the explicit numerical scheme.
Author Contributions: All authors contributed equally to this work. All authors have read and agreed to the published version of the manuscript.