Improving stability of moving particle semi-implicit method by source terms based on time-scale correction of particle-level impulses
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: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: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 = P/ρgHF, 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)
- et al.
A class of renormalised meshless Laplacians for boundary value problems
J Comput Phys
(2018) - et al.
Two-dimensional fluid-structure impact computations with regularization
Comput Methods Appl Mech Eng
(1981) - et al.
Multi-phase SPH modelling of air effect on the dynamic flooding of a damaged cabin
Comput Fluids
(2018) - et al.
Improving stability of MPS method by a computational scheme based on conceptual particles
Comput Methods Appl Mech Eng
(2014) - et al.
An SPH projection method
J Comput Phys
(1999) - et al.
A multiphase MPS solver for modeling multi-fluid interaction with free surface and its application in oil spill
Comput Methods Appl Mech Eng
(2017) - et al.
The truncation and stabilization error in multiphase moving particle semi-implicit method based on corrective matrix: which is dominant?
Comput Fluids
(2019) - et al.
An improved high-order ISPH method for simulation of free-surface flows and convection heat transfer
Powder Technol
(2020) - et al.
Numerical simulation of Rayleigh-Bénard convection and three-phase Rayleigh-Taylor instability using a modified MPS method
Eng Anal Bound Elem
(2021) On enhancement of Incompressible SPH method for simulation of violent sloshing flows
Applied Ocean Research
(2014)