TRACKING METHODS FOR FREE SURFACE AND SIMULATION OF A LIQUID DROPLET IMPACTING ON A SOLID SURFACE BASED ON SPH*

LI Da-ming, XU Ya-nan, LI Ling-ling, LU Hui-jiao, BAI Ling

State Key Laboratory of Hydraulic Engineering Simulation and Safety, School of Civil Engineering, Tianjin University, Tianjin 300072, China, E-mail: lidaming@tju.edu.cn

TRACKING METHODS FOR FREE SURFACE AND SIMULATION OF A LIQUID DROPLET IMPACTING ON A SOLID SURFACE BASED ON SPH*

LI Da-ming, XU Ya-nan, LI Ling-ling, LU Hui-jiao, BAI Ling

State Key Laboratory of Hydraulic Engineering Simulation and Safety, School of Civil Engineering, Tianjin University, Tianjin 300072, China, E-mail: lidaming@tju.edu.cn

With some popular tracking methods for free surface, simulations of several typical examples are carried out under various flow field conditions. The results show that the Smoothed Particle Hydrodynamics (SPH) method is very suitable in simulating the flow problems with a free surface. A viscous liquid droplet with an initial velocity impacting on a solid surface is simulated based on the SPH method, and the surface tension is considered by searching the free surface particles, the initial impact effect is considered by using the artificial viscosity method, boundary virtual particles and image virtual particles are introduced to deal with the boundary problem, and the boundary defect can be identified quite well. The comparisons of simulated results and experimental photographs show that the SPH method can not only exactly simulate the spreading process and the rebound process of a liquid droplet impacting on a solid surface but also accurately track the free surface particles, simulate the free-surface flow and determine the shape of the free surface due to its particle nature.

Smoothed Particle Hydrodynamics (SPH) method, viscous liquid droplet, impacting on a solid surface, free surface

Introduction

The phenomenon of a liquid droplet impacting on a solid surface can be found extensively in nature, being studied in materials science, chemical engineering, metallurgy, aerospace, manufacturing sector, agriculture and other fields. It is an important issue in production and modern technology.

Wu[1]predicted the size distribution of secondary ejected droplets by the crown splashing of droplets impinging on a solid wall based on experiments, Guo and Dai[2]established a heat and mass transfer model for the thin liquid film on the hot solid surface cooled by the impinging small droplets, with consideration of the effect of the droplet impact, surface tension, thermocapillary, evaporation/condensation, and van der Waals attraction. Li et al.[3]conducted a visual experimental study on droplets being impacted onto a hori-zontal solid surface. Ma et al.[4]carried out elastic / plastic impact simulations of water jet using the SPH and the finite element method. He et al.[5]established a droplet-wall collision model based on Euler-Lagrange two-phase flow model; Zhang et al.[6]simulated the breaking waves impacting on a building in a complex boundary with the Smoothed Particle Hydrodynamics (SPH) method. Various methods were used in studies of impacting problems related with free surface, and the key of a liquid droplet impacting on a solid surface is to track the free surface. Among many movement interface tracking methods, the most popular methods include the Volume Of Fluid (VOF)[7]methods and the Level Set method[8]based on grids, and the SPH method – a meshfree particle method.

The SPH is a new meshfree Lagrangian particle computational technique[9], becoming popular in recent years. It was introduced by Lucy and Monaghan in 1977, was applied in astrophysics at first and then widely used in studying problems of hydrodynamics and other related fields. The SPH method uses a number of discrete particles to replace the whole continuous medium and uses the particle set and the interpolation function to estimate the space function and its derivatives, thus overcomes manydifficulties in the process of fluid simulation based on grids. Therefore, it is very convenient to use the SPH method to simulate fluid flow problems with large deformation and free surface.

The focus of this article is to simulate some typical examples by six popular movement interface tracking methods in various flow field conditions and based on comparisons among the tracking results, to simulate the motion of a viscous liquid droplet with an initial velocity impacting on a solid surface by the SPH method. The surface tension is considered by searching the free surface particles, the initial impact effect is considered by using the artificial viscosity, the boundary virtual particles and the image virtual particles are both introduced to deal with the boundary problem, and the boundary irregularity is taken care of quite well. A comparison between the simulated results of a viscous liquid droplet impacting on a solid surface and experimental photographs is made at last.

1. Basic theory of SPH method

1.1Basic principle of SPH

The SPH method is a meshfree particle method, using a number of discrete particles to replace the problem domain, approximating “the kernel function estimation” by an integral kernel named “the kernel function”, and approximating the field function by an integral form, which is called the kernel approximation method as shown in Eq.(1), and on these bases, the particles are applied to approximate the kernel approximation equation further, the method of superimposing the corresponding value of the adjacent particles in subdomain to replace the integral forms of the function and its derivative is called the particle approximation method as shown in Eq.(2), in this way, the particle approximation equations of the hydrodynamics equations are obtained as

wherefdenotes the field function,Ωthe calculating domain,xthe coordinate vector,hthe smoothing length andWthe kernel function.

whereρis the density, andmthe mass.

1.2Fluid governing equations of SPH method

The basic equations of the SPH method can be obtained by solving the Lagrange two-dimensional N-S equations, as follows:

Mass conservation equation

1.3Searching free surface particles and calculating surface tension

It is difficult to search the free surface particles in the numerical calculation. The Koshizuka method is always applied to identify the free surface particles in the SPH method: ifρi＜βρ0,iis considered as a free surface particle, whereβis a parameter of less than 1, whose value is always taken to be about 0.99. However, a liquid droplet impacting on a solid surface is a strong impacting problem, it is better to use the continuous density method to calculate the density of the particles rather than the Koshizuka method.

For a liquid droplet impacting on a solid surface without splash, an improved searching method for free surface particles is adopted in this article. The main idea of the method is as follows: a part of the particles inside are removed by the Koshizuka method firstly, then the calculating domain is divided intoKparts inxdirection with a certain spacing, according to the initial spacing of particles, in every subdomain, all particles are sorted from small to big according to the value ofy, and the value Δyof every two sorted particles is calculated, then a comparison is made between Δyand a certain limit, defined by trial computations, which usually takes a value 1.5 times the initial spacing of particles, if the value of Δyis bigger than that value, then this place is considered as a discontinuity surface and the two particles will be marked in the surface particle document, else if there is no discontinuity surface in the subdomain, the particle with the biggest value ofywill be marked in the surface particle document, all the surface particles in subdomains marked in the surface particle document and sorted from left to right inxdirection, form the free surface of the whole fluid.

The surface tension in the free surface is calculated by the surface tension equation

wherep0denotes the atmosphere pressure,κthe surface tension coefficient, andζthe curvature.

The curvature of the surface particles is calculated by the cubic spline curve fitting when the surface particles are determined, then the pressure of the surface particles can be calculated according to Eq.(6). Because the surface tension takes effect near the liquid level and the dividing spacing of particles is usually very small, it can be assumed that the pressure caused by the surface tension should not only be added to the free surface particles but also to all the particles in the influential domain of the surface particles.

1.4Artificial viscosity

For the problem of a liquid droplet impacting on a solid surface, the artificial viscosity[10]is adopted in the simulation in order to avoid nonphysical oscillations in calculation. The Monaghan artificial viscosity was used extensively in the SPH method, which not only provides for the dissipation that the shock wave surface requires but also prevents the nonphysical penetration when the particles move towards each other. The form of the Monaghan artificial viscosity is as followα∏andβ∏are the standard constants, whose values are always taken to be about 1.0,φ=0.1hij, which is used to solve the numerical divergence when the particles move towards each other,candvare the sound velocity and the velocity vector of particles,α∏is related to the bulk viscosity andβ∏is related to von Neumann-Richtmyer artificial viscosity, which are used to avoid mutual penetration of particles.

1.5Boundary treatments

Two kinds of boundary conditions are considered in this article, the free surface boundary and the rigid boundary. The SPH method can treat the free surface boundary condition naturally, but when it comes to a rigid boundary condition, a boundary truncation will take place if the particles are to be integrated on the boundary or near the boundary, which may affect the gradient calculation of variables near the boundary, where a careful treatment is required. The detailed treatment scheme can be found in Ref.[9].

Fig.1 Initial data in constant flow field

Two kinds of virtual particles are adopted in the rigid boundary condition. One kind of particles are distributed on the boundary directly, and are called the boundary virtual particles, to be used to generate a strong repulsive force on the inner real particles near the boundary, and in this way, to avoid the nonphysical penetration of real particles to the boundary. The expression of the repulsive force is

where the parametersn1andn2usually take the values of 12 and 4, respectively, the value ofDis fixed according to practical problems, which usually

Fig.2 Tracking results of free surface in constant flow field

takes a standard value of the square of the utmost velocity,r0is an approximate value of the initial particle spacing.

The second kind of virtual particles are called the image virtual particles, and are formed in the following manner: to a real particle with a distance to the boundary withinκhi, a corresponding virtual particle is symmetrically distributed outside the boundary, with the same density, pressure and smoothing length, an equal but opposite velocity as the corresponding real particle. Thus, the boundary defects in the SPH method are remedied and the boundary condition fixed in this way is in conformity with the generalized boundary conditions proposed in Ref.[9]. Because the boundary penetration of real particles can not be prevented completely, the first kind of virtual particles must also be used.

1.6Artificial compression ratio

The liquid in the simulation is incompressible, so the real state equations would set limits on the length of the time step. In order to calculate the pressure in the momentum equation effectively, slightly compressible fluid is used to simulate incompressible fluid in the SPH method. The artificial compressible method proposed by Monaghan (1994) is adopted in this article, and the form is as follow

whereγis a constant,γ=7,ρ0denotes the reference density,ρ0=1000kg/m3,patmis the atmosphere pressure,csis the sound velocity, with a value usually 10 times of the utmost flow velocity, but it can also be calculated ascs=Vmax/(Δρ/ρ)1/2, whereVmaxis the predicted utmost flow velocity.

2. Comparative analysis of tracking methods

Simulations of several typical examples based on the VOF methods (Hirt and Nichols method[7], FLAIR method[11], Youngs method[12], FCT method[13]), the level set method and the SPH method are carried out in various flow field conditions.

2.1Constant flow field

For the constant flow field, we assume that (u,v)=(0.5,?0.5) as shown in Fig.1(a), the initial state is a square with sides of 0.2 m as shown in Fig.1(b), and the whole calculating domain is [0,1]× [0,1].

With the method of VOF, the value of the functioncinside the square area is 1.0, and is 0.0 elsewhere, with grids 100×100, spatial step Δx=0.01, time step Δt=0.1Δx=0.001.

With the level set method, the contour?=0 of the level set function is the boundary of the square, and?＜0 inside the square,?＞0 outside the square. Because the finite volume method is used to discrete equations, the level set functions need to be placed at the center of cells. In order to compare with the VOF method, the computing field is extended to[?0.005,1.005]×[?0.005,1.005], with grids 101×101, spatial step Δx=0.01, and time step Δt=0.1.

With the SPH method, the rectangular arranging node method with spacing Δx=Δy=0.01 is used in this simulation as shown in Fig.1(c), the particles on the boundary of the square are marked as shown in Fig.1(d), all the particles are assigned a valuec,c=0 when the particles are on the boundary, andc=1 elsewhere; the free surface is traced by tracking the marked particles.

The simulated results att=0.5s andt=1s by the methods mentioned above in the constant flow field is shown in Fig.2. It can be seen that the exact translation location of the square can be obtained by all methods, but there is some shape deformation in varying degrees. There is a gradient change on top and bottom sides of the square as obtained by the Hirt and Nichols VOF method, the square obtained by the FLAIR VOF method is in a high degree of smoothness and similar to a circle, the squares obtained by the Youngs VOF method and the level set method are both smooth at corners, almost without deformation, the squares obtained by the FCT VOF method and the SPH method are the best.

Fig.3 Initial data in rotational flow field

u(x,y)=?π(y?0.5),v(x,y)=π(x?0.5), as shown in Fig.3(a), the initial state is a circle with a notch as shown in Fig.3(b), its radius is 0.4 m, the centre is at (0.5,0.5), the width of the notch is 0.15 m, the length is 0.393 m, which forms the famous Zalesak Disk Problem[14]. The whole calculating domain is [0,1]× [0,1].

With the VOF method and the Level Set method, the grids, spatial steps and time steps are same as those in the constant flow field.

With Δx=0.01 and by encircling the centre of the circle by 1o, 40 particles are arranged by the circle arranging node method, the rectangular arranging node method with spacing Δx=Δy=0.01 is used in the notch as shown in Fig.3(c). The particles on the disk boundary are marked as shown in Fig.3(d), with all the particles being assigned a valuec,c=0 when the particles are on the boundary, andc=1 elsewhere. The free surface is identified by tracking the marked particles.

The sim ulated resu lts att=0.5sandt=1 s by themethodsmentionedaboveintherotationalflow

2.2Rotational flow field

Fig.4 Tracking results of free surface in rotational flow field

For a rotational flow field, we assume that field is shown in Fig.4. The rotational flow field is different from the constant flow field, and there is no preferable direction in its velocity field, and the Zalesak Disk Problem concerns a circle with a notch, so it is difficult to keep a high resolution at the sharp corner of the notch. The calculation results show that: the tracking disks obtained by the Hirt and Nichols VOF methods are quite coarse, with serrated edges in both circle and notch, the circles obtained by the FCT VOF, FLAIR VOF and Youngs VOF methods are smooth, but with deformation in the notch, especially, in corners, the result obtained by the Youngs method is better, with a little smoothness in the notch. With no deformation and a high resolution at corners, the Zalesak Disk obtained by the SPH method is the best.

2.3Shear flow field

For a shear flow field, it is assumed that

as shown in Fig.5(a), the velocity field is repeated every second, as follows:

The initial state is a circle as shown in Fig.5(b), with a radius of 0.2 m, and the centre at (0.5,0.3). The whole calculating domain is [0,1]×[0,1].

Fig.5 Initial data in shear flow field

With the VOF and the Level Set methods, the grids, spatial steps and time steps are the same as those in the constant flow field and the rotational flow field.

With Δx=0.01 and by encircling the centre of the circle by 1o, 20 particles are arranged by the circle arranging node method in the circle as shown in Fig.5(c). The particles on the circle boundary are marked as shown in Fig.5(d), with all particles being assigned a valuec,c=0 when the particles are on the boundary, andc=1 elsewhere. The free surface is identified by tracking the marked particles.

There is no deformation in the movement interface in the constant flow field and the rotational flow field, but a real fluid flow field includes tension, shear, collision, broken and so on. So the simulation in a shear flow field is more reflective of the real fluid flow and is better in testing the performance of the tracking methods. The simulated results att=1sandt=2s in the shear flow field are shown in Fig.6. It can be seen that the simulated results by the VOF methods are far from perfect, the circles att=2s do not return to the initial state completely, and there are bubbles in the sheared figures, the simulated results obtained by the Level Set method are better, and those obtained by SPH method are the best.

2.4Analysis of comparative results

(1) The movement interface built by the Hirt and Nichols VOF method is only of the first order accuracy, so parts of the movement interface see a stepped distribution, the normal of the movement interface in the grids has to be calculated in the Youngs VOF method, which makes the movement interface more smooth.

(2) A form of the first order accuracy is adopted by the VOF methods for the force conservation of the mass locally, so the numerical dispersion in results can not be controlled, which is responsible for the presence of bubbles.

(3) The outline of the interface obtained by the Level Set method is more clear than those obtained by the VOF methods, when the movement of the free interface is complex, the VOF methods produce results not only dispersive but also dissipative.

(4) Because in the Level Set method, a high precision discrete form is adopted in solving equations, it can very well simulate the fluid motion with large deformation. But the VOF methods are better than the Level Set method for simple fluid motions, such as translation and rotation.

(5) Because of the particle nature of the SPH method, it can better simulate the free surface motion, track the free surface position and describe the free surface shape than other methods.αΠ=1.0 andβΠ=2.0, the viscosity coefficient is 1.0×10–3Pa·s, the coefficient of surface tension is 7.3×10–4N/m,Vmaxis 2.01 m/s, the fluid compression is 0.01, the atmosphere pressure isPatm= 0.0Pa, the initial density is calculated by formula (9), the time step is 1.0×10–5s, which is determined by trials, in searching the free surface particles, the density coefficientβis 0.99, the value to be compared with Δyis deter- mined by trials and is 1.5 times the initial spacing of particles.

3. Simulation of a liquid droplet impacting on a solid surface

3.1Computing model and parameters

A liquid droplet impacting on a solid surface is simulated, with the length of the solid surface of 0.06 m, the initial height of the droplet center of 0.045 m, the initial velocity of 2.0 m/s, and the radius of the droplet of 0.005 m, The rectangular arranging node method is adopted in the simulation, with 556 real particles in the droplet and 239 virtual particles in the solid surface. The complete pairing searching method is used to search adjacent particles. The smoothing length is 1.5 times of the initial spacing of particles, the artificial viscosity coefficients are

Fig.6 Tracking results of free surface in shear flow field

3.2Simulated results and analyses

Figures 7(a), 7(b) and 7(c) are the particle distributions, the outlines of the free surface and the experimental photos of a viscous liquid droplet impacting on a solid surface at different instants, which show that the motion of a liquid droplet impacting on a solid can be mainly divided into two stages: the spreading process and the rebound process, and the inertial force and the surface tension play an important role in the two processes, respectively. In the spreading process, the droplet spreads in the horizontal radial direction mainly du e to the inertial force after the collision, as showninfiguresfromt=2.1×10?2s tot=2.4× 10?2s. Beforet=2.4×10?2s , the base radius of the droplet increases gradually, with a shape high in centre but low on edges, whent=2.4×10?2s, the base radius of the droplet reaches the largest value with a shape slightly low in centre but slightly high on edge. The rebound process takes place after the horizontal radius of the droplet reaches the largest value, and the base radius of the droplet decreases but the vertical height increases mainly due to the surface tension as shown in figures fromt=2.6×10?2s tot= 3.3×10?2s. When the height of droplet reaches a certain value, the droplet will repeat the process ofspreading and rebound till an equilibrium is reached. The results agree with the experimental results proposed by Kim and Chun[15]as shown in Fig.7(c). The free surface outlines of the droplet at the corresponding time can be obtained from the calculated results, as the objective of using SPH method.

Fig.7 The comparison between simulated results and experimental results

4. Conclusions

(1) The SPH method can overcome many shortcomings of the VOF methods and the Level Set method under various flow field conditions, can track the free surface effectively and is applicable in simulating the flow problems with a strong nonlinear free surface.

(2) The motion of a liquid droplet impacting on a solid surface can be mainly divided into two stages: the spreading process and the rebound process, the inertial force plays a main role in the motion state of the droplet in the spreading process, and the surface tension plays a main role in the motion state of the droplet in the rebound process.

(3) Because of its particle nature, the SPH method can accurately simulate the free surface motion, track the free surface position and describe the free surface shape. The numerical simulated results of a liquid droplet impacting on a solid surface are basically in agreement with the experimental results.

[1] WU Zi-niu. Prediction of the size distribution of secondary ejected droplets by crown splashing of droplets impinging on a solid wall[J].Probabilistic Engineering Mechanics,2003, 18(3): 241-249 .

[2] GUO Jia-hong, DAI Shi-qiang. Research on stability of liquid film on hot solid surface impinged by small droplets[J].Journal of Hydrodynamics, Ser. B,2007, 19(3): 264-271.

[3] LI Wei-zhong, ZHU Wei-ying and QUAN Sheng-lin et al. Visual experimental study on droplet impacted onto horizontal solid surface[J].Journal of Thermal Science and Technology,2008, 7(2): 155-160(in Chinese).

[4] MA Li, TAO Wei-ming and GUO Yi-mu et al. Elastic/ plastic impact simulation of water jet using smoothed particle hydrodynamic and finite element method[J].Journal of Zhejiang University (Engineering Science),2008, 42(2): 259-263(in Chinese).

[5] HE Zheng, GAO Ye and GU Xuan et al. Investigating a droplet-wall collision model[J].Journal of HarbinEngineering University,2009, 30(3): 267-270(in Chinese).

[6] ZHANG Xue-ying, PAN Zhong-jian and MA Xiao-fei. Numerical simulation for breaking waves impacting on a building in complex boundary with the SPH method[J].Chinese Journal of Hydrodynamics,2010, 25(3): 332-336(in Chinese).

[7] HIRT C. W., NICHOLS D. B. Volume of fluid method for the dynamics of free boundaries[J].Journal of Computational Physics,1981, 39(1): 201-225

[8] OSHER S., SETHIAN J. A. Fronts propagating with curvature-dependent speed: Algorithms based on Hamilton-Jacobi formulations[J].Journal of Computational Physics,1988, 79(l): 12-49

[9] GONG Kai, LIU Hua and WANG Ben-long. Water entry of a wedge based on SPH model with an improved boundary treatment[J].Journal of Hydrodynamics,2009, 21(6): 750-757.

[10] ZHENG Kun, SUN Zhao-chen and SUN Jia-wen. Numerical simulations of water wave dynamics based on SPH methods[J].Journal of Hydrodynamics,2009, 21(6): 843-850.

[11] LI Da-ming, CHEN Hai-zhou and ZHANG Jian-wei. The comparison of two different load methods of external force in simulating liquid sloshing by SPH[J].Engineering Mechanics,2010, 27(2): 241-249(in Chinese).

[12] DU Xiao-tao, WU Wei and GONG Kai et al. Two dimensional SPH simulation of water waves generated by underwater landslide[J].Journal of Hydrodynamics, Ser. A,2006, 21(5): 579-586(in Chinese).

[13] HU Lei, GUO Jia-hong and WANG Xiao-yong. Numerical simulation of two droplets impact on the liquid film at the same time using lattice Boltzmann method[J].Chinese Journal of Hydrodynamics,2011, 26(1): 11-18(in Chinese).

[14] VINCENT S., CALTAGIRONE J. P. Efficient solving method for unsteady incompressible interfacial flow problems[J].International Journal for Numerical Methods In Fluids,1999, 30(6): 795-811.

[15] KIM H.-Y., CHUN J.-H. The recoiling of liquid droplets upon collision with solid surfaces[J].Physics of Fluids,2001, 13(3): 643-658.

December 16, 2010, Revised April 4, 2011)

* Project supported by the National Natural Science Foundation of China (Grant No. 51079095), the Science Fund for Creative Research Groups of the National Natural Science Foundation of China (Grant No. 51021004).

Biography: LI Da-ming (1957-), Male, Ph. D., Professor

2011,23(4):447-456

10.1016/S1001-6058(10)60135-7