Enhancement of stabilization of MPS to arbitrary geometries with a generic wall boundary condition

ARTICLE IN PRESS

JID: CAF

[m5G;October 8, 2018;9:14]

Computers and Fluids 0 0 0 (2018) 1–25

Contents lists available at ScienceDirect

Computers and Fluids journal homepage: www.elsevier.com/locate/compfluid

Enhancement of stabilization of MPS to arbitrary geometries with a generic wall boundary condition Tiangang Zhang a,b, Seiichi Koshizuka c, Ping Xuan d,∗, Jinbao Li d,∗, Cheng Gong b a

Electronic Engineering College, Heilongjiang University, Harbin 150080, China School of Mathematical Science, Heilongjiang University, Harbin 150080, China c Graduate School of Engineering, The University of Tokyo, 7-3-1, Hongo, Bunkyo-ku, Tokyo 113-8656, Japan d School of Computer Science and Technology, Heilongjiang University, Harbin 150080, China b

a r t i c l e

i n f o

Article history: Received 11 March 2018 Revised 4 August 2018 Accepted 12 September 2018 Available online xxx Keywords: Moving particle semi-implicit Boundary conditions Numerical oscillation Gradient model Incompressible fluid

a b s t r a c t As a Lagrangian meshless method, moving particle semi-implicit (MPS) method has been proven useful in the analysis of the free-surface flow, especially accompanied by the large deformation and fragmentation of fluids. The improvement of pressure distribution in three dimensions is an important aspect to verify the effectiveness of MPS. The accurate representation of 3-D geometries especially complex geometries is the premise of obtaining convincing pressure distribution. However, most of MPS applications cannot accurately represent complex wall geometries, which highly affects the reliability of MPS. For this reason, the triangle meshes are used to accurately represent the complex wall geometries in this research. The polygon wall boundary condition (PW) is adopted to enforce the wall boundary condition to the triangle meshes. The pressure of the wall boundaries is derived from the Neumann boundary condition to improve the velocity distribution of fluid particles near the wall boundaries. A first-order gradient model is presented to improve the accuracy and stability of the PW. Our approach can enhance the numerical stabilization to arbitrary geometries. We simulate several 3-D examples such as the classic hydrostatic simulation and the complex 3-D geometries with sharp angles and curved surfaces to demonstrate the general applicability of our new model. © 2018 Elsevier Ltd. All rights reserved.

1. Introduction Moving particle semi-implicit (MPS) [1] method is a Lagrangian, mesh-free method developed by Koshizuka and Oka in 1996 to simulate incompressible fluids with free surface. Since the MPS is based on the fully Lagrangian description, the governing equations in MPS are transformed into the interactions between particles, and the convective term instigating the numerical instabilities disappears. The free surface of fluids is also be simply judged by the concentration extent of the particles. Thus, the large deformation of fluids can be easily tracked, which is an advantage over the grid-based methods. Another mesh-free method that is analogous to MPS is the smoothed particle hydrodynamics (SPH) [2,3]. SPH was originally proposed to simulate the problem of astronomy. Afterwards, it was employed in the simulation of the compressible fluids. Monaghan [4] extended SPH to the incompressible fluids so-called weakly compressible SPH (WCSPH). Then Shao and Lo [5] developed the

∗

Corresponding authors. E-mail address: [email protected] (J. Li).

fully incompressible SPH (ISPH) to simulate the incompressible fluids. The wall boundary conditions are important in MPS and SPH. Since the fluid and wall boundaries usually have different characteristic, the coupled accuracy between them is important to improve the overall numerical stability. Generally, five major wall boundary conditions are available in MPS and SPH: •

•

•

•

Arranging some layers of particles that is supposed to have the same characteristic as fluid particles to represent the wall boundaries so-called the dummy particles [1,6]. Using the mirror reflection technique to reflect the fluid particles to the outside of the wall boundaries to supplement the deficiency of the particle number density of fluid particles, i.e. the mirror particles [7]. The boundary particles are not necessary in this boundary condition. The repulsive force method [4]. Employing Lennard-Jones potential to enforce the repulsive force from the boundary particles to the fluid particles. A single layer of boundary particles is used to represent the wall boundaries. The unified semi-analytical wall boundary condition (USAW) [8]. The USAW calculates the contribution of the wall

https://doi.org/10.1016/j.compfluid.2018.09.008 0045-7930/© 2018 Elsevier Ltd. All rights reserved.

Please cite this article as: T. Zhang et al., Enhancement of stabilization of MPS to arbitrary geometries with a generic wall boundary condition, Computers and Fluids (2018), https://doi.org/10.1016/j.compfluid.2018.09.008

JID: CAF 2

ARTICLE IN PRESS

[m5G;October 8, 2018;9:14]

T. Zhang et al. / Computers and Fluids 000 (2018) 1–25

•

Fig. 1. Distance relationship between ith fluid particle and wall boundaries.(For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 2. Initial set-up of the hydrostatic simulation.

Fig. 3. Pressure profile of MPS-PW-IS and PMPS at different to the theoretical values at t = 2.0 (s).(For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

boundaries with the boundary integrals, and it also uses a single layer of boundary particles to represent the wall boundaries. The polygon wall boundary condition (PW) [9]. The geometries are represented by triangle or quadrangle meshes to accurately represent the wall boundaries. The boundary calculations are transformed to the distance function to improve the efficiency of simulations. The boundary particles are also not necessary.

The most widely used wall boundary condition is the dummy particles that was firstly proposed by Koshizuka and Oka [1] when they invented the MPS method. Since sharing the same features with the fluid particles, the dummy particles can be calculated in the same manner as the fluid particles. Due to this merit, a serial of new methods have been proposed to improve the pressure distribution based on the dummy particles. Khayyer and Gotoh [10] derived a symmetric gradient model to enhance the conservation of the linear and angular momentums and suppress the pressure oscillation. Tanaka and Masunaga [11] presented a new source term that coupled the variation of the particle number density with the divergence of velocity so-called the quasi-compressible model to suppress the pressure oscillation. This method can dramatically improve the pressure distribution. Kondo and Koshizuka [12] proposed a general source term that is composed of three parts including one main and two error compensation terms to suppress the pressure oscillation. Khayyer and Gotoh [13] developed a higher order Laplacian model to improve the accuracy of the pressure calculation. Khayyer et al. [14] simulated the interaction between the fluid particles and elastic body with the dummy particles. Another improvement of MPS is presented by Lee et al. [15] who tested different combinations of the collision coefficients and found the proper range of the collision coefficients to improve the spatial stability. Gotoh and Khayyer [16] systematically reviewed the latest advancements of the particle methods from accuracy, stability and conservation properties. The dummy particles are also a convenient wall boundary condition in SPH method. In WCSPH, Morris et al. [6] used the dummy particles to represent 2-D curved surfaces. Adam et al. [17] presented a simple formulation to improve the pressure distribution of WCSPH based on the balance of the boundary forces. In this method, the slip and no-slip boundary conditions can be enforced. In ISPH, Hosseini and Feng [18] presented a new pressure Poisson equation to suppress the pressure oscillation by deriving a rotation projection scheme which can impose non-homogeneous Neumann boundary condition. Hu and Adam [19] enforced the densityinvariance and divergence free conditions during the calculation of the pressure Poisson equation and obtained stable pressure distribution. However, this method is a little time-consuming. Asai et al. [20] adopted the similar scheme as Tanaka and Masunaga [11] to suppress pressure oscillation of ISPH and acquired stable pressure field. Khayyer et al. [21] investigated the energy conservation property in the different MPS and ISPH models and indicated that the corrected Taylor series consistent pressure gradient models provided best energy conservation. Aly et al. [22] imposed the Neumann boundary condition to the dummy particles by solving the linear system of the pressure Poisson equation in ISPH. Xu et al. [23] proposed a particle adjustment method that can slightly change the position of the fluid particles to avoid non-uniform particle distributions. This method is simple and effective to improve the pressure distribution. Lind et al. [24] generalized the particle adjustment method and it can also obtain smooth pressure distribution. Khayyer et al. [25] compared two particle shifting techniques and proposed an optimized particle shifting scheme. Since the dummy particles can be calculated in the same manner as the fluid particles, many improvements of MPS and SPH are based on this wall boundary condition. However, despite these merits, Li et al. [26] compared the variation of the particle number density

Please cite this article as: T. Zhang et al., Enhancement of stabilization of MPS to arbitrary geometries with a generic wall boundary condition, Computers and Fluids (2018), https://doi.org/10.1016/j.compfluid.2018.09.008

JID: CAF

ARTICLE IN PRESS

[m5G;October 8, 2018;9:14]

T. Zhang et al. / Computers and Fluids 000 (2018) 1–25

3

Fig. 4. Pressure time histories of MPS-PW-IS and PMPS against theoretical values at point A.

of the fluid particles at different positions near the wall boundaries composed of the dummy particles and pointed out that the dummy particles produced lots of discontinuous “bumps” along the wall boundaries due to the shape of the boundary particles. These “bumps” causes the non-physical movements of the fluid particles close to the wall boundaries. Furthermore, it is difficult for the dummy particles to precisely represent the complex wall boundaries such as the sharp angles and curved surfaces, which affects the reliability of the simulations with the dummy particles. Since the dummy particles must be the same as the fluid particles, the calculation of the wall boundaries is also expensive, especially for large scale simulations. Due to these restrictions, the dummy particles are only suitable for 2-D simulations and 3-D simulations with simple geometries. Mirror particles are another popular wall boundary condition. Since the reflective process does not affect the representation of the wall boundaries, the physical boundaries can be accurately represented by the polygonal meshes. Akimoto [27] simulated 2-D fall wedge on the free surface with the mirror particles. The reflective process near the corners is simplified by using the dividing line of the angles to avoid overlarge wall contributions to the fluid particles. However, this simplification can only be used in two dimensions. Park and Jeun [28] used mirror particles in 2-D problems to reduce the computing load by the dummy particles. Lee et al. [15] also employed mirror particles to enforce slip and no-slip boundary conditions to improve the pressure distribution. Yildiz et al. [29] proposed a multiple boundary tangent method in SPH to simulate two-dimensional cylinder obstacle and obtained stable results. However, this technique is so time consuming that it is difficult to be employed in three dimensions. Marrone et al. [30] used the concept of the mirror particles to generate layers of the dummy particles and enforced the Neumann boundary condition. This method can obtain smooth pressure distribution. However, the generation process must be carefully designed according to the angle of the geometries to avoid excessive or deficient contribution of the wall boundaries. Thus, this method can only be extended to simple 3-D problems. Manuel et al. [31] presented a new open boundary condition for ISPH with mirror particles to simulate two-dimensional micro flow and obtained good pressure distribution. In sum, mirror particles can maintain the physical shape of the geometries and obtain stable pressure distribution. In addition, the slip and no-slip boundary conditions can be easily enforced. However, due to the complexity of mirror reflective technique, the mirror particles can only be used in 2-D and simple 3-D geometries.

The repulsive force method [4] was proposed by Monaghan. It is a convenient wall boundary condition that can treat 3-D complex geometries. The wall boundaries are discretized by a single layer of boundary particles and the interactions between the fluid and boundary particles are imposed by the Lennard-Jones potential. Since only one layer of boudary particles are employed, the complex geometries can be simply represented. However, the interaction between the fluid and boundary particles is separated, which causes the pressure oscillation in the simulations of complex geometries. To address this issue, Monaghan and Kajtar [32] proposed the boundary force integrals method to suppress the spurious behaviors caused by the repulsive force. However, this method is lack of the unified formulations, which restricts the range of applications. The USAW was proposed by Ferrand et al. [33] to accurately represent the wall boundaries and improve the pressure distributions in 2-D. Leroy et al. [34] developed a turbulent model and then proposed a new open boundary condition [35] based on the USAW. These two methods can both obtain smooth pressure distributions. Aly et al. [22] used the analytical kernel function to accurately compute the renormalization factor and its gradient term. Mayrhofer et al. [36] extended the USAW to three dimensions. However, the accurate boundary integrals are substituted by the approximate equations to improve the computational efficiency, which leads to the accuracy of the USAW dramatically decreases. Thus, the USAW can attain smooth pressure distributions in 2-D. However, it is difficult to be employed in 3-D simulations. The PW was proposed by Harada et al. [9] in 2008 to solve the problem that the dummy particles cannot accurately represent the complex geometries. The geometries are represented by the triangle meshes in PW, which is the same as the mirror particles. A staggered background mesh is used in PW. The shortest distance from a single grid node to the triangle meshes is calculated and stored. The distance from a fluid particle to the wall boundaries is calculated by the linear interpolation of the distance stored at grid nodes around this fluid particle. In 2-D, four and in 3-D, eight surrounding grid nodes are used to interpolate the distance of the fluid particles. The Laplacian and gradient models between the fluid particles and wall boundaries are turned into the distance function. Since the boundary particles are not involved in the calculation of the Laplacian and gradient models, the PW can effectively improve the computational efficiency in large scale simulations. However, the wall contributions to the fluid particles near the corners are not accurate due to the elimination of the boundary particles. To address this issue, Zhang et al. [37]

Please cite this article as: T. Zhang et al., Enhancement of stabilization of MPS to arbitrary geometries with a generic wall boundary condition, Computers and Fluids (2018), https://doi.org/10.1016/j.compfluid.2018.09.008

JID: CAF 4

ARTICLE IN PRESS

[m5G;October 8, 2018;9:14]

T. Zhang et al. / Computers and Fluids 000 (2018) 1–25

Fig. 5. Snapshots of hydrostatic simulations of MPS-PW-IS and PMPS at t = 0.23, 0.3, 0.4, 0.5, 0.7 and 2.0 (s).

proposed a simple formulation to calculate the wall contribution to the fluid particles and suppress the strong pressure oscillation. Although this method can improve the particle number density to a certain extent, it is still not accurate near the corners. Then Zhang et al. [38] presented a boundary particle arrangement technique to improve the particle number density near the corners. To suppress the pressure oscillation, Zhang et al. [39] further proposed a new source term for the PW, namely the MPS-PW-IS method to suppress the pressure oscillation. This method can improve the pressure distribution to arbitrary geometries. However, the pressure calculation near the wall boundaries is still not accurate. To the simulations of the complex geometries, the arbitrary Lagrangian Eulerian (ALE) [40] approach can be used to transform the grid points to conform to the complex wall boundaries. However, it is limited to moderate body deformation because the mesh is required to conform to the body all the times. To the non-boundary conforming methods, the Cartesian (or cut cell)

methods [41,42] can be used in two and three dimensional problems. Another grid based method is the immersed boundary conditions [43,44], which uses the fixed Cartesian grid to simulate the complex wall boundaries and moving bodies. Although the Cartesian method and immersed boundary conditions can simulate complex wall boundaries, the efficiency of these two methods is not high, which can be hardly used in the large scale industry simulations with the complex wall geometries. In the particles methods, the simulations of the complex wall boundaries are restricted by the enforcement of the wall boundary conditions. Most of the researches adopted the dummy particles to treat the complex wall geometries due to the convenient imposition of the wall boundary conditions. Shibata et al. [45] simulated three dimensional shipping water on the deck comparing to the experiments. The smooth vessel body is represented by the dummy particles. Tartakovsky and Meakin [46] used SPH to mimic the phenomenon of the flow into the fracture and complex

Please cite this article as: T. Zhang et al., Enhancement of stabilization of MPS to arbitrary geometries with a generic wall boundary condition, Computers and Fluids (2018), https://doi.org/10.1016/j.compfluid.2018.09.008

JID: CAF

ARTICLE IN PRESS

[m5G;October 8, 2018;9:14]

T. Zhang et al. / Computers and Fluids 000 (2018) 1–25

5

Fig. 5. Continued

apertures represented by the dummy particles. Cleary et al. [47] applied SPH to modelling the die casting process with the dummy particles. Ulrich et al. [48] outlined the application of SPH to the interaction between floating bodies, water/oil and ship propeller flows. Although the dummy particles are employed in the simulations of the complex geometries, the error between the physical boundaries and the numerical boundaries still exists, which influences the accuracy of the simulations. To represent complex wall geometries, Monaghan and Kos [49] employed the repulsive force method to represent the slope of the beach and simulated the process of a run-up of a solitary wave onto the beach. However, the pressure oscillation is strong near the slope. To treat complex geometries in industry simulations with efficiency, the PW gradually becomes a convenient boundary condition to simulate the interactions between the fluid and complex wall boundaries in recent years. Mitsume et al. [50] used the PW to simulate the interaction between the fluid and structure. Regmi

et al. [51] mimicked the filling and solidification process of the box cavity composed of the aluminum alloy with the PW. Yuhashi et al. [52] employed the PW to simulate the stirring of the fluid in the rotational cam-shaft comparing to the experiments. The PW is imposed to the cam-shaft constituted by the triangle meshes. Tanaka et al. [53] proposed a multi-resolution MPS method with the PW to simulate the complex wall geometries. In this work, the PW is adopted to guarantee that 3-D geometries can be accurately represented. The problem of the gradient model used in MPS-PW-IS is elaborated. To improve the accuracy of the gradient model, we derive the pressure of the wall boundaries from the Neumann boundary condition. The first order gradient model proposed by Suzuki [54] is introduced to the PW to represent the interaction between fluid particles. The interaction between the fluid particles and wall boundaries in the gradient model is derived from the first order gradient model [54] and the pressure of the wall boundaries. The new gradient model can

Please cite this article as: T. Zhang et al., Enhancement of stabilization of MPS to arbitrary geometries with a generic wall boundary condition, Computers and Fluids (2018), https://doi.org/10.1016/j.compfluid.2018.09.008

JID: CAF 6

ARTICLE IN PRESS

[m5G;October 8, 2018;9:14]

T. Zhang et al. / Computers and Fluids 000 (2018) 1–25

Fig. 5. Continued

improve the velocity calculations of the fluid particles near the wall boundaries and further improve the pressure distributions. The new gradient model coupled with the improved source term [39] forms a new method so-called PMPS method. The PMPS is simple and can improve the pressure distribution to arbitrary geometries. We test PMPS with several 3-D cases. At first, we validate the hydrostatic simulation with a rectangular and a cylindrical water tank to show the improvement of the pressure calculation of PMPS. Then, the classic dam break problem is simulated to show that PMPS can effectively suppress the pressure oscillation. Subsequently, the simulation results are compared with the experiment of the dam break with an obstacle. At last, two geometries with sharp and obtuse angles, and a complex geometry are simulated to show the improvement of PMPS to arbitrary geometries. 2. MPS-PW-IS method In this section, MPS-PW-IS method is briefly introduced. The detailed descriptions are illustrated in [39]. The governing equations for the incompressible viscous flows are the mass conservation and Navier–Stokes equations

Dρ = 0, Dt

(1)

1 Du = − ∇ P + ν∇ 2 u + g, Dt ρ

(2)

whereρ is the density, t is the time, u is the velocity vector, P is the pressure, ν is the kinematic viscosity and g is the acceleration due to gravity. The geometry shapes are composed of the triangle meshes and a staggered background mesh is utilized in the MPS-PW-IS. The shortest distance and normal vector of each grid node are calculated and stored in the memory before the main loop. If the grid point is in the flow region, the shortest distance is designated as positive. If it is on the wall boundaries, the distance is set as zero. Otherwise, it is negative. The normal vector is calculated by the central difference of the shortest distance stored at the grid points. The distance and normal vector are only calculated for one time. The shortest distance of the fluid particle is calculated by the trilinear interpolation of the distances stored at the grid points around this fluid particle. Three layers of boundary particles are generated by the boundary particle arrangement (BPA) [38] technique to supplement the particle number density of the fluid particles. The pressure Poisson equation is represented by the variation of

Please cite this article as: T. Zhang et al., Enhancement of stabilization of MPS to arbitrary geometries with a generic wall boundary condition, Computers and Fluids (2018), https://doi.org/10.1016/j.compfluid.2018.09.008

JID: CAF

ARTICLE IN PRESS

[m5G;October 8, 2018;9:14]

T. Zhang et al. / Computers and Fluids 000 (2018) 1–25

(a) Time histories of free surface against theoretical values.

(b) Heights of free surface at different X position at t=2 (s).

Fig. 6. Variation of heights of free surface of MPS-PW-IS and PMPS.

Fig. 7. (a) Initial configuration of cylinder composed of triangle meshes, (b) Initial arrangement of fluid particles in cylinder.

Fig. 8. Pressure time histories of free falling of fluid particles in cylinder at t = 0.3, 0.5 and 0.7 (s).

Please cite this article as: T. Zhang et al., Enhancement of stabilization of MPS to arbitrary geometries with a generic wall boundary condition, Computers and Fluids (2018), https://doi.org/10.1016/j.compfluid.2018.09.008

7

JID: CAF 8

ARTICLE IN PRESS

[m5G;October 8, 2018;9:14]

T. Zhang et al. / Computers and Fluids 000 (2018) 1–25

Fig. 9. (a) Simulation result of PMPS shown by black-edged circles at t = 3.0 (s), (b) Simulation result of PMPS shown by areas at t = 3.0 (s).

the particle number densities

 2  n∗i f ρ nk − n∗ ρ n0 − nski }, ∇ P i = ∗ { (1 − γ ) 2 i 0 i + γ 2 ni t n t n0

(3)

where nki is the particle number density of ith fluid particle at kth time step, n0 is the particle number density constant, n∗i is the particle number density of ith fluid particle at intermediate step, n∗i f is the fluid contributions to the particle number density. Parameter γ is chosen in the range of 0.01 < γ < 0.05. nski is the weighted average of the particle number density of ith particle at kth time step Fig. 10. Initial set-up of dam break simulation.

 

 nski

j=i, j∈ f luid

= 

nkj w(ri j  )

j=i, j∈ f luid

  .

w(ri j  )

(4)

where w is the weight function defined as

Fig. 11. Comparison of pressure time histories between MPS-PW-IS and PMPS at points A and B.

Please cite this article as: T. Zhang et al., Enhancement of stabilization of MPS to arbitrary geometries with a generic wall boundary condition, Computers and Fluids (2018), https://doi.org/10.1016/j.compfluid.2018.09.008

ARTICLE IN PRESS

JID: CAF

[m5G;October 8, 2018;9:14]

T. Zhang et al. / Computers and Fluids 000 (2018) 1–25

9

Fig. 12. Snapshots of dam break simulation with MPS-PW-IS and PMPS from the side view at t = 0.43 and 0.7 (s).

⎧   ⎨  re  − 1 (0 ≤ ri j  < re )       ri j w ( ri j ) = ,   ⎩0 (re ≤ ri j  )

(5)

ith particle. The wall part becomes a distance function

where |rij | is the distance between ith and jth particles, re is the effective radius. The fluid particles are initially arranged with the equal distance L0 that is the diameter of the fluid particles. In this research, re is chosen as 3.1L0 in all the simulations. The gradient term is divided into the fluid and wall parts. The fluid part uses the strictly derived gradient model proposed by Khayyer and Gotoh [10]

∇ Pi f

d = 0 n j=i



Pi + Pj − 2Pˆi

  r j − ri 2

   (r j − ri )w(r j − ri  ) ,

where Pi , Pj are the pressure of ith and jth fluid particles, respectively, Pˆi is the minimum pressure within the effective radius re of

(6)

∇ Piw =

ρ|dr| riw , t 2 |riw |

(7)

where riw = rw − ri is the shortest distance vector from ith fluid particle to the wall boundaries, |dr|=L0 − |riw |. Fig. 1 shows the position relationship between riw and dr. dr is the distance vector represented by the blue arrow pointing to the balance position of fluid particle to keep the distance from ith fluid particle to the wall boundaries equal to L0 . riw is the green arrow pointing to the wall boundaries. The black circle w on the wall boundaries is the closest distance point of ith fluid particle to the wall boundaries. Since the closest distance point is only related to the position of the coun-

Please cite this article as: T. Zhang et al., Enhancement of stabilization of MPS to arbitrary geometries with a generic wall boundary condition, Computers and Fluids (2018), https://doi.org/10.1016/j.compfluid.2018.09.008

JID: CAF 10

ARTICLE IN PRESS

[m5G;October 8, 2018;9:14]

T. Zhang et al. / Computers and Fluids 000 (2018) 1–25

Fig. 13. Pressure distributions of the dam break simulation with MPS-PW-IS and PMPS from the bottom view at t = 0.5, 0.6, 0.72 and 0.8 (s).

terpart fluid particle, we do not need to calculate the position of this point. riw can be simply obtained by multiplying the normal vector by the closest distance. If |riw | is greater than or equal to L0 , |dr| is zero. Then ∇ Piw is zero, which means the wall boundaries do not enforce repulsive force to the fluid particle. If |riw | is less

than L0 , ∇ Piw takes effect and the wall boundaries will push the fluid particles back to the balance position. The accuracy of |riw | is affected by the grid size. Generally, the fine grid can obtain more accurate interpolation distance than the coarse grid. However, fine grid will also slow down the computational efficiency of the algo-

Please cite this article as: T. Zhang et al., Enhancement of stabilization of MPS to arbitrary geometries with a generic wall boundary condition, Computers and Fluids (2018), https://doi.org/10.1016/j.compfluid.2018.09.008

JID: CAF

ARTICLE IN PRESS

[m5G;October 8, 2018;9:14]

T. Zhang et al. / Computers and Fluids 000 (2018) 1–25

11

Fig. 14. Comparison of particle distributions between MPS-PW-IS and PMPS at t = 2.2 (s). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

the number of neighboring particles and formulated as

Ni =



ws (ri j ),

(9)

j=i

wherews is formulated as



w ( ri j ) = s

1 ( ri j ≤ re ) . 0 ( ri j > re )

(10)

To improve the spatial stability, the collision coefficients are chosen the same as [15]. 3. Current challenge in MPS-PW-IS Although the pressure distribution can be improved by MPSPW-IS, the gradient model still encounters a problem i.e. the wall part of the gradient model is invariable no matter what gradient model is selected. To illustrate this problem, we derive the wall part of the gradient model from the beginning. In MPS [1], the lefthand side of Navier–Stokes equation is divided into two parts

Fig. 15. Schematic description of leading edge.

rithm at the initial procedure. For balance, the grid size is chosen equal to L0 in all simulations in this research. The surface detection uses the method proposed by Tanaka and Masunaga [11]

Ni < β N0 ,

(8)

where N0 is the maximum number of neighboring particles within effective radius for fully submerged particles in the initial distribution, β is the experience coefficient which is assigned as 0.8, Ni is

uk+1 − u∗ u∗ − uk 1 + =− ∇P t t ρ

k+1

+



k ν∇ 2 u + gk ,

(11)

whereu∗ is the temporary velocity, superscript k and k + 1 represent kth and (k + 1)th time step. Split (11) into two equations

u∗ = uk + t u = −

 2 k ν∇ u + g ,

t k+1 ∇P , ρ

(12)

(13)

Please cite this article as: T. Zhang et al., Enhancement of stabilization of MPS to arbitrary geometries with a generic wall boundary condition, Computers and Fluids (2018), https://doi.org/10.1016/j.compfluid.2018.09.008

ARTICLE IN PRESS

JID: CAF 12

[m5G;October 8, 2018;9:14]

T. Zhang et al. / Computers and Fluids 000 (2018) 1–25

Since dr and riw satisfy dr/|dr| = −riw /|riw | as shown in Fig. 1, (18) can be further formulated as

Pw =

ρ n0 |dr||riw | + 2Pˆi − Pi . dt 2 n∗iw

(19)

Substitute (19) back into (16) and then the wall part of the gradient model namely (7) can be obtained. When different fluid part of the gradient models is chosen, uif should be different. The wall part of the gradient model should also produce corresponding uiw so as to obtain consistent velocity distribution near the wall boundaries. However, the wall part of the gradient model is invariable to any gradient model, which instigates the pressure oscillation near the wall boundaries in MPS-PW-IS. 4. PMPS method To treat the issue of the gradient model in MPS-PW-IS, we present a new gradient model in this part. The pressure of the wall boundaries is derived from the Neumann boundary condition and the wall part of the gradient model can change with the fluid part of the gradient model. The Neumann boundary condition takes the form

∂ y (x ) = f (x ) ∂n

(20)

where n denotes the normal vector to the wall boundaries; y(x) is a smooth function and f(x) is a scalar function. The normal derivative on the left side of (20) is defined as

∂ y (x ) = ∇ y (x )  n ∂n

Fig. 16. Comparison of position change of leading edge of collapsed water column.

where u =uk+1 − u∗ is the corrective velocity. (13) can be further divided into the fluid and wall parts

ui f = −

t ∇ Pi f , ρ

(14)

t ∇ Piw , ρ

(15)

uiw = −

∇ Piw =

d Pi + Pw − 2Pˆi ∗ niw riw , n0 |riw |2

(16)

where Pw is the pressure of the shortest distance point of ith fluid particle, n∗iw is the wall contribution of ith fluid particle to the particle number density. To calculate Pw , substitute (16) into (15)

uiw = −

dt Pi + Pj − 2Pˆi ∗ niw riw . ρ n0 |riw |2

(17)

Multiply t on both side of (17), dr can be derived as following

dt Pi + Pw − 2Pˆi ∗ dr = − niw riw . ρ n0 |riw |2

where ∇ y(x) represents the gradient vector of y(x). In PW, the most important step is to derive the pressure of the wall boundaries, specifically, the pressure value of the shortest distance point of each fluid particle. Here we will derive the pressure of the wall boundaries not from the gradient model but from the Neumann boundary condition. Multiply t on both sides of (15), and transform uiw to the distance vector

dr = −

where uif is the velocity of ith particle caused by surrounding fluid particles, uiw is the velocity of ith particle yielded by the wall boundaries, and u = ui f + uiw . ∇ Pif can be solved in the same manner as (6). The interaction between the fluid particles and wall boundaries can be regarded as the interaction between the fluid particles and their closest distance points. Thus, the summation in (6) disappears

2

(18)

(21)

t 2 ∇ Piw . ρ

(22)

Then, dot product of (22) with n on both sides

dr · n = −

t 2 ∇ Piw · n, ρ

(23)

where n is the normal vector of the wall boundaries to ith fluid particle as shown in Fig. 1. Since every surface mesh has a normal vector, there will be many normal vectors in the effective radius of ith fluid particle. To decide the direction of n, we may regard n as the opposite direction of the strongest repulsive force of the wall boundaries to ith fluid particle. Since riw is the direction vector of the shortest distance, it can be regarded as the opposite direction of the strongest repulsive force. Thus, the direction of n can be replaced by the normalization of riw . Replace n on the left side of (23) with riw /|riw | and we can get

dr ·

riw t 2 =− ∇ Piw · n. ρ |riw |

(24)

The right side of (24) satisfies the definition of the Neumann boundary condition (21). Thus, (24) can be further transformed into

dr ·

riw

|riw |

=−

t 2 ∂ P , ρ ∂n

(25)

where∂ P/∂ n is the normal derivative of pressure. Extend ∂ P/∂ n along the normal direction of riw with the first order difference

Please cite this article as: T. Zhang et al., Enhancement of stabilization of MPS to arbitrary geometries with a generic wall boundary condition, Computers and Fluids (2018), https://doi.org/10.1016/j.compfluid.2018.09.008

ARTICLE IN PRESS

JID: CAF

[m5G;October 8, 2018;9:14]

T. Zhang et al. / Computers and Fluids 000 (2018) 1–25

13

Fig. 17. Leading edge simulation between MPS-PW-IS and PMPS from the side view at t = 0.6, 0.8, 0.9, 1.05, 1.1, 1.4 (s).

approximation

dr ·

riw

|riw |

=−

t 2 Pw − Pi . ρ |riw |

(26)

Thus, Pw can be yielded

Pw = Pi +

ρ |dr||riw |. t 2

(27)

As we can see, the pressure value of the wall boundaries is different from original wall pressure (19). Then substitute (27) into (16) and a new formulation of ∇ Piw can be attained. However,

using this wall part of the gradient model, the fluid particles penetrate the wall boundaries and flowing out. The reason is that the fluid part of the gradient model (6) used in MPS-PW-IS is lack of the accuracy comparing to Pw derived from the first order approximation. Thus, u cannot be accurately calculated near the wall boundaries, which leads to the leakage of the fluid particles. To solve this problem, the gradient model with the first order accuracy should be considered. We test different gradient models and finally choose the first order gradient model proposed by Suzuki [54] to balance the forces between the fluid particles and wall boundaries. The gradient model proposed by Suzuki [54] is

Please cite this article as: T. Zhang et al., Enhancement of stabilization of MPS to arbitrary geometries with a generic wall boundary condition, Computers and Fluids (2018), https://doi.org/10.1016/j.compfluid.2018.09.008

ARTICLE IN PRESS

JID: CAF 14

[m5G;October 8, 2018;9:14]

T. Zhang et al. / Computers and Fluids 000 (2018) 1–25

Fig. 17. Continued

written as



∇ Pi =

 1  w ( r j − ri  ) n0



r −r

 j i r j − ri 

r −r

 j i r j − ri  j=i

   1 Pj − Pi r j − ri    w(r j − ri  ) , r j − ri  r j − ri  n0 

the fluid particles. The formulation of ∇ Pif is the same as (28). ∇ Piw can also be derived by replacing j with w, w(|r j − ri | ) with w(|rw − ri | ) and removing the summation in (28). Then the new wall part of gradient model can be formulated as

−1



(28)

j=i

where symbol “” represents the tensor product. ∇ Pi can also be divided into ∇ Pif and ∇ Piw . ∇ Pif is the interaction between

∇ Piw =

riw riw  |riw | |riw |

−1 

 ρ|dr| riw . t 2 |riw |

(29)

When the distance from the fluid particles to the wall boundaries is less than L0 (29) is functioned. Comparing to (28), (riw /|riw |)(riw /|riw |) is generally irreversible. However, we may

Please cite this article as: T. Zhang et al., Enhancement of stabilization of MPS to arbitrary geometries with a generic wall boundary condition, Computers and Fluids (2018), https://doi.org/10.1016/j.compfluid.2018.09.008

JID: CAF

ARTICLE IN PRESS

[m5G;October 8, 2018;9:14]

T. Zhang et al. / Computers and Fluids 000 (2018) 1–25

15

Fig. 18. Leading edge simulation between MPS-PW-IS and PMPS from the bottom view at t = 0.6, 0.8, 0.9, 1.05, 1.1, 1.4 (s).

calculate the pseudoinverse of (riw /|riw |)(riw /|riw |) although some accuracy is lost. Fortunately, the accuracy loss is slight and the simulation results can be improved after testing different cases. The merit of (29) is that it is simple and can be extended to arbitrary geometries without special treatment. Thus, it can be eas-

ily employed in 3-D complex geometries. The purpose of inverse of the tensor product is to calculate the repulsive direction of the wall boundaries to the fluid particles to attain accurate velocity distributions and improve the pressure distributions.

Please cite this article as: T. Zhang et al., Enhancement of stabilization of MPS to arbitrary geometries with a generic wall boundary condition, Computers and Fluids (2018), https://doi.org/10.1016/j.compfluid.2018.09.008

ARTICLE IN PRESS

JID: CAF 16

[m5G;October 8, 2018;9:14]

T. Zhang et al. / Computers and Fluids 000 (2018) 1–25

By solving the pseudoinverse of tensor product, ∇ Piw can be calculated and the result is the same as (7). It manifests that (29) can be simplified to (7) when the wall boundaries are flat, which is consistent with the physical sense of the repulsive force. Thus, the inverse of the tensor product in (29) only adjusts the direction of repulsive force near the corners to avoid over strong repulsive force from the polygonal surfaces to the fluid particles. The MPS-PW-IS with the new gradient model (28) and (29) is dubbed as the polygon MPS (PMPS) method. 5. Numerical results and discussions In this section, we present four test cases including the hydrostatic, dam break, dam break with a rectangle obstacle and three simulations of complex geometries to show the improvement of PMPS against MPS-PW-IS. 5.1. Hydrostatic simulation

Fig. 19. The schematic dam break with an obstacle.

To the flat wall boundaries as shown in Fig. 1, the repulsive force from the wall boundaries to ith fluid particle should be vertical to the wall boundaries and points inside to the flow domain. Since riw /|riw | is the normal vector that is vertical to the flat wall boundaries, we can directly designate the normal vector of the wall boundary as riw /|riw | = (0, −1, 0 ) . Substitute riw /|riw | into (29) and yield

 ∇ Piw =

0 −1 0



 

0 −1 0

−1 

 ρ|dr| riw . t 2 |riw |

(30)

To illustrate the performance of PMPS, we first apply it to the hydrostatic simulation. The initial set-up of the hydrostatic simulation is shown in Fig. 2. The length, width, and height of the water tank are 1.08 (m), 1.08 (m) and 1.64 (m), respectively. The x-direction is designated to align with the length direction. The height of the water column is 1.0 (m) with 18252 fluid particles. The background mesh is arranged to coincide with the wall boundaries of the water tank since the geometry shape is simple. Three layers of boundary particles are generated uniformly along the wall boundaries with the BPA technique. The simulation conditions are shown in Table 1. The equations used in MPS-PW-IS and PMPS are listed in Table 2. For the PMPS, the coefficient γ in (3) is chosen equal to 0.02 after a serial of tests to guarantee that the pressure values of fluid particles near the bottom surface are closest to the theoretical value. To the

Fig. 20. Snapshots of PMPS compared with experiment at t = 0.4 (top) and 0.56 (bottom) (s).

Please cite this article as: T. Zhang et al., Enhancement of stabilization of MPS to arbitrary geometries with a generic wall boundary condition, Computers and Fluids (2018), https://doi.org/10.1016/j.compfluid.2018.09.008

JID: CAF

ARTICLE IN PRESS

[m5G;October 8, 2018;9:14]

T. Zhang et al. / Computers and Fluids 000 (2018) 1–25

17

Fig. 21. Pressure time histories at P1 (left) and P2 (right).

Fig. 22. Initial set-up of water tank with triangular and trapezoidal wedges.

Table 1 Calculation conditions. Initial condition

Value

Unit

Diameter of fluid particle L0 Effective radius re Density ρ Collision coefficient α 1 Collision coefficient α 2 Gravity acceleration

0.04 3.1L0 10 0 0 0.9 0.2 9.81

m m kg/m3

m/s2

Table 2 Equations used in MPS-PW-IS and PMPS.

MPS-PW-IS PMPS

Fig. 23. Initial arrangement of fluid particles in the water tank with the triangular and trapezoidal wedges.

Fluid part of gradient model

Wall part of gradient model

Eq. (6) Eq. (28)

Eq. (7) Eq. (29)

MPS-PW-IS, γ used in [39] is 0.01. Thus, in this study, we respectively test the hydrostatic pressure with γ = 0.01 and 0.02 in MPSPW-IS to find the proper coefficient. Fig. 3 shows the hydrostatic

Please cite this article as: T. Zhang et al., Enhancement of stabilization of MPS to arbitrary geometries with a generic wall boundary condition, Computers and Fluids (2018), https://doi.org/10.1016/j.compfluid.2018.09.008

JID: CAF 18

ARTICLE IN PRESS

[m5G;October 8, 2018;9:14]

T. Zhang et al. / Computers and Fluids 000 (2018) 1–25

Fig. 24. Snapshots of the geometry with the triangular and trapezoidal wedges at t = 2.0 (s), (a) Simulation result of PMPS shown by black-edged circles, (b) Simulation result of PMPS shown by areas.

Fig. 25. Initial set-up of water tank with two triangular wedges.

pressures of MPS-PW-IS and PMPS to the theoretical values at different depths at t = 2.0 (s). The black line is the theoretical values at different depths. The red spots and green triangles are the pressure values calculated by MPS-PW-IS with γ = 0.01 and 0.02, respectively. The blue squares are the pressure values calculated by PMPS with γ = 0.02. The pressure values calculated by PMPS are in good agreements with the theoretical values. The pressure value calculated by MPS-PW-IS with γ = 0.01 is almost the same as the theoretical values at 0.2 (m). However, with the depths increasing, the errors gradually increase. Although the error decreases at 1.0 (m), the pressure value is still less than PMPS and theoretical value. When γ is chosen as 0.02 in MPS-PW-IS, the pressure values are a little larger than those with PMPS at 0.2 and 0.4 (m), and almost the same as those with PMPS at 0.6 and 0.8 (m). However, at 1.0 (m), the pressure value is far greater than the theoretical value and that of PMPS. Comparing to MPS-PW-IS, the PMPS can attain better pressure calculations. Although MPS-PW-IS with

γ = 0.02 has a better agreement with theoretical values than that with γ = 0.01 at most depths, the pressure near the wall boundaries with γ = 0.02 is overlarge than the theoretical value, which dramatically affects the accuracy of the pressure calculations near the wall boundaries and induces pressure oscillations. Thus, we choose γ = 0.01 in MPS-PW-IS in all simulations to avoid overlarge pressure calculations near the boundaries. Fig. 4 shows the time histories of pressure between MPS-PWIS with γ = 0.01 and PMPS at a fixed measuring point A which is 0.04 (m) above the bottom of the water tank. When simulation begins, the fluid particles move downwards due to gravity. Since the distance between the fluid particles and bottom surface decreases, the repulsive force of the bottom surface is imposed to the fluid particles and pushes the fluid particles to move upwards. Thus, the pressure at point A increases at first and then decreases. After 0.7 (s), the fluid particles gradually reach the balance status. As shown in Fig. 4, the pressure calculated by PMPS smoothly

Please cite this article as: T. Zhang et al., Enhancement of stabilization of MPS to arbitrary geometries with a generic wall boundary condition, Computers and Fluids (2018), https://doi.org/10.1016/j.compfluid.2018.09.008

JID: CAF

ARTICLE IN PRESS T. Zhang et al. / Computers and Fluids 000 (2018) 1–25

Fig. 26. Initial arrangement of fluid particles in the water tank with two triangular wedges.

increases and decreases at first 0.4 (s), whereas the pressure by MPS-PW-IS increases smoothly but decreases with the obvious pressure oscillation. After 0.4 (s), the pressure values calculated by PMPS are closer to the theoretical values than those by MPS-PW-IS as shown in the enlarged portion of Fig. 4. Thus, PMPS can obtain more accurate pressure calculations than MPS-PW-IS. Fig. 5 shows the snapshots of the hydrostatic simulation results of MPS-PW-IS with γ = 0.01 and PMPS at different time steps. The color represents the pressure values. At first 0.7 (s), the obvious pressure oscillation can be seen near the all boundaries in MPSPW-IS due to the inaccurate repulsive force calculated by the wall part of the gradient model. After 0.7 (s), the stable pressure distribution can be acquired. However, the fluid particles are clustered in MPS-PW-IS at t = 2.0 (s). The distance between each pair of fluid particles is less than L0 as shown in the magnified portion of MPSPW-IS. To the PMPS, the pressure distributions are still smooth at each time step and the clustering phenomenon does not occur. It illustrates that the wall part of the gradient model in PMPS are consistent with the fluid part of the gradient model, which improve the velocity field near the boundaries and overall pressure distributions.

[m5G;October 8, 2018;9:14] 19

Fig. 6 shows the variation of heights of the free surface simulated by PMPS and MPS-PW-IS with γ = 0.01 to validate the conservation of the fluid. Fig. 6(a) shows the time histories of peaks of free surface of MPS-PW-IS and PMPS comparing to the theoretical values. The free surface of PMPS becomes stable after 1.0 (s). Then the heights of the free surface gradually increase and become almost the same as the theoretical values after 1.25 (s). The free surface of MPS-PW-IS becomes stable after 0.7 (s) and the heights of it is much lower than the theoretical values. At each time step, the volume conservation of PMPS is better than that of MPS-PW-IS. Fig. 6(b) shows the heights of the free surface of MPSPW-IS and PMPS at different x-coordinate at 2.0 (s). The heights of the free surface by PMPS are closer to the theoretical values than those by MPS-PW-IS. The root mean squares (RMS) of MPS-PW-IS and PMPS with respect to hi -1 where hi is the depth of the surface particles are 0.064 and 0.0198, respectively. The error of PMPS is much smaller than that of MPS-PW-IS. Thus, the overall volume conservation of PMPS is better than that of MPS-PW-IS. Since MPSPW-IS causes the clustering of the fluid particles, the conservation of the volume of MPS-PW-IS must be worse than that of PMPS. To show the improvement of PMPS to the curved surfaces, the hydrostatic simulation in a cylinder is tested. Fig. 7(a) depicts the configuration of the cylinder which is composed of 4016 triangle meshes. The height of the cylinder is 1.6 (m). Since the cylinder is complex, the background mesh cannot align with the complex wall boundaries. Practically, the background mesh is randomly arranged without considering the position relationship to the cylinder. Fig. 7(b) shows the initial arrangement of the fluid particles in the cylinder and 11894 fluid particles are used in the simulation. Due to the restriction of the curved wall boundaries, the fluid particles can also not be arranged along the wall boundaries with the equal distance L0 to the wall boundaries. In this simulation, the fluid particles are arranged in the center of the cylinder hovering above the bottom surfaces. The distance from the fluid particles to all wall boundaries are greater than 2L0 at the initial procedure. Fig. 8 depicts the pressure distributions of the fluid particles in cylinder at three different time steps. When the simulation begins, the fluid particles descend and collide with the wall boundaries, and then fill the empty space. On account of the repulsive force of the wall boundaries, the fluid particles then move upwards. Fig. 8(a) exhibits the pressure increase process while the fluid particles falling on the bottom surfaces. Fig. 8(b) shows the

Fig. 27. Snapshots of the geometry with two triangular wedges at t = 2.0 (s), (a) Simulation result of PMPS shown by black-edged circles, (b) Simulation result of PMPS shown by areas.

Please cite this article as: T. Zhang et al., Enhancement of stabilization of MPS to arbitrary geometries with a generic wall boundary condition, Computers and Fluids (2018), https://doi.org/10.1016/j.compfluid.2018.09.008

JID: CAF 20

ARTICLE IN PRESS

[m5G;October 8, 2018;9:14]

T. Zhang et al. / Computers and Fluids 000 (2018) 1–25

Fig. 28. Initial configuration of complex geometry.

Fig. 29. Initial arrangement of fluid particles in the complex geometry.

moving up process of the fluid particles pushed by the wall boundaries. After 0.7 (s), the fluid particles gradually become stable. The pressure distributions at three time steps are all smooth even the drastic collision of the fluid particles on the bottom wall boundaries at 0.3 (s). Fig. 9 shows the pressure distribution of PMPS at 3.0 (s). The color represents the pressure values. The fluid is represented by the particles and area in Fig. 7(a) and (b), respectively. The pressure distribution of PMPS is smooth, which illustrates that the PMPS can obtain stable pressure distribution to the geometries with the curved surfaces. 5.2. Dam break simulation In this part, the dam break simulation is conducted to show the improvement of PMPS to suppress the pressure oscillation. The simulation conditions are shown in Table 3. The initial distance between particles is 0.04 (m). Fig. 10 shows the initial set-up of a 3-D dam break example. The length, width, and height of the water tank are 2.4 (m), 2.4 (m) and 1.6 (m), respectively. The height of water column is 1.6 (m) and the total number of fluid particles

Table 3 Calculation conditions. Initial condition

Value

Unit

Effective radius re Density ρ Collision coefficient α 1 Collision coefficient α 2 Gravity acceleration γ in MPS-PW-IS γ in PMPS

3.1L0 10 0 0 0.9 0.2 9.81 0.01 0.02

m kg/m3

m/s2

used in the simulation is 160 0 0. Points A and B are chosen to test the pressure time histories of MPS-PW-IS and PMPS. Point A is at the center of the right edge of the bottom surface and B is vertical to A. The distance between A and B is 0.04 (m). The background mesh is also aligned with the wall boundaries. Fig. 11 shows the pressure time histories of MPS-PW-IS and PMPS calculated at points A and B. At first 1.0 (s), the pressure values of PMPS are on the top of those of MPS-PW-IS at both points A and B because the pressure values calculated by MPS-PW-IS is

Please cite this article as: T. Zhang et al., Enhancement of stabilization of MPS to arbitrary geometries with a generic wall boundary condition, Computers and Fluids (2018), https://doi.org/10.1016/j.compfluid.2018.09.008

JID: CAF

ARTICLE IN PRESS T. Zhang et al. / Computers and Fluids 000 (2018) 1–25

[m5G;October 8, 2018;9:14] 21

Fig. 30. Pressure distributions of the complex geometry at t = 1.0, 1.5, 4.0, 5.0 and 9.0 (s).

less than those by PMPS against the theoretical values as depicted in Fig. 3. However, the pressure fluctuation scope of MPS-PW-IS is greater than that of PMPS at first 1.0 (s) and the pressure oscillation of MPS-PW-IS is much stronger than PMPS in the whole 3 (s). Thus, PMPS can better suppress the pressure oscillation than MPS-PW-IS. Fig. 12 shows the snapshots of the dam break simulation between MPS-PW-IS and PMPS at 0.43 and 0.7 (s). From the side view, both methods can acquire smooth pressure distribution. However, comparing to MPS-PW-IS, PMPS can obtain better pressure distributions especially near the wall boundaries as shown in the amplified portion at 0.43 (s) and another rectangle box with black edges. At 0.7 (s), the pressure oscillation of MPS-PWIS can also be seen in the enlarged portion and another rectan-

gle box. The PMPS can obtain better pressure distributions than MPS-PW-IS. Since PMPS mainly improves the pressure distributions near the wall boundaries, the accuracy of the gradient model can be more clearly shown under large pressure. Thus, we compare the pressure distribution of the dam break simulation from the bottom view. Fig. 13 shows the pressure distributions of PMPS and MPSPW-IS from the bottom view at different time steps. In MPS-PW-IS, the pressure oscillation can be seen at each time step. The pressure oscillations are almost all happened near the corners, which is caused by the inaccurate repulsive direction of the non-planar wall boundaries. On the contrary, PMPS can acquire smooth pressure distributions especially near the corners at all time steps. To more explicitly illustrate the improvement of PMPS against MPS-PW-IS, we compare the particle distributions at t = 2.2 (s) as

Please cite this article as: T. Zhang et al., Enhancement of stabilization of MPS to arbitrary geometries with a generic wall boundary condition, Computers and Fluids (2018), https://doi.org/10.1016/j.compfluid.2018.09.008

JID: CAF 22

ARTICLE IN PRESS

[m5G;October 8, 2018;9:14]

T. Zhang et al. / Computers and Fluids 000 (2018) 1–25

Fig. 30. Continued

shown in Fig. 14. The fluid particles are cut apart along the x direction and the back part is reserved. In MPS-PW-IS, the fluid particles near the bottom surface form a boundary layer and the distance gap can be clearly seen in the amplified part of black edge rectangle box. This phenomenon also appears at the free surface as shown in the enlarged part of the red edge rectangle box. Single layer of fluid particles sticks to the back surface of the water tank without falling down by gravity. At the front surface of the water tank, the same problem arises. This is caused by the inaccurate direction of the repulsive force. However, in PMPS, this problem does not appear as shown in the enlarged parts of the red and black rectangle boxes. Thus, PMPS can effectively adjust the direction of the repulsive force of the wall boundaries to improve the pressure distributions near the wall boundaries and further improve the overall pressure calculations.

Fig. 15 shows the initial set up of position change of leading edge. The length, width and the height of the water tank are 8 (m), 2 (m) and 8 (m), respectively. The length and height of the fluid are 2 (m) and 4 (m), respectively. The total number of the fluid particles is 2450 0 0. Fig. 16 shows the position changes of leading edge of MPS, MPS-PW-IS and PMPS. The horizontal axis indicates the nondimensional time t(2 g/L)0.5 , and the vertical axis represents the non-dimensional position of water column’s edge. The simulation results are compared with the experimental results by Martin and Moyce [55]. The MPS-PW-IS has a better agreement with the experimental results than MPS because the wall boundaries represented by the triangle meshes is smoother than those with dummy particles. The positions of the leading edge of PMPS are closer to the experimental results than those of MPS-PW-IS, which

Please cite this article as: T. Zhang et al., Enhancement of stabilization of MPS to arbitrary geometries with a generic wall boundary condition, Computers and Fluids (2018), https://doi.org/10.1016/j.compfluid.2018.09.008

JID: CAF

ARTICLE IN PRESS

[m5G;October 8, 2018;9:14]

T. Zhang et al. / Computers and Fluids 000 (2018) 1–25

23

Fig. 30. Continued

illustrates that PMPS can achieve more accurate velocity distribution than MPS-PW-IS. Fig. 17 shows the snapshots of the leading edge simulation between MPS-PW-IS and PMPS. The color represents the pressure. The pressure distributions of MPS-PW-IS and PMPS are both smooth at 0.6 (s). At other time steps, the slight pressure oscillation can be seen near the corners in MPS-PW-IS. The PMPS can obtain smooth pressure distributions at all times steps. Fig. 18 compares the pressure distributions of the leading edge between MPS-PW-IS and PMPS from the bottom view at the same time steps as those in Fig. 17. The pressure oscillation can be clearly seen in MPS-PW-IS at each time step except t = 0.6 (s). The PMPS acquires better pressure distributions near the wall boundaries than MPS-PW-IS at all time steps. The simulation results manifest that the new gradient model in PMPS can effectively improve the pressure distributions near the polygonal wall boundaries.

5.3. Dam break flow impinging on an obstacle The 3-D dam break experiment was originally performed by Kleefman [56] to evaluate the green water loading on the deck equipment. Here it is adopted to test the performance the PMPS against MPS-PW-IS. The water tank involves a rectangle. The geometry is depicted in Fig. 19. The height of the water tank in the experiment is 1.0 (m) and we use 2.0 (m) in this simulation. The simulation conditions are listed in Table 3 and the initial distance between fluid particles are 0.02 (m). The initial background mesh does not align with the geometry. The locations of the pressure transducers on the obstacle are also illustrated. Two point P1 and P2 is on the right side of the obstacle. The simulation uses 85080 fluid particles. The snapshots of the simulation compared with experiment at t = 0.4 and 0.56 (s) are shown Fig. 20 and good agreements are achieved on the arrival time of the flow front. The height of the flow is lower than the experiment at 0.56 (s) because the height of the water tank is two times higher than that used in

Please cite this article as: T. Zhang et al., Enhancement of stabilization of MPS to arbitrary geometries with a generic wall boundary condition, Computers and Fluids (2018), https://doi.org/10.1016/j.compfluid.2018.09.008

JID: CAF 24

ARTICLE IN PRESS

[m5G;October 8, 2018;9:14]

T. Zhang et al. / Computers and Fluids 000 (2018) 1–25

the experiment. The contours of the pressure distributions at two time steps are also smooth. In Fig. 21, the pressure time histories at P1 and P2 simulated with MPS-PW-IS and PMPS are compared with the experimental data [56]. The pressure values in MPS-PW-IS at P1 from 3.0 to 4.7 (s) are greater than the experimental values and the peak value is less than the experimental data. Although the peak pressure value of PMPS is also less than the experimental data, it is closer to the experimental data than that of MPS-PW-IS. Except the peak value, the pressure time histories of PMS at P1 agree well with the experimental data. At P2 , the PMPS and MPS-PW-IS both have the same trend as the experimental data. The PMPS agrees well with experimental data at peak value. However, the pressure values of PMPS are greater than the experimental data since 1.2 (s). This is because the pressure of P2 is measured by weighted average of pressure, the pressure values are affected by the particle resolution. If the particle resolution is improved, the error should be decreased. Comparing to PMPS, the peak value of MPS-PW-IS is far less than the experimental data and the error of pressure values to the experimental data from 1.2 (s) is also greater than PMPS.

gles, convex and concave surfaces is composed of 8988 triangle meshes. The length, width and height of the geometry are 2 (m), 3.1 (m) and 3.1 (m), respectively. There are 82274 fluid particles arranged in the complex geometry at the initial time as shown in Fig. 29. Since the geometry is very complex, fluid particles cannot be well arranged along the wall boundaries. The position of the background mesh is arranged randomly and has no relationship with the complex geometry. Fig. 30 shows the simulation results at different time steps. The color represents the pressure values. The violent flow of the fluid particles can be seen when the fluid particles are bounced back by the wall boundaries at 1.0 (s). Then the fluid particles fall into the flow at 1.5 (s). Although the flow is violent at 1.0 and 1.5 (s), the pressure distributions are both smooth. At 1.5 (s), slight pressure oscillation can be seen from the side view because the fluid particles fall into the flow and affect the pressure distribution. The free surface does not become stable until 9.0 (s). The pressure distributions are all smooth at each time step especially near the sharp angles. Thus, PMPS can improve the pressure distribution to arbitrary complex geometries with the curved surface.

5.4. Simulation of complex geometries

5. Conclusions

In this part, three complex geometries are simulated to show the improvement of the pressure distribution of PMPS to arbitrary geometries. The first two examples are simulated to show PMPS can improve the pressure distribution to the concave and convex geometries. The third example is to illustrate the improvement of PMPS to enhance the pressure distribution to the complex geometries with curved surfaces and sharp angles. The initial particle distance is 0.04 (m) and the simulation conditions are shown in Table 3. Fig. 22 shows the initial set up of a complex geometry with two wedges. There is one triangular wedge in the right surface of the geometry and one trapezoidal wedge in the bottom surface of the geometry. The length, width, and height of geometry are 2.0 (m), respectively. The geometry composed of 39 surfaces uses 93408 fluid particles in the simulation. The background mesh is arranged to coincide with the left surface of the geometry. The grid nodes do not coincide with the triangle meshes near the wedges. Fig. 23 shows the initial arrangement of the fluid particles in the geometry. Since the geometry is complex, the fluid particles cannot be accurately arranged around the wedges. When the simulation starts, the fluid particles will fall freely by gravity and then fill the space around the wedges. Fig. 24 shows the snapshot of simulation result at 2.0 (s). The color represents the pressure. Although the free surface is still not stable, the pressure distribution of PMPS is smooth even near the two wedges. Fig. 25 describes the geometry composed of 35 triangle meshes with two sharp angles. The length, width, and height of the geometry are also 2 (m), respectively. There are two triangular wedges in the geometry. One is the right triangle in the right surface of the geometry and the other is the isosceles triangle on the bottom surface. There are 90768 fluid particles used in the simulation. The background is also arranged to coincide with the left vertical surface of the geometry. The grid nodes do not coincide with the triangle meshes near the two wedges. Fig. 26 shows the initial arrangement of the fluid particles in the geometry. There are still space existed around the triangular wedges. Fig. 27 shows the snapshots of simulation result at 2.0 (s). The color also represents the pressure. The free surface is still not stable. However, the smooth pressure distribution can also be achieved near sharp angles. Fig. 28 shows the initial set-up of the geometry with complex curved surfaces. The complex geometry including sharp an-

In this research, the polygon wall boundary condition (PW) is introduced to MPS to accurately represent 3-D wall boundaries. The CAD data can be directly utilized without revising, and the numerical boundaries are unified with the physical boundaries. Since the wall boundaries have no relationship with the fluid particles, the larger fluid particles can be used and thereby, the efficiency of simulations can be improved. The wall part of the gradient model of MPS-PW-IS is proved to be invariable no matter what gradient model is adopted. Thus, the corrective velocity of fluid particles caused by the wall part of the gradient model is unchangeable to any gradient model, which disturbs the pressure distribution. To address this issue, we derive the pressure of the wall boundaries from the Neumann boundary condition and a new gradient model for the PW is presented. The new gradient model can attain consistent velocity distribution between the fluid and the wall boundaries and further improve the pressure distributions near the wall boundaries. The new method PMPS can improve the pressure distribution to arbitrary geometries with simple boundary treatments. Thus, it has immense potential in the industry simulations. The PMPS method is validated by the hydrostatic simulation with the straight and curved wall geometries to show the improvement of pressure calculations. The classic dam break simulation shows PMPS can suppress effectively the pressure oscillations. The PMPS agrees better with the experiment than MPS-PW-IS by comparing with the experiment of the dam break with an obstacle. Finally, the simulations of three complex geometries show that PMPS can obtain smooth pressure distributions to arbitrary 3-D geometries. Although the velocity and pressure distributions have been improved in the PW, the velocity boundary condition such as the slip and no-slip boundary conditions does not be imposed due to the shape restriction of the 3-D arbitrary geometries. Thus, next step, we will consider imposing the velocity boundary conditions to the polygon wall boundary condition to further enhance the accuracy of the boundary calculations. References [1] Koshizuka S, Oka Y. Moving-particle semi-implicit method for fragmentation of incompressible fluid. Nucl Sci Eng 1996;123:421–34. [2] Gingold RA, Monaghan JJ. Smoothed particle hydrodynamics: theory and application to non-spherical stars. Mon Not R Astron Soc 1977;181:375–89. [3] Lucy LB. A numerical approach to the testing of the fission hypothesis. Astron J 1977;82:1013–24.

Please cite this article as: T. Zhang et al., Enhancement of stabilization of MPS to arbitrary geometries with a generic wall boundary condition, Computers and Fluids (2018), https://doi.org/10.1016/j.compfluid.2018.09.008

JID: CAF

ARTICLE IN PRESS T. Zhang et al. / Computers and Fluids 000 (2018) 1–25

[4] Monaghan JJ. Simulating free surface flows with SPH. J Comput Phys 1994;110:399–406. [5] Shao S, Lo EY. Incompressible SPH method for simulating Newtonian and non-Newtonian flows with a free surface. Adv Water Resour 2003;26:787–800. [6] Morris JP, Fox PJ, Zhu Y. Modeling low Reynolds number incompressible flows using SPH. J Comput Phys 1997;136:214–26. [7] Cummins SJ, Rudman M. An SPH projection method. J Comput Phys 1999;152:584–607. [8] Ferrand M, Laurence D, Rogers BD, Violeau D, Kassiotis C. Unified semi-analytical wall boundary conditions for inviscid, laminar or turbulent flows in the meshless SPH method. Int J Numer Methods Fluids 2013;71:446–72. [9] Harada T, Koshizuka S, Shimazaki K. Improvement of wall boundary calculation model for MPS method. Trans Jpn Soc Comput Eng Sci 2008 (20080006).(in Japanese). [10] Khayyer A, Gotoh H. Development of CMPS method for accurate water-surface tracking in breaking waves. Coastal Eng J 2008;50:179–207. [11] Tanaka M, Masunaga T. Stabilization and smoothing of pressure on MPS method by Quasi-compressibility. Trans JSCES 2008 20080025. [12] Kondo M, Koshizuka S. Improvement of stability in moving particle semi-implicit method. Int J Numer Methods Fluids 2011;65:638–54. [13] Khayyer A, Gotoh H. A higher order Laplacian model for enhancement and stabilization of pressure calculation by the MPS method. Appl Ocean Res 2010;32:124–31. [14] Khayyer A, Gotoh H, Falahaty H, Shimizu Y, Nishijima Y. Towards development of a reliable fully-Lagrangian MPS-based FSI solver for simulation of 2D hydroelastic slamming. Ocean Syst Eng Int J 2017;7:299–318. [15] Lee B-H, Park J-C, Kim M-H, Hwang S-C. Step-by-step improvement of MPS method in simulating violent free-surface motions and impact-loads. Comput Methods Appl Mech Eng 2011;200:1113–25. [16] Gotoh H, Khayyer A. On the state-of-the-art of particle methods for coastal and ocean engineering. Coastal Eng J 2018;60:79–103. [17] Adami S, Hu X, Adams N. A generalized wall boundary condition for smoothed particle hydrodynamics. J Comput Phys 2012;231:7057–75. [18] Hosseini SM, Feng JJ. Pressure boundary conditions for computing incompressible flows with SPH. J Comput Phys 2011;230:7473–87. [19] Hu X, Adams NA. An incompressible multi-phase SPH method. J Comput Phys 2007;227:264–78. [20] Asai M, Aly AM, Sonoda Y, Sakai Y. A stabilized incompressible SPH method by relaxing the density invariance condition. J Appl Math 2012;2012:1–24. [21] Khayyer A, Gotoh H, Shimizu Y, Gotoh K. On enhancement of energy conservation properties of projection-based particle methods. Eur J Mech Fluids 2017;66:20–37. [22] Nguyen TM, Aly AM, Lee S-W. Improved wall boundary conditions in the incompressible smoothed particle hydrodynamics method. Int J Num Methods Heat Fluid Flow 2018;28:704–25. [23] Xu R, Stansby P, Laurence D. Accuracy and stability in incompressible SPH (ISPH) based on the projection method and a new approach. J Comput Phys 2009;228:6703–25. [24] Lind S, Xu R, Stansby P, Rogers B. Incompressible smoothed particle hydrodynamics for free-surface flows: a generalised diffusion-based algorithm for stability and validations for impulsive flows and propagating waves. J Comput Phys 2012;231:1499–523. [25] Khayyer A, Gotoh H, Shimizu Y. Comparative study on accuracy and conservation properties of two particle regularization schemes and proposal of an optimized particle shifting scheme in ISPH context. J Comput Phys 2017;332:236–56. [26] Li D, Sun Z, Chen X, Xi G, Liu L. Analysis of wall boundary in moving particle semi-implicit method and a novel model of fluid–wall interaction. Int J Comput Fluid Dyn 2015:1–16. [27] Akimoto H. Numerical simulation of the flow around a planing body by MPS method. Ocean Eng 2013;64:72–9. [28] Park S, Jeun G. Coupling of rigid body dynamics and moving particle semi-implicit method for simulating isothermal multi-phase fluid interactions. Comput Methods Appl Mech Eng 2011;200:130–40. [29] Yildiz M, Rook R, Suleman A. SPH with the multiple boundary tangent method. Int J Num Methods Eng 2009;77:1416–38. [30] Marrone S, Antuono M, Colagrossi A, Colicchio G, Le Touzé D, Graziani G. δ -SPH model for simulating violent impact flows. Comput Methods Appl Mech Eng 2011;200:1526–42.

[m5G;October 8, 2018;9:14] 25

[31] Hirschler M, Kunz P, Huber M, Hahn F, Nieken U. Open boundary conditions for ISPH and their application to micro-flow. J Comput Phys 2016;307:614–33. [32] Monaghan JJ, Kajtar JB. SPH particle boundary forces for arbitrary boundaries. Comput Phys Commun 2009;180:1811–20. [33] Ferrand M, Laurence D, Rogers B, Violeau D, Kassiotis C. Unified semi-analytical wall boundary conditions for inviscid, laminar or turbulent flows in the meshless SPH method. Int J Num Methods Fluids 2013;71:446–72. [34] Leroy A, Violeau D, Ferrand M, Kassiotis C. Unified semi-analytical wall boundary conditions applied to 2-D incompressible SPH. J Comput Phys 2014;261:106–29. [35] Leroy A, Violeau D, Ferrand M, Fratter L, Joly A. A new open boundary formulation for incompressible SPH. Comput Math Appl 2016;72:2417–32. [36] Mayrhofer A, Ferrand M, Kassiotis C, Violeau D, Morel F-X. Unified semi-analytical wall boundary conditions in SPH: analytical extension to 3-D. Numer Algorithms 2015;68:15–34. [37] Zhang T, Koshizuka S, Shibata K, Murotani K, Ishii E. Improved wall weight function with polygon boundary in moving particle semi-implicit method. Trans Jap Soc Comput Eng Sci 2015 2015. [38] Zhang T, Koshizuka S, Murotani K, Shibata K, Ishii E, Ishikawa M. Improvement of boundary conditions for non-planar boundaries represented by polygons with an initial particle arrangement technique. Int J Comput Fluid Dyn 2016:1–21. [39] Zhang T, Koshizuka S, Murotani K, Shibata K, Ishii E. Improvement of pressure distribution to arbitrary geometry with boundary condition represented by polygons in particle method. Int J Num Methods Eng 2017. [40] Hirt C, Amsden AA, Cook J. An arbitrary Lagrangian-Eulerian computing method for all flow speeds. J Comput Phys 1974;14:227–53. [41] Noh WF. CEL: A time-dependent, two-space-dimensional, coupled Eulerian-Lagrange code. Lawrence radiation lab. Livermore: University of California; 1963. [42] Almgren AS, Bell JB, Colella P, Marthaler T. A Cartesian grid projection method for the incompressible Euler equations in complex geometries. SIAM J Sci Comput 1997;18:1289–309. [43] Viecelli JA. A method for including arbitrary external boundaries in the MAC incompressible fluid computing technique. J Comput Phys 1969;4:543–51. [44] Viecelli JA. A computing method for incompressible flows bounded by moving walls. J Comput Phys 1971;8:119–43. [45] Shibata K, Koshizuka S, Tanizawa K. Three-dimensional numerical analysis of shipping water onto a moving ship using a particle method. J Marine Sci Technol 2009;14:214–27. [46] Tartakovsky AM, Meakin P. Simulation of unsaturated flow in complex fractures using smoothed particle hydrodynamics. Vadose Zone J 2005;4:848–55. [47] Cleary P, Prakash M, Ha J. Novel applications of smoothed particle hydrodynamics (SPH) in metal forming. J Mater Process Technol 2006;177:41–8. [48] Ulrich C, Leonardi M, Rung T. Multi-physics SPH simulation of complex marine-engineering hydrodynamic problems. Ocean Eng 2013;64:109–21. [49] Monaghan J, Kos A. Solitary waves on a Cretan beach. J Waterway Port Coastal Ocean Eng 1999;125:145–55. [50] Mitsume N, Yoshimura S, Murotani K, Yamada T. Improved MPS-FE fluid-structure interaction coupled method with MPS polygon wall boundary model. Comput Model Eng Sci 2014;101:229–47. [51] Regmi A, Shintaku H, Sasaki T, Koshizuka S. Flow simulation and solidification phenomena of AC4CH aluminum alloy in semi-solid forging process by explicit MPS method. Comput Particle Mech 2015;2:223–32. [52] Yuhashi N, Matsuda I, Koshizuka S. Calculation and validation of stirring resistance in cam-shaft rotation using the moving particle semi-implicit method. J Fluid Sci Technol 2016;11 JFST0018-JFST0018. [53] Tanaka M, Cardoso R, Bahai H. Multi-resolution MPS method. J Comput Phys 2018;359:106–36. [54] Suzuki Y. (2007) Ph.D. thesis, The University of Tokyo, Japan. [55] Martin J, Moyce W. Part IV. An experimental study of the collapse of liquid columns on a rigid horizontal plane, Philos Trans R Soc London A, 244 (1952) 312–324. [56] Kleefsman K, Fekken G, Veldman A, Iwanowski B, Buchner B. A volume-of-fluid based simulation method for wave impact problems. J Comput Phys 2005;206:363–93.

Please cite this article as: T. Zhang et al., Enhancement of stabilization of MPS to arbitrary geometries with a generic wall boundary condition, Computers and Fluids (2018), https://doi.org/10.1016/j.compfluid.2018.09.008