An Improved CSPM Approach for Accurate Second-Derivative Approximations with SPH

We compare several approximations for second derivatives with Smoothed Particle Hydrodynamics (SPH). A first-order consistent approximation, derived from the zeroth-order consistent Corrective Smoothed Particle Method (CSPM), is proposed. The accuracy of the new method (ICSPM) is similar to that of the Finite Particle Method (FPM) and Modified Smoothed Particle Hydrodynamics (MSPH), but it is computationally less expensive. We demonstrate the accuracy of our method by studying heat conduction in a slab with discontinuous conductivity coefficients. We use both uniformly and pseudo-randomly distributed particles.


Introduction
Smoothed Particle Hydrodynamics (SPH) is a spatial discretisation method to solve partial differential equations.It is a mesh-free, Lagrangian method in which the system is represented by a finite set of particles.SPH was originally developed to solve astrophysical problems [1] [2], but many advantages of the method compared to conventional methods for specific types of problems made it attract attention in other areas, such as fluid and solid mechanics.For instance, with SPH, free surfaces are very easy to deal with [3].Also, SPH is quite straightforwardly applicable to multiphase flows [4] [5] [6], flows with fluidstructure interaction [7] [8], and so on.
As opposed to the astrophysical problems that SPH was initially applied to, most fluid and solid mechanics problems have solid/physical domain boundaries.
As a consequence, the support domain of particles close to these boundaries overlaps with the empty area behind it, which has a major negative influence on the accuracy of the SPH approximations.Consistency is also lost if particles are not uniformly distributed, which is the case in most simulations.Restoring the consistency for functions and first-derivative approximations is fairly straightforward and can be achieved by using correction terms as in, e.g.[9].An extensive comparison between correction methods can be found in [10].Correcting second-derivative approximations, however, is more complicated.
Second derivatives (or the Laplacian) appear in viscous terms, conduction equations and in the pressure Poisson equation that is used in incompressible SPH.
Flebbe et al. [11] straightforwardly computed the second-derivative terms by first computing the gradient of the unknown variable and then computing the divergence of the result.This sometimes led to non-physical oscillations in the solution [12].Originally, approximations for second derivatives included second derivatives of the kernel function, which are sensitive to particle disorder [13] [14].Therefore, several ways to restore the consistency have been proposed in the literature.Many researchers have suggested to use expressions based on first derivatives of the kernel function rather than the second derivatives, because they are less sensitive to the particle distribution [13] [14].However, these expressions do not solve the consistency problem near boundaries.Some researchers have therefore proposed to include boundary terms in the approximation [15] [16].Also, there is quite some literature available on Reproducing Kernel Particle Methods (RKPM) and similar approaches [17] [18] [19], which are designed to give approximations up to a desired order of consistency.They can lead, however, to partially negative, non-symmetric and non-monotonically decreasing smoothing functions and are therefore not preferred for hydro-dynamic simulations [20].A different approach was followed by Chen and Beraun [21], who present a method in which an estimate for the leading error term is subtracted from the approximation.Other researchers use methods in which a function and all of its desired derivatives are solved simultaneously [22] [23] [24].These methods have higher accuracy, but are also computationally expensive.
In this work we only consider methods that use conventional smoothing functions.This excludes the previously mentioned RKPM.We describe the original SPH method to approximate second derivatives, as well as the methods by Chen and Beraun [21] and Zhang and Batra [23].We also propose an improvement to the estimate by Chen and Beraun.Our method is computationally only slightly more expensive than the method by Chen and Beraun, but its accuracy is similar to that of the method by Zhang and Batra.We compare numerical results in one and two dimensions for both uniformly and non-uniformly distributed particles, but are especially interested in the latter, because it most accurately resembles particle distributions from actual simulations.Initial results have been presented in [25] [26].

Smoothed Particle Hydrodynamics
There are two approaches to derive the SPH equations for fluid flows in the lite-rature.The first one is primarily concerned with choosing a density estimate.Substituting this estimate into the Lagrangian then leads to a system of equations conserving both linear and angular momentum.This approach is explained in detail by Price [14].The second approach considers SPH to be a spatial discretisation scheme which can be applied to any set of equations, e.g. the Navier-Stokes equations.In this paper we adopt the latter view on SPH.A short derivation of the SPH equations will follow.For more details we refer to, e.g.[27].

The Kernel Approximation
Consider a function f , defined on a domain Ω .Then the value of f at an arbitrary point x ∈ Ω can be written as: Here, ( ) δ ⋅ is the Dirac delta distribution: In SPH, ( ) δ ⋅ is replaced by a continuous function ( ) smoothing parameter.K is called the smoothing or kernel function and should converge to the Dirac delta distribution as h goes to zero.Preferably it is also radially symmetric, has compact support and satisfies the unity condition for all x ∈ Ω : ( ) Not satisfying the unity condition-to be more precise, its discrete equivalent-is the main reason for inaccurate approximations.This will be explained in more detail later.Replacing the Dirac delta function by K leads to the fol- lowing approximation: which is called the kernel approximation of f .Note that it is a weighted average over the values of f at all points x ∈ Ω  , including x itself.

The Particle Approximation
To get a numerically useful expression, a quadrature rule is applied to (4).This results in a partitioning of the (computational) domain Ω by a finite number of so-called particles.The summation for a particular particle is limited to those particles that are within the support domain of the particle.This is illustrated in Figure 1.Here, H is the support radius.Defining i S as the set of particles within the support domain of particle i , (4) is approximated by: where j V is the volume of particle j , j f denotes ( ) . The approximation in (5) is called the particle approxima- tion of f .When applied to a physical problem, the volume of a particle is usually replaced by the ratio of its mass and density, j j j V m ρ = , but since this is not necessary for our derivations we will stick to the notation in (5).

Derivatives in SPH
Approximations for derivatives are found by replacing f by its derivatives in (4).For instance, substituting f ∇ leads to: Integrating by parts gives: ( where B is a function depending on f and K and ∂Ω denotes the boundary of Ω .The first term on the right-hand side of ( 7) is a boundary term that is usually neglected, since it is zero for particles sufficiently far away from the boundaries.Hence, only the second term on the right-hand side remains and, after discretising, we find: : .
Similarly, we can derive an approximation for the Laplacian by substituting 2 f ∇ for f in (4).Integrating by parts twice, neglecting boundary terms and discretising leads to: : .
Note that in both the approximation for the gradient and the Laplacian the derivative operator has switched from f to K .Since K is known, these ex- pressions can be computed.Substituting a Taylor series expansion for j f around i x in (9) reveals that the leading error term of this estimate can be subtracted to find the more accurate estimate:
Here, the subscript "SPH" was added to indicate that this is the standard SPH equation used in Section 5.In practice, ( 10) is rarely used.The reason is that for most kernels the second derivatives have very steep gradients, making the ap-proximation very sensitive to irregularities in the particle distribution [13] [14].An alternative expression is given by Brookshaw [13]: ( ) where − .This approximation uses first derivatives of the kernel func- tion instead of second derivatives, which makes it less sensitive to changes in the particle distribution [13].However, it has accuracy issues that are explained in the next section.

Error Analysis
The simple expressions in (10) and ( 11) have a downside: their accuracy decreases for particles close to boundaries or when the particle distribution is irregular.This can be shown by substituting a Taylor series expansion for j f around i x .For (10), for instance, this gives: where f H denotes the Hessian matrix of f .The term that we are approxi- mating, the Laplacian, is contained in the second term on the right-hand side of (12).The first term on the right-hand side of ( 12) is an ( ) -error term that vanishes in the ideal conditions of a uniform particle distribution and no boundaries, but it does not if a particle is close to a boundary or is surrounded by irregularly located particles.A similar ( ) -error term can be found for (11).Furthermore, it is impossible to distill the exact Laplacian from the second term, because the separate second-derivative terms have different coefficients.This leads to an additional ( ) 1  -error, which also holds for (11).In the next section we discuss several methods that were designed to give more accurate approximations.

Second Derivative Approximations with Higher Accuracy
This chapter describes two methods that approximate the Laplacian with higher accuracy than the conventional estimate in (11).We discuss the Corrective Smoothed Particle Method (CSPM) and the Modified Smoothed Particle Hydrodynamics (MSPH) method.

Corrective Smoothed Particle Method
This method, described in various papers by Chen and his colleagues [21] [28] [29] [30] [31] starts from the Taylor series expansion of j f around i x : ( ) In all derivations to come, we assume Ω is a two-dimensional domain and x and y refer to the two independent spatial directions.We stress, however, that it is possible to extend everything in Sections 3 and 4 to three dimensions.To compute second derivatives, the first term on the right-hand side of ( 13) is first subtracted from both sides of the equation.The result is multiplied by the vector ij K h , with h defined as: Multiplying ( 13) by j V and summing over all the particles leads to: In (15) and in the rest of this paper, the derivatives in ij K h are taken with re- spect to j x .The first term on the right-hand side of ( 15) is a ( ) -term.We wish to subtract it from both sides of the equation, but since i f ∇ is unknown we are forced to use an approximation for that as well.A derivation similar to the one that led to ( 15), but with where is the normalisation matrix defined by: T , : .
Multiplying ( 16) by the inverse of leads to a first-order accurate approximation for the gradient: .
The gradient approximation in ( 19) is substituted for i f ∇ in (15).Subtracting the first term on the right-hand side gives: It can be shown that the first term on the right-hand side of ( 20) can be written as , where Γ is the normalisation matrix defined by: with T ij ξ the following vector: T : , 2 , . .
By using an approximation for i f ∇ in (15), an error was introduced, as indi- cated by the approximation symbol in (20).As a result, normalising the first term on the right-hand side of (20) does not lead to the desired first-order accurate approximation.Instead, the first-order error of the gradient estimate makes that the approximation in ( 23) is zeroth-order accurate, i.e.: In Section 4 we show how this approximation can be improved to be first-order accurate.

Modified Smoothed Particle Hydrodynamics
Zhang and Batra [23] [32] describe a different approach to compute second derivatives.In their method, a vector i ϕ is computed that consists of an approx- imation of f itself, all its first-order derivatives and all its second-order deriv- atives: ( ) This is achieved by multiplying (13 h , leading to six equ- ations for six unknowns (in two dimensions) for each particle.Because all unknowns are put in one single vector, this method requires the inversion of 6 6 ×matrices.This is computationally expensive, but it leads to more accurate results than with the method described in Section 3.1.In fact, if we isolate i f h from i ϕ , we find: It is possible to achieve similar accuracy without using these larger matrices.For that we need only one adjustment to CSPM.This is described in Section 4.

Improving the Normalisation Matrix
In Section 3.1 we showed that the CSPM approximation for second derivatives is zeroth-order consistent.The crucial step in improving the accuracy is keeping the first-order terms in ( 16) separate from the second and higher-order terms.This gives: so that, instead of (19), we find: Accounting for the extra term in (28) leads to the adjustment of (20) to: with , i h Γ and ij ξ as defined in ( 21) and ( 22), respectively.Note that the nor- malisation matrix , i h Γ is the normalisation matrix used in CSPM.We propose to use the following approximation for second derivatives instead of ( 23): with , i h Γ  the normalisation matrix: , which we designate as improved CSPM (ICSPM).This method is slightly more expensive than CSPM, but it requires the same effort regarding matrix inver- sions.Yet, it follows directly from the definition of , i h Γ  and ( 29) that this ap- proximation is first-order accurate: In Section 5 we explore the behaviour of ICSPM and compare it with the standard SPH, CSPM and MSPH approximations described in Sections 3.1, 3.2 and 2.3.

Accuracy versus Computational Cost
The methods described in Sections 3.1, 3.2 and 4 are more accurate, but also computationally more expensive than (11).Therefore, if the available computation time is limited, (11) is the approximation to go for.Also, if the domain has no boundaries, particles are distributed regularly, or if accuracy is simply not of great importance, ( 11) is a perfectly good estimate for the Laplacian.However, if accuracy is important and a somewhat higher computational cost is not an issue, CSPM, MSPH and ICSPM are attractive alternatives.Moreover, these methods give approximations for all second derivative-terms instead of just the Laplacian.
The MSPH method by Zhang and Batra [23] [32] requires the inversion of 6 × 6matrices (in two dimensions) and 10 × 10-matrices (in three dimensions).This makes the method computationally more expensive than CSPM, which splits the problem into smaller sets of equations, so that it only requires the inversion of 2 × 2 and 3 × 3-matrices (in two dimensions) and 3 × 3 and 6 × 6-matrices (in three dimensions).MSPH, however, is more accurate, as it is first-order accurate (26), whereas CSPM is only zeroth-order accurate (24).
With the improved normalisation step in Section 4, we have found a method that is computationally similar to CSPM-it only additionally requires the expli-cit computation of the sums in (31)-but which has the same order of accuracy as MSPH.We will verify this with numerical examples in Section 5.
In principle it is possible to increase the accuracy to any order.This might not be obvious with CSPM, but it is quite straightforward with MSPH.However, every increment in order introduces a number of extra equations to solve (simultaneously), so that it is not worth the effort for any order higher than one.Hence, if accurate approximations for second derivatives are required, we recommend ICSPM (30).This method is more accurate than CSPM, while the extra computational effort is negligible.MSPH has a similar accuracy, but is computationally more expensive.

Numerical Experiments
This section describes three numerical case studies to compare the various methods that were discussed previously.The first one yields the computation of a one-dimensional second derivative.In the second and third case study we consider a time-dependent problem on a two-dimensional domain.In all three experiments we use the Wendland 2 C kernel, which, in one-dimension, is given by: ( ) ( ) ( ) where ( ) ( ) and H is the radius of the support domain.As sug- gested by Dehnen and Aly [33], we choose the smoothing length equal to twice the standard deviation, so that ( ) . In two dimensions, the 2 C Wendland kernel has a different form: with ( ) In our test cases we use these kernels with R the ratio between particle distance and H .

Second Derivative on a One-Dimensional Domain
We start with a one-dimensional domain , on which we compute the second derivative of the given function ( ) ( ) . This test case was also performed by Zhang and Batra [23].The analytical solution is available and is equal to . Therefore, we can calculate the exact error, for which we use the infinity norm: where N is the total number of particles.We use 1.5 h x = ∆ , with x ∆ the (uniform) inter-particle distance, equal to that of Zhang and Batra [23].The two left panels of Figure 2 show the results with  We also calculate the errors according to (35) for several values of N .These are shown in the two right panels of Figure 2.They clearly show the ( ) -convergence of MSPH and ICSPM.Note that the error plots for MSPH and ICSPM coincide.

Heat Conduction on a Two-Dimensional Domain
Next, we compute the time-dependent temperature distribution on a two-dimensional domain , a test case by Cleary and Monaghan [34].The equation governing the conduction process is:  analytical solution of this problem is known and reads: , , sin π sin π , where Since our main interest is the comparison between the various spatial discretisations, it is sufficient to use the explicit Euler time stepping scheme.The time step used depends on the spatial discretisation distance and is chosen as In this case study we use two expressions for the error.The first one is similar to (35), but looks at the temperature itself rather than the Laplacian: The second one considers the relative errors at the positions of the particles instead of the absolute errors: We choose the smoothing length to be 1.2 xy h = ∆ , equal to that of Cleary and Monaghan [34].We use pseudo-random particle distributions, so as to compare the methods with distributions that resemble realistic particle configurations.
Starting with a uniform distribution, every particle is randomly shifted in both the horizontal and vertical direction, either to the left or to the right and up or down, with a maximum displacement of 2 5 xy ∆ in either direction.Boundary particles are only shifted in one direction, such that they stay on the boundaries, while corner particles are entirely excluded from the shifting process.
The results in the left panel of Figure 4 show the temperature distributions at larger N the CSPM method is less accurate than both MSPH and ICSPM.
The results in the right panel in Figure 5, for which we used the error in (40), show even more clearly that there is a significant difference in behaviour between CSPM and ICSPM (and MSPH).Clearly, the zeroth-order consistency of CSPM results in a zeroth-order convergence for the temperature as well, in contrast to ICSPM and MSPH, which are second-order convergent.MSPH is slightly more stable and accurate than ICSPM.

Heat Conduction in Slabs with Discontinuous Parameters
In the final test case we study the effect of discontinuous parameters.Cleary and Monaghan [34] describe the heat conduction in a slab of unit width which is periodic in the y -direction.The slab consists of two different materials, each with their own set of parameters, touching each other at the interface at .
When the temperature variation of the outermost points is small, the analytical solution can be approximated by [34] [35]: We start with a discontinuity only in the conductivity.For the left half it is set to 10  , while at the right 1 t = .As time progresses, the errors become smaller.This is confirmed by the errors in Table 1, which shows the maximum relative errors (40) for the Brookshaw method used Table 1.Maximum relative errors for the particles for which  2. The Brookshaw solution decreases only marginally when going from 40 x N = to 80 x N = .For the case 80 x N = , the error in the left table is one order lower for the ICSPM for these particular parameter choices.This behaviour seems to be coincidental, however, since the Brookshaw approximation and the ICSPM usually perform similarly in the infinite slab test case.
Finally, we also take the density to be discontinuous at the interface.We  2. Just like in the previous cases, the temperature distribution is captured very well and the errors become smaller for increasing numbers of particles.Since the particles are uniformly distributed and there are no boundaries causing problems, both methods behave similarly.The test cases in this section show, however, that ICSPM can also handle discontinuous parameters.

Summary
We studied approximations of second-derivative terms (or the Laplacian) in Smoothed Particle Hydrodynamics (SPH).Second derivatives appear in viscosity and conduction terms, and in Poisson equations.Traditionally, second derivatives have been approximated with the second derivatives of the kernel function, but these showed to be sensitive to particle disorder.Therefore various other methods and improvements have been suggested in the literature.
We proposed an improvement to the Corrective Smoothed Particle Method (CSPM) that increased the consistency of the estimates from zeroth-order to first-order.This method, called improved CSPM (ICSPM), was-together with several other methods-applied to one and two-dimensional test cases.The onedimensional test case clearly showed the high accuracy of ICSPM compared to CSPM.In the first two-dimensional test case we solved the heat equation and compared errors of the temperature itself, rather than the Laplacian.With the absolute error, the difference between CSPM and ICSPM was already clearly visible, but when we considered the relative error the difference between the two methods was even more obvious.In the third test case we computed the temperature distribution for an infinite slab consisting of two different materials.The results showed that our improved approximation is perfectly capable of handling discontinuous parameters.
ICSPM is computationally only slightly more expensive than CSPM.Yet, its accuracy is similar to that of more expensive methods, such as the Finite Particle Method (FMP) and Modified Smoothed Particle Hydrodynamics (MSPH), although MSPH proved to be a bit more stable in specific cases.

Figure 1 .
Figure 1.Support domain of a particle illustrated in two dimensions.
the left-hand side of (20) by the inverse of , i h Γ gives the final CSPM approximation for second derivatives: particles are distributed both uniformly (top) and pseudo-randomly (bottom), where in the latter case the boundary particles are kept on the boundary.As we can see in the figure, the particle distribution has no great effect on the results, except for the standard SPH discretisation.
sin π sin π , T x y x y = (37) which is shown in Figure 3.We assume four isothermal walls with 0 T = .The
with ref V the reference vo- lume, equal to the volume of an interior particle.

Figure 4 .Figure 5 .
Figure 4.The temperature profile at 0.5 t = for the particles for which x y ≈ and the particle distribution showing the volumes of particles, which form a Voronoi tessellation.(a) Temperature profile; (b) Particle distribution.
itially, the left half of the slab has zero temperature, 0 T =  , while the right half has temperature 1 r T = .Since κ is no longer constant on the whole domain, in the approximations it is replaced with: 2 : both halves of the slab.We uniformly distribute 40 particles in the x - direction and 20 in the y -direction and we apply an explicit Euler time stepping scheme with time step size 2 t x ∆ =∆ .
Figure 6 shows the temperature dis- tribution along the x -direction for the ICSPM derived in this work.The left panel shows the temperature at 0.2 t =

Figure 6 .
Figure 6.Heat distribution along the x -direction.The conductivity of the left half of the slab is 10 κ = 

5 t
Monaghan and the ICSPM.To avoid dividing by small numbers we only consider particles for which0.4x > .Next, we consider the case where both the conductivity and the specific heat are discontinuous.We choose 1 c = is shown in the left panel of Figure7.Again, we find that the temperature distribution is captured very well.The maximum relative errors are shown on the left in Table

2 t
= is shown in the right panel in Figure 7.The relative errors are shown on the right in Table

Figure 7 .
Figure 7. Heat distribution along the x -direction (a) at 5 t = when 3 r r c κ = =

Table 2 .
Maximum relative errors for particles of particles.(a) Errors; (b) Errors.