Fengze Xie,Weiwen Zhao and Decheng Wan

Abstract With the significant development of computer hardware, many advanced numerical techniques have been proposed to investigate complex hydrodynamic problems.This article aims to provide a detailed review of moving particle semi-implicit(MPS)techniques and their application in ocean and coastal engineering.The achievements of the MPS method in stability and accuracy,boundary conditions,and acceleration techniques are discussed.The applications of the MPS method,which are classified into two main categories, namely, multiphase flows and fluid-structure interactions, are introduced. Finally,the prospects and conclusions are highlighted.The MPS method has the potential to solve practical problems.

Keywords MPS technique; Ocean engineering; Coastal engineering; Stability;Accuracy; Boundary conditions;Acceleration techniques;Multiphase flows;Fluid-structure interactions

1 Introduction

The ocean environment is closely tied to human life and health. Products are transported from one region to another by sea, and the cost of transport by sea is the lowest,which promotes economic globalization. With the further exploitation of marine resources, industrialization has been further accelerated. Although human beings benefit from marine resources,they also suffer from marine disasters. Tsunamis induced by earthquakes destroy inshore facilities and threaten the safety of people (Hayashi and Hughes 2013). To better adapt to the ocean environment and fully utilize marine resources,the hydrodynamic problems in the field of ocean and coastal engineering were investigated by researchers. In the past decades, computational fluid dynamics techniques have advanced by leaps and bounds. With the assistance of numerical tools, more comprehensive information on the flow field can be obtained at a lower cost.

Particle-based methods exhibit their superiority in three main aspects. First, there are not topological relationships among Lagrangian particles. Therefore, particle-based methods can better capture the free surface and interface with large deformation.Second,in the simulations of moving boundaries with mesh-based methods,errors will be introduced because of the deformation of the mesh.Remeshing is needed to avoid ill-shaped meshes, which is timeconsuming (Sun et al. 2019). Particle-based methods can simulate moving boundaries without additional processing. Third, based on the Lagrangian view, the momentum equation of the particle-based method has no convection term.Accordingly, the numerical diffusion induced by the calculation of the advection term is avoided. Smoothed particle hydrodynamics (SPH) and moving particle semiimplicit (MPS) techniques are two mainstream particlebased methods. SPH was first proposed to solve astrophysical problems(Gingold and Monaghan 1977;Lucy 1977)and was subsequently applied to the hydrodynamic(Monaghan 1994) and structural dynamics (Antoci 2007) problems.Based on the hypothesis that the fluid is weakly compressible, the pressure distribution of the flow field can be obtained by explicitly solving the equation of state.The SPH can be easily parallelized with less efficiency loss (Ferrari et al. 2009). However, because of the pressure fluctuation,the original SPH cannot be further extended to complex problems. MPS was first developed to solve the hydrodynamic problems of nuclear engineering (Koshizuka and Oka 1996).The pressure is calculated by solving the pressure Poisson equation (PPE), which consumes the most time during the MPS simulation. The parallelization of MPS, particularly of PPE, is more difficult than that of SPH,and its parallelization efficiency is not as high as that of SPH. Therefore, further enlarging the simulation scale of MPS is difficult. Because of the prediction-correction semi-implicit scheme, the MPS performs better than the original SPH in terms of stability (Luo et al. 2021). In the past years,many numerical techniques have been proposed and introduced to SPH and MPS, and their stability and computational efficiency have been improved dramatically.

Many excellent review articles about particle-based methods have been published. Luo et al. (2021) comprehensively introduced the advances of particle-based methods and their applications in ocean and coastal engineering.The review published by Li et al. (2020) focused on the MPS method and its application in nuclear engineering.Gotch et al.(2021)reviewed the particle-based methods for fluid-structure interactions (FSI) in ocean engineering in three aspects, namely, reliability, adaptivity, and generality.This paper aims to provide a review of the latest advances in the MPS method and their application in the field of ocean and coastal engineering. First, the achievements of the MPS method in stability and accuracy,boundary conditions, and acceleration techniques are briefly summarized.Then, the applications of the MPS method in the field of ocean and coastal engineering are presented.Finally,future works for further improvement and potential application of MPS are discussed.A brief conclusion is drawn as well.

2 Review of the MPS method

2.1 Original MPS method

The governing equations for fluid include the continuity and Navier-Stokes (N-S) equations, which can be written as follows:

whereU,ρ,P,andνare the velocity,density,pressure,and kinematic viscosity of the fluid,respectively,tis the physical time,andgis the gravity vector.

The particle interaction models,which contain the gradient, divergence, and Laplacian models, are used to discretize the governing equations (Koshizuka and Oka 1996),which can be expressed as follows:

whereiandjare the subscripts of fluid particles,dis the dimensions of the space,ϕis the physical value of the MPS particle,n0is initial particle density,ris the position vector relative to origin andλis a parameter used to compensate for the errors resulting from the limited range of the kernel function(Koshizuka and Oka 1996).

The pressure information of MPS particles is obtained by solving the PPE as follows:

wherePk+1is the pressure of the stepk+ 1 and Δtis the time step.

Based on the difference in the particle density between free-surface and internal particles, the original criteria of free-surface detection can be expressed as follows:

2.2 Stability and accuracy

2.2.1 Modification of gradient model

The original gradient model has several drawbacks.When the interparticle interaction strength increases as the particles approach, the tensile instability may occur because of the presence of attractive interparticle forces(Khayyer and Gotoh 2011). To improve the stability of simulations, Koshizuka et al. (1998) proposed an improved gradient model(Eq.(8)),in which the pressure of particleiwas altered by the minimum pressure among the neighboring particles so that the interparticle forces were repulsion.A gradient model (Eq. (9)) with exact conservation of linear momentum and significantly enhanced preservation of angular momentum was proposed by Khayyer and Gotoh(2008). The interparticle forces of two neighboring particles are equal and opposite.Another conservative formulation of the gradient model(Eq.(10))was proposed by Tanaka and Masunaga (2010). If the MPS particles distribute uniformly, then the aforementioned gradient models are equivalent to the original gradient model.If the MPS particles distribute nonuniformly, then the artificial repulsive force drives the MPS particles from high to low concentrations; its basic concept is similar to that of the particle shifting(PS)technique.

According to the analysis by Tsuruta (2013), the artificial repulsive force rather than the original gradient force was predominant in Eqs. (8) and (9), which might lead to the unphysical motion and perturbations of fluid particles.To solve this problem, a meticulously adequate stabilizing force was introduced to the original gradient model,as expressed in Eq.(11):

The error term of the original gradient model is zero-order convergence (Duan et al. 2019). The discretization error can be reduced,and the accuracy of the gradient model can be improved after high-order schemes are adopted.Khayyer and Gotoh (2011) proposed a corrective matrix based on the Taylor series expansion for a more accurate approximation of pressure gradient,as expressed in the following equations:

Wang et al. (2017) introduced the corrective matrix to Eq.(8)and obtained a stabilized gradient model as follows:

If the particle distribution is anisotropic, then the zeroorder discretization error will be introduced. Hence, dummy particles were introduced to the neighboring field to compensate for the zero-order discretization error of the gradient model (Liu et al. 2018), as shown in Figure 1.The more accurate and stable gradient model can be written as follows

where the subscriptdrepresents the dummy particles.

Figure 1 Estimated position of a dummy neighboring particle for particle i(Liu et al.2018)

Duan et al. (2019) developed a second-order corrective matrix (SCM) to further reduce the truncation error of the gradient model. Compared with the first-order corrective matrix,the SCM significantly reduces the discretization error.

2.2.2 Modification of the PPE

The pressure distribution of the flow field is obtained by solving the PPE. Therefore, the construction and discretization of the equation are relative to the accuracy of the pressure calculation.The PPE consists of the Laplacian operator(left)and source term(right).

The original Laplacian model exhibits minus first order convergence (Duan et al. 2019). Khayyer and Gotch(2010; 2012) developed high-order Laplacian models(Eqs. (21) and (22)) based on the original MPS kernel function for 2D and 3D simulations.Duan et al.(2018)introduced a corrected vector to the Laplacian model to compensate for anisotropy error. The corrected model was proven to exhibit convergence and reduce unphysical fluctuation.Liu et al.(2018)used the dummy neighboring particle for the discretization of the Laplacian model (Eq. (25));its basic principle was similar to that of the corrected vector (Duan et al. 2018). Duan et al. (2019) introduced the SCM to the Laplacian model (Eq. (26)), which exhibited first-order convergence.

The construction of the PPE source term affects the accuracy and smoothness of the numerical flow field. In the original MPS, incompressibility is influenced by the particle number density (PND). However, the PND at each time step would not be exactly equal to the initial particle density,and the numerical errors would accumulate,which may result in considerable pressure fluctuation. Khayyer and Gotoh (2009) developed a high-order source term by calculating the time differentiation of PND based on the original kernel function. To balance the accuracy and continuity, many multiterm formulations for PPE are proposed.Kondo and Koshizuka(2011)derived a more generic formulation,including one main part and two error compensating parts,which could prevent density error accumulation.The coefficientsαandβare obtained by hydrostatic pressure calculation, which would not be necessarily appropriate for other simulation cases, particularly for violent flows.Khayyer and Gotoh(2011)proposed the coefficients based on the features of instantaneous flow and applied the high-order schemes for the main part and error compensating parts. The incompressibility of the fluid can also be satisfied by the divergence-free condition. Tanaka and Masunaga (2010) proposed a novel PPE, which had both the advantages of the PND condition that the error was not stored and the advantages of the divergence-free condition that smoother pressure was obtained. The divergence-free term varies smoothly with time, which is similar to the first and second terms in Eqs.(28)and(29).

2.2.3 PS schemes

The PS scheme can rearrange seriously distorted particle distribution and consider the transporting effects of shifts, which can also significantly reduce stability errors(Duan et al. 2019). Xu et al. (2009) first proposed the PS scheme based on Fick’s law of diffusion, which moved particles from high to low concentrations. The PS scheme was integrated into the incompressible SPH(ISPH)framework and applied to the simulation of internal flows. Because of the incomplete kernel support, the free-surface particles’ distribution may be distorted by PS. Lind et al.(2012) suggested reducing or eliminating the normal shifting of free-surface particles and applied the PS to the freesurface flows.However,this treatment may lead to numerical inconsistencies, including the unphysical diffusion of free-surface particles,increased level of perturbations,and/or an unphysical gap between free-surface and nearby particles (Khayyer et al. 2017).To overcome those drawbacks,Khayyer et al. (2017) proposed the optimized PS (OPS)and classified the shifting particles as splashing, freesurface, vicinity, and inner particles. The normal shifting of free-surface and vicinity particles is excluded. Compared with the PS (Lind et al. 2012), the OPS does not include any adjusting parameter and obtains more accurate numerical results. Duan et al. (2018) applied the OPS and conservative gradient model to produce the tangent shifting and normal adjustments of free-surface particles, respectively.

2.3 Boundary conditions

2.3.1 Free surface

In the MPS and other projection-based methods, the pressure exerted by the detected free-surface particles is set to a constant value, which is the Dirichlet boundary condition of PPE. However, if the inner particles are misjudged as the free-surface particles, then their calculated pressure and that of their neighbors would be incorrect,which may lead to pressure fluctuation. Therefore, a highprecision free-surface detection model can enhance the stability and accuracy of the simulation.

The free-surface detection method can be classified into four types of approaches (Luo et al. 2021). The first approach is focused on the concentration of the particles. In the original MPS (Koshizuka and Oka 1996), the PND was used to detect the free-surface particles. Lee et al.(2008)utilized the divergence of particle positions to track the free-surface particles.Although the first approach is efficient, some free-surface particles may be undetected,

whereas some inner particles may be incorrectly identified as free-surface particles(Khayyer and Gotch 2008;Khayyer et al.2008).The nonsymmetric distribution of the neighboring particles of free-surface particles is the basic concept of the second approach. Khayyer et al. (2009) used the summation of thexorycoordinates of neighbors to detect the free-surface particles. Zhang et al. (2014) defined a vector function to detect the free-surface particles.For the third approach,the free-surface particles are detected based on the geometric features of their relative positions.Marrone et al.(2010)defined an umbrella-shaped area and used the conditions expressed in Eq.(31)to evaluate whether the target particleibelongs to the free surface or not(as shown in Figure 2):

whereniandτiare the local unit normal and tangential vectors of the free surface, respectively. Surface detection with the third approach is accurate but time-consuming.Therefore, many scholars tried to utilize the advantages of the aforementioned approaches and developed several hybrid approaches, i.e., the fourth category of approaches.Wang et al. (2019) coupled the first approach (Lee et al.2008) and the third approach (Marrone et al. 2010) to track the free surface accurately with low computational cost. Several scholars tried to enforce the Dirichlet boundary condition without free-surface detection and avoid misjudging the free-surface particle. Sun et al. (2014) proposed a novel discretization of the PPE and a new framework of the MPS method (NSD-MPS) using conceptual particles. The pressure fluctuation can be suppressed and the cluster of free-surface particles can be minimized.

Figure 2 Sketch of the umbrella-shaped region for free-surface detection(Marrone et al.2010)

2.3.2 Wall boundary condition

The use of dummy particles to discretize solid walls was the most common approach for boundary conditions,which was first proposed by Koshizuka and Oka (1996).This approach can prevent the penetration of fluid particles and provide sufficient kernel support for them. Lee et al. (2011) applied nonslip boundary conditions for simulation by introducing the moving dummy particles.Duan et al.(2020) imposed the Neumann boundary condition on the dummy particles and obtained satisfactory numerical results in the boundary value problem. However, simulating thin-walled structures using the original dummy particles is difficult.If the particles are large,then the fluid particles on both sides of thin-walled structures may unphysically interact with each other, which may affect the accuracy of the simulation. By contrast, the simulation with small particles is time-consuming. To overcome this drawback, He et al. (2018; 2019) proposed a multilayer dummy particle method. The dummy particles in each layer provide support for the fluid particles on their corresponding side, as shown in Figure 3.Even if this boundary condition is satisfactory for the simplified structures, it is difficult for dummy particles to adhere to structures with complex geometries with uniform particle spacing.

Figure 3 Sketch of a thin-walled structure using the multilayer DBP method(He et al.2019)

The mirror particle method is another widely used boundary condition.A corresponding mirror particle will be generated for the fluid particle near the solid boundaries,which can prevent the fluid particle from penetrating the wall effectively. Colagrossi and Landrini (2003) proposed the mirror particle method, with its physical values deduced from their corresponding fluid particles near the solid boundaries. Park and Jeun (2011) arranged one-line wall particles before the simulation, and the mirror particles reflected from the fluid particles generated at every time step. The mirror particles were only generated when there were fluid particles near the wall.Although the memory was reduced, the computational cost was significantly increased.Akimoto (2013) pointed out that the fluid particles near the corner had multiple imaginary particles associated with each boundary plane, which resulted in the inaccuracy of PND near the corner.A simple split line was used to divide the corner region to prevent the excessive and deficient generation of mirror particles. Because of the complex topological relationships among the fluid particles,solid walls, and mirror particles, extending this approach to complex geometries is difficult.

The repulsive function method can be used to replace the dummy and mirror particles. In this method, the solid boundary is modeled with one-line particles and provides the repulsive force to the fluid particles.The magnitude of the repulsive force is determined by the distance between the fluid particles and the solid boundary. Monaghan(1994) first applied the artificial repulsive forces in the form of Lennard-Jones to the fluid particles near the solid wall. The influence of boundary conditions on the tangential motion of fluid particles is ignored (Violeau and Rogers 2016). Monaghan and Kajtar (2009) used radial forces to obtain suitable boundary forces and the Poisson summation to derive the repulsive function. However, the incomplete kernel support for the fluid particle was not considered in the aforementioned methods.Sun et al.(2021)proposed a generic smoothed wall (GSW) boundary and applied the contact force to repel the fluid particles away from the solid boundary,which is based on the discrete element method (DEM) for solid-solid contact. Moreover, a boundary weight function was applied in the GSW boundary model to compensate for the PND of fluid particles near the wall, which has the same effect on the dummy particles.

Wall boundaries can also be represented by the polygon mesh, which is suitable for solid wall boundaries with complex geometries. Harada et al. (2008) first proposed the polygon wall (PW) boundary conditions without using the wall particles. The boundary contributes to the fluid particles through the distance function. However, severe pressure fluctuation was induced.Zhang et al.(2015)modified the weight function for nonplanar boundaries, which significantly suppressed the pressure fluctuation. Zhang et al. (2016) proposed the boundary particle arrangement technique to make the PND near the nonplanar wall boundaries more accurate. Then, Zhang et al. (2017) proposed a new PPE source term for the PW and obtained a stable pressure distribution. To further improve the accuracy of the simulation, Zhang et al. (2019b) determined the pressure of PW through the Neumann boundary condition and introduced a high-order gradient model to the PW.

2.3.3 Open boundary condition

The open boundary condition (OBC) is widely used in the simulation of ocean engineering problems. Different from the Euler methods,the flow field is discretized as Lagrangian particles. Therefore, how to deal with the fluid particles near the inlet and outlet is a key issue that needs to be considered.The most commonly used OBC was proposed by Lastiwka et al. (2009). The inflow and outflow regions are defined to provide complete kernel support for the particles in the fluid domain. One particle will be generated in the inflow region once one particle is transferred from the inflow region to the fluid domain. The particles will be eliminated when they move out of the flow region.Shakibaeinia and Jin(2009)proposed a similar OBC associated with the particle recycling strategy.Another type of particle called storage particles was defined. The generation and elimination of fluid particles by Lastiwka et al.(2009) corresponds to the addition and subtraction from storage particles by Shakibaeinia and Jin (2009). Wave generation can also be implemented using the OBC. Ni et al. (2018) developed a numerical wave flume (NWF)based on the OBC.The waves were generated by changing the number of fluid particles in the vertical direction in the buffer region. This NWF was able to generate multiple types of waves,including solitary,linear,and second-order regular waves. To further improve the efficiency of NWF,Ni et al. (2020) incorporated the SWE-SPH and NS-SPH.The information was exchanged through the OBC.Tsuruta et al. (2021) proposed a wavy interface model for MPS,which is similar to the NWF.The prescribed wave surface elevation in the buffer region was obtained through the generation and removal of fluid particles, as shown in Figure 4.

2.4 Acceleration technique

In the field of ocean and coastal engineering, the 3D effect cannot be ignored.Therefore,2D MPS needs to be extended to 3D MPS.However,the computational cost would increase dramatically. Many techniques have been proposed to improve the efficiency of MPS.Those techniques can be classified as multi-resolution,parallel, and graphics processing unit(GPU)acceleration techniques.

2.4.1 Multiresolution technique

The basic concept of the multiresolution technique is that the particles with small spacing are only arranged at the region that needs to be considered, which reduces the computational load and meets the requirements of local accuracy.Introducing the multiresolution technique to the incompressible MPS is difficult because the PPE is solved implicitly. The overlapping particle technique (OPT) was first proposed for MPS by Shibata et al.(2012)and extended to three dimensions by Tang et al.(2016).In the refined domain, the coarse particles can overlap with each other.The PPEs for coarse and fine particles were solved separately, which has only a slight influence on the stability of simulations. However, the interaction between MPS particles of different sizes is a one-way process, and only the coarse particles can affect the fine particles but not the other way around.Shibata et al.(2017)further developed a bilateral OPT and introduced several newly developed techniques to reduce the error in the total fluid mass. Tanaka et al. (2009) proposed a multiresolution technique to generate MPS particles of different sizes before the simulation.Tang et al. (2016) extended the multiresolution technique(Tanaka et al. 2009) to three dimensions and applied the technique to the simulations of the dam break and water entry. Tanaka et al. (2018) further introduced splitting or merging algorithms to the multiresolution technique to change the spatial resolution dynamically. Liu and Zhang(2021)first introduced the multiresolution technique to the least square MPS for the simulation of multiphase flows.The particle resolution near the interface was adaptive.

The FSI phenomena can be frequently observed in ocean and coastal engineering.The plates of marine structures tend to be thin, and the corresponding wall particles tend to be small. Therefore, the total number of particles will be large with a uniform resolution, which results in a heavy computational load. The multiresolution techniques were also introduced to the particle-based FSI solvers to improve their performance.Multiresolution techniques can be classified into two classes, namely, domain and structure refinement methods. Sun et al. (2019; 2021) coupled the SPH method with adaptive particle refinement (APR)to investigate the FSI problems. The APR region was set around the elastic structure.The coarse particles were split into fine particles when they entered the refinement region,and the fine particles were merged into large particles when they left the refinement region.Khayyer et al.(2019;2021) first introduced a multiresolution technique to MPS and ISPH-SPH to solve hydroelastic problems.The elastic structures were made of fine particles, whereas the flow field was composed of coarse particles. Sun et al. (2021)applied an adaptive high-resolution region to the MPSDEM coupled method(Sun et al.2019)to improve computational efficiency.

Figure 4 Sketch of the generation and removal of fluid particles in the buffer zone(Tsuruta et al.2021)

2.4.2 Central processing units(CPU)parallel acceleration technique

The MPS method cannot be applied to large-scale simulation only with multiresolution techniques. Scholars tried to assign the heavy task to multiple central processing units(CPUs)to increase computational efficiency.Two main parallel strategies, namely, open multiprocessing (OpenMP)and message passing interface (MPI), were utilized. For OpenMP, different processors share the memory, which can be easily coded.However,OpenMP is only applicable to a single host with multiple CPUs,and its speedup is limited because of the overhead for data communication(Luo and Koh 2017). Compared with OpenMP, MPI is suitable for largescale computing and can be applied to computer clusters.

The task division approach for particle-based methods can be classified as particle and domain decomposition methods.For the particle decomposition method,the particles are assigned to the processors according to their numbers. Massive data communication is required because the relative position between particles is not fixed (Iribe et al.2006). To overcome this drawback, Iribe et al. (2010) introduced the renumbering and communication list creation techniques to the particle decomposition method.The communication was significantly reduced. For the domain decomposition method, the calculation domain is divided into several subdomains, and only the particles near the boundary of subdomains exchange the information. Ikari and Gotch (2008) compared the aforementioned decomposition methods and pointed out that the domain decomposition method performed better than the particle decomposition method. Fernandes et al. (2015) proposed a new nongeometric strategy to parallelize the MPS method for fully distributed computing in clusters.The new strategy is based on the domain decomposition method and has been successfully applied to simulations with hundreds of millions of particles.Marrone et al.(2012)proposed a parallel strategy associated with background grids for SPH, which is also suitable for other particle-based methods.The particles are allocated to the background grids that belong to different processors,and the load balance is maintained by redistributing the grids to the processors dynamically. The parallel SPH was applied to the 3D simulation of a ship in steady forward motion. Xie et al. (2021a) introduced this parallel strategy to the MPS-DEM coupled method, which was used to solve the FSI problems.

2.4.3 GPU acceleration technique

Compared with the CPU, one chip area of the GPU has more arithmetic logic units.Therefore,the GPU can achieve high floating-point operations per second and deal with large amounts of data at the same time (Chen and Wan 2019a).CUDA is a general parallel computing platform released by NVIDIA, which enables the GPU to solve complex computing problems.The CUDA platform consists of the host and device parts.The CPU is the host part,whereas the GPU is the device part.

Recently, the GPU acceleration technique has been applied in many fields,including medical imaging and intelligent driving. The GPU acceleration technique has also been used to improve the efficiency of MPS code. Hori et al. (2011) developed a GPU-accelerated MPS method,which was three to seven times faster than the code with a single CPU. Chen et al. (2019) developed the in-house solver MPSGPU-SJTU, which was used to calculate the 3D dam-break problem with 2 678 595 particles. Compared with the computation time of the CPU, the total speedup of GPU code was up to 33.53.Zhang et al.(2021)used the GPU-based MPS to simulate the impinging atomization of two jets,and the speedup is 11.75.To further improve the computational efficiency, Gou et al. (2019) first introduced the multiple GPUs to the incompressible MPS.The domain decomposition method was used to divide tasks and the data communication between GPUs through the MPI. The multiple-GPU-based MPS solver was applied to the simulation of the 3D dam break with 1 million particles.

3 Applications in ocean engineering

3.1 Multiphase flows

3.1.1 Liquid-liquid and liquid-gas flows

The major challenges in extending the MPS method to liquid-liquid/liquid-gas flow problems focus on density and viscosity jump near the interface, surface tension, and compressibility.

Several approaches have been proposed to solve the interface problems. Koshizuka et al. (1999) proposed a multiphase MPS (MMPS) method and solved the pressure fields of two fluids separately to maintain numerical stability. To calculate the pressure head exactly, a buoyancy model was introduced to the MMPS method(Shirakawa et al. 2002; Kim and Kim 2014). However, the adjusting parameter in the buoyancy model was determined through several simple tests, which may be unsuitable for complex multiphase problems (Duan et al. 2017a). Several hybrid methods have also been developed to solve the multiphase problems. Liu et al. (2005) incorporated the MPS and finite volume method(FVM).The heavy phase was simulated by MPS, whereas the light phase was calculated by FVM. The information was exchanged between particles and meshes through extrapolation,which avoided the problems caused by the discontinuity of viscosity and density in the interface.Nevertheless,errors might be induced during the extrapolation process. Several researchers tried to solve the pressure fields of different fluids in one framework. Shakibaeinia and Jin (2012) developed the WCMPS solver, which was suitable for multiphase flow with a low density ratio. The multidensity and multiviscosity models were proposed to handle the interface problem. Khayyer and Gotoh (2013) proposed a high-order density smoothing scheme and extended the incompressible MPS to solve multiphase problems with a high density ratio. Shimizu et al.(2018)dealt with the interface with a high density jump by incorporating the space potential particle concept.Khayyer et al. (2019) further introduced the OPS to the MMPS to maintain the uniform distribution of the particles at the interface. Duan et al. (2017a) proposed two approaches,namely,MMPS-HD and MMPS-CA,to stabilize the multiphase simulations. In the MMPS-HD, harmonic mean density,which was adopted in mesh methods to handle high density ratios, was used to avoid large acceleration at the interface. In the MMPS-CA, new multiphase formulations are proposed based on the local weighted average of particle acceleration to keep the acceleration continuous at the interface.

The surface tension model is important in the simulation of multiphase flows. A continuum surface force (CSF)model was first proposed for mesh-based methods (Brackbill et al.1992).In particle-based methods,the surface tension is converted into the volume forces exerted on the interface particles. The surface tension is relevant to the local curvature at the interface. Nomura et al. (2001) calculated the curvature through the fan-shaped distribution of the particles near the interface, as shown in Figure 5.This approach of curvature approximation was also used by Rong and Chen (2010) to simulate Taylor bubble formation in microchannels.Based on the classical CSF,Duan et al.(2015)proposed the contoured continuum surface force(CCSF) model, which approximates the phase interface by the contours of the smoothed color function.The CCSF was robust enough and verified through several cases by Wen et al.(2021b).

Figure 5 Fan-shaped distribution of the particles near the interface

In the multiphase simulation with a high density ratio,the compressible characteristics cannot be ignored. Ikeda et al.(2001)assumed that density was correlated with pressure and related to sound velocity. Based on this assumption, a compressible term was added to the PPE source term. To reduce the pressure fluctuation, Khayyer et al.(2009) introduced a stabilizing term to the PPE source term, which enabled the slight compressibility of singlephase flows. Then, Khayyer and Gotoh (2016) further developed the compressible-incompressible MMPS method and proposed the compressible-incompressible error compensating source of PPE. The MMPS was used to investigate the cushioning (compressibility) effect of water slamming.Duan et al.(2020)proposed an incompressible-compressible particle method by incorporating incompressible MPS and weakly compressible SPH.

The MMPS has been applied to multiphase problems in ocean and coastal engineering assisted by the advanced techniques mentioned previously.

3.1.1.1 Bubbly flows

Bubbly flow is one of the research highlights in the field of ocean engineering. Tian et al. (2009) adopted MPS to investigate single bubble collapse with noncondensable gas and successfully captured the rebound phenomenon.Guo et al. (2018) introduced the high-order scheme and ECS to the MMPS-CA. The accuracy and stability of the improved MMPS-CA were verified by simulating bubble pairs rising and coalescence. Then, the improved MMPSCA was employed to investigate the evolution of bubbly flows in a simplified closed air-water separator.Wen et al.(2021b) investigated the violent interactions between a set of 3D bubbles with MPS, as shown in Figure 6.They concluded that bubble interactions in the 3D simulations became more violent,and more complex interfaces were captured by the MMPS method, particularly in cases with a large Bond number, when compared with 2D results. Pipeline transport is an important part of the offshore oil and gas exploitation process. Therefore, bubble motion in the pipeline needs to be investigated. Duan et al. (2017a) adopted MMPS-CA to simulate bubble motion in the rightangle pipe, as shown in Figure 7. The complex breakups or coalescence of the bubbles could be observed clearly,and the breakup of the gas phase frequently occurred in the right-angle region because of the existence of vortices.Rong and Chen(2010)simulated the generation of the Taylor bubble in a T-shaped microchannel using the MMPS method. Then, the influence of junction shape on bubble generation was analyzed.

3.1.1.2 Violent flows

Figure 6 Snapshots of the rising of 3D bubbly flows involving a set of bubbles(Wen et al.2021b)

Figure 7 Snapshots of simulated bubble rising in a right-angle tube with complex breakup behavior(Duan et al.2017a)

The phenomena of air entrapment mainly occur in the problems of violent flows, which affect the characteristics of the local fluid(Luo et al.2021).Therefore,the air phase should be considered in the simulation of violent flows.Wang and Zhang(2019)proposed a novel MMPS method.The spatial difference form of the governing equation of the original MPS method was replaced by the integral equation based on the theory of peridynamics, which can reduce the numerical instability induced by discontinuity of the physical quantities near the interface. Furthermore,a corrective matrix derived from the shape function in peridynamics was introduced to the gradient and divergence equation to eliminate the error induced by the uneven distribution of particles. The MMPS was applied to simulate the 2D dam break, as shown in Figure 8. The complex interface was successfully captured, and the numerical results are consistent with the experimental data. Wen et al.(2021a) simulated the 3D dam break benefiting from the GPU acceleration technique. Snapshots of 3D two-phase dam-break flows are presented in Figure 9. The numerical results indicated that the air cavity experienced severe pressure and volume oscillations, and the variation trends of pressure and volume were opposite. For multiphase sloshing cases, Shimizu et al. (2018) simulated liquid sloshing with a high filling ratio and small excitation amplitude.A smooth phase interface was captured,and no unphysical gap was observed. However, this method has not been applied to simulate more violent sloshing flows,which can be attributed to the fact that the pressures of different phases were solved separately and the existence of air entrapment that led to the discontinuity of the interface,thereby increasing the difficultly of simulation. Khayyer and Gotoh(2013)utilized the improved MPS with the density smoothing scheme and the Taylor series consistent pressure gradient model to simulate violent sloshing flows.The numerical results presented a smooth interface and uniform distribution of particles. Moreover, coalescence of sloshing-induced secondary jets,entrainment of gas particles, and fragmentation of water particles can be captured clearly. Wen et al. (2022) employed both the singlephase MPS and MMPS to simulate violent sloshing flows,as shown in Figure 10. The comparisons showed that, as the intensity of sloshing increased, the effect of the air cushion became nonnegligible, and the accuracy of the single-phase simulation significantly decreased, whereas the simulation results obtained by the MMPS method were consistent with the experimental results in all cases. Wen et al.(2021c)further investigated the violent multilayer liquid sloshing using the MMPS, as shown in Figure 11 and Figure 12. The phenomenon of internal wave breaking was well captured. Furthermore, the vertical baffle was able to significantly reduce the sloshing intensity and limit the interfacial wave elevations. Water entry is also a significant issue in the field of ocean engineering (Lyu et al.2022).The multiphase numerical models were also used to predict the pressure distribution and trajectory of structures accurately. Fang et al. (2022) simulated the water entry of the ship panel using a multiphase Riemann-SPH method.The air cushion and slamming load were investigated in detail. Lyu et al. (2021) investigated the water entry with different density ratios using the multiphase SPH model.However, the MPS has not been applied to simulate the multiphase water entry. In the future, the MMPS method should be further modified and employed to solve more complex problems.

Figure 9 Snapshots of 3D two-phase dam-break flows (Wen et al.2021a)

Figure 10 Comparison of air entrapment and cavity evolution captured in the experiment (Rafiee et al. 2010), MMPS simulation,and single-phase MPS simulation(Wen et al.2022)

Figure 11 Comparison of phase interfaces in the experiment(Audiffren et al.2012)and MMPS simulations(Wen et al.2021c)

Figure 12 Comparison of the time histories of interfacial wave elevations in the experiment (Audiffren et al. 2012) and MMPS simulations(Wen et al.2021c)

3.1.1.3 Wave problem

The original single-phase flow MPS method has been widely used in wave problems, which is naturally capable of handling strongly nonlinear free surfaces. In recent years,researchers have tried to employ the MMPS to solve wave problems related to multiphase flows. Khayyer et al.(2019) investigated multiphase standing gravity waves.The time histories of water surface elevation obtained by MMPS incorporated with OPS were consistent with the theoretical solution compared with other numerical results.Duan et al. (2017b) investigated the effect of waves and continuous spills on the oil spreading, as shown in Figure 13. The numerical results showed that the spreading of thick oil slicks could be accelerated when the wave height increased or wave length decreased.Wen et al.(2021a)analyzed internal solitary waves, as shown in Figure 14. The internal solitary wave was able to propagate at a long distance with the waveform maintaining its stability and the wave profile simulated by MMPS being consistent with the experimental data.

3.1.2 Liquid-solid flows

The solid phase of liquid-solid flows is considered the continuum field of discrete particles (Feng and Yu 2004).Considering the computational cost, the former is more suitable for large-scale problems. Shakibaeinia and Jin(2011) developed a multiphase non-Newtonian WCMPS to solve the mobile-bed dam-break problem. The liquid and solid phases were solved using a single set of multidensity and multiviscosity equations.Moreover,the generalized viscoplastic fluid model was used to predict the behavior of the solid phase.Fu and Jin(2015)introduced the Herschel-Bulkley rheological model to the WCMPS to simulate the motion of the solid phase, which was regarded as the non-Newtonian fluid. The deformation of landslides could be presented accurately. Tajnesaie et al. (2018)employed the viscoplastic rheological model to simulate the behavior of granular flows. Their numerical results showed that the WCMPS with viscoplastic rheological model performed better than the Herschel-Bulkley rheological model.Jandaghian et al.(2021)developed a 3D liquid-solid two-phase model based on the WCMPS with a regularized visco-inertial rheological model. A new formulation for the effective pressure of solid skeleton was proposed.The WCMPS was applied to simulate the immersed granular collapse,as shown in Figure 15.

Although the approach with the solid phase being a continuum is more efficient,some important information is neglected, such as the forces acting on a single particle and the motion behavior of a single particle. By contrast, the discrete method emphasizes the interactions among solid particles, and the simulation is more realistic.The DEM is a Lagrangian method based on Newton’s second law and is widely employed in simulations of granular flows. In recent years, MPS-DEM coupled methods have been developed for liquid-solid flows. The coupling strategies of MPS and DEM can be classified as resolved and unresolved coupling. For resolved coupling, the fluid force exerted on DEM particles is in the form of pressure and viscous force. Duan et al. (2022) developed a MPS-DEM model to simulate debris remelting in the field of nuclear engineering. Zhou et al. (2021) incorporated SPH and DEM. Then, the SPH-DEM model was applied to represent the tsunamis induced by the 1963 Vajont landslide by Xu et al. (2021). To the best of the authors’ knowledge,the MPS-DEM method with resolved coupling has not been applied to investigate the hydrodynamic problems in ocean and coastal engineering.

Figure 13 Snapshots of the distributions of oil slicks in the waves(Duan et al.2017b)

Figure 14 Snapshots of the propagation of internal solitary waves(Wen et al.2021a)

Figure 15 Snapshots of immersed granular collapse(Jandaghian et al.2021)

For unresolved coupling, Sakai et al. (2012) developed a 2D MPS-DEM model with the formulation of the fluid governing equations being Model A, which assumed that the pressure drop is shared between fluid and solid phases(Feng and Yu 2004).Sun et al.(2014)further coupled WCMPS and DEM with the same coupling scheme. Then, Li et al. (2019) extended the WCMPS-DEM coupled method to non-Newtonian solid-liquid flows.Xie et al.(2021b)incorporated the incompressible MPS with DEM with Model B, which assumes that pressure drop is applied to the fluid phase only (Feng and Yu 2004). The authors pointed out that the instability of DEM particles induced by fluid might be avoided because the pressure gradient force was replaced by a more stable term in Model B. The unresolved MPS-DEM model has been used to solve coastal engineering problems. Harada et al. (2019) simulated the swash beach process and investigated the effect of seepage flows on sediment transport. Tazaki et al. (2021) analyzed the vortical sorting process in the swash region using both experimental and MPS-DEM methods.Their numerical results were consistent with the experimental data.Harada et al. (2021) investigated the ripple formation process generated on a movable bed surface with 3D MPS-DEM model.Tazaki et al. (2022) assessed the effect of plunging waves on sediment transport,as shown in Figure 16.

3.2 Interaction between fluid and structures

3.2.1 Elastic structures

The fluid dynamics of severe flows may lead to large deformations and even fractures of offshore structures.Therefore,fluid-elastic structure interaction problems need to be solved.No fixed topological relationships among Lagrangian particles were detected, and information exchange was not restricted to specific nodes.Hence,particle-based methods have the potential to handle severe FSI problems with large deformations.

Particle-based methods can be directly used to simulate elastic structures. Hwang et al. (2014) developed a structural dynamic response model based on MPS.The governing equations of the structures were derived from Hook’s law for isotropic linear elastic solids.Khayyer et al.(2018a)employed the projection-based ISPH method for the simulation of fluid and developed a novel SPH-based structure model.The fluid-structure acceleration-based and pressure integration coupling schemes were compared in detail.Khayyer et al. (2018a) proposed Hamiltonian MPS and Hamiltonian SPH structural models based on the application of variational principles for elastodynamics.Then, the Hamiltonian structural models were extended to 3D(Khayyer et al. 2021a) and composite structures (Khayyer et al.2021c). Particle-based methods can also be coupled with other structural dynamics methods to solve the FSI problems. The finite element method (FEM), as a mature approach, has been widely applied to solve the dynamic problems of complex structures.Lee et al.(2007)first proposed the MPS-FEM coupled method, and the MITC4 shell element was adopted in the FEM analysis of the structure.Fourey et al.(2017)compared the conventional parallel staggered (CPS) and conventional sequential staggered(CCS) algorithms. The CSS algorithm was more stable than the CPS algorithm and was widely used as the coupling algorithm between MPS and other methods. Zhang et al. (2018; 2019) incorporated the incompressible MPS and FEM. The shape-function-based interpolation technique and kernel-function-based interpolation technique were proposed for the information exchange between fluid and structures. Sun et al. (2019) first coupled MPS and DEM.Then,Xie et al.(2021a)extended the MPS-DEM coupled method to the FSI problems with fracture.

Figure 16 Comparison of wave profile, sediment shape, and solid-phase velocity between experiments with PIV analysis and numerical simulations(Tazaki et al.2022)

The MPS and coupling methods have been applied to solve various FSI problems with deformation. Zha et al.(2021) employed the high-order MPS to simulate the 2D water entry of elastic wedges with different deadrise angles and drop heights.The 3D water entry of elastic wedges was simulated by Khayyer et al.(2021a),and their numerical results with different resolutions were consistent with the experimental data,as shown in Figure 17.Zhang et al.(2016a)analyzed the mitigating effects of the elastic baffle on sloshing using the MPS-FEM. Zhang and Wan (2018a) further investigated the phenomenon in the 2D elastic tank with an elastic wall.By varying Young’s modulus of the tank walls,interesting characteristics of the evolution of the free surface, variation of the impact pressures, and dynamic responses of the structures in both time and frequency domains were observed. Rao and Wan (2018) simulated the interactions between solitary waves and elastic plates.The effects of structural flexibility on wave-induced forces were analyzed. Zhang et al. (2022a) used the MPS-FEM to assess the hydroelastic responses of floating structures in waves,as shown in Figure 18 and Figure 19.A series of qualitative analyses were conducted, and the FSI solver could provide reasonable predictions for wave-ship interaction.

Figure 17 Snapshots of the 3D water entry of elastic structures(Khayyer et al.2021a)

Figure 18 Snapshots of the interaction between waves and floating structures(Zhang et al.2022a)

3.2.2 Floating structures

The floating structures may be damaged in the harsh marine environment,which could threaten the normal production and life of people. In recent years, many researchers tried to use MPS to investigate the interactions between waves and floating structures.

Figure 19 Time histories of impact pressure and ship deflection(Zhang et al.2022a)

Shibata et al. (2012) used MPS to investigate the ship’s motion in waves. Their numerical results showed that the MPS had the potential to simulate the ship’s motion under conditions with high wave height.Zhang et al.(2017)analyzed the roll motion of a rectangular floating body in regular waves with different frequencies. The numerical results of response amplitude operators for roll motions were consistent with the experimental data. Zhang et al. (2020)investigated the effect of the internal baffles on the flooding and motion characteristics of the damaged ship section. Hashimoto et al. (2022) compared MPS and SPH for the simulation of flooding ship sections, as shown in Figure 20. Both the MPS and SPH were able to simulate complex flooding flows.Amaro et al.(2021)developed an MPS-DEM coupled method, which was used to investigate the dynamic behavior of the wave-ice floes interaction.

Figure 20 Comparison of the simulation snapshots of flooding flows in ships by experiment,MPS,and SPH(Hashimoto et al.2022)

3.2.3 Porous structures

Porous structures can dissipate wave energy and are widely used to protect shorelines and bridge piers (Ren et al.2016).

Porous structure models have been introduced to particle-based methods by many researchers. A macroscopic model is an efficient approach that can address the demand for high accuracy,in which the N-S equation was reformulated based on the volume averaging theory (Xu and Jin 2019)or mixture theory(Khayyer et al.2018b).Compared with SPH, only a few studies of the MPS porous model have been conducted. Pahar and Dhar (2017) developed a model of flow-porous structure interaction-based incompressible MPS.A concept of apparent density was used to account for increased representative volume by multiplying a coefficient. The numerical solver was applied to solve free-surface flow passing uniform and nonuniform porous block cases. Xu et al. (2019) proposed a transit region between porous and nonporous regions for the smooth transition of porosity. The effects of porosity and mean diameter on the permeability, impact pressure, and reflected wave height were investigated.

3.2.4 Liquid sloshing

Because of the uneven distribution of energy in different regions,a large amount of energy needs to be transported from one area to another every year. Liquid sloshing is a significant issue in the transportation of liquefied natural gas, oil, and liquefied petroleum gas. The liquid inside a partially filled tank will oscillate violently and exert large impact pressure on the tank under external excitation, particularly when the excitation frequency is close to the natural frequency of sloshing.Particle-based methods can easily capture complex free surfaces with large deformation in tanks.In recent years,MPS has been widely applied to liquid sloshing problems.

MPS can be used to investigate the flow mechanics of sloshing. Pan et al. (2008) simulated 2D sloshing with MPS and investigated the effect of kernel functions on the numerical results. Pan et al. (2012) further extended MPS to large eddy simulation(LES)by incorporating it with the subparticle-scale turbulence model.The large impact pressure on the rolling tank wall was predicted accurately by MPS-LES. Jena and Biswal (2017) developed a modified MPS to simulate liquid sloshing in the 2D tank.The wave amplitude and pressure on the tank walls were calculated accurately.Xie et al.(2020)simulated the 3D liquid sloshing in a rectangular tank using the MPS methods assisted by GPU acceleration techniques. Their numerical results showed that the occurrence of the 3D effect was closely related to the excitation frequency, dimensions of the tank,and filling ratio.Zhang et al.(2022b)used MPS to investigated the mechanism of Faraday waves in a vertically excited square tank, as shown in Figure 21. The effects of forcing amplitude and excitation frequency on the formation of Faraday waves were evaluated.

Figure 21 Comparison of waveforms between experimental results(Frandsen and Peng 2006)and numerical results(Zhang et al.2022b)

MPS can also be employed to determine whether certain devices can mitigate sloshing. Tsukamoto et al. (2011)used MPS to analyze the influence of an elastically linked moving body on liquid sloshing. Their numerical results showed that the device can mitigate sloshing in the tank.Tsukamoto et al. (2020) simulated liquid sloshing with stiffeners on the bottom or side wall. For the filling ratio of 5%, the bottom stiffeners could significantly suppress sloshing,whereas the sidewall stiffeners have only a slight effect on sloshing. For the filling ratios of 10% to 20%,the sidewall stiffeners led to significant load mitigation.Bellezi et al.(2022)investigated the effects of floating baffles on liquid sloshing.The heave motion of the baffle was unfixed, and the baffle could significantly reduce liquid sloshing in a wide range of filling conditions.

4 Future perspectives and concluding remarks

The MPS method has been considerably improved and gradually applied to solve the hydrodynamic problems in ocean and coastal engineering since it was proposed. For stability and accuracy, the researchers have proposed many improved schemes for the gradient model,Laplacian model, and PPE source term. However, spurious pressure fluctuations still can be observed in the simulation by MPS, which affects the force exerted on offshore structures.Therefore,the reduction of spurious pressure fluctuations is a significant issue in future work.Furthermore,the stabilizing methods may introduce stabilization errors to the simulation (Duan et al. 2019). Hence, how to balance the accuracy and stability should also be considered.

Another major challenge of the MPS is the boundary condition. The shape of real offshore structures is complex, and slender structures are common. Future boundary conditions should provide adequate support and create accurate physical relationships for the fluids near the complex structures.The MPS method is more suitable for violent flows. However, the computational efficiency of MPS methods is relatively low compared with other methods,such as the FVM and high order spectral method.In the future, the MPS can be coupled with other methods to improve the computational efficiency of the simulation. The region where the waves interact with structures is simulated by MPS, and the wave information is transferred from other methods through OBC. Moreover, the conservation of mass and momentum should be considered during the information transfer process.

In terms of computational efficiency, the GPU technique can significantly improve the performance of MPS.Gou et al. (2019) conducted a simulation with multiple GPUs and determined that the total number of particles was approximately 1.24 million. The advanced task decomposition approaches can be proposed and introduced to fully utilize each GPU.With the further improvement of device memory and parallel strategy, the MPS simulation with particle numbers reaching 10 million should be tried in the future.

Regarding its application in ocean and coastal engineering, MPS was rarely used to solve complex problems. For the FSI problems,the safety of structures needs to be evaluated. MPS and coupled methods can be used to predict the responses of ships and bridges under extreme sea conditions. The stress distribution of overall structures should be presented accurately, which is conducive to the design of offshore structures. Furthermore, the fracture problems should be considered.Accurate fracture models and boundary conditions of fracture surfaces need to be proposed.The wave-ice interaction is also a hot issue, and MPS can be extended to model ice, such as SPH (Zhang et al.2019a). Both ice and water can be simulated by MPS,which may be beneficial for the information exchange between solid and fluid.As for liquid-solid flows,to simplify the numerical model, the shapes of solid particles are the same. However, granular flows in nature contain multiple nonspherical particles.Therefore, MPS simulation of solid particles with different shapes should be tried in the future.Moreover, the abrasion of solid particles is common. Several novel abrasion models can be introduced to the numerical methods to make the simulation more realistic.

The development of MPS and its application in ocean and coastal engineering are reviewed in detail. The latest advances in stability and accuracy, boundary conditions,and acceleration techniques of MPS methods are presented.The applications of MPS in ocean and coastal engineering,including multiphase flows and fluid-structure interaction, are discussed. Finally, the development direction of MPS methods is predicted. In summary, the MPS method has the potential to solve practical problems. New schemes and techniques need to be proposed to further promote the development of MPS.

AcknowledgementSupported by the National Key Research and Development Program of China (2019YFB1704200), and the National Natural Science Foundation of China(51879159,52131102).

Open AccessThis article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing,adaptation, distribution and reproduction in any medium or format,as long as you give appropriate credit to the original author(s) and thesource, provide a link to the Creative Commons licence, and indicateif changes were made. The images or other third party material in thisarticle are included in the article's Creative Commons licence, unlessindicated otherwise in a credit line to the material. If material is notincluded in the article's Creative Commons licence and your intendeduse is not permitted by statutory regulation or exceeds the permitteduse, you will need to obtain permission directly from the copyrightholder. To view a copy of this licence,visit http:/creativecommons.org/licenses/by/4.0/.