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.

场地全部钻孔（CK1～CK10，ZK1～ZK21）均有揭露，揭露层厚 1.30～4.50m，平均厚度2.66m；层顶高程3.40～3.83m，层顶深度0.00～0.00m。

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]

fl ow tongue enter the splitting region,they are refined into smaller particles.From the pressure time history as shown in Fig.11,one may find that the results of the partially refined case almost coincide with those of the fully refined case and moreover,the pressure evolution is even more stable.Comparing the results of the partially refined cases between γ=0.7 and γ=0.9,the latter case shows better agreement regarding the pressure magnitude when compared with the fully refined case.The dissipation of the total mechanical energy is plotted on the right side of Fig.11,which shows that the results of the partially refined case are between the results of unrefined and fully refined cases.Besides,comparing the results between γ =0.7 and γ=0.9,it is noted that after aboutthe result of γ=0.9 dissipates more mechanical energy than that of γ=0.7.The reason is that,the magnitude of the artificial viscosity is proportion alto the smoothing length,and therefore the latter one dissipates more kinetic energy due to larger viscous dissipation.From this case,it has been demonstrated that,the partially refined cases give stable pressures in local regions and the mechanical dissipation is also smaller than the unrefined case.The numerical result implies that the particle splitting technique can be a potential way to more accurately and efficiently evaluate the local pressure load in some violent FSI problems.

对于无法避免的含有害物质的材料和零部件，可采用适当加热的方法促进有害物质提前释放，或者加大采购量，在自然通风条件下放置一段时间，使有毒有害物质充分释放后再组装。

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

where r G denotes 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.

来一场说走就走的旅行，不是任性，而是实属无奈。原本手头上的事情太多，十一长假也准备好好“奋斗”一番，谁想一通电话的相约，孩子欢呼雀跃吵着要走，如一枚石子骤然敲破平镜的湖面。计划乱了，但心之所动，已飘向了旅途的远方。

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

至今仍然有很多老师对于题海战术训练乐此不疲，长期以往，取得的效果却微乎其微、甚至相反.而随着信息科技的高速发展，初中数学习题课教学绝不应受到应试教育的影响而仅仅局限于“老师一张嘴巴讲到尾”或者“掐头去尾烧中段”的题海战术套路，而更应该教会学生学以致用、培养学生勇于创新的精神，于是变式教学应运而生，它以精心设计问题、摒弃题海战术、引导探索发现、展现知识形成过程、注重知识建构、培养创新能力为主要目标，并由变式问题拓展延伸出其他知识点，促进学生知识的迁移.

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.〈∇ρ〉L is 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].

where p max is the expected maximum pressure and U max is the expected maximum velocity.

where andis the physical viscosity of the fluid.K equals to 2(d+2)where d is the spatial dimensions as stated before.Equation(7)converges to the∇ ·V defined in Eq.(3)when the fluid tends to be incompressible,i.e.∇·u→ 0,see more in Ref.[35].

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 μarti fi cial is an artificial viscosity.Kμarti fi cial can 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

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.

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 IG are 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 F F−S and torque T F−S applied 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

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

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 the4th order Runge-Kutta scheme,which allows a relatively large time step as[5]

In Eq.(13),h m in is the minimum smoothing length inside the whole flow.

轮子直径66mm，光电码盘齿数为20，轮子周长207mm=20.7cm，轮子转一周的脉冲信号计数次数为40，一个计数变化表示轮子跑过的距离为20.7/40=0.5175cm。

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,u k.According to the conservation of total mass,we assume the mass of the daughter particle as mk=0.25mn,where mn is the mass of the mother particle.The velocities of the daughter particles are equal to that of the mother particle(i.e.u n=u k)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(r i)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 r i is defined as

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

试验组45例患者中，显效26例（57.8%），有效16例（35.6%），无效3例（6.7%），临床总有效率为93.3%（42/45）；对照组45例患者中，显效18例（40.0%），有效16例（35.6%），无效11例（24.4%），临床总有效率为75.6%（34/45），组间对比，差异具有统计学意义（χ2=5.413 5，P=0.019 9＜0.05）。

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 h coarse=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Ω(r i).In order to show the distribution of variation of eΩ(r i)more clearly,ln(eΩ(r i))is shown in Fig.2,from which one may find that,as γ increases,ln(eΩ(r i))is decreasing.For κ =1.23,the minimum ln(eΩ(r i))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Ω(r i))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=ρ0 gh d where h d is the fluid depth.In all the test cases of the present work,the reference density ρ0 is set as the density of pure water,i.e.ρ0=1000 kg/m3 and 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

A square splitting region with the side length of L=0.5m is defined as shown in Fig.3.Once a coarse particle enters the splitting region,it is split into 4 daughter particles and the coefficient of relative location β is set to be 0.5.Regarding γ,both γ =0.7 and γ =1.0 are tested.In Fig.4,a comparison of the splitting interfaces is shown.One may find that,for γ=0.7,at the two right angles in the bottom of the splitting region,the refined particles start to penetrate into the coarse particles atand finally atthe splitting interface has been quited is order;while forγ=1.0,the splitting region is quite stable,which is similar to the results introduced in Ref.[28]by the Riemann-SPH method,which demonstrates that for the hydrostatic cases,γ=1.0 is necessary since it allows smaller splitting error.Besides,since the Wendland kernel is applied here,there is no particle clumping which is usually excited by a large ratio between the smoothing length and the particle spacing.

4.2 A water discharging case

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

where

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

①寒武系石牌组（Є2s）：灰色页岩，呈页片状，岩石由粘土质矿物及少量陆屑石英组成。厚 40～50 m。

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

Here we use three different particle resolutions:H =50Δx,H=100Δx and H=200Δx to test the convergence of the present SPH model.The smoothing length is set as h=1.23Δx at the initial stage.The results are shown in Fig.9.From the pressure-time history one may find that,as the particle resolution is refined,the pressure noise decreases and the result of the higher particle resolution shows better agreement with the experimental data[41].For the particle resolutions of H = 100Δx and H=200Δx,the pressure-time curves almost coincide.The ratiowhere E pot isthe total potential energy,E mech is the total kinetic energy and denote the potential energy at the initial stage,is defined to measure the dissipation of the total mechanical energy.From the curves on the right side of Fig.9,one may find that as the particle resolution is refined,the energy dissipation is reduced.

上述公式中，锁外消耗时长表示为tp，其占比重表示为q；解及加锁消耗时长表示为tι，锁内消耗时长表示为ts。当q是整数的情况下，意为首次申请才会形成锁竞争，若q是小数，那么会重复形成锁竞争。运算并行下载时间的运算公式如下：

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

where R=(∇u+∇u T)/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]

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/m3 and 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 W s=1.4 m and the bottom of the splitting area is located at h s=hG−0.4 m where hG is 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.

盘山石结成，松峪独有土。 此中即沃壤，开田作场圃。 黄犊晓耕耘，绿衰春 带雨。 遥听布谷鸣，杏花满村坞。[3]16

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.

姚正安先生《熟人生处》讲，因是熟人而不讲礼节，招致对方反感。告诫人们不管是熟人还是陌生人都得以礼相待，揭示的也是待人接物的人生哲理。

The drag and lift force coefficients,defined as CD=respectively,aremeasured on the cylinder surface and the force coefficients during tU/D∈(35,75)have been plotted in Fig.17.The time evolution of the forces implies that in this time range,a periodical vortex shedding stage is reached.The mean drag and r.m.s.lift force coefficients are compared against the solutions from the literature in Table 1.All the solutions are comparable apart from some small discrepancies which may be due to the different mesh/particle resolution adopted in each simulation.As shown in Table 2,the computational costs are also compared between the cases with fully refined particle resolution(the whole flow is discretised with the particle spacing Δx=D/100)and partially refined particle resolution(using the particle splitting technique).In Table 2,T0 denotes the computational time within one time step of simulation with fully refined resolution.The computational cost of the partially refined one is only 36.7%of the fully refined one.In other words,the reduction rate reaches about 63.3%.

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

Fully refined Partially refined Computational cost T0 36.7%T0

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.

Acknowledgements The 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).

计算机软件在实际应用中有着比较突出的优点：（1）在社会应用中比较普遍，得到了较大的认可，同时在很大程度上也在进一步推动计算机软件产业的不断快速发展和进步，经济效益呈现出了多元化的状态。（2）软件开发工作实际上属于一种较强的系统性工作，工作极为细致和复杂，这都需要消耗大量的人力、物力、财力，但是一旦开发并推广使用，对社会发展的推动作用就不可忽视了[1]。

References

1.Marrone,S.,Bouscasse,B.,Colagrossi,A.,et al.:Study of ship wave breaking patterns using 3D parallel SPH simulations.Comput.Fluids 69,54–66(2012)

2.Sun,P.,Ming,F.,Zhang,A.:Numerical simulation of interactions between free surface and rigid body using a robust SPH method.Ocean Eng.98,32–49(2015)

3.Zhang,A.-M.,Sun,P.-N.,Ming,F.-R.,et al.:Smoothed particle hydrodynamics and its applications in fluid–structure interactions.J.Hydrodyn.Ser.B 29,187–216(2017)

4.Shao,J.R.,Li,H.Q.,Liu,G.R.,et al.:An improved SPH method for modeling liquid sloshing dynam ics.Comput.Struct.100–101,18–26(2012)

5.Bouscasse,B.,Colagrossi,A.,Marrone,S.,et al.:Nonlinear water wave interaction with floating bodies in SPH.J.Fluids Struct.42,112–129(2013)

6.Jiang,T.,Lu,L.G.,Lu,W.G.:The numerical investigation of spreading process of two viscoelastic droplets impact problem by using an improved SPH scheme.Comput.Mech.53,977–999(2014)

7.Jiang,T.,Ren,J.L.,Lu,W.G.,et al.:A corrected particle method with high-order Taylor expansion for solving the viscoelastic fluid lf ow.Acta Mech.Sin.33,1–20(2017)

8.Ming,F.R.,Sun,P.N.,Zhang,A.M.:Numerical investigation of rising bubbles bursting at a free surface through a multiphase SPH model.Meccanica 52,2665–2684(2017)

9.Zhang,A.,Sun,P.,Ming,F.:An SPH modeling of bubble rising and coalescing in threedimensions.Comput.Methods Appl.Mech.Eng.294,189–209(2015)

10.Ming,F.R.,Zhang,A.M.,Xue,Y.Z.,et al.:Damage characteristics of ship structures subjected to shockwaves of underwater contact explosions.Ocean Eng.117,359–382(2016)

11.Zhang,A.M.,Yang,W.S.,Huang,C.,et al.:Numerical simulation of column charge underwater explosion based on SPH and BEM combination.Comput.Fluids 71,169–178(2013)

12.Zhang,Z.,Wang,L.,Silberschmidt,V.V.,et al.:SPH–FEM simulation of shaped-charge jet penetration into double hull:a comparison study for steel and SPS.Compos.Struct.155,135–144(2016)

13.Liu,M.B.,Liu,G.R.:Smoothed particle hydrodynamics(SPH):an overview and recent developments.Arch.Comput.Methods Eng.17,25–76(2010)

14.Ren,J.,Jiang,T.,Lu,W.,et al.:An improved parallel SPH approach to solve 3D transient generalized Newtonian free surface flows.Comput.Phys.Commun.205,87–105(2016)

15.Gan,B.S.,Nguyen,D.K.,Han,A.L.,et al.:Proposal for fast calculation of particle interactions in SPH simulations.Comput.Fluids 104,20–29(2014)

16.Sun,P.,Colagrossi,A.,Marrone,S.,et al.:The δplus-SPH model:simple procedures for a further improvement of the SPH scheme.Comput.Methods Appl.Mech.Eng.315,25–49(2017)

17.Antuono,M.,Colagrossi,A.,Marrone,S.:Numerical diffusive terms in weakly-compressible SPH schemes.Comput.Phys.Commun.183,2570–2580(2012)

18.Barcarolo,D.A.Improvement of the precision and the efficiency of the SPH method:theoretical and numerical study.[Ph.D.Thesis],Centrale Nantes,France(2013)

19.Liu,G.R.,Liu,M.B.:Smoothed Particle Hydrodynamics:A Meshfree Particle Method.World Scientific,Singapore(2003)

20.Marsh,A.,Oger,G.,Touzé,D.l.,et al.:Validation of a conservative variable-resolution SPH scheme including∇h terms.In:The 6th International SPHERIC Workshop,Hambourg,Germany(2011)

21.Koukouvinis,P.K.,Anagnostopoulos,J.S.,Papantonis,D.E.:Simulation of 2D wedge impacts on water using the SPH–ALE method.Acta Mech.Sin.224,2559–2575(2013)

22.Feldman,J.,Bonet,J.:Dynamic refinement and boundary contact forces in SPH with applications in fluid flow problems.Int.J.Numer.Methods Eng.72,295–324(2010)

23.Omidvar,P.,Stansby,P.K.,Rogers,B.D.:Wave body interaction in 2D using smoothed particle hydrodynamics(SPH)with variable particle mass.Int.J.Numer.Methods Fluids 68,686–705(2012)

24.López,Y.R.,Roose,D.,Morfa,C.R.:Dynamic particle refinement in SPH:application to free surface flow and non-cohesive soil simulations.Comput.Mech.51,731–741(2013)

25.Antuono,M.,Bouscasse,B.,Colagrossi,A.,et al.:A measure of spatial disorder in particle methods.Comput.Phys.Commun.185,2609–2621(2014)

26.Vacondio,R.,Rogers,B.D.,Stansby,P.K.,et al.:Variable resolution for SPH:a dynamic particle coalescing and splitting scheme.Comput.Methods Appl.Mech.Eng.256,132–148(2013)

27.Omidvar,P.,Stansby,P.K.,Rogers,B.D.:SPH for3D floating bodies using variable mass particle distribution.Int.J.Numer.Methods Fluids 72,427–452(2013)

28.Barcarolo,D.A.,Oger,G.,Vuyst,F.D.:Adaptive particle refinement and derefinement applied to the smoothed particle hydrodynamics method.J.Comput.Phys.273,640–657(2014)

29.Monaghan,J.J.:Simulating free surface flows with SPH.J.Comput.Phys.110,399–406(1994)

30.Colagrossi,A.,Landrini,M.:Numerical simulation of interfacial flows by smoothed particle hydrodynamics.J.Comput.Phys.191,448–475(2003)

31.Marrone,S.,Antuono,M.,Colagrossi,A.,et al.:δ-SPH model for simulating violent impact flows.Comput.Methods Appl.Mech.Eng.200,1526–1542(2011)

32.Dehnen,W.,Aly,H.:Improving convergence in smoothed particle hydrodynamics simulations without pairing instability.Mon.Not.R.Astron.Soc.425,1068–1082(2012)

33.Adami,S.,Hu,X.Y.,Adams,N.A.:A generalized wall boundary condition for smoothed particle hydrodynamics.J.Comput.Phys.231,7057–7075(2012)

34.Zhang,A.M.,Cao,X.Y.,Ming,F.R.,et al.:Investigation on a damaged ship model sinking into water based on three dimensional SPH method.Appl.Ocean Res.42,24–31(2013)

35.Colagrossi,A.,Soutoiglesias,A.,Antuono,M.,et al.:Smoothed particle hydrodynamics modeling of dissipation mechanisms in gravity waves.Phys.Rev.E Stat.Nonlin.Soft Matter Phys.87,023302(2013)

36.Marrone,S.,Colagrossi,A.,Antuono,M.,et al.:An accurate SPH modeling of viscous flows around bodies at low and moderate Reynolds numbers.J.Comput.Phys.245,456–475(2013)

37.Guo,K.,Sun,P.-N.,Cao,X.-Y.,et al.:A 3-D SPH model for simulating water flooding of a damaged floating structure.J.Hydrodyn.Ser.B 29,831–844(2017)

38.Monaghan,J.J.,Gingold,R.A.:Shock simulation by the particle method SPH.J.Comput.Phys.52,374–389(1983)

39.Colagrossi,A.,Bouscasse,B.,Antuono,M.,et al.:Particle packing algorithm for SPH schemes.Comput.Phys.Commun.183,1641–1653(2012)

40.Greco,M.:A two-dimensional study of green-water loading.[Ph.D.Thesis],Norwegian University of Science and Technology,Norway(2001)

41.Buchner,B.Green water on ship-type offshore structures.[Ph.D.Thesis],Technische Universiteit Delft,Netherlands(2002)

42.Yettou,E.M.,Desrochers,A.,Champoux,Y.:Experimental study on the water impact of a symmetrical wedge.Fluid Dyn.Res.38,47–66(2006)

43.Federico,I.,Marrone,S.,Colagrossi,A.,et al.:Simulating 2D open-channel flows through an SPH model.Eur.J.Mech.34,35–46(2012)

44.Sun,P.,Colagrossi,A.,Marrone,S.,et al.:Detection of Lagrangian coherent structures in the SPH framework.Comput.Methods Appl.Mech.Eng.305,849–868(2016)

45.Lind,S.,Xu,R.,Stansby,P.,et al.:Incompressible smoothed particle hydrodynamics for free-surface flows:a generalised diffusion based algorithm for stability and validations for impulsive flows and propagating waves.J.Comput.Phys.231,1499–1523(2012)

46.Huang,C.,Zhang,D.H.,Shi,Y.X.,et al.:Coupled finite particle method with a modified particle shifting technology.Int.J.Numeri.Methods Eng.113,197–207(2018).https://doi.org/10.1002/nme.5608

47.Pinelli,A.,Naqavi,I.Z.,Piomelli,U.,et al.:Immersed-boundary methods for general finite-difference and finite-volume Navier–Stokes solvers.J.Comput.Phys.229,9073–9091(2010)

48.Vanella,M.,Balaras,E.:A moving-least-squares reconstruction for embedded-boundary formulations.J.Comput.Phys.228,6617–6628(2009)

49.Cai,S.G.,Ouahsine,A.,Favier,J.,et al.:Improved implicit immersed boundary method via operator splitting.In:Ibrahimbegovic,A.(ed.)Computational Methods for Solids and Fluids.Computational Methods in Applied Sciences,vol.41.Springer,Cham(2016)

50.Constant,E.,Favier,J.,Meldi,M.,et al.:An immersed boundary method in Open FOAM:verification and validation.Comput.Fluids 157,55–72(2017)