Numerical Simulation of Droplet Impacting and Sliding on Hydrophobic Granular Surfaces

Bao, Qing; Kang, Hengyi

doi:https://doi.org/10.1155/2021/5534003

Mathematical Problems in Engineering

On this page

Abstract Introduction Results and Discussion Conclusions Data Availability Conflicts of Interest Acknowledgments References Copyright Related Articles

Research Article Corrigendum

!

A Corrigendum for this article has been published. To view the article details, please click the ‘Corrigendum’ tab above.

Research Article | Open Access

Volume 2021 | Article ID 5534003 | https://doi.org/10.1155/2021/5534003

Numerical Simulation of Droplet Impacting and Sliding on Hydrophobic Granular Surfaces

Qing Bao¹and Hengyi Kang²

Academic Editor: Debao Zhou

Received11 Jan 2021

Revised14 Mar 2021

Accepted31 Mar 2021

Published14 Apr 2021

Abstract

Droplet sliding naturally happens with practical significance in developing artificial self-cleaning surfaces or impermeable barriers. On water-repellent soil surfaces, such processes evolve at very small scales, typically at the particle level. To address this, this paper presents a two-dimensional Lattice Boltzmann (LB) study on the droplet sliding dynamics on a layer of regularly arranged particles with varying size and contact angle (CA) aimed at mimicking conditions comparable to those of real soils. The numerical droplet is initialized above the inclined granular surface with different lifting distances and deposited by gravity. The droplet hits the surface with different impacting velocities and subsequently slides down the slope. Four droplet-sliding behaviors were observed: a droplet sticks to the granular surface, a droplet moves by pinning and depinning of its interface (“stick-slip”), a droplet undergoes periodic elongation and shortening during sliding, and a droplet lifts off the granular surface and may be ruptured. For a droplet that displays the “stick-slip” behavior, the sliding velocity reaches a converged terminal velocity, which increases with a higher CA, a more inclined slope, and a smaller particle size. However, nonunique terminal velocities were identified to be affected by the impacting velocities, but their correlation is not continuous and may not be positive. Finally, we propose to quantify the rotational or translational movement by effective kinematic ratio (EKR), which is defined as the translational kinematic energy divided by the total kinematic energy. The unique relation between the EKR and the terminal velocity is suggested to be one practical indicator to intrinsically characterize the water repellency at the particle level.

1. Introduction

Droplet sliding on an inclined surface has practical significance in the development of artificial synthesized self-cleaning surfaces [1] such as antifouling clothes and kitchens [2]. In nature, wildfires [3] or soil microbiological activities [4] may also generate organic hydrophobic coatings on the soil particles, reducing the chemical affinity between the soil and the water, and thus making the soil waterproof. The waterproof property may be utilized for reducing the amount of water infiltration to improve the stability of slopes or landfills during rainfall. From the microscopic perspective, a droplet falling on an inclined soil surface will travel downslope, wherein the dynamics reveals the resistance of the soil to water infiltration. To better understand the fundamental mechanism behind, the paper utilizes an advanced numerical tool that can capture the evolution of droplet interfaces attached to one or several particles, and provides an insight into the droplet interaction and sliding dynamics with granular soil surfaces.

On an idealized, smooth, and homogeneous surface with a nonzero slope angle, the droplet slides down with a steady translational movement [5]. However, a realistic surface has microscale bulges or indentations, i.e., the rough substrate. Particularly, for granular materials, such as soils, the surface morphology is related to the size, shape, and the alignment of the particles [6]. In contrast to the idealized smooth surface, previous studies [7–9] indicate that a droplet is stuck on a rough surface with an inclined angle lower than a critical value, i.e., the sliding angle [8]. If the inclined angle exceeds the sliding angle, the droplet slides by regularly pinning and depinning of its interfaces on the substrate [10]. Similar stick-slip behaviors are also observed for a droplet sliding on a geometrically smooth surface with patterned varying wettability [11–13]. The stick-slip behavior triggers the oscillation of the droplet interfaces, resulting in a viscous dissipation. At the same time, gravitational potential energy contributes to the downslope movement. It is found that the stick-slip movement can reach a steady state with a constant sliding velocity [14], indicating that a balance is achieved between the viscous dissipation and the gravitational energy input. The constant sliding velocity, usually termed as the terminal velocity, is the focus of the studies on the droplet sliding dynamics.

Hyväluoma et al. [5] simulated the sliding dynamics for a droplet in the Wenzel state on a hydrophobic rough surface patterned by ridges. If a droplet maintains nearly a rounded shape during sliding, the linear relationship between the capillary number , and the Bond number proposed in [15] is verified, written as Here, is the slope angle and is the gravitational acceleration. and indicate the dimensionless values of sliding velocity and the downslope gravitational driven force, respectively. The cutoff value reflects the minimum driven force in the tangential direction to the slope that can trigger the droplet sliding. Effects of wettability are also investigated. For a water droplet on a flat surface, the intersection angle between the interface and the base plane, i.e., the contact angle (CA) quantifies the wettability. A surface with a low CA smaller than 90° for the water droplet is considered as hydrophilic or wettable, while a high CA larger than 90° is considered as hydrophobic or water-repellent. In [5], by increasing the CA, decreases, and the ratio increases, indicating that a higher CA can facilitate the droplet sliding. They also found that it is easier for a droplet to move parallel than perpendicular to the patterned ridges. Wang et al. [16] numerically studied the effects of the solid fraction and the width of the micropillar for a droplet sliding in the Cassie-Baxter state on a patterned surface. They found that increasing the spacing between the micropillars (i.e., decreasing the solid fraction) is beneficial for a droplet sliding down. This is consistent with the experimental study in [8] that a larger slope angle is needed to initiate the droplet sliding for the surface with a higher solid fraction. However, further increasing the spacing between the micropillars may trigger the droplet configuration to transit from the Cassie-Baxter state to the Wenzel state [17, 18]. In other words, the droplet penetrates the microscale indentations and the sliding is impeded. This indicates that there exists an optimal solid fraction for the patterned surface, on which the droplet achieves the maximum terminal velocity. Varagnolo et al. [13] used a high-speed camera to trace a droplet sliding on a geometrically smooth surface with alternating hydrophilic and hydrophobic strips. They found that both the CA and the surface pattern of the domains with varying wettability affect the terminal velocity. The droplet undergoes a stick-slip motion and the terminal velocity is larger with a higher apparent CA. Zhao et al. [18] and Sbragaglia et al. [11] found that terminal velocity increases monotonically with the fraction of the hydrophobic regions. Concerning the detailed patterns, Varagnolo et al. [19] found that the terminal velocity depends mainly on the apparent CA and whether actual surface pattern is in the shape of strips, squares, or triangles has very limited effect.

The above studies on the droplet sliding initiated the downslope movement by gradually increasing the slope angle of the tested surface, on which a droplet shall be deposited beforehand. In the real world, droplets may fall on an inclined surface with an impacting velocity. However, researches on the droplet impacting and rebounding dynamics only deal with the droplet distortion and the possible rupture [20, 21] but ignore the subsequent sliding behavior. Up to now, the dependency of the terminal velocity on the impacting velocity remains unclear for a droplet, which can be clarified by conducting pore-scale numerical analysis. In the present work, the phase-field Lattice Boltzmann (LB) method is utilized as the numerical tool, due to its capability of capturing the evolution of numerous liquid-vapor interfaces interacting with irregular solid boundaries. On the contrary, traditional multiphase computational fluid dynamics (CFD) methods, e.g., the volume of fluid (VoF) method [22, 23] and the level set method [24, 25], can only trace the dynamics of a few predefined free surfaces in hydraulic problems, such as the open-channel flow. Another advantage is that surface tension and CA are accurately incorporated into the phase-field LB method by a thermodynamic approach. The feature is vital for accurately solving the dynamics of microfluids, where the capillary forces play an important role. By simulating the droplet impacting and sliding on inclined hydrophobic granular surfaces with varying CA and particle, various modes of droplet dynamics are identified with the sliding velocity monitored. In addition, we further propose an energy-based method to quantitatively characterize the rotational or translational type of droplet movement, by taking the advantages of the LB method that provides not only the global movement of the droplet but also the flow field inside the droplet.

2. Phase-Field Lattice Boltzmann Method

The numerical simulations of droplet sliding dynamics on a layer of particles need to deal with the movement of air-water interfaces on granular surfaces with controlled wettability. Such problems can be handled by adopting different multiphase LB methods, e.g., pseudopotential LB methods [26, 27], free energy models [28–30], and phase-field LB method [31, 32], noting that the LB method, originally developed from gas kinematics, is well-known for the efficient parallelization of the algorithm and flexibility in dealing with irregular boundaries. Therefore, multiphase LB methods have been widely applied in simulating pore-scale flows, such as unsaturated permeability [33], wetting-drying in porous media [34], capillary rise [35], Pickering emulsion [36], and colloid hydrodynamics [37].

In the present paper, the simulation is performed with the phase-field LB method [32]. The system of binary fluids adopts the one-fluid formulation where a phase index field C is defined to describe the phase that occupies each computational cell. The phase index C is around one for the liquid phase and is around zero for the gas phase. The value gradually changes from zero to one in the interface. The evolution of the phase field is governed by the minimization of the free energy given bywhere is the total free energy, and the subscripts b and s represent the contributions from the bulk fluids and the solid boundaries, respectively. Ψ_s is nonzero only for the fluid nodes adjacent to the solid nodes and C_s is the composition at a solid surface and ϕ_s (i.e., ϕ₁, ϕ₂, ϕ₃, etc.) are constant coefficients. consists of two parts, namely, the local bulk energy and the interface-induced energy The constants β and κ can be related to the surface tension σ and interface thickness ξ by

Minimizing the bulk free-energy Ψ_b for the fluid phase leads to the Cahn–Hilliard (CH) equation:where u is the velocity, µ is the chemical potential, and M is a numerical parameter controlling the rate of diffusion. The chemical potential is derived from the free energy as follows:where µ₀ is the local part of the chemical potential, and is the interface-induced chemical potential. Since E₀ takes the form of the double-well potential function, C will be close to either 0 or 1 in the pure single-phase region. The negative sign in the nonlocal term helps build a stable diffusive gas-liquid interface. The two competing terms mean that the phase index C is stable either completely separated as the pure single phase or well mixed as the interface. Therefore, the evolution of the CH equation results in the phase separation while the total free energy in the system is decreasing.

For the fluid nodes adjacent to the solid surface, minimizing the Ψ_b + Ψ_s leads to a partially wetting boundary condition, where the CA can be defined at the fluid-solid interface by specifying the gradient of C [38]. The boundary condition is established aswhere θ^eq is the equilibrium CA at the solid surface.

Solving the CH equation (equation (4)) provides the spatial distributions of density ρ, the kinematic viscosity ν, and the surface tension force F_s. Here, density and viscosity are both necessary material properties for solving the incompressible Navier–Stokes (NS) equations. They can be calculated from the index C simply through linear mapping, which are written aswhere the subscripts l and indicate the liquid phase and the gas phase, respectively. The surface tension force F_s applied at the interfacial region is proportional to the gradient of the interface-induced chemical potential, given by

The above equation is based on the chemical potential is called the “potential form” [31]. Together with the standard isotropic spatial discretion, the nonphysical velocities at the interfacial region in the multiphase LB method (also termed as “spurious currents”) can be minimized [39].

With the information provided by the CH equation, the LB method computes the evolution of the mesoscale particle probability functions, from which the macroscale velocity and pressure are obtained. The standard discrete Boltzmann equation (DBE) with the body force term for macroscopic motion of fluid flows iswhere is the distribution function, e_α is the discrete velocity for the direction α, c_s is the lattice speed of sound, and λ is the relaxation time, which is related to the kinematic viscosity by . F is the body force term , including the contributions from surface tension force and gravity. , and is the equilibrium distribution function defined bywhere t_α is the weight for the direction α.

The LB formulation has been proved to be equivalent to the artificial compression method in solving the incompressible NS equation by the Chapman–Enskog analysis [40]. However, the original formulation stores the base density with the hydrodynamic pressure as a small perturbation in the distribution functions and applies only to the multiphase problem with a nonuniform density field. The problem is solved by defining a new distribution function to split the density [41]. Similarly, the is transformed as . The DBE for the new distribution function is given as follows:from which the evolution of pressure and velocity in the macroscopic NS equations can be recovered as

By integrating the DBE (equation (12)) with a trapezoidal approximation over a small time step δt [32], the Lattice Boltzmann equation for numerical implementation iswhere . With the temporal discretion, the truncation error is O(∆t²). The velocity and the macroscopic pressure are obtained from the distribution functions by

The same temporal discretion can also be applied to solve the CH equation [31, 32]. However, the LB discretion can only recover the advection term in equation (4), while the nonlinear diffusion term cannot be recovered and related to the relaxation time λ. The advantages of the LB scheme to recover both the advection term and the linear diffusion as in solving the NS equations cannot be exploited. Instead, the shall be included as an additional term, which still relies on the finite difference approximation with first-order accuracy. Therefore, the CH equation can be more efficiently solved by a full finite difference (FD) scheme [42] and then coupled with the NS equations being solved by LB. We also adopt this hybrid LB-FD formulation. The advective term is discretized by a second-order GVC (Group Velocity Control) scheme. In the Cartesian coordinate system with a uniform grid, the nonlinear diffusion term is approximated with the isotropic central difference scheme involving all the neighborhood nodes [42]. The third-order TVD (Total Variation Diminishing) Runge–Kutta scheme is used for the temporal discretion [43]. The computational details for solving the CH equation can be referred to computational fluid dynamics textbooks.

3. Numerical Test

The dynamic behavior of millimeter-sized droplet impacting and subsequently sliding on inclined granular surfaces is simulated with the LB method described in the previous section. The numerical tool was validated by both the benchmark numerical tests and experimental tests, which can be referred to our previous work [44]. As it takes a significantly large number of time steps for a droplet to travel a long distance before reaching a steady state, a two-dimensional analysis is performed, so that the computation can be finished within a reasonable time. The schematic diagram of the numerical setup is illustrated in Figure 1. The droplet diameter is set to 96 in LB units to mimic a 10-uL droplet, which is commonly used in the sessile drop method for quantifying the wettability of a granular surface [45]. The substrate is a layer of hydrophobic particles with varying CA, fixed on an inclined surface at different slope angles from 5° to 75°. The lattice nodes occupied by the particles are marked as solid nodes, implemented with the bounce-back boundary condition. The particle diameters are 24, 32, and 40 in LB units, corresponding to the physical values 0.67, 0.89, and 1.15 mm, respectively. The sizes are all within the range of coarse-grained sands, where the ratios D_par/D_drop are 1 : 4, 1 : 3, and 1 : 2.4. The relatively larger values of D_par/D_drop indicate that the details of the substrate geometry shall have a significant effect on the droplet dynamics, and hence a pore-scale numerical study is required. To facilitate pore fluid flows [46, 47] in a two-dimensional analysis, artificial spacings are created to separate the particles. In all the simulations, the spacing ratio ϵ₀, defined as the gap divided by the distance between two adjacent particles, is fixed as 0.25. The gap between two adjacent particles then approximately equals the equivalent circular diameter of the pores in a square packing. The CAs used in this work are 104.5°, 113.6°, 123.4°, and 134.4°, with the cosine values of −0.25, −0.4, −0.55, and −0.7, respectively. Only the hydrophobic particles with CA > 90° are studied, which ensures that the droplet will not infiltrate into the granular surface.

For each particle layer with different slope angles, the droplet hits the granular surface with multiple initial velocities. The impacting velocities are controlled by releasing the droplet from different distances above the granular surface so that the droplet accelerates under the gravitation force. The gravity is set as = 4.27 × 10⁻⁷ in the LB unit, corresponding to 9.81 m/s² in the physical unit. The falling distances are chosen as 1/6, 2/3, 1, 2, 3, 4, 5 times the droplet diameter , covering a wide range within the low-velocity impacts. Neither splitting nor splashing shall be induced. For the largest falling distance, the Weber number W, defined as the ratio of the inertial force over the surface tension , is 9.84. The value is large to induce severe droplet deformation [24, 25]. On the other hand, the smallest Weber number is only 0.328 for , indicating a negligible inertial effect. The droplet almost remains circular during its deposition on the granular surface.

In each test, the liquid density is ρ_l = 1, the liquid kinematic viscosity is ν_l = 3.33 × 10⁻³, and the surface tension is σ = 4 × 10⁻³ in LB units, corresponding to 10³ kg/m³, 2.34 × 10⁻⁶ m² s⁻¹ and 72.8 mN/m, respectively, in physical units. To better illustrate the relationship between the LB units and the physical units, the input parameters are summarized in Table 1. Here, the kinematic viscosity is close to the value of 20% concentration saline water, which is 2.33 times that of pure water (1.006 × 10⁻⁶ m² s⁻¹). Further reducing the lattice viscosity ν_l would trigger numerical instability, so a smaller physical viscosity can only be achieved by enlarging the computational domain, which would demand the exponential growth of the computational power. Despite the compromise, in terms of relative magnitude of inertial and capillary forces to viscous forces, the dimensionless Laplace number is sufficiently large so that the droplet is considered as having low viscosity [48]. Such a large value has been seldom numerically investigated in the literature.

The simulation involves a preparation stage and a gravity deposition stage. First, the droplet is initiated in a circular shape above the particle layer with a certain distance to control the impacting velocity. The phase index C is predefined as 1 within the droplet and is predefined as 0 in the surrounding air. At the interface, a smooth transition is realized using a tanh function [32]. However, the distribution of the phase index C is not in the equilibrium state. The correct pressure field that fits the surface tension forces has not yet been built up. Hence, the simulation runs with no gravity until the system reaches equilibrium. Then, gravity can be applied to the liquid phase so that the droplet moves downward and impacts the granular surface. At specific time intervals, the subsequent evolution process is recorded in the form of two-dimensional image files, showing the information of both the phase distribution and the velocity field. During the simulation, pressure boundary conditions are applied at the top and the bottom boundaries to maintain the ambient air pressure as constant. Periodic boundary conditions are imposed on the other two sides. With a limited computational domain, it is possible for a droplet to travel a sufficiently long distance to reach a steady state. As denoted by the dashed line in Figure 1, we notice that a smaller slope angle and a larger falling height require a larger computational domain. Hence, the width in the y dimension n_y is determined by , plus the space to contain the particle layer and the gaps at the top and the bottom boundaries. The length n_x in the sliding direction is defined as 480, which is 5 times the droplet diameter. The height n_y is controlled by both the slope angle and the droplet falling height. The largest size is 480 × 720 for a 5° slope with falling height, while smallest size is 480 × 250 for a 75° slope with falling height. For each configuration of CA and particle size, a complete sequence of 105 numerical tests with 15 slope angles and 7 impacting velocities require around 96 hours on an AMD R9 390 graphics processing unit (GPU) to finish the computation.

4. Results and Discussion

4.1. Droplet Sliding Dynamics on Granular Surfaces

For clarity, a test is labelled by a combination of three segments. The first segment is the particle size and the intrinsic CA. The particle layer with the diameters 0.67 mm, 0.89 mm, and 1.15 mm is denoted by S, M, and L, respectively. The second segment, which begins with “sl,” is the slope angle. The last segment, which begins with “h,” is the falling height divided by the droplet diameter . For example, “h1/6” means that the falling height is one-sixth of the droplet diameter and “h5” means 5 times the droplet diameter. Therefore, “M113_sl60_h5” indicates that the droplet is falling from 5 times its diameter on a layer of mono-sized 0.89 mm spheres with an intrinsic CA of 113°, which is inclined at the slope angle 60°.

In the present paper, four modes for the droplet sliding on the inclined hydrophobic granular surfaces are identified: a droplet sticks to the granular surface (mode I), a droplet moves by pinning and depinning of its interface (mode II), a droplet undergoes periodic elongation and shortening during sliding (mode III), and a droplet lifts off the granular surface and may be ruptured (mode IV). Figure 2 illustrates the temporal evolution of the sliding velocity of the droplet for L104_sl60_h5, M113_sl60_h5, S134_sl60_h5, and M134_sl65_h5, which are the representative cases for modes I to mode IV, respectively. In addition, for mode I and II, the evolution of the droplet shape in continuous sequence showing the propagation of the droplet interfaces are plotted in Figure 3. For modes III and IV, since the droplet travels across the periodical boundaries for several times within a short period, we plot the snapshots in Figures 4 and 5 to show the droplet shape at selected moments.

Mode I occurs if the droplet finally gets stuck on the slope after impacting and travelling a short distance. As shown in Figure 3, after impacting at 53.9 ms, the droplet spreads asymmetrically on the surface, developing a bulge at the downstream side and a thin tail at the upstream side. From 53.9 ms to 137.5 ms, inertia drives the droplet to move downslope. The interface undergoes pinning and depinning, from which the capillary wave (defined as a wave travelling along the phase boundary of a fluid) generates and propagates through the droplet. Severe distortion and oscillation are observed and hence, the kinematic energy is dissipated quickly and the sliding velocity slows down. Finally, the gravitational force, which drives the droplet down the slope, can be counterbalanced by pairs of surface tensions at the liquid-air, liquid-solid, and solid-air interfaces. From 137.5 ms, the droplet interfaces are fixed and the droplet is eventually stabilized.

Mode II occurs if the droplet sliding finally reaches a steady state featured by periodically pinning and depinning of its interfaces. Besides, the sliding velocity is constant, which is usually termed as the terminal velocity. The snapshots from 137.5 ms to 188.1 ms in Figure 3 demonstrate that regular pinning and depinning do not cause significant distortion of the droplet, and the interfacial areas for the liquid-air, liquid-solid, and solid-air interfaces remain dynamically constant at the steady state. In terms of its shape during the sliding, the style is called “stick-slip” movement [17]. At this moment, both the kinematic energy and the surface energy are dynamically stable, which indicates that the viscous dissipation is exactly compensated with the gravitational energy input.

Mode III in Figure 4 shows that the droplet undergoes periodic elongation and shortening, in terms of its shape during the sliding. The behavior is the transitional mode from mode II “stick-slip” to mode IV featured by the droplet completely lifting off from the granular surface. Compared to mode II, the gravitational energy input will not be counterbalanced by the viscous dissipation if the droplet is in a stick-slip movement. As a response to facilitate a higher rate of energy dissipation, the droplet undergoes elongation and shortening starting from 229.9 ms. For example, the droplet is shortened at 269.5 ms and elongated at 304.7 ms. Due to the difference in the surface energy, the sliding velocity is severely fluctuating.

Mode IV occurs in the extreme condition of a high CA with a high slope angle. The droplet lifts off the granular surface and jumps down over the slope. When the droplet is flying with no retardation from the granular surface, the droplet is accelerating at a high rate. Each time the droplet falls back on the granular surface (see ① and ③), Figure 5 shows that the higher impacting velocity will drive the droplet to be more and more distorted (see ② and ④). Finally, the droplet is ruptured to halves at 502.7 ms with no terminal velocity being achieved.

4.2. Effect of Impacting Velocity

For the first time, the dependence of the terminal velocity on the impacting velocity is confirmed for a droplet that displays the stick-slip sliding dynamics. The terminal velocity is not unique, as indicated by the temporal evolution of the sliding velocity in Figures 6–9. Here, we select the representative cases S113_sl55, S123_sl30, M113_sl40, and S113_sl35, which all belong to mode II and converge to terminal velocities. It can be seen that most cases display multiple terminal velocities, although some converge to one single value. Figures 6 and 7 indicate that two terminal velocities are observed for the cases S113_sl55 and S123_sl30. Figure 8 shows that for M113_sl40, three terminal velocities are possible. In contrast, there can be one single value irrespective of the impacting velocity, for example, in S113_sl35. However, it should be noted that the terminal velocity does not have a continuous positive-correlated relationship with the impacting velocity. With as many as seven different impacting velocities, only a few discrete values of the terminal velocities can be obtained.

The phenomenon of a droplet sliding in multiple discrete terminal velocities can be understood through an energy approach. The internal energy related to the droplet sliding dynamics consists of the kinematic energy and the surface energy. In the steady state, the constant terminal velocity indicates that the rate of change in the kinematic energy is zero. The rate of change in the surface energy depends on the magnitude of the surface tensions , , and , multiplying the rate of change in the surface areas , , and , for the liquid-vapor, solid-liquid, and solid-vapor interfaces, respectively. The is written as

On average, the segments of the main liquid-vapor droplet interface and the particle surface attached to liquid or vapor phases remain dynamically constant, indicating that is also zero.

Next, the energy input and loss are analyzed. The gravitational energy is the only input energy source, written aswhere and are the sliding velocity and the area of the droplet, respectively. The viscous dissipation rate is

Equations (19) and (20) indicate that the higher sliding velocity at a higher energy level results in a higher gravitational energy input rate, and also higher viscous dissipation rate, which increases with the shear rate of the liquid inside the droplet. The opposite is true for a lower energy level. Despite various energy levels, since and are both zero, the gravitational energy input and the viscous dissipation shall always cancel each other out. Otherwise, the sliding state is unstable, triggering the droplet acceleration or deceleration through its interaction with the granular surface. When the sliding is in such an unstable state, even a small discrepancy of the initial condition may have great influence on the converged energy level. Figure 6 shows that the sliding velocities are very close to each other during the period from 148.5 ms to 336.6 ms, for the cases S113_sl55_h1/6 and S113_sl55_h2/3. However, the small difference causes their divergence at 336.6 ms and S113_sl55_h2/3 turns into a higher energy level by a sudden acceleration. It is worth noting that the droplet may be able to run in an unstable state without significant changes in the velocity, for a relatively longer duration. For instance, Figure 7 shows that the velocity is almost constant from 337.7 ms to 608.3 ms in the case S123_sl30. Finally, the unstable state breaks up and the sliding velocity transits into a higher energy level.

Another important finding is that the terminal velocity may not be positively correlated with the impacting velocity at relatively low slope angles. To better clarify this finding, Figure 10 plots the terminal velocities against various impacting velocities for the two illustrative series of M113 and M134. For example, in the series of M113_sl40 with three energy levels, “h2/3” and “h1” correspond to the low energy level; “h1/6” and “h3” correspond to the middle energy level; and “h2”, “h4,” and “h5” correspond to the high energy level. We see that a higher impacting velocity may not necessarily lead to a higher energy level for a relatively lower slope angle. This is because the impact is more developed in the normal direction to the slope as the spreading and recoiling processes on the complex granular surface. The relatively small downslope component of the impacting velocity does not efficiently drive the droplet to slide down at a higher velocity. Moreover, the recoiling process results in a larger viscous dissipation rate, which retards the downslope movement of the droplet. However, in the case of a higher slope angle, for example, M113_sl60, “h3”, “h4,” and “h5” slide at the higher energy level, while all the other cases with lower impacting velocities correspond to the low energy level. The relatively larger downslope component of the impacting velocity can overcome the energy barrier, so that the droplet tends to slide at a higher energy level. A similar phenomenon can also be observed in the series M134 with a higher CA in the figure. However, the threshold slope angle for the terminal velocity to be positively correlated with the impacting velocity is only 15°, compared with the value of 60° in M113.

(a)

(b)

4.3. Terminal Sliding Velocity

The terminal sliding velocities are summarized in Figure 11 for all the combinations of CA and particle size. At each slope angle, since multiple terminal velocities may be achieved for the simulations conducted with seven different impacting velocities, the values are averaged together with the error bars indicating the upper and lower bounds. The figure shows that the terminal sliding velocity increases with the slope angle since the energy input is higher. For example, in S134, the average terminal velocity increases from 0.0178 m/s on a 5° slope to as large as 0.534 m/s on a 45° slope. The CA is found to have a significant positive effect on the sliding velocities. A larger CA causes less distortion and oscillation every time the depinning of the droplet interface happens at the upstream, leading to less viscous damping and a higher velocity. Take the 0.67 mm particle layer inclined at 40° as an example; the droplet sticks to the granular surface for the CA 104°, irrespective of the impacting velocities. With larger CAs 113°, 123°, 134°, the average terminal velocities are 0.0619 m/s, 0.277 m/s, and 0.492 m/s, respectively. The CA also affects the static behavior of critical slope angle to trigger the downslope sliding. As the depinning of the droplet interface is more difficult, an inclined surface with a smaller CA can retain the droplet statically at a larger slope angle. For example, the minimum slope angles to induce the droplet sliding for the particle size 0.67 mm are 5°, 20°, 35°, and 60° for the CAs 134°, 123°, 113°, and 104°, respectively.

The effect of particle size in terminal velocity is small but not negligible, within the relatively narrow range of coarse sandy sizes. Figure 11 shows that at the same CA, the terminal velocity is generally smaller for a larger particle size, which may be interpreted as the larger geometric roughness and hence result in a higher energy dissipation rate. This is consistent with Wang et al. [20] where increasing the width of micropillars slows down the droplet sliding. Next, the larger geometric roughness for the larger particle creates a larger energy barrier between the higher and the lower energy levels. In other words, it is possible for the droplet to run in significantly different terminal velocities on the particle layer formed by larger particles, as indicated by the larger error bars. Figure 11 also shows the deviation between the upper and lower bounds of the multiple terminal velocities for each configuration of CA and particle size. The deviation is generally larger for the larger particle, except for L104 with only a few samples. Take CA = 113° for instance, the maximum deviations are 0.0816 m/s at the slope angle of 50°, 0.0934 m/s at the slope angle of 40°, and 0.121 m/s at the slope angle of 65°, for the particle sizes 0.67 mm, 0.89 mm, and 1.15 mm, respectively.

4.4. Quantification of the Rotational or Translational Behavior

In this subsection, we propose a new method to quantify the motion style to describe whether a droplet is undergoing rotational or translational movement. For a droplet sliding at a lower velocity, the droplet appears to rotate over the particles in a more rounded shape. With a higher sliding velocity, the droplet is more elongated and slides in a translational motion. Previous studies describing the droplet rotation or translation need to assume a droplet radius so that the velocity can be decomposed into the rotational component and the translational component [18, 49]. However, this approach only applies to a droplet that maintains its rounded shape during sliding.

From the perspective of energy, we propose to use the ratio of the translational kinematic energy in the total kinematic energy to quantitatively describe whether a droplet is rotating or translating. We define it as the effective kinematic ratio EKR, which is written as and are the sliding velocity and the kinematic energy unit area of the droplet as a whole. They are averaged over the region that is occupied by the liquid phase (C > 0.5), based on the Euler description of the fields of phase index C and the velocity u. Such a kinematic ratio from the energy perspective can be applied to a droplet with an arbitrary shape. The value ranges between 0 and 1, from rotation to translation.

Figure 12 plots the effective kinematic ratio EKR against the sliding velocity for CA = 104°, 113°, 123°, and 134°. EKR increases with the sliding velocity in a high rate at low velocities but finally converges to a constant value. Hence, a hyperbolic equation is adopted to empirically fit the relationship. The CA = 104° is not fitted since it does not have enough data points. The particle size has a negligible effect on the EKR. This may due to the fact that the selected particle sizes are all within a narrow range of coarse sandy sizes. Hence, for each CA, the results from all the particle sizes are analyzed together. The table in Figure 12 shows that for the CAs of 113°, 123°, and 134°, the values of a are 0.05978 m/s, 0.04316 m/s, and 0.02227 m/s, respectively, while the values of b are 1.231, 1.218, and 1.187, respectively. Considering the geometric interpretation, the inverse of parameter a means the rate of effective kinematic ratio increasing with the terminal velocity at the point of origin, while the inverse of parameter b means the maximum effective kinematic ratio that the droplet is converging to. The results indicate that the droplet sliding is getting closer to the rotational motion with a lower CA. This is because the depinning of the droplet interface is more difficult for a lower CA, and therefore additional kinematic energy is required to drive the droplet to rotate over the particles. Compared to the multiple terminal velocities that depend on the impacting velocity, we demonstrate that the unique relationship between the EKR and the terminal velocity intrinsically characterizes the particle-level water repellency.

5. Conclusions

The present work investigated the interaction between an inclined hydrophobic granular surface and a millimeter-sized droplet. The droplet dynamics, including the impacting and the subsequent downslope sliding, was simulated with the LB method. Granular surfaces with different geometric configurations and water repellency were constructed. To elucidate the mechanisms controlling the interaction of droplets and hydrophobic granular slopes, a series of two-dimensional numerical droplet sliding tests were carried out on a granular surface formed by a single layer of spherical particles with different sizes and intrinsic CAs. Particularly, the droplet was impacting the granular surface with different velocities by controlling the falling height. The major findings can be drawn as follows:(1)Under different conditions, the droplet dynamics can be described within four modes in the downslope sliding: (I) a droplet sticks to the granular surface, (II) a droplet moves by pinning and depinning of its interface, (III) a droplet undergoes periodic elongation and shortening during sliding, and (IV) a droplet lifts off the granular surface and may be ruptured. A higher CA and a larger slope angle favor the downslope sliding from mode I to mode IV. The higher slope angle means a higher gravitational energy input, and a larger CA causes less energy dissipation, since the easier depinning process triggers less distortion and oscillation.(2)Mode II, a droplet moving steadily by pinning and depinning of its interface (stick-slip) can reach a converged terminal velocity through droplet-particle interaction, which means the viscous dissipation is exactly balanced with gravitational energy input.(3)The dependence of the terminal velocity on the impacting velocity is confirmed, and the relationship is not continuous. We first demonstrate that the droplet has multiple discrete choices of its terminal velocity under different impacting velocities. Moreover, the terminal velocity is positively correlated with the impacting velocity only at relatively high slope angles, since the large downslope component of the impacting velocity can overcome the energy barrier, and the droplet tends to slide at a higher energy level. At low slope angles, the convergence to a higher or lower terminal velocity is unpredictable, since the impact develops in the normal direction to the slope, featured by the spreading and recoiling.(4)The terminal velocity increases with both the slope angle and the CA, which, respectively, indicates a higher gravitational energy input and a lower viscous energy dissipation.(5)The larger geometric roughness induced by a larger particle size results in a higher energy dissipation rate during the droplet-particle interaction, and therefore a smaller terminal velocity. The larger roughness also creates a larger energy barrier, so the droplet may converge to significantly different terminal velocities on granular surfaces formed with larger particles.(6)The effective kinematic ratio, EKR, defined as the translational kinematic energy divided by the total kinematic energy, is proposed as the quantitative indicator to describe a rotational or translational type of movement. Within the “stick-slip” sliding mode, the EKR increases with the terminal velocity at a high rate at small velocities but finally converges to a maximum value. An empirical equation can be adopted to empirically fit the relationship. The fitted results capture the phenomenon that for a larger CA, the EKR increases faster with the terminal velocity and the maximum EKR that can be achieved is also larger. Besides, compared to multiple terminal velocities, which may be affected by the impacting velocity, the unique relationship between the EKR and the terminal velocity is a better indicator to intrinsically characterize the particle-level water repellency.

The current study provides a basis to investigate the factors that affect the efficient removal of a droplet under its gravity on an inclined granular surface, which is fundamental to applications of synthetic water-repellent soils in the geotechnical engineering discipline. The study can also extend to other applications involving self-cleaning rough surfaces, for e.g., water-repellent clothes and antifouling kitchenware.

Data Availability

The data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest

The authors declare that there are no conflicts of interest regarding the publication of this paper.

Acknowledgments

This work was supported by the National Natural Science Foundation of China (Grant no. 61806061).

References

C. Neinhuis and W. Barthlott, “Characterization and distribution of water-repellent, self-cleaning plant surfaces,” Annals of Botany, vol. 79, no. 6, 1997.
View at: Publisher Site | Google Scholar
P. Zhu, T. Kong, X. Tang, and L. Wang, “Well-defined porous membranes for robust omniphobic surfaces via microfluidic emulsion templating,” Nature Communications, vol. 8, p. 15823, 2017.
View at: Publisher Site | Google Scholar
L. F. Debano and J. S. Krammes, “Water repellent soils and their relation to wildfire temperatures,” International Association of Scientific Hydrology. Bulletin, vol. 11, no. 2, pp. 14–19, 1966.
View at: Publisher Site | Google Scholar
G. W. Jex, B. H. Bleakley, D. H. Hubbell, and L. L. Munro, “High humidity-induced increase in water repellency in some sandy soils,” Soil Science Society of America Journal, vol. 49, no. 5, pp. 1177–1182, 1985.
View at: Publisher Site | Google Scholar
J. Hyväluoma, A. Koponen, P. Raiskinmäki, and J. Timonen, “Droplets on inclined rough surfaces,” The European Physical Journal E, vol. 23, no. 3, pp. 289–293, 2007.
View at: Publisher Site | Google Scholar
L. Jing, C. Y. Kwok, Y. F. Leung, and Y. D. Sobral, “Characterization of base roughness for granular chute flows,” Physical Review E, vol. 94, no. 5, p. 52901, 2016.
View at: Publisher Site | Google Scholar
H.-Y. Kim, H. J. Lee, and B. H. Kang, “Sliding of liquid drops down an inclined solid surface,” Journal of Colloid and Interface Science, vol. 247, no. 2, pp. 372–380, 2002.
View at: Publisher Site | Google Scholar
C. Lv, C. Yang, P. Hao, F. He, and Q. Zheng, “Sliding of water droplets on microstructured hydrophobic surfaces,” Langmuir, vol. 26, no. 11, pp. 8704–8708, 2010.
View at: Publisher Site | Google Scholar
N. Savva and S. Kalliadasis, “Droplet motion on inclined heterogeneous substrates,” Journal of Fluid Mechanics, vol. 725, pp. 462–491, 2013.
View at: Publisher Site | Google Scholar
Y. Wang and S. Chen, “Numerical study on droplet sliding across micropillars,” Langmuir, vol. 31, no. 16, pp. 4673–4677, 2015.
View at: Publisher Site | Google Scholar
M. Sbragaglia, L. Biferale, G. Amati et al., “Sliding drops across alternating hydrophobic and hydrophilic stripes,” Physical Review E, vol. 89, no. 1, Article ID 012406, 2014.
View at: Publisher Site | Google Scholar
S. Suzuki, A. Nakajima, K. Tanaka et al., “Sliding behavior of water droplets on line-patterned hydrophobic surfaces,” Applied Surface Science, vol. 254, no. 6, pp. 1797–1805, 2008.
View at: Publisher Site | Google Scholar
S. Varagnolo, D. Ferraro, P. Fantinel et al., “Stick-slip sliding of water drops on chemically heterogeneous surfaces,” Physical Review Letters, vol. 111, no. 6, Article ID 066101, 2013.
View at: Publisher Site | Google Scholar
P. Hao, C. Lv, Z. Yao, and F. He, “Sliding behavior of water droplet on superhydrophobic surface,” EPL (Europhysics Letters), vol. 90, no. 6, Article ID 66003, 2010.
View at: Publisher Site | Google Scholar
T. Podgorski, J.-M. Flesselles, and L. Limat, “Corners, cusps, and pearls in running drops,” Physical Review Letters, vol. 87, no. 3, Article ID 036102, 2001.
View at: Publisher Site | Google Scholar
Q. Wang, Y. Li, Z. Yu, and B. Guo, “Numerical simulation of gravity-driven droplet displacement on an inclined micro-grooved surface,” in Proceedings of the ASME 2016 5th International Conference on Micro/Nanoscale Heat and Mass Transfer, American Society of Mechanical Engineers (ASME), 2016.
View at: Google Scholar
G. Hu and Y. Yao, “Movement of a droplet on a structured substrate: a dissipative particle dynamics simulation,” in APS Division of Fluid Dynamics Meeting Abstracts, 2012.
View at: Google Scholar
J. Zhao, S. Chen, and Y. Liu, “Droplets motion on chemically/topographically heterogeneous surfaces,” Molecular Simulation, vol. 42, no. 17, pp. 1452–1459, 2016.
View at: Publisher Site | Google Scholar
S. Varagnolo, V. Schiocchet, D. Ferraro et al., “Tuning drop motion by chemical patterning of surfaces,” Langmuir, vol. 30, no. 9, pp. 2401–2409, 2014.
View at: Publisher Site | Google Scholar
D. Khojasteh, M. Kazerooni, S. Salarian, and R. Kamali, “Droplet impact on superhydrophobic surfaces: a review of recent developments,” Journal of Industrial and Engineering Chemistry, vol. 42, pp. 1–14, 2016.
View at: Publisher Site | Google Scholar
Y. H. Yeong, J. Burton, E. Loth, and I. S. Bayer, “Drop impact and rebound dynamics on an inclined superhydrophobic surface,” Langmuir, vol. 30, no. 40, pp. 12027–12038, 2014.
View at: Publisher Site | Google Scholar
C. W. Hirt and B. D. Nichols, “Volume of fluid (vof) method for the dynamics of free boundaries,” Journal of Computational Physics, vol. 39, no. 1, pp. 201–225, 1981.
View at: Publisher Site | Google Scholar
W. J. Rider and D. B. Kothe, “Reconstructing volume tracking,” Journal of Computational Physics, vol. 141, no. 2, pp. 112–152, 1998.
View at: Publisher Site | Google Scholar
S. Osher and R. P. Fedkiw, “Level set methods: an overview and some recent results,” Journal of Computational Physics, vol. 169, no. 2, pp. 463–502, 2001.
View at: Publisher Site | Google Scholar
S. Osher and J. A. Sethian, “Fronts propagating with curvature-dependent speed: algorithms based on Hamilton-Jacobi formulations,” Journal of Computational Physics, vol. 79, no. 1, pp. 12–49, 1988.
View at: Publisher Site | Google Scholar
X. Shan and H. Chen, “Lattice Boltzmann model for simulating flows with multiple phases and components,” Physical Review E, vol. 47, no. 3, p. 1815, 1993.
View at: Publisher Site | Google Scholar
X. Shan and H. Chen, “Simulation of nonideal gases and liquid-gas phase transitions by the lattice Boltzmann equation,” Physical Review E, vol. 49, no. 4, p. 2941, 1994.
View at: Publisher Site | Google Scholar
M. R. Swift, W. R. Osborn, and J. M. Yeomans, “Lattice Boltzmann simulation of nonideal fluids,” Physical Review Letters, vol. 75, no. 5, pp. 830–833, 1995.
View at: Publisher Site | Google Scholar
M. R. Swift, E. Orlandini, W. R. Osborn, and J. M. Yeomans, “Lattice Boltzmann simulations of liquid-gas and binary fluid systems,” Physical Review E, vol. 54, no. 5, p. 5041, 1996.
View at: Publisher Site | Google Scholar
T. Inamuro, T. Ogata, S. Tajima, and N. Konishi, “A lattice Boltzmann method for incompressible two-phase flows with large density differences,” Journal of Computational Physics, vol. 198, no. 2, pp. 628–644, 2004.
View at: Publisher Site | Google Scholar
T. Lee and C.-L. Lin, “A stable discretization of the lattice Boltzmann equation for simulation of incompressible two-phase flows at high density ratio,” Journal of Computational Physics, vol. 206, no. 1, pp. 16–47, 2005.
View at: Publisher Site | Google Scholar
T. Lee and L. Liu, “Lattice Boltzmann simulations of micron-scale drop impact on dry surfaces,” Journal of Computational Physics, vol. 229, no. 20, pp. 8045–8063, 2010.
View at: Publisher Site | Google Scholar
L. Hao and P. Cheng, “Lattice Boltzmann simulations of water transport in gas diffusion layer of a polymer electrolyte membrane fuel cell,” Journal of Power Sources, vol. 195, no. 12, pp. 3870–3881, 2010.
View at: Publisher Site | Google Scholar
C. Pan, M. Hilpert, and C. Miller, “Lattice-Boltzmann simulation of two-phase flow in porous media,” Water Resources Research, vol. 40, no. 1, 2004.
View at: Publisher Site | Google Scholar
M. E. A. B. Amara, P. Perre, and S. B. Nasrallah, “Lattice-Boltzmann analysis of capillary rise,” Journal of Porous Media, vol. 19, no. 5, pp. 453–469, 2016.
View at: Google Scholar
F. Jansen and J. Harting, “From bijels to pickering emulsions: a lattice Boltzmann study,” Physical Review E, vol. 83, no. 4, Article ID 046707, 2011.
View at: Publisher Site | Google Scholar
M. E. Cates, J.-C. Desplat, P. Stansell et al., “Physical and computational scaling issues in lattice Boltzmann simulations of binary fluid mixtures,” Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences, vol. 363, no. 1833, pp. 1917–1935, 2005.
View at: Publisher Site | Google Scholar
T. Lee and L. Liu, “Wall boundary conditions in the lattice Boltzmann equation method for nonideal gases,” Physical Review E, vol. 78, no. 1, Article ID 017702, 2008.
View at: Publisher Site | Google Scholar
T. Kevin, “A review of spurious currents in the lattice Boltzmann method for multiphase flows,” Journal of Mechanical Science & Technology, vol. 26, no. 12, pp. 3857–3863, 2012.
View at: Google Scholar
H. Chen, S. Chen, and W. H. Matthaeus, “Recovery of the Navier-Stokes equations using a lattice-gas Boltzmann method,” Physical Review A, vol. 45, no. 8, p. R5339, 1992.
View at: Publisher Site | Google Scholar
X. He, S. Chen, and R. Zhang, “A lattice Boltzmann scheme for incompressible multiphase flow and its application in simulation of Rayleigh-taylor instability,” Journal of Computational Physics, vol. 152, no. 2, pp. 642–663, 1999.
View at: Publisher Site | Google Scholar
Y. Wang, C. Shu, H. B. Huang, and C. J. Teo, “Multiphase lattice Boltzmann flux solver for incompressible multiphase flows with large density ratio,” Journal of Computational Physics, vol. 280, pp. 404–423, 2015.
View at: Publisher Site | Google Scholar
X. Shan, “Pressure tensor calculation in a class of nonideal gas lattice Boltzmann models,” Physical Review E, vol. 77, no. 6, Article ID 066702, 2008.
View at: Publisher Site | Google Scholar
H. Kang, S. D. N. Lourenço, and W. M. Yan, “Lattice Boltzmann simulation of droplet dynamics on granular surfaces with variable wettability,” Physical Review E, vol. 98, no. 1-1, 2018.
View at: Publisher Site | Google Scholar
S. H. Y. Ng and S. D. N. Lourenço, “Conditions to induce water repellency in soils with dimethyldichlorosilane,” Géotechnique, vol. 66, no. 5, pp. 441–444, 2016.
View at: Publisher Site | Google Scholar
Y. Abdelhamid and U. El Shamy, “Pore-scale modeling of surface erosion in a particle bed,” International Journal for Numerical and Analytical Methods in Geomechanics, vol. 38, no. 2, pp. 142–166, 2014.
View at: Publisher Site | Google Scholar
D. F. Boutt, B. K. Cook, B. J. McPherson, and J. Williams, “Direct simulation of fluid-solid mechanics in porous media using the discrete element and lattice-Boltzmann methods,” Journal of Geophysical Research: Solid Earth, vol. 112, no. B10, 2007.
View at: Publisher Site | Google Scholar
K. G. Winkels, J. H. Weijs, A. Eddi, and J. H. Snoeijer, “Initial spreading of low-viscosity drops on partially wetting surfaces,” Physical Review E, vol. 85, no. 5, Article ID 055301, 2012.
View at: Publisher Site | Google Scholar
S. P. Thampi, R. Adhikari, and R. Govindarajan, “Do liquid drops roll or slide on inclined surfaces?” Langmuir, vol. 29, no. 10, pp. 3339–3346, 2013.
View at: Publisher Site | Google Scholar

Copyright

Copyright © 2021 Qing Bao and Hengyi Kang. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

PDF Download Citation

Download other formats

Order printed copies

Views

363

Downloads

464

Citations