Application of particle splitting method for both hydrostatic and hydrodynamic cases in SPH

W.T.Liu·P.N.Sun·F.R.Ming·A.M.Zhang

Abstract Smoothed particle hydrodynamics(SPH)method with numerical diffusive terms shows satisfactory stability and accuracy in some violent fluid–solid interaction problems.However,in most simulations,uniform particle distributions are used and the multi-resolution,which can obviously improve the local accuracy and the overall computational efficiency,has seldom been applied.In this paper,a dynamic particle splitting method is applied and it allows for the simulation of both hydrostatic and hydrodynamic problems.The splitting algorithm is that,when a coarse(mother)particle enters the splitting region,it will be split into four daughter particles,which inherit the physical parameters of the mother particle.In the particle splitting process,conservations of mass,momentum and energy are ensured.Based on the error analysis,the splitting technique is designed to allow the optimal accuracy at the interface between the coarse and refined particles and this is particularly important in the simulation of hydrostatic cases.Finally,the scheme is validated by five basic cases,which demonstrate that the present SPH model with a particle splitting technique is of high accuracy and efficiency and is capable for the simulation of a wide range of hydrodynamic problems.

Keywords Smoothed particle hydrodynamics·Particle splitting ·Particle refinement·Solid wall boundary ·Fluid–structure interaction

1 Introduction

Recently,smoothed particle hydrodynamics(SPH)method has been widely applied in the simulation of some violent fluid-structure interaction(FSI)phenomena in the field of the naval architecture and ocean engineering[1–3].Thanks to its Lagrangian characteristic,SPH is very robust in the simulation of rolling,breaking,splashing phenomena of the free surface,e.g.the rolling bow wave[1],sloshing[4],nonlinear wave-body interactions[5],viscoelastic fluid flow[6,7],etc.Meanwhile,SPH has also been widely applied in the simulation of some multiphase flows,e.g.the dynamic behavior of rising bubbles[8,9],underwater explosions[10,11],etc.Besides,SPH has also been introduced to solve some FSI problems by coupling with other numerical solvers,e.g. finite element method(FEM)[12].

On one hand,the research on SPH method is conducted focusing on improving the approximation accuracy and the numerical stability[13,14];on the other hand,researchers are trying to improve the computational efficiency[15,16].Indeed,these aspects,accuracy,stability and efficiency,are correlated to each other.For example,in weakly compressible SPH(WCSPH),due to the assumption of the fluid to be weakly compressible,the time step is restricted to the stability reasons,which makes a low efficiency of WCSPH in the practical applications if no parallel computing methods are used;In incompressible SPH(ISPH),the time step is free of the stability conditions,therefore the amplitude of the time step may be larger than that of WCSPH[17].However,as pointed out in Ref.[18],the increase of the time step in ISPH reduces the accuracy of the numerical results,so the time step in ISPH is actually constrained by accuracy reasons,even not stability reasons.Besides,due to the solving of Poisson equation and the capturing of the free surface,the time consuming of each calculating step in ISPH is much longer than that in WCSPH.As a result,we get the conclusion that the efficiencies are comparative between ISPH and WCSPH,and they are both actually constrained by the either stability reasons or accuracy reasons,which makes the efficiency of SPH method not-so-satisfying when compared to some traditional mesh based methods,e.g.boundary element method(BEM),FEM,etc.

Along with the development of the SPH theory and applications,more and more published papers are also aimed to address the efficiency problems.One way is to use advanced computational technologies,e.g.adopting parallel computing method based on high performance computers;another way is to reduce the amount of computation.For example,similar to the adaptive mesh refinement(AMR)method in some mesh based methods,one can use particle splitting method in SPH to refine the particle resolution in some local regions,while in the rest part of the flow,coarse particles are still distributed.The latter method decreases the total amount of computation,besides that,it also increases the local accuracy in the local regions.Due to this reason,the particle splitting has become very attractive recently.

In the early stage,SPH with the variable smoothing length is often applied in a way that the particle resolution is smoothly decreased from the region of interest up to the borders of the computational domain,see e.g.Liu et al.[19,20].In a similar way,Koukouvinis et al.[21]refined the particle resolutions in some areas of the flow at the initial stage,but during the simulation,no more particles are refined.The above methods are easy to be conducted since the total particle number is conserved.However,these methods cannot dynamically refine the particle resolution since in some cases we only need to refine the particle resolution in a certain area when the flow evolves to a certain stage.Moreover,the interface between the coarse particles and the refined particles(named as splitting interface hereinafter)are usually difficult to treat.Numerical stability problems often take place,for example,the penetration and mixture observed at the splitting interface,and,therefore,the particle resolution needs to be varied very smoothly from the flow center to the border,which makes the arrangement of the particles very difficult especially for some simulations characterized by complex boundaries.

In Ref.[22],a dynamic particle splitting was proposed and it was very attractive since only the mother particle which enters the region of interest will be split.Two splitting patterns were proposed in that work,characterized by a hexagon and triangle,and they allow for satisfying accuracy at the splitting interface.However,both the two kinds of patterns cannot refine the ghost particles very well in the solid wall boundary(e.g. fixed ghost particle)since the ghost particles are usually distributed in a squared lattice and after the refinement,the distribution of the daughter particles will be quite disordered which worsens the accuracy there.Omidvar et al.[23]simplified Feldman’s method with a squared particle splitting pattern in two dimensions and a cubic pattern in three dimensions,however,in their work,particle clumping phenomenon caused by tension instability exists,which was addressed by applying a punishment force.López et al.[24]gave an improved method for the splitting error analysis in their work,where the splitting error analysis based on?f(r)(f(r)is a generic variable at the particle location denoted by vector r)is applied.However,since the smoothing lengths of the daughter particles are constrained in order to avoid particle clumping,an partially overlapped distribution of daughter particles is adopted in Ref.[24]to decrease the splitting error.In that way,the distribution of the daughter particle becomes quite disordered.However,as pointed in Ref.[25],the rate of convergence of the SPH simulation is tightly related to the spatial disorder of the particle distribution.For the simulation without free surface,Vacondio et al.[26]applied a position adjustment method to avoid the particle clumping,however,that method is difficult to be applied in free surface flows.As an alternative,Omidvar et al.[27]applied a punishment force to avoid particle clumping.Recently,based on the Riemann-SPH,Barcarolo[18] successfully applied López’s scheme to simulate both hydrostatic and hydrodynamic problems.In Ref.[28],an adaptive particle refinement(APR)with overlapping particles are also proposed.Recently APR is also combined with the well-known δ-SPH in Ref.[16].One limitation of APR with overlapping particles is that it is difficult to ensure strict conservation of the momentum in the flow since it uses a buffer zone to exchange the physical information between different particle resolution levels,see more in Ref.[28].In addition to that,the overlapping of particles of different resolution levels enlarges the number of neighbor particles which may reduce the efficiency from the programming point of view.Therefore,in the present work,we seek to use only particle splitting and test its performance in both hydrostatic and hydrodynamic problems.

It can be seen from the literature reviewed above that,the researches regarding the particle splitting in SPH have been widely conducted,however there still exists the following problems:(1)The weakly compressible SPH with numerical diffusive terms(also known as δ-SPH in the literature),which has been well applied in many practical problems,has still not been well combined with the particle splitting method;(2)Only Riemann-SPH has been proven to give good results in the hydrostatic case which is very significant in some practical ocean engineering problems,for example in the development of the numerical water wave tank;however the performance of this particle splitting technique in other SPH variants are not tested;(3)As the smoothing length of the daughter particle is increased,the particle clumping is usually excited,which constrains the improving of accuracy on the splitting interface since it will be shown in the following part of this paper that nearly identical smoothing lengths on the interface supplies the highest accuracy;(4)Apart from the repulsive force boundary[29]and mirror-ing ghost boundary[30],a robust solid wall boundary namely fixed ghost particle boundary(see Ref.[31])deserves to be tested in the violent FSI problems when a particle splitting technique is implemented since we find that the last boundary treatment gives the most satisfying accuracy.

In the present work,based on the standard SPH method with numerical diffusive terms,a squared splitting pattern(similar to Ref.[24])is applied,while the splitting error is analyzed in a broader varying range of smoothing length,and,therefore,the splitting error is further decreased.Meanwhile,Wendland kernel is applied to avoid particle clumping[32].Both mirroring ghost boundary[30]and fixed ghost particle boundary[33]are applied since the former shows good consistency and the latter is robust for freely moving solid bodies with curved boundary shapes.In addition,the particle splitting method is combined with OpenMP parallel computing method,which has been proven to obviously accelerate the calculating speed[34].In the numerical results,a hydrostatic test case is validated firstly,and then a water discharge case is tested to show the improvement of the splitting accuracy.Followed by that,a classic dam breaking benchmark is tested to show the increase of accuracy and the decrease of spurious pressure noise in the splitting region.A wedge water entry case is also tested to show its potential applications in ocean engineering,where the numerical results agree well with the experimental data.Finally,viscous flow around rigid bluff body is simulated and validated.In the last case,besides the validation of the numerical accuracy of the SPH scheme,the improving of the efficiency is also highlighted.

2 Theoretical background

2.1 Governing equations for the fluid

The WCSPH method is applied in this paper to solve hydrostatic and hydrodynamic problems.The fluid is assumed to be barotropic and the evolution of the inner energy can be ignored[17].The governing equations in Lagrangian form are written as follows[35]

In Eq.(1),the symbol D/D t stands for the Lagrangian derivative.ρ,u,p and f denote density,velocity,pressure and the acceleration due to the body force,respectively.T denotes the stress tensor of Newtonian fluid,which can be expanded as follows[35]

where R=(?u+?uT)/2 denotes the tensor of strain rate and one can derive tr R= ? ·u.λ and μ denote the bulk viscosity and dynamic viscosity respectively,and λ = ?2μ/3 is usually adopted according to the Stokes hypothesis.I is the identity matrix with the order of d×d where d denotes the spatial dimensions.In the momentum equation,the term?·T can be expanded as follows[35]

where?·V=(λ +μ)? (? ·u)+μ?2u is the part of viscous force and V is the viscous stress tensor.

In Eq.(1),an equation of state is applied to uncouple the momentum equation and the continuum equation[17],that is

In Eq.(4),ρ0is the reference density when the pressure is zero.c0is the artificial sound speed,which is usually set as[16]

where pmaxis the expected maximum pressure and Umaxis the expected maximum velocity.

The idea of SPH is to discretize the continuum fluid into freely moving particles.These particles carry the generic variables like p,ρ,u,etc.and they are transported in a Lagrangian way under the pressure gradient and viscous force.After using kernel and particle approximations[13,19]and adding the density diffusive term[17],the governing equations can be written as follows

In Eq.(6),the density diffusive term added in the continuum equation is aimed to tackle the spurious pressure noise that often occurs in the WCPSH results.δ is a parameter independent on the specific problem and is usually set as 0.1.〈?ρ〉Lis the gradient of the density evaluated using the renormalized gradient equation[17].

In the SPH modeling of different hydrodynamics problems,the physical viscous force or artificial viscous force can be considered in the momentum equation. For the simulation of viscous flows where the boundary layer plays an important role, physical viscous force is often considered (see e.g. Refs.[5,36]), while in the modeling of FSI problems where the boundary layer is not relevant, artificial viscous force is usually applied to stabilize the velocity field, see e.g.Refs. [3,37].

The physical viscous force(the second term on the right hand side of the momentum equation)is similar to the one proposed in Ref.[38].This term can be written as

For the modeling of inviscid flows,the physical viscous force term can be replaces by an artificial viscous term(see more in Ref.[38])as

where μartificialis an artificial viscosity.Kμartificialcan be replaced by αc0ρ0h where the coefficient α = 0.05 is adopted in the present work.

The kernel function applied in this paper is the C2 Wendland function which has been proven to be free of particle clumping even when the smoothing length is large,see Ref.[32].The kernel function is written as follows

2.2 Governing equations for the motion of the solid body

For the freely moving solid bodies,their motions are controlled by the following equations

In Eq.(10),the subscripts F and S represent the set of fluid particles and the set of solid particles,respectively.M and IGare the mass and the moment of inertia of the body around the gravity center.U andΩ denote the translational and rotational velocity respectively.The force FF?Sand torque TF?Sapplied on the moving solid body can be derived according to the New ton’s third law.That is the resultant of forces imposed on the solid body is equal to the summation of those forces applied to the fluid particles.Based on the momentum equation in Eq.(6),one has the resultant of forces on the solid body as

where in a particle pair,i belongs to set of fluid particles and j belongs to the set of solid body particles.The torque on the center of the rigid body is written as

where rGdenotes the gravity center of the solid body.The viscosity μ in Eqs.(11)and(12)can be the physical one(see Eq.(7))or the artificial one(see Eq.(8)),depending on the specific cases.

2.3 Boundary conditions and time stepping

Mirroring ghost boundary and dummy particle boundary are applied in the present work,for more details,we recommend the readers to Refs.[30,31,33].Note that,in this paper,we have extended the dummy particle boundary to freely moving walls,which play an important role in FSI problems with irregular shapes in ocean engineering,see e.g.Ref.[2].The governing equations are explicitly integrated by the4thorder Runge-Kutta scheme,which allows a relatively large time step as[5]

In Eq.(13),hminis the minimum smoothing length inside the whole flow.

Fig.1 Illustration of the particle splitting process

3 Particle splitting technique

In this section,a particle splitting method in two dimensional spaces is introduced.When a particle denoted by n,namely the mother particle,enters the splitting region,it will be split into 4 daughter particles located on the 4 vertices of a square whose side length is denoted by βΔx,β ∈ (0,1),where β is a coefficient of the distance between two daughter particles and Δx is the distance between two mother particles,see Fig.1.We use k to denote the index of the daughter particles and the mass,density,pressure,velocity of the daughter particles can be denoted as mk,ρk,pk,uk.According to the conservation of total mass,we assume the mass of the daughter particle as mk=0.25mn,where mnis the mass of the mother particle.The velocities of the daughter particles are equal to that of the mother particle(i.e.un=uk)since it contributes to an exact conservation of momentum[22].Regarding the density,a Shepard kernel interpolation is adopted in Ref.[24]to evaluate the density from the mother particles,however,on one hand,it enlarges the amount of calculation since an additional program loop is needed,on the other hand,it may lead to the non-conservation of the total volume,therefore in a more convenient way,we set the densities of the daughter particles equal to that of the mother particle,i.e.ρn= ρk.A key factor that is relevant to the numerical accuracy is the smoothing length ratio between the mother and the daughter particle.Here the latter is defined as γ =hk/hn,γ ∈ (0,1).

In order to achieve the optimal accuracy,an error analysis has been conducted to determine the optimal β and γ .Since the derivative?f(r)is the mostly used in the governing equations,the error of?f(r)deserves more attention than that of f(r)as pointed out in Ref.[24].If one particle n in the compactly supported domain of particle i is split,then we have the new approximation of?f(ri)as follows

Therefore,one can obtain the square of the error as

More specifically,substituting the continuity equation into Eq.(15),one can obtain the splitting error of the density approximation as

The general splitting error mainly due to the following term as

Here we consider an extreme case that the coarse particles around one coarse particle are all split,therefore,the total splitting error at riis defined as

eΩ(ri)is named as the splitting error hereinafter.If we vary the coefficient β and γ and calculate each corresponding eΩ(ri),the optimal β and γ can be obtained when eΩ(ri)has the minimum value.

Fig.2 The distribution of the splitting error ln(eΩ(r i))when β varies from 0.2 to 0.8 and γ varies from 0.1 to 1

Similar to the numerical analysis in Ref.[24],in the initial condition,all the particles are coarse particles and the initial smoothing length is set as hcoarse=1 and therefore Δx=1/κ.Here two frequently-used κ are tested,i.e.κ =1.23 and κ=2 since in different cases,one may adopt different κ to achieve a balance between the numerical accuracy and efficiency,see Ref.[35].In our scheme, we choose one coarse particle and split all the other particles adjacent to it.We split the particle with β varying from 0.2 to0.8andγ varying from 0.1 to1.Based on this testing matrix,one can easily obtain the optimal β and γ which make the minimum value of eΩ(ri).In order to show the distribution of variation of eΩ(ri)more clearly,ln(eΩ(ri))is shown in Fig.2,from which one may find that,as γ increases,ln(eΩ(ri))is decreasing.For κ =1.23,the minimum ln(eΩ(ri))locates at β =0.4andγ =0.9 or β =0.5 and γ =0.9 or β =0.6 and γ =0.8 or β =0.7 andγ =0.7.Forκ =2.0,the minimum ln(eΩ(ri))locates at β=0.7 and γ=0.9.Note that,the above parameters only contribute to the minimum error at the splitting interface;while in the rest of the refined region,choosingw ill make the daughter particles in a disordered distribution.Based on our numerical tests,even we choose,after several steps,the particles tends to a uniform distribution(i.e.β=0.5),which is due to the particle packing mechanism,see Ref.[39].For the particle splitting occurring near the fixed ghost particle boundary,both the fluid particles and the ghost particles are split,if,the ghost particles would be non-uniform forever since they are fixed while the fluid particles will be regularized during the SPH simulation.In that case,the particle disorder takes place near the boundary and worsens the numerical accuracy.Therefore,β=0.5 should be maintained.Then regarding γ,from Fig.2,wefind that,γ＞0.85 shows relatively small error for both κ=1.23 and κ=2.Note that,since the Wendland kernel is applied,the particle pairing due to the large smoothing length ratio can be avoided,see Ref.[32].

4 Numerical results and discussions

4.1 A hydrostatic case

Fig.3 Illustration of the initial condition for the hydrostatic case

In the SPH modeling of water wave propagation problems in ocean engineering,the dynamic effects are not so apparent below the fluid surface and most particles are in quasi hydrostatic state,therefore small errors in the flow may lead to significant spurious particle motions.Besides,the simulating duration is usually quite long in these cases and therefore small errors may grow up and worsen the SPH results.Therefore if particle splitting technique is implemented,the splitting error should be quite minor in order to keep the splitting interface stable.Similar to Ref.[28],a simple hydrostatic test case is applied in this part to test the ability of the present SPH method in maintaining the hydrostatic pressure. It needs to be noted that the present test case is simple but quite difficult for SPH models.To the knowledge of us,only Barcarolo et al.[5]give good results regarding this problem in the framework of the Riemann-SPH method.

The initial condition for the hydrostatic test is shown in Fig.3.The depth of the stationary fluid is H=1 m and the width of the water tank is W=2 m.The fluid is discretised with the particles pacing Δx=0.01 m at the initial stage and the smoothing length is set as h=2Δx.The initial pressure inside the flow is p=ρ0ghdwhere hdis the fluid depth.In all the test cases of the present work,the reference density ρ0is set as the density of pure water,i.e.ρ0=1000 kg/m3and the gravity acceleration is g=9.8 m/s2.After the pressure p is calculated,we can determine the density ρ based on the equation of state.Regarding the boundary method,the mirroring ghost boundary is applied in this test case since it allows the highest consistency near the solid wall boundary.

Fig.4 The comparison of the numerical results between γ =0.7 and γ =1.0 at different time instants for the hydrostatic test case

4.2 A water discharging case

Fig.5 Illustration of the initial condition for the waterd is charging case

In this part,similar to Ref.[24],a water discharging case is tested.The difference is that,in López’s work,a partly overlapping particle distribution for the daughter particles is applied which was proven to significantly improve the numerical accuracy;while in the present work,we maintain the daughter particles to be uniform but improve their smoothing lengths,which also contributes to a higher accuracy.The initial condition for this case is shown in Fig.5 where the depth of the undisturbed water is H=2 m and the width is W=1 m.The height of the gate at the right bottom of the tank is H′=0.4 m.After releasing the gate at t=0,the water flows out under the gravity force.An out flow condition is applied at the gate through which any particle flows out is transferred to a certain location with all the physical variables set to zero.

Fig.6 The discharging of the fluid at different time instants:on the left is the fully refined case with Δx=0.02 m and on the right is the partially refined case with β =0.5,γ =0.9

Firstly the simulations are carried out with a single particle resolution(without particle splitting),two particle resolutions are used.Δx=0.04 m is defined as the unrefined particle resolution and Δx=0.02 m is defined as the fully refined particle resolution.The smoothing length is set as h=1.23Δx at the initial stage.Dummy particle boundary is applied to treat the solid wall boundary.And then in the partially refined simulations,based on the case of Δx=0.04 m,we split the particles when they enter the splitting region with the parameters of β =0.5,γ=0.7 and β =0.5,γ =0.9,and compare the results between them.One can imagine that since the height of the gate is restricted,the smaller the particle is,the faster the fluid surface drops.

The comparisons between the flows of the full refined case and the partially refined case of β =0.5,γ =0.9 at different time instants are shown in Fig.6,from which one may find that as the coarse particles enter the splitting region,they are split into daughter particles and then flow out.Therefore,the out flow speed of the partially refined cases should be similar to that of the fully refined case.In order to more clearly show the effects of the particle resolution,we compare the time evolution of the water height in Fig.7,from which one may find that,for the case of single resolution,the water height of the unrefined case drops slower than that of the fully refined case and the results of the partially refined case are between them.One should also notice the difference between the results of γ =0.9 and γ =0.7 that the former is more close to the fully refined result,which verifies the conclusion in Sect.3 that,with β =0.5,the splitting error of γ =0.9 is smaller than that of γ=0.7.

4.3 Dam breaking test

Dam breaking test is a classic benchmark which has been widely adopted by SPH researchers to validate their numerical schemes,e.g.Refs.[30,31].The initial condition is illustrated in Fig.8,in which all the parameters for the size of the fluid and the water tank are labeled.The initial pressure distribution is obtained by the solution based on an incompressible hypothesis,see more details in Ref.[40].Time evolution of the pressure load on the pressure transducer at the height of H′=0.266 m will be measured and compared against the experimental data[41].

Fig.7 The time evolution of the water height under different particle resolutions

Fig.8 The initial condition for the dam-breaking benchmark

Fig.9 On the left shows the time evolution of the pressure measured on the transducer;SPH results are compared with the experimental data in Ref.[41];on the right side shows the time evolution of the total mechanical energy inside the flow

Fig.10 The flow features of the partially refined dam-breaking case at different time instants

In the practical applications,the pressure load on the solid wall makes sense for the evaluation of the structure strength.The particle resolution of H=100Δx is shown to be enough for the pressure approximation.Therefore it is named as the fully refined case and H=50Δx is defined as the unrefined case.If we split the particles which enter the splitting domain shown in Fig.8,we can on one hand achieve a better pressure time history curve on the solid wall and on the other hand improve the computational efficiency.This case with the particle splitting is named as partially refined case.Similar to Sect.4.2,two groups of splitting parameters β =0.5,γ =0.7 and β =0.5,γ=0.9 are tested and the results will be compared against each other.

The flow features of the partially refined case are shown in Fig.10,from which we can find that,after the particles on the

Fig.11 On the left shows the time evolution of the pressure measured on the transducer;SPH results are compared with the experimental data in Ref.[41];on the right side shows the time evolution of the total mechanical energy inside the flow

Fig.12 The initial condition for the wedge water entry case

4.4 The wedge water entry

In this part,a wedge water entry case is tested through the comparison of SPH results against the experiment data in Ref.[42].The initial condition for this problem is depicted in Fig.12.The mass of the wedge is M=94 kg.The fluid density is 1000 kg/m3and the gravity acceleration is 9.8 m/s2.The depth H and width W of the tank is determined as 1.1 m and 2.5 m respectively.The width of the upper panel of the wedge is W′=1.2 m and the dead rise angle is θ=25°.The wedge enters the undisturbed water surface with the initial velocity as u0=5.05 m/s.The splitting region is defined with the width of Ws=1.4 m and the bottom of the splitting area is located at hs=hG?0.4 m where hGis the height between the gravity center G of the wedge and the tank bottom.Therefore,as the wedge falls,the splitting region also falls following the wedge,and a dynamic particle refinement is achieved.

The unrefined case is defined with Δx=10 mm and the fully refined case is defined with Δx=5 mm.In the partially refined case,we split the particles in the splitting region from the particle spacing ofΔxn=10 mm into Δxk=5 mm.The smoothing length is set as h=1.23Δx for the unrefined particles.The parameters used for the particle splitting technique isβ=0.5andγ=0.9.Snapshots of the water entry process of the wedge are depicted in Fig.13 for three time instants.The results show that the splitting region moves following the motion of the wedge,which contributes to a dynamic particle refinement in SPH.

In order to more clearly show the effects of the particle resolution on the final numerical results,comparisons for the wedge penetration velocities of different particle resolutions are made in Fig.14.As the particle resolution is refined,the penetrating velocity converges to the experimental data[42].In addition,the result of the partially refined case is very close to that of the fully refined case.

4.5 Viscous flow around a circular cylinder

In this section,viscous flow around a circular cylinder is tested using the SPH model with the particle splitting technique to reduce the overall particle number and therefore accelerate the computational speed.Different to the previous cases,a physical viscosity is adopted in this case.The Reynolds number,defined as Re= ρUD/μ where U is the free stream velocity and D is the cylinder diameter,is set as 185.

Fig.13 The dropping of the wedge,the translational motion of the splitting region and the pressure distribution inside the flow domain at different time instants

Fig.14 Time evolutions of the penetrating velocity of the wedge.SPH results with different particle resolutions are compared against the experimental data in Ref.[42]

Fig.15 Illustration for the set-up the SPH simulation of the viscous flow around a fixed circular cylinder at Re=185

In the simulation,the center of the circular cylinder is located at the origin of the frame of reference.The in flow boundary(see Refs.[43,44])is set at x=?8D and outflow boundary at x=24D,while the two lateral boundaries imposed by free-slip conditions are arranged at y=±8D(see Fig.15).Under this condition,the flow is wide enough to avoid the blockage effect for this problem,see Ref.[36].The particle spacing in the in flow buffer zone is Δxn=D/50 and the smoothing length is h=1.23Δxn.After these particles enter the splitting domain which covers the fixed circular cylinder,they are split into daughter particles with the spacing of Δxk= Δxn/2,i.e.Δxk=D/100 and these daughter particles can freely exit the splitting domain.The smoothing length parameter for the daughter particle is γ=0.9.A particle shifting technique has been applied for all the particles in order to prevent the tensile instability,see more in Refs.[16,45,46].

Fig.16 The mass distribution at time instant tU/D=38.88 in the viscous flow around a circular cylinder at Re=185

Fig.17 Time evolution of the drag and lift force coefficients measured on the cylinder for the viscous flow over a cylinder at Re=185

In order to rapidly obtain a periodical vortex shedding,when tU/D≤5.0,the upper half of the cylinder is imposed by a free-slip condition while the lower half is imposed by a no-slip condition and when tU/D＞5.0,the whole cylinder surface is imposed by a no-slip condition.The particle mass distribution at the time instant tU/D=38.88,when the vortex shedding reaches a steady stage, is depicted in Fig. 16.The mass of the daughter particle is about a quarter of the one of mother particle.It is observed that the daughter particles are distributed along the vortex structure in the wake with the shape resembling the von Karman vortex street.

Table 1 Comparison of the mean drag and the r.m.s.lift force coefficients for the viscous flow over a cylinder at Re=185

Table 2 Comparison of the computational cost between the cases of fully refined and partially refined particle resolutions

5 Conclusion

In the present work,a further extended particle splitting method is introduced and tested.The parameter β is set to be 0.5 in order to allow for a uniform distribution of the daughter particles.Regarding the smoothing length parameter γ,as γ increases,the splitting error is decreasing,for both h=1.23Δx and h=2Δx,γ ＞ 0.85 has to be maintained for a satisfactory accuracy.

We test the particle splitting technique through both hydrostatic and hydrodynamic cases.The hydrostatic case is challenging when a particle splitting is imposed;while the present splitting method is demonstrated to be able to give a stable multi-resolution interface.Three hydrodynamic test cases are also tested,respectively water discharging,dambreaking and wedge water entry,through which the effects of the splitting parameters are tested and compared.The last test is the viscous flow around a fixed circular cylinder at Re=185.Since a free stream flow is modeled,a large flow domain is requested and consequently the particle splitting in a local region contributes to a significant reduction of the total computational cost.The reduction rate reaches about 63.3%.

In conclusion,in the present work a preliminary attempt has been made to apply the particle splitting technique in the SPH modeling of different hydrodynamic problems.The results show that,after using the present particle splitting method,the numerical results of partially refined cases are similar to those of the fully refined results,while the computational cost is reduced.

AcknowledgementsThe project was supported by the National Natural Science Foundation of China (Grant 51609049), the Science Foundation of Heilongjiang Province(Grant QC2016061),and the Fundamental Research Funds for the Central Universities(Grants HEUGIP201701,HEUCFJ170109).