A projection-based particle method with optimized particle shifting for multiphase flows with large density ratios and discontinuous density fields

Accepted Manuscript

A Projection-Based Particle Method with Optimized Particle Shifting for Multiphase Flows with Large Density Ratios and Discontinuous Density Fields Abbas Khayyer , Hitoshi Gotoh , Yuma Shimizu PII: DOI: Reference:

S0045-7930(18)30784-9 https://doi.org/10.1016/j.compfluid.2018.10.018 CAF 4035

To appear in:

Computers and Fluids

Received date: Revised date: Accepted date:

28 February 2018 5 October 2018 15 October 2018

Please cite this article as: Abbas Khayyer , Hitoshi Gotoh , Yuma Shimizu , A ProjectionBased Particle Method with Optimized Particle Shifting for Multiphase Flows with Large Density Ratios and Discontinuous Density Fields, Computers and Fluids (2018), doi: https://doi.org/10.1016/j.compfluid.2018.10.018

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

ACCEPTED MANUSCRIPT

Highlights 

A multiphase projection-based particle method for flows of large density ratios and discontinuous density fields.



A computational algorithm on the basis of optimized particle shifting is



CR IP T

incorporated for interface stabilization. Validations and comparisons are conducted through several multiphase benchmark tests.

An extended version of an iterative shifting scheme is implemented and tested.



The present study portrays the robustness of particle shifting concept for

AN US



AC

CE

PT

ED

M

multiphase particle-based simulations.

ACCEPTED MANUSCRIPT

A Projection-Based Particle Method with Optimized Particle Shifting for Multiphase Flows with Large

CR IP T

Density Ratios and Discontinuous Density Fields

Abbas Khayyer1, Hitoshi Gotoh, Yuma Shimizu

Department of Civil and Earth Resources Engineering, Kyoto University, Katsura

AN US

Campus, Nishikyo-ku, Kyoto 615-8540, Japan

Abstract

A novel projection-based particle method is presented for simulation of multiphase flows

M

characterized by large density ratios and discontinuous density fields at the phase interface. The method considers a multi-fluid continuous system and comprises of a specific

ED

computational algorithm utilizing the recently developed Optimized Particle Shifting (OPS [1]) scheme to maintain the regularity of particles at the phase interface and free-surface. The

PT

method is founded on an improved version of Moving Particle Semi-implicit (MPS [2]) as a projection-based particle method. A set of previously developed improved schemes are also

CE

adopted and hence the proposed method is referred to as improved MPS + OPS. Validations

AC

are made both qualitatively and quantitatively in terms of accuracy, energy conservation properties as well as convergence properties by consideration of several benchmark tests.

Keywords: multiphase flows, particle method, projection method, MPS method, particle shifting, optimized particle shifting

1

Corresponding Author, Email: [email protected], Phone: +81-75-380-3180

ACCEPTED MANUSCRIPT 1. INTRODUCTION

Multiphase flows are widely encountered in engineering and industrial applications. On the other hand, numerical simulations of multiphase flows, especially those characterized by large density ratios, are known as one of the most challenging issues in computational physics due to presence of large and abrupt density drop at the phase interface that would lead to a

CR IP T

mathematical discontinuity of density and accordingly a discontinuous pressure gradient field. The fully-Lagrangian particle methods appear to be well suited for simulation of multiphase flows (e.g. [3-6]) due to their inherent potential robustness and their distinct advantage in tracking moving boundaries [7-10]. Up to now, several attempts have been made to simulate

AN US

multiphase flows in the context of two well-known particle methods, namely MPS (Moving Particle Semi-implicit) method [2] and SPH (Smoothed Particle Hydrodynamics) method [11]. To tackle the density discontinuity issue, several stabilizing approaches are applied. These approaches include density smoothing schemes [12-15], density averaging [16], background

M

pressure [17-19], artificial repulsive pressure force [20] or artificial surface tension force [17,18]. Indeed, such numerical treatments may adversely affect the accuracy as well as the

ED

overall consistency and conservation properties of the numerical method and lead to imprecise reproductions, as in case of reproduced entrapped air in water slamming [21]. Thus, it is

PT

preferable to develop stable and accurate multi-phase numerical methods that would keep the

CE

material discontinuity at the interface and at the same time provide stable and accurate results. Lind et al. [22] proposed a two-phase incompressible-compressible SPH with discontinuous

AC

density at the interface through coupling a projection-based SPH, namely, Incompressible SPH (ISPH) for water and compressible SPH for the air phase. In their approach, pressure was solved explicitly for the air phase using an equation of state providing a pressure boundary condition for the surface particles in the incompressible water phase. The calculation of water phase was conducted via a semi-implicit procedure providing a velocity boundary condition for the compressible air phase particles located in the vicinity of water-air interface. In order to enable stable/accurate multi-phase calculations with a discontinuous material interface, Lind et

ACCEPTED MANUSCRIPT al. [22] applied a well-known and potentially robust particle regularization scheme, namely, Particle Shifting (PS) scheme [23] for the water phase particles. The PS scheme was originally proposed by Xu et al. [24] and then presented in a generalized form by Lind et al. [23] allowing for extended applications or interfacial flows. This scheme is founded on Fick's diffusion law and a Taylor-series expansion for evaluation of particle quantities in new

CR IP T

positions and has been applied in several developed SPH-based methods, including both single-phase simulations (e.g. [25-27]) and multi-phase ones (e.g. [8,17,22,28-30]).

A challenge for application of PS scheme to interfacial flows (including multi-phase or free-surface flows) corresponds to achievement of a precise and physically-consistent

AN US

elimination of shifting normal to the interface (as highlighted by Lind et al. [23]; Khayyer et al. [31]). In order to achieve an accurate and consistent implementation of particle shifting for free-surface flows, Khayyer et al. [1] proposed a so-called Optimized Particle Shifting (OPS) scheme through a careful implementation of only tangential shifting for free-surface and free-

M

surface vicinity particles. In contrast to PS, the OPS is free from any tuning parameters for free-surface and was shown to consistently result in perfect elimination of shifting normal to

ED

an interface, providing accurate results and highly regular particle distributions.

PT

This paper aims at presenting a novel projection-based particle method for reliable simulation of multiphase flows characterized by large density ratios and discontinuous density

CE

fields at the phase interface. The proposed method considers a multi-fluid integrated system and comprises of a specific computational algorithm utilizing the OPS scheme both at the free-

AC

surface and phase interface. The considered projection-based particle method corresponds to an improved version of MPS method which benefits from several enhanced schemes including HS (Higher order Source term of Poisson Pressure Equation), HL (Higher order Laplacian), ECS (Error Compensating Source term), GC (Gradient Correction) and DS (Dynamic Stabilization, merely for PPP method). Detailed information on the mentioned schemes can be found in recent papers by Khayyer et al. [31] as well as Gotoh and Khayyer [32]. The proposed multiphase improved MPS method will be simply referred to as improved MPS + OPS, as it

ACCEPTED MANUSCRIPT benefits from several improved schemes as well as particle regularization at the phase interface or free-surface through a specific computational algorithm utilizing the OPS scheme. As shown by Khayyer et al. [1], the OPS scheme is capable of regularizing the free-surface and its vicinity particles as well as minimizing the unphysical particle perturbations in a more accurate, conservative and physically consistent manner with respect to the PS scheme [23].

CR IP T

Indeed, a crucial issue in particle-based simulations of multiphase flows characterized by large density ratios and discontinuous density fields corresponds to preservation of the regularity of particle distributions at the phase interface in a physically consistent manner such that overall conservation features and reproduced physics are not adversely affected. The presented

AN US

improved MPS + OPS method takes advantage of OPS scheme for regularization at the phase interface with careful implementation of a two-step shifting algorithm such that regular particle distributions are maintained more consistently and conservatively, with respect to previously developed multiphase shifting approaches [8,17,22,28-30] founded on the concept of PS

M

scheme [23]. In addition, another distinct feature of improved MPS + OPS is that the method is free of any tuning parameters for interface particle shifting.

ED

The paper is organized as follows. In section 2, the multiphase improved MPS + OPS will

PT

be described. This section comprises of descriptions of PS and OPS schemes in the MPS framework as well as explanation of the algorithm of multiphase improved MPS with interface

CE

regularization on the basis of the OPS scheme. In section 3, first, a simple validation is conducted for the considered MPS-based expressions of PS and OPS schemes for a single-

AC

phase fluid flow to portray the robustness of both PS and OPS schemes and attained enhancements by the OPS scheme in the MPS framework, since PS or OPS have not yet been implemented in the MPS context, to the best knowledge of authors. Then, a set of validations and comparative investigations will be presented through a set of multiphase benchmark tests. In the forth section, an extended version of shifting scheme, the so-called iterative shifting scheme, originally proposed by Vacondio and Rogers [33], will be implemented for the

ACCEPTED MANUSCRIPT improved MPS + OPS, resulting in improved MPS + iterative OPS. The fifth section of this paper is dedicated to concluding remarks.

2. IMPROVED MPS + OPS

CR IP T

2.1. PS and OPS schemes in MPS context The PS (Particle Shifting) scheme is proposed by Lind et al. [23] as an effective particle regularization scheme in SPH context. Recently, Khayyer et al. [1] suggested an enhanced version of PS, referred to as OPS (Optimized Particle Shifting), for accurate particle shifting at and in vicinity of free-surface in free-surface fluid flows. The OPS was proposed in the

AN US

context of ISPH (Incompressible SPH). In this section, the OPS scheme is formulated in MPS context, through considerations of concept of particle number density n, and radius of influence re.

M

The PS scheme is founded on Fick's law of diffusion to shift particles from regions of high concentration to those of low concentration. The Fick's law is expressed as: (1)

ED

J  DC

PT

where J is the flux, C denotes the concentration gradient, and D stands for a diffusion coefficient. Accordingly, a particle shifting displacement vector is written as: (2)

CE

 ri   DCi

AC

where D (= D Δt) signifies the shifting coefficient and Ci is the particle concentration at target particle i, which is defined as follows: Ci 

1 n0

w

ij

j i

; Ci 

1 n0

w

ij

j i

;n

w

ij

(3)

j i

where j represents a typical neighboring particle, n0 stands for the initial particle number density (or the particle number density corresponding to perfectly regular particle distributions for a particle with full compact support) and wij denotes the kernel function.

ACCEPTED MANUSCRIPT Considering the von Neumann stability analysis of advection-diffusion equation, the following equation is derived (Lind et al. [23]; Khayyer et al. [1]):

 ri  Cshift d 02 Ci

(4)

where r denotes position vector, ri signifies the particle shifting displacement vector for target particle i, Cshift ≤ 0.5 is a constant number and d0 stands for the particle diameter. The

CR IP T

upper limit of shifting distance is set as 0.1re (re = radius of influence) equivalent to 0.2h in ISPH context as considered by Lind et al. [23] and Khayyer et al. [1].

In order to achieve a reliable physical particle shifting at and in vicinity of free-surface

AN US

particles, the particle shifting displacement vector for these particles is modified as follows [1]:

~ n ~ C  ri  Cshift d 02 I  n i i i

(5)

follows for target particle i:

M

where n~ is the corrected local normal unit vector to the free-surface, which is defined as

1

ED

 ~   Bi  Ci ; B   1 n rij  wij   i i  Bi  Ci  n0 j i 

(6)

PT

The above modification ensures accurate particle regularizations and eliminates shifting

CE

normal to the interface. After shifting particles, the velocity field is updated by considering the following Taylor-series based approximation:

AC

ui  ui  ui   rii

(7)

where u symbolizes the velocity vector, i and i’ correspond to position of target particle i before and after shifting, respectively, and rii’ denotes the relative position vector between post-shifting and pre-shifting positions of particle i.

ACCEPTED MANUSCRIPT 2.2. Multiphase shifting algorithm This section presents details of implementation of shifting schemes for multiphase improved MPS + OPS method. As previously stated, throughout the manuscript the multiphase improved MPS with phase interface/free-surface particle regularization on the basis of OPS will be simply referred to as improved MPS + OPS.

CR IP T

The algorithm of multiphase OPS is illustrated in Fig. 1. At each time step, MPS-based calculation is conducted for the whole domain. Then, particle shiftings are performed for a consistent and accurate particle regularization to ensure the stability and accuracy. The particle shifting process is composed of a two-step procedure. First, particle shifting is carried

AN US

out only for the heavy phase while the light phase is ignored. In this case, the OPS scheme is applied at the phase interface since it is being treated as a free-surface considering that the light phase particles are temporarily ignored. Then, the shifting of light phase is conducted and at this stage heavy phase particles are also considered. Hence, light phase particle shifting

M

at the phase interface is simply performed by applying the PS scheme. The described two-step particle shifting process ensures that there will be no shifting normal to the phase interface

ED

and at the same time guarantees an almost perfect regularized particle distribution in this region. Fig. 2 presents an illustration of the described two-step shifting process. Up to now, a

PT

number of multiphase shifting methods have been proposed [8,17,22,28-30] and in particular

CE

some similarities are observed in between the multiphase shifting method SA2 by Mokos et al. [17] and the method proposed in this paper. However, the proposed improved MPS + OPS

AC

method is characterized by a more careful and physically consistent implementation of particle shifting at the phase interface by incorporation of the OPS scheme and is free of tuning parameters for shifting at interfaces, as highlighted in the introduction section. Similar to [1], the divergence of position vector,   r , is considered for assessment of particle type, i.e. 1) free-surface (or phase interface), 2) free-surface or phase interface vicinity, 3) splash and 4) inner particles. In MPS context, the expression of   r ,

ACCEPTED MANUSCRIPT corresponding to that used in ISPH context by Lee et al. [34], Lind et al. [23] and Khayyer et al. [1], is expressed as:   ri 

1 n0

r

ij

j i

 wij

(8)

In this study, a normalized quintic Wendland kernel is applied and   r is equal to 3 for a

normalized quintic Wendland kernel is expressed as:   　　 for  

rij  re

for

rij  re

AN US

4    r rij ij    1  r  1  4 r wij   e e    0　　　　　　　　　　 

CR IP T

particle with a full compact support and regularly distributed neighbors. The applied

(9)

It should be noted that for heavy phase shifting since the light phase is temporarily ignored,

1 n0

  ri 

w

kK

1 n0

ik

; Ci 

r

kK

ik

 wik

1 n0

 w

kK

ik

ED

Ci 

K  k | k  i and rki  re and k  Ωh 

PT

(10)

M

(3) and (8) are modified as follows.

(11)

(12)

CE

where Ωh signifies the heavy phase. Fig. 3 portrays the detailed implementation of particle shifting in multiphase OPS method. The figure clarifies the particle types considered in

AC

implementation of shifting, their corresponding shifting displacement vector as well as details related to assessment of particle type. It should be noted that the   r would be equal to 3.0 for a particle with full compact support and regularly distributed neighbors in the considered framework (Eq. 11), while the corresponding number in SPH/ISPH frameworks [1,23,34] would be obtained as 2.0. Considering the fact that the related value in free surface criteria is set as 1.5 in the SPH/ISPH contexts [1,23,34], a corresponding value of 2.25 is considered in this study.

ACCEPTED MANUSCRIPT

3. VALIDATIONS AND COMPARISONS

The accuracy of the proposed multiphase method, referred to as improved MPS + OPS, is validated comparatively in this section. For comparisons two other versions of multiphase MPS-based methods are considered, namely, multiphase improved MPS implementing the PS

CR IP T

scheme [23] for interface regularization (referred to as improved MPS + PS) and multiphase improved MPS implementing a First-order accurate Taylor series-based Density Smoothing [12] for interface stabilization (referred to as improved MPS + FDS). All the improved MPS schemes, namely, HS, HL, ECS, GC and DS, applied for improved MPS + OPS are also

AN US

applied for other two considered multiphase methods.

The improved MPS + PS applies a two-stage shifting process similar to that described in Fig. 2, with only one difference, i.e. PS being applied in place of OPS for heavy phase particle shifting. During the heavy phase shifting the same treatment for interface described by Lind et

 C  C   s     n   s  n    

ED

 ri  Cshift d 02 

M

al. [23] is considered (Eq. 27 in [23]), i.e.:

(13)

PT

with s and n denoting the tangent and normal vectors to the free-surface, respectively. Both

CE

parameters  and  in Eq. (13) are set as zero. As for improved MPS + FDS, instead of particle shifting-based regularization, a density

AC

smoothing scheme is applied. At each time step, calculated densities are smoothed by using Eq. (14), as presented by Khayyer and Gotoh [12].

i 



    1 wij

jJ

jJ

i



  i  xij  i yij  wij xij yij 

(14)

where xij (= xj  xi) and yij (= yj  yi) are components of relative position vector rij and spatial variations of density at target particle i are evaluated as follows:

ACCEPTED MANUSCRIPT  i  i rij  i xij  i  i rij  i yij   ;   xij rij xij rij rij yij rij yij rij rij

; rij  rij

(15)

For validations and comparisons four benchmark tests, namely, a simple multiphase perturbation test, oscillating multi-fluid concentric elliptical regions [20], evolution of a multifluid square patch of fluid [12,35] and multi-phase standing graving wave [22] are considered.

CR IP T

In addition, prior to conducting these multiphase benchmark tests, a single-phase test case, namely, evolution of a square patch of fluid [1,35], is performed in section 3.1 since PS or OPS have not yet been implemented in the MPS context, to the best knowledge of authors. In all performed simulations, the calculation time step is set according to CFL stability condition

AN US

(e.g. [34]) and a maximum allowable time step.

3.1 Validation of single-phase improved MPS + OPS

The OPS scheme is originally developed in ISPH context and to the best of our knowledge,

M

there has not been any applications of either PS or OPS in MPS framework. Hence, in this

ED

section a simple validation is made to illustrate the effectiveness of both PS and OPS schemes for a single-phase MPS-based simulation and to portray enhanced performance of OPS with

PT

respect to PS. The evolution of a weightless, initially square patch of inviscid/incompressible single-phase fluid subjected to a pure rigid rotation [1,35] is considered. The length of the

CE

fluid, L, is set as 1.0 m (L = 1.0 m). The initial velocity field is given as:

(16)

AC

u0   y  0    u0      0  v0   x

where angular velocity, ω, is considered as 1.0 rad/s. The particles are 1.0E-2 m in diameter (d0 = 1.0E-2 m). The maximum allowable time step is set as Δtmax = 1.0E-3 s. The density of fluid is set as 1000 kg/m3. The shifting coefficient is set as 0.5 (Cshift = 0.5) [23]. Fig. 4 shows a qualitative comparison in between single-phase improved MPS implementing either OPS or PS scheme. The figure portrays snapshots of particles together

ACCEPTED MANUSCRIPT with normalized pressure field at a normalized time of tω = 2.0, illustrating the effectiveness of both OPS and PS in providing stable and qualitatively accurate pressure distribution. On the other hand, from the enlarged portions of the figure, the free-surface reproduced by OPS is smoother and characterized by more regular particle distribution compared with that by PS. Fig. 5 presents time histories of normalized pressure at the center of square patch by single-

CR IP T

phase improved MPS implementing either OPS or PS. Both time variations of pressure are in relatively good agreement with that of BEM by Colagrossi [35] until t = 2.0, while OPS provides superior performance after t = 2.0. The difference is caused by growth of irregularities in particle distributions at and in vicinity of the free-surface in PS simulation, as

AN US

the simulation proceeds. Accordingly, the global solution of PPE (Poisson Pressure Equation) is affected by free-surface particle irregularities and inaccurate detection of free-surface particles in imposition of free-surface boundary conditions.

M

3.2. A simple multiphase perturbation test

ED

A simple multiphase perturbation test is conducted to verify the effectiveness of OPS for interface particle regularization in presence of perturbations in particle motions. This test was

PT

originally conducted in a single-phase simulation by Khayyer et al. [1]. Fig. 6 shows a schematic sketch corresponding to the initial setup of the calculation. A density ratio of 1 to

CE

1000 is considered, i.e. h /l = 1000. The diameter of particles and the maximum allowable

AC

time step are set as d0 = 1.0E-2 m and Δtmax = 1.0E-4, respectively. An interface heavy phase particle at the center of the tank is perturbed 0.1d0 horizontally at the beginning of calculation at t = 0.

The shifting coefficient is set as 0.01 (Cshift = 0.01) for simulations conducted in section 3.2. The considered shifting coefficient is considerably smaller than the upper limit of 0.5 obtained from a von Neumann stability analysis of the advection-diffusion equation [23]. It should be noted that by combining this estimate with a local CFL condition, Skillen et al. [36] proposed

ACCEPTED MANUSCRIPT an updated shifting coefficient ( Cshift  t ui / d0 ). For the conducted simulations of this study, the computational time step is not only governed by the CFL condition, but also a maximum allowable time step which is set quite small (about one order of magnitude smaller than that obtained from the CFL condition) for multiphase simulations characterized by large density ratios (e.g. [12,15]). Therefore, the shifting coefficients considered for conducted multiphase

CR IP T

flow simulations are about one order of magnitude smaller than that considered for singlephase flow simulations.

Fig. 7 portrays the particle distributions by improved MPS + OPS and improved MPS + PS at t = 30.0 s. The figure illustrates the enhancing effect of OPS in regularization of particle

AN US

distributions in a multi-fluid system of large density ratio. Fig. 8 shows the pressure distribution at the center of the tank (x = 1.0 m) at t = 30.0 s, where the reproduced pressure is shown to be in a good agreement with the theoretical solution. Fig. 9 presents time variation of the x-component of a so-called randomness index  [12] for the perturbed particle. This index

M

was presented by Khayyer and Gotoh [12] as a simple indicator for particle regularity to show

ED

relative enhancements in attaining regular particle distributions. The quantitative value of the index is not normalized and simply approaches to zero as neighbors of a target particle are

PT

more regularly distributed. More rigorous quantifications of randomness (e.g. [37]) should be made in the MPS context for precise evaluations of accuracy, consistency, convergence and

CE

conservation properties. The presented figure illustrates the effectiveness of OPS-based regularization in providing an almost regular particle distribution throughout the simulation of

AC

a simple multi-fluid system of large density ratio.

3.3. Oscillating multi-fluid concentric elliptical regions A two dimensional drop evolving under the action of a central conservative force field is simulated. This test is originally proposed by Monaghan and Rafiee [20] and has been reproduced by Khayyer and Gotoh [12]. The conservative force field corresponds to a potential field of 2r2/2 and is thus equal to 2 r , with  being a dimensional parameter. The

ACCEPTED MANUSCRIPT concentric elliptical regions are of different densities where the lighter fluid surrounds the heavier one. The initial ellipses are circles with radii of 0.5 and 1.0. Both fluids are considered to be incompressible and inviscid. The fluids are subjected to an initial velocity field as:

u0   (0) x  v0   (0) y

(17)

CR IP T

In this study, δ(0) and  are set equal to 0.5 and 1.0, respectively. A density ratio of 1 to 1000 is considered (ρh /ρl = 1000). The particles are 1.0E-2 m in diameter (d0 = 1.0E-2 m). The calculation time step is set according to CFL stability condition and a maximum allowable time step of Δtmax = 1.0E-4 s. The period of the oscillating drop is calculated as T = 2 / 

AN US

4.44 s. The shifting coefficient, Cshift , is set equal to 0.01. Fig. 10 presents a schematic illustration of this test.

Fig. 11 illustrates a set of snapshots of particles together with density field at t = 6.0 s corresponding to the oscillating elliptical regions test by improved MPS methods incorporating

M

either OPS, PS or FDS [12]. From Fig. 11, utilization of OPS has resulted in stable and

ED

qualitatively accurate simulations with almost regular particle distributions at the phase interface without the presence of an unphysical gap. The snapshot by PS is characterized by

PT

irregularities at the phase interface as well as a distinctive unphysical gap in that region. Due to smoothing of density at the phase interface, an unphysical gap does not appear in FDS

CE

snapshot. However, presence of irregularities in particle distributions is evident in FDS

AC

snapshot.

Fig. 12(a) portrays the time histories of the semi-major axis of the outer ellipse. Compared

with FDS and PS schemes, OPS has provided more accurate time histories. The time histories of semi-minor axis of the inner ellipse are shown in Fig. 12(b), where the discrepancies for both PS and FDS results seem to have been intensified. Indeed, application of density smoothing at phase interface has resulted in a clear decrease of the amplitude of oscillations that suggests an unphysical decay of energy by FDS. Table 1 presents the RMSE (Root Mean

ACCEPTED MANUSCRIPT Square Error) and NRMSE (Normalized RMSE) corresponding to the numerical results presented in Fig. 12. To have a more quantitative investigation on energy conservation properties of considered multiphase MPS methods, the time histories of Kinetic Energy (KE) and Potential Energy (PE) are considered and are portrayed in Fig. 13. From this figure, application of FDS has resulted

CR IP T

in clear dissipations of both kinetic and potential Energies. The results by improved MPS + PS scheme, on the contrary, show unphysical increases of both kinetic and potential energies. The improved MPS + OPS, however, portray consistent and acceptable time histories of kinetic and potential energies.

AN US

Fig. 14 provides a comparison of time variations of normalized total energy for the three considered methods and illustrates the enhanced energy conservation properties of improved MPS + OPS. The enhanced energy conservation properties of improved MPS + OPS can be attributed to implementation of a careful and consistent particle regularization at the phase

M

interface. Indeed, further enhancement of the implemented particle regularization at the phase

ED

interface is expected to result in further improved preservation of mechanical energy. Fig. 15 presents the time histories of the semi-major axis of the outer ellipse simulated by

PT

improved MPS + OPS with a set of three different diameters of d0 = 0.008 m, d0 = 0.010 m and d0 = 0.020 m. The RMSE and NRMSE corresponding to the results shown in Fig. 15 are

CE

presented in Table 2. Fig. 15 and Table 2 indicate the possession of an acceptable convergence property by improved MPS + OPS.

AC

A qualitative comparison between improved MPS + OPS and ICSPH or Incompressible-

Compressible SPH method of Lind et al. [22] is presented in Fig. 16. The ICSPH is founded on coupling of ISPH (heavy phase) and Weakly Compressible SPH (light phase), and it benefits from PS implementation for the heavy phase. Both methods have provided stable results with almost negligible noise in pressure field. However, the results by improved MPS + OPS seem to be superior in terms of volume continuity at the phase interface as well as absence of perturbations/noise in velocity field.

ACCEPTED MANUSCRIPT

3.4. Evolution of a multi-fluid square patch Evolution of a weightless multi-fluid square patch subjected to a pure rigid rotation is considered as the second benchmark test. This test is originally proposed by Colagrossi [35] and further extended to multiphase flow simulation by Khayyer and Gotoh [12]. The initial

CR IP T

velocity field is given as:

u0   y  0    u0      0  v0   x

(18)

where angular velocity, ω, is considered as 1.0 rad/s. Here a lighter square-shaped fluid is

AN US

placed inside a heavier one concentrically. The length of the inner fluid and that of the outer one are set as 0.5 m and 1.0 m, respectively. A density ratio of 10 to 1000 is considered (ρh/ρl = 100). The particles are 5.0E-3 m in diameter (d0 = 5.0E-3 m). The maximum allowable time step is set as Δtmax = 2.0E-4 s. The shifting coefficient Cshift is set as 0.02. Both fluids are

M

considered to be inviscid and incompressible. Fig. 17 illustrates a schematic sketch of this

ED

numerical test together with considered computational conditions. Fig. 18 portrays a qualitative comparison in between FDS, PS and OPS, in reproducing the

PT

multi-fluid square patch test. The figure shows snapshots of particles together with density field at a normalized time of tω = 0.44. The snapshot by FDS is characterized by anisotropic

CE

particle distributions due to absence of particle shifting in a purely Lagrangian simulation. The superiority of implementation of OPS in place of PS scheme in a multiphase flow simulation is

AC

clearly seen at the corners of the inner fluid. Fig. 19 illustrates a set of snapshots of particles together with density and pressure fields by

improved MPS + OPS and improved MPS + PS methods at tω = 2.04. The figure highlights the superiority of improved MPS + OPS in providing continuous pressure and volume fields as well as a smooth phase interface with minimized unphysical perturbations. It should be highlighted here that physically consistent minimization of unphysical perturbations at the

ACCEPTED MANUSCRIPT phase interface would have a significant contribution in enhancing the stability and accuracy of multiphase projection-based particle methods. Similar to the previous benchmark tests, the snapshot by multiphase PS is characterized by presence of particle irregularities and more distinctively an unphysical one-particle-gap, at the phase interface similar to those seen in single-phase PS results of free-surface flows [23]. This gap is seen in heavy phase since this

CR IP T

phase is considered in the first stage of two-stage shifting process and shifting normal to the first layer of heavy phase is set as zero through the application of Eq. 13 (or Eq. 27 in [23]).

Fig. 20 presents the time histories of normalized kinetic energy of the patch by improved MPS

implementing

either

OPS,

PS

or

FDS

for

phase

interface

particle

AN US

regularization/stabilization. From the presented figure, the improved MPS + OPS has resulted in a smooth numerical dissipation of kinetic energy of the patch. The results by other two considered methods are inferior to that by improved MPS + OPS in terms of showing more

M

dissipations as well as unsmooth variations.

ED

3.5. Multi-phase standing gravity wave

An oscillating standing gravity wave is simulated based on the conditions described in the

PT

paper by Lind et al. [22]. Fig. 21 presents a schematic sketch of this test together with the considered computational conditions. Here, the system comprises of a heavy fluid (water, h =

CE

1000 kg/m3) positioned beneath a light one (air, l = 1 kg/m3). The water depth h and the bottom width are set 0.5 m and 1.0 m, respectively. The periodic boundary condition is

AC

enforced in the horizontal direction (at the side walls) and no penetration conditions are imposed at top and bottom walls. Initial profile of water surface is given as:

η0 ( x)  Acos kx

(19)

where η0 is the initial surface elevation, A signifies wave amplitude ( = 0.05 m), k represents wave number ( = 2π/λ), and λ denotes wave length ( = 1.0 m). The diameter of particles is set

ACCEPTED MANUSCRIPT as d0 = 5.0E-3, and the maximum allowable time step is considered as Δtmax = 5.0E-5 s, respectively. The shifting coefficient is set as 0.025 (Cshift = 0.025). Fig. 22 illustrates the snapshots of particles together with the density field at t = 2.0 s by improved MPS methods incorporating OPS, PS and FDS schemes. The improved MPS + OPS method has resulted in a clear and almost symmetric air-water interface characterized by

CR IP T

relatively regular particle distributions. The snapshot by improved MPS + PS indicates the presence of an almost one-particle-thick gap at the water phase. The worst snapshot seems to correspond to improved MPS + FDS, in terms of absence of a clear air-water interface and irregularities at both phase interface and inside the air phase. Fig. 22 portrays the effectiveness

AN US

of particle shifting concept in achieving enhanced multiphase flow simulation results.

Fig. 23 shows the time histories of water surface elevation measured at the center of the tank by the three considered multiphase improved MPS methods. In this figure, the produced water surface elevations are compared with the second-order theory of Wu and Taylor [38].

M

The time history provided by improved MPS + OPS is shown to be in a better agreement with the theoretical solution compared with those by improved MPS + PS and improved MPS +

ED

FDS.

PT

Fig. 24 presents the time histories of water surface elevation at the center of the tank by improved MPS + OPS with a set of three different particle diameters of d0 = 1.0E-2 m, d0 =

CE

5.0E-3 m and d0 = 4.0E-3 m. The RMSE and NRMSE corresponding to the results shown in Fig. 24 are presented in Table 3. Fig. 24 and Table 3 reconfirm the possession of an acceptable

AC

convergence property by improved MPS + OPS. Besides, with respect to Table 3, the relation between particle diameter and obtained NRMSE is plotted in Fig. 25. From this figure, the convergence rate of improved MPS + OPS is shown to be almost first order.

ACCEPTED MANUSCRIPT 4. IMPROVED MPS + ITERATIVE OPS

The proposed improved MPS + OPS method comprises of a two-step shifting process, where heavy phase particles are shifted first (ignoring the light phase) and then shifting of light phases are carried out. Indeed, from Fig. 2, it is evident that a consistent shifting of light phase requires a precise and consistent shifting of the heavy phase at the first step of two-step shifting

CR IP T

process. Hence, the overall accuracy of the particle shifting process implemented in improved MPS + OPS significantly depends on the accuracy of heavy phase shifting by OPS.

Recently, Vacondio and Rogers [33] proposed an iterative particle shifting scheme, which employs an iterative algorithm for particle shifting in order to enhance achievement of particle

AN US

regularity. In this algorithm, the maximum absolute value of particle concentration gradient corresponding to typical target particles in the whole computational domain, i.e. max Ci , is considered as a criterion for iterating the particle shifting process. The particle shifting is

M

carried out until the value, max Ci , becomes A times smaller than a pre-defined threshold, Ct , where A = 0.1 [33]. Selection of an appropriate value of Ct leads to an adequate

PT

calculation domain.

ED

number of iterative shifting to guarantee a regularized particle distribution in the whole

In this appendix, the aforementioned iterative shifting scheme is extended for free-

CE

surface/multiphase flows and employed for improved MPS + OPS, resulting in improved MPS + iterative OPS. Comparisons would be carried out in between improved MPS + OPS

AC

and improved MPS + iterative OPS. In case of OPS, due to existence of free-surface or freesurface vicinity particle shifting, Ci , or the absolute value of concentration gradient of a typical target particle i, is newly defined as follows:

 CSi  Ci   Ci  0 

i   F   FV  i  I ; CSi  I  n~i  n~i Ci i   SP

(20)

ACCEPTED MANUSCRIPT where ΩF, ΩFV, ΩI and ΩSP signify Free-surface, Free-surface Vicinity, Inner and Splash Particles as described in Fig. 3. In Eq. (20), CSi symbolizes the concentration gradient calculated for either free-surface or free-surface vicinity. In this study, the constant Ct is set as 1.0, which leads to a relatively small threshold value to guarantee an almost regular particle

incorporated in Improved MPS + iterative OPS.

CR IP T

distribution (A Ct = 0.1). Fig. 26 presents a pseudocode for iterative shifting scheme

The effect of iterative shifting is investigated in simulation of three benchmark tests. Fig. 27 corresponds to multiphase perturbation test and depicts the time histories of x-direction randomness index  x [12] for the initially perturbed particle (Fig. 6) by multiphase improved

AN US

MPS incorporating either OPS or iterative OPS. As the figure shows implementation of iterative shifting has resulted in minimization of randomness index, which is a measure for particle disorder.

M

Figs. 28 and 29 correspond to the standing gravity wave simulations. Fig. 28 portrays snapshots of particles together with density field by improved MPS + OPS and improved MPS

ED

+ iterative OPS. The particle distributions for both OPS and iterative OPS snapshots are almost regular, while iterative shifting has apparently resulted in a slightly more regular particle

PT

distribution, e.g. close to the interface. Fig. 29 presents the time histories of water surface elevation at the center of the tank, where both improved MPS + OPS and improved MPS +

CE

iterative OPS seem to have resulted in almost similar results. This would suggest that OPS and

AC

iterative OPS seem to possess almost similar energy conservation properties, while iterative OPS tends to provide at least slightly more regular particle distributions. Fig. 30 corresponds to the multi-fluid square patch test. The figure portrays snapshots of

particles together with density field at t = 2.04 by improved MPS + OPS and improved MPS + iterative OPS. The effectiveness of iterative shifting is distinguishable from several enlarged portions of the reproduced fluid patch. For instance, the iterative OPS has been effective in reproducing a smoother free-surface. In addition, for inner fluid particles, particle distribution

ACCEPTED MANUSCRIPT by improved MPS implementing iterative OPS seems to be more regular and isotropic. Indeed, in our future works, we will further examine and validate the iterative OPS through consideration of various benchmark tests. Tables 4 and 5 present the required CPU time in simulation of multi-fluid square patch until t = 2.04 corresponding to improved MPS + OPS and improved MPS + iterative OPS,

CR IP T

respectively. In the presented tables, prediction and correction correspond to calculations conducted in prediction (viscous/gravitational acceleration calculations, corresponding neighbor search, etc.) and correction (solution of PPE, calculation of pressure gradient acceleration, corresponding neighbor search, etc.) subroutines except for all required

AN US

calculations corresponding to particle shifting as well as others (data storage and output). From the presented tables, application of iterative OPS in place of OPS has increased the overall required CPU time by about 66% for multi-fluid square patch simulation. The increase in total

ED

5. CONCLUDING REMARKS

M

CPU time for the standing gravity wave simulation is about 50% (Tables 6 and 7).

A novel projection-based [39] particle method is presented for simulation of multiphase

PT

flows characterized by high density ratios. The distinct feature of the proposed method corresponds to its capability of preserving material discontinuity at the interface along with

CE

volume/pressure continuity, thanks to a careful implementation of particle shifting [23] concept. The method considers a multi-fluid system within the framework of MPS (Moving

AC

Particle Semi-implicit) [2] and applies a two-step shifting process at the phase interface. At the first step, shifting of the heavy phase particles are conducted based on the OPS (Optimized Particle Shifting) scheme [1] with ignoring the light phase particles, such that a consistent shifting of heavy phase particles is achieved without any unphysical penetrations into the lighter phase. At the second step both heavy and light phase particles are considered and thus, shifting is simply made based on PS (Particle Shifting) [23] scheme. The proposed multiphase method benefits from a set of previously developed improved schemes to achieve more

ACCEPTED MANUSCRIPT accurate MPS-based simulations and comprises of a two-step shifting process distinctively on the basis of OPS, hence, it is simply referred to as improved MPS + OPS. Validations are made both qualitatively and quantitatively in terms of accuracy, energy conservation properties as well as convergence properties by consideration of four benchmark tests. Comparisons are also made by considering multiphase improved MPS methods that

CR IP T

apply PS [23] or FDS (First-order accurate Density Smoothing) [12] schemes. The improved MPS + PS possesses exactly the same computational algorithm and schemes as improved MPS + OPS, the only difference is related to application of PS [23] in place of OPS in heavy phase particle shifting. The improved MPS + FDS utilizes a First-order accurate Taylor series-based

AN US

Density Smoothing scheme for interface stabilization instead of particle shifting.

The results of this study highlight the effectiveness of particle shifting concept to achieve consistent and accurate results in multiphase particle-based simulations. At the same time, the results underline the importance of precise and consistent implementation of particle shifting to

M

achieve enhanced results in terms of accuracy, energy conservation and convergence. Particle shifting evidently outperformed density smoothing in all conducted multiphase simulations of

ED

this study.

PT

The iterative shifting scheme of Vacondio and Rogers [33] is extended for freesurface/multiphase flows and is employed for improved MPS + OPS, resulting in improved

CE

MPS + iterative OPS. The iterative OPS has been shown to be effective in providing further regularized particle distributions while maintaining the overall accuracy/conservation property.

AC

Future works include further validations as well as further development of the proposed

method for simulation of violent multi-phase flows through more precise implementation of particle shifting scheme as well as incorporation of accurate models for physical stabilizing terms such as interfacial surface tension and viscosity. In this regard, more challenging numerical tests characterized by high density ratios as well as low Bond numbers will be targeted. In addition, in our future works, detailed investigations on the performance of iterative shifting approach will be conducted by considering violent test cases. The proposed

ACCEPTED MANUSCRIPT method is also expected to be extended for more practical engineering/industrial applications, e.g. hydroelastic FSI (Fluid-Structure Interaction) [40] comprising of two-phase air-water flows.

Acknowledgements

CR IP T

The authors would like to express their gratitude to Dr. Steven Lind, at the University of Manchester for kindly providing them with a figure related to their ICSPH [22] considered for

AN US

comparison in Fig. 16 of this paper.

References

[1] Khayyer A., Gotoh H. and Shimizu Y., “Comparative study on accuracy and conservation properties of two particle regularization schemes and proposal of an optimized particle

M

shifting scheme in ISPH Context,” Journal of Computational Physics, 332, 236-256, 2017.

ED

[2] Koshizuka S. and Oka Y., “Moving particle semi-implicit method for fragmentation of

PT

incompressible fluid,” Nuclear Science and Engineering, 123, 421-434, 1996. [3] Ran Q., Tong J., Shao S., Fu X. and Xu Y., “Incompressible SPH scour model for

CE

movable bed dam break flows,” Advances in Water Resources, 82, 39–50, 2015. [4] Shao S., “Incompressible smoothed particle hydrodynamics simulation of multifluid

AC

flows,” International Journal for Numerical Methods in Fluids, 69(11), 1715–1735, 2012.

[5] Gong K., Shao S., Liu H., Wang B. and Tan S.-K., “Two-phase SPH simulation of fluid– structure interactions,” Journal of Fluids and Structures, 65, 155-179, 2016.

[6] Shimizu Y., Gotoh H. and Khayyer A., “An MPS-based particle method for simulation of multiphase flows characterized by high density ratios by incorporation of space potential particle concept,” Computers & Mathematics with Applications, 76(5), 1108-1129, 2018.

ACCEPTED MANUSCRIPT [7] Lind S.J., Stansby P.K. and Rogers B.D., “Fixed and moored bodies in steep and breaking waves using SPH with the Froude–Krylov approximation,” Journal of Ocean Engineering and Marine Energy, 2(3), 331-354, 2016. [8] Xenakis A.M., Lind S.J., Stansby P.K. and Rogers B.D., “Landslides and tsunamis predicted by incompressible smoothed particle hydrodynamics (SPH) with application to

CR IP T

the 1958 Lituya Bay event and idealized experiment,” Proceedings of The Royal Society A Mathematical Physical and Engineering Sciences, 473(2199), 20160674, 2017.

[9] Gotoh H. and Khayyer A., “On the state-of-the-art of particle methods for coastal and ocean engineering,” Coastal Engineering Journal, 60(1), 79-103, 2018.

Khayyer A., Gotoh H., Shimizu Y., Gotoh K., Falahaty H. and Shao S.,

AN US

[10]

“Development of a projection-based SPH method for numerical wave flume with porous media of variable porosity,” Coastal Engineering, 140, 1-22, 2018. [11]

Monaghan J.J., “Simulating free surface flows with SPH,” Journal of Computational

[12]

M

Physics, 110, 399-406, 1994.

Khayyer A. and Gotoh H., “Enhancement of performance and stability of MPS

ED

meshfree particle method for multiphase flows characterized by high density ratios,”

[13]

PT

Journal of Computational Physics, 242, 211-233, 2013. Shakibaeinia A. and Jin Y.C., “MPS mesh-free particle method for multi-phase

CE

flows,” Computer Methods in Applied Mechanics and Engineering, 229-232, 13-26, 2012.

Rezavand M., Taeibi-Rahni M. and Rauch W., “An ISPH scheme for numerical

AC

[14]

simulation of multiphase flows with complex interfaces and high density ratios,” Computers and Mathematics with Applications, 75(8), 2658-2677, 2018..

[15]

Duan G., Chen B., Koshizuka S. and Xiang H., “Stable multiphase moving particle

semi-implicit method for incompressible interfacial flow,” Computer Methods in Applied Mechanics and Engineering, 318, 636-666, 2017.

ACCEPTED MANUSCRIPT [16]

Luo M., Koh C.G., Gao M. and Bai W., “A particle method for two-phase flows with

large density difference,” International Journal for Numerical Methods in Engineering, 103(4), 235-255, 2015. [17]

Mokos A., Rogers B.D. and Stansby P.K., “A multi-phase particle shifting algorithm

for SPH simulations of violent hydrodynamics with a large number of particles,” Journal

[18]

CR IP T

of Hydraulic Reseearch, 55(2), 143-162, 2017. Cao X.Y., Ming F.R., Zhang A.M. and Tao L., “Multi-phase SPH modelling of air

effect on the dynamic flooding of a damaged cabin,” Computers and Fluids, 163, 7-19, 2018.

Krimi A., Rezoug M., Khelladi S., Nogueira X., Deligant M. and Ramírez L.,

AN US

[19]

“Smoothed Particle Hydrodynamics: A consistent model for interfacial multiphase fluid flow simulations,” Journal of Computational Physics 358, 53-87, 2018. [20]

Monaghan J.J. and Rafiee A., “A simple SPH algorithm for multi-fluid flow with

M

high density ratios,” International Journal for Numerical Methods in Fluids, 71, 537-561, 2013.

Khayyer A. and Gotoh H., “A multiphase compressible-incompressible particle

ED

[21]

PT

method for water slamming,” International Journal of Offshore and Polar Engineering, 26(1), 20-25, 2016.

Lind S.J., Stansby P.K. and Rogers B.D., “Incompressible–compressible flows with a

CE

[22]

transient discontinuous interface using smoothed particle hydrodynamics (SPH),” Journal

AC

of Computational Physics, 309, 129-147, 2016.

[23]

Lind S.J., Xu R., Stansby P.K. and Rogers B.D., “Incompressible smoothed particle

hydrodynamics for free-surface flows: a generalised diffusion-based algorithm for stability and validations for impulsive flows and propagating waves,” Journal of Computational Physics, 231(4), 1499-1523, 2012.

ACCEPTED MANUSCRIPT [24]

Xu R., Stansby P.K. and Laurence D., “Accuracy and stability in incompressible SPH

(ISPH) based on the projection method and a new approach,” Journal of Computational Physics, 228, 6703–6725, 2009. [25]

Leroy A., Violeau D., Ferrand M. and Joly A., “Buoyancy modelling with

incompressible SPH for laminar and turbulent flows,” International Journal for Numerical

[26]

CR IP T

Methods in Fluids, 78(8), 455-474, 2015. Sun P.N., Colagrossi A., Marrone S. and Zhang A.M., “The δplus-SPH model:

Simple procedures for a further improvement of the SPH scheme,” 315, 25-49, 2017. [27]

Fourtakas G., Stansby P.K., Rogers B.D. and Lind S.J., “An Eulerian–Lagrangian

AN US

incompressible SPH formulation (ELI-SPH) connected with a sharp interface,” Computer Methods in Applied Mechanics and Engineering, 329, 532-552, 2018. [28]

Fourtakas G. and Rogers B.D., “Modelling multi-phase liquid-sediment scour and

resuspension induced by rapid flows using Smoothed Particle Hydrodynamics (SPH)

M

accelerated with a Graphics Processing Unit (GPU),” Advances in Water Resources, 92, 186-199, 2016.

Lind S.J., Stansby P.K., Rogers B.D. and Lloyd P.M., “Numerical predictions of wave

slam

using

incompressible–compressible

smoothed

particle

PT

water–air

ED

[29]

hydrodynamics,” Applied Ocean Research, 49, 57-71, 2015. Lind S.J., Fang Q., Stansby P.K., Rogers B.D. and Fourtakas G., “A Two-Phase

CE

[30]

Incompressible-Compressible (Water-Air) Smoothed Particle Hydrodynamics (ICSPH)

AC

Method Applied to Focused Wave Slam on Decks,” The 27th International Ocean and Polar Engineering Conference, California, USA, June, 2017.

[31]

Khayyer A., Gotoh H., Shimizu Y. and Gotoh K., “On enhancement of energy

conservation properties of projection-based particle methods,” European Journal of Mechanics - B/Fluids, 66, 20-37, 2017.

ACCEPTED MANUSCRIPT [32]

Gotoh H. and Khayyer A., “Current achievements and future perspectives for

projection-based particle methods with applications in ocean engineering,” Journal of Ocean Engineering and Marine Energy, 2(3), 251-278, 2016. [33]

Vacondio R. and Rogers B.D., “Consistent Iterative shifting for SPH methods”,

Proceedings of 12th international SPHERIC workshop, Ourense, Spain, June, 2017. Lee E.S., Moulinec C., Xu R., Violeau D., Laurence D. and Stansby P.K.,

CR IP T

[34]

“Comparisons of weakly compressible and truly incompressible algorithms for the SPH mesh free particle method, “ Journal of Computational Physics, 227, 8417–8436, 2008. [35]

Colagrossi A., “A meshless Lagrangian Method for Free-Surface and Interface Flows

[36]

AN US

with Fragmentation,” PhD Thesis, Universita di Roma, La Sapienza, 2003.

Skillen A., Lind S.J., Stansby P.K. and Rogers B.D., “Incompressible smoothed

particle hydrodynamics (SPH) with reduced temporal noise and generalised fickian smoothing applied to body–water slam and efficient wave–body interaction”, Comput.

[37]

M

Methods Appl. Mech. Eng. 265, 163-173, 2013.

Antuono A., Bouscasse B., Colagrossi A. and Marrone S., “A measure of spatial

Wu G.X., and Taylor E.R., “Finite element analysis of two-dimensional nonlinear

PT

[38]

ED

disorder in particle methods,” Computer Physics Communications, 185, 2609-2621, 2014.

transient water waves,” Applied Ocean Research, 16(6), 363–372, 1994. Chorin A.J., “Numerical solution of the Navier-Stokes equations,” Mathematics of

CE

[39]

Computation, 22, 745-762, 1968. Khayyer A., Gotoh H., Falahaty H. and Shimizu Y., “An enhanced ISPH–SPH

AC

[40]

coupled method for simulation of incompressible fluid–elastic structure interactions,” Computer Physics Communications, 232, 139-164, 2018.

M

AN US

CR IP T

ACCEPTED MANUSCRIPT

AC

CE

PT

ED

Figure 1. Computational flowchart of multiphase improved MPS + OPS method

CR IP T

ACCEPTED MANUSCRIPT

AC

CE

PT

ED

M

AN US

Figure 2. Two-step shifting process in multiphase improved MPS + OPS method

CE

PT

ED

M

AN US

CR IP T

ACCEPTED MANUSCRIPT

AC

Figure 3. Details of particle shifting process in multiphase improved MPS + OPS method

CR IP T

ACCEPTED MANUSCRIPT

AN US

Figure 4. Effectiveness of PS and OPS schemes for a single-phase simulation by improved MPS – snapshots of particles together with normalized pressure field at t =

AC

CE

PT

ED

M

2.0 in a square patch of fluid test (Colagrossi [35])

CR IP T

ACCEPTED MANUSCRIPT

Figure 5. Time histories of normalized pressure at the center of the patch – a (single-

AC

CE

PT

ED

M

AN US

phase) square patch of fluid test (Colagrossi [35])

ACCEPTED MANUSCRIPT Figure 6. Schematic illustration of a simple multiphase perturbation test - a particle at

CE

PT

ED

M

AN US

CR IP T

the phase interface is perturbed 0.1d0 horizontally as shown in the figure

AC

Figure 7. Snapshots of particles together with density field at t = 30.0 s by multiphase improved MPS incorporating either PS or OPS – a simple multiphase perturbation test

M

AN US

CR IP T

ACCEPTED MANUSCRIPT

Figure 8. Pressure profile at x = 1.0m / t = 30.0 s by multiphase improved MPS + OPS

AC

CE

PT

ED

– a simple multiphase perturbation test

AN US

CR IP T

ACCEPTED MANUSCRIPT

Figure 9. Time variations of x-direction randomness index  x [12] for the initially

AC

CE

PT

ED

M

perturbed particle – a simple multiphase perturbation test

AN US

CR IP T

ACCEPTED MANUSCRIPT

AC

CE

PT

ED

M

Figure 10. Schematic illustration of oscillating concentric elliptical regions test

M

AN US

CR IP T

ACCEPTED MANUSCRIPT

Figure 11. Snapshots of particles together with density field at t = 6.0 s by multiphase

AC

CE

PT

ED

improved MPS implementing either OPS, PS or FDS – oscillating concentric elliptical regions test

CE

PT

ED

M

AN US

CR IP T

ACCEPTED MANUSCRIPT

AC

Figure 12. Time histories of the semi-major axis of the outer ellipse (A) and semiminor axis of the inner ellipse (b) – oscillating concentric elliptical regions

AC

CE

PT

ED

M

AN US

CR IP T

ACCEPTED MANUSCRIPT

Figure 13. Time histories of kinetic and potential energies – oscillating concentric elliptical regions

CR IP T

ACCEPTED MANUSCRIPT

AN US

Figure 14. Time histories of normalized total energies – oscillating concentric elliptical

AC

CE

PT

ED

M

regions

AN US

CR IP T

ACCEPTED MANUSCRIPT

M

Figure 15. Time histories of the semi-major axis of the outer ellipse (A) by multiphase improved MPS + OPS for simulations with a set of three different particle diameters,

AC

CE

PT

ED

d0 = 0.008 m, 0.010 m and 0.020 m – oscillating concentric elliptical regions test

ED

M

AN US

CR IP T

ACCEPTED MANUSCRIPT

PT

Figure 16. Snapshots of particles together with pressure and density fields by multiphase improved MPS + OPS and Incompressible-Compressible SPH (ICSPH)

AC

CE

results by Lind et al. [22] at t = 2.2 s – oscillating concentric elliptical regions test

AN US

CR IP T

ACCEPTED MANUSCRIPT

AC

CE

PT

ED

M

Figure 17. Schematic illustration of multi-fluid square patch test

AC

CE

PT

ED

M

AN US

CR IP T

ACCEPTED MANUSCRIPT

Figure 18. Snapshots of particles and density field in a multi-fluid square patch test at t = 0.44 – results by multiphase improved MPS incorporating the OPS scheme (improved MPS + OPS)

PT

ED

M

AN US

CR IP T

ACCEPTED MANUSCRIPT

Figure 19. Snapshots of particles as well as density and pressure fields at t = 2.04 by

AC

test

CE

multiphase improved MPS implementing either PS or OPS – multi-fluid square patch

CR IP T

ACCEPTED MANUSCRIPT

Figure 20. Time histories of normalized kinetic energy by multiphase improved MPS

AC

CE

PT

ED

M

AN US

implementing either OPS, PS or FDS – multi-fluid square patch test

Figure 21. Schematic illustration of standing gravity wave test

CE

PT

ED

M

AN US

CR IP T

ACCEPTED MANUSCRIPT

Figure 22. Snapshots of particles together with density field at t = 2.0 s – standing

AC

gravity wave

AN US

CR IP T

ACCEPTED MANUSCRIPT

Figure 23. Time histories of water surface elevation at the tank center by multiphase

AC

CE

PT

ED

M

improved MPS incorporating either OPS, PS or FDS – standing gravity wave

AN US

CR IP T

ACCEPTED MANUSCRIPT

Figure 24. Time histories of water surface elevation at the tank center by multiphase improved MPS incorporating the OPS with three different sets of particle diameter of

AC

CE

PT

ED

M

d0 = 0.004 m, 0.005 m and 0.010 m – standing gravity wave

CR IP T

ACCEPTED MANUSCRIPT

AN US

Figure 25. Quantitative profile of particle diameter versus NRMSE corresponding to the data in Table 3 (Fitted curve: log y = 0.9609 log x + 0.9173; x and y correspond to

AC

CE

PT

ED

M

particle diameter and NRMSE, respectively)

AC

CE

PT

ED

M

AN US

CR IP T

ACCEPTED MANUSCRIPT

Figure 26. Pseudocode for iterative shifting scheme incorporated in Improved MPS + iterative OPS

CR IP T

ACCEPTED MANUSCRIPT

AN US

Figure 27. Time histories of x-direction randomness index  x for the initially perturbed particle (Fig. 6) by multiphase improved MPS incorporating either OPS or iterative

AC

CE

PT

ED

M

OPS – a simple multiphase perturbation test

ACCEPTED MANUSCRIPT Figure 28. Snapshots of particles together with density field at t = 2.0 s by multiphase

ED

M

AN US

CR IP T

improved MPS incorporating either OPS or iterative OPS – standing gravity wave

PT

Figure 29. Time histories of water surface elevation at the center of the tank by multiphase improved MPS incorporating either OPS or iterative OPS – standing

AC

CE

gravity wave

M

AN US

CR IP T

ACCEPTED MANUSCRIPT

Figure 30. Snapshots of particles together with density field at t = 2.04 by multiphase

ED

improved MPS implementing either OPS or iterative OPS in a multi-fluid square patch

AC

CE

PT

simulation

ACCEPTED MANUSCRIPT

Table 1. Root Mean Square Error (RMSE) and Normalized RMSE (NRMSE) corresponding to numerical data in Fig. 12 - oscillating concentric elliptical regions test PS

FDS

RMSE (A)

0.01268

0.03593

0.018013

NRMSE (A)

0.01799

0.05097

0.025553

RMSE (b)

0.00383

0.01830

0.00691

NRMSE (b)

0.01095

0.05230

0.01976

CR IP T

OPS

AN US

Table 2. Root Mean Square Error (RMSE) and Normalized RMSE (NRMSE) corresponding to numerical data in Fig. 15 - oscillating concentric elliptical regions test d0 = 8.0E-3 0.00992

NRMSE

0.01407

d0 = 2.0E-2

0.01268

0.03541

0.01799

0.05023

ED

M

RMSE

d0 = 1.0E-2

Table 3. Root Mean Square Error (RMSE) and Normalized RMSE (NRMSE) corresponding

PT

to numerical data in Fig. 24 - standing gravity wave d0 = 5.0E-3

d0 = 1.0E-2

RMSE

0.004909

0.005603

0.011559

NRMSE

0.042516

0.048526

0.100113

AC

CE

d0 = 4.0E-3

Table 4. Required CPU time for improved MPS + OPS in multi-fluid square patch simulation until t = 2.04

ACCEPTED MANUSCRIPT

Correction

Others

OPS

Overall

1 step (s)

0.086511

0.5845632

0.07722

0.226417185

0.974716

Total (s)

882.4991

5963.129

787.766

2309.6817

9943.0783

8.9%

60.0%

7.9%

23.2%

100%

Percentage (%)

CR IP T

Prediction

Table 5. Required CPU time for improved MPS + iterative OPS in multi-fluid square patch simulation until t = 2.04 Correction

1 step (s)

0.0822863

0.5684876

0.07809

0.887001216

1.6158661

Total (s)

839.4028

5799.1424

796.598

9048.2994

16483.45

5.1%

35.2%

4.8%

54.9%

100%

Percentage

Iterative OPS

Overall

AC

CE

PT

ED

M

(%)

Others

AN US

Prediction

Table 6. Required CPU time for improved MPS + OPS in standing gravity wave simulation until t= 1 s

ACCEPTED MANUSCRIPT

Prediction

Correction

Others

OPS

Overall

1 step (s)

0.2075293

0.8206522

0.0784

0.462077316

1.5686578

Total (s)

2117.0061

8371.4734

799.744

4713.6507

16001.878

13.2%

52.3%

5.0%

29.5%

100%

(%)

CR IP T

Percentage

Table 7. Required CPU time for improved MPS + iterative OPS in standing gravity wave simulation until t= 1 s Correction

1 step (s)

0.207639

0.8085539

0.08037

1.265718371

2.3622806

Total (s)

2118.1256

8248.0585

819.844

12911.5931

24097.624

8.8%

34.2%

3.4%

53.6%

100%

Percentage

AC

CE

PT

ED

M

(%)

Others

Iterative OPS

AN US

Prediction

Overall