{'fraction_non_alphanumeric': 0.03237924865831843, 'fraction_numerical': 0.012641621943947525, 'mean_word_length': 4.936637168141593, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Bifurcation analysis of liquid films over low wettability surfaces

Thin liquid layer evolution over a solid substrate and film instability phenomena are involved in a number of engineering applications: in chemical absorption through structured packing, the corrugated sheets are covered by the liquid solvent, offering an enhanced interface surface between the solvent and the gas solute; in coating process, the liquid pattern influences the resulting coating quality; in condensation over finned dehumidifier, heat transfer performances are influenced by the evolving liquid layer, which may arrange as a droplets population or an ensemble of rivulets. Here, the evolution of a liquid layer flowing down an inclined plate bounded by lateral walls, which is the simplest configuration describing hydrodynamics inside structured packing, is numerically investigated. An in-house code, previously developed and largely validated in case of film instability and rivulet buildup, is used in order to solve governing lubrication equations. The full implementation of capillary pressure allows to simulate contact angles up to 60°. Film break is observed due to instability induced by lateral walls, if the imposed liquid flow rate exceeds a critical value, leading to the formation of a rivulet pattern. Fixing the size of the investigated physical domain, the number of observed rivulets, which strongly influences the resulting wetted area, is traced as a function of the flow characteristics (identified by the Bond number), the substrate wettability and the liquid properties; the corresponding bifurcation diagram is presented.


Introduction
The dynamics of driven liquid films is a complex phenomenon, which involves physical discontinuities due to moving contact lines and sudden transition between different patterns. Indeed, multiple rivulets and moving dry patches may form from a continuous film, induced by fingering instability, or neighbor droplets may coalesce and drain, forming a stable rivulet. Such film instability phenomena are largely investigated in literature [1][2][3][4][5][6][7][8]. In [1], bifurcation analysis is numerically conducted in order to analyse heterogeneous dewetting on a patterned substrate: the rupture of 1D film driven by capillary forces and viscous dissipation is investigated, with the solid substrate wettability modeled through disjoining pressure, and all the possible configurations are presented as a function of amplitude (multiplying disjoining pressure) and periodicity (pattern length). Numerical bifurcation analysis is also carried out in [2] for a 1D falling film in case of non-negligible inertia: the influence of two parameters identifying the flow condition and the wavenumber of the periodic perturbations applied to the physical system is investigated.
Knowing the evolution of a thin liquid layer over a solid substrate is crucial for many engineering applications. For example, in carbon dioxide capture through chemical absorption in structured packing, gas CO 2 flows up, while a liquid solvent falls down a collection of corrugated layers composing the packed column. Thus, the prediction of the hydrodynamics through structured packing is crucial for the design of packed columns. In fact, enhancing liquid-gas interface area ensures a more efficient mass transfer. However, increasing the solvent flow rate in order to ensure a continuous film pattern leads to high pressure drops on the contercurrent gas flow and promotes the occurrence of flooding condition (liquid plug inside packed column). Different packing configurations and solvent flow rates were experimentally investigated in [9] and the structured packing performances estimated in terms of both the height equivalent to a theoretical plate (HETP) and the pressure drop on the gas side. A sensible influence of the liquid flow rate on the pressure drop and the flooding condition occurrence was detected in [9]. Different packing configurations were also experimentally investigated in [10] and the effect of the liquid flow rate on the pressure drop presented. Due to difficulty to physically observe what happens inside the packed column, a number of numerical studies on hydrodynamics over structured packing geometries are also available. A fully 3D approach is usually adopted [11][12][13], but the inherently multiscale physics just allows to focus on a facet of the problem. In fact, packed columns may have a diameter of 10 m, while the characteristic length scale of the packing layer is 20 cm and the film thickness may be less than 1 mm. A simplified configuration, defined by an inclined plate bounded by lateral walls, was investigated through 3D computations in [11,12] and the effect of both the liquid properties and the flow characteristics was presented in terms of plate wetted area. A larger portion of packing layer was simulated by [13] using the VoF method incorporated in Ansys Fluent R and running intensive parallel computations in clusters. However, the computational costs deriving from a fully 3D model implementation allowed the authors to investigate just one packing configuration. A 2D mathematical approach was proposed by the authors [14], who numerically solved the lubrication equations. The small slope approximation, usually adopted with lubrication theory [5,[15][16][17][18], was abandoned in favor of the full implementation of capillary pressure, allowing to investigate relatively low wettable surfaces with equilibrium contact angles up to θ = 60 • [14,19]. Furthermore, a great reduction of computational costs compared to a fully 3D model was observed, without sensible loss of information about liquid pattern behavior. Here, the same approach as [14] is adopted. The in-house solver, previously developed and largely validated in case of film instability, rivulet build up, moving contact lines and droplet coalescence [14,[19][20][21][22], is parallelized through OpenMP library and used to investigate the evolution of a gravity driven film over an inclined plate bounded by lateral walls. Different plate dimensions and flow conditions (i.e. liquid properties and flow rate) are investigated. Results are presented in nondimensional form for the sake of generalization. The influence of the investigated parameters on the number of observed rivulets, which in turn affects the normalized liquid-gas interface area, is discussed. In particular, the rivulet bifurcation diagram is presented, looking for configurations leading to enhanced liquid-gas free surface area (and, thus, allowing for optimization of the absorption process through structured packing).

Mathematical model
Consider a thin film flowing down an inclined plate, driven by gravity. Let α be the plate inclination, h 0 and u 0 the undisturbed film thickness and film velocity, calculated according to Nusselt theory: Lubrication theory allows to integrate the continuity equation along the plate normal direction. Retaining the viscous dissipation terms arising from velocity gradient across the film thickness and neglecting liquid inertia, the momentum equation gives a parabolic velocity profile. Averaging the film velocity and substituting in the continuity equation gives the well known lubrication equation, with h being the film thickness. The local film pressure p is given by the hydrostatic, capillary and disjoining contributions, with d h 0 being the precursor film thickness and θ the equilibrium contact angle. The small slope approximation, strictly valid for low free surface slope and usually adopted in literature, is here abandoned and the free surface curvature is estimated according to [14,19,21].
Eq. (6) was largely validated in [14] for contact angles up to 60 • . Defining the following dimensionless quantities and introducing the Bond number the governing equations, Eqs. (2) and (3), are recast in a nondimensional form: Eqs. (9) and (10) are numerically solved on a structured, orthogonal grid. An in-house solver, developed in FORTRAN by the authors and largely validated [14,19,21], is used. Space and time discretizations are explained in detail in [14,23]. In order to speed up simulations, the source code is parallelized through OpenMP library for a shared memory machine, allowing for a sensible reduction of computational costs. In particular, the two algebraic, pentadiagonal systems to be solved at each time step and deriving from ADI approximate factorization, are respectively decomposed in n x and n y independent subsystems, that are solved in parallel, with n x and n y being the number of elements of the implemented structured grid along Cartesian directions.

Problem setup
The dynamics of a gravity driven film flowing down an inclined plate of nondimensional size L X × L Y , bounded by lateral walls, is investigated. However, only half of the plate is simulated along Y direction due to symmetry of the problem, in order to reduce computational cost. The influence of both the Bond number, Eq. (8), and the dimensionless plate width L X on the flow pattern evolution is investigated. The undisturbed film flow rate is supposed to enter the computational domain through the top section of the with Q Q Q being the nondimensional film flux, Q Q Q = − ∇P 3 H 3 . Fully developed flow condition, ∇P ·n n n = 0, ∇H ·n n n = 0, is imposed through the bottomest section of the plate, Y = L Y . Symmetry condition, Q Q Q ·n n n = 0, ∇H ·n n n = 0, is imposed through the middle section of the plate, X = 0, while wall boundary condition is implemented through X = L X /2, leading to: The vertical plate is initially dry. Thus, the precursor film thickness covers the whole computational domain, H| T =0 = δ , at the first stage of the computation. The dry initial condition allows us to study contact line dynamics and phenomena related to finger instability of a falling film.

Result
Distillation and absorption processes through structured packing are very complex, multiphysics phenomena, with many length scales involved. However, a gravity driven film down an inclined plate is the simplest configuration for better understanding the hydrodynamics inside packed columns [11,12]. Thus, a film flowing down a vertical plate, α = 90 • , is chosen as a representative setup. The instability induced by the lateral walls determines the formation of a rivulets array, with dry regions observed over the computational domain [12,14]. However, different configurations may occur, in terms of rivulets number and, thus, interface area, as observed in [14]. Thus, a parametric analysis is conducted, in order to trace the flow regime transitions. In particular, the influence of both the Bond number, Bo = ρ g h 0 2 /σ , which is the ratio between gravitational and surface tension forces, and the nondimensional plate width L X is analyzed. However, results are presented in terms of Bo and L X Bo 1/2 , the latter being the ratio between the dimensional plate width L x and the capillary length, L c = σ /(ρ g): Physically, fixing L X Bo 1/2 and varying the Bond number can be interpreted as fixing both the physical domain dimension and the liquid properties, while varying the liquid flow rate through the plate with Bo. However, the following results can be extended to arbitrary fluid properties and physical plate dimensions and, thus, are valid for a wide range of applications, since dimensionless parameters are used in the analysis. The equilibrium contact angle, which also affects the contact line dynamics, is set to a fixed value of θ = 60 • , giving a solid surface that can be considered low wettable in the context of lubrication theory, while the precursor film thickness and the disjoining parameters are set to δ = 5 × 10 −2 and n = 3, m = 2 according to [5,14,15]. The contact angle hysteresis, which is implicitly modeled by disjoining pressure, is checked running preliminary computations of a 1D falling film. Reduction to 1D mathematical problem is obtained through imposition of symmetry condition at X = 0 and X = L X . Film rupture (occurring when the falling 1D film breaks forming a train of self similar droplets, as reported in Fig. 1) is always observed when Bo ≤ 5.01 × 10 −2 . The contact angle hysteresis, estimated for Bo = 5.01 × 10 −2 , is given by receding and advancing contact angles of a self similar droplet: θ rec = 41.3 • and θ adv = 65.5 • . Moreover, it was verified that a spatial discretization step equal to ∆X = 10 −1 ensures grid independency, since refining the mesh to ∆X = 5 × 10 −2 leads to a less than 1% uncertainty in the computed θ rec , θ adv . Thus, ∆X, ∆Y ≤ 10 −1 is set for all the 2D computations. and shed away. Thus, the corresponding rivulets, observed in the transient film distribution, shed as well and disappear in the stationary solution. Note that residual droplets may form during shedding of unstable dry patches, Fig. 3(b); a similar phenomenon is experimentally and numerically reported in literature [14,24]. This further bifurcation, related to dry patch stability, explains the trend of the number of stationary rivulets, that may both increase or decrease with a small variation in the Bond number.  Fig. 6(a); however, 2 rivulets disappear from transient to stationary solution due to instability when Bo < 0.12 (thus, 3 stationary rivulets are observed instead of 5), as reported in Fig. 6(b).
Since the maximization of the interface area between liquid solvent and gas mixture is crucial in absorption process through structured packing, the normalized interface area, resulting from the stationary film distribution, is also computed over the wetted domain (i.e. where H > δ ) and plotted in Fig. 7 as a function of the Bond number. Note that the normalized interface area increases for increasing Bo, until the bifurcation point is reached; in fact, a discontinuity, with a sudden decrease of ξ I , is observed. Thus, local maximum points are observed, at different values of the Bond number, for each plate width. The maximum value of the normalized interface area, ξ I 0.49, is obtained in Fig. 7(a) for Bo = 0.16 and L X Bo 1/2 = 25.0, corresponding to the 5 rivulets solution family. Thus, an optimal configuration in terms of plate width and Bond number can be identified, in order to maximize the liquid-gas interface area. Physically, this means that, for a given liquid solvent, both the optimal corrugation frequency of the structured packing geometry and the optimal liquid flow rate can be determined, in order to enhance the absorption process. It is important to point out that further increasing the Bond number (i.e. the solvent flow rate) leads to the fully wetted condition. However, this may not correspond to the optimal configuration, due to: higher pressure drop of the gas flow, induced by a reduction of the available void volume, through which the gas mixture flows; reduction of the loading and flooding gas velocities.

Conclusion
Instability of a falling film over a vertical plate bounded by lateral walls was numerically investigated. The effect of nondimensional plate dimension and Bond number on the resulting pattern was analysed in terms of number of observed rivulets and normalized free surface area. Two different bifurcations were observed: depending on the imposed parameters, namely L X and Bo, from 3 to 5 rivulets were observed during the transient solution; however, if some of the generated dry patches are not stable, they sheds away and the number of stationary rivulets reduces. The influence of bifurcations on the interface area were also discussed, showing that, fixing the physical plate dimension and the liquid properties (fixed L X Bo 1/2 ), a local maximum of the free surface area can be obtained at low flow rate (at low Bond number). Since the investigated test case is the simplest configuration describing hydrodynamics in structured packing, the presented results can be used as a novel approach to test the packing geometry and enhance absorption process, that requires increased liquid-gas interface area at low solvent flow rate in order to promote mass transfer and reduce the pressure drops of the gas flow. Equilibrium contact angle equal to 60 • was successfully investigated due to the full implementation of the capillary pressure, while it was demonstrated in [25] that small slope approximation leads to inaccurate results for θ ≥ 30 • . Furthermore, the effect of lateral walls on the film stability was never investigated in literature involving lubrication theory. As a future work, the solver should be coupled with the core flow solver, in order to investigate the effect of shear, applied at the free surface by the gas flow. The effect of liquid inertia should be also included in the mathematical model.