Min Chai,Kun Luo*,Changxiao Shao,Song Chen,Jianren Fan

State Key Laboratory of Clean Energy Utilization,Zhejiang University,Hangzhou 310027,China

1.Introduction

As ubiquitous phenomena encountered in a range of technical applications and natural sciences,drop impacts on to a thin preex is ting liquid film have been studied for more than a century for a deeper knowledge of the gas–liquid–wall interaction.According to Rio booet al.'sinvestigation[1],the evolution phenomena after an impact can be classified into three concepts,i.e.deposition,crown for mationands plash.Gener ally in many situations,a sheet jetting flow almost immediately emanates from the neck region connecting the drop and the preexisting liquid film.The jet then grows into a crown and propagates in a short time.The edge of the crown is unstable and will eventually breaks up into small droplets,which is a typical splash process.

With the rapid development of high speed photograph techniques,a great advance on the topological evolutions mentioned previously has been made.Many researchers[2–4]for example V and erWalet al.systematically studied the affected parameters on after-impact evolutions. A conclusion is generally accepted that impact velocity,drop diameter, film thickness and liquid properties,i.e.density,viscosity and surface tension,are the most important parameters affecting after-impact evolutions.

However,the scale limitations in experimental study and the complex and mutable physics phenomena after drop impacts make experimental study difficult.High- fidelity,detailed numerical simulations,which can represent the internal field information visually and truly,provide a possibility to overcome this difficult task,thereby supplement experimental diagnostics.Xieet al.[5]numerically analyzed the influence of drop velocity and film thickness on single drop impact with moving particle semi-implicit(MPS)method,proving the square-root correctionxc=c(t-t0)0.5originally proposed by Weiss and Yarin[6].Leeet al.[7]investigated the splashing and spreading results using level set(LS)method.This paper also pointed out that 2D simulation shows good agreement with experimental data in incipient stage,which is one of the bases of the present work.Other numerical methods,e.g.VOF and CLSVOF,can be easily found in many other articles[8–11].In consideration of the vast disparity of length and time scales present in the impact simulations,high-resolution is required as computational domains must be large enough to resolve the most energetic scales on the order of the drop diameter,while also providing sufficient resolution to capture the smallest structures which are important to reveal the intrinsic mechanisms.Thus direct numerical simulation(DNS)is preferred.

In summary,almost all these studies reviewed above pay their attention to late time drop dynamics while incipient stage researches are still not thorough,and disputed.Josser andet al.[12]once indicated that jet formation and shape evolution of the emerging jet in neck regions are interesting points remaining to be investigated.Similar point can also be found in Yarin's paper[13].But in spite of the numerous achievements in impact morphology,the intrinsic mechanisms especially for oblique impacts are still on their road up to now.What's more,studies on oblique impacts,e.g.Refs.[14]and[15],are relatively rare compared to those on normal impacts.It's obvious that in the neck regions exits large pressure gradient,which is referred as the reason of jet formation.But Thoravalet al.[16]also proposed a mechanism to explain splashing by destabilizing the liquid sheet base through vortex shedding from the free surface.Spontaneously,here comes an interesting question,which is the starting point of the present work,that is pressure gradient still the reason for jet formation when the impact becomes oblique?So the present paper,using DNS coupled with an accurate level set method(ACLS),concentrates on the early time of a drop impacting onto a thin preexisting liquid film,highlighting the analysis of jet formation and evolution when the film is turned from flat to inclined(from a normal impact to an oblique impact,in other words).Note that,deposition or prompt splash without any observable emitting jets is not taken into consideration in the present work.

2.Numerical Methods

2.1.Governing equations

To describe the immiscible gas–liquid flow,the incompressible Navier–Stokes equations are introduced.The continuity and momentum equations can separately be written as:

where ρ is the density,u is the velocity vector,pis the pressure,μ is the dynamic viscosity and g is the gravity acceleration.

The present paper assumed that the physical properties in each phases are constant,but jump at the gas–liquid interface,which can be represented as[ρ]Γ=ρl-ρg,[μ]Γ=μl-μg.The velocity crossing the interface remains the same,[u]Γ=0,while the pressure at the interface is discontinuous,

where σ is the surface tension coefficient,κ is the curvature and n is the unit normal vector of the interface.

2.2.Accurate conservation level set method

To deal with the complex topological evolution in drop impact,an accurate interface capturing method is needed.There are two main categories,VOF method and LS method.The latter implicitly represents the interface as the zero level set of a continuous function,making it convenient to calculate the interface normal and curvature.This function is usually defined as signed distance function:

regarding?(x,t)＞0 as liquid phase,?(x,t)＜0 as gas phase and ?(x,t)=0 as the interface.Because the gradient of this function is numerically equal to one,|??(x,t)|=1,it's convenient to calculate the interface normal and curvature,which can be written as

However,the signed distance function cannot guarantee mass conservation of each phase even with a re-initialization step.The volume of liquid phase may change significantly as time advance,bringing numerical errors.To overcome this problem,an ACLS method is adopted following the works of Desjardinset al.[17–19]and Shaoet al.[20]using a hyperbolic tangent function:where ε is the thickness of the pro file.Phase interface is located at the is osurface of ψ=0.5.Actually,the hyperbolic tangent function is almost the same as a Heavyside function when ε goes to zero.

The transport equation of the ACLS method is

To modify the deterioration of levelset pro file,are-initializationstep is necessary:

where τψis a pseudo time.

2.3.Solving Navier–Stokes equations

Discretization of the Navier–Stokes equations is based on staggered uniform grid.The spatial discretization of the Navier–Stokes equations is performed using second-order finite central-difference scheme.A second-order semi-implicit iterative procedure[21]for time integration is utilized.

The computations of velocity and pressure fields are decoupled by using the projection method.The jump conditions for the pressure gradient in the Poisson equation as well as for the viscous terms are taken into account by the ghost fluid method(GFM)[17].Consider an interface Γ located atxΓbetween the two grid locationsxiandxi+1,wherexi+1is within the liquid phase.The pressure jump in the pressure Poisson equation is then written as

where ρ*=ρgζ+ρl(1-ζ)and ζ=(xΓ-xi)/Δx.The continuum surface force(CSF)approach[22]is used to deal with the viscous term in Eq.(3).

The full solution procedure is summarized here:

(1)Advance the gas–liquid interface implicitly to determine the corresponding phase properties using ACLS method according to Eqs.(5),(7)and(8).

(2)Advance the velocity field by solving Eq.(2)without pressure gradient.

(3)Solve the Poisson equation making use of GFM.

(4)Correct the velocity field using the pressure gradient with GFM.

2.4.Simulation conditions and initialization

The present work performs 2D simulations of single drop impacting onto a liquid film with an inclined angle θ with the horizontal in ambient air.As shown in Fig.1,the drop is initially assumed as a sphere with a diameterD0.The initial thickness of the film ish.The timetis set to zero when the drop meets the film.The velocity at this moment is impact velocityU.Note thatUand g are always vertical.For convenience,four dimensionless parameters are adopted:

After an impact happens,the jets in the front neck region and back neck region may involute in different ways as shown in Fig.1b.So it's natural to distinguish the radii of the jet base,i.e.the distance from the middle of the base to the impact point,in the front or back neck regions,denoted asxcfandxcbrespectively.The physical properties of liquid and gas are listed in Table 1.

Fig.1.Schematic of a drop impacting onto a liquid film(a)just before and(b)after the impact.

Table 1The physical properties of liquid and gas

The computational domain is set to 8D0×3.6D0as a rectangle,and straggled uniform Cartesian grid is used.Once the grid size is determined,the grid can be successfully generated using linear partition.Then the coordinate variables,e.g.index of the cell and position of the cell in domain,are correspondingly determined and the value of signed distance function for every cell is set.No-slip wall boundary condition is applied to the bottom wall while periodic condition is used for the other boundaries.

3.Results and Discussion

The fundamental normal impact case(i.e.θ=0°)is the same to that reported in Cossaliet al.'s experiment[23]withWe=297,Oh=0.0019 and δ=0.29.Equivalently,this means that the initial drop diameter is 3.82 mm,the impact velocity is 2.389 m·s-1and the film thickness is 1.1078 mm.For the inclined cases,remain these parameters but change the inclined angle from 0°to 15°,30°,45°and 60°.The grid size is set to 1024×512 and the ACLS method is used.As Weiss and Yarin[6]predicted,the jet happens at the initial moments of strong impacts(We≥40)att～10-6to 10-5s,resulting in the serious lack of experimental data.To validate the accuracy of ACLS and grid-independence,we set two normal impact cases with modifications of using standard LS method and 2048×1024 grid size respectively as contrastive examples.The setups of single drop impacts are listed in Table 2.

Table 2The setups of single drop impacts

3.1.Numerical validations

Before any further analysis of numerical results,several validations about the method adopted are necessary.As pointed out in Section 2,LS-based methods have great advantages in handling interface normal and curvature if their biggest problem about mass loss is solved.In Fig.2,the instantaneous liquid mass is presented.The value is dimensionless with the theoretical one which can be easily calculated according to their simulation conditions.All these cases have an overestimation of1%–2%.But only in the case using LS method the overestimation grows as time evolves while the others almost remain in the same.Compare to the case with fine grid,the other cases have bigger numerical errors,which means that the overestimation is mainly due to mesh size rather than numerical method.A convincing conclusion can be drawn that the ACLS method has a significant advantage than the standard LS method in mass conservation and accuracy.What's more,the error from adopted mesh size 1024×512 is acceptable.

Fig.2.Mass conservation of drop impacting onto a liquid film.

The numerical result of normal impact neck jet in Case 1 is qualitatively contrasted to that experimentally investigated by Thoroddsen[24]in Fig.3 as Ref.[11]did.These two results are similar even though their impact conditions are quite different.For a more quantitative validation,the square-root correctionxc=c(t-t0)1/2proposed for normal impact is used.It is rewritten to dimensionless type(xc/D0)2=C2(τ-τ0)for a visualized linear relationship as shown in Fig.4.

For a single impact on a pre-existing liquid layer of thicknessh,base positions can be given by[13]

Fig.3.Qualitatively comparison of neck jet between(a)experimental result in Ref.[24]and(b)numerical result.

The coefficientCis a univariate function of dimensionless liquid film thickness.The prediction of Eq.(11)C=1.231 for Case 1 is quite the same to that calculated from our simulationC=1.232 with a 0.05%disagreement.

In this section,we validate the grid-independence as well as illustrate the superiority of ACLS method than standard LS method in mass conservation,thereby in impact morphology prediction.

3.2.Square-root correction modification

In Fig.4,a linear relationship is observed not only in normal impacts but also in oblique impacts.This means that the correctionxc=c(tt0)1/2is also suitable for oblique impacts though the coefficientsCin the front and back necks are discrepant.To ascertain it,we introduce the theoretically study in Ref.[25],which regards the position of crown base after a normal impact onto a moving liquid film as a circle with moving center point in early stage.The radius of this circle is inferred to be equal to that of normal impact onto a still film if other condition is kept on.Note that the so-called oblique impact in Ref.[25]is different in condition setups from that defined in the present paper.But if we turn the coordinate system to move with the liquid drop with its tangential velocityUsinθ,the oblique impact condition can also be equivalently deemed as a normal impact onto a moving film.So the center point can be derived as

following Ref.[25]when shifting time is relatively small,and thereby simplifications are acceptable.

Meanwhile a good prediction of the radius is also necessary to improve accuracy.Eq.(11)is imperfect and often has an overestimation about 10%–15%attributed partly to the effect of the rather thick liquid film which gives rise to a significant velocity component normal to the wall[6]or of the rather thin liquid film where the liquid is damped by the viscous forces.As derived in Ref.[6],the crown base position is

where τf=2πft,f=U/D0for single drop normal impacts when surface tension is totally dominated by inertia.Re-derive to dimensionless parameters,Eq.(13)becomes

The coefficientCis a univariate function of Reynolds number.This correction bases on an assumptionh≈(υ/f)1/2,which is accurate only in certain circumstances.Both of Eqs.(11)and(14)are incomplete since the coefficientCis a univariate function.To relieve it,this paper proposes a modified equation to take more influence factors into account as

wherep∈[0,1].We recommendp=1 ifh≈(υ/f)1/2andp=0 ifh?(υ/f)1/2.In the present work,pis adopted as zero.Note that Eq.(15)is derived for normal impacts,i.e.,Umeans the component perpendicular to the liquid film.

Fig.4.Evolution of dimensionless radius squared with dimensionless time.

In summary,combine Eqs.(12)and(15)to predict the front and back neck base positions for both normal and oblique impacts as

where non-dimensionalization is based on the velocity magnitudeU.The agreement of Eq.(16)and numerical results,shown in Table 3,is quite good except those of the back base in large inclined angle.In these cases the back neck region phenomenon is more like a wave rather than a jet in early stage,making Eq.(16)unsuitable.What's more,the base position measurement is also inaccurate with large personal errors.

Table 3coefficient C in different cases

3.3.Dynamics of jet formation and evolution

After an impact happens,a neck jet flow arises immediately.Due to the kinetic energy transferring from the drop to the film,the jet gradually evolves into a crown-like structure,and then grows outward and upward.The evolutions of the preset impacts are illustrated in a sequence of snapshots in Fig.5.For a detailed presentation,enlarged views of neck regions are drawn in those snapshots needed.The interfaces are presented by different colors,i.e.,red,black,pink,blue and purple for inclined angle θ =0°,15°,30°,45°,60°.As shown in Fig.5,the ACLS method can conceivably capture complex small scale structures such as cavitation and splashing droplets.Qualitatively,the inclined angle has a significant effect on impact dynamics.The front jet flows are obviously sharp,thin and long while those in the back are eased,thick and short in oblique impacts.When θ increases,the front jets become firstly enhanced and then weakened while the back jets are monolithically weakened.In the present work,θ=15°has the biggest enhancement for the front jets.Radius and height of the crown in the back are decreased when θ increased.In cases where θ≥45°the topology in the back neck region is more like a wave rather than a jet.These phenomena can be explained by kinetic energy transmission.Slightly increasing the inclined angle can concentrate more energy in the front neck region,making the front jet fiercer but the back jet slower.But if the angle becomes too large,the impact is weakened since the normal component of the impact velocity becomes too small.Another interesting point is the angle α between the jet and the film.At first,α in the front is bigger than that in the back with a weak dependence on θ.As time advances,bigger θ makes smaller α in the front neck while α in the back tends to 90°independent on θ.

3.3.1.Normal impacts

Thanks to the advantages of numerical simulation in internal information visualization,which is especially important at incipient stage,a qualitative analysis is feasible.In Fig.6,a max pressure gradient is observed in the impact neck region in Case 1.At early time(t=0.15 ms),the pressure gradient is large because the kinetic energy is firstly concentrated in the neck region.As time advances,a part of the energy is gradually absorbed by the liquid film.Meanwhile,the base of the jet flow becomes thicker.Both of these two reasons make the pressure gradient smaller at later time(t=0.5 ms).The energy transmission can also be observed in Fig.7,which illustrates the maximum pressure in liquid phase as a square root function of time.In normal impacts,the direction of the pressure gradient is firstly the same to the jet direction,and then turns normal to the film,seen in Fig.6.This change leads to a “?！眘hape intermediate topology.

Fig.5.Evolution of the interface for drop impacting onto liquid film((a)–(f)for time=0.15,0.25,0.35,0.5,0.75 and 1 ms.Line color red,black,pink,blue and purple for Cases 1–5.).

Fig.6.Pressure field near the neck region at(a)t=0.15 ms and(b)t=0.5 ms for Case 1.

Fig.7.Maximum pressure as a function of time for Case 1.

Qualitatively,we set five probes,whose positions form a“+”shape,in the front neck region and collect information in five time steps starting att=0.2 ms andt=0.4 ms.If we only consider the effects of pressure gradient and ignore the others,a force equilibrium equation is easily deduced as

Discrete Eq.(17)inXandYdimensions respectively.The equivalent velocity increment owing to pressure gradient per time step(dt=1.0×10-6s),dvip,can be written as

whereAis the cross area,Vis the volume,andaPis the equivalent acceleration owing to pressure gradient.The results are shown in Table 4.At 0.2 ms,the directions of the pressure gradient are–Xand+Y,promoting the neck jet outward and upward.Then as time advances,the acceleration inXdimension becomes small.Whent=0.4 ms,the pressure gradient even impedes the jet from growing outward.But the effects inYdimension are still positive. These observations support the mechanism of pressure gradient.

Eventually,we analyze the influence of vortex.The vorticity field ofQcriterion in liquid phase and the velocity vector at 0.15 ms are shown in Fig.8.Vorticity and direction change of velocity can be found near the neck region,making contribution to the neck jet flow.However,in normal impacts the influence of vorticity is relatively small compared to that of pressure gradient.

細節(jié)描寫就是對具有典型意義的細小情節(jié)或事物的某一細微特征進行細致描寫，如人物的外貌、語言、行為、心理、環(huán)境或場面等都可以根據(jù)作品的需要進行細致描繪。細節(jié)描寫能以小見大，增強作文的表現(xiàn)力，具有較強的藝術(shù)感染力，能突出人物的個性特征，揭示文章的主題思想，起到“秤砣小壓千斤”的作用?！墩Z文課程標(biāo)準(zhǔn)》要求“能根據(jù)需要，運用常見的表達方式寫作，發(fā)展書面語言運用能力”，在小學(xué)語文教材中有許多名篇佳作，其中的寫作方法值得教師深入挖掘，借此指導(dǎo)學(xué)生在習(xí)作中模仿和借鑒。下面筆者將談?wù)勅绾我龑?dǎo)學(xué)生學(xué)習(xí)和運用細節(jié)描寫，增強作文的表現(xiàn)力和感染力，使文章引人入勝。

Fig.8.(a)Vorticity field and(b)velocity vector at 0.15 ms for Case 1.

3.3.2.Oblique impacts

As for oblique impacts in Cases 2–5,the evolutions are roughly similar to those in normal impacts.The pressure fields in liquid phase for Case 3 and Case 4 are shown in Fig.9.It's obvious that the front neck region has extremely large pressure gradients while in the back region the pressure gradient is small and will decrease if θ is increased.This means that the mechanisms of jet formation can be different in the front and back regions in oblique impacts.We purposely present the vorticity fields ofQcriterion and velocity vector at 0.15 ms for θ =30°in Fig.10.The vortex is clearly concentrated near neck regions and forms vortex rings to expel the liquid toward the jet.What's more,the vortex intensities in the front and back regions are similar in spite of the difference in jet flow.A change of velocity reversal can be also obviously found in both neck regions.

In summary,we find that the pressure gradient and vortex are both the direct acting factors of jet formation.They are coexisting and competitive.When θ increases,the dominant mechanisms of jet formation becomes different in the front and back regions,i.e.,the pressure gradient is still the dominant in the front while the effect of vortex becomes relatively larger in the back.As time advances,the jet evolves into a “?！眘hape inter mediate topology and then a crown-like structure due to the direction change of pressure gradient.

Table 4Equivalent velocity increment owing to pressure gradient per time step

Fig.9.Pressure fields in liquid phase for(a)Case 3 and(b)Case 4.

Fig.10.Vorticity field and velocity vector in the(a)front and(b)back neck regions.

Fig.11.PDF distributions at 0.2 ms near the(a)front and(b)back necks.

3.4.Mathematical analysis of competitive mechanisms

In this section we statistically analyze the competitive effects of the pressure gradient and vorticity from the insight of mathematics.The analysis regions in this section are all concentrated on the jet bases in liquid phase.Methods of probability density function(PDF)distribution and correlation coefficient are adopted.In the present work the basic definition of vorticity can be rewritten as

3.4.1.PDF distribution

To describe the interaction between the phase interface and pressure or vorticity,a definition ofλ,i.e.cosine of two vectors,is introduced.However,in a two-dimensional simulation where the velocity is assumed to have no z component,the vorticity vector is always parallel to thezaxis.This inevitable nature makes λ unsuitable for vorticity in the present paper since the vorticity vector will be perpendicular to the interface gradient all the time.So we use the velocity to indirectly present the effect of vorticity.As a result,the interaction between the phase interface and pressure or vorticity can be respectively presented as

The specific PDF distributions of λpand λvare shown in Fig.11 and Fig.12.Owing to space reasons only typical 0°,30°,and 45°cases at 0.2ms and0.6ms,respectively presented as jet formation and evolution stages,are considered.The front and back necks are calculated apart,expressed by symbols filled or not.As Fig.11 shows,a peakcan be obviously found in PDF of λpwhatever the inclined angle is,which means the pressure gradient is the driving force for jet formation.In this stage,λvvalues are concentrated at the left end, revealing that the velocity reversals in neck bases are mainly accordant to interface changes,i.e.the effect of vorticity is not big.However,in the jet evolution stage ast=0.6 ms the inclined angle has an important effect.Neither the interface topology nor the driving forces are alike in the front and back necks.When θ increases,the pressure gradient plays a more significant part in the front neck while its effect in the back neck becomes even smaller.In this stage the λvdistribution tends to be uniform especially in the back necks,hence reveals the existence and influence of vorticity.

3.4.2.Correlation coefficient

Fig.12.PDF distributions at 0.6 ms near the(a)front and(b)back necks.

Fig.13.Correction coefficients of(a)Case 1 and(b)Case 4.

whereiis a scale factor from 0 to 1,corr denotes a correlation calculation algorithm andCIis the correlation coefficient,mathematically representing the combined influence on interface.Correction coefficients for various impact angles at different times are shown in Figs.13 and 14.Wheni=0,i.e.considering only the effect of vorticity,CIis nearly zero due to the 2D simulation nature mentioned above.Conversely wheni=1,only the pressure gradient is considered.A compromise,e.g.seen in Fig.13(a),in the com petition bet ween pressure gradient and vorticity can be identifiably observed wheniincreases since a linear change of scaling factor will not alterCIvalues if there is no compromise exits.The effect of pressure gradient can refer to the final values ofCIwhile that of vorticity can be relatively estimated from the slopes in these figures.The coexisting and competitive effects are rather complex in oblique impacts.As time advances,the pressure gradient has an increasing large influence in the front neck while a stable or smaller one in the back as shown in Fig.13(b).In Fig.14,the effect of impact angles on competition and compromise is presented.As θ increases,the effect of pressure gradient firstly decreases and then increases in the front neck while remains relatively small in the back,consistent with those conclusions drawn before.To estimate the influence of vorticity,slopes are summarized in Table 5 where smaller values mean larger relative effect of vorticity.The effect of vorticity is small in the front neck and will decrease as time advances while it is large in the back.Its accurate and specific fraction value is hard to reach,but there are reasonable grounds to believe that more than a quarter relative fraction can be acceptable in seriously oblique impacts.

Table 5Slopes for various angles and times

Fig.14.Correction coefficients at(a)0.2 ms and(b)0.6 ms.

In this section,PDF distribution and correction calculation are introduced to quantitatively describe the jet formation and evolution dynamics.A compromise in the competition bet ween pressure gradient and vorticity is observed.Their relative effects are complex related to impact angle and time.These mathematical methods will be more suitable for 3D analysis.

4.Conclusions

The present paper,using 2D direct numerical simulation coupled with an accurate level set method,concentrates on the incipient time of sing ledrop impacting on to a th in preexisting liquid film.The adopted ACLS method is validated to be efficient with perfect mass conservation in both normal and oblique impacts.A square-root correction for neck bases is modified in accuracy as well as scope of application.

Processes of jet formation and evolution have been studied to reveal internal dynamics in drop impacts.The pressure gradient and vortex are both the direct acting factors of jet formation while the inclined angle intuitively has a significant effect on them.When θ increases,the dominant mechanisms of jet formation becomes different in the front and back regions,i.e.,the pressure gradient is still the dominant in the front while the vortex effect becomes larger in the back.As time advances,the jet evolves into a “?！眘hape intermediate topology and then a crown-like structure mainly due to the direction change of pressure gradient.To quantitatively analyze the competitive mechanisms,PDF distribution and correction calculation are introduced.A compromise in the competition between pressure gradient and vorticity is observed.

The research shown in this work is not sufficient to reveal the detailed mechanism of drop impacts.But the methodology of ACLS and mathematical analysis presented here shall be a promising tool for further study.

Nomenclature

Cdimensionless square-root coefficient

csquare-root coefficient

OhOhnesorge number

t0shifting time,s

WeWeber number

xcneck base radius(e.g.xcfandxcb),m

δ dimensionless film thickness

τ dimensionless time

τ0dimensionless shifting time

υ kinematic viscosity,m2·s-1

Subscripts

g gas

lliquid

Γ interface

[1]R.Rioboo,C.Bauthier,J.Conti,M.Vou,J.De Coninck,Experimental investigation of splash and crown formation during single drop impact on wetted surfaces,Exp.Fluids35(6)(2003)648–652.

[2]G.E.Cossali,A.Coghe,M.Marengo,Theimpactof a single drop on a wetted solid surface,Exp.Fluids22(6)(1997)463–472.

[3]R.L.V.Wal,G.M.Berger,S.D.Mozes,The splash/non-splash boundary upon a dry surface and thin fluid film,Exp.Fluids40(1)(2006)53–59.

[4]R.L.Vander Wal,G.M.Berger,S.D.Mozes,Droplets splashing upon films of the same fluid of various depths,Exp.Fluids40(1)(2006)33–52.

[5]H.Xie,S.O.Y.Koshizuka,Modelling of a single drop impact onto liquid lmusing particle method,Int.J.Numer.Methods Fluids45(9)(2004)1009–1023.

[6]A.L.Yarin,D.A.Weiss,Impact of drops on solid surfaces:Self-similar capillary waves,and splashing as a new type of kinematic discontinuity,J.Fluid Mech.283(1995)141–173.

[7]S.H.Lee,N.Hur,S.Kang,A numerical analysis of drop impact on liquid film by using a level set method,J.Mech.Sci.Technol.25(10)(2011)2567–2572.

[8]H.Liu,S.C.Wang,M.Z.Xie,M.Jia,J.R.Shi,Y.F.Shi,Numerical simulation of secondary atomization on impacting ofa droplet upon a thin liquid film,J.Combust.Sci.Technol.18(1)(2012)38–43(in Chinese).

[9]G.T.Liang,S.Q.Shen,Y.Yang,CLSVOF simulation for splashing of single drop impact on flat liquid film,J.Therm.Sci.Technol.11(1)(2012)8–12(in Chinese).

[10]M.R.Davidson,Spreading of an inviscid drop impacting on a liquid film,Chem.Eng.Sci.57(2002)3639–3647(PII S0009-2509(02)00266-X17).

[11]G.T.Liang,Y.L.Guo,S.Q.Shen,Analysis of liquid sheet and jet flow mechanism after droplet impinging onto liquid film,Acta Phys.Sin.Chin.Ed.02(2013)415–421(in Chinese).

[12]C.Josserand,S.Zaleski,Droplet splashing on a thin liquid film,Phys.Fluids15(6)(2003)1650.

[13]A.L.Yarin,Drop impact dynamics:Splashing,spreading,receding,bouncing…,Annu.Rev.Fluid Mech.38(2006)159–192.

[14]G.T.Liang,S.Q.Shen,Y.L.Guo,J.X.Chen,H.Yu,Y.Q.Li,Special phenomena of droplet impact on an inclined wetted surface with experimental observation,Acta Phys.Sin.Chin.Ed.62(8)(2013)374–380(in Chinese).

[15]J.F.Dai,X.P.Fan,B.Meng,J.F.Liu,A coupled level-set and volume-of- fluid simulation for splashing of single droplet impact on an inclined liquid film,Acta Phys.Sin.Chin.Ed.64(9)(2015)(in Chinese).

[16]M.J.Thoraval,K.Takehara,T.G.Etoh,S.Popinet,P.Ray,C.Josserand,S.Zaleski,S.T.Thoroddsen,Von Kármán vortex street within an impacting drop,Phys.Rev.Lett.108(26)(2012)264506.

[17]O.Desjardins,V.Moureau,H.Pitsch,An accurate conservative level set/ghost fluid method for simulating turbulent atomization,J.Comput.Phys.227(18)(2008)8395–8416.

[18]E.Olsson,G.Kreiss,A conservative level set method for two phase flow,J.Comput.Phys.210(1)(2005)225–246.

[19]E.Olsson,G.Kreiss,S.Zahedi,A conservative level set method for two phase flow II,J.Comput.Phys.225(1)(2007)785–807.

[20]C.Shao,K.Luo,J.Yang,S.Chen,J.Fan,Accurate level set method for simulations of liquid atomization,Chin.J.Chem.Eng.23(4)(2015)597–604.

[21]C.D.Pierce,P.Moin,Progress-variable approach for large eddy simulation of turbulent combustion,Technical report TF80, flow physics and computation division,Department of Mechanical Engineering Stanford University,2001.

[22]J.U.Brackbill,D.B.Kothe,C.Zemach,A continuum method for modeling surface tension,J.Comput.Phys.100(2)(1992)335–354.

[23]G.E.Cossali,M.Marengo,A.Coghe,S.Zhdanov,The role of time in single drop splash on thin film,Exp.Fluids36(6)(2004)888–900.

[24]S.T.Thoroddsen,The ejecta sheet generated by the impact of a drop,J.Fluid Mech.451(2002).

[25]I.V.Roisman,C.Tropea,Impact of a drop onto a wetted wall:Description of crown formation and propagation,J.Fluid Mech.472(2002).