An improved predictive-corrective incompressible smoothed particle hydrodynamics method for fluid flow modelling ＊

Chong Peng , Christoph Bauinger Kamil Szewc Wei Wu, Hui Cao

1. ESS Engineering Software Steyr GmbH, Berggasse 35, 4400, Steyr, Austria 2. Institut für Geotechnik, Universität für Bodenkultur, Feistmantelstrasse 4, 1180, Vienna, Austria

Abstract: Predictive-corrective incompressible smoothed particle hydrodynamics (PCISPH) is an efficient SPH variation originally developed for computer graphics. Its application in modeling physics-focused fluid flows has not yet been reported. In this work, an improved PCISPH method is presented for physics-focused fluid flow modeling. Different from traditional weakly-compressible SPH(WCSPH) and incompressible SPH (ISPH), PCISPH satisfies the incompressibility requirement at the particle level without a predefined equation of state. The pressure is obtained using an iterative predictive-corrective scheme at each individual particle. The presented PCISPH allows much larger time steps compared with WCSPH. It also avoids the time-consuming solution of Pressure Poisson Equation (PPE) in ISPH. Consequently, the PCISPH has high computational efficiency even with millions of computational particles. To ensure physically correct modeling of fluid flows, we employ several techniques to enhance the accuracy and stability of the PCISPH: (1) the continuity equation is used to predict density, replacing the mass summation approach used in the original PCISPH,(2) numerical diffusion and pressure smoothing are introduced to improve the pressure computation, (3) a generalized boundary treatment approach which can handle arbitrarily complex geometries is employed, (4) an adaptive time-stepping algorithm is used,allowing efficient simulation as well as ensuring stability. The performance of the improved PCISPH is systematically investigated using three standard SPH validation cases. Comparisons between the improved PCISPH and the state-of-the-art δ-SPH are presented.It is found that the improved PCISPH gives numerical results as accurate as δ-SPH, except for having moderate temporal pressure oscillations. However, the numerical results show that the improved PCISPH is approximately five times faster than δ-SPH. The improved PCISPH method shows to be a promising tool for large-scale three-dimensional fluid flow modeling.

Key words: Smoothed particle hydrodynamics, incompressibility, predictive-corrective scheme, high efficiency

Introduction

Smoothed particle hydrodynamics (SPH) is a meshless method developed for astrophysical research.Since its invention, its application has been expanded to many areas, e.g., computational fluid dynamics(CFD)[1], fracture simulation[2], and geophysics[3].Among them, SPH probably gains most popularity in CFD, especially in modeling free-surface flows[1]. The reasons are mainly twofold: (1) SPH uses Lagrangian particles, which conserve mass exactly and model free-surfaces naturally, (2) its meshless property makes the preparation of numerical models simple compared with mesh-based numerical methods.

In CFD, enforcement of incompressibility constraints is crucial. Based on how to ensure incompressibility, most of the SPH variations can be grouped into two categories: weakly-compressible SPH (WCSPH) and incompressible SPH (ISPH). In WCSPH, incompressible flows are modelled by treating fluids as weakly compressible with an equation of state, which relates pressure to density variation. WCSPH has simple formulations, and complex physics can be modeled[4]. Therefore,WCSPH finds applications not only in CFD, but also in other areas such as solid mechanics[4], high speed impacts[5], explosions[6], and granular flows[7], where materials are treated in a similar weakly-compressible manner. Despite the wide application, WCSPH has two weaknesses. Firstly, the time step in WCSPH is limited by the Courant-Friedrichs-Lewy (CFL)condition. According to the CFL condition, the time step in WCSPH is mainly controlled by the speed of sound. This usually gives rise to very small time steps,leading to low efficiency. Secondly, in WCSPH, the pressure results suffer from spurious oscillations,leading to unreliable results and numerical instability.

The ISPH method is proposed as an alternative to WCSPH[8]. In ISPH, the pressure is obtained by solving a Poisson equation that ensures a diver- gence-free velocity field. ISPH allows larger time steps and produces smooth pressure fields. However, the large time step does not necessarily guarantee ISPH higher efficiency, because the solution of the linear system resulted from the Poisson equation is rather time-consuming. For simulations with a very large number of particles, ISPH can be even less efficient than WCSPH[9].

In computer graphics, an efficient SPH variation was proposed by Solenthaler and Pajarola[10]. In this SPH variation, the incompressibility is enforced by iteratively predicting and correcting the density. The pressure is determined in a way that the pressure forces can reduce the density fluctuation. The predictive-corrective incompressible SPH (PCISPH)allows larger time step size like ISPH. Moreover,PCISPH solves the incompressibility at the particle level, so there is no need to solve Poisson equations.Therefore, its computational cost is lower than ISPH.Like WCSPH, PCISPH is also suitable for parallelization. It combines the advantages of WCSPH and ISPH while avoiding their shortcomings.

Although PCISPH is used in computer graphics,its application in CFD focusing on flow physics has not yet been reported. In computer graphics, one focuses only on flow dynamics. However, the CFD community can embrace PCISPH as a reliable tool only if it delivers physically correct results. In this paper, we propose several improvements to enhance the accuracy and stability of PCISPH. First, instead of using the mass summation, we employ the continuity equation for the density prediction. Second, numerical diffusion of density and pressure smoothing are employed to enhance pressure accuracy and numerical stability. Third, a generalized boundary treatment method based on the Adami boundary condition[11]is introduced in PCISPH. Besides, a dynamic time stepping method is employed, which achieves large time steps while keeping simulations accurate and stable. For the first time, standard SPH validation cases are simulated using PCISPH, the numerical results are compared with experimental data and numerical results in literature. The accuracy, stability and efficiency of PCISPH in modeling realistic flows are investigated by comparing it with the state-ofthe-art WCSPH schemeδ-SPH.

1. Smoothed particle hydrodynamics fundamentals

The governing equations for viscous flows include the continuity equation and the momentum equation read as follows:

whereρis the fluid density,uis the velocity,pthe pressure,νthe kinematic viscosity andextbis the acceleration induced by external forces.

In standard WCSPH, the governing equations have the following SPH forms

whereWis the kernel function, ∇Wis the gradient of the kernel function,ηis a regularizing constant usually taken as 0.1h, withhthe soothing length. In Eq. (4), the first and second terms are accelerations induced by pressure gradient and viscous effect, respectively. The above SPH formulation is suitable for laminar flows[12]. An alternative form employing the artificial viscosity can be written as[4]

whereαis a tuning parameter,csis the speed of sound. The artificial viscosity can stabilize computation of flows with violent impacts and shocks.

An equation of state (EOS)p=p(ρ) is employed to close the governing equations where0ρdenotes the reference density,γis a parameter often chosen as 7. Usually, the speed of sound is at least ten times greater than the magnitude of the maximum fluid velocity to ensure that the density fluctuation is lower than 1%.

2. Predictive-corrective SPH

2.1 Predictive-corrective procedure with splitting

In WCSPH, a given speed of soundcsis used to achieve a certain compressibility. The speed of soundcsis like a stiffness coefficient that scales the pressure and pressure forces. A higher speed of sound results in a stricter compressibility constraint but leads to smaller time steps. ISPH computes the pressure field by solving a Poisson equation which ensures divergence-free velocity field. Although ISPH admits larger time steps, it needs higher computational effort to solve the Poisson equation. The main concept of PCISPH is to replace the global Poisson equation in ISPH with a local predictive-corrective procedure;thus, PCISPH avoids the time-consuming solution of linear systems. Like in ISPH, a splitting scheme is employed. Firstly, the intermediate velocities and positions are predicted using non-pressure forces

The intermediate densityρ*can be obtained based on the predicted intermediate velocityand position. In the original work by Solenthaler and Pajarola[10], a mass summation approach is used

and their formulations of PCISPH are based on this density prediction. However, this mass summation approach leads to underestimated densities near boundaries and free-surfaces. We find that this approach leads to unstable behaviors in violent free-surface flows. Therefore, we base our density prediction and therefore the PCISPH formulations on the more popular continuity equation, which reads

In the predictive step, the incompressibility is not satisfied. The intermediate density deviatesfrom the reference density0ρby an error

To minimize the density error, a correction of the current pressure Δpiis computed as

whereKiis a pressure coefficient. This pressure coefficient plays a crucial role in PCISPH and its specific formulation will be derived in Section 2.2.

At the beginning of a computational step, the pressure is initially set as zero. For each iterative correction, the pressure is accumulated as

The pressure forces resulted from the corrected pressure field are computed according to Eq. (4). This pressure forces are then used to update the intermediate velocities, positions and densities. Because the pressure corrections are formulated to minimize density variation, incompressibility is achieved in iterations. This procedure is very similar to the Jacobi iteration used in solving linear systems.

The convergence criterion for PCISPH is

whereemaxis a tolerance of compressibility violation.Once convergence is reached, the solutions at timet+Δtare computed as follows:

Table 1 Comparison of computational times

The complete algorithm of PCISPH is given in Table 1. Note that in the predictive-corrective procedure, there is no need to solve any linear system. As in standard WCSPH, all the computations take place at the particle level. Therefore, PCISPH is suitable for massive parallel computing.

2.2 Derivation of pressure coefficient

In this section, we proceed to derive the pressure coefficientKi. The density changes during time span Δtcan be approximated using the continuity equation

where Δx=uΔtis used. To obtain the relation between pressure correction and density variation, we focus on the density change caused by the pressure forces. The particle position change Δxresulted from the pressure forceis computed as

However, the pressure corrections are unknown and need to be determined. Besides, the standard approach for pressure force computation at one particle requires pressure information from neighboring particles. This dependency between unknowns is undesirable for developing a local predictive-corrective procedure.Following Solenthaler and Pajarola[10], we assume that all particles in the neighborhood of particleihave the same unknown pressure correctionp~i, and all particles have the reference density0ρowing to incompressibility. Based on these two assumptions,we can compute the pressure force at particleias

Substituting Eq. (18) into Eq. (17), we have

Similarly, we can compute Δxj. Note that our purpose is to derive the relation between pressure correction and density variation at the central particlei. To this end, we only consider the position changecaused by the pressure force between particleiandj. Owing to force symmetry, the particlejgets the following force from the particlei

The corresponding position change of particlejis

With the position changes at particlei,j, the Eq.(16) becomes

where

Solving for the unknown pressure correctionp~i,we have

The meaning of Eq. (23) is that a pressure correctioncan result in a density variation of Δρiin the time duration Δtat particlei.

In the predictive step, we predict the intermediate densityand obtain the density errorAccording to Eq. (23), a pressure correctionresulting in a density change ofcan be applied to reverse this error

Note that several assumptions are made in the derivation of the pressure coefficientKi, therefore, a single corrective step cannot guarantee incompressibility at all particles. In PCISPH, multiple corrective iterations are used to ensure that the maximum density fluctuation is lower than a predefined tolerance. In the original PCISPH[10], the pressure coefficient is a constant calculated based on a prototype particle distribution to avoid instability caused by particle deficiency at free-surfaces and boundaries. In our modeling, the pressure coefficientKiis calculated independently for each particle, and fast convergence is observed.

2.3 Numerical diffusion and pressure smoothing

PCISPH is not a truly incompressible approach but an iterative pseudo-EOS method. From the derivations, it can be found that like WCSPH,PCISPH also has an EOS-like equation, i.e., Eq. (12).However, unlike WCSPH, it makes no use of a predefined stiffness. Instead, the stiffness is calculated automatically during the simulation. It is found that simulations using the original PCISPH exhibits spurious pressure oscillations which are even stronger than pressure oscillations observed in classical WCSPH. The reasons are twofold. First, it uses an EOS like WCSPH; therefore, it also suffers the so-called short-length-scale-noise[13]. Second, much larger time step is used in PCISPH simulations, which makes the pressure field even more oscillating. It may give nice animations in computer graphics, but its accuracy needs to be improved before it can be applied to CFD focusing on physics. In this work, we introduce two improvements for PCISPH to enhance its accuracy.

Recently, several numerical diffusion terms are proposed to improve the pressure accuracy in WCSPH[14]. These numerical diffusion terms are added to the continuity equation to reduce the numerical noise. Because the pressure field in WCSPH is calculated based on the density field, the diffusion terms can significantly reduce the oscillations thus give rise to smooth pressure field. Although PCISPH does not calculate the pressure field directly from the density field, a smooth density field can obviously help to reduce numerical noise. To this end, the numerical diffusion term proposed by Molteni and Colagrossi[14]is added into the continuity equation when predicting the densities. With the numerical diffusion term, the new equation for predicting the density can be written as

whereδis a parameter usually taken as 0.1. Like in the formulation of artificial viscosity, the speed of soundcsis a numerical coefficient which does not affect the time step of PCISPH.

To further reduce the numerical noises, the original pressure is replaced by a smoothed pressure when computing the pressure forces

whereχis a parameter taken as 0.5 in this work.This smoothed pressure field can avoid extremely large pressure forces caused by singular particles; thus,it helps to stabilize the computation. Numerical examples in Section 4 also show that it is also helpful in producing smooth pressure field and has negligible influence on the flow kinematics.

3. Numerical implementation

3.1 Generalized boundary condition

In PCISPH, boundary treatment is not only used to prevent particle penetration, but also related to convergence. An inappropriate boundary treatment method may result in strong density fluctuations near boundaries. Consequently, it may lead to slow convergence, small time steps or even unexpected breakdowns of computations. Generally, all boundary treatment methods used in WCSPH can be applied to PCISPH. In this work, the generalized boundary method proposed by Adami et al.[11]is employed with some modifications.

Fig. 1 Sketch of the generalized boundary treatment of solid walls

As shown in Fig. 1, in the presented method the solid boundaries are represented using boundary particles, which are either fixed or moved with prescribed velocities. For simplicity, the particles are placed on a regular lattice with particle spacing Δr,the same as the initial fluid particle spacing. The mass of the boundary particles is initially computed using the particle spacing Δrand the reference density0ρ.The boundary particles ensure that the kernel support domains near solids are complete, thus solve the so-called particle deficiency problem. As illustrated in Fig. 1, when performing SPH computations, the support domains of fluid particlefinclude neighboring boundary particles, which participate in PCISPH computations like normal fluid particles.However, the densities and pressures at boundary particles are not computed in the same predictivecorrective manner for fluid particles. Instead, they are extrapolated from the fluid field using the SPH kernel function

where subscriptwandfdenote boundary and fluid particles, respectively. Following Adami et al.[11],the pressure at boundary particle is extrapolated as

wheregis the gravity,xwfis the vector from boundary particle to fluid particle. The gravity correction is included in the pressure extrapolation to reproduce a correct pressure gradient. Note that in this work, the densities at boundary particles are directly obtained using extrapolation, not computed from pressure using an EOS as in Adami et al.[11].

Using Eqs. (27), (28), the updating of densities and pressures at boundary particles is performed once in each PCISPH iteration just after new densities and pressures are obtained in the fluid region.Consequently, the extrapolated densities and pressures are always consistent with the physical quantities in the fluid field. Both non-slip and free-slip boundary conditions can be considered using the boundary treatment method. For non-slip boundary condition,the viscous forces between wall and fluid particles are taken into consideration; while these viscous forces are neglected in the case of free-slip boundary condition. The employed boundary treatment method is accurate, giving rise to smooth density and pressure fields near boundaries. Furthermore, its numerical implementation is convenient as it only involves one additional extrapolation process. It can handle boundaries with arbitrarily complex geometries and motions without difficulty.

3.2 Dynamic time stepping

In SPH simulations, dynamic time stepping is usually employed to achieve higher efficiency. In WCSPH, time step is controlled by the CFL condition and viscous condition. For low viscosity fluids, a popular form of variable time step is[1]

whereλCFLis a constant called CFL number in the order of 0.1,aiindicates the acceleration at particlei, and the speed of soundcsis 10 times bigger than the maximum fluid velocity. Except for problems involving extremely violent impact, the dominant factor restricting the time step in WCSPH is the speed of soundcs.

The above variable time step only works for WCSPH. In PCISPH, spurious density and pressure oscillations can be observed if Eq. (29) is used,because the pressure coefficientK, hence the pressure, are highly dependent on the time step Δt.From Eq. (24), it is obvious that a smaller time step results in a larger pressure coefficientK, and vice versa. If the time step changes too much between computational steps, the convergence is hard to achieve. Ihmsen et al.[15]observed that if the difference in time step between two sequential steps is less than 0.2%, PCISPH converges and produces results without extreme density and pressure oscillations.

The variable time step algorithm used in this work is based on that proposed by Ihmsen et al.[15].The time step is automatically adjusted according to the maximum density error, the average density erro r(Nis the total particle number), the maximum velocityumax=and the maximum accelerationamax= Initially, the time step Δtiniis estimated as

whereUmaxis an estimated maximum velocity in the whole simulation. It can be found that the initial time step is quite large compared with the time step used in WCSPH, becauseUmaxis approximately 10 times smaller thancs. To ensure smooth change in simulations, the time step Δtis not obtained from the CFL conditions directly but calculated by slightly increasing or reducing the time step used in the last computational step. During simulations, the time step is increased by 0.2% if all the following conditions are satisfied

On the other hand, the time step should be reduced by 0.2% if any of the following conditions holds

The constants in Eqs. (31), (32) are similar to the empirical values used by Ihmsen et al.[15]. However,our choices are more conservative, resulting in smaller time step. Moreover, we set an upper limit Δtmaxto the time step

Due to the conservative choices of constants and the upper limit, generally our approach gives smaller time step compared with that in Ihmsen et al.[15].However, our approach has two benefits in fluid flow simulations focusing on physics. Firstly, the smaller time step gives rise to more accurate pressure and force results. More importantly, limiting the maximum time step can solve the difficulties caused by strong numerical shocks. In case of strong numerical shocks, the maximum density error can be so big that reducing the time step by 0.2% is insufficient if the current time step is already very large. Therefore, handling of strong numerical shocks is a challenging task in dynamic time stepping of PCISPH. Ihmsen et al.[15]proposed to use some criterion to detect strong numerical shocks. If a strong numerical shock is detected, their algorithm goes two computational steps backwards and restarts the simulation with smaller time step. Their approach is more difficult to implement and demands more memory. It does not satisfy the requirement of smooth time step changing, leading to large pressure oscillations. On the other hand, as the maximum time step is limited in our approach, convergence can always be achieved even in the case of strong numerical shocks.

Comparing Eqs. (29), (31) and (32), it can be found that the time step in PCISPH is of one order larger than that in WCSPH. With the proposed dynamic time step, PCISPH usually converges in a few iterations. Consequently, it is expected that PCISPH is much faster than WCSPH. The actual performance comparison between PCISPH and WCSPH will be discussed in Section 4.

4. Numerical examples and discussion

Three standard validation cases suggested by the SPH European Research Interest Community(SPHERIC, sphericsph.org) are simulated in this section. Experimental and numerical results of these cases are available; therefore, investigations on the performance of PCISPH and the comparison between PCISPH and WCSPH can be carried out. All the PCISPH and WCSPH simulations are performed using a GeForce GTX 1080Ti graphics card. The GPU implementation of PCISPH follows the algorithm given in Table 1. For other details on GPU-accelerated SPH, we refer the readers to Crespo et al.[16]. In all simulations, the Wendland kernel is used with a smoothing lengthh=1.35Δr, where Δris the initial particle spacing. The radius of the support domain is taken as 2h. Sinceδ-SPH is arguably the state-of-the-art WCSPH scheme, it is employed in all the WCSPH simulations withδ=0.1. The high-accuracyδ-SPH is used as a reference method for checking the performance of PCISPH. In all the PCISPH simulations, the speed of soundcs, which is only used in artificial viscosity and numerical diffusion, is taken as 30 m/s. In theδ-SPH simulations,the speed of sound is taken 10 times larger than the maximum velocity, giving rise to a density fluctuation of 1%. The convergence criterion of PCISPH is alsoemax=1%.

Fig. 2 The geometry and boundary conditions of the 3-D liddriven cavity flow

4.1 3-D lid-driven cavity

The 3-D lid-driven cavity problem consists of a cubic cavity filled with fluid which is covered by a horizontal moving lid. As shown in Fig. 2, the cubic cavity is of side lengthL, the covering lid is moving with a constant velocityuwin thexdirection. A case of Reynolds numberis simulated using PCISPH andδ-SPH with four resolutions, i.e.,L/Δr=25, 50, 75 and 100. In this paperL/Δr=1m is used, resulting to particle size of respectively Δr=0.04 m, 0.02 m, 0.0133 m and 0.01 m in the four simulations. In the simulations, the boundary particles in the lid move with a velocity ofuw=1m/s . For this Reynolds number, the flow is in the laminar regime, so the laminar viscosity treatment method in Eq. (4) is used. The total physical time is 100 s. Gravity is not considered. Despite the simple boundary conditions, the lid-driven cavity problem has no analytical solutions. Therefore, numerical results are compared with data obtained using a finite volume method (FVM) solver by Yang et al.[17].

Fig. 3 (Color online) The steady-state velocity fields on the middle plane /=0.5 y L

The steady-state velocity fields on the middle planey=L/2 obtained with resolution /=50LΔrare shown in Fig. 3 for both PCISPH andδ-SPH.Generally, the two velocity fields are quite similar,showing a central main vortex and two secondary vortices in the lower corners. This observation is consistent to other numerical results obtained using FVM[17]and SPH methods[9]. For further comparison,the velocitiesuxalongx=L/2 anduzalongz=L/2 on the planey=L/2 are shown in Figs. 4,5, respectively.

Overall, the velocity profiles obtained using PCISPH,δ-SPH and FVM have the same trend.However, the SPH results show lower peak velocities in both thexandzdirections. Although the liddriven cavity is simple in geometry and boundary condition, previous studies show that accurate SPH simulations of this problem requires much attention in boundary treatment and SPH formulations[9,18]. The current generalized boundary treatment method is convenient but does not impose exact non-slip boundary condition on solid walls, which might be the cause of the discrepancy between SPH and FVM results. Nevertheless, generally the results from PCISPH andδ-SPH simulations show high degree of agreement. It is observed that both PCISPH andδ-SPH show similar convergence behaviors as the initial particle spacing decreases. It can be concluded that in the modeling of the 3-D lid-driven cavity problem, the presented improved PCISPH scheme can achieve the same level of accuracy asδ-SPH.

Fig. 4 (Color online) Profiles of uz on the middle plane y=L/2 along x=L/2

Fig. 5 (Color online) Profiles of uz on the middle plane y=L/2 along z=L/2

Table 2 shows the comparison of numerical efficiency between PCISPH andδ-SPH. As expected, tt is found that PCISPH is much faster. The speed-up is about 5-6 compared withδ-SPH regardless of the particle number used in simulation.

Table 2 Comparison of efficiency between PCISPH and δ-SPH in the 3-D lid-driven cavity problem, where Δr is the initial particle spacing, N is the total particle number, FPS is the computational step per second of wall-clock time, Δt is the average time step and Tsim is the total wall-clock time used for simulation

4.2 Dam-breaking with an obstacle

In this section, we investigate PCISPH's capability of simulating violent free surface flow and impact, which are the main application field of SPH methods. A three-dimensional dam-breaking problem with a box obstacle on the bed of tank is considered,which was experimentally and numerically studied using a Eulerian volume of fluid (VOF) method by Kleefsman et al.[19]. The experimental set-up is shown in Fig. 6. A tank 3.22 m long and 1 m wide is used in the experiment. In the right part of the tank, a body of waterH=0.55 m deep is confined by a gate. After pulling up the gate, the water flows into the tank and has violent impact with the box obstacle and walls of the tank. Four height probesH1-H4and eight pressure sensorsP1-P4are installed in the tank to measure the flow.

Fig. 6 (Color onli ne) Experimental s et-up and measurement positionsforthedam-breakingproblem. H1-H4 arethe positions of height probes, and P1-P8 are the positions of pressure sensor. The blue region represents the water body initially confined by the gate

Both the improved PCISPH andδ-SPH are employed to model this problem. The results are also compared with numerical results data from Kleefsman et al.[19]. As the flow is highly turbulent, artificial viscosity given in Eq. (5) is used to stabilize computation. The artificial viscosity coefficient is taken as 0.05 for both PCISPH andδ-SPH.Simulations with four resolutions, i.e., Δr=H/30,H/45,H/60 andH/75, are performed, withHthe initial water depth. As a result, the resulted particle sizes areΔr=0.0183 m, 0.0122 m, 0.0092 m and 0.0073 m, respectively.

Fig. 7 (Color online) Free-surface profiles and pressure results obtained using three different SPH methods

Initially, the fluid is in hydrostatic state, where the pressure is computed from the height of the overburden. The total simulated physical time is 6 s.For post-processing, the height of the free surface at a certain position is determined using the free-surface detection algorithm proposed by Marrone et al.[20],while the pressure is obtained by extrapolating the pressure field of the fluid region.

Figure 7 shows the free-surface profiles and pressure distributions obtained using three different SPH schemes with computational resolution Δr=H/45. The three SPH schemes are the original PCISPH[10], the improved PCISPH presented in this paper andδ-SPH. It is observed the three SPH schemes give rise to similar free-surface profiles. The obtained profiles are in good agreement with the experimental observations[19], and numerical results obtained using grid-based methods[20]and SPH methods[21-22]. Violent impact between the flow and the obstacle is observed. Despite the high degree of similarity in the results, PCISPH gives rise to simulations with slightly more fragments and fluid drops, which is more realistic when compared to the test pictures. The spatial distribution of pressure obtained using the original PCISPH suffers spurious oscillations, especially in the area with much particle rearrangement. On the other hand, smooth pressure fields are obtained using the improved PCISPH andδ-SPH.

Fig. 9 (Color online) History of the free-surface height recorded at H2

The evolutions of free-surface height at different locations are shown in Figs. 8, 9. The free-surfaces obtained using PCISPH andδ-SPH agree well with the test data. The pressures at1P,3Pare shown in Figs. 10, 11, respectively. Generally, the numerical results agree well with the experimental data. The results from PCISPH show higher temporal fluctuations than those obtained usingδ-SPH. In PCISPH, pressure is computed by accumulating the pressure corrections, which is obtained using pressure coefficientKand density errorρerr,*. Between two adjacent steps, the pressure coefficientKvaries in response to the changes in time step and particle arrangement. Therefore, it is unsurprising to observe some temporal pressure oscillations in PCISPH results.Nevertheless, the temporal oscillations of pressure are within a low level. The accuracy of the pressure results is like that of ISPH, and better than classical WCSPH without the new improvements such as density diffusion.

The comparison of efficiency between PCISPH andδ-SPH is given in Table 3. The speed-up of PCISPH overδ-SPH is between 4.57 and 6.45. The speed-up in this dam-breaking case is slightly lower than that in the lid-driven cavity case. This reason is that the violent impact makes the convergence more difficult to achieve, which increase the computational cost at each step and reduce the time step. This can be confirmed by comparing the average time step in the two cases. In the lid-driven cavity problem, the average time step of PCISPH is more than 10 times bigger than that ofδ-SPH, while in the dambreaking case, the ratio between time steps of PCISPH andδ-SPH is less than 10. Given the fact that PCISPH andδ-SPH has similar accuracy, PCISPH proves to be an appealing choice.

Table 3 Comparison of efficiency between PCISPH and δ-SPH in the 3-D dam-breaking problem

Fig. 10 (Color online) Pressure evolution recorded at 1P

Fig. 11 (Color online) Pressure evolution recorded at 3P

Fig. 12 (Color online) Set-up and measurement positions for the horizontal sloshing problem

4.3 Horizontal sloshing

The problem of a partially filled tank subjected to horizontal movement is simulated. The horizontal sloshing involves free-surface flow and moving boundary. Experimental results and analytical solutions of this problem can be found in Faltinsen et al.[23]. In this work, three-dimensional simulations using PCISPH andδ-SPH are performed. Figure 12 shows the geometry of the simulated problem, which is identical to the experimental set-up[23]. The size of the tank is 1.73 m in length, 0.2 m in width and 1.0 m in height.The filled water height isH=0.6 m . The tank moves following a functionS=Acos(2πt/T), whereSis the horizontal displacement of the tank,A=0.032 m ,T=1.5s are amplitude and period of the external excitation. Eight simulations using PCISPH andδ-SPH with varying resolutions are carried out. Artificial viscosity is employed in both PCISPH andδ-SPH w ith v iscosit y co efficientα=0.05. Thechangesofwaveheightandpressure are measured atH1,P1, respectively.

Figure 13 shows the fluid flow and pressure distribution obtained using PCISPH with Δr=H/90.It is found that the spatial distribution of pressure is very smooth. The employment of numerical diffusion and pressure smoothing greatly enhances the accuracy of pressure and prevent spatial fluctuations.

Fig. 14 (Color online) Change of wave height recorded at H1

Fig. 15 (Color online) Change of pressure recorded at 1P

Figure 14 shows the wave height atH1. Both PCISPH andδ-SPH results are well corroborated by the analytical solutions. The PCISPH results slightly underestimate the amplitude, while theδ-SPH simulations overestimate the amplitude. Generally, the simulations show that PCISPH yields slightly better results thanδ-SPH in terms of the flow motion.Note that the peaks given byδ-SPH are more and more delayed, indicating excessive energy dissipation.This phenomenon is also observed in PCISPH results,but to a lesser extent.

Table 4 Comparison of efficiency between PCISPH and δ-SPH in the horizontal sloshing problem

The pressure history atP1is shown in Fig. 15.The pressure changes smoothly in the results obtained usingδ-SPH regardless of the resolution. As for PCISPH, the pressure history shows more temporal oscillations, especially for simulations with low resolutions Δr=H/30,H/60. The reason for this temporal pressure oscillations in PCISPH simulations are discussed in Section 4.2. Again, the pressure oscillations are within an acceptable level. With high resolution, smooth temporal pressure can be achieved.The comparison of efficiency is shown in Table 4.Again, around five times of speed-up is observed.

Fig. 16 (Color online) Summary of the speed-up in the threenumerical examples with different particle numbers

4.4 Discussion on the efficiency of PCISPH

Through three examples we investigate the performance of the improved PCISPH. It is found that the speed-up of PCISPH overδ-SPH is from 4.23 to 6.77. Except the influence of implementation, the speed of PCISPH depends on many factors, including the convergence criterion, the dynamic time-stepping parameters and the characteristics of the problem,among others. On the other hand, th e e fficienc y ofδ-SPH also varies according to speedofsound,CFL number and many other factors. Therefore, it is very difficult to give a strict and universal speed-up value when comparing PCISPH andδ-SPH.

Nevertheless, owing to the large time step allowed in simulations, PCISPH is undoubtedly much faster than WCSPH. Moreover, like WCSPH, all computations in PCISPH are in the particle level, the solution of a linear system is not required. Therefore,PCISPH is very suitable for parallel computing. This can be confirmed by Fig. 16, where the speed-up in all the three numerical simulations are plotted against the particle number. The speed-up oscillates with the particle number, no clear trend of decreasing efficiency is observed. Even with millions of particles,PCISPH is still 4-5 times faster thanδ-SPH. This is a desirable feature, as in 3-D simulations, the particle number increases very fast with higher resolution and larger computational domain.

It is worth noting that the speed-up of PCISPH over WCSPH is owing to that the time step in WCSPH is restricted by the speed of sound. However,from the CFL condition, we can find that there are other factors determining the time step. In cases where the speed of sound is not the dominant factor, PCISPH may perform poorer than WCSPH. Examples of such cases are flows with very high viscosity. In these cases, the time step is controlled by the viscous effect.As a result, PCISPH may be slower due to the iterative process. Nevertheless, as shown in this work,in common CFD problems PCISPH bears much higher efficiency than WCSPH.

5. Conclusions

PCISPH is the SPH variant originally developed for computer graphics. It features high numerical efficiency compared with the traditional SPH methods such as WCSPH and ISPH. In this paper, we employed PCISPH in physics-focusing fluid flow modelling for the first time. An improved PCISPH method was presented, where the density prediction is based on the continuity equation instead of the mass summation used in the original PCISPH. Several enhancements, including numerical diffusion, pressure smoothing, generalized solid boundary condition and dynamic time stepping algorithm, were introduced to improve the accuracy and stability of this method.

The performance of PCISPH in modelling physics-focusing fluid flows was investigated. Three standard SPH validation cases were simulated, the results were compared to experimental data, analytical solutions and numerical results obtained using other methods. Comparisons between the improved PCISPH and the state-of-the-artδ-SPH were presented. It is found that the accuracy of the improved PCISPH is almost of the same level asδ-SPH. It provides simulation results with accurate flow motion and smooth pressure field. However, PCISPH gives more temporal pressure oscillations thanδ-SPH, and till now this is the only weakness we found in PCISPH when comparing it toδ-SPH. This temporal pressure oscillations are caused by the unique way of calculating pressure in the predictivecorrective procedure. How to alleviate this temporal pressure oscillations still requires further study.

The most appealing feature of PCISPH is its high efficiency. In our three numerical examples, the speed-up of PCISPH overδ-SPH is in the range of 4.23-6.77 despite that very conservative choices of the time stepping parameters are used. Like in WCSPH,PCISPH computations are at the particle level,solutions of global linear system are not required. This means that in simulations using a very large number of particles its advantage of high efficiency still holds.This fact renders PCISPH particularly attractive to massive parallel computing. It has great potential in large-scale three-dimensional modelling of fluid dynamics problems.