A new family of second order convergent weakly-compressible SPH schemes

by   Pawan Negi, et al.
IIT Bombay

Despite the many advances in the use of weakly-compressible smoothed particle hydrodynamics (SPH) for the simulation of incompressible fluid flow, it is still challenging to obtain second-order convergence numerically. In this paper we perform a systematic numerical study of convergence and accuracy of kernel-based approximation, discretization operators, and weakly-compressible SPH (WCSPH) schemes. We explore the origins of the errors and issues preventing second-order convergence. Based on the study, we propose several new variations of the basic WCSPH scheme that are all second-order accurate. Additionally, we investigate the linear and angular momentum conservation property of the WCSPH schemes. Our results show that one may construct accurate WCSPH schemes that demonstrate second-order convergence through a judicious choice of kernel, smoothing length, and discretization operators in the discretization of the governing equations.


page 40

page 41


Implicit and semi-implicit second-order time stepping methods for the Richards equation

This study concerns numerical methods for efficiently solving the Richar...

Efficient second-order semi-implicit finite element method for fourth-order nonlinear diffusion equations

We focus here on a class of fourth-order parabolic equations that can be...

Second Order Operators in the NASA Astrophysics Data System

Second Order Operators (SOOs) are database functions which form secondar...

A Consistent Quasi-Second Order Staggered Scheme for the Two-Dimensional Shallow Water Equations

A quasi-second order scheme is developed to obtain approximate solutions...

A Deep Dive into the Distribution Function: Understanding Phase Space Dynamics with Continuum Vlasov-Maxwell Simulations

In collisionless and weakly collisional plasmas, the particle distributi...

Convergence of the likelihood ratio method for linear response of non-equilibrium stationary states

We consider numerical schemes for computing the linear response of stead...

Frame Field Operators

Differential operators are widely used in geometry processing for proble...

1 Introduction

Smoothed Particle Hydrodynamics has been used to simulate weakly-compressible fluids since the pioneering work of Monaghan (1994). Many variations of the basic method have been proposed to create an entire class of weakly-compressible SPH schemes (WCSPH). One particularly difficult challenge has been the poor convergence displayed by the WCSPH methods making it one of the SPH grand-challenge problems Vacondio et al. (2020).

The SPH method works by using a smoothing kernel to approximate a function wherein the choice of the kernel influences the accuracy of the method. The length scale of the smoothing kernel is often termed the support radius or smoothing length, . A variety of kernels are used in the literature and the smoothing length may either be fixed in space/time or varying. One can show that for a symmetric kernel, the SPH kernel approximation is spatially second-order accurate in . However, the particle discretization of this approximation seldom achieves this and even first-order convergence often requires care and tuning of the smoothing length. Hernquist and Katz (1989) proposed that the support radius, be increased such that in three dimensions where is the local inter-particle separation. Subsequently, Quinlan et al. (2006)

derived error estimates for the standard SPH discretization and found that the ratio

must increase as the value is reduced to attain convergence; this is because of error terms of the form , where is a measure of the smoothness of the kernel at the edge of its support. This is an issue because as increases, the number of neighbors for each particle increases resulting in a prohibitive increase of computational effort. Furthermore, increasing the smoothing radius also reduces the accuracy of the method. This is the approach used in the work of Zhu et al. (2015) who proposed that the number of neighbors , where is the number of particles, in order to get convergence using SPH kernels.

Kiara et al. (2013a, b) shows that when the particles are distributed uniformly it is possible to obtain second-order convergence. The results of Quinlan et al. (2006) show that when using sufficiently smooth kernels (where is large or infinite), one can obtain second-order convergence. Indeed, Lind and Stansby (2016) demonstrate that for particle distributions on a Cartesian mesh one can obtain higher order convergence using higher order kernels.

However, for kernels that are normally used in SPH, the SPH approximations of derivatives become inaccurate even on a uniform grid unless a very large smoothing radius is used. Many methods have been proposed to correct the gradient approximation Bonet and Lok (1999); Liu and Liu (2006); Rosswog (2015); Huang et al. (2019). These typically ensure that the derivative approximation of a linear function is exact. This linear consistency is achieved by inverting a small matrix for each particle and using this to correct the computed gradients. This makes the derivative approximation second-order accurate but increases the computational cost of the gradient computation two-fold.

In the context of incompressible fluid flows, the governing equations involve the divergence, gradient, and Laplacian operators. These operators must be discretized and used in the context of Lagrangian particles. The divergence operator is encountered in the continuity equation and the discretization proposed by Monaghan (1994) is widely used. Rather than using a continuity equation, some authors Adami et al. (2013) prefer to use the summation density formulation proposed in Monaghan and Gingold (1983) to directly evaluate density. The gradient operator is encountered in the momentum equation. Many authors prefer using a discretized form that manifestly preserves linear momentum and as a result employ the symmetric form of the gradient operator Monaghan (1994); Bonet and Lok (1999). The symmetrization can be done in two different ways and Violeau (2012) shows that the selection of one form dictates the form to be used for divergence discretization in order to conserve volume (energy) in phase space.

Given the inaccuracy of the SPH approximation in computing derivatives accurately, the kernel corrections of Bonet and Lok (1999); Liu and Liu (2006) maybe applied to obtain linear consistency of the gradient and divergence operators. Unfortunately, the use of the corrections implies that linear momentum is no longer manifestly conserved. Frontiere et al. (2017) propose a symmetrization of the corrected kernel as originally suggested by Dilts (1999, 2000) to conserve linear momentum but the symmetrization implies that the operator is not first order consistent. Thus, the unfortunate consequence of demanding linear consistency is lack of conservation and vice-versa. A conservative and linear consistent gradient operator is currently not available.

The Laplacian is a challenging operator in the context of SPH. The simplest method is the one where the double derivative of the kernel is employed. However, the double derivatives of the kernel are very sensitive to any particle disorder. Chen and Beraun (2000) propose an approach by considering the inner product with each of the double derivatives and taking into account the leading order error terms. Zhang and Batra (2004) propose using the inner product with all the derivatives of the kernel lower and equal to the required derivative. This generates a system of 10 equations in two-dimensions. Korzilius et al. (2017) propose an improvement over the method of Chen and Beraun (2000) to evaluate the correction term. All of these methods require the computation of higher order kernel derivatives. Many authors Schwaiger (2008); Macià et al. (2012) proposed methods to correct the Laplacian near the boundary. In all of these formulations linear momentum is not manifestly conserved.

The Laplacian may also be discretized using the first derivative of the kernel using an integral approximation of the Laplacian. This was first suggested by Brookshaw (1985) and has been improved by Morris et al. (1997); Cleary and Monaghan (1999). They employ a finite difference approximation to evaluate the first order derivative and then convolve this with the kernel derivative. This formulation was structured such that it conserves linear momentum. However, these approximations do not converge as the resolution increases especially in the context of irregular particle distributions. Fatehi and Manzari (2011) propose an improved formulation by accounting for the leading error term; this makes the method accurate and convergent but makes the approximations non-conservative.

Another method to discretize the Laplacian is the repeated use of a first derivative and this has been used by Bonet and Lok (1999), and Nugent and Posch (2000). The formulation is generally not popular since it shows high frequency numerical oscillations when the initial condition is discontinuous. Recently, Biriukov and Price (2019) show that these oscillations can be removed by employing smoothing near the discontinuity.

Various SPH schemes have been proposed that use the above methods for discretization of the different operators. The simplest of the schemes is the original weakly compressible SPH (WCSPH) method Monaghan (1994); Gomez-Gesteria et al. (2010). This method is devised such that it conserves linear momentum as well as the Hamiltonian of the system. However, as the particles move they become highly disorganized and this significantly reduces the accuracy of the method. Many particle regularization methods popularly known as particle shifting techniques (PST) have been proposed which can be incorporated into WCSPH schemes Xu et al. (2009); Lind et al. (2012a); Sun et al. (2019a); Huang et al. (2019). These methods ensure that the particles are distributed more uniformly. Instead of displacing the particles directly, Adami et al. (2013) propose to use a transport velocity instead of the particle velocity to ensure a uniform particle distribution. This approach is also framed in the context of Arbitrary Lagrangian Eulerian (ALE) SPH schemes by Oger et al. (2016). A similar approach is used by Sun et al. (2019a) to incorporate the shifting velocity in the momentum equation with the -SPH scheme Marrone et al. (2011); Antuono et al. (2010). Ramachandran and Puri (2019) also employ a transport velocity formulation and additionally propose using the EDAC scheme Clausen (2013) in the context of SPH, which removes the need for an equation of state (EOS). The resulting method is accurate but does not converge with an increase of resolution. An alternative approach to ensure particle homogeneity is the approach of remeshing proposed by Chaniotis et al. (2002)

where the particles are periodically interpolated into a regular Cartesian mesh. The method can be accurate but the remeshing can be diffusive and makes the method reliant on a Cartesian mesh. In a subsequent development,

Hieber and Koumoutsakos (2008) employ remeshing but couple it with an immersed boundary method to deal with complex solid bodies. Recently, Nasar et al. (2019) modify the method introduced by Lind et al. (2012a) to devise an Eulerian WCSPH scheme that also uses ideas from immersed boundary methods to handle complex geometry.

To summarize the discussion in the context of convergence, some authors Quinlan et al. (2006); Schwaiger (2008); Macià et al. (2012) demonstrate numerical convergence for the derivative and function approximation. Many authors Kiara et al. (2013b); Huang et al. (2019); Adami et al. (2013); Chen and Beraun (2000); Zhang and Batra (2004); Sun et al. (2019b); Marrone et al. (2011) only show convergence in the form of plots that approach an exact solution with increasing resolution without formally computing the order of convergence. Some authors demonstrate second order convergence for simpler problems with a fixed particle configuration like the heat conduction equation Cleary and Monaghan (1999); Schwaiger (2008); Fatehi and Manzari (2011); Korzilius et al. (2017), the Poisson equation Macià et al. (2012), and the evolution of an acoustic wave Frontiere et al. (2017). Second order convergence has also been demonstrated for Eulerian SPH methods where the particles are held fixed Lind and Stansby (2016); Nasar et al. (2019) or where the particles are re-meshed Hieber and Koumoutsakos (2008). Some authors Rosswog (2015); Frontiere et al. (2017); Xu et al. (2009); Lind et al. (2012b); Ramachandran and Puri (2019); Dehnen and Aly (2012) show first order convergence for Lagrangian SPH schemes but this does not persist as the resolution is increased. Therefore, to the best of our knowledge, none of the contemporary Lagrangian SPH schemes appear to demonstrate a formal second order convergence for simple fluid mechanics problems like the Taylor-Green vortex problem for which an exact solution is known.

In this paper, we carefully construct a family of Lagrangian SPH schemes that demonstrate second order numerical convergence for the classic Taylor-Green vortex problem. We first study several commonly used SPH kernels in the context of function and derivative approximation using particles that are either in a Cartesian arrangement or in an irregular but packed configuration of particles encountered when employing some form of a particle shifting technique. We choose a suitable correction scheme that produces second order approximations. We then select a suitable kernel and smoothing radius based on this study. We then systematically study the various discretization operators along with suitable corrections. Our investigations are in two-dimensions although the results are applicable in three dimensions as well. Our numerical investigation covers a wide range of resolutions with our highest resolution using a quarter million particles with , where is the length of the domain. Once we have identified suitable second order convergent operators we carefully construct SPH schemes that display a second order convergence (SOC). We use the Taylor-Green vortex problem to demonstrate this. We also compare our results with those of several established SPH methods that are currently used. We study the accuracy, convergence, and also investigate the computational effort required. We construct both Lagrangian and Eulerian schemes that are fully second order convergent. We provide schemes that use either an artificial compressibility in the form of an equation state or using a pressure evolution equation.

Once we have demonstrated second order convergence for the Taylor-Green vortex problem we proceed to investigate the Gresho-Chan vortex Gresho and Chan (1990) problem as well as an incompressible shear layer problem Di et al. (2005) and look at how the lack of manifest conservation impacts the conservation of linear and angular momentum. In the interest of reproducibility, all the results shown in the paper are automatically generated through the use of an automation framework Ramachandran (2018), and the source code for the paper is available at https://gitlab.com/pypr/convergence_sph. In the next section we discuss the SPH method briefly and then proceed to look at the SPH kernel interpolation.

2 The SPH method

We define the SPH approximation of any scalar (vector) field

() in a domain by


where , is the kernel function, and is the support radius of the kernel. It is well known Monaghan (2005); Violeau (2012) that for a suitably chosen quadrature rule and a symmetric kernel which satisfies that,


Some of the widely used kernels are Gaussian Monaghan (2005), cubic spline Monaghan (1994), quintic spline, and Wendland quintic Wendland (1995). We note that in this work we take to be a constant.

When the kernel support is completely inside the domain boundary then we can evaluate the gradient of a function by taking the gradient of the kernel inside the integral. This approximation is also second-order in  Violeau (2012). We discretize the domain using particles having mass , and density . The discretization of the domain into particles introduces additional error in the approximation and is discussed in Quinlan et al. (2006). We can approximate the gradient of as,


where , is a measure of the volume of the particle, is the smoothness of the kernel at the edge of its support, and the sum is taken over all the particles under the support of the kernel. The value of is defined as the smallest order of derivative of the kernel at the edge of its support that is non-zero. For example, for cubic spline kernel, and for quintic spline kernel.

We note that the numerical volume is an approximation and solely depends upon the spatial distribution of the particles. The density may be computed for a particle using the summation density as,


Therefore, for constant mass we may write the volume as,


Since the function and its derivative approximation depends on the kernel , support radius , and the scaling factor , we perform a numerical study of the effect of the kernel in the next section.

3 Comparison of kernels

Figure 1: The unperturbed periodic particle and perturbed periodic particles.

In this section, we compare various SPH kernels for their accuracy and order of convergence in a discrete domain. We evaluate the error in function and its derivative approximation in a two-dimensional periodic domain. We simulate periodicity by copying the appropriate particles and their properties near the boundary such that the boundary particles have full support Randles and Libersky (1996). The particles are either placed in a uniform mesh (unperturbed) or in a packed arrangement referred as unperturbed periodic (UP) or perturbed periodic (PP), respectively. In order to obtain the packed configuration, the particles are slightly perturbed from a uniform mesh and their positions are moved and allowed to settle into a distribution with a nearly constant density using a particle packing algorithm Negi and Ramachandran (2021). The algorithm effectively ensures that the particles are not clustered and have minimal density variations. It mimics the effect of many recent particle shifting algorithms Colagrossi et al. (2012); Lind et al. (2012a); Huang et al. (2019). In the fig. 1, we show both the domains.

Name Radius Remark
- Gaussian Monaghan (2005) 3 0 Truncated for low
- Quintic spline Ramachandran and Puri (2019) 3 3 Tensile instability
- Cubic spline Monaghan (2005) 2 5 Paring and tensile instability
- Wendland Wendland (1995) 2 5 No tensile or pairing instability
- Wendland Wendland (1995) 2 8 Produces higher accuracy
- Wendland Wendland (1995) 2 11 Produces higher accuracy
Table 1: Kernels and their properties

In this paper, we compare the set of kernels listed in table 1. The set of kernels considered covers a wide range (high order, kernels having tensile instability and pairing instability Dehnen and Aly (2012); Monaghan (2000)). In order to assess the effect of for a kernel, we perform the numerical experiment proposed by Dehnen and Aly (2012). We evaluate particle density using eq. 4 for increasing the number of neighbors , for each of the kernels. The increase in corresponds to the scaling of the smoothing kernel using the parameter. In this numerical experiment, we change both the resolution and .

Figure 2: The particle density for different kernels with varying number of neighbors

In the fig. 2, we plot the absolute error in the particle density of one particle in an UP domain for different kernels with the change in under the kernel support. Clearly, the Wendland class of kernel shows a monotonic decrease in error with the increasing . However, in the case of the and kernels, the errors are an order less at a lower compared to Wendland class of kernels. The error in the kernel does not change significantly with the change in the compared to others. It is because we truncate the kernel to have compact support. In the , the error is lower than the in the entire plot. Therefore, we drop and in the subsequent investigations since it reaches the order of accuracy of when is approximately . High results in higher computational cost.

Figure 3: Order of convergence for the approximation of a function for different values in an UP domain. The dashed line shows the second order rate.
Figure 4: Order of convergence for derivative approximation for different values in an UP domain. The dashed line shows the second order rate.
Figure 5: Order of convergence for function approximation for different values in a PP domain. The dashed line shows the second order rate.
Figure 6: Order of convergence for the derivative approximation for different values in a PP domain. The dashed line shows the second order rate. The eq. 8 is used for the approximation.

We compare the four kernels , , and for convergence of function and its gradient approximation. We consider the field,


Given a function and its approximation , we evaluate the error using,


where is the total number of particles in the domain. Since the and kernels have support radius of whereas, the and kernel have support radius of , we set the such that the is same in an UP domain. Therefore, when for (or ), we take for the (or ). For the convergence study, in this paper, we consider , , , , , and resolutions for all the test cases unless stated otherwise.

3.1 Unperturbed periodic domain

In fig. 3, we plot the in the function approximation as a function of the resolution for different values of in a UP domain. We observe similar error values for all the kernels except . We obtain second-order convergence (SOC) in an UP domain upto a considerably high resolution of as expected, for all the kernels Kiara et al. (2013a).

In fig. 4, we plot the error in the derivative approximation of the function in eq. 6 in a UP domain. The and kernels show a better convergence rate compared to and for lower . The kernel does not show SOC even at , since we use a truncated Gaussian. The and kernel shows SOC only when . The kernel shows SOC at however, a reasonable convergence can be seen for as well.

3.2 Perturbed Periodic domain

In an SPH simulation, the particles advect with different velocities, and thus the distribution of particles is no longer uniform. Some particle shifting techniques (PST) can be used to make the particle distribution uniform Xu et al. (2009); Lind et al. (2012a). Thus, it is essential to observe the convergence rate in the PP domain as well.

In fig. 5, we plot the error of the approximation of the field given in eq. 6 as a function of resolution in a PP domain for different values. The convergence rates tend to zero for higher resolution for low value of for all the kernels. The kernel performs worse than the kernel at lower values however, the errors are significantly lower in when the value increase. On comparing and , the error plot looks exactly same except when .

Figure 7: The convergence of the derivative approximation with different kernels when the gradient correction of Bonet and Lok (1999) is used on a PP domain. Here we use .

The SPH approximation of the gradient of a function is not even zero order accurate in a perturbed domain Nguyen et al. (2008); Kiara et al. (2013a). The derivatives diverge when we evaluate it using eq. 3. We use a zero order consistent method proposed by Monaghan (1994) to compare the kernels. We write this approximation as,


In fig. 6, we plot the error in the function derivative approximation as a function of resolution using eq. 8 in a PP domain for different values and kernels. Clearly, the approximation for all the kernels shows at least zero-order convergence. The kernel does not show SOC for high which is the same as observed in the case of the UP domain. The accuracy in the case of and oscillates when going from lower to higher values. Zhu et al. (2015) suggest that one should increase the as one increases the resolution but given the inconsistent behavior of the and kernels; these may not be suitable for that approach. The zero-order convergence rate occurs due to dominance of discretization error (the term in eq. 3) when the resolution increases in the PP domain.

In order to observe the effect of kernel gradient correction, we apply the correction proposed by Bonet and Lok (1999) for all the kernels. In fig. 7, we plot the error in the derivative approximation as a function of resolution. Clearly, all the kernels , , , and show more or less the same behavior. Thus, we can choose any of these kernels for our convergence study of the WCSPH schemes. In the figure, it can be seen that the and kernels do not sustain the second-order behavior. Therefore in this work, we choose the kernel with for all the test cases henceforth.

4 Comparison of SPH discretization operators

The SPH method is widely used to solve fluid flow problems. In this work, we focus on weakly-compressible SPH schemes that are used to simulate incompressible fluid flows. We write the Navier-Stokes equation for a weakly compressible flow along with the equation of state (EOS) as,


where , , and are the fluid density, velocity, and pressure, respectively, is the kinematic viscosity, is the reference fluid density, and is the artificial speed of sound in the fluid.


We note that the fluid density, is independent of the summation density . Normally, in SPH simulations these two are treated as the same and we discuss the reasons behind this choice in section 5. The property, does not depend upon particle configuration and should be prescribed as an initial condition.

For a scheme to be SOC, the SPH discretization of the RHS of all the equations in eq. 9 must be . Thus, in this section, we perform the comparison of the schemes by comparing the convergence of different approximations for the gradient, divergence, and Laplacian operators.

There are many different ways to discretize the momentum equation as can be seen from Gomez-Gesteria et al. (2010); Sun et al. (2018); Adami et al. (2013); Violeau (2012); Frontiere et al. (2017); Sun et al. (2017). One of the key features of the discretization of the momentum equation is to ensure linear momentum conservation. However, some researchers trade conservation for better accuracy Sun et al. (2017) and others use a conservative form but chose to loose accuracy Frontiere et al. (2017). In view of this, we consider both conservative and non-conservative discretization in the present study.

4.1 Comparison of approximation

Name Expression Used in
sym1 SPH Sun et al. (2018)
sym2 , WCSPH Gomez-Gesteria et al. (2010),
TVF Adami et al. (2013); Zhang et al. (2017); Ramachandran and Puri (2019), ISPH Shao and Lo (2003)
asym WCSPH Monaghan (2005)
Table 2: Different gradient approximations for . The column “expression” is assumed to be summed over the index over all the neighbor particles inside the kernel support. The term for gradient comparison. is the density averaged pressure.
Figure 8: The rate of convergence in UP (left) and PP(right) domains for various pressure gradient listed in table 2. The dashed line shows the SOC rate. The _bc and _lc suffixes represent the corresponding form with Bonet correction Bonet and Lok (1999) and Liu correction Liu and Liu (2006). The sym_sl is the sym1 formulation with symmetrization of kernel proposed in Dilts (1999)
Name error Order
asym_bc 1.01e-02 1.97 1.11e-04 1.99
asym 3.10e-04 1.00 1.01e-03 0.98
sym1_bc 1.08e-02 1.87 7.65e-02 -0.75
sym1_lc 1.01e-02 2.35 1.11e-04 1.99
sym1 -2.53e-11 1.05 7.65e-02 -0.76
sym2_bc 1.10e-02 1.87 7.65e-02 -0.75
sym2_lc -1.90e-11 2.43 7.65e-02 -0.76
sym2 -1.87e-11 1.01 7.65e-02 -0.76
sym_sl -1.73e-11 2.50 6.46e-02 -0.72
Table 3: The ratio showing the total force in the system due to lack of conservation in the approximation, the time taken relative to the asym formulation, the error for particle in a PP domain, and last column shows the order of convergence for all the kind of formulations considered in the first column.

In this section, we compare various pressure gradient approximations. In the table 2, we list the gradient approximations considered in this study. The sym1 and sym2 are the symmetric, conservative form of the gradient approximation. We note that conservative forms have . The asym is the asymmetric form. Since the SPH kernel gradient does not show SOC in a perturbed domain Quinlan et al. (2006), we also consider the kernel correction employed to each of the approximation. In this paper, we refer to the correction proposed by Bonet and Lok (1999) as Bonet correction and the one proposed by Liu and Liu (2006) as Liu correction. We add the suffix _bc, and _lc respectively in the plots and tables to indicate these corrections. The application of corrections renders the symmetric forms non-conservative, we use the method of symmetrization of the kernel proposed by Dilts (1999) to again make it conservative. We refer to this formulation as sym_sl which we write as


where is the Liu correction applied to the kernel gradient 111We select the sym1 formulation over sym2 as the latter does not perform well with the linear correction c.f A for details.. This formulation is used in the scheme proposed by Frontiere et al. (2017).

In order to compare the convergence, we consider a pressure field, . We determine the error using eq. 7, where is the pressure gradient evaluated using the approximation and is the exact pressure gradient. The exact pressure gradient, . We compare only the x-component of the results. In fig. 8, we plot the error in the various gradient approximations discussed above in both an UP and PP domain. In the UP domain, barring the sym2_lc, all the corrected gradient approximations behave the same, whereas the uncorrected gradients do not display SOC. The corrected versions retains SOC even at high resolution since it reduces the discretization error in the approximation Fatehi and Manzari (2011); Quinlan et al. (2006). We also observe that with the correction the second term involving the term is zero in an UP domain leading to the same expression.

In the case of the PP domain, we observe that both sym1 and sym2 and their corresponding _bc versions overlap. The symmetric formulations show an increase in the error in the approximation with increasing resolution as suggested in Fatehi and Manzari (2011). Furthermore, the Bonet correction does not correct the symmetric formulations. Clearly, the asym formulation shows better convergence, and the Bonet correction version shows SOC. Therefore, the Bonet correction can be applied only when an asymmetric formulation is employed. Moreover, the Liu correction only corrects the symmetric form sym1, which suggests that the sym2 cannot be corrected using traditional correction techniques. We investigate this in more detail in A. Finally, the sym_sl method has a slightly lower error but looses SOC behavior due to the symmetrization of the kernel gradient. Frontiere et al. (2017) reported the similar behavior.

We also compare the linear momentum conservation and time taken to evaluate the gradient for the case with particles. As shown in Bonet and Lok (1999), linear momentum is conserved when the total force, , where the sum is taken over all the particles and . In table 3, we tabulate the ratio of total force to the maximum force (), the time taken to evaluate the gradient scaled by the minimum time taken by all the methods, and the error with the order of convergence 222

In this paper, we report order of convergence by fitting a linear regression line and finding its slope.

, for all the formulations plotted in fig. 8. As expected, all the symmetric forms of approximation have zero total force. The asymmetric formulation has a very small total force. Clearly, the use of Bonet corrections increases the total force and slows down the computation by a factor of , whereas the Liu correction makes it times slower. The sym_sl formulation shows zero residual force as expected. Using the table 3, we can see that asym_bc and sym1_lc show SOC and have a very low total force which makes them a suitable candidate for a scheme with SOC.

4.2 Comparison of approximation

Figure 9: The rate of convergence UP (left) and PP(right) domains for velocity divergence in eq. 11. The dashed line shows the SOC rate. The suffix _bc represents the corresponding form with Bonet correction.

A zero order consistent SPH approximation for the divergence operator Monaghan (2005) is,


We refer to the approximation given in eq. 11 as div. We apply the Bonet correction as done in the case of gradient approximation for a first-order consistent approximation. We refer to the corrected form as div_bc.

We consider the velocity field, . The divergence of the velocity is given by, . We evaluate the error in the approximation using eq. 7. In fig. 9, we plot the error in the divergence approximation in an UP and PP domain. The uncorrected approximation does not display SOC since the discretization error dominates as we approach higher resolutions. Clearly, the corrected form shows SOC even in the case of a PP domain.

In order to evaluate the accuracy of the approximation in a divergence-free field, we consider the velocity field,

Figure 10: The rate of convergence UP (left) and PP(right) domains for velocity divergence for a divergence-free field. The dashed line shows the SOC rate. The suffix _bc represents the corresponding form with Bonet correction.

In the fig. 10, we plot the error using eq. 7 in the divergence computation as a function of resolution for the UP and PP domain. Clearly, the divergence is zero in a UP domain owing to the symmetry of the particles. However, the error in PP domain remains about the same order as seen in the case of general field in fig. 9 333We discuss the source of this behavior in B. Clearly, the Bonet correction does not correct this issue. We observe the implication of this behavior when we compare the schemes in section 5.4.

The continuity equation corresponds to the mass conservation of the system, since mass of each particle is kept constant, we satisfy the global conservation of mass implicitly.

4.3 Comparison of approximation

Name Expression Used in
Cleary WCSPH Cleary and Monaghan (1999); Gomez-Gesteria et al. (2010)
Fatehi modified WCSPH Fatehi and Manzari (2011)
tvf TVF Adami et al. (2013), EDAC Ramachandran and Puri (2019)
coupled Bonet and Lok (1999)
Table 4: The various approximations of . The column “expression” is assumed to be summed over the index over all the neighbor particles inside the kernel support. The term are calculated using first-order consistent formulation i.e. asym_bc
Figure 11: The rate of convergence UP (left) and PP(right) domains for various approximations of the Laplacian operator in table 4. The dashed line shows the SOC rate. The suffixes _bc and _lc represent the corresponding form with the Bonet correction and Liu correction, respectively. The Fatehi_c refers to the fatehi formulation with the correction proposed by Fatehi and Manzari (2011) (See C).
Name error Order
Cleary_bc -1.28e+00 1.95 6.55e-02 -0.83
Cleary_lc -2.15e-01 2.43 4.59e-02 -0.79
Cleary -1.08e-10 1.12 6.54e-02 -0.84
Fatehi_c 1.69e+00 4.10 2.88e-04 1.30
Fatehi 1.30e+00 2.70 8.43e-04 0.68
Korzilius 1.61e+00 4.66 2.49e-04 1.70
coupled_bc 1.61e+00 3.05 2.54e-04 1.95
coupled 1.30e+00 2.59 8.38e-04 1.46
tvf_bc -1.26e+00 2.07 6.80e-02 -0.66
tvf_lc -2.16e-01 2.35 5.01e-02 -0.56
tvf -1.06e-10 1.00 6.77e-02 -0.67
Table 5: The ratio showing the total force in the system due to lack of conservation of the approximation and the time taken, relative to the tvf formulation, and error for particle case in a PP domain. The last column shows order of convergence for all the methods listed in first column.

In this section, we compare various approximations for the Laplacian operator listed in table 4. We refer to the symmetric formulations of Cleary and Monaghan (1999); Brookshaw (1985) as Cleary, and those of Adami et al. (2013) as tvf. These ensure that linear momentum is conserved. We also consider the coupled formulation used by Bonet and Lok (1999); Nugent and Posch (2000) and refer to these as coupled. This formulation shows oscillations in the approximation when the initial condition is discontinuous. However, to remedy this, one can perform a first-order accurate approximation near the discontinuity and then perform this approximation as shown in Biriukov and Price (2019). We consider the improved formulation proposed by Fatehi and Manzari (2011) referred as Fatehi. Both the coupled and Fatehi formulations are asymmetric. These formulations are performed in two steps where the first step involves computation of velocity gradient for each particle. We also consider the kernel correction applied to each of these formulations. In the case of Cleary, tvf, and coupled methods, we use the standard Bonet and Liu corrections. However, in the case of Fatehi

, we use the correction tensor proposed by

Fatehi and Manzari (2011) given by


where the subscripts are SPH summation indices and superscripts are in tensor notation 444The term is widely used in SPH literature and so we use tensor indices in the superscript to derive eq. 13 in C.. We refer to this corrected formulation as Fatehi_c. Additionally, the second derivative can also be obtained by taking the double derivative of the kernel Chen and Beraun (2000); Zhang and Batra (2004); Korzilius et al. (2017). Therefore, we also consider the method proposed by Korzilius et al. (2017) that remedies the deficiencies in earlier approaches where the second derivative was employed. We obtain the second derivative of a scalar field using,


where is the operator, . The gradient is approximated using the asym_bc formulation. The correction is given by


where and is the Bonet correction matrix (see (40)). We refer to this formulation as Korzilius.

In the fig. 11, we plot the rate of convergence for the various formulations discussed above in both UP and PP domains. In the UP domain, all the methods at least show zeroth-order convergence. All methods without corrections suffer from high discretization error that dominates at higher resolutions Fatehi and Manzari (2011). When either Bonet or Liu corrections are employed, Cleary, coupled, Fatehi_c, and Korzilius methods show SOC. The coupled method is approximately half an order less accurate as compared to Cleary and Fatehi. The accuracy of Korzilius method is in between the coupled and Fatehi method. The tvf method is very inaccurate as the discretization error increases due to the introduction of .

It is important to note that in the PP domain, the symmetric methods diverge due to discretization error of , where is the deviation from the regular particle arrangement Fatehi and Manzari (2011). Only the coupled, Korzilius, and Fatehi methods show a positive convergence rate. On applying the corresponding correction, the coupled, Korzilius, and Fatehi methods improve. The accuracy for coupled, Korzilius, Fatehi is maintained as observed in the case of UP domain.

In the table 5, we tabulate the total force as a result of the approximation, the time taken for the approximation, and the error on a PP domain consisting of particles with the order of convergence in the last column for each method plotted in fig. 11. We observe a similar increase in computational time due to the Bonet and Liu corrections as seen in the case of gradient approximation. The coupled Korzilius, and Fatehi formulation have even higher computational cost due to the additional step of velocity gradient computation. The Korzilius method requires additional time since the double derivative of the kernel is involved. The Fatehi_c method has an additional step where we compute the second-order tensor in eq. 13 for each particle resulting in a further increase in computation time. We observe a similar increase in total force when an asymmetric version of the formulation is employed, as seen in the case of gradient approximation. Clearly, both the coupled, Korzilius, and Fatehi formulation results in an equal amount of total force resulting in a lack of conservation of linear momentum. In order to get a SOC approximation, we can use either of coupled_bc, Korzilius, or Fatehi_c formulations for viscous force estimation.

5 Comparison of schemes

In this section, we compare the convergence of several SPH schemes. However, given the results of the previous section, we first propose a SOC scheme using the right combination of discretizations.

5.1 A SOC WCSPH scheme

In this section, we discuss a SOC scheme for weakly compressible fluid flows. We consider the fluid density as a property, carried by a particle. The numerical density, , and volume, are a function of the surrounding particle distribution. The mass, of the particle satisfy where is the physical volume of the particle occupied by the particle and is the numerical volume used for integration. Thus, we can approximate the fluid density using the standard SPH approximation given by


In case of weakly compressible SPH, the requirement of linear momentum conservation condition may be relaxed and is only satisfied approximately Oger et al. (2007); Bonet and Lok (1999); Sun et al. (2017). Therefore, we use the SOC approximations that are non-conservative discussed in section 4. In table 6, we list all the discretizations that we can employ to obtain a SOC WCSPH scheme.

Operators Possible discretization for SOC
Gradient asym_c, sym1_l
Divergence div_c
Laplacian coupled_c, Fatehi_c, Korzilius
Table 6: The operators and their discretization suitable for a SOC scheme.

Furthermore, one can solve the fluid flow equations by assuming the Lagrangian nature as well as the Eulerian nature. We discuss the scheme for both the context in the following sections.

5.1.1 Lagrangian method

In the Lagrangian description, the continuity equation and the momentum equation are given by


In order to evaluate the RHS of the above equations, one may employ any method listed in table 6. The pressure is evaluated using an equation of state given by


where , is the reference density, and is the reduced speed of sound. We note that the linear equation of state where in eq. 18, , works equally well. We integrate the particles in time using a Runge-Kutta order integrator.

Since we use an asymmetric form of the pressure gradient approximation, particles tend to clump together due to absence of a redistributing background pressure Sun et al. (2018). We use the iterative particle shifting proposed by Huang et al. (2019) after every few iterations to redistribute the particles. We compute the shifting vector for the iteration (of the shifting iterations) using


where . The new particle position,


is computed. The particles are shifted until the criterion,


is satisfied up to a maximum of 10 iterations, where ,

is the value for uniform distribution, and

is an adjustable parameter. In order to keep the approximation of the particle accurate, we update the particle properties after shifting by,


where is the final position after iterative shifting, is the property to be updated, and is the gradient of the property on the last position computed with the Bonet correction. In a variation of the above scheme discussed in section 5.4, we observe that usage of non-iterative PST proposed by Sun et al. (2019b) results in slightly higher errors but still retains its SOC. We refer to the scheme discussed above as L-IPST-C (Lagrangian with iterative PST and coupled_c viscosity formulation). Similarly the method using Fatehi_c and Korzilius formulation are referred as L-IPST-F and L-IPST-K respectively. We note that we only perform the IPST step every timesteps rather than at every timestep.

5.1.2 Eulerian method

In the Eulerian description, the continuity equation is written as


Since the fluid density is not same as the particle density we do not ignore this term, this is unlike what is done by Nasar et al. (2019). The momentum equation is written as


We discretize all the terms using SOC operators listed in table 6. We perform time integration using the RK2 integrator; however, we note that the positions of the particles are not updated.

5.2 Comparison with other SPH schemes

As discussed in the section 1, many weakly-compressible SPH schemes have been proposed. In this section, we compare the following schemes:

  1. TVF : Transport velocity formulation proposed by Adami et al. (2013).

  2. -SPH : The improved -SPH formulation proposed by Sun et al. (2019a).

  3. EDAC : Entropically Damped artificial compressibility SPH formulation proposed by Ramachandran and Puri (2019).

  4. EWCSPH : The Eulerian SPH method proposed by Nasar et al. (2019).

  5. L-IPST-C : The Lagrangian method with iterative PST and coupled_c viscosity formulation discussed in section 5.1.1.

In order to compare these schemes, we consider the Taylor-Green vortex problem. We choose this problem, since it is periodic, has no solid boundaries, and admits an exact solution. The solution of the Taylor-Green problem is given by


where , where is the Reynolds number of the flow. We consider and . For the Lagrangian schemes, we consider a perturbed periodic (PP) arrangement of particles shown in the fig. 1 for different resolutions. At we initialize the pressure and velocity using eq. 25 for all the schemes. Since the fluid density is a function of pressure, we initialize density inverting eq. 18. In the case of the EWCSPH scheme we consider an unperturbed periodic (UP) arrangement of particles and initialize the using the prescribed pressure. We compute the error in pressure and velocity by


where, is the smoothing length of the kernel, is the timestep, is the total number of particles in the domain, is either pressure or velocity, and is the exact value obtained using eq. 25. The particle spacing, is set according to the resolution. We consider resolutions of to particles in a periodic domain. In order to isolate the effect of spatial approximations on the convergence, we set timestep, , where is set corresponding to highest resolution, , for all the simulations. We run all the simulations for timestep and observe convergence. We choose one timestep since most of the scheme considered diverge.

Figure 12: Convergence of error in pressure (left) and velocity (right) with the change in resolution. The , , and only timestep taken.
2.19e-05 1.49 5.69e-05(0.00) 2.03e-01(0.00)
EDAC 1.88e-07 1.32 1.68e-05(-1.24) 7.00e-05(-0.07)
EWCSPH 7.03e-15 2.25 3.13e-07(0.38) 2.10e-05(0.00)
L-IPST-C 1.34e-05 3.36 1.18e-07(1.42) 6.41e-05(0.00)
TVF 3.70e-16 1.00 9.45e-04(-0.88) 2.88e-01(-0.07)
Table 7: Table showing total force w.r.t. the maximum force in the domain and the time taken for iteration w.r.t. the TVF scheme for all the schemes.

In the fig. 12, we plot the error evaluated using eq. 26 for pressure and velocity in the domain for different schemes. Clearly, none of the schemes show convergence in pressure. It is because the initial velocity is divergence-free, so there is no change in density and thereby pressure. We observe that the EDAC, EWCSPH, and L-IPST-C schemes are almost four orders more accurate than the TVF and SPH schemes. In case of both the TVF and SPH schemes, we link the pressure with particle density , which is a function of the particle configuration. The particle positions are a result of the particle shifting, and therefore, the pressure is incorrectly captured. On the other hand, the other schemes either use a pressure evolution equation (EDAC) or a fluid density to evaluate pressure. In the case of the EWCSPH scheme, we initialize density using the pressure values in the eq. 18 which results in better accuracy.

The error in velocity diverges in the case of the TVF and EDAC schemes since these use a symmetric form of type sym2 in table 2 to discretize the momentum equation. Whereas, in the case of the SPH scheme, sym1 type of discretization is employed leading to less errors. Moreover, the SPH scheme uses a consistent formulation and both TVF and EDAC schemes are inconsistent when the shifting (transport) velocity is added to the momentum equation Sun et al. (2019a). The EWCSPH, and L-IPST-C formulations show convergence (not second-order) as expected. We observe that in the velocity convergence a constant leading error term dominates resulting in flattening at higher resolutions. Since, we use second-order accurate formulations in L-IPST-C and EWCSPH 555It is second-order accurate since a uniform stationary grid is used. formulations, the only equation which is not converging with resolution is the equation of state (EOS).

In this section, we focus on highlighting the effect of using fluid density different than the numerical density . It allows a better convergence rate and independence of density from particle positions. In contrast to this, the use of numerical density as a function of particle position is consistent with the volume used for the SPH approximation. In TVF, -SPH, EDAC and, EWCSPH schemes, we make no such distinction, and use the fluid density in the numerical volume . The poor convergence for these schemes show that it is important to treat the fluid and numerical densities differently.

In the table 7, we tabulate the total force and the time taken by the scheme for one timestep with the errors and order of convergence in pressure and velocity for the resolution case. It is clear that the TVF and EWCSPH schemes conserve linear momentum and take the least amount of time. The EDAC and the -SPH scheme do not conserve linear momentum exactly. In the case of the EDAC scheme the use of average pressure in the pressure gradient evaluation results in lack of conservation. Whereas, in the case of -SPH the asymmetry of the shifting velocity divergence causes lack of conservation. The L-IPST-C scheme is known to be non-conservative; however, the value is comparable to other schemes. The time taken by the L-IPST-C scheme is significantly higher due to the evaluation of correction matrices.

5.3 The effect of

In the schemes discussed in the previous section, except the EDAC scheme 666The EDAC scheme uses a pressure evolution equation., we impose artificial compressibility (AC) using the EOS, which is accurate Chorin (1967); Clausen (2013), where is the Mach number of the flow. Chorin (1967) originally proposed this method to obtain steady-state solutions of an incompressible flow. Some authors have used artificial compressibility with dual-time stepping to achieve truly incompressible time-accurate results Fatehi et al. (2019); Ramachandran et al. (2021). We achieve the incompressibility limit when .

Figure 13: Convergence rates for pressure (left) and velocity (right). The L-IPST-C and L-IPST-F methods are compared for different values of .
L-IPST-C 7.13e-05 1.00 2.03e-04(1.38) 3.66e-03(0.93)
L-IPST-C 9.31e-05 2.15 7.94e-05(1.78) 2.09e-03(1.60)
L-IPST-C 1.44e-04 3.76 5.32e-05(1.93) 2.98e-03(1.85)
L-IPST-F 7.11e-05 1.23 1.80e-04(0.85) 3.77e-03(0.73)
L-IPST-F 5.99e-05 2.84 4.81e-05(1.44) 2.54e-03(1.31)
L-IPST-F 1.66e-04 5.46 1.36e-05(1.98) 3.52e-03(1.65)
Table 8: Comparison of total force, time taken relative to the L-IPST-C with , error in velocity and order for different values of .

This section compares the convergence of the L-IPST-F and L-IPST-C scheme with change in AC parameter i.e. the speed of sound . We consider the Taylor-Green problem; however, we run the simulation for . In fig. 13, we plot the error in the pressure and velocity for both the schemes and different . We observe that both L-IPST-C and L-IPST-F methods are affected by the change in value significantly as expected. In case of pressure, with the increase in the value, the error in the lower resolutions increases; however, for higher resolution it improves. Clearly, we attain the increase in the order of convergence in case of the pressure due to increased error scales at lower resolutions. The increase in error is attributed to the inability of the SPH operators to correctly capture a divergence free velocity field as discussed in section 4.2. However, on looking at the velocity convergence, both the schemes attain SOC even at higher resolutions. We observe, though, at lower resolutions, the error in the pressure increases with the value. As observed in the case of the Laplace operator comparison in section 4.3, the use of Fatehi_c discretization offers better accuracy.

In the table 8, we tabulate the total force, relative time, and the error in pressure and velocity with the order of convergence for particles. We observe that at a higher value, the total force is less compared to the simulation when values are lower. From the table, we can see that the use of L-IPST-F scheme offers better accuracy at the cost of the extra time taken. We also note that one can choose to use lower values of at lower resolutions and increase the value as the resolution increase to get better accuracy in pressure. Doing this does not affect the convergence in velocity.

5.4 Study of other SOC variants

In this section, we show that the scheme presented in the section 5.1.1 (L-IPST-C) can be easily converted into other forms for improved accuracy and ease of calculation. We note that regardless of the set of governing equations employed, the discretizations from table 6 must be used to achieve SOC.

In order to remove high frequency oscillations, one could modify the continuity equation given by,


where, is the damping constant, where . This corresponds to the -SPH scheme Antuono et al. (2010). In this case we also use the linear equation of state to evaluate given by


The following are different variations of the basic scheme:

  1. Using different PST : One could use either IPST proposed by Huang et al. (2019) or the non-iterative PST proposed by Sun et al. (2019a). The properties like and need to be updated using first order Taylor expansions given by


    where is the desired property. We use the coupled formulation for the viscosity and non-iterative PST. We refer to this method as L-PST-C.

  2. Using pressure evolution: On taking the derivative of EOS in eq. 28 w.r.t. time and using the eq. 27, we get the pressure evolution equation given by


    This is very similar to the EDAC pressure evolution equationRamachandran and Puri (2019). The value of can be evaluated from the EOS in eq. 28 given by


    We employ the coupled formulation for viscosity and use IPST for regularization. We refer to this method as PE-IPST-C.

  3. Using remeshing for regularization: The regularization step performed using PST in the L-IPST-C method can be replaced with the remeshing procedure of Hieber and Koumoutsakos (2008). The remeshing is performed using the kernel given by,


    where , where is the initial particle spacing. The properties on the regular grid are computed using


    where are points on a regular Cartesian mesh. The remeshing procedure can be performed every few steps; however, we perform remeshing after every timestep. We use the coupled formulation for viscosity. We refer to this method as L-RR-C.

  4. Including regularization in the form of shifting velocity: the methods of Adami et al. (2013); Oger et al. (2016); Sun et al. (2019a) use shifting by perturbing the velocity of the particles and adding corrections to the momentum equation. Thus the particles are advected using the transport velocity, and the displacement is given by,


    The new continuity and momentum equations are given by


    where . In this method, we employ SOC approximations mentioned in table 6. We use the PST proposed in Sun et al. (2019a) for this scheme along with the coupled formulation for viscosity. We refer to this method as TV-C.

  5. Eulerian method: The Eulerian method solves the equation of motion on a stationary grid. This can be derived using the TV-C method by setting . This substitution makes the transport velocity in the eq. 34 equal to zero, thus the particle does not move. The modified equation on setting in eq. 35, we get


    Therefore, we recover the governing equation for the Eulerian method. We note that unlike Nasar et al. (2019), we retain the last term in the continuity equation. We use the coupled formulation to discretize viscous term. We refer to this method as E-C (Eulerian method with coupled_c viscosity formulation) method.

Figure 14: Convergence rates for pressure (left) and velocity (right) of different variants of SOC schemes.
E-C 8.13e-15 1.19 5.11e-05(1.94) 1.50e-03(0.71)
EWCSPH 1.07e-14 1.00 3.25e-05(1.57) 1.13e-02(-0.20)
L-IPST-C 1.44e-04 1.53 5.32e-05(1.93) 2.98e-03(1.85)
L-PST-C 2.64e-06 2.02 7.31e-05(1.84) 8.93e-03(1.68)
L-RR-C 1.51e-18 1.41 3.74e-05(1.92) 9.11e-04(0.65)
PE-IPST-C 6.17e-05 1.65 5.05e-05(1.96) 3.01e-03(1.84)
TV-C -5.01e-06 1.88 1.06e-04(2.18) 1.51e-02(2.14)
Table 9: Comparison of total force, time taken relative to L-IPST-C with , error in velocity and pressure for variation of schemes with .

We simulate the Taylor-Green problem using for a duration of with different resolutions. In addition, we study the performance of the EWCSPH scheme since it is computationally efficient and accurate. In the fig. 14, we plot the error in pressure and velocity for all the schemes. In the table 9, we tabulate the total force, relative time, error in pressure and velocity at resolution, and the order of convergence.

The L-IPST-C and PE-IPST-C overlap in both the pressure and velocity convergence plots, and these are both approximately second-order. Compared to L-IPST-C, the L-PST-C shows lower convergence rate, and TV-C shows higher order of convergence, whereas E-C and L-RR-C show very poor convergence rates in pressure; however, the L-RR-C method shows very low errors in pressure. The EWCSPH has a negative convergence rate in pressure. While the TV-C shows a high convergence rate, it has much larger errors than all the other schemes considered for both pressure and velocity. The E-C, L-IPST-C, L-RR-C, and PE-IPST-C show high convergence rates in velocity as expected. The L-PST-C shows slightly high error and convergence rate. The EWCSPH shows a lower convergence rate of but is the most accurate of all the schemes regarding the velocity error.

The TV-C scheme shows low accuracy since we perform the shifting using an additional term in the momentum equation compared to the PE-IPST-C and L-IPST-C. This decrease in accuracy is also visible in the case of velocity. The L-PST-C scheme show higher error suggesting that the non-iterative PST does not perform the required amount of regularization. Both L-RR-C and E-C are comparable and most accurate. These schemes have lower error since the particles are fixed on a cartesian grid resulting in accurate computation of divergence as discussed in section 4.2. The pressure convergence flattens since it reaches the limit of accuracy possible with this value of and further accuracy may be seen by increasing this further.

Clearly, the total force in case of E-C, L-RR-C, and EWCSPH scheme is zero since we compute the acceleration on a uniform Cartesian grid of particles. However, the total force in other schemes are accurate to order . The times taken shows that L-PST-C is the highest since we apply the PST at every timestep, the TV-C involves many terms in the equations and therefore takes a lot of time. The E-C, and EWCSPH take the least time since they do not use a PST 777These methods can be made even faster since the neighbors need not be updated, and the correction matrices can also be computed once and saved.. The L-RR-C, L-IPST-C, and PE-IPST-C take a similar amount of time.

5.5 Conservation of SOC schemes

Thus far we, have looked at the convergence of the various schemes. In this section, we look at the schemes listed in table 9 from the perspective of conservation of linear and angular momentum.

5.5.1 The Grehso vortex

We consider the Gresho vortex problem Gresho and Chan (1990), which is an inviscid incompressible flow problem having the pressure and velocity fields given by,

Figure 15: The velocity of particles with the distance from the center of the vortex (left) and the -component of the total linear momentum (right) for all the schemes.
Figure 16: The velocity of particles with the distance from the center of the vortex (left) and the the -component of the total linear momentum (right) for the variation of SOC scheme.

We consider an unperturbed periodic domain of size with the center at . We set the kinematic viscosity, , and the time step and other properties as done in the Taylor-Green problem. The problem is simulated until . Since the problem is inviscid, we expect the scheme to retain the velocity and pressure field. We do not use artificial viscosity in the simulations for any of the schemes. We perform the simulation of all the schemes listed in table 7 except the EWCSPH scheme 888We discuss The failed simulations are discussed in the E.. We note that using an initial perturbed particle configuration results in very diffused results for all schemes except the L-IPST-C.

In the fig. 15, we plot the velocity of the particles with the distance, from the center (on left) and the -component of the total linear momentum with time for a particle simulation. The L-IPST-C scheme retains the velocity profile very well. The -SPH, EDAC, and TVF schemes show diffusion due to inaccuracy in the pressure gradient evaluation. Except for the TVF scheme, the rest show a finite increase in the momentum bounded at . Clearly, approximate linear momentum conservation is sufficient to obtain accurate results in the case of weakly compressible flows.

We also perform the simulations with different versions of the SOC scheme listed in table 9 999The L-RR-C, TV-C, and E-C scheme fail to complete the simulation, and these are discussed in the E. In the fig. 16, we plot the velocity with the distance from the center and the -component of the linear momentum with time for particle simulation. Clearly, all the schemes are accurate and approximately conserve linear momentum as expected.

Figure 17: The angular momentum variation with time for Gresho-Chan vortex for different schemes.

In fig. 17, we show the angular momentum variation with time for different schemes. None of the schemes conserve angular momentum, but for the SOC schemes, the variations are very small and .

5.5.2 The incompressible shear layer

The incompressible shear layer simulates the Kelvin-Helmholtz instability in an incompressible flow. This test case produces non-physical vortices for the schemes where the operators are under resolved even when the scheme is convergent Di et al. (2005). The initial condition for the velocity in x direction is given by


where . In order to begin the instability, a small velocity is given in y direction,


where . We consider a small viscosity . We simulate the problem using all the schemes listed in table 7. In fig. 18 and fig. 19, we plot the vorticity field for the schemes 101010The L-RR-C method failed to run due to discontinuity in the initial velocity field. discussed in this paper. We note that unlike the inviscid problem of Gresho-Chan vortex, the scheme EWCSPH, TV-C and E-C shows results matching other SOC schemes. In fig. 18, we observe that TVF scheme and -SPH scheme show high frequency oscillations and while the EDAC scheme is much better; However, it shows some undesired artifacts.

Figure 18: Vorticity contour plot for resolution for all the schemes.
Figure 19: Vorticity contour plot for resolution for all the schemes.

6 Cost of computation

Figure 20: The error in velocity with respect to the time taken to evaluate 5000 timesteps for all the schemes discussed in the previous sections.

In this section, we compare the cost of computation of all the schemes considered in this study. We simulate the Taylor-Green problem for timesteps with , , and resolutions for all the schemes. We use Intel(R) Xeon(R) CPU E5-2650 v3 processor and ran all the simulations in serial. In fig. 20, we plot the error in velocity computed using the eq. 26 as a function of time taken for the simulation. Clearly, all the SOC schemes are close to each other in terms of errors. The E-C and EWCSPH scheme takes very less time and are very accurate; however, EWCSPH is not convergent in pressure as shown in section 5.2. The EDAC scheme has lower error comparable to the SOC schemes; however, its convergence rate reduces with increase in resolution. We show that despite having higher time taken by the SOC schemes, they achieve higher accuracy with a fewer number of particle also, for some schemes, these accuracy levels are not achievable at all.

We implement the schemes using the open source PySPH 

Ramachandran et al. (2020) framework and automate the generation of all the figures presented in this manuscript using the automan framework Ramachandran (2018). The source code is available at https://gitlab.com/pypr/convergence_sph.

7 Discussions and Conclusions

In this paper, we have performed a numerical study of the accuracy and convergence of a variety of SPH schemes in the context of weakly-compressible fluids. Based on the numerical study performed in the previous sections, we summarize the key findings below.

7.1 Choice of smoothing kernel

We first considered the SPH approximation of a function and its derivative using different kernels. All the kernels considered here show second-order convergence when the support radius is high. The accuracy is marginally effected by the change in type of kernel. The smoothing error of an SPH approximation scales as and this necessitates that the smoothing length of the kernel be as small as possible. This implies that be small. As is well known, the discretization errors scale as and this necessitates that the smoothing radius be larger. These two requirements are contradictory. We find that by using a modest , along with the kernel corrections of Bonet and Lok (1999) or Liu and Liu (2006) we are able to obtain close to second-order convergence for the kernels considered in this work. It holds up to a resolution of , where which appears to be among the highest resolutions we have seen in the literature concerning the convergence of SPH methods. In the literature, we find kernels like the cubic spline to demonstrate pairing instabilities Dehnen and Aly (2012). We can avoid this instability by using a particle shifting technique (PST).

7.2 Choice of suitable operators

The SPH approximation of operators like the gradient, divergence, and Laplacian must be chosen carefully. In this paper, we recommend two methods for gradient approximation and three methods for viscous term approximation that ensure second-order convergence. The approximations which ensure pair-wise linear momentum conservation are always divergent. In the future, one could explore pair-wise linear momentum conserving and second-order convergent SPH approximation in a perturbed domain using SPH. Furthermore, the widely used artificial compressibility assumption makes the scheme accurate. We recommend using high values or a dual-time stepping criteria to achieve convergence. Solving the pressure using the pressure Poisson equation may also provide SOC, although those have not been studied in this work.

7.3 Particle density and fluid density

We recommend that one employ the fluid density in the governing differential equation as a property that convects with the particle. The approximation of the SPH operators should not be a function of a property of the fluid i.e. density. We obtain the integration volume by eq. 5 where the mass and kernel support radius of particles is kept constant. In the future, it would be important to explore the convergence of SPH operators when the mass as well as the support radius of the particles are varying as required by an adaptive SPH algorithm.

7.4 New scheme and variations

We demonstrate Eulerian as well as Lagrangian SPH schemes that are second-order convergent. We show that the Eulerian schemes captures the divergence accurately due to symmetry in the particle distribution resulting in better accuracy in pressure. We derive a pressure evolution equation using the continuity equation that resembles the EDAC SPH scheme in literature. We show that the PST step in the Lagrangian method can be replaced by a remeshing step which is another moment conserving regularization. However, remeshing is not stable in the presence of jumps in the properties as observed in the case of the Gresho vortex (See

E) and incompressible shear layer. The PST step can be included in the momentum equation resulting in the -SPH scheme. From the -SPH method, one can obtain the Eulerian form of WCSPH method by setting the shifting velocity to . All these schemes are SOC when we use a second-order convergent approximation for the operators. We show that even though the schemes are non-conservative in the absolute sense, approximate conservation also produce accurate results in the case of incompressible flows.

Thus, by a judicious choice of discretization, particle shifting, and a separation of the fluid and particle densities we have shown that second-order convergence is possible using the SPH method for weakly-compressible flows. We do observe that the SPH discretization of the divergence operator introduces errors for divergence-free fields which are noticeably absent in the case of an Eulerian method due to the symmetry of the particle distribution. This introduces significant errors into the pressure; it would be valuable to develop more accurate divergence operators for the Lagrangian case.

Given that the proposed schemes are second-order, it would be important to study the boundary conditions employed in the SPH to see how they affect the accuracy and order of convergence. In the future a similar analysis in the context of variable smoothing length, and mass would be very useful in light of many recent developments of an adaptive SPH method. One concern of note is the increased computational effort required to maintain second-order convergence and future developments in this area would be important for practical simulation using the SPH method.

Appendix A The kernel corrections

We obtain the kernel correction by Bonet and Lok (1999) such that the new kernel, satisfies,


where is the correction matrix for the kernel gradient, . The eq. 40 represents the first moment of the kernel gradient tends to unity as . The first order Taylor series of a function, defined at about is,


integrating both sides with


Using the SPH approximation of the integral on both the sides, we get


The eq. 43 is accurate Fatehi and Manzari (2011). Therefore, the kernel gradient correction by Bonet and Lok (1999) imply the use of asym form of the gradient approximation for an accurate approximation. On the other hand, the correction proposed by Liu and Liu (2006), originates by convolving the Taylor series by and , and solving all the equation simultaneously. The matrix form is given by


where, is the destination particle index, is the neighbor particle index, for is the kernel gradient component in the direction. All the terms containing are summed over all the neighbor particles. On solving the eq. 44, we obtain a first order consistent gradient Liu and Liu (2006). For a constant field, this method ensures that we satisfy and . Therefore, with this correction in both the sym1 and asym forms the second term becomes zero, and we get the SOC approximation. Whereas, in sym2 the term does not become zero, thus even this correction fails to correct the approximation.

Appendix B The divergence-free velocity field

In this section, we explain the reason behind the astonishing accuracy obtained when evaluating the divergence of a divergence-free field using a fixed Cartesian mesh of particles as discussed in section 4.2. We consider a 2-dimensional velocity field. We write the error , in the divergence evaluation as