Effect of Roughness and Contact Coefficients on the Dynamic Behavior of Droplet Impact onto a Solid Surface

In this paper, a lattice Boltzmann model (LBM) is applied to simulate the dynamic behavior of liquid droplets on a solid substrate with variable wettabilities and different micro-topographical patterns. The relationship among fluid/solid interaction parameter, contact angle and surface wettability is obtained by the linear fitting of the pre-simulation results. In the case of 0.314 < |Gt| < 0.563, 0.149 < |Gt| < 0.314 and 0.066 < |Gt| < 0.149, the contact angle will be θ < 90°, 90° < θ < 150° and 150° < θ < 180°, namely the material is the hydrophilic, hydrophobic and superhydrophobic, respectively. Flowing state of droplets on different surfaces are presented. For the hydrophilic rough surface, the droplet fills up the grooves among roughness bumps, while for the hydrophobic one, liquid drop suspends and sits on the tops of the posts with pockets of air beneath it. Most of all, the effects of roughness on droplet behavior are studied and some comprehensive results have been obtained. When the surface is hydrophilic and super hydrophobic, changing the structure of rough units has less impact on the steady-state flow pattern of droplet. And for hydrophobic surface, the effect of surface roughness on the impact characteristics is significant. Only when the structures of rough units meet that the height of the rough unit is greater than the depth the droplets dipping into grooves, droplets can suspend on the top of the rough units with no contact with the bottom wall of the groove.


Introduction
In various natural and industrial processes, such as spray combustion, pesticide application, internal combustion engines, and printing circuit boards, the impact process of liquid droplet is a basic component.Understanding the accompanying physical phenomena is of significance for many technical applications.(Worthington, 1876) performed the relevant systematical research work firstly.From then on, a large number of theoretical, experimental studies and numerical simulations were performed.
From the point of research interest, the studies have experienced two distinctive stages: 1) the early studies focused on the tracking the drop' shape during the deformation process and the influences of impact velocity, drop size, impact direction relative to the surface and the properties (its density, viscosity, viscoelasticity and so on).The salient details associated with the impact process can be found in the comprehensive reviews reported by (Rein, 1993) and (Yarin, 2005) and the references therein.2) the recent researches laid more emphases upon the effects of surface characteristics, such as chemical properties and topographical characteristics, on the dynamic of drops, it is because that the advent of microfabrication technique leads to the appearance of the substrates whose chemical or topographical properties can be controlled.
As to the latter stage researches, a plenty of experimental investigations have been conducted.For example, different types of micro-patterned substrates have been fabricated as the targeted surfaces and the drops deformation processes have been recorded by means of photographic approach.The typical experiments were reported by (Zhu et al., 2005), (Shinoda et al., 2007(Shinoda et al., , 2009)), (Kannan and Sivakumar, 2014), (McDonald et al., 2008) and so on.The results indicated that the deformation process of impacting liquid drops was dramatically changed by the groove structures.At the same time, the impaction against the hydrophilic or hydrophobic surfaces has been carried out to investigate the wettability characteristics, such as the (Chen et al., 2005), (Bormashenko et al., 2007), (Vaikuntanathan et al, 2014)., and (Lee et al., 2013).In addition, Wang and Huang tried to analyze the relationship for the evolution of dimensionless droplet height and wetting diameter during the initial spreading stage of droplet impingement.However, even with the ingenious experimentations, it is still difficult to prepare so many solid surfaces controlled structures and chemistry to carry out extensive experiments.As a result, a clear basic mechanism and a systemic discipline on the effects of surface wettability and roughness have not been presented.
Therefore, the numerical simulations, which are more powerful and less expensive, have turned out to be very useful for modeling liquid drops on well-dined substrates.At present, several numerical studies were also carried out to investigate this problem.(Dupuis & Yeomans, 2005), (Shinoda et al., 2009), and (Semprebon et al., 2013) employed different mathematic methods to simulate this impaction process.Despite the number of previous studies, it is still not clear the effects of wettability and roughness of surface on liquid drop impacts on solid surface.Especially, the systematic study is remains to be conducted.
In this paper, a new numerical approach, the Shan-Chen Multiphase LBM, is presented to describe the dynamic behavior of liquid droplets on a solid substrate with variable chemical and micro-topographical patterns.The main objective of this work is to obtain some comprehensive results about the effects of wettability and roughness on impacts.The paper is organized as follows.In Section 2, the LBM theory and LBM model for multiphase flow will be briefly reviewed.In Section 3, the model will be validated by representing the Laplace law quite well; the connections between the contact angle and the fluid/solid interaction parameter and between the droplet's surface tension and fluid/fluid interaction parameter will be obtained.In Section 4, numerical simulations will be presented and the effects of wettability and roughness will be discussed.In Section 5, conclusions will be drawn.

Basics of lattice-Boltzmann Method
The idea of the LB approach originates from the kinetic theory of gases.The integro-differential Boltzmann equation describes the relation of the molecular distribution function ( , ) f t x , representing the number of molecules found in the phase space volume element, is given as: (1) where ( , ) is the collision operator representing the changing rate of the particle distribution due to collision.The collision operator is generally simplified by means of the single time relaxation approximation (Bhatnagar et al., 1954 andQian, 1988) Then, the discrete-velocity Boltzmann equation model is written in the following form: ] where ( , ) i f t x is the particle distribution function at position x at time t ; the relaxation time, τ , is calculated by the kinematic viscosity of the fluid . i e is the particle velocity, where the subscript i denotes the ith velocity vector.For the D2Q9 model, the velocity i e is given as (0, 0) i = e for 0 i = ; (cos ,sin ) , and 2 (cos ,sin ) .
Here the lattice speed can be set as  eq i f t x , the corresponding equilibrium particle distribution function, is expressed as: 2 2 eq eq eq eq eq i i where the weights i w is 4/9 for the rest particles ( 0 i = ), 1/9 for 1, 2, 3, 4 i = , and 1/36 for 5, 6, 7,8 i = and 3 s c c = is defined as the lattice speed of sound.

Multiphase (MP) Lattice Boltzmann Model
At the early stage of the development of the LBM, considerable efforts were directed to develop multiphase models.Early examples of lattice gas multiphase model were proposed by (Rothman, 1988) and (Appert and Zaleski, 1990).The lattice Boltzmann implementation of these models began with (Shan & Chen, 1993).There are also so-called ''free energy'' approaches proposed by (Swift et al., 1996) , and ''finite density" models that use the Enskog equation for dense gases have also proposed an approach based on tracking an energy (temperature) component method (Luo, 2000 andHe &Doolen, 2002).Each one has its own advantages and application limits.How to choose it depends on the limit range of the correlation and the working conditions involved in the simulation working conditions.In the present study, the Shan-Chen model (Shan & Chen, 1994), also known as pseudo-potential, is chosen because it is easy to handle fluids with different densities, separate the interface automatically, and implement different wettability conditions.
Long-range intermolecular interaction is introduced into the above system by defining a lattice version of the interaction potential between two neighboring sites x and ′ x ( i t δ ′ = + x x e ).In single-component case, it can be written as: x is a Green's function and ( ) ψ x is the "effective mass" representing the inter-particle interaction.
In this study, the two-phase (gas and liquid) flow model will be adopted and a solid surface will be considered.Thus, the surface tension between one fluid (gas) and the other one (liquid) as well as the interaction force between fluid and solid should be taken into account.According to Eq.5, ( ) F σ x , the fluid/fluid interaction can be expressed in the following form: The interaction with nearest neighbor particles will be sufficient to simulate the basic phenomena of multiphase fluid.Thus, for the D2Q9 model, the Green function ( for other cases, hence, the interaction force between fluid and fluid are simplified as: The value of G controls the strength of the interaction force between fluid components.By varying G , different surface tensions between fluids can be obtained.ψ is the interaction potential, according to (Shan &  Chen, 1994), it must be monotonically increasing and bounded.One scheme of the interaction potential ψ commonly used is Raiskinmäki, 2000, 2002and Hyväluoma, 2004).Some other kinds of schemes are ( ) as proposed by (Martys & Chen, 1996) and (Pan et al., 2004) as well as suggested by (Qian et al., 1988).The first one which has been used extensively is employed in this paper.
In addition, the interaction between the fluid and the solid surface, ( ) t F x , is determined by the presence of the solid in the nearest and next nearest neighbours that surround a fluid node.It has a form as: where s is an indicator function of a solid phase, 1 s = or 0 s = is for fluid or solid, respectively.By adjusting t G , surface wetting characteristic can be realized.
Besides, the gravity, ( ) g F x , which is also considered in the present simulation.It is given by the following form: where g is the gravitational acceleration.
At each lattice, all the forces will be incorporated into governing equation by shifting the eq u in Eq. 4 as: ( , ) ( , ) ( , ) / ( , ) is the total force acting on the fluid, including fluid/fluid interaction F σ , fluid/solid interaction t F as well as the gravity effect g F .
Properties of fluid, such as fluid density ρ , and fluid velocity u , can be obtained from: ( , ) ( , ) / ( , ) As shown by (Shan & Doolen, 1996), with the above definition of the interaction potential, the equation of state for D2Q9 lattice-Boltzmann model can be written as: where the first term on the right-hand side is a kinetic contribution, while the second is a contribution due to inter-particle interaction.As long as the interaction potential is chosen properly, any equation of state can be established.

Surface Tension and Fluid/Fluid Interaction Parameter
The surface tension is obtained by conducting a series of droplet tests according to the Laplace's law which gives the relationship between the final radii of the droplet and pressure difference as / p R σ Δ = .In this research, circular droplets with different radii are generated under different G .Initially, a liquid droplet with the radius of 10 lattices is put in the center of the computational domain (200×200 lattice system), and it is surrounded by vapor.The density ratio is set as 10.And then it will evolve towards a steady state.No external force, such as gravity, is applied.Periodic boundary condition is imposed on each direction of the domain without consideration of wettability.The density ratio of liquid to gas is set to be 10 lattices and 0.85 τ = is employed for all simulations.The relationship between p Δ and 1/ R is presented in Figure2.The data points suggest that the simulation results and the lines are linear regression.

Surface Wettability and fluid/Solid Interaction Parameter
The surface wettability is measured by contact angle.If the contact angle of a surface with water is less than 90° and the surface is tending to wetting, and the surface is then usually classified as hydrophilic.Otherwise, if this angle is greater than 90° and the surface is tending to non-wetting, it is usually classified as hydrophobic.Moreover, if a surface whose contact angle is greater than 150°, it is regarded as super-hydrophobic.Contact angle of water on a surface is determined mainly by the characteristics of the solid surface.
In order to get the values of t G certain pre-simulations must be conducted.A circular liquid droplet with the radius of 10 lattices is placed upon the bottom surface of the computational domain (200×200 lattice system).For both the left-and right-hand sides, periodic boundary condition is employed and for the top and bottom, the bounce-back boundary condition is adopted.The density ratio and the relax time are the same as those used in previous simulations.t G < 0.149, the contact angle between liquid droplet and smooth solid surface will be θ < 90°, 90° < θ < 150° and 150° < θ < 180°, namely the material is the hydrophilic, hydrophobic and superhydrophobic, respectively.These results are concluded in Table 2 and will be used in the following simulations.

Results and Discussions
The impact process of a liquid droplet onto a solid surface will be investigated numerically in this section and the effects of surface roughness and surface wettability will be discussed in details.
A domain with 300×200 lattices is used, within which gas filled.The bottom surface is a rough one, distributed with the rectangle roughness elements.The geometric parameters (height h, width w and separation distance s) will be changed in our simulations.In initial condition, a liquid droplet surrounded by gas is placed upon the roughness elements without any contacting and then collides with the bottom surface at the speed of V0, as presented in Figure5.As to the boundary conditions, the periodic conditions are adopted for the left and right while the solid surface conditions with the bounce-back boundary condition considering the fluid/solid interaction are used for the top and bottom ones as well as the roughness elements' boundaries.Surface roughness, often shortened to roughness, is a component of surface texture.As to the micromachining surface, the micro geometric characteristics, especially the distributions of the nanometer-sized or micrometer-sized peaks and valleys, need to be lain more emphasis because these multiscale protuberances change the actual area of the surface and the contact area between the liquid droplet and the surface.Primarily, the roughness coefficient ( r ), which is defined as the ratio of actual area to apparent area, is applied to describe the changes of the actual area of the surface.The actual area includes the grooves' area and the bumps' area, while the apparent area is one of the ideal surface which is considered to be smooth completely.Thus, it could be calculated as the follow: 2 w s h r w s In addition, our research focuses on the impact process of liquid droplet upon the hydrophobic surfaces, therefore, another coefficient should be introduced to describe the changes of the contact area between the droplet and the hydrophobic surface.(Cassie & Baxter, 1944) studied many kinds of hydrophobic surfaces in the nature and presented a gas cavity model, also called as the C-B model, for hydrophobic surface wettability, as shown in Figure 6.When a liquid droplet gets contacted with a solid surface decorated with the micro-grooves, the cavities among the roughness bumps could not fill with liquid, due to the effects of liquid surface tension, instead, they are full of air.Therefore, the contact area between a liquid droplet and a solid surface consists of two parts, the contact area between the droplet and the small micro-bumps of the solid surface and the area between the droplet and the air in the micro-grooves of the solid surface.In this paper, the ratio of contact area between the liquid droplet and the top of bumps to apparent area is used and defined as the contact coefficient ( s f ), whose calculation formula is given as follow: Figure 6.The gas cavity model proposed by Cassie and Baxter (C-B model)

Analysis of Impact Process
Primarily, the geometric parameters are set as constants (height h = 20, width w =10 and separation distance s =10) and the influences of surface wettabilities on the impact process of droplets are studied.The flow patterns of droplets on the rough surface of hydrophilic and super hydrophobic material are given in the Figure7 and Figure8, respectively.Two surfaces have the same roughness coefficient ( 3 r = ) and contact coefficient ( 0.5 ), the droplet spread out on surface and dip into the grooves between rough units, until reaching the stable state.As shown in Figure7, after 1000 steps, most of the droplet is immersed inside the groove, and a few stay at the top of the rough unit wall forming a thin liquid membrane with a small contact Angle.When the surface material presents super hydrophobic ( | | t G =0.1), the flow form of liquid droplets is more complicated.It can be divided into four processes: spreading, wetting, rebound, fallback and receding.In the initial stage of the impaction, the droplets spread out on the wall along the radius direction under the action of kinetic energy and gravity.At the same time, some of the liquid dips into the groove.It forms the crown shape when t = 200.When immersed liquid touch the bottom wall of the groove, it will rebound because of the strong repulsive force for liquid of the super hydrophobic materials and the surface tension of droplet.At the same time, the whole droplet will also rebound from the top of the rough unit.When t = 1000, the middle of the droplet will lose contact with the solid wall completely.Only on the edge of the both sides will stay small contact.When the droplet rebounds to the highest point, it will fall down under the action of gravity with thinner middle part and thicker edge.When t = 1500, it forms a shape similar to the dumbbell.At the same time, under the action of surface tension, the droplet retracts inward along the radius direction.Finally, the retraction liquid polymerizes in the center part forming a shape close to the spherical.Only a tiny part of the liquid dips inside the groove right now, the vast majority sits on top of the rough unit remaining spherical shape.In this case, the contact angle of the droplet is larger, and the contact area with solid surface is very small.

Effects of Roughness and Contact Coefficients
Moreover, the effects of roughness and contact coefficients on the droplet steady-state shapes are also investigated in this paper.The geometric parameters of roughness units are changed and their values employed in our simulated are listed in Table 3.The simulation results are presented in Figure9 -12.), changing the structure of rough units has less impact on the steady-state flow pattern of droplet, as shown in Figure9.Due to the hydrophilic material, most of the droplet is immersed inside the groove, and a few stay at the top of the rough unit wall forming a thin liquid membrane with a small contact angle.When the wall material is hydrophobic ( | | t G =0.2, intrinsic contact angle is 129.97°), the structure of rough units has a significant impact on the final shape of droplet.Two rules could be observed from Figure8: a) the rough unit is higher, the contact angle under steady state of droplet is larger, when it has the same width and separation distance.As shown in Figure10 (f), (g) and (b), the shape of the droplet changes significantly with the difference of height.When liquid dips into the groove, it presents flat in the case of h = 5 but bulged in the case of h = 20.And b) the smaller the separation distance is, the larger the contact angle is, when the height and width are the same.As shown in Figure10 (d), (b) and (e), in the case of s = 20, most of the droplet dips into the groove, while in the case of s = 5, droplet presents coronary shape sitting on the top of the rough unit.Among these simulation results, only one situation, as shown in Figure 10 (d), could be considered as the compound contact of infiltration in C-B model.
Another hydrophilic surface is used in our simulation ( | | t G =0.15, intrinsic contact angle is 148.82°) and the results are provided in Figure11.In the case of 1, 2, 3, and 4, as shown in Figure11 (a), (b), (c) and (d),only a small amount of liquid immerses inside grooves, but does not contact with the wall at the bottom of the groove.The vast majority of liquid suspends on the top of the rough unit keeping spherical shape with slight difference.These situations could be considered as the compound contact of infiltration in C-B model.On the contrary, in the case of 5, 6, and 7, as shown in Figure11 (e), (f), and (g), the droplet dips into grooves with no contact with the bottom of the trench wall and the compound contact does not be formed.
It could be illustrated from Figure10 and 11 that the structure and separation distance of rough units have a direct impact on the way of infiltration for hydrophobic material.On the basis of studying conscientiously, the following conclusion could be get: when the structure and the size of rough units meet that the height of the rough unit is greater than the depth the droplets dipping into grooves, the droplets can suspend on the top of the rough units with no contact with the bottom wall.The depth of droplets dipping into the groove is related to many factors, such as the Weber number, the Reynolds number and the wettability of mural face material.For the material surface with intrinsic contact angle close to 180°, the arrangement, structure and separation distance of rough units have less impact on changing the final flow state of droplets.But according to C-B model, when the intrinsic contact angle invariant, the contact angle will become lager with the contact coefficient becoming smaller.It could be revealed that, accord to our simulation results, the C-B model does not apply to super hydrophobic material with a intrinsic contact angle nearly 180°.
purpose, c is set as unity under the assumption that lattice spacing x δ and time step t δ are unity as well.The schematic diagram of the discrete velocity set is illustrated in Figure1.

Figure 1 .
Figure 1.Lattice geometry and velocity vectors of the D2Q9 model

Figure 2 .
Figure 2. Relationship between 1/ R and p Δ with different G Several values of | | t G (within the range of 0-0.5) representing the different surface wettabilities are used to simulated.Under the interaction force between fluid and solid phase action, the droplets evolved and spread on the surface with different contact angles.The final state of the process of droplets' evolution are presented in Figure3, where | | G = 0.7 is adopted.It can be found that the contact angle increases with the | | t G decreases.

Figure 5 .
Figure 5.Initial state and boundary condition for droplet impact onto rough surface

Figure 9 .
Figure 9. Flow states of droplet on solid surface with different roughness and contact factor (| | 0.4 t G = )

Figure10.
Figure10.Flow states of droplet on solid surface with different roughness and contact factor ( | | 0.2 t G = )

Figure 12 .
Figure 12.Flow states of droplet on solid surface with different roughness and contact factor ( | | 0.1 t G = )

Table 1 .
The values of σ , b and 2

Table 2 .
The relationship among fluid/solid interaction parameter, contact angle and surface wettability

Table 3 .
Structural parameters for rough elements t G =0.4, intrinsic contact angle is 61.13° parameter