A review on approaches to solving Poisson’s equation in projection-based meshless methods for modelling strongly nonlinear water waves

Abstract

Three meshless methods, including incompressible smooth particle hydrodynamic (ISPH), moving particle semi-implicit (MPS) and meshless local Petrov–Galerkin method based on Rankine source solution (MLPG_R) methods, are often employed to model nonlinear or violent water waves and their interaction with marine structures. They are all based on the projection procedure, in which solving Poisson’s equation about pressure at each time step is a major task. There are three different approaches to solving Poisson’s equation, i.e. (1) discretizing Laplacian directly by approximating the second-order derivatives, (2) transferring Poisson’s equation into a weak form containing only gradient of pressure and (3) transferring Poisson’s equation into a weak form that does not contain any derivatives of functions to be solved. The first approach is often adopted in ISPH and MPS, while the third one is implemented by the MLPG_R method. This paper attempts to review the most popular, though not all, approaches available in literature for solving the equation.

1 Introduction

Marine structures are widely used in ocean transportation, exploitation and exploration of offshore oil and gas, utilization of marine renewable energy and so on. All these are vulnerable to harsh weather and so to very violent waves. Under action of violent waves, they may suffer from serious damages. Therefore, it is crucial to be able to model the interaction between violent waves and structures for designing safe and cost-effective marine structures. The available numerical models for strongly nonlinear interactions between water waves and marine structures are mainly based on solving either the fully nonlinear potential flow theory (FNPT) or the Navier–Stokes (NS) equations. For dealing with the problems associated with violent waves, the NS model should be employed.

The NS model may be solved by either mesh-based methods or meshless methods. The former is usually based on the Eulerian formulation, but the latter on the Lagrangian formulation. In the meshless methods, the fluid particles are largely followed and so the methods are also referred to as particle methods. The mesh-based methods have been developed for several decades and mainly based on finite volume and finite different methods (Greaves 2010; Causon et al. 2010; Chen et al. 2010; Zhu et al. 2013). The meshless (or particle) methods are of relative new development, but have been recognized as promising alternative methods in recent years, particularly for modelling violent waves and their interaction with structures owing to their advantages that meshes are not required and numerical diffusion associated with convection terms is eliminated in contrast to mesh-based methods. Extensive review of all the methods would divert the focus of this paper. A brief overview for meshless methods is given below, as this paper is concerned only on topics related to them. For more information about mesh-based methods, the readers are referred to other publications, such as Causon et al. (2010) and Zhu et al. (2013).

1.1 Overview of meshless methods

Many meshless methods have been developed and reported in literature, such as the moving particle semi-implicit method (MPS) (e.g. Koshizuka 1996; Gotoh and Sakai 2006; Khayyer and Gotoh 2010), the smooth particle hydrodynamic method (SPH) (e.g. Monaghan 1994; Shao et al. 2006; Khayyer et al. 2008; Lind et al. 2012), the finite point method (e.g. Onate et al. 1996), the element free Galerkin method (e.g. Belytschko et al. 1994), the diffusion element method (Nayroles et al. 1992), the method of fundamental solution (e.g., Wu et al. 2006), the meshless local Petrov–Galerkin method based on Rankine source solution (MLPG_R method) (e.g. Ma 2008) and so on. Among them, the MPS, SPH and MLPG_R methods have been used to simulate violent wave problems.

When the meshless methods are applied to model strongly nonlinear or violent waves, two formulations are employed. One is based on the assumption that the fluid can be weakly compressed, while the other just assumes that the fluid is incompressible. The first one is mainly adopted for SPH, e.g. Monaghan (1994), Dalrymple and Rogers (2006), Gomez-gesteira et al. (2010) and so on. More references can be found in Violeau and Rogers (2016). The second formulation has been implemented in SPH, MPS and MLPG_R methods. The SPH based on incompressible assumption is called as incompressible smooth particle hydrodynamic (ISPH). Most of the publications that employ the three meshless methods for modelling incompressible flow are based on the projection scheme developed by Chorin (1968). One of the main tasks associated with the projection-based meshless methods is to find the pressure through solving Poisson’s equation. Various SPH methods have been reviewed very recently by Violeau and Rogers (2016). All the aspects of ISPH and MPS have also been discussed by Gotoh and Khayyer (2016). This paper tries to only review the approaches of solving Poisson’s equation in the meshless methods for incompressible flow.

1.2 Mathematical formulation of projection-based meshless methods

For completeness, the mathematical formulation and numerical procedure of projection-based meshless methods are summarized in this subsection. The incompressible Navier–Stokes equation (referred to as NS equation) and continuity equation together with proper boundary conditions including the free surface one are applied. In the fluid domain,

$$\begin{aligned} \nabla \cdot \vec {u}=0, \end{aligned}$$

(1)

$$\begin{aligned} \frac{D\vec {u}}{Dt}=-\frac{1}{\rho }\nabla p+\vec {g}+\upsilon \nabla ^2\vec {u}, \end{aligned}$$

(2)

where $\vec {g}$ is the gravitational acceleration; $\vec {u}$ is the fluid velocity; $\rho$ and $\upsilon$ are the density and the kinematic viscosity of fluid, respectively; and p is the pressure. On a rigid boundary, the velocity and pressure satisfy

$$\begin{aligned} \vec {n}\cdot \vec {u}=\vec {n}\cdot \vec {U}, \end{aligned}$$

(3a)

$$\begin{aligned} \vec {n}\cdot \nabla p=\rho ( {\vec {n}\cdot \vec {g}-\vec {n}\cdot \vec {\dot{U}}+\upsilon \,\vec {n}\cdot \nabla ^2\vec {u}}), \end{aligned}$$

(3b)

where $\vec {n}$ is the unit vector normal to the rigid boundary; $\vec {U}$ and $\vec {\dot{U}}$ are the velocity and acceleration of a rigid boundary which may be a part of the structures. When the structures are floating and freely responding to waves, $\vec {U}$ and $\vec {\dot{U}}$ need to be found by solving floating dynamic equations which will not be discussed here. On the free surface, a dynamic condition must be imposed, i.e.

$$\begin{aligned} p=0. \end{aligned}$$

(4)

If the multiphase flow is considered, the condition on the fluid–fluid interface needs to be considered. For more details, readers may refer to, e.g. Shao (2012) and Hu and Adams (2007, (2009).

In the projection-based meshless methods, the above equations are solved using the following time-split procedure (Chorin 1968). In the procedure, when or after the velocity, pressure and the location for each particle at $n\mathrm{th}$ time step ($t=t_{n})$ are known, one uses the following steps to find the corresponding variables at $(n+1)\mathrm{th}$ time step.

• (1) Calculate the intermediate velocity ($\vec {u}^*)$ and position ($\vec {r}^*)$ of particles using

  $$\begin{aligned} \vec {u}^*=\vec {u}^n+\vec {g}\Delta t+\upsilon \nabla ^2\vec {u}^n\Delta t, \end{aligned}$$

  (5)
  $$\begin{aligned} \vec {r}^*=\vec {r}^n+\vec {u}^*\Delta t, \end{aligned}$$

  (6)

  where $\vec {r}$ is the position vector of particles and the superscript n represents the n-th time step; $\Delta t$ is the increment of the time step.

• (2) Evaluate the pressure $p^{n+1}$ using

  $$\begin{aligned} {\nabla ^2}{p^{n + 1}} = \Lambda \frac{{{\rho ^{n + 1}} - \rho ^ *}}{{\Delta {t^2}}} + (1 - \Lambda )\frac{\rho }{{\Delta t}}\nabla .\mathop u\limits ^{ \rightarrow *}, \end{aligned}$$

  (7)

  where $\Lambda$ is a coefficient taking a value between 0 and 1. $\rho ^{n+1}$ and $\rho ^*$ are the fluid densities at $(n+1)$ time step and intermediate fluid density, respectively.

• (3) Calculate the fluid velocity and update the position of the particles using

  $$\begin{aligned} \vec {u}^{**}=-\frac{\Delta t}{\rho }\nabla p^{n+1}, \end{aligned}$$

  (8a)
  $$\begin{aligned} \vec {u}^{n+1}=u^*+\vec {u}^{**}=\vec {u}^*-\frac{\Delta t}{\rho }\nabla p^{n+1}, \end{aligned}$$

  (8b)
  $$\begin{aligned} \vec {r}^{n+1}=\vec {r}^n+( {1-\beta })\vec {u}^{n+1}\Delta t+\beta \vec {u}^n\Delta t, \end{aligned}$$

  (9)

  where $\beta$ is often taken as 0 (such as in Ma and Zhou 2009) or 0.5 (such as Cummins and Rudman 1999).

• (4) Go to (1) for the next time step.

The above procedure is followed by all the three meshless methods (ISPH, MPS and MLPG_R) discussed here. These methods are different in several aspects including estimating velocities and finding solution for pressure. In this paper, our attention is focused on discussing how different they are in solving Eq. (7) for pressure. That is because solving the equation dominates the cost of the computational time, and also because the accuracy of the solution for pressure determines the accuracy of the overall solution for wave dynamic problems.

It is noted that multiple sub-steps in each time step in the above procedure may be applied as in Hu and Adams (2007). No matter how many sub-steps are used, the solution of Poisson’s equation is always concerned in the projection-based meshless methods. Again, as our focus here is on the approaches for solving Poisson’s equation, readers who are interested in multiple sub-steps procedure may refer to relevant publications such as Hu and Adams (2007).

The right hand side of Eq. (7) is the source term of the Poisson’s equation combining the terms of density invariant and velocity divergence. The appropriate choice of the $\Lambda$ value has been discussed for achieving relatively more ordered particle distribution and more smoothing pressure field by, e.g. Ma and Zhou (2009), Gui et al. (2014, (2015). In addition, attempts are also made by improving the source term. Khayyer et al. (2009) replaced the source term by a higher-order source term, while Kondo and Koshizuka (2011), Khayyer and Gotoh (2011), Khayyer and Gotoh (2013), Gotoh and Khayyer (2016) and Gotoh et al. (2014) introduced an error-compensating term (including a high-order main term and two error-mitigating terms multiplied by dynamic coefficients). The higher-order source and the error-compensating terms help to enhance the pressure field calculation, volume conservation and uniform particle distributions throughout the simulation that minimizes the perturbations in particle motions. As the work related to improving the formulation and evaluation of the source term has been well covered by the cited papers, further details will not be given in this paper. This review hereafter focuses on the ways to deal with the Laplacian on the left hand side of Eq. (7).

Another issue in solving Eq. (7) is related to the boundary conditions satisfied by pressure on the free surface, rigid (fixed or moving) wall and arbitrary boundaries (also called in/outlets) introduced for computation purpose. To numerically implement the boundary condition on the rigid wall, several approaches have been suggested, including, for example, addition of dummy particles (e.g. Lo and Shao 2002; Gotoh and Sakai 2006) and unified semi-analytical wall boundary condition (e.g. Leroy et al. 2014). To numerically implement the boundary condition on the free surface, the key issue is how to identify the particles on it. There are several approaches for doing so, such as detecting if the density (particle number density) is smaller than a specified value (e.g. Lo and Shao 2002; Gotoh and Sakai 2006), mixed particle number density and auxiliary function method (Ma and Zhou 2009) and an auxiliary condition proposed by Khayyer et al. (2009). For implementing in/outlet conditions and other more details about the treatment of boundary conditions, the readers are referred to the recent review papers by Violeau and Rogers (2016) and Gotoh and Khayyer (2016).

2 Approaches of ISPH in solving Poisson’s equation

As far as we know, most publications based on the ISPH method adopt an approach that is to discretize Poisson’s equation directly. In such an approach, discretization of Laplacian is a key. Various different formulations of Laplacian discretization for the ISPH method are discussed in this section.

Use of the ISPH method appears to start in Cummins and Rudman (1999), which just gave the results for 2D problems irrelevant to water waves. In that paper, they employed the following approach to approximate the Laplacian in Eq. (7), i.e.

LP-SPH01:

$$\begin{aligned} \nabla \cdot {\left( {\frac{{\nabla p}}{\rho }} \right) _i} = {\sum \limits _j}\frac{{{m_j}}}{{{\rho _j}}}\frac{4}{{\rho i + \rho j}}\frac{{{p_{ij}}\mathop {\vec r}_{ij}.\nabla {\mathrm{{w}}_{ij}}}}{{{{\left| {{{{\vec r}_{ij}}}} \right| }^2} + {\eta ^2}}}, \end{aligned}$$

(10)

where $r_{ij} =\sqrt{( {x_i -x_j })^2+( {y_i -y_j })^2+( {z_i -z_j })^2}$, $p_{ij} =p_i -p_j$, $m_j$ is the mass of a particle and $\mathrm{W}_{ij} =W( {{\vec r}}_{ij})$ is the weight function or kernel function. One of the typical definitions for the kernel function is

$$\begin{aligned} w(| {{{\vec r}_{ij}}}|) = \left\{ \begin{array}{ll}1 - 6\left( {\frac{{\left| {{{\vec r}_{ij}}} \right| }}{h}} \right) ^{2} + 8{{\left( {\frac{{\left| {{{\vec r}_{ij}}} \right| }}{h}} \right) }^3} &{}\\ \quad - 3\left( {\frac{{\left| {{{\vec r}_{ij}}} \right| }}{h}}\right) ^4 &{} \quad 0 \le \,\,{{\frac{{\left| {{{\vec r}_{ij}}} \right| }}{h}}} \le 1\\ 0 &{} \quad {{\frac{{\left| {{{\vec r}_{ij}}} \right| }}{h} \le 1}} \end{array} \right. ,\nonumber \\ \end{aligned}$$

(11)

where h is the smoothing length. In the above equation, $\eta$ is a small number introduced to keep the denominator non-zero and often taken as h/10 (e.g. Lo and Shao 2002). LP-SPH01 may be obtained using the expression for estimating viscous stresses in Morris et al. (1997).

The above formulation was followed by Lo and Shao (2002), which considered the water waves propagating near shore. In their work, the Laplacian in Eq. (7) was approximated by

LP-SPH 02:

$$\begin{aligned} \nabla \cdot {\left( {\frac{{\nabla p}}{\rho }} \right) _i} = {\sum \limits _j}{m_j}\frac{8}{{{{(\rho _i + \rho j)}^2}}}\frac{{{p_{ij}}{r_{ij}}.\nabla {\mathrm{{w}}_{ij}}}}{{{{\left| {{{\vec r}_{ij}}} \right| }^2} + {\eta ^2}}}. \end{aligned}$$

(12)

This approximation was also employed by many other researchers, e.g. Shao et al. (2006), Rafiee et al. (2007), Ataie-Ashtiani and Shobeyri (2008), Ataie-Ashtiani et al. (2008) and Khayyer et al. (2008). It is noted here that the density in Eqs. (10) and (12) may be numerically estimated at the intermediate step even for incompressible fluids (e.g. Lo and Shao 2002; Gui et al. 2015). In these cases, $\rho _i \ne \rho _j$, though they may be quite close to each other except on the boundary. The density may also be specified as the physical density (e.g. Asai et al. 2012; Lind et al. 2012; Leroy et al. 2014) and thus should be the same at all particles, i.e. $\rho _i =\rho _j$. If it is, the above two expressions become exactly the same. Actually, if the fact is taken into account, both approximations are equivalent to

LP-SPH03:

$$\begin{aligned} {({\nabla ^2}p)_i} = \frac{2}{\rho }{\sum \limits _j}{m_j}\frac{{{p_{ij}}{{\vec r}_{ij}}.\nabla {\mathrm{{w}}_{ij}}}}{{{{\left| {{{\vec r}_{ij}}} \right| }^2} + {\eta ^2}}}, \end{aligned}$$

(13)

which was employed by Lee et al. (2008) and Xu et al. (2009). This one is the same as that given by Jubelgas et al. (2004) if $\eta =0$, which was derived using the idea employed by Brookshaw (1985). In their derivation, Jubelgas et al. (2004) expanded a function into a Taylor series ignoring all third- or higher-order terms. Following the same line, Schwaiger (2008) showed that

$$\begin{aligned} {\nabla ^2}p({r_i})\approx & {} {\Gamma _{\beta \gamma }}p{,_{\beta \gamma }}({\vec r_i}) \approx 2{\int _{\Omega } }(p(\vec r) \nonumber \\&- p(\vec r_i))\frac{{(\vec r_ij){\alpha ^W},a}}{{{{\left| {{\vec r}_{ij}} \right| }^2}}}\mathrm{d}\Omega \nonumber \\&- 2{\int _\Omega }\,W{,_\alpha }p({\vec r_i})\mathrm{d}\Omega , \end{aligned}$$

(14)

where the subscripts indicate the components of coordinates; the repeated subscripts such as $\alpha$ denote summation over it; $p_{,\alpha } ( {{\vec r}}_i)$ is the partial derivative of the function with respect to a coordinate and ${\Gamma _{\beta \gamma }}\, = {\int _\Omega }{({\vec r_{ij}})_\beta }{({\vec r_{ij}})_\gamma }\frac{{({{\vec r}_{ij}}){\alpha ^{w,\alpha }}}}{{{{\left| {{{\vec r}_{ij}}} \right| }^2}}}\mathrm{d}\Omega$. In his paper, the equation was given in terms of a general function. Here, it is written specifically for pressure to be consistent with other equations. As he indicated, if the weight function is symmetric and the support domain is entire, ${\Gamma }_{\beta \gamma } =\delta _{\beta \gamma }$ ($\delta _{\beta \gamma } =1\,\mathrm{if}\,\beta =\gamma ;\mathrm{otherwsie}\,\mathrm{zero})$ and ${\int _\Omega }\,W{,_\alpha }{f_{,\alpha }}({\vec r_i})\mathrm{d}\Omega = 0$, and so Eq. (14) will be the same as that given by Jubelgas et al. (2004) and its discretized form will be LP-SPH03 with $\eta =0$. Corresponding to his formulation, Schwaiger (2008) gave the following discrete Laplacian:

LP-SPH04:

$$\begin{aligned} {\nabla ^2}p({\vec r_i})\approx & {} \frac{{{\Gamma ^{{ - 1}}}\beta \beta }}{\kappa }\nonumber \\&\times \left\{ 2{\sum \nolimits _j}\frac{{{m_j}}}{{{p_j}}} (p(\vec r_j) - p(\vec r_i)) \frac{{(\vec r_ij){\alpha ^W},a}}{{{{\left| {{{\vec r}{_i}_j}} \right| }^2}}}\right. ,\nonumber \\&-\left. 2p,\alpha ({{\vec r}_i}).\left[ {{\sum \nolimits _j}\frac{{{m_j}}}{{{p_j}}}W,\alpha } \right] \phantom {\left\{ 2{\sum \nolimits _j}\frac{{{m_j}}}{{{p_j}}} (p(\vec r_j) - p(\vec r_i)) \frac{{(\vec r_ij){\alpha ^W},a}}{{{{\left| {{{\vec r}{_i}_j}} \right| }^2}}}\right. }\right\} , \end{aligned}$$

(15a)

$$\begin{aligned} {\nabla ^2}p({\vec r_i})\approx & {} {\sum \limits _j}\frac{{{m_j}}}{{{p_j}}}(p(\vec r_j) - p(\vec r_i)){C_{\alpha \beta }}W{,_\beta } \end{aligned}$$

(15b)

$$\begin{aligned} {C_{\alpha \beta }}= & {} {\left[ {{\int _\Omega }{{({{\vec r}_{ij}})}_\alpha }W{,_\beta }\,\,\mathrm{d}\Omega } \right] ^{ - 1}}, \end{aligned}$$

(15c)

where $\kappa$ is the number of dimensions, e.g. $\kappa$ = 2 for 2D cases. The sizes of ${\Gamma }_{\alpha \beta }^{-1}$ and $C_{\alpha \beta }$ are both 2 $\times$ 2 for 2D cases and 3 $\times$ 3 for 3D cases, and only the trace of ${\Gamma }_{\alpha \beta }^{-1}$ is required. Compared to others, this formulation requires inverse matrixes ${\Gamma }_{\beta \beta }^{-1}$ and $C_{\alpha \beta }$ at each particles and so may bear extra computational costs. It was employed for solving thermal diffusion problem without a free surface in Schwaiger (2008), but extended and tested by Lind et al. (2012) to solve water wave problems.

Hu and Adams (2007) suggested the following approximation by considering particle-averaged spatial derivative:

LP-SPH05:

$$\begin{aligned} \nabla \cdot {\left( {\frac{{\nabla p}}{\rho }} \right) _i} = 2{\sigma _i}{\sum \limits _j}\frac{1}{{\left| {{{\vec r}_{ij}}} \right| }}\frac{{\partial {\mathrm{{w}}_{ij}}}}{{\partial \left| {{{\vec r}_{ij}}} \right| }}\Bigg [ {\frac{1}{{\sigma _{i}^{2} }} + \frac{1}{{\sigma _{j}^{2}}}} \Bigg ]\frac{{{p_{ij}}}}{{\rho i + \rho j}},\nonumber \\ \end{aligned}$$

(16)

where $\sigma _i =\mathop \sum \nolimits _j \mathrm{W}_{ij}$ and $\rho _i =m_i \sigma _i$. If $\rho _i =\rho _j$ and $\sigma _i =\sigma _j$ with the entire support domain, LP-SPH05 becomes equivalent to LP-SPH03 (with $\eta =$ 0) as discussed above. That means that they are largely in the same order of accuracy.

Hu and Adams (2009) suggested another approximation with double summations:

LP-SPH06:

$$\begin{aligned} \nabla \cdot \left( {\frac{\nabla p}{\rho }}\right) _i= & {} -\frac{1}{\sigma _i }\mathop \sum \limits _j \nabla \text{ W }_{{ij}} \cdot \left[ \frac{1}{m_i }\mathop \sum \limits _k \nabla \text{ W }_{{ik}} \left( {\frac{p_i }{\sigma _i^2 }+\frac{p_k }{\sigma _k^2 }}\right) \right. \nonumber \\&\left. -\frac{1}{m_j }l \sum \nolimits _l \nabla \text{ W }_{jl} \left( {\frac{p_j }{\sigma _j^2 }+\frac{p_l }{\sigma _l^2 }}\right) \right] . \end{aligned}$$

(17)

In Hu and Adams (2009), LP-SPH06 was only used for the calculation of intermediate pressure for particle density correction. For full step velocity updating, they still used the pressure obtained by LP-SPH05. As far as we know, the use of LP-SPH06 has not been found in other publications. It is not clear if it can be employed alone.

Hosseini and Feng (2011) used the following approximation, which was derived to ensure that the gradient of a linear function is accurately evaluated as proposed by Oger et al. (2007).

LP-SPH07:

$$\begin{aligned} {({\nabla ^2}p)_i} = \frac{2}{{{\rho _{}}}}{\sum \limits _j}\,{m_j}\frac{{pi{j^{\varvec{r}}}ij.C({{\vec r}_{ij}}).\nabla {\mathrm{{W}}_{ij}}}}{{{{\left| {{{\vec r}_{ij}}} \right| }^2}}} \end{aligned}$$

(18a)

$$\begin{aligned} C(\vec r_{ij}) = \left( \begin{matrix} {{\sum {V_j}({x_j} - {x_i}){\textstyle {{\partial {W_{ij}}} \over {\partial x}}}\,\sum {V_j}({x_j} - {x_i}){\textstyle {{\partial {W_{ij}}} \over {\partial y}}}}}\\ {{\sum {V_j}({y_j} - {y_i}){\textstyle {{\partial {W_{ij}}} \over {\partial x}}}\,\sum {V_j}({y_j} - {y_i}){\textstyle {{\partial {W_{ij}}} \over {\partial y}}}}}\end{matrix}\right) ^{ - 1}.\nonumber \\ \end{aligned}$$

(18b)

where $V_j =\frac{m_j }{\rho _j }$ denotes the volume occupied by a particle. The expression of Matrix C here, the discrete form of Eq. (15c), is only for two-dimensional (2D) problems. For three-dimensional (3D) problems, it will be a 3 $\times$ 3 matrix (see, e.g. Schwaiger 2008). Khayyer et al. (2008) suggested a similar correction to the kernel gradient and applying it to estimating internal viscous force calculation to preserve both linear and angular momentum, but not for solving Poisson’s equation.

Gotoh et al. (2014) derived the following expression using the divergence of the pressure gradient,

LP-SPH08:

$$\begin{aligned} {({\nabla ^2}p)_i} = {\sum \limits _{j \ne i}}\frac{{{m_j}{p_{ij}}}}{{{\rho _j}}}\left( {\frac{1}{{\left| {{{\vec r}_{ij}}} \right| }}\frac{{\partial {w_{ij}}}}{{\partial \left| {{{\vec r}_{ij}}} \right| }} - \frac{{{\partial ^2}{w_{ij}}}}{{\partial {{\left| {{{\vec r}_{ij}}} \right| }^2}}}} \right) . \end{aligned}$$

(19)

The definition of $p_{ij}$ here is slightly different from that in Gotoh et al. (2014) and so the equation appears to be different but it is actually the same.

Apart from these forms described above, another discrete Laplacian was formed by Chen et al. (1999, (2001). In the scheme, all the second derivatives are found by solving the set of following equations:

$$\begin{aligned} {B}_{\eta \xi } F_\xi =\varPhi _\eta \end{aligned}$$

(20a)

$$\begin{aligned} {\Phi _\eta }= & {} {\sum \limits _j}\,{V_j}\,\left( {p({{\vec r}_j}) - p({{\vec r}_i})} \right) {W_{,\lambda \gamma }} \nonumber \\&- {p_{,\alpha }}({\vec r_i}){\sum \limits _j}\,{V_j}\,{({\vec r_{ij}})_\alpha }{W_{,\lambda \gamma }}, \end{aligned}$$

(20b)

$$\begin{aligned} {\mathrm{{B}}_\eta }_\xi = (1 - {\delta _{\alpha \beta }}/2){\sum \limits _j}\,{V_j}\,{({\vec r_{ij}})_\alpha }{({\vec r_{ij}})_\beta }{W_{,\lambda \gamma }}, \end{aligned}$$

(20c)

where $F_\xi =p_{,\alpha \beta } =\partial p/\partial r^\alpha \partial r^\beta$ with correspondence between $\xi$ and $\alpha \beta$ being $1\leftrightarrow 11,2\leftrightarrow 22,3\leftrightarrow 33$, $4\leftrightarrow 12$, $5\leftrightarrow 23$, $6\leftrightarrow 13$, as between subscripts $\eta$ and $\lambda \gamma$. After solving for all the second derivatives of pressure, the discrete Laplacian can be formed by summing up of $p_{,\alpha \alpha }$. This Laplacian discretization is named as corrective smoothed particle method (CSPM) following Schwaiger (2008). Fatehi and Manzari (2011) derived a new scheme (they called it as Scheme 4) using error analysis. Careful examination reveals that their new scheme is almost the same as the CSPM. The only difference is that $\text{ B }_{\eta \xi }$ in their scheme contains a correction to the leading error caused by approximation to the gradient. The correction may not play a very significant role if the gradient used in Eq. (20) is accurately estimated. Therefore, the scheme by Fatehi and Manzari (2011) may be considered as one with the similar accuracy as CSPM. As indicated by Chen et al. (1999), the solution of the equation theoretically gives the exact value of the second derivatives for any particle distribution if the pressure is a constant, linear or parabolic field and if $p_{,\alpha } ( {x_i })$ equals the exact value of the first derivatives of the pressure. None of other approximations (LP-SPH01 to LP-SPH08) have such a good property. In view of this fact, this formulation can be considered as the most accurate one among all those discussed above. However, when this approach is employed for solving Poisson’s equation about pressure, ${\Phi }_\eta$ is not evaluable, as it contains the pressure itself. One must form the matrix and work out its inversion before it is used for discretising Poisson’s equation. Clearly, it is the most time-consuming one (Schwaiger 2008), as it requires finding the inversion of two matrixes for every particle. One is ${B}_{\eta \xi }$, which is 3 $\times$ 3 for 2D cases and 6 $\times$ 6 for 3D cases, and the other is matrix C, which is 2 $\times$ 2 for 2D cases and 3 $\times$ 3 for 3D cases. In addition, the property of matrixes is sensitive to distribution of particles and to the number of particles falling in the region characterized by the smoothing length. More discussions about this will be given in the section about the patch tests.

3 Approaches of MPS in solving Poisson’s equation

Moving-particle semi-implicit (MPS) method was proposed by Koshizuka et al. (1995) and Koshizuka (1996). In this method, Poisson’s equation (Eq. 7) is solved also by directly approximating the Laplacian. In the cited papers, the Laplacian was approximated by

LP-MPS01:

$$\begin{aligned} {({\nabla ^2}p)_i} = \frac{{2\kappa }}{{{\lambda _0}{\sigma _0}}}{\sum \limits _{i \ne j}}({p_j} - {p_i}){w_{ij}}, \end{aligned}$$

(21)

where $\kappa$ is the same as before, $\sigma _0$ is the initial value of $\sigma _i ,\lambda _0$ are defined by

$$\begin{aligned} {\lambda _0} = {\int _V}\,{\left| {{{\vec r}_{ij}}} \right| ^2}\,W({\vec r_{ij}})\mathrm{d}V/{\int _V}\,W({\vec r_{ij}})\mathrm{d}V. \end{aligned}$$

To stabilize the pressure calculation, an improved Laplacian discretization was proposed by Khayyer and Gotoh (2010).

LP-MPS02:

$$\begin{aligned} {({\nabla ^2}p)_i} = \frac{1}{{{\sigma _0}}}{\sum \limits _{i \ne j}}{p_{ij}}\left( {\frac{1}{{\left| {{{\vec r}_{ij}}} \right| }}\frac{{\partial {w_{ij}}}}{{\partial \left| {{{\vec r}_{ij}}} \right| }} - \frac{{{\partial ^2}{w_{ij}}}}{{\partial {{\left| {{{\vec r}_{ij}}} \right| }^2}}}} \right) . \end{aligned}$$

(22)

This equation was derived by applying the same principle as that for LP-SPH08 and is applicable for 2D simulations. The extension to 3D has also been developed by Khayyer and Gotoh (2012), in which the first term in the bracket of Eq. (22) disappears. Here, there is no term of $\frac{m_j }{\rho _j }$ or $V_j$ in the summation like in the SPH formulations. However, if $\frac{m_j }{\rho _j }$ or $V_j$ is constant and taken as their initial value, i.e. $1/\sigma _0$, LP-SPH08 becomes exactly the same as LP-MPS02. As the authors of the cited papers indicated, the LP-MPS02 gave better results than LP-MPS01. Nevertheless, simple tests show that LP-MPS02 cannot give an exact value of Laplacian even if the pressure is a simple function like $p=x^{2}$ and the particles are uniformly distributed, while LP-MPS01 can give a right value in such a situation. This is probably because the normalization ($\lambda _0 )$ is conducted for LP-MPS01, but not for LP-MPS02 which leads to 0th order consistence of LP-MPS02 as pointed out by Tamai et al. (2016). If it is the case, it is easily rectified. An alternative formation was proposed by Ikari et al. (2015) with an aim to improve LP-MPS02 and given by

LP-MPS03:

$$\begin{aligned} {l} {\left\langle {{\nabla ^2}p} \right\rangle _i}= & {} \sum \limits _{j \ne i} {p_{ij}}\,\left\{ \frac{2}{{{{\left| {{{\vec r}_{ij}}} \right| }^3}}}\frac{{\partial {w_{ij}}}}{{\partial \left| {{{\vec r}_{ij}}} \right| }}\,({C_{11}}{x_{ij}^2} \right. \nonumber \\&+\,({C_{12}} + {C_{21}}){x_{ij}}{y_{ij}} + {C_{22}}{y_{ij}^2}) \nonumber \\&- \left[ {\left( {\frac{1}{{{{\left| {{{\vec r}_{ij}}} \right| }^2}}}\frac{{{\partial ^2}{w_{ij}}}}{{\partial {{\left| {{{\vec r}_{ij}}} \right| }^2}}} - \frac{1}{{{{\left| {{{\vec r}_{ij}}} \right| }^3}}}\frac{{\partial {w_{ij}}}}{{\partial \left| {{{\vec r}_{ij}}} \right| }}\,} \right) } \right. \nonumber \\&\times \left\{ {{C_{11}}{x_{ij}^2} + ({C_{12}} + {C_{21}}){x_{ij}}{y_{ij}} + {C_{22}}{y_{ij}^2}} \right\} \nonumber \\&+ \frac{1}{{\left| {{{\vec r}_{ij}}} \right| }}\frac{{\partial {w_{ij}}}}{{\partial \left| {{{\vec r}_{ij}}} \right| }}({C_{11}} + {C_{22}}) - \frac{{{x_{ij}}}}{{{{\left| {{{\vec r}_{ij}}} \right| }^3}}}\frac{{\partial {w_{ij}}}}{{\partial \left| {{{\vec r}_{ij}}} \right| }}\nonumber \\&\times ({C_{11}}{x_{ij}} + {C_{21}}{y_{ij}}) \nonumber \\&\left. {\left. { + \frac{1}{{{{\left| {{{\vec r}_{ij}}} \right| }^3}}}\frac{{\partial {w_{ij}}}}{{\partial \left| {{{\vec r}_{ij}}} \right| }}\,({C_{12}}{x_{ij}} + {C_{22}}{y_{ij}})} \right] } \right\} , \end{aligned}$$

(23a)

where

$$\begin{aligned} C = {\left( \begin{array}{l} - \frac{1}{{{\sigma _0}}}{\sum \limits _{j \ne i}}\frac{{\partial {w_{ij}}}}{{\partial \left| {{{\vec r}_{ij}}} \right| }}\frac{{{x^2}_{ij}}}{{\left| {{{\vec r}_{ij}}} \right| }} - \frac{1}{{{\sigma _0}}}{\sum \limits _{j \ne i}}\frac{{\partial {w_{ij}}}}{{\partial \left| {{{\vec r}_{ij}}} \right| }}\frac{{{x_{ij}}{y_{ij}}}}{{\left| {{{\vec r}_{ij}}} \right| }}\\ - \frac{1}{{{\sigma _0}}}{\sum \limits _{j \ne i}}\frac{{\partial {w_{ij}}}}{{\partial \left| {{{\vec r}_{ij}}} \right| }}\frac{{{x_{ij}}{y_{ij}}}}{{\left| {{{\vec r}_{ij}}} \right| }} - \frac{1}{{{\sigma _0}}}{\sum \limits _{j \ne i}}\frac{{\partial {w_{ij}}}}{{\partial \left| {{{\vec r}_{ij}}} \right| }}\frac{{{y^2}_{ij}}}{{\left| {{{\vec r}_{ij}}} \right| }} \end{array} \right) ^{ - 1}} \end{aligned}$$

(23b)

The expression of C in Eq. (23b) is written only for 2D problems, though it can be straightforwardly extended to 3D problems. Theoretically, LP-MPS03 should be reduced to LP-MPS02 if C is a unit matrix, but actually it is not. The reason is perhaps attributed to the approximation adopted when deriving the LP-MPS03. Interested readers can find more details about this from the cited papers. Very recently, Tamai et al. (2016) proposed another formation given by

LP-MPS04:

$$\begin{aligned} {\nabla ^2}p({\vec r_i})\approx & {} tr\left\{ 2{M^{ - 1}}{\sum \limits _j}\left( p({{\vec r}_j}) - p({{\vec r}_i})\right. \right. \nonumber \\&-\left. \left. p({{\vec r}_{ij}})_\alpha p{,_\alpha }({{\vec r}_i}) \right) {w_{ij}}{q_{ij}} \phantom {\left\{ 2{M^{ - 1}}{\sum \limits _j}\left( p({{\vec r}_j}) - p({{\vec r}_i})\right. \right. }\right\} , \end{aligned}$$

(24)

where $q_{ij}$ is a vector of order 3 for 2D cases and order 6 for 3D cases, requiring inversion of another matrix associated with first-order derivatives. M is a matrix based on $q_{ij},$ also with a size of 3 $\times$ 3 for 2D cases and 6 $\times$ 6 for 3D cases. For the detailed definition of the matrixes, readers are referred to the cited paper. The derivation of the formulation is analogous to CSPM (Eq. 20) and Fatehi and Manzari (2011). However, the content of matrix M is different from Eq. (20c), in that a term associated with the first derivative is involved in M, like in Fatehi and Manzari (2011).

Apart from these, Tamai and Koshizuka (2014) proposed a scheme based on a least square method, but Tamai et al. (2016) pointed out that this scheme needed inversion of a larger size matrix (an order of 5 for 2D cases and 9 for 3D cases) and also a larger support domain (or smoothing length) to keep the matrix invertible. More discussions can be found in the cited papers.

4 Patch tests on different discrete Laplacians

To investigate the behaviours of different forms of Laplacian discretization, a few papers carried out patch tests. In some patch tests, a Laplacian discretization is applied to estimate the value of Laplacian for a specified function, which is defined on a specified domain, giving the exact evaluation of the error. This section will summarize these tests available in published papers.

Schwaiger (2008) investigated several discrete Laplacians, including CSPM, LP-SPH04 and LP-SPH03 with $\eta =0$ (unless the $r_{ij} \sim 0$, a small value of $\eta$ in LP-SPH03 does not play a significant role). He considered functions of $x^m+y^m$ and $x^my^m$ defined on a 2D domain of $2<x<3$ and $2<y<3$ with m $=$ 2, 3, 4, 5 and 6, and calculated the value of discrete Laplacians, respectively. Two particle configurations were considered, one with uniform distribution at a distance S between particles and the other with random perturbation of $\left| \epsilon \right| \le 0.4S$ (i.e. the distance between particles is determined by S + $\epsilon )$ in x- and y-directions on the basis of uniform distribution. They found that for uniform particle distribution, the results of LP-SPH04 and CSPM were very similar at the interior particles away from boundaries with the error at a level of the machine error. LP-SPH03 can also give quite a good estimation at these particles, though its error is larger. However, at the particles close to boundaries, all approximations can produce large errors, though the relative errors of LP-SPH04 and CSPM are smaller in most cases. Even for them, the relative error close to the boundaries can reach the level of 70 % as observed in Fig. 1 of Schwaiger (2008).

Fig. 1

Variation of errors with changes of particle distances (originally presented in Zheng et al. 2014): a mean error; b maximum error

Fig. 2

Variation of errors with changes of randomness (originally presented in Zheng et al. 2014): a mean error; b maximum error

For irregular or disorderly particle distributions, the behaviour of discrete Laplacians also depend on how to choose the smoothing length. Schwaiger (2008) studied two options, one is h = 1.2S and the other is h $=$ $0.268\sqrt{S}$. Based on the relative errors in the region 2.25 $< x<$ 2.75 and 2.25 $< y<$ 2.75 without accounting for the particles near the boundaries, they found that for $h =$ 1.2S, both LP-SPH04 and CSPM did not show a fully convergent behaviour, but remained at a fairly constant relative error,reducing the average particle distance. They also found that LP-SPH03 became divergent, i.e. the relative error increasing with reducing the particle distance. CSPM should give converged results even for irregular particle distribution. The reason it did not do so is perhaps because there were no sufficient number of particles within the region of size $h = 1.2S,$ due to irregular shifting of particle positions, yielding that the property of matrixes involved in the CSPM became worse, and so leading to non-convergent results. For $h = 0.268\sqrt{S}$, they showed in Fig. 5 of their paper that all approximations exhibited convergent behaviour and that LP-SPH04 and CSPM results converged much faster, with the rate being near the second order, while LP-SPH03 results converged much slower with its convergent rate being less than first order. The reason for LP-SPH04 and CSPM results to be in second-order convergent rate in this case is perhaps because $0.268\sqrt{S}$ is much larger than 1.2S, and so there were always sufficient number of particles involved in the cases studied. Schwaiger (2008) mentioned that the smoothing length h was often set proportional to S, but the divergence behaviour corresponding to the case is perhaps troubling.

Lind et al. (2012) carried out similar investigations by comparing LP-SPH03 with LP-SPH04 for the functions of $x^m+y^m$ with $m =$ 1, 2 and 3 defined on the same domain as that by Schwaiger (2008). They just confirmed that the relative error of LP-SPH04 could reach 70 %, while LP-SPH03 yielded an error of 4000 % on the boundary. At the row next to the boundary, the relative error of LP-SPH04 reduced to 4 %, while that of LP-SPH03 remained to be 500 %.

Lind et al. (2012) also carried out investigations by solving the equation of $( {\nabla ^2p})_i =1$ for 1D problem with the boundary conditions of dp/dy $=$ 10 at y $=$ 0 and $p =$ 1 at $y =$ 1 using LP-SPH04. For this purpose, they employed both uniform and non-uniform particle configurations. The latter was produced by specifying different small random perturbations of $(\pm 0.1 \sim \pm 0.5)S$ to the particle distance for uniform distribution. They particularly indicated that the relative error of the solution became larger with increased random perturbation: 1.5 % corresponding to $(\pm 0.1)S,$ but 17 % to $(\pm 0.5)S$. They also demonstrated that the LP-SPH04 may lead to results with a convergent rate of 1.2–1.3 (less than 2 as shown for $h=0.268\sqrt{S}$ by Schwaiger 2008) with the particle shifting scheme to maintain the particle orderliness.

Zheng et al. (2014) performed similar tests, but used the function of f(x,y) $=$ cos(4$\pi x$ + 8$\pi$y) defined in the region of 2 $\le$ x $\le$ 3 and 2 $\le$ y $\le$ 3. This function is closer to the real pressure in water waves than $x^m+y^m$ and $x^my^m$. The discrete Laplacian they considered also included the LP-SPH03 and LP-SPH04. In their tests, the domain was first divided into small squared elements with $\Delta x=\Delta y=S$. The particles were then redistributed according to $\Delta {x}'\, \text{ and }\, \Delta {y}'$ determined by $S[1+k(Rn-0.5)]$, where Rn is a random number between 0 and 1.0 and different for $\Delta {x}'\, \text{ and }\, \Delta {y}'$, and k is a constant. Clearly, $k=0$ leads to regular distribution of particles. $k>0$ makes the distribution of particles irregular or disorderly. As k increases, the disorderliness increases. The accuracy of the Laplacian approximations is quantified in a similar way to that in Schwaiger (2008), by evaluating the average relative errors

$$\begin{aligned} \mathrm{Er}=\sqrt{\sum \limits _{i=1}^N {\left( \frac{\nabla ^2f_{i,c} -\nabla ^2f_{i,a} }{\nabla ^2f_{i,a,m} }\right) ^2V_i } } , \end{aligned}$$

(25a)

where $\nabla ^2f_{i,a}$ is the analytical value of Laplacian with $\nabla ^2f_{i,a,m}$ being its magnitude, e.g. $\nabla ^2f_{i,a,m} =80\pi ^2$, for $f(x,y) =$ cos(4 $\pi$ x + 8 $\pi$ y); and $\nabla ^2f_{i,c}$ is the values of discrete Laplacian. When estimating the error, only the particles within the region of $2.2\le x\le 2.8$ and $2.2\le y\le 2.8$ are considered as in Schwaiger (2008). The accuracy of the Laplacian approximations is also quantified by estimating their maximum relative errors given as

$$\begin{aligned} \mathrm{Er}_{\max } =\max \left( \left| {\frac{\nabla ^2f_{i,c} -\nabla ^2f_{i,a} }{\nabla ^2f_{i,a,m} }} \right| \right) , \quad i=1,2,3\ldots N. \end{aligned}$$

(25b)

In their tests, $S=0.1,0.05,0.02,0.0125,0.01,0.08$ and $k=0,0.2,0.4,0.8,1.0,1.2$ were considered. Some of their results are reproduced in Figs. 1, 2, and 3. Figure 1a presents the average relative errors for different values of S with a value of k being fixed to be 0.8, i.e. with the random shift up to ±0.4S, the same as that in Schwaiger (2008). From the figure, one can see that the average error of LP-SPH04 is consistently reduced with reduction of S. This trend is similar to the results of Schwaiger (2008) for a function of $x^m+y^m$ obtained using $h=0.268\sqrt{S} ,$ but different from those of Schwaiger (2008) obtained using $h=1.2S$ which is shown to be constant with the reduction of S in their papers. The reason is perhaps because the smooth length used in Zheng et al. (2014) was larger, though it was still proportional to S. The average errors of LP-SPH03 can increase with the reduction of S, which is a divergent behaviour. Figure 1b demonstrates that the maximum error of LP-SPH04 consistently decreases until $S=0.0125$ or Log$(S)\approx -1.9,$ but increase with increasing the resolution of the particles after that. In addition, the smallest value of the error is Log$(\mathrm{Er}_\mathrm{max})>-0.6$, corresponding to $\mathrm{Er}_\mathrm{max} =25~{\% }$, which is considerably larger than the average errors for the same case (Fig. 1a) and may be considered to be significant as the error occurs inside the domain. Again, the maximum error of LP-SPH03 shows a divergent behaviour when $S<0.05$ (Log$(S)<-1.3)$. It is noted that overall, the accuracy of numerical methods are controlled by the maximum error, and not the average error.

Figure 2 plots the average and maximum errors for different values of k with $S=0.01$. One can see from the figures that the errors of both approximations (LP-SPH03 and LP-SPH04) increase with the increase of k values, i.e. with particle being more disorderly, which is consistent with observation of Lind et al. (2012). Furthermore, the maximum error inside the domain can become very large, for example, Log$(\mathrm{Er}_\mathrm{max} )>-0.4$, corresponding to $\mathrm{Er}_\mathrm{max} >40~\%$, at $k>0.8$ even for LP-SPH04.

Fig. 3

The number of particles with an error larger than a certain values for $S=0.01$ and $k=1.2$ (originally presented in Zheng et al. 2014)

To demonstrate if there is a significant number of particles with a large error, Zheng et al. (2014) plotted a figure similar to Fig. 3. In this figure, the horizontal axis shows the different ranges of relative error, e.g. [20, 30 %], while the vertical axis shows the number of particles whose error lies in a range. For example, in the range of [20, 30 %], there are about 230 particles for the LP-SPH04. The relative error at each individual particle used in this figure is estimated by $\mathrm{Er}_i =\left| {\frac{\nabla ^2f_{i,c} -\nabla ^2f_{i,a} }{\nabla ^2f_{i,a,m} }} \right| \quad (i=1,2,3\ldots N)$. This figure demonstrates that a quite large relative error (>20 %) can happen at a considerable number of particles for the approximations even when they are applied to computing the Laplacian of the quite simple function, though the number for the LP-SPH04 is much smaller than for the LP-SPH03.

In most of the above tests (except for some cases in Schwaiger 2008), the value of S / h is fixed with the smoothing length varying and sometimes with different randomness. Quinlan et al. (2006) discussed the theoretical convergence of approximating the gradient of a function used in SPH. They showed that the error caused by numerical approximations to the gradient did not only depend on the smoothing length and randomness (non-uniformity), but also on the ratio S / h. Specifically speaking, the error increases with the larger randomness and can be proportional to 1 / h if S / h is not small enough, which is consistent with that observed in the above results. When S / h is small enough, the convergent behaviour of the approximation to the gradient can be improved. Graham and Hughes (2007) particularly investigated the behaviour of LP-SPH03 with $\eta =0$ by varying the value of S / h. They studied the pressure-driven flow between parallel plates with a constant pressure gradient with the diffusion term estimated by LP-SPH03 ($\eta$ = 0) for three values (1, 1/1.25 and 1/1.5) of S / h. They showed that the method was not convergent in several cases they studied and that random particle configurations could have a dramatic effect on the accuracy of the SPH approximations. More specially, the results are divergent if their random factor is larger than 0.25, and their best results are these obtained by using $S/h =$ 1/1.5 with the particles fixed, among which the error reduces at a rate less than first order when their random factor is relatively small. Fatehi and Manzari (2011) also carried out tests by varying S / h from 1/1.5 to 1/3.5 on a scheme similar to LP-SPH03 with $\eta =0$ and their new scheme which is similar to CSPM (discussed above) by solving a thermal diffusivity problem defined on a unit square $0\le x\le 1$ and $0\le y\le 1$, which has a similar equation to the problem with the zero pressure gradient considered by Graham and Hughes (2007) . In their tests, regular and irregular particle distributions were considered, and the relative errors of numerical results to the analytical ones at a time near steady state were presented in their paper. The random perturbation they employed was $\left| \epsilon \right| \le 0.05S$ or $\left| \epsilon \right| \le 0.1S$, much less than $\left| \epsilon \right| \le 0.4S$ used by Schwaiger (2008). Their results showed that the scheme LP-SPH03 with $\eta =0$ had a convergent rate of first order at the best, and that their new scheme similar to CSPM had a convergent rate of second order. However, they indicated that the scheme did not work when smoothing length was 1.5S, consistent with the analysis of Quinlan et al. (2006). This is perhaps because the number of neighbouring particle is not sufficient, which may make the matrixes involved in CSPM invertible. It is not sure if the convergent rate would maintain when the random perturbation is larger.

Gotoh et al. (2014) presented some convergent test results on LP-SPH08. For this purpose, the approximation was used together with their error-compensation term to simulate a pressure field caused by a modified gravitational acceleration. Their results showed that for irregular particle distributions (the initial distribution randomly altered and half of the fluid particles displaced by ±0.02S), the normalized root mean square error reduced with decrease of the initial particle distance, a convergent behaviour. According to the cited paper, the errors are 0.108, 0.068 and 0.065 corresponding to S $=$ 0.004, 0.003 and 0.002, respectively, which gives an average convergent rate at about 0.7, though it is about 1.6 from 0.004 to 0.003.

Ikari et al. (2015) tested the discrete Laplacians (LP-MPS02 and LP-MPS03) for the MPS method. Their results are summarized here. The first case they presented was about the computation of a pressure field due to a sinusoidal disturbance to gravitational acceleration. The particles were randomly shifted by $\pm 0.05S$ on the basis of uniform distribution. As they indicated, the results of LP-MPS03 were better than those of LP-MPS02. They also showed that there were some spurious fluctuations in the pressure time histories from LP-MPS02 on reducing the particle distance. Their second case was similar to their first case except for a difference that the sinusoidal disturbance was multiplied by an exponential growing factor. The results for this case also showed the outperformance of LP-MPS03 compared to LP-MPS02. They pointed out that the clear convergence of results from LP-MPS03 was not observed in terms of root mean square error of numerical results relative to the analytical solution for the case. The third case they investigated was about a 2D diffusion problem on a square domain. For this case, the performances of both LP-MPS02 and LP-MPS03 were satisfactory, though LP-MPS03 was slightly better. The convergent trend was not, however, exhibited. For example, the root mean square error of LP-MPS03 is 19.2850, 30.8089 and 24.2292 corresponding to the mean particle distance of 10, 5 and 2.5 mm. The convergent property of LP-MPS02 with a higher-order source term on the right hand side of Eq. (7) was also examined by Khayyer and Gotoh (2012), showing an improved and more stabilized (without fast fluctuation) pressure for the similar case (but in 3D here) to that in Gotoh et al. (2014) discussed above. In this test, the initial distribution of particles is randomly altered and half of the fluid particles are displaced by ${\mp } 0.05 S$, similar to that in Ikari et al. (2015). The results demonstrated that the normalized root mean square error reduced with decrease of the initial particle distance. In other words, convergent behaviour was observed. The specific information is that the errors are 0.241, 0.228 and 0.192 corresponding to S $=$ 0.012, 0.010 and 0.008, respectively. The average convergent rate is near 0.8.

Tamai et al. (2016) carried out tests using the discrete Laplacians to estimate the values of Laplacian for a given function on a square domain $0\le x \le 1 \text { and } 0 \le y \le 1$. This is similar to the method used by Schwaiger (2008) and Zheng et al. (2014) discussed above, but Tamai et al. (2016) used a sum of four exponential functions. In their tests, the random distribution of particles was also achieved by randomly shifting the particle position on the basis of uniform distribution. The random disturbance was given by a normal distribution with zero expectation and standard deviation of 0.1. The degree of randomness was higher than Gotoh et al. (2014) and Khayyer and Gotoh (2012), in the same level as in Fatehi and Manzari (2011), but not as large as in Schwaiger (2008). The smoothing length they used was not less than 2.7S, quite large compared to the tests mentioned above. The results of the tests in Tamai et al. (2016) indicated that (a) the maximum errors of LP-SPH03, LP-MPS01 and LP-MPS02 grew with reduction of mean particle distance (S), i.e. showing a divergent behaviour, similar to Fig. 1b given by Zheng et al. (2014) and (b) the convergent rate of LP- MPS04 is about 2, which is similar to the observation on CSPM by Schwaiger (2008) and Fatehi and Manzari (2011).

Table 1 Summary of the main features of different discrete Laplacians

In summary, the above tests are clearly not extensive to cover all Laplacian approximations, but they indeed cover some best approaches available so far in literature. Their main features and typical behaviours observed in the tests described above are outlined in Table 1. In the table, the approximations are classified into three types for the convenience of discussion here. Type 1 includes those without the need of matrix inversion, such as LP-SPH03, Type 2 includes those with one matrix inversion, while Type 3 refers to those with the need of two matrix inversions. According to the results, one may find that the schemes can be improved in the following aspects.

• The discretization of Laplacians can be a notable issue at particles near a boundary, especially for disordered particle distributions, such as water surface, without additional appropriate treatment. This may not be a big issue in some applications where the solution near the boundary is not mainly concerned, but would be a critical issue for modelling water waves and their interaction with structures in marine or coastal engineering, in which the accuracy of pressure near the water and body surface is important. The corrected Laplacian operator in integral formulation (Souto-Iglesias 2013) improves the Laplacian evaluation near the boundary and gives convergent solutions of Poisson’s equation with boundary conditions which is also applied to evaluate the curvature in Khayyer et al. (2014). This may suggest that the discretization schemes of Laplacians discussed above should be employed together with the correction to improve their behaviour near boundaries.

• The error of discrete Laplacians can become larger when the degree of particle disorderliness is higher or results converge slower even inside computational domains. In the cases for violent water waves, the particle distribution always becomes highly disordered even though they are uniformly and regularly located initially. More effort may be required to make them less sensitive to particle disorderliness.

• It is observed that Type 1 Laplacian approximations may not converge for a high degree of particle distribution randomness (or disorderliness), but they may converge at a rate less than first order for a low degree of particle distribution randomness (or disorderliness). That means that the results obtained from approximations may become worse with reduction of particle distance or increase of the number of particles used when particle distribution randomness level is high. Type 2 may have similar problem, though it may be more accurate for the same number of particles.

• Type 3 has a convergent rate of 2nd order if the random level of particle distribution is not very high, but the rate may become lower with the increase of the random level. The computation costs of the type are high compared with others. In addition, the number of neighbouring particles must always be high enough to ensure the matrixes to be invertible. This is not necessarily guaranteed when modelling violent water waves as the configuration of particles can dramatically and dynamically vary during simulation, which is not a priori predictable.

These issues associated with larger errors and divergent behaviours observed in the tests on Type 1 and 2 Laplacian approximations as indicated in Table 1 do not mean that the methods using the approximations could not give acceptable results. Actually, a large body of literature has provided abundant evidence that the ISPH (e.g. Lo and Shao 2002; Lind et al. 2012) and MPS (e.g. Khayyer and Gotoh 2010, 2011, 2012, 2013) based on the schemes are successful in many applications. These issues mainly imply that researchers need to make much effort in determining the right number of particles to be used to achieve acceptable results. If their number of particles is not right, their solutions bear a large error. It may also imply that to achieve results with a specified accuracy requires a large number of particles.

These results demonstrate that although a lot of effort has been undertaken to develop better approximations to the Laplacian discretization, it still needs to be improved, in particular for modelling violent water waves where the water particles can become severely disordered. To improve the accuracy of computation, one may adjust the distribution of particles as done by Xu et al. (2009) to reduce the disorderliness of particles. The other way to circumvent the issues is by avoiding the direct discretization of Laplacian. This approach will be discussed in the next section.

5 Approaches of MLPG_R in solving Poisson’s equation

Bonet and Kulasegaram (2000, (2002) adopted a variational formulation to solve Possion’s equation for their ISPH method. In their basic formulation, only the gradient was involved, and second-order derivatives were avoided entirely. The discrete gradient in the formulation was involved and so the inversion of a matrix had to be performed for each of all particles, similar to that in LP-SPH07 and LP-MPS03. They, however, showed that the analytical gradient of a linear function with a constant gradient could not be guaranteed to be correctly evaluated even with a correction in their basic variational formulation. They then introduced an integration correction factor that is actually a vector. The factor was indeed effective to overcome the problem associated with basic variational formulations. However, to estimate the factor, they needed iterations which involve summation and multiplication of matrixes at each of all particles, requiring significant extra computational costs. Another problem with their variational formulation arises also from evaluating function gradients, so-called ‘spurious modes’, i.e. the nonzero gradient of a function being perhaps estimated numerically as zero, which exists even with the use of the integration correction factor. To eliminate the spurious modes, they introduced a least-square stabilization method by adding a stabilization potential to the variational formulation. The stabilization procedure requires the evaluation of the Laplacian (i.e. the second-order derivatives) of the solution function. As they discussed, a special correction had to be applied to ensure that the evaluation of Laplacian was correct for linear or quadratic functions. One can see that the use of the variational formulation does not only require extra computational costs, but also actually still require dealing with the Laplacian. That would be perhaps the reason why it has not become popular in the SPH community.

Fig. 4

Illustration of integration and support domains for MLPG_R method

Ma (2005a) started to employ the MLPG method to solve water wave problems. In this method, the Poisson’s equation (Eq. 7 with $\Lambda =0$ without loss of generality in this paper) is first integrated over a circle in 2D (Fig. 4) and a sphere in 3D cases to give

$$\begin{aligned} \int \limits _{\Omega _I } {\left[ {\nabla ^2p-\frac{\rho }{\Delta t}\nabla \cdot \vec {u}^{( *)}} \right] \varphi \mathrm{d}\Omega } =0, \end{aligned}$$

(26)

where $\Omega _I$ is the integration domain centred at node I and $\varphi$ is a test function, which can be arbitrarily chosen. Ma (2005a) employed a Heaviside step function as the test function and arrived at MLPGR-01:

$$\begin{aligned} \int \limits _{\partial \Omega _I } {\vec {n}\cdot \nabla p\mathrm{d}S} =\frac{\rho }{\Delta t}\int \limits _{\partial \Omega _I } {\vec {n}\cdot \vec {u}^*\mathrm{d}S}, \end{aligned}$$

(27)

where $\partial \Omega _I$ is the boundary of $\Omega _I$ and $\vec {n}$ is the normal vector of $\partial \Omega _I ,$ pointing out of the integration domain. In this formulation, the pressure and its gradient are estimated by

$$\begin{aligned} p( {\vec {x}})\approx \sum \limits _{J=1}^N {\Phi _J ( {\vec {x}}){\hat{p}}_J } \,\,\, \text {and} \,\,\, \nabla p( {\vec {x}})\approx \sum \limits _{J=1}^N {\nabla \Phi _J ( {\vec {x}}){\hat{p}}_J }, \end{aligned}$$

(28)

where $\Phi _J ( {\vec {x}})$ is a shape function formulated by moving the least square (MLS) method and using a local weight function defined on a support domain (Fig. 4), similar to the one used in SPH or MPS, e.g. Eq. (11); ${\hat{p}}_J$ is the presentative pressure, which may not be equal to the real pressure at the point concerned. This formulation is similar to that of Bonet and Kulasegaram (2000) in the sense that only gradients, no second-order derivatives, are directly involved.

This formulation was soon enhanced into the MLPG method based on rankine source solution (shortened as MLPG_R method) by Ma (2005b). In the MLPG_R method, the solution of the Rankine source is taken as the test function, i.e. the function $\varphi$ satisfies $\nabla ^2\varphi =0$ in $\Omega _I$ except for the centre and $\varphi =0$ on $\partial \Omega _I$ with a radius of $R_I$. The expression of the solution for Rankine source is proposed to be

$$\begin{aligned} \varphi =\frac{1}{4\pi }(1-R_I /r) \text {for 3D cases} \end{aligned}$$

(29a)

$$\begin{aligned} \varphi =\frac{1}{2\pi }\ln (r/R_I )\text {for 2D cases}, \end{aligned}$$

(29b)

where r is the distance between the concerned point and the centre of $\Omega _I$. Based on this test function, the integration of Eq. (26) can be changed into

MLPGR-02:

$$\begin{aligned} \int \limits _{\partial \Omega _I } {\vec {n}\cdot ( {p\nabla \varphi })\mathrm{d}S-R_I ^{\kappa -2}p_I } =\int \limits _{\Omega _I } {\frac{\rho }{\Delta t}\vec {u}^{( *)}\cdot \nabla \varphi \, \mathrm{d}\Omega }, \end{aligned}$$

(30)

where $\kappa$ is the number of dimensions, $\kappa$ = 2 for 2D cases and $\kappa$ = 3 for 3D cases, as defined in previous sections. The MLPGR-02 formulation is similar to MLPGR-01 in the sense that only the boundary integral on the left hand side is involved. However, it is distinct from the latter and also from the variational formulation of Bonet and Kulasegaram (2000) in the sense that it is the pressure, rather than pressure gradient, that is dealt with in the MLPGR-02. Therefore, there are no issues associated with evaluating the pressure gradient, such as integration correction factor and spurious modes as discussed by Bonet and Kulasegaram (2000). Apart from this, there is of course no issue related to Laplacian discretization.

Fig. 5

Variation of mean error with changes of particle distance

It is noted that the right hand side of MLPGR-02 is a volumetric integration, i.e. the integration domain is a sphere in 3D and a circle in 2D cases, rather than a surface integration as in MLPGR-01. Use of normal numerical integration techniques, like Gaussian quadrature, may need a significant amount of computational time for estimating the term. To overcome this problem, Ma (2005b) and Zhou et al. (2010) developed semi-analytical techniques for 2D and 3D cases. With use of these techniques, the computational costs for evaluating the right hand term in MLPGR-02 is similar to that for evaluating the right hand term in MLPGR-01.

It is also noted that a special interpolation technique was developed by Ma (2008) for discretizing pressure and velocity in Eq. (30). For completeness, it is given below in terms of the pressure.

$$\begin{aligned} p( {\vec {r}_0 })=\sum \limits _{J=1}^N {\Phi _J ( {\vec {r}_0 ;\vec {r}_I })} p( {\vec {r}_J }), \end{aligned}$$

(31a)

$$\begin{aligned} \Phi _J ( {\vec {r}_0 ;\vec {r}_I })= & {} \frac{w( {\left| {\vec {r}_J -\vec {r}_0 } \right| })}{\sum \nolimits _J^N {w( {\left| {\vec {r}_J -\vec {r}_0 } \right| })} }-( {1-\delta _{IJ} })B_{0,J} ( {\vec {r}_I }) \nonumber \\&+\, \delta _{IJ} \sum \limits _{J\ne I}^N {B_{0,J} ( {\vec {r}_I })}, \end{aligned}$$

(31b)

$$\begin{aligned} B_{0,J} ( {\vec {r}_I })=\vec {R}_0 \cdot \vec {A}_{IJ}, \end{aligned}$$

(31c)

$$\begin{aligned} \vec {A}_{IJ}= & {} \left[ \frac{n_{I,y} B_{IJ,x} -\bar{n}_{I,xy} B_{IJ,y} }{n_{I,x} n_{I,y} -\bar{n}_{I,xy} ^2},\nonumber \right. \\&\left. =\frac{n_{I,x} B_{IJ,y} -\bar{n}_{I,xy} B_{IJ,x} }{n_{I,x} n_{I,y} -\bar{n}_{I,xy} ^2} \right] ^T, \end{aligned}$$

(31d)

$$\begin{aligned} \vec {R}_0 =\frac{\sum \nolimits _J^N {( {\vec {r}_J -\vec {r}_0 })w( {\left| {\vec {r}_J -\vec {r}_0 } \right| })} }{\sum \nolimits _J^N {w( {\left| {\vec {r}_J -\vec {r}_0 } \right| })} }, \end{aligned}$$

(31e)

$$\begin{aligned} B_{IJ,x_m } =\frac{( {\vec {r}_{J,x_m } -\vec {r}_{I,x_m } })}{\left| {\vec {r}_J -\vec {r}_I } \right| ^2}w( {\left| {\vec {r}_J -\vec {r}_I }, \right| }) \end{aligned}$$

(31f)

$$\begin{aligned} n_{I,x_m }= & {} \sum \limits _{J=1,J\ne I}^N {\frac{( {\vec {r}_{J,x_m } -\vec {r}_{I,x_m } })^2}{\left| {\vec {r}_J -\vec {r}_I } \right| ^2}W( {\left| {\vec {r}_J -\vec {r}_I } \right| })} \nonumber \\&\quad (m=1\,\,\mathrm{or}\,\,2), \end{aligned}$$

(31g)

$$\begin{aligned} \bar{n}_{I,xy}= & {} \sum \limits _{J=1,J\ne I}^N {\frac{( {\vec {r}_{J,x_m } -\vec {r}_{I,x_m } })( {\vec {r}_{J,x_k } -\vec {r}_{I,x_k } })}{\left| {\vec {r}_J -\vec {r}_0 } \right| ^2}w( {\left| {\vec {r}_J -\vec {r}_I } \right| })} \nonumber \\&\quad (m=1\,\,\mathrm{or}\,\,2,k=1\,\,\mathrm{or}\,\,2,m\ne k), \end{aligned}$$

(31h)

Fig. 6

Sketch of sloshing tank ($d=0.5$ m, $L=2$ d , $a_0 =0.001\,d)$

Fig. 7

Wave time histories and convergent behaviour of different methods, originally presented in Zheng et al. (2014). a Wave time histories on the left wall obtained by different methods, b error of numerical results

$$\begin{aligned} \delta _{IJ} =\left\{ {\begin{array}{l} 1 \qquad I=J \\ 0 \qquad I\ne J \\ \end{array}} \right. , \end{aligned}$$

(31i)

where $\vec {r}_{I,x_m }$ is the component of $\vec {r}\,\,$ in $x_m$ (m $=$ 1, 2, or 3) direction. The set of equations is obtained by replacing Eq. (12) in Ma (2008) for estimating the gradient with Eq. (15) for 2D cases in that paper. It can be straightforwardly extended to 3D cases using Eq. (15) for 3D cases in Ma (2008). When working out the gradient, one needs inversion of a matrix which has a size of 2 $\times$ 2 for 2D cases and 3 $\times$ 3 for 3D cases similar to that for Type 2 discrete Laplacians. In other words, the MLPG_R scheme based on Eq. (31) in this paper needs about the same level of computational costs of that based on Type 2 discrete Laplacians, but the inversion of such a higher order matrix will be much more efficient than Type 3 discrete Laplacians as the latter requires the version of two matrixes, the larger one with a size of 3 $\times$ 3 for 2D cases and 6 $\times$ 6 for 3D cases, respectively.

6 Patch tests of MLPGR-02 and comparative studies

Two cases are considered in this section. One is that the MLPGR-02 scheme is applied to solve Poisson’s Equation about a simple problem defined by $\nabla ^2p=0$ with $p(0,y)=0, p(1,y)=0, p(x,0)=0, p(x,1)=\mathrm{sin}(\pi x)$, solely for this paper. The analytical solution for this case is $p(x,y)=\mathrm{sin}h(\pi y)\mathrm{sin}(\pi x)$. It is the steady-state solution of the case studied by Fatehi and Manzari (2011). To solve the case, the domain is first divided into small squared elements ($\Delta x\times \Delta y$ with $\Delta x=\Delta y=S)$. The particles are then redistributed according to $\Delta {x}'\, \text{ and }\, \Delta {y}'$ determined by $[1+k\vartheta ]S$, where $\vartheta$ is a random number between −0.5 $\sim$ 0.5 and different for $\Delta {x}'\, \text{ and }\, \Delta {y}'$, and k is a constant factor, in the same way as in Zheng et al. (2014). Clearly, as k increases, the disorderliness becomes higher as discussed by Zheng et al. (2014). Particularly, the random factor $\left| \epsilon \right| \le 0.1S$ employed by Fatehi and Manzari (2011) corresponds to $k=0.2$.

Fig. 8

CPU times used by three numerical methods corresponding to different number of particles, originally presented in Zheng et al. (2014)

Fig. 9

Pressure distribution of violent sloshing at two instants of time [CISPH (upper row) and ISPH_R (lower row), originally presented in Zheng et al. (2014)]. a $\tilde{t}$ = 20.8, b $\tilde{t}$ = 29.6

Fig. 10

Comparison of pressure time histories obtained by using different methods a SPH, b CISPH (CISPH2 in Zheng et al. 2014) and c ISPH_R, originally presented in Zheng et al. (2014)

The errors estimated by $\mathrm{Er}_p ={\sqrt{\sum \nolimits _{i=1}^{N_t } {\left| {p_i -p_{i,a} } \right| ^2} } } / {\sqrt{\sum \nolimits _{i=1}^{N_t } {\left| {p_{i,a} } \right| ^2} } }$ for $k=0.2$ and 0.3 are shown in Fig. 5, where $p_i$ is the numerical results, while $p_{i,a}$ is the analytical solution. From this figure, one can see that the convergent rate is close to the second order. Compared with the methods using Type 1 and Type 2 Laplacians which converge at a lower rate, the convergent properties of the MLPGR-02 scheme is much better. Compared with these adopting the Type 3 Laplacian, the computational efficiency of the MLPGR-02 scheme is higher as the inversion of a matrix with a size of 3x3 for 2D cases and 6x6 for 3D cases is required by the Type 3 Laplacian, but the inversion of such a higher order matrix is not needed in the MLPGR-02 scheme.

The second case was considered by Zheng et al. (2014) who developed a hybrid method. In the hybrid method, pressure is solved using MLPGR-02 and all others are the same as ISPH. They named the hybrid method as incompressible smoothed particle hydrodynamics based on Rankin source solution (ISPH_R) method. They applied the ISPH_R method together with the ISPH method based on LP-SPH04 (named as CISPH) and with the traditional weakly compressible SPH (named as SPH) to simulate the sloshing waves in a tank (Fig. 6) that is subject to the motion $X_s =a_0 \Omega /\sqrt{gd} (1-\cos \Omega t)$. The error of the numerical solution against the analytical solution is evaluated by $Er_\eta ={\sqrt{\sum \nolimits _{i=1}^{N_t } {\left| {\eta _i -\eta _{i,a} } \right| ^2} } } / {\sqrt{\sum \nolimits _{i=1}^{N_t } {\left| {\eta _{i,a} } \right| ^2} } },$ where $\eta _i$ is the numerical result at the time instant, $N_t$ is the total time steps in the simulation duration of $\tilde{t}=t\sqrt{L/g} =50.0$, and $\eta _{i,a}$ is the analytical solution at $i\mathrm{th}$ time step, given by Faltinsen (1976).

To simulate the case, the particles are uniformly distributed at start, but can become disorderly during simulation as they are moving together with waves, though the disordered level is not very high as the motion of waves is not very big in this case. Figure 7a depicts the wave time histories on the left wall obtained by different methods and compared with the analytical solution of Faltinsen (1976). It can be seen from Fig. 7a that the results of the SPH has a good agreement with the analytical solution at the first three periods, but in the later stage numerical dissipation in the wave amplitude becomes evident. In addition, in the time range (such as t/ $\surd$ g/$L=$ 35–40) of short and small waves, the traditional SPH method cannot catch the details correctly. In contrast, the results from the CISPH and ISPH_R method can well catch the details of short and small waves and do not show visible dissipation. Figure 7b presents the errors of numerical results from the SPH, CISPH and ISPH_R methods relative to the analytical solution corresponding to different numbers of particles employed, originally presented in Zheng et al. (2014). It clearly shows that the convergent rate of results from the ISPH_R method is about second order and the error is considerably smaller than those of SPH and CISPH methods. In other words, to achieve any specified accuracy, the ISPH_R method needs much less number of particles (or larger particle sizes) than others. For example, corresponding to $\mathrm{Log}(\mathrm{Er}_\eta )\approx -3.55$, the particle size required by ISPH_R and CISPH are $\mathrm{Log}(S)\approx$ −1.70 and −2.08, or $S=0.02$ and 0.008, respectively. In addition, the ISPH_R can lead to a low level of error, such as $\mathrm{Log}(\mathrm{Er}_\eta )=-4,$ but CISPH cannot yield a result with such a low error.

Fig. 11

Qualitative comparison of pressure field obtained by CISPH-HS (left), experiment (center) and CISPH-HS-HL-ECS (right) at t $=$ 0.1, 0.2, 0.3 and 0.4 $T_{0}$, originally presented in Gotoh et al. (2014)

Fig. 12

Time histories of pressure obtained by CISPH-HS, experiment and CISPH-HS-HL-ECS, originally presented in Gotoh et al. (2014)

To explore the properties of the methods in another way, Fig. 8 depicts the CPU time spent by all the methods corresponding to different numerical errors on the same computer, which is also originally presented in Zheng et al. (2014). One can see from Fig. 8 that the ISPH_R method needs much less CPU time to achieve the same level of accuracy. For example, corresponding to $\mathrm{Log}(\mathrm{Er}_\eta )\approx -3.545$, the CPU times spent by the ISPH_R and CISPH are $\mathrm{Log}(\mathrm{CPU}\_t)\approx$ 2.97 and 4.02, corresponding to $\mathrm{CPU}\_t\approx 925$ and 10, 531 (about 11 times of the former) seconds, respectively. It is noted that the CPU time for running a case may depend on the choice of solver and preconditioner for solving the system of linear algebraic equations resulting from the discretized equations. As far as we know, the results in Fig. 8 were obtained by a solver combining the GMRES with Gauss–Seidel method. At each time step, they firstly run the Gauss–Seidel procedure for a specified number of iterations and then run the GMRES if necessary. If different procedure would have been used, the CPU time would be different.

To further show the performance of different methods in cases involving violent waves, Fig. 9 presents the pressure distribution at two time instants for sloshing in a rectangular tank shown in Fig. 6, but with the parameters of $L =$ 0.6 m, $d = 0.12$ m $= 0.2$ L, $a_0 (\text{ moiton } \text{ amplitude })=0.05$ m and $T_0 (\text{ moiton } \text{ peirod })=1.5$ s. Figure 10 depicts the pressure time histories recorded at a point on the left wall with a height of 0.1667L, resulting from three approaches (SPH, ISPH_R and CISPH defined above). The results are also compared with Kishev and Kashiwagi (2006).

Gotoh et al. (2014) presented some results for the same cases as in Figs. 9 and 10, which are produced using CISPH-HS and CISPH-HS-HL-ECS. LP-SPH02 together with a higher order source term was used in CISPH-HS. LP-SPH08 was employed in CISPH-HS-HL-ECS together with the error-compensating source (ECS) term. Their original figures for pressure fields and pressure time histories at the same point as in Fig. 10 are duplicated in Figs. 11 and 12, respectively. As they indicated, the pressure trace by CISPH-HS is characterized by frequent and, relatively, large-amplitude unphysical oscillations, while the results of CISPH-HS-HL-ECS are much smoother. If comparing the results in Figs. 10 and 12, one may find that the pressure time histories produced by ISPH_R and CISPH-HS-HL-ECS have a similar level of smoothness and agreement with the experimental data. This indicated that the ECS term is quite effective, because there are still visible unphysical oscillations if it is not applied as shown in Fig. 4 of Gotoh et al. (2014). It may be interesting to see more comparisons of different approaches to improve their performance further.

7 Conclusions

This paper has reviewed the approaches to solve Poisson’s equation for pressure involved in incompressible smoothed particle hydrodynamic (ISPH), moving particle semi-implicit (MPS) and meshless local Petrov–Galerkin method, based on Rankine source solution (MLPG_R) methods for simulating nonlinear or violent water waves assuming fluids are the incompressible. As summarized in Table 2, there are three different approaches, i.e. discretizing Laplacian directly (DLD) by approximating the second-order derivatives, transferring Poisson’s equation into a weak form containing only gradient of pressure (WCG) and transferring Poisson’s equation into a weak form that does not contain any derivatives of functions to be solved (WCF). The first approach DLD has been employed by most publications related to ISPH and MPS, while the third approach WCF has been implemented by the MLPG_R method for modelling water waves.

For effectively implementing the first approach DLD, three types of discrete Laplacians have been proposed as summarized in Table 1. Type 1 does not need inversion of matrix and so is relatively computationally efficient. However, the patch tests available have shown that this type of discrete Laplacians may converge at a rate less than first order for random (or disorderly) particle distribution. Type 3 is the most accurate one, but it requires the inversion of two matrices (one of them with a size of 3 $\times$ 3 for 2D cases and 6 $\times$ 6 for 3D cases) for each particle and so is relatively computationally inefficient. Patch tests discussed above have demonstrated that the convergent rate of this type can reach to second order.

Table 2 Three approaches for solving Poisson’s equation

Type 1 discrete Laplacians in the DLD approach is relatively easier to implement and have been the most popular one so far, particularly in the community which employs the ISPH and MPS methods. Their performance can be improved by applying the error-compensating term on the right hand side of Poisson’s equation, by reducing the randomness of particle distribution and by adding a correction term near the boundaries. More efforts may be made to improve their performance by enhancing the convergent rate.

The third approach (WCF) adopted by the MLPG_R method does not need to deal with any derivative and so has no issue related to discretizing derivatives when solving Poisson’s equation for pressure. Using relatively simple approximation (requiring only inversion of one matrix with a size of 2 $\times$ 2 for 2D cases and 3 $\times$ 3 for 3D cases) to the pressure, one can achieve the second-order convergent rate in patch pests, as that achieved by using Type 3 discrete Laplacians in the first approach, and also in solving the sloshing waves with small amplitudes. In addition, limited tests in literature available demonstrate that the ISPH based on the third approach (WCF) requires less CPU time to achieve the results with the same accuracy compared to ISPH based on the first approach (DLD) for simulating sloshing waves. Nevertheless, the WCF approach is currently less popular than others, perhaps because it is relatively new as it has just started to be used since 2005.

More comparative studies are encouraged, in particular for applying all the schemes to the same cases. The studies should compare convergent rate, accuracy and computational efficiency of different approaches. Such studies help fully understand the behaviours of the different approaches and select the best for modelling water waves in general cases.