Lattice Boltzmann investigation of droplet inertial spreading on various porous surfaces

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers. L’archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d’enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.


I. INTRODUCTION
The spreading of liquid drops on solid surfaces is of both fundamental and industrial interest. In many applications, surfaces are porous or covered with a thin porous layer. The presence of a porous layer modifies the wettability of the substrate and, hence, the spreading patterns [1,2]. The spreading of liquid drops on porous media is a wide-spread phenomenon [3] involved in natural and industrial processes such as ink jet printing [4], adhesion, coating [5,6], lubrication, detergency, plant treatment, composite manufacturing, painting, and oil recovery. This topic has received much attention from researchers [7][8][9][10][11][12][13]. To the best of our knowledge, the understanding of the underlying mechanisms of drop spreading on a porous surface remains limited, in particular with respect to spreading on a smooth surface. Earlier investigations were mainly devoted to viscous wetting, and most modeling and simulation approaches were performed within the framework of a continuous porous material, pore-scale approaches being quite rare [14][15][16]. Long-term spreading dynamics, arising from a competition between capillarity and viscosity both inside and outside the porous material, are preceded by an initially faster spreading. The spherical cap-shaped droplet approximation fails to describe this short-term spreading.
Recently, both the inertial spreading dynamics and capillary invasion of pores were investigated by means of numerical simulations [15]. However, the dependence of spreading dynamics on the porosity remains quite qualitative in this study. The present paper aims to address this point with the * xavier.frank@supagro.inra.fr help of lattice Boltzmann (LB) numerical experiments. We will focus on simple porous substrate models, in which pore space is made of parallel holes of infinite length.
In the present study, gas-liquid flows are simulated with the help of the well-known Shan-Chen pseudopotential approach [42]. The velocity used to compute equilibrium distributions is modified according to a local fluid-fluid force that emerges from interactions between neighboring fluid particles. A liquid phase and a vapor phase emerge spontaneously within this framework. Liquid density ρ L , vapor density ρ V , and surface tension σ depend on the magnitude of fluid attraction in the Chen-Shan model. An interface is located at the middle of the transition from ρ G to ρ L , where density ρ reaches the value ρ M = (ρ G + ρ L )/2.
In this study, Martys and Chen's fluid-solid pseudopotential [43] describes interactions between fluid particles and the solid. The equilibrium contact angle θ eq s depends on the magnitude of both fluid-fluid and fluid-solid local interactions in the Chen-Shan model. Please see the previous works [24,42,43] for more details about the numerical approach.

III. NUMERICAL SETUP
We employ symmetries through the xz and yz planes to reduce the demand for numerical resources (Fig. 1, panel A). As a consequence, the initial position of the droplet center has to be placed on the z axis. The initial z position of the droplet center of mass is fixed to place the bottom of the liquid phase at themiddleofz range. Simulation box dimensions are fixed to 180 × 180 × 450, and the initial droplet diameter is set to 190. Simulations are conducted without gravity, and the wettability of the solid is defined using the solid equilibrium contact angle θ eq s . First, the droplet is equilibrated during 10 000 LB iterations without any solid node in the simulation box. The equilibrated droplet radius R, which is slightly different from the initial radius, has to be measured at the final state of this equilibrating phase (Fig. 1, panel B1). Then, the solid nodes are introduced in the simulation box, with a solid position adjusted in such a way that the droplet interface slightly touches the solid phase. To compute the radius of the droplet's footprint r (Fig. 1, panel B2), a density profile is extracted from a z,x = const. line, just above the solid surface. This specific position is chosen to extract the density profile only above the solid, and not above the pores. Wetting simulation runs during 15 000 iterations, and the spreading radius r is computed and saved at each time step.

IV. SPREADING ON FLAT SURFACES
First, numerical spreading experiments are conducted with flat surfaces. As spreading dynamics can depend significantly on liquid viscosity μ L [44], it is crucial to compare the magnitude of inertial and capillary forces with that of viscous forces. In the present work, the value of the Laplace number La = ρ L σR μ L 2 is La = 135 (in other words, the Ohnesorge number is Oh = 1/ √ La = 8.6 × 10 −2 ), compared with the La values in the Bird et al. study [45], which ranged from 130 to 72 000. As La ≫ 1 we can then consider that the droplet is of low viscosity in our simulations. In such a case, the three following steps can be identified within droplet spreading dynamics [46]: a first step, independent of wettability, in which the spreading radius scales as r ≈ t 1/2 , a second step, in which the spreading dynamics depend on the solid equilibrium contact angle [45][46][47], and, finally, a slow viscous step [48]. The first step is not observed in our numerical simulations, as a consequence of grid resolution limits, and only the intermediate regime is studied in the present work. Under such conditions, the inertial dynamics of droplet spreading is dominated by the propagation of a capillary wave from the bottom to the top of the droplet interface [45]. As a consequence, the time scale τ of inertial spreading is τ = ρ L R 3 /σ . The normalized spreading radius r/R is deduced from r using R as the length scale, and it is plotted as a function of t/τ for 5 values of the equilibrium contact angle θ eq s (Fig. 2). Strictly speaking, only the first inertial regime ( t τ 0.04) can be characterized by a power law [46], and the intermediate regime droplet spreading dynamics in this intermediate phase [45]. This power law approximation is expressed as where α is the power-law exponent and C is the prefactor, which is in good agreement with the experimental results [45]. The evolution of the power law prefactor C(θ eq s ) and exponent α(θ eq s ) with θ eq s , shown in Fig. 2 as insets, has to be underlined and compared with earlier experimental results. Both decrease with θ eq s , which is in agreement with the experiments from Bird et al. [45]. The duration of inertial spreading [49] can be easily calculated from τ as T ≈ (ρ L σR/μ L 2 ) 1/8 τ .B i r d et al. showed that the duration of the inertial regime does not depend significantly on the equilibrium contact angle [45]. In the present study, the evaluation leads to T ≈ 1.84τ , which is in good agreement with our numerical simulations (Fig. 2).

V. SPREADING ON POROUS SURFACES
We consider now surfaces perforated by a regular pattern (square or centered square lattice) of holes of infinite length. Between neighboring cylindrical pores, solid walls are effective to avoid any liquid communication. Surface porosity is defined as ǫ = S pores /λ 2 ; S pores is the sum of the transversal sections of all pores belonging to one lattice cell, and λ is the lattice period. The higher is the porosity, the lower is the driving force of spreading. As a consequence, the spreading dynamics are expected to slow down as ǫ increases. However, the porosity of a surface can be distributed in various ways: small pores with high numerical density or larger pores with lower numerical density, or various pore lattices and various pore section shapes. Do these details matter or not? If this is not the case, the spreading dynamics will depend only on ǫ regardless of the pore section shape, the pore lattice unit cell pattern, or the pore lattice period. To clarify this point, porous surfaces are built up using various geometrical parameters for both square and centered square lattices of circle-sectioned pores, and the square lattice of square-sectioned pores (Fig. 3).
For each kind of lattice, the geometry of the porous surface is completely defined by the pore size d pore (diameter for circlesectioned ones and side for square-sectioned ones) and spatial period λ. Porosity ǫ is easily computed from d pore and λ. Surface parameters are summarized in Table I. As an example, a perspective view of one solid phase is shown in Fig. 4.
Numerical spreading experiments are conducted in the same way as for flat surfaces. The surface porosity leads to TABLE I. Parameters of pore lattices: pore lattice type, pore size d pore , and pore period λ, which is deduced the surface porosity ǫ.Both d pore and λ in lattice units (l.u.).  (Fig. 5). Liquid spreads on the solid surface and crosses the pores. Imbibition can start after the pore opening has been covered by the liquid interface. As expected, pores are invaded when θ eq s < 90 • , not when θ eq s > 90 • (Fig. 6). Both droplet height evolution and droplet volume variation emerge from a competing mechanism between spreading on the solid phase of the surface and capillary imbibition inside the porous medium. Such a point was detailed in an earlier paper [15] and is not developed here.
The simplest approach to consider both the contact angle and morphology of a complex surface in a single law is to assume that the surface behaves as an equivalent smooth surface. Such a description leads to the famous Cassie-Baxter law [50]. Considering a composite surface involving two materials with two different equilibrium contact angles, the surface is supposed to exhibit an effective homogeneous equilibrium contact angle θ eq eff , depending on both the intrinsic equilibrium contact angle and the proportion of basic materials. In the present case, we have a solid phase with an equilibrium contact angle θ eq s and a surface fraction (1 − ǫ), and a pore phase with equilibrium contact angle θ eq pore and surface fraction ǫ. Using these assumptions, an effective equilibrium contact angle θ eq eff (θ eq s ,ǫ) could be deduced from the Cassie-Baxter law as follows: cos θ eq eff = (1 − ǫ) cos θ eq s + ǫ cos θ eq pore .
Effective power-law parameters can stem from the computed values of θ eq eff , as if the porous surface is an equivalent smooth solid surface with equilibrium contact angle θ eq eff : α eff θ eq s ,ǫ = α θ eq eff .
The application of the Cassie-Baxter law to our case requires θ eq pore to be fixed. A simple hypothesis can then be proposed: For each of the 12 pore lattices, 5 values of θ eq s were used (cf. Fig. 2), which means a total of 60 curves. as long as a pore is not invaded by liquid, it remains an obstacle to global spreading. Thus, θ eq pore = 180 • seemstobea straightforward hypothesis.
To compare our effective law with the numerical results, we plot the ratio r/R as a function of C eff (t/τ) α eff , t/τ being a curve parameter (Fig. 8). Concretely, for each simulation, θ eq eff is computed as a function of both θ eq s and ǫ, and effective power-law parameters C eff and α eff are deduced from linear fits of values extracted from the numerical simulation of droplet spreading on a smooth surface (Fig. 2). Despite oscillations, it is clear that the spreading radius dynamics from numerical simulations are similar to a power-law whose prefactor and exponent are predicted by the Cassie-Baxter law. Moreover, various curves fit this simple effective law well, although the proposed spreading law does not explicitly take into account the distributed pattern of porosity on the surface. Only the pore fraction matters whatever its distribution. The global behavior of the droplet is very close to that on an equivalent smooth surface.
Recently, Stapelbroek et al. conducted droplet spreading experiments with various complex solid surfaces, such as microtextured and chemically striped substrates [51]. In this study, the data collapse of the crossover time between the early inertial spreading and intermediate regime could be achieved using only a single parameter, the effective contact angle, regardless of the details of the substrate. Such result corroborates our findings. To overcome the energy barriers in the presence of topographic or chemical patterns that lead to a spatial variation of the local surface energy, an additional source of energy is required in the vicinity of the contact line. This energy must come from the kinetic energy of the flow inside the droplet. As a consequence, it could be reasonably argued that the balance between the kinetic energy of the spreading drop and the surface energies is a relevant mechanism, and that the details of a complex surface play a minor role in the inertial spreading global dynamics. Of course, such a picture is no longer valid in the viscous regime [47].

VI. CONCLUSIONS
In summary, we have proposed a droplet spreading law on porous surfaces, including both the intrinsic solid contact angle and surface porosity. Lattice Boltzmann numerical simulations were performed for smooth surfaces with various equilibrium contact angles and the predicted power law for spreading dynamics compares favorably with published experimental results. Simulations with various porous surfaces and solid equilibrium contact angles show that increasing porosity slows down spreading in the same way as increasing the equilibrium contact angle for smooth surfaces. Assuming that pores are perfectly nonwetting patches, we deduced an effective equilibrium contact angle of the surface, and both the power law prefactor and exponent from an effective equilibrium contact angle are in satisfactory agreement with the simulation.