An accurate and stable alternating directional moving particle semi-implicit method for incompressible flow simulation

Abstract

As a Lagrangian mesh free method, Moving Particle Semi-implicit (MPS) method can easily handle complex incompressible flow with free surface. However, some deficiencies of MPS method such as inaccurate results, unphysical pressure oscillation and particle thrust near free surface still required to be further resolved. In this paper, an improved MPS method based on Crank-Nicolson time scheme, named Alternating Direction Moving Particle Semi-implicit (ADMPS) method, is proposed. In addition, a time integral source term (TIS) is developed to suppress unphysical pressure oscillation and a free-surface revise (FR) technique is proposed to avoid particle thrust near free surface. Enhancement of accuracy and stability by these improvements is proved by hydrostatic and dam breaking benchmarks, and the numerical results of the ADMPS-FR-TIS model show good agreement with the hydrostatic analytical solution and the dam breaking experiment data.

Introduction

The Moving Particle Semi-implicit (MPS) method was firstly proposed by Koshizuka et al. ¹ for simulating incompressible flow with free surface. As a mesh-free Lagrangian particle method, the MPS method can easily handle with simulating problems with complex flow, such as ship-wave interaction in rough seas, ² hemodynamics of small vessels,^3,4 tsunami-structure interaction, ⁵ hydrodynamic of multiphase flow,^6,7 liquid sloshing,^8,9 bubble rising and coalescence,^10,11 flow and solidification behavior of molten metals, ¹² rotary atomization, ¹³ fluid induced deformation and destruction. ¹⁴ Despite the success in simulating complex flows, however, the MPS method still faces the challenges of accuracy and stability.

In order to alleviate pressure oscillation, several kinds of improved source term of Poison pressure equation (PPE) have been proposed. Khayyer and Gotoh^15,16 proposed a derivative form of source term which is called the Higher order Source (HS) term. For reducing the pressure oscillation due to particle number density, mixed DI/DF terms by adopting weighted form to combined DI (Density Invariant) and DF (Divergence-Free) terms were proposed.^17–19 In order to enhance the volume conservation, some multi-terms with error compensating parts were developed to keep the relative positions of particles close to the initial particle spacing.^20–23 Considering the compressibility of fluid, weak compressible source term ¹⁵ and quasi-compressibility source term ²⁴ were proposed. Duan et al. ¹² proposed a mixed DI/DF compressible term. Guo et al. ²⁵ proposed a compressible multi-term with compensating parts.

The original MPS pressure gradient operator is simple, and has several deficiencies such as lack of accuracy, non-conservation of momentum and insufficient of repulsive force. In order to enhance repulsive force and avoid particle aggregation, Koshizuka and Oka ²⁶ proposed a repulsive gradient operator by replacing the pressure of the target particle with the minimum pressure in the neighbor particle. Toyota et al. ²⁷ proposed a conservative gradient model to guarantee momentum conservation. Combined the advantages of the above two gradient operators, Khayyer and Gotoh ¹⁵ developed a conservative and repulsive gradient operator, which is similar to that in the CSPH method developed by Bonet and Lok. ²⁸ In view of the lack of accuracy of gradient operators, Khayyer and Gotoh ²⁰ proposed a correction gradient operator to ensure the consistency of the first order Taylor expansion. Duan et al. ¹² proposed a repulsive correction gradient operator, ensure both the accuracy and sufficient repulsive force. Wang et al. ²⁹ proposed another different repulsive correction gradient operator by directly introducing a corrective matrix into the repulsive correction gradient operator. By adding a dummy particle correction term into the gradient operator proposed by Wang et al., ²⁹ Liu et al. ³⁰ developed a high-order stabilized gradient model for further enhancing stability of MPS calculation.

One way to enhancement of MPS pressure calculation is to apply a more accurate Laplacian operator for PPE. By calculating the divergence of original gradient operator with the standard kernel function, Khayyer and Gotoh ³¹ obtained the high order Laplacian (HL) operator. Chen et al. ³² proposed a Laplacian model with virtual particle. Xu and Jin ³³ developed a new Laplacian model (abbr. as LMN) by calculating the divergence of original gradient operator with a kernel function which is different from the standard one. Liu et al. ³⁰ proposed Laplacian models with dummy particle. Liu et al. ³⁴ further developed a renormalized Laplacian model with dummy particles.

The judgment of free surface is very important for the accuracy and stability of pressure calculation. Once the internal particles are misjudged as free surface particles, the calculated pressure field may be seriously distorted, leading excessive particle displacement and causing particle thrust. The original method proposed by Koshizuka et al. ¹ judged free-surface particles just by their number density. Since particles are not always uniformly distributed, only judging free surface by particle number density might cause misjudge. Therefore, Tanaka and Masunaga ²⁴ made judgment by the number of neighboring particles. These simple one-judge methods are very simple but easy to cause misjudge. In order to further reduce misjudge, hybrid strategies were developed, such as twice-judge methods proposed by Lee et al., ¹⁸  Pan et al., ³⁵ and Zhang and Wan. ³⁶ There are also some methods to judge the free surface by geometric arrangement of neighboring particles, such as methods proposed by Sun et al., ¹⁹ Shibata et al., ³⁷ and Sun et al. ³⁸

In traditional MPS method, partial or differential equations and boundary conditions are discretized directly without using any integral form, and finally a set of strong-form discrete equations are obtained. However, strong-form mesh free method could be lack of accuracy and instability compared to the weak-form one. Therefore, a weak-form scheme called hybrid FVM-MPS method were developed recently by Liu et al. ³⁹

Despite the success of the efforts for improving the MPS method, few of them considered the effect of intermediate velocity, one of the most important factors contributing to pressure calculation. In MPS method, the intermediate velocity obtained by the explicit step has no prediction component, which may affect the accuracy of the pressure calculation. Therefore, in this paper, we proposed an improved MPS framework based on the Crank-Nicolson time scheme, called “Alternating Direction Moving Particle Semi-implicit (ADMPS) method,” to enhance the pressure calculation of the MPS method. In this method, each time step is divided into two equal parts, a directional step and a semi-implicit step. In the directional step, particle velocity is predicted explicitly by the pressure Poisson equation of incompressible flow, which is obtain directly under the velocity divergence free condition; while in the semi-implicit step, MPS algorithm is adopt to recover particle number density and calculate the displacement and velocity of next time step. In addition, in order to further suppress unphysical pressure oscillation, a time integral source term of pressure Poison equation is developed.

In MPS method, Dirichlet boundary condition is applied directly to free surface particles. According to the original gradient operator ²⁶ and improved gradient operators,^15,27,31,40 there is no interacted repulsion force between free-surface particles no matter how close they are. Thus, particle thrust might happen on the free surface, and would finally cause instability of MPS calculation. In order to prevent particle thrust on the free surface, modified PPEs and pressure gradient operators were proposed by an assumption that virtual (dummy) particles were over the free surface.^30,32,34,37 Under the virtual particle assumption, interacted repulsion force between free-surface particles could be exist in the calculation, and finally prevent particle thrust happened on the free surface. However, under this assumption, pressure of free surface particles could no longer be uniform (e.g. all equal to 0 Pa), and the initial particle number density on free surface is difficult to define. Hence, in the paper, a revise model without virtual particle assumption for free surface particles is also developed for preventing particle thrust.

After the Introduction, MPS method is simply introduced in Section 2. Then, the proposed improvements, including the ADMPS algorithm, the time integral source term (TIS) and the free-surface revise (FR) technique, are presented in Section 3. In Section 4, verification analyses are conducted to prove enhancement of accuracy and stability of the proposed measurements. Main conclusions are pointed out in Section 5.

MPS method

Governing equations

The governing equations of MPS method include the momentum equation (equation (1)) and the mass conservation equation (equation (2)).

$$\frac{d\overrightarrow{u}}{dt}=-\frac{1}{\rho }\nabla P+\nu {\nabla }^2\overrightarrow{u}+\overrightarrow{g}$$

(1)

$$\frac{d\rho }{\mathrm{dt}}+\rho \nabla ·\overrightarrow{u}=0$$

(2)

In equations (1) and (2), $\overrightarrow{u}$ , $t$ , $\rho$ , $P$ , $\nu$ , $\overrightarrow{g}$ represent velocity vector, density, pressure, viscous coefficient and acceleration of gravity respectively.

Gradient operator and Laplacian operator

The MPS gradient operator can be represented as:

$$<\nabla f{>}_i=\frac{d}{n_0}{\sum }_{j\ne i}[\frac{f_j-f_i}{{r_{\mathrm{ij}}}^2}{\stackrel{\rightharpoonup }{r}}_{\mathrm{ij}}w(r_{\mathrm{ij}})]$$

(3)

where, $f$ denotes the physical quantity; $d$ represents dimension of the calculation; $n_0$ denotes initial particle number density; ${\stackrel{\rightharpoonup }{r}}_{\mathrm{ij}}$ represents the displacement of particle $i$ to particle $j$ . $w(r)$ represents kernel function. In standard MPS, $w\left(r\right)$ follows the below equation:

$$w\left(r\right)=\{\begin{array}{cc}\frac{r_e}{r}-1,\hfill & r\le r_e;\hfill \\ \multicolumn{1}{c}{0,}& \mathrm{others}.\hfill \end{array}$$

(4)

where, $\mathrm{re}$ represents the searching radius. Usually, $r_e=2.1l_0$ ( $l_0$ denotes the particle resolution).

Particle number density follows below equation:

$$n={\sum }_{j\ne i}w\left(r_{\mathrm{ij}}\right)$$

(5)

When calculating the pressure gradient, the following pressure gradient (equation (6)) is usually used in order to prevent particle aggregation in MPS methods:

$$<\nabla P{>}_i=\frac{d}{n_0}{\sum }_{j\ne i}[\frac{P_j-{\hat{P}}_i}{{r_{\mathrm{ij}}}^2}{\stackrel{\rightharpoonup }{r}}_{\mathrm{ij}}w(r_{\mathrm{ij}})]$$

(6)

where ${\hat{P}}_i$ is the minimum pressure in the search domain of particle $i$ .

Laplacian operator follows the below expression:

$$<{\nabla }^{2}f{>}_i=\frac{2d}{{\lambda }_0{n}_0}{\sum }_{j\ne i}(f_j-f_i)w(r_{\mathrm{ij}})$$

(7)

$${\lambda }_0=\frac{{\sum }_{j\ne i}{{r_{\mathrm{ij}}}^0}^{2}w({r_{\mathrm{ij}}}^0)}{{\sum }_{j\ne i}w({r_{\mathrm{ij}}}^0)}$$

(8)

where, the superscript “0” indicates the standard arrangement of particles.

Algorithm of MPS method for incompressible flow

Normally, the MPS algorithm consists of two steps (an explicit step and an implicit correction step), and belongs to the operator-splitting algorithm.

In the explicit step, an intermediate velocity ${\overrightarrow{\mathrm{u}}}^{*}$ can be obtained by equation (9):

$$\frac{{\overrightarrow{u}}^{*}-{\overrightarrow{u}}^n}{\mathrm{\Delta }t}=\nu {\nabla }^2{\overrightarrow{u}}^n+\overrightarrow{g}$$

(9)

Then in the implicit step, velocity of the $(\mathrm{n}+1{)}^{\mathrm{th}}$ time step ${\overrightarrow{\mathrm{u}}}^{\mathrm{n}+1}$ can be calculated by equation (10):

$$\frac{{\overrightarrow{u}}^{n+1}-{\overrightarrow{u}}^{*}}{\mathrm{\Delta }t}=-\frac{1}{\rho }\nabla P^{n+1}$$

(10)

Added equations (9) and (10), the momentum equation is obtained. Taking the divergence of both sides of equation (10),

$${\nabla }^2{P}^{n+1}=\frac{\rho }{\mathrm{\Delta }t}(\nabla ·{\overrightarrow{\mathrm{u}}}^{*}-\nabla ·{\overrightarrow{\mathrm{u}}}^{\mathrm{n}+1})$$

(11)

For incompressible flows, ${\overrightarrow{\mathrm{u}}}^{\mathrm{n}+1}$ satisfies the incompressible continuity equation:

$$\nabla ·{\overrightarrow{u}}^{n+1}=0$$

(12)

Substitute equation (12) into equation (11), a DF-type PPE is obtained:

$${\nabla }^2{P}^{n+1}=\frac{\rho }{\mathrm{\Delta }t}\nabla ·{\overrightarrow{u}}^{*}$$

(13)

It can be found in equation (13) that, the value of the source term is determined by $\nabla ·{\overrightarrow{u}}^{*}$ when $\rho$ and $\mathrm{\Delta }t$ are known before. If the support domain of the calculation of $\nabla ·{\overrightarrow{u}}^{*}$ includes only fluid particles, according to equation (9),

$$\begin{array}{c}\multicolumn{1}{c}{\nabla ·{\overrightarrow{u}}^{*}=\nabla ·({\overrightarrow{u}}^n+(\nu {\nabla }^2{\overrightarrow{u}}^n+\overrightarrow{g})\mathrm{\Delta }t)}\\ =\nabla ·{\overrightarrow{u}}^n+\nu \nabla ·\left({\nabla }^2{\overrightarrow{u}}^n\mathrm{\Delta }t\right)\hfill \end{array}$$

(14)

It should be noted that, because the intermediate velocity of wall boundary particles is not determined by equation (9), the equation (14) is not applicable to the case where the support domain contains wall boundary particles.

If the viscous term is ignored in equation (14), then

$$\frac{\rho }{\mathrm{\Delta }t}\nabla ·{\overrightarrow{u}}^{*}\approx \frac{\rho }{\mathrm{\Delta }t}\nabla ·{\overrightarrow{u}}^n$$

(15)

It can be seen from equation (15) that in the MPS method, the value of the source term is almost determined by the velocity divergence at the previous time, and does not contain any predicted component.

According to equation (13) and the mass conservation equation of equation (2),

$$\frac{\rho -{\rho }^{*}}{\rho \mathrm{\Delta }t}=\nabla ·{\overrightarrow{u}}^{*}$$

(16)

where ${\rho }^{*}$ denotes intermediate density.

Replacing $\rho$ and ${\rho }^{*}$ by $n_0$ and intermediate number density $n^{*}$ in the left side of equation (16):

$$\frac{n^{*}-n_0}{n_0·t}=-\nabla ·{\overrightarrow{u}}^{*}$$

(17)

Simultaneous equation (13) and (17), a DI-type PPE is obtained (equation (18)):

$${\nabla }^2{P}^{n+1}=-\frac{\rho }{\mathrm{\Delta }{t}^2}\frac{n^{*}-n_0}{n_0}$$

(18)

Directly substituted density by particle number density can lead to errors and cause pressure oscillation. Thus, an empirical parameter $\gamma$ was introduced to equation (18) for suppressing pressure oscillation and improving accuracy. ³⁷  Then, a modified pressure Poison Equation (equation (19)) is obtained.

$${\nabla }^2{P}^{n+1}=-\gamma \frac{\rho }{\mathrm{\Delta }{t}^2}\frac{n^{*}-n_0}{n_0}$$

(19)

Shibata et al. ³⁷ gave a guideline for setting $\gamma$ (as equation (20) shows). In equation (20), $n^{**}$ is the permissible upper limit of intermediate number density, and ${|\frac{1}{\rho }\nabla P|}_{\max }$ denotes the expected maximum acceleration due to pressure gradient term on wall boundaries. Generally speaking, $\gamma =0.01$ is recommended.

$$1\ge \gamma \ge {\left|\frac{n^{**}-n_0}{n_0}\right|}^{-1}\frac{{\left(\mathrm{\Delta }t\right)}^{2}d}{\sqrt{{\lambda }_0}}{|\frac{1}{\rho }\nabla P|}_{\max }>0$$

(20)

ADMPS method with time integral source term

In this section, the proposed ADMPS method is presented. First, the PPE of incompressible flow is introduced in Section 3.1. Then the ADMPS algorithm and a new improve time integral source term are presented in Section 3.2 and 3.3 respectively. In Section 3.4, we expound the reason why particle thrust happens on free surface, and derive a revised model to suppress particle thrust. The free surface particle detection algorithm proposed by Pan et al. ³⁵ is presented in section 3.5.

The PPE of continuous incompressible flow

For incompressible flow, the continuity equation can be expressed by

$$\nabla ·\overrightarrow{u}\equiv 0$$

(21)

Dispersing both sides of the momentum equation (equation (1)):

$$\nabla ·\left(\frac{d\overrightarrow{u}}{\mathrm{dt}}\right)=-\frac{1}{\rho }\nabla ·(\nabla P)+\nu \nabla ·\left({\nabla }^2\overrightarrow{u}\right)+\nabla ·\left(\overrightarrow{g}\right)$$

(22)

For 2D problems, let $\overrightarrow{u}={(u,v)}^T$ , where $u$ and $v$ denote the horizontal and vertical components of $\overrightarrow{u}$ respectively. In equation (22),

$$\begin{array}{c}\multicolumn{1}{c}{\nabla ·\left(\frac{d\overrightarrow{u}}{\mathrm{dt}}\right)=\frac{\partial (\nabla ·\overrightarrow{u})}{\partial t}+{\left(\frac{\partial u}{\partial x}\right)}^2+{\left(\frac{\partial v}{\partial y}\right)}^2}\\ +2\frac{\partial v}{\partial x}\frac{\partial u}{\partial y}+u\frac{\partial (\nabla ·\overrightarrow{u})}{\partial x}+v\frac{\partial (\nabla ·\overrightarrow{u})}{\partial y}\hfill \end{array}$$

(23)

$$\nabla ·\left({\nabla }^2\overrightarrow{u}\right)={\nabla }^2(\nabla ·\overrightarrow{u})$$

(24)

$$\nabla ·\left(\overrightarrow{g}\right)=0$$

(25)

Substitute equation (21) into equations (23) and (24) respectively,

$$\nabla ·\left(\frac{d\overrightarrow{u}}{\mathrm{dt}}\right)={\left(\frac{\partial u}{\partial x}\right)}^2+{\left(\frac{\partial v}{\partial y}\right)}^2+2\frac{\partial v}{\partial x}\frac{\partial u}{\partial y}$$

(26)

$$\nabla ·\left({\nabla }^2\overrightarrow{u}\right)=0$$

(27)

Substitute equations (25)–(27) into equation (22), the pressure Poison Equation of incompressible flow (equation (23)) is obtained.

$${\nabla }^{2}P=\mathrm{SQ}$$

(28)

where,

$$\mathrm{SQ}=-\rho {[{\left(\frac{\partial u}{\partial x}\right)}^2+{\left(\frac{\partial v}{\partial y}\right)}^2+2\frac{\partial v}{\partial x}\frac{\partial u}{\partial y}]}^n$$

When the motion is hindered, fluid element will be compressed, resulting in velocity divergence not being 0 or the density changing. However, under the assumption of incompressible flow, the compressed fluid element will then recover instantaneously, that is, the velocity divergence will return to 0 and the density will return to the initial value instantaneously. Because equation (28) is derived under the divergence free condition, it can be adopted for predicting velocity explicitly.

ADMPS algorithm

In this part, the ADMPS algorithm based on the Crank-Nicolson time scheme is introduced. The Crank-Nicolson time scheme had been adopted to construct a MPS algorithm for highly viscous fluid by Sun et al. ⁴¹ Though the ADMPS algorithm and the MPS algorithm proposed by Sun et al. are both based on the Crank-Nicolson time scheme, the two algorithms are different. The MPS algorithm proposed by Sun et al. is constructed by applying the Crank-Nicolson scheme to the viscous term. While the ADMPS algorithm is constructed by applying the Crank-Nicolson time scheme to the entire momentum equation directly. The following is the construction process of the ADMPS method.

After applying the Crank-Nicolson time scheme to the Momentum Equation, equations (29) and (30) are obtained. The calculation of ${\overrightarrow{u}}^{\prime }$ through equation (29) is called the directional step, while the calculation of ${\overrightarrow{u}}^{n+1}$ through equation (30) is called the semi-implicit step.

$${\overrightarrow{u}}^{\prime }-{\overrightarrow{u}}^n=(-\frac{1}{\rho }\nabla P^{\prime }+\nu {\nabla }^2{\overrightarrow{u}}^n+\overrightarrow{g})\mathrm{\Delta }t/2$$

(29)

$${\overrightarrow{u}}^{n+1}-{\overrightarrow{u}}^{\prime }=(-\frac{1}{\rho }\nabla P^{n+1}+\nu {\nabla }^2{\overrightarrow{u}}^{\prime }+\overrightarrow{g})\mathrm{\Delta }t/2$$

(30)

In equation (29), $P^{\prime }$ is obtained by equation (31).

$${\nabla }^2{P}^{\prime }=S{Q}^{\prime }$$

(31)

where,

$$S{Q}^{\prime }=\{\begin{array}{cc}0,\hfill & \mathrm{SQ}>0;\hfill \\ \multicolumn{1}{c}{0,}& \mathrm{SQ}<0\mathrm{and}\nabla ·{\overrightarrow{u}}^n>0;\hfill \\ \multicolumn{1}{c}{\mathrm{SQ},}& \mathrm{others}.\hfill \end{array}$$

Equation (31) is a modified form of equation (28). SQ > 0 indicates that the fluid is under stretching. However, it is particularly notorious that fluid is not stretchable. Hence, when $\mathrm{SQ}>0$ , we set the value of the source term to 0. Besides, if $\mathrm{SQ}<0$ when $\nabla ·{\overrightarrow{u}}^n>0$ , velocity divergence would further increase. In this situation, the source term is reset to zero for reducing numerical instability.

In the semi-implicit step, the equation (30) is divided into equation (32) and equation (33):

$$\frac{{\overrightarrow{u}}^{*}-{\overrightarrow{u}}^{\prime }}{\mathrm{\Delta }t/2}=\nu {\nabla }^2{\overrightarrow{u}}^{\prime }+\overrightarrow{g}$$

(32)

$$\frac{{\overrightarrow{u}}^{n+1}-{\overrightarrow{u}}^{*}}{\mathrm{\Delta }t/2}=-\frac{1}{\rho }\nabla P^{n+1}$$

(33)

In equation (32), ${\overrightarrow{u}}^{*}$ denotes the temporary velocity. And in equation (33), pressure $P^{n+1}$ can be calculated by using the PPE of MPS method (equation (19)). The entire simulation process is shown in Figure 1.

Figure 1.

Simulation process of ADMPS method.

Combining equations (32) and (29),

$$\begin{array}{c}\multicolumn{1}{c}{{\overrightarrow{u}}^{*}={\overrightarrow{u}}^n+(-\frac{1}{\rho }\nabla P^{\prime }+\nu {\nabla }^2{\overrightarrow{u}}^n+\overrightarrow{g})\mathrm{\Delta }t/2}\\ +(\nu {\nabla }^2{\overrightarrow{u}}^{\prime }+\overrightarrow{g})\mathrm{\Delta }t/2\hfill \end{array}$$

(34)

If the viscous term is ignored, then

$$\frac{\rho }{\mathrm{\Delta }t}\nabla ·{\overrightarrow{u}}^{*}\approx \frac{\rho }{\mathrm{\Delta }t}\nabla ·{\overrightarrow{u}}^n-\frac{1}{2}\nabla ·(\nabla P^{\prime })$$

(35)

From equations (15) and (35), we can find that the source term of the ADMPS method contains an extra prediction term (the second term in the right side of equation (35)) comparing to that of the MPS method.

Improved time integral source term

In equation (19), the factor γ is introduced into the source term to alleviate pressure oscillation. Essentially, the pressure oscillation can be alleviated because γ reduces the absolute value of the source term. Generally speaking, as γ decreases, the pressure oscillation decreases. However, if the value of γ is excessive small, the value of the source term would tend to zero, causing the slamming feature to be erased. In order to further stabilize the pressure on the premise of ensuring the slamming feature, a new source term is constructed by performing time integration on the pressure Poison equation in this section.

Performing time integration on both sides of equation (19), equation (36) is obtained.

$${\int }_{t^m}^{t^n}{\nabla }^{2}P(t)\mathrm{dt}=-\frac{\gamma \rho }{\mathrm{\Delta }{t}^2}{\int }_{t^m}^{t^n}\frac{n^{*}\left(t\right)-n_0}{n_0}\mathrm{dt}$$

(36)

In equation (36), $n^{*}\left(t\right)$ represents the intermediate particle density at time $t$ . $t^n$ and $t^m$ represent the time of the $n^{\mathrm{th}}$ and $m^{\mathrm{th}}$ time step respectively. If $(t^n-t^m)$ is small enough, ${\int }_{t^m}^{t^n}{\nabla }^{2}P(t)\mathrm{dt}$ can be approximated replaced by ${\nabla }^2{P}^n(t^n-t^m)$ :

$${\nabla }^2{P}^n(t^n-t^m)\approx {\int }_{t^m}^{t^n}{\nabla }^{2}P(t)\mathrm{dt}$$

(37)

Substitute equation (37) into (36), a time integral source term of pressure Poison equation (equation (38)) is obtained.

$${\nabla }^2{P}^n=-\frac{1}{t^n-t^m}\frac{\gamma \rho }{\mathrm{\Delta }{t}^2}{\int }_{t^m}^{t^n}\frac{n^{*}\left(t\right)-n_0}{n_0}\mathrm{dt}$$

(38)

The discrete form of the time integral source term is given in equation (39):

$$\begin{array}{c}\multicolumn{1}{c}{{\nabla }^2{P}^n=-\frac{1}{t^n-t^m}\frac{\gamma \rho }{\mathrm{\Delta }{t}^2}{\sum }_{m+1\le i\le n}}\\ \left[\frac{n^{*}\left(t^i\right)-n_0}{n_0}(t^i-t^{i-1})\right]\hfill \end{array}$$

(39)

where, integer $i$ represents the $i^{\mathrm{th}}$ time step. $(t^n-t^m)$ should satisfy the following inequality:

$$0<(t^n-t^m)\le \mathrm{\Delta }{T}_{\max }$$

(40)

where, $\mathrm{\Delta }{T}_{\max }$ is the maximum acceptable value of $(t^n-t^m)$ . In this paper, $\mathrm{\Delta }{T}_{\max }=0.01s$ is recommended.

Revise model for free surface particles

According to the MPS gradient operator in equation (3),

$$<\nabla P{>}_i=\frac{d}{n_0}{\sum }_{j\ne i}[\frac{P_j-P_i}{{r_{\mathrm{ij}}}^2}{\stackrel{\rightharpoonup }{r}}_{\mathrm{ij}}w(r_{\mathrm{ij}})]$$

(41)

Let,

$$<\nabla P{>}_i={\sum }_{j\ne i}<\nabla P{>}_{\mathrm{ij}}$$

(42)

where,

$$<\nabla P{>}_{\mathrm{ij}}=\frac{d}{n_0}\frac{P_j-P_i}{r_{\mathrm{ij}}^2}{\stackrel{\rightharpoonup }{r}}_{\mathrm{ij}}w(r_{\mathrm{ij}})$$

(43)

In equation (42), $<\nabla P{>}_{\mathrm{ij}}$ represents the interaction force between particle $i$ and particle $j$ .

In traditional MPS methods, the pressure of free surface particles is set to the same value (usually 0 Pa), thus for free surface particles

$$\begin{array}{cc}{P}_j=P_i\hfill & i,jϵ\mathrm{free}\mathrm{surface}\hfill \end{array}$$

(44)

According to equations (42) and (44)

$$\begin{array}{cc}<\nabla P{>}_{\mathrm{ij}}\equiv 0\hfill & i,jϵ\mathrm{free}\mathrm{surface}\hfill \end{array}$$

(45)

Equation (45) indicates that there is no interacted repulsion force between free surface particles no matter how close they are. This may lead to particle overlapping and particle thrust on the free surface, and finally cause numerical instability. In order to keep a reasonable distance between free surface particles to suppress particle thrust, we derive a displacement revise model for free surface particles to modify their displacement.

According to the MPS Laplacian operator,

$$<{\nabla }^{2}P{>}_i=\frac{2d}{{\lambda }_0{n}_0}{\sum }_{j\ne i}(P_j-P_i)w(r_{\mathrm{ij}})$$

(46)

According to equation (46), we can found that $<{\nabla }^{2}P{>}_i$ is a linear combination of $(P_j-P_i)w\left(r_{\mathrm{ij}}\right)(j\ne i)$ . Let

$$<{\nabla }^{2}P{>}_i={\sum }_{j\ne i}<{\nabla }^{2}P{>}_{\mathrm{ij}}$$

(47)

where

$$<{\nabla }^{2}P{>}_{\mathrm{ij}}=\frac{2d}{{\lambda }^0{n}^0}(P_j-P_i)w(r_{\mathrm{ij}})$$

(48)

Substitute equation (19) into equation (5):

$$\begin{array}{c}\multicolumn{1}{c}{<{\nabla }^{2}P{>}_i=-\gamma \frac{\rho }{\mathrm{\Delta }{t}^2}\frac{n^{*}-n_0}{n_0}}\\ \multicolumn{1}{c}{=-\gamma \frac{\rho }{\mathrm{\Delta }{t}^2}\frac{\sum w(r_{\mathrm{ij}})-\sum w({r_{\mathrm{ij}}}^0)}{n_0}}\end{array}$$

(49)

Hence,

$$<{\nabla }^{2}P{>}_{\mathrm{ij}}=-\gamma \frac{\rho }{\mathrm{\Delta }{t}^2}\frac{w(r_{\mathrm{ij}})-w({r_{\mathrm{ij}}}^0)}{n_0}$$

(50)

Simultaneous equations (48) and (50):

$$\frac{2d}{{\lambda }^0{n}^0}(P_j-P_i)w\left(r_{\mathrm{ij}}\right)=-\gamma \frac{\rho }{\mathrm{\Delta }{t}^2}\frac{w(r_{\mathrm{ij}})-w({r_{\mathrm{ij}}}^0)}{n_0}$$

(51)

According to equation (51),

$$\begin{array}{c}\multicolumn{1}{c}{(P_j-P_i)=-\gamma \frac{{\lambda }^0{n}^0\rho }{2d\mathrm{\Delta }{t}^2}\frac{w(r_{\mathrm{ij}})-w({r_{\mathrm{ij}}}^0)}{w\left(r_{\mathrm{ij}}\right)n_0}}\\ \multicolumn{1}{c}{=-\gamma \frac{{\lambda }^0\rho }{2d\mathrm{\Delta }{t}^2}\frac{w(r_{\mathrm{ij}})-w({r_{\mathrm{ij}}}^0)}{w\left(r_{\mathrm{ij}}\right)}}\end{array}$$

(52)

In this paper, it is assumed that the interaction force arises only when the distance between two free surface particles is less than $l_0$ . Substitute $\varnothing$ for P to denote the interaction force between the free surface particles in equation (52):

$$({\varnothing }_j-{\varnothing }_i)=\{\begin{array}{cc}-\gamma \frac{{\lambda }^0\rho }{2d\mathrm{\Delta }{t}^2}\frac{w(r_{\mathrm{ij}})-w(l_0)}{w\left(r_{\mathrm{ij}}\right)},\hfill & r_{\mathrm{ij}}<l_0;\hfill \\ \multicolumn{1}{c}{0,}& \mathrm{others}.\hfill \end{array}$$

(53)

Substitute equation (53) into (48), we get:

$$<\nabla \varnothing {>}_{\mathrm{ij}}=\{\begin{array}{cc}-\gamma \frac{{\lambda }^0\rho }{2d{n}_0\mathrm{\Delta }{t}^2}\frac{w\left(r_{\mathrm{ij}}\right)-w\left(l_0\right)}{r_{\mathrm{ij}}^2},\hfill & r_{\mathrm{ij}}<l_0;\hfill \\ \multicolumn{1}{c}{0,}& \mathrm{others}.\hfill \end{array}$$

(54)

If the standard kernel function is adopted, and the searching radius is 2.1 $l_0$ , then:

$$<\nabla \varnothing {>}_{\mathrm{ij}}=\{\begin{array}{cc}-\gamma \frac{{\lambda }^0\rho }{2n_0\mathrm{\Delta }{t}^2}\frac{w\left(r_{\mathrm{ij}}\right)-1.1}{r_{\mathrm{ij}}^2}{\stackrel{\rightharpoonup }{r}}_{\mathrm{ij}},\hfill & r_{\mathrm{ij}}<l_0;\hfill \\ \multicolumn{1}{c}{0,}& \mathrm{others}.\hfill \end{array}$$

(55)

Therefore,

$$\mathrm{\Delta }{\stackrel{\rightharpoonup }{x}}_{\mathrm{ij}}^{\prime }=-\frac{\mathrm{\Delta }{t}^2}{\rho }<\nabla \varnothing {>}_{\mathrm{ij}}=\{\begin{array}{cc}\gamma \frac{{\lambda }^0}{2n_0}\frac{w(r_{\mathrm{ij}})-1.1}{r_{\mathrm{ij}}^2}{\stackrel{\rightharpoonup }{r}}_{\mathrm{ij}},\hfill & r_{\mathrm{ij}}<l_0;\hfill \\ \multicolumn{1}{c}{0,}& \mathrm{others}.\hfill \end{array}$$

(56)

where, $\mathrm{\Delta }{\stackrel{\rightharpoonup }{x}}_{\mathrm{ij}}^{\prime }$ represents the corrected displacement component caused by particle $j$ .

Hence, the corrected displacement of particle $i$ is

$$\begin{array}{cc}\mathrm{\Delta }{\stackrel{\rightharpoonup }{x}}_i^{\prime }=-\frac{\mathrm{\Delta }{t}^2}{\rho }{\sum }_{j\in \mathrm{free}\mathrm{surface}}<\nabla \varnothing {>}_{\mathrm{ij}},\hfill & i\in \mathrm{free}\mathrm{surface}\hfill \end{array}$$

(57)

Considering that the distance between ordinary particles and free surface particles may also be less than $l_0$ , particle j should include ordinary particles to maintain a reasonable distance between a free surface particle and its neighboring particles, thus,

$$\begin{array}{l}\mathrm{\Delta }\stackrel{\rightharpoonup }{x}{}_i^{\prime }=-\frac{\mathrm{\Delta }{t}^2}{\rho }{\sum }_{j\in theneighboringparticlesofi}<\nabla \varnothing {>}_{ij},\\ i\in freesurface\end{array}$$

(58)

It can be seen in equation (58) that the calculation of the modified displacement does not change the pressure of the free surface particles, nor the initial particle number density on the free surface.

Free surface particle detection algorithm

In this paper, a twice-judgment algorithm proposed by Pan et al. ³⁵ is adopted to detect free surface particle in the ADMPS method.

If a condition

$$\frac{n_i^{*}}{n_0}\ge {\beta }_1$$

(59)

is satisfied, then particle $i$ is recognized as a normal particle; if

$$\frac{n_i^{*}}{n_0}\le {\beta }_2$$

(60)

then, particle $i$ is recognized as a free surface particle; when a condition

$${\beta }_2<\frac{n_i^{*}}{n_0}<{\beta }_1$$

(61)

is satisfied, the secondary judgment is required. ${\beta }_1$ , ${\beta }_2$ are user-defined tuning parameters. In this paper, we chose ${\beta }_2=0.85$ and ${\beta }_1=0.98$ .

For the particle that requires secondary judgment, its neighboring region between 1.0 $l_0$ and 2.1 $l_0$ is divided into eight parts (as Figure 2 shows). If there are not any neighboring particles in two adjacent parts, then the particle is recognized as a free surface one.

Figure 2.

Schematic diagram of the secondary judgment.

Results and discussion

In this section, several numerical tests are conducted to show the advantage of the proposed methods. Table 1 presents a brief description of the computational models for verification analysis. As shown in Table 1, five models including MPS, MPS-FR, ADMPS, ADMPS-FR, ADMPS-FR-TIS are involved. The Free surface particle detection algorithm proposed by Pan et al. ³⁵ is adopted in all these models. A hydrostatic problem and a dam breaking problem are simulated by these models. The hydrostatic problem is for the verification of low-pressure status, while the dam breaking one is for the verification of high-pressure status.

Table 1.

Computational models for verification analysis.

Model name	Description
MPS	MPS scheme with equation (19)
MPS-FR	MPS modified by equation (58)
ADMPS	ADMPS scheme with equation (19)
ADMPS-FR	ADMPS modified by equation (58)
ADMPS-FR-TIS	ADMPS modified by equations (38) and (58)

Hydrostatic pressure problem

Simulation condition

In this section, a hydrostatic pressure problem is simulated by the proposed method. Figure 3 shows the simulation domain. It’s a $0.36m$ long rectangular pool with water depth $H=0.48m$ . The spatial resolution $l_0=0.015m$ , and the time step length $\mathrm{\Delta }t=0.0002s$ .

Figure 3.

Simulation domain of the hydrostatic calculation.

Simulation results

Figure 4 shows the particle and pressure distribution of hydrostatic simulation t = 4.0 s. It can be seen from the figure that, the particle distributions of these five models are all uniform, as well as the pressure distributions.

Figure 4.

Fluid behavior and pressure distribution of hydrostatic simulation at t = 4.0 s.

Figure 5 shows the time histories of pressure on P1 calculated by the MPS, MPS-FR, ADMPS, ADMPS-FR, and ADMPS-FR-TIS models. All these pressure curves follow a similar variation law. After a short oscillation in the initial stage, all curves tend to be stable gradually, and finally in good agreement with the analytical solution (the differences are all less than 100 Pa). The steady-state pressures of the ADMPS models are almost the same as those of the MPS models, because the predicted component of the intermediate velocity can be ignored in the hydrostatic problem. Compared with these models, it can be found that the stability of the models with FR technology is better than that of the models without FR technology. In addition, the oscillation frequency of the ADMPS models is smaller than that of the MPS models. We also found that, the oscillation amplitude of ADMPS-FR-TIS is much smaller than that of the other four models. This means that the time integral source term improves the stability of hydrostatic pressure calculation and makes the pressure curve more smooth.

Figure 5.

Time histories of pressure on P1.

Dam breaking problem

The dam breaking experiment conducted by Lobovský et al. ⁴² is simulated by the five models. The experiment was conducted in a glass tank with 0.600 m height and 1.610 m length (as shown in Figure 6). The initial height of the water column H = 0.600 m and the initial width B = 0.600 m as well. A monitoring point P2 was set to observe the pressure variation. As shown in Figure 6, P2 is located at the left corner of pool and is 3 mm from the bottom. Obviously, P2 is in the area where strong slamming would occur. The simulations are performed in the spatial resolution $l_0=0.0100\mathrm{m}$ . The time step length $\mathrm{\Delta }t=0.0002\mathrm{s}$ .

Figure 6.

Schematic diagram of pool size.

Figure 7 shows the particle distributions of dam breaking simulation at t = 0.90 s. From the simulation results of the MPS and the ADMPS models, we can found that, the phenomenon of particle aggregation occurs near free surface, and some free surface particles are even almost overlapped. By contrast, according to the results of the MPS-FR and the ADMPS-FR models, the phenomenon of particle aggregation does not appear on free surface, and the distributions of particles are more uniform than that of the MPS or ADMPS model. The results indicate that the FR technology is capable to provide a reasonable particle distance on free surface.

Figure 7.

Particle distributions of dam breaking simulations at t = 0.9 s.

In Figure 8, fluid flows at t = 0.316, 0.413, and 0.463 s simulated by MPS, MPS-FR, ADMPS, ADMPS-FR, and ADMPS-FR-TIS are compared to the experiment conducted by Lobovský et al. ⁴² In general, the shapes of fluid flow simulated by these models are in good agreement with the experimental results, especially at 0.316 s (before the slamming on vertical wall) and 0.413 s (the slamming just happened on the vertical wall). However, the raising water columns of these models at 0.463 s are slightly different. Water splash of the models with FR technology is more remarkable compared to that of the models without FR technology. Compared with the other three models, the water column shapes of ADMPS-FR and ADMPS-FR-TIS are more consistent with the experiment result. It can also be found in Figure 8 that the ADMPS-FR-TIS model provides a more uniform pressure distribution than other models.

Figure 8.

Snap shot of dam breaking fluid flows of experiment and simulations at (a) t = 0.316 s, (b) t = 0.413 s, and (c) t = 0.463 s.

Figure 9 shows the time histories of pressure calculated by MPS and MPS-FR at point P2. As shown in the figure, the median values of these two curves are significantly smaller than the experimental values most of the time after slamming (at dimensionless time = 1.5). The two curves have significant oscillations between dimensionless time of 1.5 and 2.6, and the amplitude of oscillation gradually decreases with time. Comparing these two curves, it can also be found that the oscillation amplitude of the MPS-FR curve is slightly smaller than that of the MPS curve generally, and the maximum value of the MPS-FR curve (about 5.2) is smaller than that of MPS curve (greater than 6). Figure 10 shows the time histories of pressure calculated by ADMPS, ADMPS-FR and ADMPS-FR-TIS at point P2. In this figure, the median values of the ADMPS curves are more consistent with the experimental values compared with that of the MPS curves in Figure 9. The oscillation law of the ADMPS curves is similar to that of MPS curves, that is the oscillation amplitude gradually decreases with time during 1.5–2.6. Comparing these curves, it can be found that the oscillation amplitude of the ADMPS-FR-TIS curve is significantly smaller than that of the other two curves, and the oscillation amplitude of the ADMPS-FR curve is slightly smaller than that of the ADMPS curve. In addition, the maximum value of the ADMPS-FR-TIS curve (about 3.6) is much smaller than other curves as well.

Figure 9.

Pressure curves on P2 calculated by MPS and MPS-FR.

Figure 10.

Pressure curves on P2 calculated by ADMPS, ADMPS-FR, and ADMPS-FR-TIS.

Conclusion

In the paper, an accurate and stable improved MPS method based on Crank-Nicolson time scheme for incompressible flow with free surfaces is proposed. The present methodology includes a Lagrangian alternating direction moving particle semi-implicit (ADMPS) algorithm to improve intermediate velocity, a time integral source term (TIS) to suppress unphysical pressure oscillation and a free-surface revise (FR) model to avoid particle thrust near free surface. Hydrostatic and dam breaking examples are applied to verify these techniques. Results of five numerical models including MPS, MPS-FR, ADMPS, ADMPS-FR, ADMPS-FR-TIS are compared in the paper. By discussing and comparing the simulation’s results of these methods, we draw the following conclusions:

1. The proposed FR model is proved to be effective for preventing particle thrust on free surface;

2. The slamming pressure curves of the ADMPS-FR models are more accurate than those of the MPS models when applying to simulate dam breaking process.

3. The proposed TIS model is demonstrated to be effective for suppressing unphysical pressure oscillation;

4. Remarkable enhancements of stability and accuracy is demonstrated when using the ADMPS-FR-TIS model.