Study on divergence approximation formula for pressure calculation in particle method

Abstract

The moving particle semi-implicit method is a meshless particle method for incompressible fluid and has proven useful in a wide variety of engineering applications of free-surface flows. Despite its wide applicability, the moving particle semi-implicit method has the defects of spurious unphysical pressure oscillation. Three various divergence approximation formulas, including basic divergence approximation formula, difference divergence approximation formula, and symmetric divergence approximation formula are proposed in this paper. The proposed three divergence approximation formulas are then applied for discretization of source term in pressure Poisson equation. Two numerical tests, including hydrostatic pressure problem and dam-breaking problem, are carried out to assess the performance of different formulas in enhancing and stabilizing the pressure calculation. The results demonstrate that the pressure calculated by basic divergence approximation formula and difference divergence approximation formula fluctuates severely. However, application of symmetric divergence approximation formula can result in a more accurate and stabilized pressure.

Keywords

Moving particle semi-implicit method pressure Poisson equation divergence approximation formula

Introduction

The moving particle semi-implicit (MPS) method is a Lagrangian gridless method originally proposed by Koshizuka¹ for incompressible free-surface flows. The MPS method is applied in a wide variety of engineering fields such as nuclear engineering,^2–7 ocean engineering,^8–10 bioengineering,¹¹ etc. Despite the wide application in various engineering fields, the MPS method suffers the drawbacks that may significantly affect the stability and accuracy of pressure calculation. Many attempts have been made to resolve this problem. Sueyoshi^12,13 studied the nonlinear phenomena in naval architecture and ocean engineering and proposed the time average method to control the pressure oscillations. Hibi and Yabushita¹⁴ investigated the unphysical pressure fluctuations and proposed the modified MPS method with postprocessing for incompressible fluid flow. Khayyer and Gotoh¹⁵ proposed a modified MPS method and improved the performance of the prediction of wave impact pressure. Khayyer and Gotoh¹⁶ further improved the error-compensating parts in the multiterm source, and a corrective matrix for pressure gradient is proposed. The proposed improvements are shown to stabilize and enhance the performance of MPS method. Kondo and Koshizuka¹⁷ proposed a new formulation for the source term of pressure Poisson equation which consisted of three parts, one main part and two error-compensating parts. Lee et al.¹⁸ illustrated that there are several defects including nonoptimal source term, gradient and collision models, and detection of free-surface particle, which caused the unphysical pressure fluctuations. These defects are remedied by step-by-step improvements. Taimai and Koshizuka¹⁹ introduced a new methodology including arbitrary high-order accurate mesh-free spatial discretization schemes, and this new methodology showed drastic improvement of accuracy and stability. Daneshvar et al.²⁰ proposed two modified pressure gradient models for discrete pressure gradients and improved the pressure calculation.

In this paper, three various divergence approximation formulas, including basic divergence approximation formula (BDAF), difference divergence approximation formula (DDAF), and symmetric divergence approximation formula (SDAF) are proposed and used to discretize the source term in pressure Poisson equation. We applied these formulas to two problems: hydrostatic pressure problem and dam-breaking problem. Performance in pressure calculation by different divergence approximation formulas is discussed.

MPS method

Governing equations

MPS method is an approximation method for particle systems; the mass and momentum conservation equations of incompressible flows in MPS method are as follows

\frac{D ρ}{D t} = 0

(1)

ρ \frac{D u}{D t} = - \nabla P + μ \nabla^{2} u + G

(2)

where

D / D t

means the Lagrangian derivative involving convection terms,

u

represents the velocity of fluid,

P

represents the pressure of fluid,

ρ

represents the density of fluids,

μ

represents the viscosity of fluid, and

G

represents the gravity acceleration.

Particle interaction models

The governing equations written with partial differential operators are transformed to the equation of particle interactions. The particle interactions in the MPS are based on the weight function. In the traditional MPS method, the following weight function is employed

w (r) = {\begin{array}{l} \frac{r_{e}}{r} - 1, & 0 < r < r_{e} \\ 0 & r \geq r_{e} \end{array}

(3)

The distance between two particles is $r$ and $r_{e}$ represents the effective range of particle interactions. Since the area covered with this weight function is bounded, a particle interacts with a finite number of neighboring particles. The nearer the distance between two particles, the larger the weight of interactions. If the distance between two particles is quite long, their interactions can be neglected. Summation of the weight functions is called particle number density, $n_{i}$ , which is used as a normalization factor for averaging

n_{i} = \sum_{j \neq i} w (| r_{j} - r_{i} |)

(4)

The particle number density represents its fluid density. It should be constant for incompressible flows: $n_{i} = n^{0}$ , where $n^{0}$ is dependent on the initial arrangement of particles.

Gradient model

The gradient operator in the governing equations is transformed to equivalent particle interactions. If $Φ$ is an arbitrary scalar, particle interaction models for differential operators are expressed as

{〈 \nabla Φ 〉}_{i} = \frac{d}{n^{0}} \sum_{j \neq i} [\frac{Φ_{j} - {\hat{Φ}}_{i}}{{| r_{j} - r_{i} |}^{2}} (r_{j} - r_{i}) w (| r_{j} - r_{i} |)]

(5)

where

d

is the number of space dimensions,

n^{0}

is the particle number density fixed for incompressibility of the initial condition of particle arrangement, and

{\hat{Φ}}_{i}

is the minimum value in neighbor particles. The gradient model is calculated as the average of gradient vectors of the neighbors of the

i

th particle.

Laplacian model

In the traditional MPS method, the Laplacian operator model is derived from the physical concept of diffusion. The Laplacian operator in the governing equations is transformed to equivalent particle interactions. If $Φ$ is an arbitrary scalar, particle interaction models for differential operators are expressed as

{〈 \nabla^{2} Φ 〉}_{i} = \frac{2 d}{λ n_{0}} \sum_{j \neq i} (φ_{j} - φ_{i}) w (| r_{j} - r_{i} |)

(6)

where

d

is the space dimension number and

λ

is defined as

λ = \frac{\sum_{j \neq i} {| r_{j} - r_{i} |}^{2} w (| r_{j} - r_{i} |)}{\sum_{j \neq i} w (| r_{j} - r_{i} |)}

(7)

Parameter $λ$ is introduced to keep the same variance increase as that of the analytical solution. The Laplacian operator model is derived from the physical concept of diffusion. Substituting the above particle interaction models into governing equations, we can obtain the particle dynamics to simulate fluid flows.

However, this formulation of discretization scheme in the traditional MPS method has some disadvantages, including no foundation of exact mathematical derivation, only one order of precision, and unphysical pressure fluctuation effects. In this paper, a brand new second derivative for pressure Poisson equation is proposed which is based on the idea of smoothed particle hydrodynamics (SPH). SPH is a particle method which is similar to MPS method. SPH was proposed by Lucy²¹ and Gingold²² independently for astrophysical problems. We first briefly introduce the basis of interpolation integral in the SPH method. In the SPH method, a quantity, which is a function of the spatial coordinates, is based on the integral interpolant

A (r) = \int A (r^{'}) W (r - r^{'}) d r^{'}

(8)

where the function

W (r - r^{'})

is the kernel function and

d r^{'}

is a differential volume element. The integral can then be approximated by a summation over the mass elements. This gives the summation interpolant

A (r) = \sum_{b} m_{b} \frac{A_{b}}{ρ_{b}} W (r - r_{b})

(9)

where the summation is over all the neighbor particles. Brookshaw²³ solved heat diffusion in three-dimensional particle simulations by a new second-derivative formulation. The following new second derivative for Laplacian model is analogous to the approximation used by Cummins and Rudman²⁴ for thermal diffusion and the viscous diffusion term. The derivation process is as follows.

For an arbitrary scalar $f (r^{'})$ , by Taylor series expansion of $f (r^{'})$ about $r^{'} = r$

\int \frac{[f (r^{'}) - f (r)]}{(r - r^{'})} \frac{\partial W (r - r^{'})}{r} d r^{'} = \frac{1}{2} \frac{d^{2} f}{d r^{2}} + o ({(Δ r)}^{2})

(10)

discretization form of integral approximation follows

{(\frac{d^{2} f}{d r^{2}})}_{i} = 2 \sum_{j \neq i} \frac{m_{j}}{ρ_{j}} \frac{(f_{i} - f_{j})}{(r_{i} - r_{j})} \frac{\partial W_{i j}}{\partial r_{i}} + o ({(Δ r)}^{2})

(11)

where

W_{i j} = W (| {\overset{⇀}{r}}_{i} - {\overset{⇀}{r}}_{j} |) .

The approximation of the second-derivative integral interpolant can be expressed as follows

\begin{array}{l} {(\frac{d^{2} f}{d r^{2}})}_{i} \approx 2 \sum_{j \neq i} \frac{m_{j}}{ρ_{j}} \frac{(f_{i} - f_{j})}{(r_{i} - r_{j})} \frac{\partial W_{i j}}{\partial r_{i}} \\ = 2 \sum_{j \neq i} \frac{m_{j}}{ρ_{j}} \frac{(f_{i} - f_{j}) (r_{i} - r_{j})}{{(r_{i} - r_{j})}^{2}} \frac{\partial W_{i j}}{\partial r_{i}} \end{array}

(12)

To keep the stability of simulation, a formula can be written as follows

{(\frac{d^{2} f}{d r^{2}})}_{i} \approx 2 \sum_{j \neq i} \frac{m_{j}}{ρ_{j}} \frac{(f_{i} - f_{j}) (r_{i} - r_{j})}{{(r_{i} - r_{j})}^{2} + {(0.1 l_{0})}^{2}} \frac{\partial W_{i j}}{\partial r_{i}}

(13)

${(0.1 l_{0})}^{2}$ is a small parameter included to ensure that the denominator remains nonzero. Applying the above approximation discretization form to pressure Poisson equation

{〈 \nabla^{2} P^{n + 1} 〉}_{i} = 2 \sum_{j \neq i} \frac{m_{j}}{ρ_{j}} \frac{{\overset{⇀}{r}}_{i j} \cdot \nabla_{i} W_{i j}}{{| {\overset{⇀}{r}}_{i j} |}^{2} + {(0.1 l_{0})}^{2}} (P_{i} - P_{j})

(14)

where

{\overset{⇀}{r}}_{i j} = {\overset{⇀}{r}}_{i} - {\overset{⇀}{r}}_{j}

\nabla_{i} W_{i j} = \partial W_{i j} / \partial r

. Equation (31) is the new second derivative for Laplacian model. Our new second derivative is an improvement of Cummins’s²⁴ and Brookshaw’s proposal. It has a simpler formulation than the approximate projection operator used by Cummins. Consequently, it saves the computational time and enhances the convergence of simultaneous linear equations. Meanwhile, this formula is insensitive to particle disorder and stabilizes the pressure simulation which the original formulation proposed by Cummins. At the same time, the kernel function used in this paper also differs from the kernel function utilized by Cummins and Brookshaw.

Kernel function

In the traditional MPS method, the kernel function represents the weight of physical quantities. At the same time, it does not require for the derivatives of kernel function. The following weight function is employed in the traditional MPS method

w (r) = {\begin{matrix} \frac{r_{e}}{r} - 1, & 0 < r < r_{e} \\ 0 & r \geq r_{e} \end{matrix}

(15)

The distance between two particles is $r$ and $r_{e}$ represents the effective range of particle interactions. To employ the new second derivative for PPE in equation (14), a new smoothing kernel function which matches the SPH frame should be put forward. The smoothing kernel function should meet the requirement of normalization condition in the compact support domain; the following function is the criterion of normalization condition

\int_{Ω} W (r_{i j}, h) d r = 1

(16)

where the

W

means kernel function,

r_{i j}

represents the distance between two particles, and

h

is the smooth length. More generally, in order to ensure the second-order accuracy of the approximations, the kernel function must verify the following conditions which is mentioned in Koshizuka¹

{\begin{array}{l} \int_{Ω} W (x - x^{'}, h) d x^{'} = 1 \\ \int_{Ω} (x - x^{'}) W (x - x^{'}, h) d x^{'} = 0 \\ \int_{Ω} {(x - x^{'})}^{2} W (x - x^{'}, h) d x^{'} = 0 \end{array}

(17)

The quadratic smoothing function has advantages in compressive stability and computational efficiency. The quadratic smoothing function and the derivative is as follows

\begin{matrix} W (s, h) = α_{d} (\frac{3}{16} s^{2} - \frac{3}{4} s + \frac{3}{4}), 0 \leq s \leq 2 \end{matrix}

(18)

\begin{array}{l} \nabla_{i} W_{i j} = \frac{{\overset{⇀}{r}}_{i j}}{r_{i j}} \frac{\partial W_{i j}}{\partial r_{i j}} \\ = α_{d} \frac{1}{h} \frac{{\overset{⇀}{r}}_{i j}}{r_{i j}} (\frac{3}{8} s - \frac{3}{4}), 0 \leq s \leq 2 \end{array}

(19)

where

s = r_{i j} / h

r_{i j}

represents the distance between two particles and

h

is the smooth length. In the MPS method, the smooth length

h

equals to

1.05 l_{0}

in particle number density model and gradient model. In Laplacian model, the smooth length

h

equals to

2.0 l_{0}

. The parameter

α_{d}

1 / h

2 / π h^{2}

, and

5 / 4 π h^{3}

in one-dimensional, two-dimensional, and three-dimensional space, respectively. The quadratic smoothing function is employed in this paper.

Pressure Poisson equation

The divergence free velocity field is employed as the source of pressure Poisson equation which is found on the Helmholtz–Hodge decomposition of an intermediate velocity field. According to the Helmholtz–Hodge decomposition, any vector field, $V$ , can be decomposed into a divergence-free component, $V^{d}$ , and the gradient of some scalar, $\nabla \emptyset$ (a curl-free component). Thus

V = V^{d} + \nabla φ

(20)

For the N–S equation, it can be written as follows

\frac{u^{n + 1} - u^{n}}{Δ t} = - \frac{1}{ρ} \nabla P + \frac{μ}{ρ} \nabla^{2} u + \frac{1}{ρ} G

(21)

And the left side of the equation (17) can be written as

\begin{matrix} \frac{u^{n + 1} - u^{n}}{Δ t} = \frac{u^{n + 1} - u^{*} + u^{*} - u^{n}}{Δ t} \\ = \frac{u^{n + 1} - u^{*}}{Δ t} + \frac{u^{*} - u^{n}}{Δ t} \end{matrix}

(22)

In the tradition MPS method, the equation can be decomposed to equations as follows

\frac{u^{*} - u^{n}}{Δ t} = \frac{μ}{ρ} \nabla^{2} u + \frac{1}{ρ} G

(23)

\frac{u^{n + 1} - u^{*}}{Δ t} = - \frac{1}{ρ} \nabla P

(24)

Equation (24) can be written as follows

u^{*} = u^{n + 1} + \frac{Δ t}{ρ} \nabla P

(25)

Take the divergence on both sides of equation (25)

\nabla \cdot u^{*} = \nabla \cdot u^{n + 1} + \frac{1}{Δ t} \cdot \nabla \cdot (\frac{1}{ρ} \nabla P)

(26)

The divergence of velocity field at $n + 1$ step equals to zero for the incompressible fluid, which satisfies the Helmholtz–Hodge decomposition. As a result, equation (26) can be expressed as

\nabla^{2} P = \frac{ρ}{Δ t} \nabla \cdot u^{*}

(27)

On the right side of equation (27), the source term of the pressure Poisson equation is the divergence-free velocity.

Choosing the divergence approximation formula of source term

In the SPH frame, the derivative of a function can be obtained by using the formula of equation (9). The following formula is the divergence of the velocity field obtained in this way

{〈 \nabla \cdot u 〉}_{i} = \sum_{j \neq i} \frac{m_{j}}{ρ_{j}} u_{j} \cdot \nabla W (| {\overset{⇀}{r}}_{j} - {\overset{⇀}{r}}_{i} |)

(28)

In the present paper, the formula in equation (28) is called the BDAF. Equation (28) is based on the basic approximation idea of SPH method and has a simple discretization scheme.

Apart from the BDAF, the divergence value for the velocity can be calculated by other ways. The second formula is derived by utilizing the derivative of a product. The process of derivation is as follows, first we have

\nabla \cdot (ρ u) = ρ \nabla \cdot u + u \nabla \cdot ρ

(29)

which can be rewritten as

ρ \nabla \cdot u = \nabla \cdot (ρ u) - u \nabla \cdot ρ

(30)

This leads to the following expression

ρ_{i} \nabla \cdot u_{i} = \sum_{j} m_{j} (u_{j} - u_{i}) \cdot \nabla W (| {\overset{⇀}{r}}_{j} - {\overset{⇀}{r}}_{i} |)

(31)

Thus, the divergence of the velocity can be approximated as

\nabla \cdot u_{i} = (\frac{1}{ρ_{i}}) \sum_{j} m_{j} (u_{j} - u_{i}) \cdot \nabla W (| {\overset{⇀}{r}}_{j} - {\overset{⇀}{r}}_{i} |)

(32)

We called equation (32) the DDAF in this paper.

In order to satisfy Newton’s second law with two particles, some researchers adopt the derivative of a quotient function to discrete the divergence. The derivative of a quotient of functions is expressed by

\nabla \cdot (\frac{u}{ρ}) = \frac{ρ \nabla \cdot u - u \nabla \cdot ρ}{ρ^{2}}

(33)

The above equation can be written as follows

\frac{\nabla \cdot u}{ρ} = \nabla \cdot (\frac{u}{ρ}) + \frac{u \nabla \cdot ρ}{ρ^{2}}

(34)

Thus, the following formula can be obtained

\frac{\nabla \cdot u_{i}}{ρ_{i}} = \sum_{j} m_{j} (\frac{u_{i}}{ρ_{i}^{2}} + \frac{u_{j}}{ρ_{j}^{2}}) \cdot \nabla W (| {\overset{⇀}{r}}_{j} - {\overset{⇀}{r}}_{i} |)

(35)

Thus, the divergence of the velocity can be approximated as

\nabla \cdot u_{i} = ρ_{i} \sum_{j} m_{j} (\frac{u_{i}}{ρ_{i}^{2}} + \frac{u_{j}}{ρ_{j}^{2}}) \cdot \nabla W (| {\overset{⇀}{r}}_{j} - {\overset{⇀}{r}}_{i} |)

(36)

We called equation (36) the SDAF in this paper.

Hydrostatic pressure problem

The hydrostatic pressure problem is chosen to verify pressure calculation as the first problem. In order to understand the differences in the divergence approximation of BDAF, DDAF, and SDAF, the hydrostatic pressure problem is simulated to compare the different divergence formula for simulating stable fluid flows. Figure 1 shows the simulation domain. The height of the water is 0.48 m and the width of the water is 0.36 m. The fluid mass density is $1000 kg / m^{3}$ and the kinematic viscosity coefficient is $1.0 \times 10^{- 6} m^{2} / s$ . The acceleration of gravity is $9.8 m / s^{2}$ . The calculation time step for the verified problem is adopted based on Courant–Friedrichs–Lewy condition. $Δ t = C_{C F L} \cdot l_{0} / u_{\max}$ , where $u_{\max}$ is the maximum velocity of fluid particle. The initial particle distance is $0.004 m$ and the total particle number is 12,392. The time variation of hydrostatic pressure at point P is recorded. The analytical solution is calculated as follows

P = ρ g H

(37)

where

ρ

is the fluid mass density,

g

is the acceleration of gravity, and

H

is the height of the fluid.

Figure 1.

Schematic view of hydrostatic pressure problem.

Accuracy of divergence formula

In this section, three different divergence formulas are explicitly evaluated by the hydrostatic pressure problem. Figure 2 shows the divergence approximation results of point P by BDAF, DDAF, and SDAF at different time. In our simulations, however, the three different forms of divergence approximation behave differently, as presented in Figure 2. SDAF acquired the least error result and yielded a more accurate solution. While the other two formulas had much greater errors and are not as accurate as the SDAF performs. We also observe that BDAF had a larger error than DDAF.

Figure 2.

Divergence value of hydrostatic pressure problem by different formulas. BDAF: basic divergence approximation formula; DDAF: difference divergence approximation formula; SDAF: symmetric divergence approximation formula.

Comparison of pressure field

The time histories of pressure at measuring point P are shown in Figure 3. From this figure, the BDAF has resulted in largest amplitude of pressure similar to the amplitude of divergence seen in Figure 3. The amplitude of pressure of DDAF is second largest among the three divergence formulas. Because solving the pressure Poisson equation equals to solving the large sparse matrix

A x = b

(38)

Figure 3.

Comparison of pressure field of different divergence formulas. BDAF: basic divergence approximation formula; DDAF: difference divergence approximation formula; SDAF: symmetric divergence approximation formula.

The system of equations actually solved is as follows

(A + δ A) x = (b + δ b)

(39)

If the || $δ A$ || and || $δ b$ || are small, however || $δ x$ || is large during the calculation, we call the equation set $A x = b$ ill-conditioned equation set. On the contrary, if the || $δ A$ || and || $δ b$ || are large, however || $δ x$ || is small during the simulation, we call the equation set $A x = b$ well-conditioned equation set. The ill-conditioned equation set will have numerical oscillation for any algorithm. As for the pressure Poisson equation, the pressure oscillation will occur. Because pressure is a direct function of time variation of source term and therefore the key issue in obtaining an accurate and stable pressure field with minimized unphysical oscillations is to achieve a smooth result of the source term. In equation (27), it is equivalent to obtain a smooth result of divergence value of velocity. Application of the SDAF has decreased the amplitude of pressure, resulting in smoothest pressure variations. As a result, SDAF appears to be the most effective formula in reducing the amplitude of pressure and the pressure result is closest to the theoretical value.

Figure 4 presents the pressure field of hydrostatic pressure problem by different divergence approximation formulas. Not surprisingly, the pressure field calculated by BDAF is characterized by large particle movement on the fluid boundary, which is deduced by the larger particle distance on the free surface by BDAF. As a result, the height of the free surface by this formula is higher than other approximation formulas. For the unstable phenomenon on the free surface, the pressure is unphysical and unstable on the free surface. Accordingly, the pressure calculation at point P is affected from the pressure oscillation on the free surface and unphysical pressure fluctuation occurs at point P. The particle vibration on the free surface by DDAF is also obvious and the height of the free surface is shorter than the result of BDAF. It means the pressure calculation on the free surface by DDAF is smoother than BDAF. Accordingly, the unphysical pressure fluctuation at point P by DDAF is smoother than BDAF. Compared to BDAF and DDAF, the particle distance is almost the same between any particles (including the particle on the free surface) by SDAF, which means the pressure calculation on the free surface is accurate and stable; as a result, it leads to smoother pressure field and less pressure fluctuation at point P.

Figure 4.

Effect of spatial resolution

In general, the stability of the pressure calculation in MPS method is related to the sensitivity to oscillation of divergence approximation formula of pressure Poisson equation. Another critical issue which may affect the stability of the pressure calculation of MPS method is the randomness of particle distribution. The spatial resolution of simulation will affect the randomness of particle distribution during the simulation. Figure 5 portrays the effect of spatial resolution (the same problem with different particle number) on the performance of different divergence approximation formulas. From the presented figure, refinement of spatial resolution from 1000 to 8000 particle numbers has considerably decreased the pressure error of hydrodynamic pressure problem. Further refinement of spatial resolution from 8000 to 16,000 particle numbers only slightly decreases the pressure error. This would signify the existence of an optimum particle size considering the desired level of computational efficiency and pressure oscillation.

Figure 5.

Average pressure error of different particle number by different divergence formulas. BDAF: basic divergence approximation formula; DDAF: difference divergence approximation formula; SDAF: symmetric divergence approximation formula.

Figure 6 shows the improvement of pressure error corresponding to resolution refinement by various divergence formulas, qualitatively. From this figure, refinement of spatial resolution has resulted in enhanced pressure fields characterized by more uniform particle distribution for all the three approximation formulas, especially the particles on the free surface. Not surprisingly, the particle distribution by the SDAG is much uniform than others with the same spatial resolution, which is consistent with the previous conclusion. For all the pressure fields by BDAF with different particle number, the particles on the surface are in severe vibration. However, the pressure by different spatial resolution is different. Because when the particle number equals to 1132, the distance between particles is bigger than the condition with larger particle numbers, such as particle number equals to 5872. When the particle moves the same proportion of initial particle distance, the divergence variation caused by low spatial resolution is larger than the high spatial resolution. As a result, the pressure error by the low resolution is larger than the high spatial resolution. When the particle number is bigger than 8000, the refinement of spatial resolution could not decrease the pressure error. It can be seen that the particle vibration on the surface still exists by BDAF and DDAF when the particle number equals to 12,392 and it could not be eliminated by the refinement of spatial resolution. The particle vibration on the surface by SDAF is also not eliminated when the particle number is larger than 8000, although it is not obvious in the last picture in Figure 6.

Figure 6.

Pressure field by refinement spatial resolution by different divergence formulas (particle numbers are 1132, 5872, and 12,392, respectively). BDAF: basic divergence approximation formula; DDAF: difference divergence approximation formula; SDAF: symmetric divergence approximation formula.

Verification by dam-breaking problem

Dam-breaking problem

Simulation of dam breaks has been studied by many researchers^20,25,26 in most particle method as a complicated bench problem. In order to understand the differences in the divergence approximation of BDAF, DDAF, and SDAF in violent flows, the dam-breaking problem is simulated to compare the different divergence formula for simulating violent fluid flows. Shibata et al.²⁶ and Hu and Kashiwagi²⁷ conducted an experimental dam-breaking problem to study the free-surface flow. In order to compare with Hu’s experimental results, the same dam-breaking model is employed in this research. As is shown in Figure 7, the height of the vertical wall is 0.22 m, the width of the tank is 1.18 m, the height and width of the water column are 0.12 and 0.68 m, respectively. In the initial condition, the water column is stable in the left side of the tank. The density of the water is $1000 kg / m^{3}$ and the kinematic viscosity coefficient is $1.0 \times 10^{- 6} m^{2} / s$ . The acceleration of gravity is $9.8 m / s^{2}$ .The pressure date of the point A is recorded every 0.01 s.

Figure 7.

Schematic view of the dam-breaking experiment.

Result and discussion

Figure 8 shows time histories of the calculated pressure at point A by different approximation formulas comparing with the experimental result by Hu. It could be seen that all the simulated results are roughly consistent with the experimental result of Hu. For the BDAF, the calculated pressure at point A fluctuates severely at some instants, including the instants of $t = 0.35 s$ , $t = 0.45 s$ , $t = 0.71 s$ , $t = 0.95 s$ , and $t = 1.05 s$ . For the pressure fluctuation at $t = 0.35 s$ , the fluid first hits the vertical wall and causes the strenuous fluid flow. While the second pressure oscillation at $t = 0.71$ s is due to the fall of water splash along the wall that forms the air entrapment. This process also increases the magnitude of unphysical oscillations at point A. For pressure fluctuation at other instants, the reason may be due to the particles’ randomness movement. For the DDAF, the condition is similar to BDAF. However, for the SDAF, it can be seen that the pressure fluctuation is alleviated and the pressure value is more close to the experimental result by Hu. It means the SDAF performs well in pressure calculation in dam-breaking problem and results in an enhanced time history of pressure.

Figure 8.

Comparison of time history of pressure of dam breaking at point A by different divergence formulas. BDAF: basic divergence approximation formula; DDAF: difference divergence approximation formula; SDAF: symmetric divergence approximation formula.

Figure 9 shows time histories of source term value by dam-breaking problem at point A by different approximation formulas. Similar to the time history of pressure in Figure 8, the fluctuations in source term value by different approximation formulas are different. The fluctuations in source term value by BDAF and DDAF are large; however, the fluctuations in source term value by SDAF are smallest. During the process of simulation, the fluid at point A has undergone through severe and unphysical expansions and compressions, which causes the variation of source term value. The SDAF tends to have a more stabilizing effect than the other two formulas, at least for this specific test. Among the applied three formulas, the most stable result belongs to SDAF. As previously stated, pressure is a direct function of time variation of source term value. Therefore, the less fluctuation in source term value, the less oscillation in pressure calculation. It can be seen in Figures 7 and 8, when the source term has a peak value at $t = 0.45 s$ , the peak value occurs in the pressure at the same time. It is consistent with the results in hydrostatic problem.

Figure 9.

Comparison of time history of source term value of dam breaking at point A by different divergence formulas. BDAF: basic divergence approximation formula; DDAF: difference divergence approximation formula; SDAF: symmetric divergence approximation formula.

Figure 10 shows the distributions of fluid particles and pressure fields at different instants for BDAF, DDAF, and SDAF. The snapshot by three different formulas is characterized by significant discrepancies in pressure field. For the BDAF, severe pressure oscillation in spatial distribution occurs at different times. It is due to the fact that the DDAF has the shortcoming of unphysical pressure fluctuations. A similar spatial distribution of pressure is proposed by DDAF. On the other hand, the smoother pressure spatial distribution is provided by SDAF. In particular, the SDAF has provided a more reliable spatial distribution of pressure at $t = 0.11 s$ . At these instants, however, the BDAF and DDAF have apparently resulted in an unphysical pressure distribution. Utilization of SDAF is expected to alleviate such unphysical pressure oscillation and results in an enhanced and accurate pressure field.

Figure 10.

Comparison of pressure field of different MPS methods. BDAF: basic divergence approximation formula; DDAF: difference divergence approximation formula; SDAF: symmetric divergence approximation formula.

Table 1 presents the comparison of CPU time for 10 steps by different divergence formulas of various particle numbers. The workstation utilized in this study is equipped with the CPU: Intel Xeon Processor E5-2680 v2*2 and the main memory: 128 GB. A single core is used without parallel computing. From Table 1, it is found that the CPU time increases with the refinement of spatial resolution. It is the result of that the calculated particle numbers increase. The CPU time of BDAF is smaller than the other two formulas for its simple discrete scheme. The CPU time of SDAF is largest because the mathematical formula of SDAF is little complicated than the other two formulas. However, the discrepancy between three approximation formulas is not significant.

Table 1.

Comparison of CPU time by different formulas.

Particle numbers	BDAF (s)	DDAF (s)	SDAF (s)
3772	1	1.12	1.18
5672	2.26	2.48	2.74
8354	4.94	5.39	6.17
12,666	11.11	12.32	13.25
15,896	17.76	19.54	21.65

BDAF: basic divergence approximation formula; CPU: central processing unit; DDAF: difference divergence approximation formula; SDAF: symmetric divergence approximation formula.

Conclusion

The MPS method is a mesh-free particle method which is popular for its flexibility. It has the shortcoming of pressure fluctuation; however, it is resolved by various improvements. In this study, three different divergence approximation formulas for source term in pressure Poisson equation are investigated. The hydrostatic pressure problem is simulated as the first case; the results show that the SDAF has a better performance than the other two formulas in suppressing pressure oscillation and obtaining a smoother pressure field. It is due to the accuracy calculation of movements on free surface particles. The dam-breaking problem is selected as the second benchmark; not surprisingly, the SDAF stabilizes the pressure calculation in dam-breaking problem and an enhanced and accurate pressure field is achieved.

Footnotes

Declaration of conflicting interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: The authors would like to acknowledge the supports from the National Natural Science Foundation of China (Grant No. 51779126,11472155) and the National High Technology Research and Development Program of China (863) (2014AA052701).

References

Koshizuka

Oka

Moving-particle semi-implicit method for fragmentation of incompressible fluid.

Nucl Sci Eng 1996; 123: 421–434.

Koshizuka

Particle method for incompressible viscous flow with fluid fragmentation. Comput Fluid Dyn 1995; 4: 29–46.

Yoon

HYKS.

A mesh-free numerical method for direct simulation of gas–liquid phase interface.

Nucl Sci Eng 1999; 133: 1–9.

Koshizuka

Oka

Application of moving particle semi-implicit method to nuclear reactor safety. Comput Fluid Dyn 2001; 9: 366–375.

Xie

Koshizuka

Oka

Simulation of drop deposition process in annular mist flow using three-dimensional particle method.

Nucl Eng Des 2005; 235: 1687–1697.

Tian

Ishiwatari

Ikejiri

et al . Numerical computation of thermally controlled steam bubble condensation using Moving Particle Semi-implicit (MPS) method. Ann Nucl Energy 2010; 37: 5–15.

Yoon

HYKS.

Direct calculation of bubble growth, departure and rise in nucleate pool boiling.

Int J Multiphase Flow 2011; 27: 277–298.

Shibata

Koshizuka

Numerical analysis of shipping water impact on a deck using a particle method.

Ocean Eng 2007; 34: 585–593.

Sueyoshi

Kashiwagi

Naito

Numerical simulation of wave-induced nonlinear motions of a two-dimensional floating body by the moving particle semi-implicit method.

J Mar Sci Technol 2008; 13: 85–94.

10.

Jena

Biswal

KC.

A numerical study of violent sloshing problems with modified MPS method.

J Hydrodyn Ser B 2017; 29: 659–667.

11.

Tsubota

KWSK.

A particle method for blood flow simulation: application to flowing red blood cells and platelets. J Earth Simul 2006; 5: 2–7.

12.

Sueyoshi

A study of nonlinear fluid phenomena with particle method (part 1). J Kansai Soc Naval Architects 2001; 236: 191–198.

13.

Sueyoshi

A study of nonlinear fluid phenomena with particle method (part 2). J Kansai Soc Naval Architects 2002; 236: 181–186.

14.

Hibi

Yabushita

A study on reduction of unusual pressure fluctuation of MPS method. J Kansai Soc Naval Architects Jpn 2004; 125–131.

15.

Khayyer

Gotoh

Modified moving particle semi-implicit methods for the prediction of 2D wave impact pressure.

Coastal Eng 2009; 56: 419–440.

16.

Khayyer

Gotoh

Enhancement of stability and accuracy of the moving particle semi-implicit method.

J Comput Phys 2011; 230: 3093–3118.

17.

Kondo

Koshizuka

Improvement of stability in moving particle semi-implicit method.

Int J Numer Meth Fluids 2011; 65: 638–654.

18.

Lee

Park

Kim

et al . Step-by-step improvement of MPS method in simulating violent free-surface motions and impact-loads. Comput Methods Appl Mech Eng 2011; 200: 1113–1125.

19.

Tamai

Koshizuka

Least squares moving particle semi-implicit method.

Comp Part Mech 2014; 1: 277–305.

20.

Daneshvar

Rakhshandehroo

Talebbeydokhti

New modified gradient models for MPS method applied to free-surface flow simulations.

Appl Ocean Res 2017; 66: 95–116.

21.

Lucy

LB.

A numerical approach to the testing of fission hypothesis.

Astronom J 1977; 82: 1013.

22.

Gingold RA and

Monaghan JJ.

Smoothed particle hydrodynamics: theory and application to non-spherical stars. Mon Not R Astron Soc 1977; 181: 375–389.

23.

Brookshaw

Method of calculating radiative heat diffusion in particle simulations.

Publ Astron Soc Aust 1985; 6: 207–210.

24.

Cummins

Rudman

An SPH projection method.

J Comput Phys 1999; 152: 584–607.

25.

Lee

Moulinec

et al . Comparisons of weakly compressible and truly incompressible algorithms for the SPH mesh free particle method. J Comput Phys 2008; 227: 8417–8436.

26.

Shibata

Masaie

Kondo

et al . Improved pressure calculation for the moving particle semi-implicit method. Comp Part Mech 2015; 2: 91–108.

27.

Kashiwagi

A CIP-based method for numerical simulations of violent free-surface flows.

J Mar Sci Technol 2004; 9: 143–157.