Improving stability of moving particle semi-implicit method by source terms based on time-scale correction of particle-level impulses

https://doi.org/10.1016/j.enganabound.2021.06.018Get rights and content

Highlights

  • Unstable pressure computation investigated from a totally new viewpoint of particle-level contact/collisions.

  • Proposition of correction based on momentum conservation of particle-level collisions.

  • New stable source terms are derived for pressure Poisson equation (PPE), mitigating pressure oscillations in a wider range of time and spatial resolutions.

  • Very simple implementation and no additional computational effort, which are desirable for practical engineering applications.

  • The approach can be combined with other strategies and incorporate into other projected-based particle methods.

  • Extensive verifications through hydrostatic and hydrodynamic benchmarks.

Abstract

The aim of this paper is to investigate the unstable nature of pressure computation focusing on incompressible flow modeling through the projection-based particle methods. A new approach from the original viewpoint of the momentum conservation regarding particle-level collisions, hereinafter refered to as time-scale correction of particle-level impulses (TCPI), is proposed to derive new source terms of pressure Poisson equation (PPE). This results in more stable computations with drastic reduction of unphysical pressure oscillations and more robust simulation with pressure magnitudes almost independent to time step. Moreover, compared to other strategies, no additional computational effort is required, its implementation is extremely simple, and the only numerical parameter is the propagation speed of the perturbations, of which the calibration is much more straightforward due to its physical meaning. Simulations were carried out using moving particle semi-implicit (MPS) method improved by the proposed approach. The comparisons of computed results with theoretical and experimental ones confirmed the effectiveness of the proposed approach.

Introduction

Many numerical methods have been proposed to predict hydrodynamic loads to solve multidisciplinary problems in hydraulic, coastal, marine, and offshore areas, and grid or mesh-based methods are generally applied for such problems. However, these methods suffer from some limitations, especially for problems involving complex or deformable domain, e.g., violent free-surface flow, and numerical treatments like free surface tracking or remeshing techniques are required, resulting in a very time-consuming activity. Due to the easy implementation and flexibility to model highly complex problems, the meshless methods have opened new perspectives in recent years, especially to solve violent free-surface flows and strong fluid-structure interaction (FSI) problems, see, e.g. Mazhar et al. [100]. An important class of meshless method is the particle-based methods, where the domain is represented by a collection of points (particles) with its own physical properties, such as mass and the internal/external forces, and of which the motions are evaluated by the interaction with the neighboring particles.

Numerical methods dealing with problems associated with strongly nonlinear interactions between water waves and structures are mainly based on solving either the fully nonlinear potential flow theory or the Navier-Stokes equations [59]. In general, particle-based methods applied to free-surface flows that solved the Navier-Stokes equations can be categorized into two groups: weakly compressible and incompressible projection-based methods. The former, such as smoothed particle hydrodynamics (SPH) [16,55] or weakly compressible moving particle semi-implicit (WC-MPS) [72] methods solve an appropriate equation of state in a fully explicit form. The latter, such as moving particle semi-implicit (MPS) [45] or incompressible smoothed particle hydrodynamics (ISPH) [10,54] methods solve a pressure Poisson equation (PPE) through a Helmholtz-Hodge decomposition [3] and application of the projection method [8,22,85]. Incompressible projection-based particle methods are generally expected to provide higher accuracy in terms of hydrodynamic pressure calculation and volume conservation [21,46][101]. Therefore, they are preferable in many applications focused on hydrodynamic load assessment. However, as the main challenges of these methods, compute-intensive and unstable pressure oscillation can be mentioned.

To achieve accurate and stable computations, substantial efforts have been carried out, as previously detailed by Khayyer and Gotoh [17], Gotoh and Khayyer [18], Ye et al. [97] and Li et al. [49]. These efforts covering at least one of the following approaches:

  • i.)

    regularization for particle distribution;

  • ii.)

    high-order and/or modified discrete differential operators;

  • iii.)

    improved free-surface particle detection;

  • iv.)

    new formulations for the source term of PPE;

  • v.)

    background mesh.

In the context of a particle regularization procedure to improve non-uniform particle distribution, several methods have been proposed. Monaghan [65] proposed the XSPH scheme which helps the particles to move with a velocity close to that of their neighboring particles, improving the smoothness of velocity field. Afterwards, Monaghan [64] highlighted the fact that the XSPH scheme does not conserve energy and proposed an implicit XSPH to resolve this issue. Xu et al. [96] proposed the particle shifting method, which slightly shifts the particles to prevent anisotropic particle structures, posteriorly extended by Lind et al. [50], allowing applications to free-surface flows. Tsuruta et al. [88] presented the so-called dynamic stabilization (DS) scheme which produces radial and anti-symmetric inter-particle forces and thus, at least, preserves both linear and angular momentum exactly. Further, Tsuruta et al. [89] introduced a potential in void space as space potential particle (SPP) to solve inconsistency in local volume conservation and to reproduce physical motions of particles near the free surface through a particle-void interaction. In the context of the ISPH method, an optimized particle shifting (OPS) scheme, free of adjusting parameters and allowing only the tangential shifting for free-surface particles, was presented by Khayyer et al. [38] and particle shifting technique combined to a collision model, named as hybrid particle shifting technique (HPST), was proposed by Garoosi and Shakibaeinia [14]. Khayyer et al. [39] also applied the OPS to an improved MPS in the context of multiphase flow. The implementation process was divided into two steps, where initially the heavy phase particles are shifted ignoring light phase particles and then both heavy and light phase particles are displaced. Despite better distribution, the updating process of particles position can be computationally expensive, which is not practical [93].

Instead of directly avoiding non-uniform particle distribution, modifications and corrections have been made on discrete differential operators considering irregular particle distribution. Discrete differential operators’ models have been proposed to assure momentum conservation [31,82], by changing the conventional operation of subtraction in the gradient operator. The so-called higher order Laplacian model was developed by Khayyer and Gotoh [33,35], by directly taking the divergence of a gradient model. Higher order approximations for differential operators were proposed based on a Taylor series expansion, giving a combination of these operators and corrective matrix [12,15,30,34,51,70,78,79]. However, generally these formulations are complex and demands additional computational cost to the compute-intensive methods, aside from the considerable coding implementation effort required.

Despite not being a definitive solution itself, the improvement of the free-surface particle detection techniques contributes significantly to more stable pressure computations by avoiding misdetection of free surface inside the fluid domain and the consequent mis imposition of the boundary condition for the pressure calculation. Among available fluid interface particle detection techniques, we can mention the techniques based on the sum of a weight function [45], amount of neighborhood particles [82], gradient vector of a property of the particles [29], relative distance between particles [19], a combination of the criteria based on the weight function and the amount of neighborhood particles [47] or the weighted distribution of neighboring particles [87], threshold angles between adjacent connecting lines of a given particle and its neighbor particles [77], or a combination of purely geometric sphere covering tests based on interval analysis with an adaptive spatial subdivision [71].

Regarding the source term of PPE, one of the considerable improvements was the combination of two error mitigating parts, the zero variation of the density and the velocity-divergence-free condition [25,82,98], which corresponds to the accumulative and instantaneous density deviations, respectively. In Hu and Adams [26], both the zero variation of density and velocity-divergence-free constraints of the incompressibility condition were enforced through the resolution of two PPEs and via application of a fractional time-step integration algorithm. Khayyer and Gotoh [32] and Khayyer et al. [37] introduced a higher order source term, which allows a slight compressibility. The authors also proposed a different formulation for the calculation of density variation by the total differentiation of the weight function. Kondo and Koshizuka [43] proposed a multi-term source consisting of one main part and two error-compensating parts. With proper selection of the coefficients multiplied by the error-compensating parts, unphysical pressure oscillation was suppressed. Khayyer and Gotoh [34] and Gotoh et al. [20] introduced a modified source term comprising of a set of high order error compensating, similar to the work of Kondo and Koshizuka [43], but with dynamic coefficients as functions of the instantaneous flow field. Following the idea of dynamic coefficients, Tamai and Koshizuka [79] and Sun et al. [76] also derived improved source terms and demonstrated that the adoption of dynamic coefficients mitigates the non-physical pressure oscillations. Zheng et al. [99] proposed a new formulation of incompressible flows based on Rankine source solution. The PPE was modified into a form that does not require any direct approximations for function derivatives, then leading to a more robust numerical method, of which the main advantage is no need to approximate second-order derivatives in the PPE. Ngo-Cong et al. [67] proposed a novel numerical approach for incompressible smoothed particle hydrodynamics by solving the PPE on a set of so-called moving integrated radial basis function networks.

Apart from aforementioned improvements, some researchers have used background meshes to mitigate problems due to the particle's maldistribution. Hwang [27] and Ng et al. [68] proposed an accurate MPS based on an embedded pressure mesh, namely moving particle pressure mesh (MPPM) and unstructured moving particle pressure mesh (UMPPM), to handle the continuity constraint. In both works, the main idea is to consider the pressure as Eulerian variable, which is in contrast with the original MPS whereby the pressure and velocity are computed on the Lagrangian particles. A background mesh scheme was proposed by Wang et al. [92] aiming to an accurate, smooth and spatially continuous source term for the PPE. The source terms are first calculated at background mesh nodes with the interpolated velocities at fixed and regular neighboring nodes. Then, the computed source terms at mesh nodes are extrapolated to the nearest neighboring particles. Different from the last two works, the imaginary background meshes are only used for the calculation of the source terms of PPE, while other calculations are performed in the Lagrangian framework. The effectiveness of the approach for the enhancement of the stability and accuracy of pressure calculation is showed through several benchmark tests. The incorporation of a background mesh into Lagrangian particle-based methods seems a promising direction to solve challenging issues related to accurate and stable computations, but it may lead to complex data interpolations.

Although several efforts have been undertaken to enhance the computation of the hydrodynamic pressure, the source terms of previous works depend on the time step and, in practice, the computed pressure peak magnitudes are sensitive to time step, with an unstable behavior of increasing magnitude as the time step decreases. In summary, since the solution does not converge as time step goes to zero, these formulations are not rigorously consistent in time domain.

The important topics of consistency and stability, i.e. convergence, in particle-based methods have been addressed mainly on the spatial discretization [1,28,52,56,66,74,75,80]. To the best knowledge of authors, only in Chen et al. [7] the pressure oscillation associated to the temporal discretization was investigated by using the so-called NSD-MPS, in which truncated compact support are supplemented by conceptual/virtual particles with pressures solved from a modified PPE. The NSD-MPS prevents the clustering of zero-pressure free-surface particles, and, as consequence, the non-physical pressure oscillations could be significantly mitigated even when the time step is reduced. Within this context and aiming to provide a practical engineering tool for more time-stable assessment of loads due to extremely nonlinear hydrodynamic phenomena, an improvement on the pressure computation of the incompressible projection-based particle methods is proposed in the present work. Based on a novel viewpoint of momentum conservation in particle-level collisions, a correction of the mismatch between numerical and physical time scales is introduced to derive new formulations for the PPE. This leads to original source terms that depend directly on the spatial discretization and independent to time step, which results in time-stable computations of the pressure. Furthermore, the proposed formulation preserves advantages of the original meshless particle-based numerical methods such as their flexibility. Moreover, it is very easy in coding implementation, requires no additional computational effort on the already compute-intensive method and can be incorporated to the above mentioned approaches to improve the stability of the pressure assessment. To demonstrate the effectiveness of the new PPE source terms derived from the proposed approach, namely time-scale correction of particle-level impulses (TCPI), an in-house code based on MPS method is adopted and four representative phenomena are considered: hydrostatic condition, water jet hitting a perpendicular flat wall, dam breaking, and liquid sloshing inside a prismatic tank. The computed pressures from MPS with original and proposed source terms are compared to show the effectiveness of the proposed improvement. The numerical results are also compared against the analytical and experimental ones available in the literature or conducted by the authors. In the Appendix A, 2D standing wave simulations are conducted to demonstrate that the conservation features (volume and energy) are not deteriorated due the inclusion of the proposed source terms.

Section snippets

Governing equations

The Lagrangian form of the conservation laws of mass and momentum for incompressible viscous flow are:DρDt=ρ·u=0,DuDt=Pρ+νk2u+f,where ρ is the fluid density, u denotes the velocity vector, P represents the pressure, νk stands for the kinematic viscosity and f is the vector of the external body force per unit mass.

In the MPS method, the differential operators of the governing equations are replaced by discrete differential operators on irregular nodes derived from a positive compactly

Time‐scale Correction of Particle‐level Impulses (TCPI)

In explicit numerical scheme, specifically for phenomena dominated by advection effects, the propagation of perturbations is limited for numerical stability by the Courant-Friedrichs-Lewy (CFL) condition [9], which establishes the relation between the propagation velocity of the perturbations cs and the numerical parameters:Cr=csΔtl0<1.0,where, Cr is the Courant number, Δt is the numerical time step and l0 is the spatial resolution.

For time discontinuous phenomena, such as impulsive load due to

Calibration of the propagation speed of the perturbations cs for TCPI‐PND and TCPI‐PND‐DF source terms

In this section, one is intended to find the range of the speed of the perturbations cs that enables stable simulations. For this purpose, the fluid volume variation and the pressure time series computed by the proposed source terms TCPI-PND and TCPI-PND-DF are analyzed considering hydrostatic condition. Figure 7 illustrates the tank dimensions of length LT = 30 × l0, initial water level height HF = 20m, and the position of the pressure sensor SH1 located at the tank bottom. Even in the

Cases of Study

In order to evaluate the performance of the source terms derived from the proposed TCPI approach, four phenomena were considered. First, a 2D hydrostatic case is simulated and the pressure of the particle on the bottom of the tank computed by the proposed source terms is compared to original ones and analytical results. After that, three simulations covering dynamic cases are considered. The first one is 2D water jet impact on a flat wall. The analytical pressure at the stagnation point and

2D hydrostatic tank

The first test case consists of a 2D hydrostatic condition with tank of length LT = 0.36m and initial column water height HF = 0.48m. The fluid properties are density of ρ = 1000kg/m3 and kinematic viscosity of νk = 1.0 x 10−6m2/s. The computed hydrostatic pressures at the sensor SH1, placed at the bottom of the tank, using different source terms were analyzed. The hydrostatic pressure is represented by the non-dimensional pressure coefficient CP = PgHF, and the dimensionless time (τ) is

Concluding Remarks

In the present study, the spurious numerical oscillations associated to the Lagrangian particle-based simulations are analyzed through the time-stability issue of the computed pressure, and an original interpretation based on the numerical modeling issues of the particle-level collision/impact is presented. By considering the momentum conservation of the particle-level collision/impact, a correction between numerical and physical duration of the impulses is proposed. Then, relaying on so-called

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Acknowledgments

This material is based upon research supported by the Office of Naval Research Global under the Award Number: N62909-16-1-2181. This study was financed in part by the Coordenação de Aperfeiçoamento de Pessoal de Nível Superior - Brasil (CAPES) - Finance Code 001. The authors are also grateful to Petrobras for financial support on the development of the MPS/TPN-USP simulation system based on MPS method and Frade Japão Petróleo Ltda for the financial support during the performance of the

References (101)

  • X.Y. Hu et al.

    An incompressible multi-phase SPH method

    J Comput Phys

    (2007)
  • A. Khayyer et al.

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

    Coastal Engineering

    (2009)
  • A. Khayyer et al.

    A higher order Laplacian model for enhancement and stabilization of pressure calculation by the MPS method

    Applied Ocean Research

    (2010)
  • A. Khayyer et al.

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

    J Comput Phys

    (2011)
  • A. Khayyer et al.

    A 3D higher order laplacian model for enhancement and stabilization of pressure calculation in 3D MPS-based simulations

    Applied Ocean Research

    (2012)
  • A. Khayyer et al.

    Enhancement of performance and stability of MPS mesh-free particle method for multiphase flows characterized by high density ratios

    J Comput Phys

    (2013)
  • A. Khayyer et al.

    Enhanced predictions of wave impact pressure by improved incompressible SPH methods

    Applied Ocean Research

    (2009)
  • A. Khayyer et al.

    Comparative study on accuracy and conservation properties of two particle regularization schemes and proposal of an optimized particle shifting schem in ISPH context

    J Comput Phys

    (2017)
  • A. Khayyer et al.

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

    Comput Fluids

    (2019)
  • A. Khayyer et al.

    On enhancement of energy conservation properties of projection-based particle methods

    European Journal of Mechanics - B/Fluids

    (2017)
  • A. Khayyer et al.

    Multi-resolution MPS for incompressible fluid-elastic structure interactions in ocean engineering

    Applied Ocean Research

    (2019)
  • E.-.S. Lee

    Comparisons of weakly compressible and truly incompressible algorithms for the SPH mesh free particle method

    J Comput Phys

    (2008)
  • B.-.H. Lee et al.

    Step-by-step improvement of MPS method in simulating violent free-surface motions and impact-loads

    Comput Methods Appl Mech Eng

    (2011)
  • S.J. Lind et al.

    Incompressible smoothed particle hydrodynamics for free-surface flows: a generalised diffusion-based algorithm for stability and validations for impulsive flows and propagating waves

    J Comput Phys

    (2012)
  • X. Liu et al.

    An advanced moving particle semi-implicit method for accurate and stable simulation of incompressible flows

    Comput Methods Appl Mech Eng

    (2018)
  • X. Liu et al.

    A stable moving particle semi-implicit method with renormalized Laplacian model improved for incompressible free-surface flows

    Comput Methods Appl Mech Eng

    (2019)
  • L. Lobovský

    Experimental investigation of dynamic pressure loads during dam break

    J Fluids Struct

    (2014)
  • S. Marrone et al.

    Prediction of energy losses in water impacts using incompressible and weakly compressible models

    J Fluids Struct

    (2015)
  • D.D. Meringolo et al.

    On the filtering of acoustic components in weakly-compressible SPH simulations

    J Fluids Struct

    (2017)
  • C. Mokrani et al.

    Conditions for peak pressure stability in VOF simulations of dam break flow impact

    J Fluids Struct

    (2016)
  • K.C. Ng et al.

    On the accuracy assessment of Laplacian models in MPS

    Comput Phys Commun

    (2014)
  • D. Ngo-Cong et al.

    Incompressible smoothed particle hydrodynamics-moving IRBFN method for viscous flow problems

    Eng Anal Bound Elem

    (2015)
  • K.C. Ng et al.

    Unstructured moving particle pressure mesh (UMPPM) method for incompressible isothermal and non-isothermal flow computation

    Comput Methods Appl Mech Eng

    (2016)
  • P.W. Randles et al.

    Smoothed particle hydrodynamics: some recent improvements and applications

    Comput Methods Appl Mech Eng

    (1996)
  • K. Shibata et al.

    The overlapping particle technique for multi-resolution simulation of particle methods

    Comput Methods Appl Mech Eng

    (2017)
  • L.D.G. Sigalotti

    A new insight into the consistency of the SPH interpolation formula

    Appl Math Comput

    (2019)
  • A. Souto-Iglesias et al.

    On the consistency of MPS

    Comput Phys Commun

    (2013)
  • Z. Sun et al.

    Modified MPS method for the 2D fluid structure interaction problem with free surface

    Comput Fluids

    (2015)
  • C. Sun et al.

    Surface treatment technique of MPS method for free surface flows

    Engineering Analysis with Boundary Elements,

    (2019)
  • T. Tamai et al.

    On the consistency and convergence of particle-based meshfree discretization schemes for the Laplace operator

    Comput Fluids

    (2017)
  • M. Tanaka et al.

    Multi-resolution MPS method

    J Comput Phys

    (2018)
  • M. Tanaka et al.

    Stabilization and smoothing of pressure in MPS method by Quasi-Compressibility

    J Comput Phys

    (2010)
  • N. Tsuruta et al.

    A short note on dynamic stabilization of moving particle semi-implicit method

    Comput Fluids

    (2013)
  • D. Violeau et al.

    Optimal time step for incompressible SPH

    J Comput Phys

    (2015)
  • L. Wang

    Enhancement of pressure calculation in projection-based particle methods by incorporation of background mesh scheme

    Applied Ocean Research

    (2019)
  • J. Wang et al.

    Modified particle pethod with integral Navier–Stokes formulation for incompressible flows

    J Comput Phys

    (2018)
  • G.X. Wu et al.

    Finite element analysis of two-dimensional non-linear transient water waves

    Applied Ocean Research

    (1994)
  • R. Xu et al.

    Accuracy and stability in incompressible SPH (ISPH) based on the projection method and a new approach

    J Comput Phys

    (2009)
  • X. Zheng et al.

    Incompressible SPH method based on Rankine source solution for violent water wave simulation

    J Comput Phys

    (2014)
  • Farrukh Mazhar et al.

    On the meshfree particle methods for fluid-structure interaction problems

    Engineering Analysis with Boundary Elements

    (2021)
  • Cited by (0)

    View full text