Smoothed particle hydrodynamics and its applications in fluid-structure interactions

187

0

2017,29(2):187-216

DOI: 10.1016/S1001-6058(16)60730-8

Smoothed particle hydrodynamics and its applications in fluid-structure interactions* A-man Zhang (张阿漫)1, Peng-nan Sun (孙鹏楠)1,2, Fu-ren Ming (明付仁)1, A. Colagrossi2 1. College of Shipbuilding Engineering, Harbin Engineering University, Harbin 150001, China, E-mail: [email protected] 2. CNR-INSEAN, Marine Technology Research Institute, Rome, Italy (Received November 24, 2016, Revised January 25, 2017) Abstract: In ocean engineering, the applications are usually related to a free surface which brings so many interesting physical phenomena (e.g. water waves, impacts, splashing jets, etc.). To model these complex free surface flows is a tough and challenging task for most computational fluid dynamics (CFD) solvers which work in the Eulerian framework. As a Lagrangian and meshless method, smoothed particle hydrodynamics (SPH) offers a convenient tracking for different complex boundaries and a straightforward satisfaction for different boundary conditions. Therefore SPH is robust in modeling complex hydrodynamic problems characterized by free surface boundaries, multiphase interfaces or material discontinuities. Along with the rapid development of the SPH theory, related numerical techniques and high-performance computing technologies, SPH has not only attracted much attention in the academic community, but also gradually gained wide applications in industrial circles. This paper is dedicated to a review of the recent developments of SPH method and its typical applications in fluid-structure interactions in ocean engineering. Different numerical techniques for improving numerical accuracy, satisfying different boundary conditions, improving computational efficiency, suppressing pressure fluctuations and preventing the tensile instability, etc., are introduced. In the numerical results, various typical fluid-structure interaction problems or multiphase problems in ocean engineering are described, modeled and validated. The prospective developments of SPH in ocean engineering are also discussed. Key words: Smoothed particle hydrodynamics, ocean engineering, fluid-structure interaction, bubble dynamics, underwater explosion, hydrodynamics

Introduction Smoothed particle hydrodynamics (SPH) method and its related variants have been developing as a new generation of computational fluid mechanics (CFD) solvers for complex hydrodynamic problems for decades[1-9] and have been widely applied in a wide range of engineering problems[10-17]. The present work is devoted to a review of the recent developments and achievements of SPH methodology and an introdu* Project supported by the National Natural Science Foundation of China (Grant Nos. U1430236, 51609049), the China Postdoctoral Science Foundation (Grant No. 2015M581432) and the China Scholarship Council (CSC, Grant No. 201506680004). Biography: A-man Zhang (1981-), Male, Ph. D., Professor

ction to its typical applications in fluid-structure interactions in ocean engineering. The distinctive features of the problems in ocean engineering are mostly linked to the existence of a free fluid surface (when the air phase is neglected) and moving or deforming wall boundaries. To deal with these kinds of problems, different numerical solvers based on Eulerian, Lagrangian or Arbitrary Lagrangian Eulerian (ALE) formulations are developed[18-21]. These numerical solvers are mesh-based or meshless. The Eulerian mesh-based methods may demand complex algorithms and huge computational costs in the tracking of the free surface, especially in simulations with drastic free surface splashing jets. While in the Lagrangian mesh-based solvers, the fluid surface and boundary are explicitly tracked. However, the problems of the Lagrangian mesh distortion and re-mesh in some violent fluid-structure interaction problems may become cumbersome. By contrast with the mesh-based

188

methods, in Lagrangian meshless methods, e.g. SPH[22] or Moving Particle Semi-implicit method (MPS)[23], inherently the complex algorithm of the free surface tracking and the mesh distortion are avoided. The Lagrangian motions of the material particles also supply a convenient way for the analysis of the flow features in a Lagrangian way, which is another merit compared to the Eulerian methods[24,25]. As SPH is a distinct representative of Lagrangian meshless methods and it has been widely applied in the field of hydrodynamics recently. Moreover, it has been proven in Souto-Iglesias et al.[26,27] and Shao and Gotho[28] that SPH and MPS are closely related, or in other words, equivalent, hereinafter our review and discussion are mainly based on SPH and its applications in fluid-structure interaction problems in ocean engineering. In the present review, the problems solved by SPH are classified into two categories namely hydrodynamic problems (e.g. water-body interactions, water impacts, water entries and exits, flow past bodies, etc.) and explosion problems (e.g. underwater explosions, contact explosions, etc.). The two categories of problems have obviously different features. For the former, the conservation of mass, momentum and energy of the numerical solver is crucial for maintaining the numerical accuracy and stability in a long-term simulation, while for the latter, the fluid should be considered as compressible and a large transient and violent deformation of the structure (fluid boundary) are usually involved. For long-term hydrodynamic problems, expensive computational costs are usually encountered. To obtain a fast and efficient simulation, on the one hand the adoption of advanced computing platforms can be crucial, recently different high performance computing techniques, e.g. OpenMP[13,29], MPI[30,31], CUDA[32-34], OpenCL[11] are introduced and applied in the framework of SPH method to accelerate the computing. On the other hand, different algorithms are also proposed to reduce the amount of computation by SPH itself. In the weakly-compressible SPH (WCSPH) method, e.g. the well-known δ - SPH scheme, a fourth-order Runge-Kutta time-integration is adopted, which allows for a relatively large CFL coefficient for the time step[35-37]. While in the Incompressible SPH (ISPH) method[38-41], the time-step is not restricted by the stability reason but by the accuracy reason[42,43]. The time step of ISPH can be 5 times larger than WCSPH, however in each time step, the ISPH requires a more expensive computational cost due to the tracking of the free surface particles[44] and the solving of the Poisson equation[45]. Finally, after a compromise between the stability and accuracy, an equivalent efficiency can be achieved between WCSPH and ISPH[43]. Another cause that restricts the efficiency of SPH is

the request for a uniform particle distribution to satisfy the conditions in the particle approximation, see the theoretical analysis in Liu and Liu[46]. As has been recently explored in Antuono et al.[47], only when the particle distribution is in a sufficiently low disorder, a convergence rate between 1 and 2 can be obtained. This characteristic means that in SPH the multi-resolution for the particle scale cannot be as straightforward as the adaptive mesh-refinement widely used in the mesh-based CFD solvers[48]. However, recently different numerical techniques are developed to address the adaptive particle refinement (APR)[49-51], in which not only the total computational cost is reduced but also the local accuracy is drastically improved. The achievements of APR in SPH make the large scale simulation, which is quite common in ocean engineering applications, become possible. The explosion problem is another important branch in the fluid-structure interaction problems in ocean engineering, e.g. the high-velocity impact[52], underwater explosion[53-56], contact explosion[16], etc. In the simulation of these problems SPH can be more reliable than other mesh-based solvers[46,56,57]. A proper equation of state is important for these problems since the fluid phases are usually regarded as compressible. Besides, SPH approximation formulations of high accuracy is also crucial since a large deformation of the material is usually involved, such as the fluid surface splashing, fragmentations, moving discontinuities, large inhomogeneity and structure fractures. In these circumstances, perhaps the particles distributions are quite disordered or the kernel functions near the boundary can be truncated. Therefore, recently SPH equations of high accuracy[58-61] are developed to resolve these problems. To model the structure responses in fluid-structure interactions, different methods have been developed. The first way is to adopt an SPH-FEM coupling[52,53,62-65]. The second way is to employ another particle based method namely shell SPH to model the structure[66,67]. The third way is to adopt an elastic SPH model to model the response of an elastic structure. Recently, this technique has been widely adopted in the sloshing or dam breaking problems considering an elastic baffle[68,69]. Finally, it is worth mentioning in Yang et al.[70-72], a coupling method of SPH and element bending group (EBG) is developed for modeling the interaction between viscous flows and flexible structures with free or fixed ends. The consideration of the air phase above the fluid surface is also highlighted in the present work. Indeed, in most of the numerical simulations for free surface flows, the air effect is not considered since it is difficult to treat and also brings more expensive computational cost. Recently, in Marrone et al.[73], it is pointed out that energy evolutions can be similar when inclu-

189

ding or not the presence of air in a fluid-structure simulation. However, from the experimental and numerical studies[74,75], the air effect can be crucial in some problems, e.g. the effect of the bubble oscillation on the load induced by liquid sloshing[76,77] or dam breaking[78], the cavitation effect on the denudation of hydraulic equipment[79,80] and rising bubble problems in the marine resource exploitations[17,81,82]. In the present work, multi-phase SPH formulations[41,83-85] are reviewed and the surface tension effect is considered[85,86]. A 3-D rising bubble benchmark case will be considered as an example to validate the multiphase SPH scheme. 1. The numerical models 1.1 Kernel and particle approximations 1.1.1 Classic SPH formulations In SPH, a generic function f (r ) and its derivatives at a general position denoted by the vector r can be approximated with the kernel approximation by integrating over the compact domain Ω as f (r ) = ∫ f (r ′)W (r − r ′, h)d r ′

(1a)

∇ r f (r ) = ∫ f (r ′)∇ rW (r − r ′, h)d r ′

(1b)

W

W

where ∇ r denotes the gradient with respect to r . r ′ denotes the positions all over Ω . W (r − r ′, h) is the kernel function where h is named as the smoothing length. The radius of Ω is equal to an integer multiple of h depending on the choice of kernel functions. The kernel function must satisfy a number of conditions, namely an even function, the normalization condition, the Dirac delta function property and the compact condition, for more details see Liu and Liu[46]. Due to the kernel properties, we have ∫ f (r ) ⋅ Ω

∇ rW (r − r ′, h)d r ′ = 0 . Therefore, the derivatives of f (r ) can be written in the following two forms:

∇ r f (r ) = ∫ [ f (r ′) − f (r )]∇ rW (r − r ′, h)d r ′

(2a)

∇ r f (r ) = ∫ [ f (r ′) + f (r )]∇ rW (r − r ′, h)d r ′

(2b)

W

W

In these Eqs.(2), we underline that the second equation is more suitable for the discretization of the momentum equation since it allows an anti-symmetric property between r and r ′ which is crucial for maintaining the momentum conservation when f (r ) represents pressure.

Similarly, the divergence of vector f ( x ) can also be written as ∇ r ⋅ f (r ) = ∫ [ f (r ′) − f (r )] ⋅ ∇ rW (r − r ′, h)d r ′

(3a)

∇ r ⋅ f (r ) = ∫ [ f (r ′) + f (r )] ⋅ ∇ rW (r − r ′, h)d r ′

(3b)

W

W

In spite of the applicability of the second Eqs.(2) to the momentum equation, the first equation is more suitable for the discretization of the continuum equation. The reason is owing to the existence of a free surface which makes the kernel function truncated, the first equation decreases the errors introduced by the kernel truncations, see more analysis in Colagrossi et al.[87]. In SPH, the domain is usually discretized into uniformly distributed particles. The particles carry all the generic variables, e.g. pressure, density, velocity, acceleration etc.. The continuous integration can be approximated with the sum of finite particle volumes. This process is named particle approximation[46]. Therefore, after a selection for these SPH approximation equations, the typical kernel approximations can be rewritten in a particle approximation form as follows:

f (r i ) = ∑ f (rj )W (ri − rj , h)V j

(4a)

∇i f (ri ) = ∑ [ f (rj ) + f (rj )]∇iW (ri − rj , h)V j

(4b)

∇i ⋅ f (ri ) = ∑ [ f (rj ) − f (rj )] ⋅ ∇iW (ri − rj , h)V j

(4c)

j

j

j

where the subscript i denotes the i - th particle, j denotes the particle within the compactly supported domain Ω of particle i . V j denotes the volume of particle j . V j can be calculated by m j / ρ j , where m and ρ denote particle mass and density respectively. The selection of kernel function directly affects the computational efficiency, accuracy and stability of SPH simulations. An improved Gaussian kernel function[84] was suggested for the simulation of free surface flows and multiphase flows[17,82]. However, as analyzed in Dehnen and Aly[88], the increase of N H which is the number of neighbor particles may lead to pairing (or clumping) instability. Dehnen and Aly have recommended a series of kernels named Wendland kernels. Regarding the determination of the smoothing length, usually a ratio κ = h / ∆x ( ∆x is the initial particle spacing) is adopted. We underline that for 2-D

190

problems, 50 particles in the neighboring domain is enough for an accurate simulation. It corresponds to κ = 2 for the Wendland C2 kernel. In 3-D simulations, 113 neighboring particles can be enough, which corresponds to κ = 1.5 for the Wendland C2 kernel.

If we multiply both sides of Eq.(5) with ∂W (r ′ − r , h) / ∂x and ∂W (r ′ − r , h) / ∂y , integrating on Ω and omitting the terms with seconder order derivatives, we have[60,61]

1.1.2 SPH formulations of higher order of accuracy In some cases, the SPH approximations introduced above are not accurate enough or the errors converge to zero slowly, especially in the cases with kernel truncations or disordered particle distributions. In this part, SPH formulations with higher order of accuracy are reviewed. Here we suggest several different variants: normalization for the divergence approximation[89], corrected smoothed particle hydrodynamics (CSPH)[60,90], finite particle method (FPM)[7,59] and kernel gradient free SPH (KGF-SPH)[91]. The similarity of these variants is they can be derived by the adoption of the Taylor series. Here we use a 2-D case for instance. The Taylor series expansion with the first order completeness is

f (r ′) d r ′ = f (r ) ∫ ∫WW ∂x

f (r ′) = f (r ) +

∂ f (r ) ∂ f (r ) ( x′ − x) + ( y ′ − y ) + O (h 2 ) ∂x ∂y (5)

where x and y are the two components in x and y directions of the vector r . (1) The corrected smoothed particle hydrodynamics (CSPH) CSPH operators have been widely applied in SPH simulations. A typical one is the evaluation of the gradient of the density used in the δ - term in the wellknown δ - SPH scheme[92]. Similar to what was done in Jiang et al.[60], if we multiply both sides of Eq.(5) with W (r ′ − r , h) , integrate on Ω and we have

∫WW

∂ f (r ) ( x′ − x)W (r ′ − r , h)d r ′ + ∂x ∫W

(6)

By omitting the terms with derivatives of f (r ) , we obtain the interpolation formulation for f (r ) with the first order accuracy as

f (r ) =

∫W f (r ′)W (r ′ − r , h)d r ′ ∫W W (r ′ − r , h)d r ′

∂W (r ′ − r , h) dr′ + ∂x

∂ f (r ) ∂W (r ′ − r , h) ( x′ − x) dr′ + ∫ W ∂x ∂x

∂ f (r ) ∂W (r ′ − r , h) ( y′ − y) d r ′ + O (h 2 ) ∫ ∂y W ∂x

(8)

and ∂W (r ′ − r , h)

f (r ′) d r ′ = f (r ) ∫ ∫WW ∂y

∂W (r ′ − r , h) dr′ + ∂y

∂ f (r ) ∂W (r ′ − r , h) ( x′ − x) dr′ + ∫ W ∂x ∂y ∂ f (r ) ∂W (r ′ − r , h) ( y′ − y) d r ′ + O (h 2 ) ∫ W ∂y ∂y

(9)

Solving Eq.(8) and Eq.(9) together we have  ∂ f (r )   ∂W (r ′ − r , h)   ∂x    ∂x   = ∫ [ f (r ′) − f (r )] M −1   dr′  ∂ f (r )  W  ∂W (r ′ − r , h)  ∂y  ∂y    (10)

where ∂W (r ′ − r , h)  dr′  ∫W ( x′ − x) ∂x M =  ( x′ − x) ∂W (r ′ − r , h) d r ′  ∫Ω ∂y

f (r ′)W (r ′ − r , h)d r ′ = f (r ) ∫ W (r ′ − r , h)d r ′ +

∂ f (r ) ( y ′ − y )W (r ′ − r , h)d r ′ + O (h 2 ) ∂y ∫W

∂W (r ′ − r , h)

(7)

Equation (7) is often used to filter the density fluctuation in WCSPH[93].

∂W (r ′ − r , h)  d r ′ ∂x  ∂W (r ′ − r , h)  dr′ ∫W ( y′ − y)  ∂y

∫W ( y′ − y)

(11)

And therefore we have corrected kernel derivatives allowing a second order accuracy as follows  ∂W C (r ′ − r , h)   ∂W (r ′ − r , h)      ∂x ∂x   = M −1   C  ∂W (r ′ − r , h)   ∂W (r ′ − r , h)      ∂y ∂y  

(12)

191

Equation (10) has the second order accuracy even for a truncated kernel function or in a disordered particle distribution. Since the Eq.(12) cannot ensure a symmetric property between r and r ′ , it is not often adopted for the discretization of the governing equations. However Eq.(10) is often employed for the calculation of the gradients, for example in vorticity, diffusive terms or deformation gradients in the evaluation of finite time lyapunov exponent (FTLE), see e.g. in Sun et al.[24]. (2) The finite particle method (FPM) FPM can be regarded as an improvement of CSPH introduced in the previous subsection[58,59]. In FPM, for a 2-D case, Eqs.(6), (8) and (9) are solved together. With three equations and three unknowns, we have      f (r )   W ( r ′ − r , h)       ∂ f (r )  −1  ∂W ( r ′ − r , h)   ∂x  = ∫W f (r ′) N   dr′ ∂x      ∂ f (r )   ∂W (r ′ − r , h)   ∂y    ∂y

(13)

where   ∫ W (r ′ − r , h)d r ′ ∫ ( x′ − x)W (r ′ − r , h)d r ′  WW ∂W (r ′ − r , h)  ∂W (r ′ − r , h) ′ N = ∫ d r ∫ ( x′ − x) dr′ WW ∂x ∂x  ∂W (r ′ − r , h)  ∂W (r ′ − r , h) d r ′ ( x′ − x ) dr′ ∫  ∫WW ∂y ∂y    ∂W (r ′ − r , h)  d r ′ ∫W ( y′ − y) ∂x  ∂W (r ′ − r , h)  dr′ ∫W ( y′ − y)  ∂y

FPM allows a second order accuracy for both f (r ) and its gradients for internal and boundary particles alike. In case only derivatives of the f (r ) are needed, the differences of accuracy and stability between Eq.(10) and Eq.(13) still need further investigations. Although here the formulations are merely derived in 2-D, they can be extended straightforwardly to 3-D. (3) The kernel gradient free SPH (KGF-SPH) method derived based on CSPH and FPM Recently, another kernel interpolation formulation has been developed[94]. In that method, the kernel derivatives are never used, which helps to make the kernel interpolation free from the choice of kernel functions since for some kernels, the derivatives are not as smooth as the kernel functions themselves. This technique was initially named symmetric smoothed particle hydrodynamics (SSPH) method, based on which Jiang et al.[60] developed a mixed corrected SSPH (MC-SSPH) and Huang et al.[91] developed a kernel gradient free (KGF) SPH method. Recently, Xu et al.[95] extended this technique to 3-D dam breaking problems. Hereinafter, as long as the gradient of the kernel function is not used, similar to Huang et al.[91], this group of SPH variants is simply called KGF-SPH. Instead of multiplying both side of Eq.(5) with ∂W (r ′ − r , h) / ∂x and ∂W (r ′ − r , h) / ∂y , in the KGFSPH, ( x′ − x)W (r ′ − r , h) and ( y ′ − y )W (r ′ − r , h) is adopted. Consequently, Eq.(8) and Eq.(9) is updated by f (r ′)( x′ − x)W (r ′ − r , h)d r ′ = f (r ) ∫ ( x′ − x) ⋅ ∫WW W (r ′ − r , h)d r ′ +

∂ f (r ) ( y ′ − y )( x′ − x)W (r ′ − r , h)d r ′ + O (h 2 ) ∂y ∫W (16)

∫W ( y′ − y)W (r ′ − r , h)d r ′

(14) and

and therefore we have the modified kernel function and its derivatives as follows      W C 2 ( r ′ − r , h)   W ( r ′ − r , h)       ∂W C 2 (r ′ − r , h)  −1  ∂W ( r ′ − r , h)   =N   ∂x ∂x     C2  ∂W (r ′ − r , h)   ∂W (r ′ − r , h)    ∂y ∂y    

∂ f (r ) ( x′ − x) 2 W (r ′ − r , h)d r ′ + ∂x ∫W

f (r ′)( y ′ − y )W (r ′ − r , h)d r ′ = f (r ) ∫ ( y ′ − y ) ⋅ ∫WW W (r ′ − r , h)d r ′ +

r , h)d r ′ +

(15)

∂ f (r ) ( x′ − x)( y ′ − y )W (r ′ − ∂x ∫W

∂ f (r ) ( y ′ − y ) 2 W (r ′ − r , h)d r ′ + O (h 2 ) ∂y ∫W (17)

Solving Eq.(16) and Eq.(17) together, we obtain the kernel gradient free CSPH equations

192

 ∂ f (r )   ∂x    = ∫ [ f (r ′) − f (r )]( M ′) −1 ⋅  ∂ f (r )  Ω  ∂y   ( x′ − x)W (r ′ − r , h)  ( y ′ − y )W (r ′ − r , h)  d r ′  

  ∫W ( y′ − y)( x′ − x)W (r ′ − r , h)d r ′  2 ∫W ( y′ − y) W (r ′ − r , h)d r ′ 

∫W ( y′ − y)W (r ′ − r , h)d r ′

(18)

 ( x′ − x)2 W (r ′ − r , h)d r ′ ∫ W  ′ M =  ( x′ − x)( y ′ − y )W (r ′ − r , h)d r ′  ∫W

2

and therefore we have the kernel function and its derivatives allowing a second order accuracy as follows    W C 2′ ( x ′ − x , h )     W ( x ′ − x , h)   ∂W C 2′ ( x ′ − x , h)   −1  ′ ′ ′ = ( ) ( x x ) W ( , h ) − − N x x     x ∂   ( y ′ − y )W ( x ′ − x , h)  C 2′  ∂W ( x ′ − x , h)    ∂y  

where

∫W ( y′ − y)( x′ − x)W (r ′ − r , h)d r ′ ∫W ( y′ − y) W (r ′ − r , h)d r ′ 

(22)

(19)

and therefore we have kernel derivatives allowing a second order accuracy for the generic gradient approximations as follows  ∂W C ′ (r ′ − r , h)    ( x′ − x)W (r ′ − r , h)  ∂x   = ( M ′)−1   (20) ′ C  ∂W (r ′ − r , h)  ( y ′ − y )W (r ′ − r , h)    ∂y  

And similarly if we solve Eq.(6), Eq.(16) and Eq.(17) together, we obtain the kernel gradient free FPM as follows    f (r )     W ( r ′ − r , h)   ∂ f (r )   −1   ∂x  = ∫W f (r ′)( N ′)  ( x′ − x)W (r ′ − r , h)  d r ′   ( y ′ − y )W (r ′ − r , h)   ∂ f (r )   ∂y 

(21)

(23)

KGF-SPH is suitable for a wide range of kernel functions[60,61]. Besides, as stated in Jiang et al.[60], the corrective matrices M ′ and N ′ are symmetric, which helps to improve the numerical stability. Moreover, the symmetric matrices help to save storage requirement as stated in Xu et al.[96]. With the kernel interpolation of second order accuracy, it is straightforward to extend these continuous equations to discretized particle approximations, as been done in Eq.(4). Due to the limitation of the paper length, for more details the readers can refer to Refs.[60,61,95,96]. In spite of the higher accuracy of the equations aforementioned, the governing equations are still discretized by classic SPH interpolating formulations in most of the literature. The reason why the classic SPH formulations are mostly adopted is their superior properties of conservation and stability. The traditional kernel and its derivatives are symmetric between the particle pair, therefore the mass, momentum are strictly conserved, for example, in standard SPH and δ SPH. The matrices used to correct the kernels can be singular when the kernel is truncated severely, which is harmful to the numerical stability.

where

1.2 Governing equations

 ∫W W (r ′ − r , h)d r ′  N ′ =  ∫ ( x′ − x)W (r ′ − r , h)d r ′  W  ( y ′ − y )W (r ′ − r , h)d r ′  ∫W

1.2.1 Governing equations for the fluid evolution (1) Classic governing equations: Form 1 Here we recall a set of governing equations which are applicable for both free surface and multiphase flows as follows:

∫W ( x′ − x)W (r ′ − r , h)d r ′ ∫W ( x′ − x) W (r ′ − r , h)d r ′

Dρi = ρi ∑ (ui − u j ) ⋅ ∇iWijV j Dt j

2

∫W

( x′ − x)( y ′ − y )W (r ′ − r , h)d r ′

Dui 1 =− ρi Dt

∑( p j

i

+ p j )∇iWijV j +

(24a)

193

α hc0

ρ0 ∑π ∇ W V + g ρi j ij i ij j

Dri = ui dt pi = p ( ρi ) , πij =

(ui − u j ) ⋅ (ri − rj )

ρi =

(24c)

D ui 1  1 = − ρi  Vi Dt

(24d)

2

ri − rj + (0.01h) 2

Equations (24) is based on the kernel approximation and particle approximation introduced in Section 1.1.1. We underline that in the continuity equation, the divergence for the velocity is calculated based on the velocity difference ui − u j between the neighboring particles while in the momentum equation, the pressure gradient is based on the pressure summation pi + p j . The reason is due to the kernel truncation on the free surface, as stated in Section 1.1.1 and more details can be found in Colagrossi et al.[87]. The second term on the right hand-side of the momentum is an artificial viscosity, in which α is the coefficient to smooth the velocity field. c0 is the artificial speed of sound, determined by the weaklycompressible hypothesis as c0 ≥ 10 max(U max , pmax / ρ0 ) where U max and pmax are the expected maximum velocity and pressure. For a viscous problem, α hc0 can be replaced by 2(d + 2)ν where d is the spatial dimension and ν is the kinetic viscosity. For an inviscid flow, the determination of α is usually set to ensure a local Reynolds number Reh =

U max h /ν  O (10)

for stability reasons[13,37]. The

linearized equation of state written as pi = c 2 ( ρi − ρ0 ) can be adopted under the weakly compressible regime. ρ0 is the reference density when the particle is on the free surface with a zero pressure. Equations (24) can be used in either the free surface flows or multiphase flows, some applications can be seen in Colagrossi and Landrini[22]. In the present work, Eqs.(24) are mainly employed for simulations of free surface flows, e.g. wave propagations, wavebody interactions, etc.. For the multiphase flow, we recall another set of formulations in the following subsection. (2) Governing equations for multiphase flows: Form 2 The second form of governing equations shown in the following is established for the simulation of multiphase flows without a free surface. The formulations read as:

Vi ≈ ∑ j χ ijV j = ∑ j

Wij

σ (ri )

Vj =

∑ WV j

ij

σ (ri )

j



1

σ (ri )

mi Vi

(24b)

,

(25a)

∑

j

( piVi 2 + p jV j2 )∇iWij + FiV +

Fi S ] + g

(25b)

γl  cl 2 ρ0l  ρi   pli =  − 1 + pb γ l  ρ 0 l   

pgi =

cg 2 ρ 0 g  ρ i  g g  ρ0 g 

FiV = ∑ j

2hih j (hi + h j )

(25c)

gg    − 1 + pb   

(Vi 2 + V j2 )

(25d)

(ri − rj ) ⋅ ∇iWij (ui − u j ) 2

ri − rj + (ε h) 2

Vi

(25e) Fi

S

∑ (nˆ − j nˆ ) ⋅ ∇ W V = βd ∑ r −r ⋅ ∇W V i

j

j

i

i

j

j

j

i

i

ij

ij

j

∇ ci

(25f)

j

In Eqs.(25), a particle interpolation for the volume is adopted and the density is updated using ρi = mi / Vi , therefore the density approximation is free from the restriction of the large density difference near the multiphase interface[83]. The momentum equation can be derived based on the principle of virtual work (PVW), as shown in Grenier et al.[84] and Zhang et al.[17]. The equation of state for the denser phase is denoted by the subscript l and that for the lighter phase is denoted by the subscript g where γ l = 7 and g g = 1.4 are usually adopted as proposed in Colagrossi et al.[22]. ρ0l and ρ0 g are the reference density when the pre-

ssure p is equal to pb which is the background pressure used to prevent the SPH specific numerical instability, namely tensile instability, and keep the particles distributing uniformly[82,85]. The viscous force FiV is calculated with the method proposed in Hu and Adams[83]. On the multiphase interface, the dynamic viscous coefficient is calculated using the harmonic mean value as ηij = S 2ηηηη can be doi j /( i + j ) . The surface tension Fi minant in some cases, e.g. the small scale bubbly flows. Here the surface tension formulation proposed in Adami et al.[86] is recommended. A variant of the normalized divergence operator similar to the CSPH introduced in Section 1.1.2 is applied for the curvature calculation of the water-air interface, but here the ma-

194

trix inverse is avoided thanks to the mathematical work done in Adami et al.[86]. In Fi S , β is the surface tension coefficient. nˆi is the unit normal vector calculated as nˆi = ∇ci / ∇ci , where the color function

cij is defined as follows: 2 ρi cij = if particle i and j belong to ρi + ρ j different phase

(26a)

cij = 0 if particle i and j belong to the same phase (26b) and ∇ci is calculated as = ∇ ci

1 Vi

∑

(Vi 2 + V j2 ) j

cii + cij ∇iWij 2

(27)

ji j was defined in Zhang et al.[17] to inverse the nor-

mal directions for different fluid phases on both sides of the interface as

ji j = −1 if particle i and j belong to different phase

(28a)

ji j = 1 if particle i and j belong to the same phase

(28b)

In Eqs.(25), a summation density is adopted, it leads to kernel truncation near the boundary. But if we arrange air particles above the fluid surface, it can still be applied for the simulation of flows with a water-gas surface[85], but that brings about more expensive computational costs. Regarding the the multi-phase model in the framework of ISPH, see more in the paper by Shao[41]. (3) Governing equations for explosion problems: Form 3 For explosion problems like underwater explosions in free field, underwater contact explosions, the continuity and momentum equations have the same form with Eqs.(24). The main differences lie in the equation of state which takes into account the inner energy. In addition, the variating smoothing length[97] due to the big compressibility of the fluid phase is also considered. The energy equation is written as Dei 1 = Dt 2 ρi

α hc0

∑( p

i

j

ρ0

2 ρi

+ p j )(ui − u j ) ⋅ ∇iWijV j −

∑π j

ij

(ui − u j ) ⋅ ∇iWijV j

(29)

During the expansion of the explosive gas, the smoothing length needs to be updated according to dh 1 h dρ =− dt d ρ dt

(30)

Different from the equation of state for free surface flows, the state equation for explosion problems should consider also the inner energy, thus real sound speed of the fluid material should be adopted. For the water modeling, Gruneisen form of the equation of state should be adopted[57]. The standard Jones-Wilkins-Lee (JWL) state equation is employed to model the products of explosion. If the material strength is considered for modeling the structure responses in contact explosions or high velocity impact (HIV) problems, the Mie-Gruneisen equation of state for solid mechanics[16,98] can be applied. The Johnson-Cook yield model can be employed to evaluate the dynamic yield stress and the von Mises yield criterion is used to determine whether the solid material yields or no[57,99]. (4) Governing equations of higher accuracy: Form 4 The main difference between form 4 and the forms 1-3 is the adoption of the kernel gradients. As reviewed in Section 1.1.2, the KGF-SPH[91] can be derived based on CSPH or FPM. The kernel gradient ∇Ci 2′Wij are used to replace the kernel gradients in Eqs.(24) to discretize the continuity equation and momentum equation. Many applications using these governing equations such as dam-breaking flows[95], 3-D transient generalized Newtonian free surface flows[61] are published recently. In summary, for governing equations, form 1 is widely used for various applications, particularly suitable for free surface flows in ocean engineering, form 2 is widely used for multiphase flows and form 3 is widely used to model HIV or underwater explosions. The form 4 allows of a higher accuracy in the simulations particular when the kernel truncation is truncated. In Jiang et al.[60], it is shown that with the MCSSPH model, the modeling of rotating square obtains a more uniform particle arrangement and for the deformation of the fluid patch, a smoother vorticity field can be achieved with respect to traditional SPH models. But since the symmetric characteristic, i.e. ∇Ci 2′Wij = −∇Cj 2′Wij , is lost, with respect to the classic kernel gradients, it is difficult for the kernel gradient ∇Ci 2′Wij to ensure a strict conservation of mass and momentum for a long term simulation with a disordered free surface deformation. In addition, the treatment for the singularity of the corrective matrix is still an open problem. Therefore, for engineering applications, the first three forms are mostly used at the present stage and form 4 still needs further investigations

195

in different benchmark cases[50]. 1.2.2 Governing equations for the rigid body motions The governing equations for the motions of six degrees of freedom (6-DOF) of the solid body are expressed as follows[100,101]:

M U x = Fxfluid −solid + g , M U y = Fyfluid −solid + g , M U z = Fzfluid −solid + g

(31)

I x Ωx − ( I y − I z )Ω y Ωz = Txfluid −solid

(32a)

I y Ω y − ( I z − I x )Ωz Ωx = Tyfluid −solid

(32b)

I z Ωz − ( I x − I y )Ωx Ω y = Tzfluid −solid

(32c)

rG = U , θ = Ω

(33)

Ω = (Ωx , Ωy , Ωz ) denote the transitional and angular velocity for the solid body in different axis components. θ denotes the vector of angle of rotation. M and I = ( I x , I y , I z ) refer to the total mass and moment of

where

U = (U x ,U y ,U z )

and

inertia with respect to the gravity center G . The function f denotes a time derivative of a generic variable f . Similar to Kajtar and Monaghan[102] and Bouscasse et al.[37], the approximations of the forces and moments on the solid body are derived based on the momentum balance between the fluid and solid body particles as follows: Ffluid-solid =

Tfluid-solid =

∑ ∑

i∈fluid j∈solid

( pi + p j -p α c0 h ρ0 ij ) ⋅ ∇iWij

mi m j

ρi ρ j (34)

mi m j  ρρ  i + j - ρG  ×  2 i∈fluid j∈solid ρi ρ j  

∑ ∑

[( pi + p j − α c0 h ρ0 pij )∇iWij ]

(35)

where rG is the position vector of the gravity center. For the coupling algorithm between the fluid and body particles can refer to Bouscasse et al.[37] when a fixed ghost particle boundary is implemented and to Guo et al.[101] when a dummy particle boundary is adopted. When an ISPH model is adopted, two different kinds of coupling algorithms can be found in Shao and Gotoh[103] and Liu et al.[38].

1.3 Techniques for preventing high frequency pressure fluctuations In WCSPH, since the fluid is assumed to be weakly compressible and an equation of state is adopted to link pressure and density, the pressure field is very sensitive to the density oscillations, which leads to a main challenging issue in the WCSPH model. In the last decades, several remedies have been proposed to smooth the pressure field. Recently, a comprehensive study is also carried out in Cercos-Pita et al.[104] to study the diffusive terms. In the present section, some widely used numerical techniques are reviewed. 1.3.1 The density filtering technique The density filtering is an efficient way to remove high frequency density oscillations. Generally, there are two widely adopted techniques respectively with the first order and second order accuracy. The former uses the Shepard kernel interpolation[93] and the latter uses the well-known MLS interpolation[22,29]. Moreover, generally for multiphase flows, the previous mentioned techniques are not applicable. However, recently a novel idea was proposed by Chen et al.[78], in which the pressure filter technique based on the Shepard kernel was modified as

ρinew =

∑ W m ∑ WV j

ij

j

j

ij

j

 j = ρ jV j , ρ j = , m

p j − pb ci2

+ ρi 0

(36)

where ρ j is the density recalculated based on the pressure of particle j but the speed of sound of particle i . The novel idea overcomes the limitation of the drastic density difference across the multiphase interface. 1.3.2 The Rusanov flux Although the density filtering technique supplies a useful way to smooth the pressure field, it has two drawbacks. The first one is that it cannot satisfy the volume conservation since the hydrostatic component can be improperly filtered[105]. The second drawback is its additional computational cost which requests another loop throughout all the particles. Considering this point, in many works, the density filtering is conducted every 20 time steps[22]. In this subsection, the Rusanov flux is recalled to smooth the density field which benefits also the pressure distribution. One merit is its value can be calculated along with the density derivatives avoiding the apparent increase of computational cost. Vila[106] and Lanson[107] used the Riemann solvers to evaluate a numerical flux between the neighboring particles, but it leads to superabundant numerical diffusion. Ferrari et al.[108] proposed to consider diffusion based on Rusanov flux only in the continuity

196

equation, which significantly reduces the numerical viscosity while still allows a smooth pressure distribution. The formulation is written as Dri = ri ∑ (ui − u j ) ⋅ ∇iWijV j + DiFerrari Dt j Ferrari i

D

= c0 ∑ ( rr j − i) j

rj − ri rj − ri

⋅ ∇iWijV j

(37a)

(37b)

The Rusanov flux is independent of any artificial parameters, which therefore makes the scheme convenient to be applied in different applications. However, in long term simulations, similar to the first order density filtering technique, the Rusanov flux cannot preserve the hydrostatic equilibrium either[35,36]. The main reason is the inconsistency due to the singularity of the density approximation close to the free surface[104]. 1.3.3 The δ - term To maintain the hydrostatic equilibrium in long term simulations, the well-known δ - SPH was proposed[35,36]. When the δ - term is taken into account, the continuity equation becomes Dρi = ρi ∑ (ui − u j ) ⋅ ∇iWijV j + Diδ Dt j Diδ = δ hc0 ∑ψ ij j

(rj − ri ) ⋅ ∇iWij rj − ri

ψ ij = 2( ρ j − ρi ) − ( ∇ρ

L j

2

(38a)

Vj

(38b)

L

(38c)

+ ∇ρ i ) ⋅ (ρρ j − i)

L

where ∇ρ is a renormalized gradient of the density using the corrected kernel gradients in Section 1.1.2. With the second term on right-hand side in ψ ij , the singularity caused by the kernel truncation on the free surface is avoided. δ is not a free altered parameter which is usually set at 0.1. The δ - term is an improved version of the so-called simple procedure to stabilize the pressure as proposed in Molteni and Colagrossi[105]. 1.3.4 Techniques adopted in ISPH or MPS models The pressure fluctuation is also a major issue in ISPH or MPS models. The error-compensating source term[109] and higher-order Laplacian formulation[110] developed by Khayyer and Gotoh in Kyoto University are widely used. Besides, the hybrid PPE source term[39,40] is also effective in reducing the pressure noise.

1.4 Boundary conditions In SPH, many efforts have also been devoted to the implementation of boundary conditions. Generally, the boundary conditions in hydrodynamics can be mainly divided into two categories: free surface boundary and solid wall boundary. 1.4.1 Free surface boundary For the free surface flows, two conditions (kinematic and dynamic conditions) have to be satisfied on the free surface[87,111]: (1) The fluid particle initially on the free surface should be always on the free surface, except for the water wave rolling and impact. (2) The pressure on the free surface is equal to the background/atmosphere pressure. The two above mentioned conditions can be implicitly satisfied in the framework of WCSPH if proper governing equations are adopted as discussed in Colagrossi et al.[87]. However, in ISPH or MPS, it is another story in which the free surface should be precisely detected and a free surface condition (dynamic condition) should be imposed. At present, several good algorithms have been proposed to detect the free surface particles[28,44,112-115]. 1.4.2 Solid wall boundary Regarding the solid wall boundary conditions, they can be divided into two categories: free slip boundary or no-slip boundary. In both the two different boundary conditions, a non-penetration condition has to be satisfied. To achieve this goal, usually two manners are adopted. The first one is to use a surface integral on the boundary to satisfy the boundary condition and the second one is to distribute ghost particles on and inside the boundary surface and extrapolate the generic variables from the fluid particles to these ghost particles. In the first category of boundary treatment, the well-known techniques are normal flux method[116] and the so called unified semi-analytical wall boundary condition[117]. It is worth noting that the normal flux method has been successfully applied to discretize the irregular ship hull in the simulation of bow breaking waves in 3-D[30]. The well-known repulsive force boundary can also be classified into this category since the repulsive force works in a similar way to the boundary integral term[118] but it lacks consistency. In the second category where ghost particles are distributed, there are plenty of different techniques published. The simplest one is the mirroring ghost boundary as proposed in Cummins and Rudman[119] and Colagrossi and Landrini[22], however it is only restricted to the applications with regular shaped boundary shapes. In order to improve the mirroring ghost boundary, fixed ghost boundary was proposed in Marrone et al.[92]. A MLS interpolation is used to interpolate the fluid properties from the fluid particles

197

to the interpolation point, which matches the fixed ghost particles one by one. The procedure for the implementation of the fixed ghost boundary is relatively tough, especially in 3-D cases. Therefore, procedures with a more convenient way in interpolating the fluid properties to the boundary are proposed. In Adami et al.[120], a generalized wall boundary condition was proposed and has been applied in many ocean engineering applications, see examples in Refs.[15]. In Liu et al.[121], beside the dummy particles (or ghost particles) inside the boundary surface, a layer of boundary particle is also distributed along the boundary surface to prevent the particle penetrations. This method was named the coupled dynamic SBT algorithm. This technique was successfully applied in problems of water entry and exit[122] and sloshing[90]. Recently, it is shown in Cercos-Pita et al.[123] the presence of some extra terms from the ghost particles does not affect the consistency of the SPH model, which further proves that using ghost particles to model the wall boundary is reliable. The difference between free-slip boundary and no-slip boundary lies in the implementation of viscous force. In order to satisfy the free-slip condition, the viscous force is imposed to be zero for a particle pair in which the two particles are respectively fluid and boundary particles[120]. While in the implementation of a no-slip boundary, the viscous force evaluation between fluid and boundary should be carefully addressed. Up to now, different techniques have been proposed concerning the accuracy, consistency, stability and operability for regular or irregular boundary shapes in 2-D or 3-D. Maciá et al.[124] carried out a theoretical analysis of the no-slip boundary condition enforcement in SPH methods. Four different extensions of the velocity from fluid to boundary particles were investigated, namely, dummy particles fixed boundary velocity (U0M)[125], ghost particles with antisymmetric mirroring (ASM)[120], the Takeda et al.[126] imaginary particles and symmetric mirroring technique (SSM)[13]. They proposed a renormalized Takeda extension permitting a consistent redefinition such that all the operators converge to the correct values. However, as was pointed out in Marrone et al.[13] the local inconsistency does not jeopardize the consistency of global quantities. This means that in practical applications, ASM and SSM are also applicable when the particle resolution is fine enough, see also the applications in Sun et al.[50]. 1.4.3 Inflow and outflow boundary Inflow and outflow boundary conditions are essential for the applications in ocean engineering. A representative example is the study of flow past solid bodies. Free stream flow past the bluff body leads to a periodical shedding of vortices which may induce vibrations being harmful to the ocean structures. In the

WCSPH framework, the inflow-outflow boundary is relatively convenient to be implemented. The technique proposed in Federico et al.[127] is the most widely adopted, see the applications in Refs.[13,24,50]. Recently in Kazemi et al.[128], another algorithm is proposed. In the SPH model where semi-analytical conditions and Riemann solver are adopted, an unsteady open boundary was proposed in Ferrand et al.[129]. Recently, in the framework of ISPH, different open boundaries are also proposed in Refs.[130-132]. It is worth noting that when the Lagrangian Coherent Structure is detected inside the flow, some modifications should be made on the inflow boundary in order to evaluate the FTLE, see more details in Sun et al.[24]. 1.4.4 Non-reflection boundary (1) Sponge layer As described in Lastiwka et al.[133], the pressure wave due to the fluid compressibility reflecting on the boundary and propagating back to the interior fluid domain is a non-physical behavior since they should propagate freely across the boundary. A sponge layer can be implemented near the stationary solid boundary wall with a certain thickness s to prevent the nonphysical pressure wave reflection. This sponge layer was originally used in Gong et al.[12] and then applied by Sun et al.[15] in water entry problems. For the fluid particles inside the sponge layer, the time derivatives of the density and velocity should be multiplied with a term as 50 λ  Dρi   Dρ  = (1 − 100−0.9 )  i   Dt   sponge  Dt 

(39a)

50 λ  Dui   Du  = (1 − 100−0.9 )  i   Dt   sponge  Dt 

(39b)

where λ is defined as λ = ( s − d ) / s , d is the distance between the fluid particles i and the solid boundary wall. (2) Improving viscosity near the wall boundary The sponge layer works well for preventing spurious pressure wave reflection, however in the simulation of the propagation of gravity water waves, the reflection of the wave can also be damped by improving the viscosity near the wall boundary. This technique has been very well adopted in Wen et al.[134]. (3)Non-reflecting Boundary condition When an inflow and outflow boundary is implemented, the sponge layer introduced previously is not applicable. Therefore in this circumstance, non-reflecting boundary condition introduced in Lastiwka et al.[133] is recommended. 1.4.5 Particle packing technique For the fluid-solid interactions with irregular

198

boundary shapes, the arrangement of the initial particle distribution close to the solid wall boundary cannot be trivial. If the particles are not very well distributed initially, spurious pressure waves can be generated due to the particle resettlements. Thanks to the particle packing technique proposed in Colagrossi et al.[135], a uniform particle distribution around the solid wall boundary can be obtained. In Fig.1, the particle distribution around a submarine shape is shown after particle packing algorithm was used. The zoom view on the right-hand side of the same figure demonstrates a uniform particle distribution around the model stern with a complex configuration.

Fig.1 The particle arrangement surrounding the submarine shape after the using of the particle packing algorithm[135]

1.5 Adaptive-particle resolution In order to well reproduce the physical phenomena in ocean engineering applications, a sufficiently fine particle resolution should be imposed. For example, in the numerical simulation of wave propagations, it was shown in Antuono et al.[136] that particle number of the order of 10 should be distributed in the wave height to avoid overdamping and achieve a convergent solution. Another example is in the simulation of flow past a bluff body. In order to obtain a convergent evaluation of the force on the body, in particular for reproducing the viscous force in the boundary layer, a particle resolution of D / Dx ≥ 80 should be distributed with D being the “projected frontal side” of the rigid body, see more in Marrone et al.[13]. However in these cases, the fluid domain is quite large in order to mimic the free stream and therefore the adoption of single particle resolution leads to very expensive computational costs. In the mesh-based methods, e.g. FDM[137], the adaptive mesh scales are usually adopted, which improves the accuracy in local areas but reduces the total computational cost due to the reduction of grid elements. In SPH, due to the Lagrangian trajectory of the particle, to achieve an adaptive-particle resolution is not as straightforward as what done in mesh-based methods. Despite the difficulty, special treatments are developed in recent years and will be reviewed in the following subsections. 1.5.1 Particle splitting and coalescing One technique is to split one big particle (namely

mother particle) into several smaller particles (namely daughter particles) to achieve a finer particle resolution in local regions. At the early stage of development of this technique, SPH with variable smoothing length was often applied, in which, the particle resolution is gradually increased from the region of interest up to the borders of the computational domain[138]. In a similar way, Koukouvinis et al.[139] refined the particles in certain regions at the initial stage, and during the simulation, no more particles are refined. Recently, Tang et al.[140] applied this technique to water entry problems in the context of MPS. The methods mentioned above are easy to be implemented since the total particle number is fixed. Starting from Feldman et al.[141], a dynamic particle splitting method has been proposed, which is very attractive since only the particles (mother particle) which enter the region of interest will be split. That means the total particle number during the simulation is always updating. It increases the complexity of the algorithm. The most challenging problem lies in the treatment of accuracy and stability on the splitting interface where the smoothing lengths and particle arrangements have a sharp transition. The error is related to the adopting of particle splitting configurations and the assignment of the smoothing lengths of the daughter particles. Feldman et al.[141] proposed two splitting patterns of hexagon and triangle, which allow satisfactory accuracy at the splitting interface. Omidvar et al.[142,143] simplified Feldman’s method with a square particle splitting pattern in two dimensions and a cubic pattern in three dimensions. An systematic investigation was carried out in López et al.[144]. An improved splitting error analysis was proposed in their work, in which the error analysis was based on the derivative of f (r ) ( f (r ) is a generic physical variable). However, since the smoothing lengths of the daughter particles are constrained in order to avoid particle clumping, a non-uniform distribution of daughter particles has to be adopted in order to decrease the splitting error. This may create problems when a particle shifting technique is adopted, see e.g. Sun et al.[50], since the latter always tends to make the particle distribution more uniform. After the daughter particles are transported out of the splitting domain, a particle merging technique should be adopted as the one proposed in Vacondio et al.[145,146]. 1.5.2 Particle overlapping technique In the particle overlapping technique, the mother particle split in a certain region is still kept but switched off and turned on when they flow out of that region. Therefore, the particle coalescing, which is complicated for programming, is avoided. Sketches of the particle splitting and overlapping are depicted in Fig.2. Figure 2(a) shows the mother particle (Level 0,

199

blue), the daughter particles after the first splitting (Level 1, green) and the daughter particles after the second splitting (Level 2, red). The splitting interface 1 and 2 means the mother particles across the interface will be split into four particles and the corresponding daughter particles flowing out of the splitting interface will be erased. Figures 2(b) and 2(c) show the enlarged views of the two square boxes by dashed lines in Fig.2(a). Figure 2(b) shows the particle distribution after the first splitting with two scales of particles overlapping. Figure 2(c) shows the particle distribution after the second splitting with three scales of particle overlapping.

Fig.2 Sketch of the particle distribution in the particle overlapping technique

For the time being, there are two different overlapping techniques developed respectively in the frameworks of WCSPH and MPS. The overlapping technique in WCSPH was first proposed in Barcarolo et al.[49] and then further extended and generalized in Sun et al.[50]. The movements of the mother particles

are based on the interpolation of the velocity from the daughter particles. While in MPS, the movements of the mother particles are obtained by solving the Poisson equation based only on the mother particles. And then the movements of the daughter particles are solved by the Poisson equation based on the quantities on the splitting boundary where the quantities are interpolated from the mother particles, see more details in Tang et al.[115]. 1.6 Preventing of tensile instability SPH method suffers from tensile instability when the pressure becomes negative. The behavior can be described as the particles tend to stay in pairs or clustered. As the tensile instability develops, numerical fluid voids caused by the fluid material fragmentation are observed in local regions. The most common tensile instability phenomena appear in the simulation of lid driven cavity[15] and downstream the flow past a bluff body[13]. To prevent this non-physical behavior, Monaghan[147] proposed to add an artificial stress added into the momentum equation to prevent the tensile instability. However, in WCSPH, the pressure is linked to the density variation, which is directly related to the evolution of the particle velocity. When an artificial stress is added into the momentum, the velocity evolution is never consistent with the N-S equations, and therefore the pressure field becomes unstable and even worse than all, the numerical results are wrong. An example is the well-known benchmark test case called rotational fluid patch which has been widely adopted to test the consistency, accuracy, stability of a numerical scheme[50,60,148-150]. In Le Touz et al.[148], it was shown that after adopting the artificial stress, the phenomenon of particle clustering was removed, while another problem of pressure fluctuation in the fluid center was observed. This is due to the inconsistency brought about by the artificial stress. Another remedy proposed to remove the tensile instability is named the dynamic stabilization (DS) in the framework of MPS, see Tsuruta et al.[150]. This technique is suitable for an incompressible scheme (MPS or ISPH) while in WCSPH schemes it has the same drawback of introducing pressure noises. For the flow without a free surface, background pressure is the most suitable and convenient way to prevent tensile instability, see the simulations of the lid driven cavity[15] and flow past body[13]. As pointed out in Adami et al.[151], a drawback of background pressure is the inducing of the spurious numerical dissipation. Therefore in that work, a transport-velocity formulation for smoothed particle hydrodynamics was proposed. In another work of the same authors[152], the difference between background pressure and transport-velocity formulation in SPH was analyzed. However, these techniques cannot be applied to free surface problems. Therefore in Oger et al.[149], a con-

200

sistent ALE formalism was implemented in the framework of Riemann solver based SPH (one of the variants of WCSPH). That scheme is consistent and can be applied to violent fluid-solid interacting problems with a free surface. Recently in Sun et al.[50], an δ + SPH scheme was proposed based on the well-known δ - SPH by introducing a particle shifting technique, which was originally proposed by Xu et al.[45] and extended by Lind et al.[153] to free surface flows. In Sun et al.[50], a novel technique was introduced to treat the free surface in order to make δ + - SPH suitable for violent free surface simulations. Moreover, the shifting was also generalized to multi particle resolutions.

correction scheme and also contributes to a better numerical stability. In either the modified predictioncorrection scheme or modified Verlet scheme, the evaluation of the general variables is conducted at the midpoint time step n +1 / 2 except the updating of the particle position which is based on the velocity evaluated on the time step n +1 . It is proven that with this modification, a larger CFL factor ( CFL ≈ 1 ) can be used. In ISPH and MPS models, the popular two-step projection solution (prediction-correction) method is usually adopted, see more in Refs.[41,109,119,125, 154].

1.7 Time-stepping scheme

2. The recent applications of SPH in fluid-structure interactions in ocean engineering

1.7.1 Runge-Kutta integration scheme For the modeling of hydrodynamic problems in ocean engineering, the well-known δ - SPH scheme[92] and its further improved variant δ + - SPH[50] can be adopted. In these schemes, the most adopted explicit time integration may be the 4th order Runge-Kutta integration[22,36,92]. The merit of such a time integration scheme is the allowance of a relatively larger time step. The time step ∆t in the 4th order Runge-Kutta integration method is constrained by the maximum acceleration amax = max( ai ) (both the fluid and body

2.1 SPH simulations of hydrodynamic problems With the rapid development of SPH schemes and numerical techniques, SPH has been widely and successfully applied in the simulation of fluid-structure interactions in ocean engineering. Particularly in the WCSPH method, thanks to the benefits of its low CPU costs, ease of implementation and parallelization, large simulations can be launched on high-performance computational machines. In this section, different SPH simulations of hydrodynamic problems are presented.

i

particles), the maximum velocity umax = max( ui ) i

(both the fluid and body particles) and artificial sound speed c0 as

 h ∆t ≤ min CFL1 , amax    h    CFL 2 min i  c0 + umax + h max φij   j  

(40)

where the CFL coefficients are heuristically found to be: CFL1 ≤ 0.25 and CFL 2 ≤ 2.0 . 1.7.2 Prediction-correction scheme In the modeling of a multiphase flow, a modified prediction-correction time stepping scheme was presented in Zhang et al.[17]. In that work, the bubble rising and bubble interaction were simulated in 3-D. It was shown that the modified prediction-correction scheme has a better numerical stability. In Molteni and Colagrossi[105], a similar time-stepping scheme was studied, and was named modified Verlet scheme, which has a form similar to the modified prediction-

2.1.1 Simulation of liquid sloshing Sloshing is a process of violent fluid motions inside a tank[155]. It has a wide application in ocean engineering, particularly important for the research and analysis of liquefied natural gas (LNG) carriers, oil tanks and other large scale liquid containers. SPH has been successfully applied to sloshing simulations by many researchers. Sloshing in different conditions are modeled, validated and investigated. To successfully model this violent fluid-structure interacting process, a robust SPH scheme with a stable pressure approximation is crucial (e.g. using techniques for preventing high frequency pressure fluctuations as presented in Section 1.3). Colagrossi et al.[87] conducted a theoretical analysis on the free surface conditions in WCSPH, explained why the surface condition is implicitly satisfied in the framework of WCSPH and a sloshing test case was also tested and validated in that paper by using the results of a BEM solver. Delorme et al.[156] carried out an experimental study to validate SPH results in the sloshing modeling. The influences of viscosity and of density re-initialization on the SPH results were discussed. After that, the proposal for δ - SPH further improved the SPH scheme allowing much more stable pressure approximations[92], removing the singularity on the free surface[36] and

201

giving a fast dissipation of the high-frequency energy in the water wave impact[157]. A good solid wall boundary is also essential. Chen et al.[158] carried out an investigation on the pressure on solid walls using SPH method. An improved boundary technique was also proposed in that work. Fixed ghost particle boundary works well in 2-D simulations and was adopted in Bouscasse et al.[155] to investigate shallow water sloshing. A 3-D simulation was made by Luo et al.[159] recently. In Cao et al.[10], sloshing with a baffle was investigated using a dummy particle boundary. Recently, in Hwang et al.[68], sloshing flows interacting with elastic baffles were modeled and investigated. Considering the coupled sloshing seakeeping effects, a SPH-FEM coupling technique was recently developed in Serván-Camas et al.[63]. After plenty of validations, SPH has been proven to be a reliable numerical tool to analyze sloshing flows in different conditions. 2.1.2 Water entry and exit of a circular cylinder (1) The water entry of a circular cylinder Water entry of a rigid body is an important phenomenon in ocean engineering operations. From releasing lifeboats, freak wave impact on ship structures to the launching of missile projectiles, the water entry can be quite often encountered. In the study of these problems, the circular cylinder water entries are the most representative ones since it is the shape mostly adopted as structural components in the construction of offshore platforms. 2-D simulations and investigations have been carried out thoroughly, but 3-D studies have mostly been conducted in experiments[160,161]. Recently, 3-D water entry problems are very well studied in Nguyen et al.[162] by using a Navier-Stokes solver. In the present subsection, 3-D water entry of a circular cylinder with the length of 0.20 m and diameter of 0.05 m is simulated. The vertical velocity when the cylinder just touches the undisturbed water surface is 6.22 m/s. The cylinder mass is 0.533 kg. The adopted SPH model is the well-known δ - SPH model. A free-slip wall boundary condition is applied on the cylinder surface and on the wall of the outside water tank whose geometrical scale is with the length of 0.625 m, width of 0.375 m and depth of 0.45 m. An experimental study was conducted for this case in Wei and Hu[161]. The numerical results are shown in Fig.3. Similar cavity shapes can be observed between the experimental and SPH results. Only half of the fluid particles are shown in order to more clearly demonstrate the cavity shape induced by the cylinder motion. The contour on the fluid particles shows the velocity variation while the contour on the cylinder surface shows the pressure distribution. In Fig.4, a good agreement is also obtained for the vertical displacements.

Fig.3 The snapshots of the water entry of a cylinder at t = 142 ms from different views

Fig.4 Time evolution of the vertical displacement Table 1 Parameters for the water exit tests Test number

No. 1

No. 2

Froude number Fr

0.4627

0.7644

Reynolds number Re

1 000

1 000

(2) The water exit of a circular cylinder Water exit, compared to water entry, has been much less investigated[163]. Beside the challenge from the dramatic free surface deformation, the water exit is also tightly related to the viscous effect which becomes particularly important when viscous flow separation occurs. While in the water entry, the viscous force can be less important since the non-viscous flow separation at the initial slamming stage is mandatory. Therefore, many works addressing slamming load during the water entry by BEM are published, see Refs.[164-167]. The water exit problems are mainly tackled in experiment or CFD solvers considering the viscous stress, see Refs.[122,168-170]. Recently, Ni et

202

Fig.5 The process for the water exit of No. 2 test

al.[171] simulated a complete water exit of a fully-submerged body through an enhanced BEM, mainly focus on the free surface deformation: free surface breaking and detachment. Wang et al.[172] adopted a viscous modification based on the results of BEM to evaluate the load on a 2-D circular cylinder during the water exit process. Taking into account the viscous effect, in this subsection, 2-D water exit of the cylinder is investigated with the recently developed δ + - SPH scheme. The numerical results are validated by the experimental data[173] and the other numerical results[163,172]. The water exit tests are characterized mainly by the Froude number Fr = U / gD and the Reynolds

number Re = UD / ν , where U is the velocity of the cylinder, g being gravity acceleration, D being cylinder diameter and ν being kinetic viscosity of the fluid. The parameters are shown in Table 1 corresponding to the experimental study[173]. Here the Reynolds number we use is 1 000, which on the one hand allows a satisfactory numerical stability for the SPH simulation, on the other hand considers the viscous effect on the flow separation. According to the experiment, the cylinder starts at Vt / R = −5.5 and after a short time acceleration, the cylinder reaches the expected Froude number at the time of Vt / R = −5 . The processes of the water exit for the test No. 2 at four time instants are plotted in Fig.5. A good

203

agreement for the cylinder position and the free surface deformation is achieved between δ + - SPH[50] and CIP[163]. Through the SPH results, we can see clearly how the vortices are detached from the cylinder at different time instants. In the cores of these vortices, the pressure is negative. The movements of the vortex affect the total drag force acting on the cylinder. Additionally, on the free surface, the pressure is always zero, as the cylinder gradually approaching the free surface, the total drag force which is equal to the pressure difference between the cylinder top and bottom is decreasing.

separation which is not taken into account in BEM. Besides, the results of BEM at Vt / D > 1 are not presented due to free surface breaking which terminates the simulation. Therefore, we can summarize that for the water exit problem, SPH can be a numerical solver more appropriate.

Fig.6 The time evolution of the coefficient of the vertical drag force acting on the cylinder during the water exit

In order to give a more clear comparison, time evolutions of the force coefficients for test No. 1 and No. 2 are plotted in Fig.6. The vertical force coefficient is evaluated by ce = Fv /(0.5 ρ V 2 D) . The SPH results are compared with some reference solutions in the literature. A fair agreement is obtained between the results of SPH and CIP since the viscous forces are both considered in the two solvers. The experimental results are accompanied with some fluctuations due to the automatic control of the hydraulic system used to enforce the test cylinder to exit the water[163]. However the overall trends of the force evolutions of the experiment and numerical results are in accordance. It is worth noting that the results of BEM is a little larger during the time range of Vt / D = (−5,1) . This is owing to the viscous flow

Fig.7 Time evolutions of the wave elevations at different horizontal positions. SPH results are validated by the experimental data[175]

Fig.8 Time evolution of the wave elevation at x / H = 7.6 with three different particle resolutions. SPH results are compared with the experimental data[175]

204

Fig.9 The snapshots of the interaction between the floating body and the freak wave at different time instants

2.1.3 Wave-body interaction Wave-body interaction plays an important role in ocean engineering. In rough seas, the wave induced motion, wave impact and green water loading can lead to a disaster for the offshore structure. Recently, SPH has been widely applied to model these kinds of interactions[15,37,103,174], particularly for the cases involving violent freak waves with structure motions presenting large amplitudes. To model such a complex interacting process, three aspects should be considered carefully. (1) The water wave generation using a wave maker. (2) The wave damping at the end of the numerical water wave tank. (3) The coupling between the fluid and floating body considering 6-DOF motions[162]. In the water wave generation, generally two types of wave makers

can be adopted: the piston wave maker or the flap one. The former is more suitable for the shallow water wave generation while the latter is better for relatively deeper water waves. In Antuono et al.[136], different water waves were studied, including standing waves, regular waves and wave packages using piston or flap wave makers, and the SPH results were validated by the numerical results through a BEM-MEL solver or the experimental data. The results demonstrated that WCSPH with the diffusive terms developed in Antuono et al.[36] is robust in modeling different water waves in particular for nonlinear waves with wave rolling and breakings. Regarding the wave damping, the non-reflection boundary introduced in Section 1.4.4 can be employed. Liu et al.[154] also introduced a non-reflection internal wave maker algorithm in the

205

framework of ISPH. The advantage of this algorithm is the wave generations are not affected by the reflected waves. Regarding the 6-DOF coupling motions between the fluid and floating body, in Bouscasse et al.[37], a technique was introduced when the fixed ghost particle boundary[92] was implemented. In Sun et al.[15], a coupling method between the floating body and the fluid based on the dummy particle boundary was presented. The governing equations introduced in Section 1.2.2 can be applied for the updating of the rigid body positions. As an example, the interaction between a 2-D floating body and a freak wave introduced in Zhao et al.[175,176] is simulated here by using a coupling technique based on the dummy particle boundary. Previously, this case was studied experimentally and numerically by a CIP method[175,176]. For more details regarding the case set-up, see Zhao et al.[176]. In the following, the SPH results are presented comparing with the experimental data.

Fig.10 Time evolutions of the motions of the floating body are plotted in the first three subplots and the wave elevation at x / H = 12.7 is plotted in the last subplot

Firstly, the freak wave generation by using a piston wave maker is validated. The time evolutions of the wave elevations at x / H = 7.6 , x / H = 12.7 , x / H = 17.5 , x / H = 22.2 and x / H = 27.4 are plotted in Fig.7. The adopted particle resolution is

H / dx = 80 where H is the undisturbed water height and dx is the initial particle spacing for discretizing the fluid domain. A good agreement between the SPH results and the experimental data is obtained in Fig.7. A convergence study for the wave elevation is also presented in Fig.8. As the particle resolution is increased, the wave elevations converge to the experimental data[177]. With the highest particle resolution of H / dx = 80 , the interaction between the freak wave and a floating body is simulated. The floating body can freely heave and pitch but the surging motion is prohibited. The snapshots of the motions of the floating body are depicted in Fig.9. The time history of the body motions and the wave elevation at x / H = 12.7 are plotted in Fig.10. A fair agreement is obtained between SPH and the experiment, which demonstrates the robustness of the present SPH scheme in modeling wave-body interactions. We underline that even though a 2-D simulation is conducted here, the SPH scheme containing a dummy particle boundary is very straightforward to be extended in 3-D cases.

Fig.11 Time evolution of the inflooding water through the opening of a floating structure

2.1.4 Simulation of the water flooding and sinking of a damaged floating structure During the lifetime of a surface ship, the crash, stranding or underwater explosion may make the ship lose its integrity. The water flooding through the

206

damaged hole of a ship and the water spreading through the openings inside the ship structure are big threats to the ship survivability[101]. Therefore a numerical model to predict the dynamic floating state of a damaged ship in different inflooding conditions is significant, which can provide a reference for the ship operations and sea rescues. To our knowledge, Idelsohn et al.[178] first proposed to model the sinking process of a damaged oil tank using particle methods. In Le Touzé et al.[179], SPH was proven to correctly simulate transient flooding behavior that occurs during and immediately after a side collision between two vessels. Recently, in Zhang et al.[29], numerical and experimental studies were carried out to study the motions of a damaged floating structure. A fairly good agreement was obtained between the results of SPH and experiment. Fig.11 shows the time evolution of the flooding process of a damaged floating structure. The structure has a length of 46 m, width of 16 m and height of 8.5 m, consisting of three water-tight cabins with the same scales. The water draft is d = 4.2 m . There is a circular opening with a diameter of 3.0 m located in the position which has the height of 1.6 m from the bottom and the distance of 2.5 m from the first bulkhead. Through the opening the water gradually floods into the cabin and a coupled motion between the fluid and structure occurs, see the snapshots at different time instants in Fig.11. 2.1.5 Flow past a circular cylinder by using adaptiveparticle resolutions Flow past a bluff body is a classic test case for CFD solvers. It also has important applications in the ocean engineering, for example, the applications in the preventing of vortex induced vibrations in flow past risers, offshore platforms and transmission pipelines. The flow features are dependent on Reynold numbers, Froude numbers and the boundary conditions, e.g. free stream condition (risers), flow past a body close to a free surface (offshore jacket platforms, transmission pipelines) and flow past a body close to seabed (subsea pipelines). SPH suffers from the so called tensile instability, which can induce a non-physical flow cavity downstream the body[13]. The newly developed δ + - SPH scheme which contains a particle shifting technique can help to avoid this numerical cavity. In this subsection, flow past a circular cylinder at the Reynold number Re = UD / ν = 1 000 is simulated firstly and then in the next subsection the cylinder is transported close to a free surface at a Froude number of Fr = U / gD = 1 to see the effect of the free surface on the vortex shedding. In the definition of the Froude number, U is the inflow velocity, D is the cylinder diameter, L is the distance between the cylinder top

and the undistributed free surface and g is the gravity acceleration. For the flow past a cylinder in free stream, the inflow boundary is imposed at x = −8D and outflow boundary at x = 24 D . The lateral boundaries are imposed at y / D = ±8D to avoid the blockage effect. Three levels of particle resolutions are adopted in this simulation. Close to the cylinder, within the radius of R2 = 1.75D , the particle spacing is D / Dx = 100 . Within the region of 1.75D < R1 < 2.25D , the particle resolution is D / Dx = 50 and in the rest of the fluid domain, D / Dx = 25 is adopted. After a short time acceleration, the inflow speed reaches the expected inflow velocity U . In order to quickly enter the vortex shedding stage, we use a free slip boundary for the upper half of the cylinder while a no-slip boundary for the lower half when tU / D < 5 . When tU / D ≥ 5 , a no-slip boundary condition is imposed on all the cylinder surface. The pressure and vorticity distributions at tU / D = 64.5 are depicted in Fig.12. A negative pressure region can be observed in the cores of the vortices. We underline that the background pressure[13] is not adopted here. Thanks to the shifting technique, the particle clustering and the non-physical flow cavitation can be avoided in the δ + - SPH when a small time step is used. The time evolution of the drag and lift force coefficients are plotted in Fig.13 along with the results by using the recently developed diffused vortex hydrodynamics (DVH) method[180]. As can be seen from the curves of force coefficient, fair agreement between the two solvers is observed.

Fig.12 The flow past a circular cylinder at Re = 1 000

2.1.6 Flow past a circular cylinder close to a free surface Differing from the flow past a body in free stream,

207

the flow past the body close to a free surface is far less investigated by hydrodynamic community. Some discussions of this topic can be found in Refs.[181,182] and recently in Refs.[183,184]. The SPH scheme validated in the previous subsection is employed here to simulate the flow past a circular cylinder close to a free surface. The Reynolds number adopted here is still Re = UD / ν = 1 000 . Owing to the benefits from the δ + - SPH scheme[50], the numerical cavitation due to the tensile instability behind the cylinder is avoided. The Froude number is set at Fr = U / gD = 1 as stated before. Initially, the cylinder center is placed at (0,0) . The inflow boundary is set at x = −6 D and outflow boundary at x = 18D . The fluid bottom is at y = −4.5D and the undisturbed free surface at y = −1.5D . It is worth underlining that to our knowledge, Re = UD / ν = 1 000 is the highest Reynolds numerical employed for the simulation of flow past cylinder positioning in the vicinity of a free surface.

Fig.13 Time evolutions of the drag and lift force coefficients on the circular cylinder at Re = 1 000 (Results of DVH from Ref.[180])

In Fig.14, the vortex structure downstream the cylinder is depicted by using the vorticity and FTLE[24] distributions. Comparing the vorticity field with the one in Fig.12, we find that the Von Karman vortex street is completely suppressed by the existing of a free surface. Beside the vortex shedding from the cylinder, it is also found that the wave breaking generates strong vortices which interact with the vortices from the cylinder. This is different from the results shown in Sun et al.[24] since here the Reynold number has been increased to 1 000 rather than 180. In such a high Reynolds number, dipoles composed of two large vortices of opposite signs are observed, see Fig.14.

Fig.14 The flow past a circular cylinder close to a free surface at t ( g / D )1/2 = 21.55 , Re = 1 000 and Fr = 1

2.1.7 Splashing bow wave of a high-speed ship model In this subsection, an example of the application of SPH to the breaking and splashing of bow waves of high speed ships are presented. Contrasted with the simulations in Sun et al.[24] and Marrone et al.[14], in order to make this case easier for the readers to reproduce, we adopt a simplified ship hull similar to the one invented in Sun et al.[185]. The water draft of the ship is one fifth of the ship height. In Marrone et al.[14] the breaking wave pattern of a fast ship has been studied in detail for different forward ship speeds. In that work a 2-D+t SPH model was used and the vortical structures due to plunging of the bow wave are tracked using passive markers which are initially positioned on the undisturbed free surface. The same technique has been also adopted in Landrini et al.[186] combining with vorticity contour plots. In this subsection, a case study is conducted to show the benefits of using the finite time Lyapunov exponent (FTLE) for capturing the bow breaking wave features[24,185]. In the present study, the water depth of the fluid domain is 5 times deeper than the ship draft. Comparing with the water depth of 2 times of the ship draught as adopted in Marrone et al.[14], the present study avoids the shallow water effect on the ship wave propagation. The ship Froude number considered here is Fr = U / gL = 0.5 with L being the ship length and U the forward ship speed. Figure 15 depicts a panorama view of the ship waves along with the simplified ship. The particles belonging to the free surface are identified through the algorithm presented in Marrone et al.[44] and plotted in blue. The other particles are plo-

208

tted with a contour plot of the FTLE indicated with λT . Only the particles having a λT larger than 0.2, which detect the Lagrangian coherent structures (LCSs) related to the two submerged vortex tubes, are plotted. Through the use of FTLE it is possible to reduce the amount of data for representing the flow feature, which is a remarkable advantage for the analysis and for the storing of 3-D SPH data[185].

rising velocity are plotted in Fig.17. As the particle resolution is increased, it is clearly shown that the SPH results converge to the reference solution.

Fig.15 Breaking wave pattern generated by the simplified ship model at Fr = U /

gL = 0.5

2.1.8 Rising bubble dynamics Rising bubble is a familiar physical phenomenon in ocean engineering from the exploring and processing of marine resources (e.g. nature gas, combustible ice) to small scale bubbles generated by ship wave breaking and propeller cavities. With a multiphase SPH scheme, the rising bubble problems can be simulated properly. Since Colagrossi et al.[22], many efforts have been devoted to modeling of rising bubbles. In Grenier et al.[84], an interface sharpness force was proposed and it plays an important role in preventing the non-physical penetration on the multiphase interface due to the implementation of the surface tension model[186,187]. In Sun et al.[82] and Zhang et al.[17], a comprehensive multiphase SPH model combining different numerical techniques from different models are employed to simulate the rising bubbles in 2-D and 3-D respectively. Here, with the multiphase SPH scheme [17,82], a benchmark of a rising bubble in 3-D is simulated as an example. The parameters for the 3-D benchmark of rising bubble can be found in Adelsberger et al.[188]. To the best of our knowledge, this may be the first attempt that the 3-D rising bubble is quantitatively validated using SPH method. The bubble shapes at different time instants are compared with the results of DROPS[188] as shown in Fig.16 in which the fluid domain is discretized in the particle number of 90×90×180. A fair agreement is obtained between SPH and the reference solution. A convergence study is carried out with three different resolutions, respectively discretizing the fluid domain into particle numbers of 30×30×60, 60×60×120 and 90×90×180, under which the time evolutions of the

Fig.16 The comparisons of the bubble shapes at different time instants

Fig.17 Time evolution of the rising velocity

2.2 SPH simulations of impact and explosion problems Underwater explosion is an important branch in the ocean engineering research. In the off shore operation, the explosion may cause structural destruction which may even lead the platform bilge to capsize. For example, in the 2010 oil spill in the Gulf of Mexico, the “Deepwater Horizon” platform sunk due to the explosion which made the structure lose its integrity. In naval engineering, the underwater explosion is a main threat that makes the warship lose its survivability[53]. The shock wave by the explosion first creates local damages and then the pressure load generated by the explosive bubble may cause the ship to lose its longitudinal strength[189-193]. SPH method has been proven to be a good numerical method to predict the shock wave load of underwater explosions[55].

209

2.2.1 Underwater explosions in a free field It is important to predict the pressure load generated by an underwater explosion since the former is an important reference for the design and evaluation of the strength and survivability of marine structures. 2-D simulations were carried out to validate the SPH results by comparing with the solutions of other numerical results[16,55]. In Zhang et al.[54] and Ming et al.[99], an axisymmetric SPH method was developed in the cylindrical coordinate system to simulate the underwater explosion when the TNT charge has the characteristic of rotational symmetry. A validation for the axisymmetric SPH on predicting the shock wave load is plotted in Fig.18. The present result is calculated with a finer particle resolution than the one used in Ming et al.[99]. Here the particle spacing is ∆x = 0.00532 m. The pressure is measured at R = 6.834 m far away from the center of a spherical TNT charge with the weight of 8 kg. A fair agreement is achieved between SPH and experimental data. The advantage of this technique is its efficiency when compared with purely 3-D simulations. However, it cannot be applied to problems with general boundary shapes.

the FEM is devoted for the elastoplasticity, damage and fracture of ductile structures. The “Glue” algorithm[53] for the fluid-structure coupling is adopted. The Johnson-Cook strength and damage models (see Table 2) are used to model the structure. An air-backed plate subjected to an underwater contact explosion is simulated since the air-backed plate is an element of the ship structure that below the water line. The geometry parameters of the plate and stiffeners are shown on the top of Fig.19. The TNT charge with the weight of 0.16 kg attached to the center of the plate is detonated. The snapshot at t = 2.0 ms is shown on the bottom of Fig.19, from which a petal cracking is observed.

Fig.18 Validation of the shock wave pressure calculated by an axisymmetric SPH model

Fig.19 Parameters and results of the contact explosion problem

Table 2 The parameters for the Johnson-Cook model

3. Future works Particle methods have been recognized as the next generation of numerical tools for multi-physics simulations. Developed for more than three decades, SPH has been proven a successful method solving violent fluid-structure interactions. Despite the great success, there are still some challenges that constrain the application of SPH to industry circles. In the future studies, a few items urgent to be solved or optimized are summarized as follows: (1) In Eulerian CFD solvers, adaptive structured or unstructured grids have been adopted successfully to simulate multiscale flows with variable resolutions[48]. However, in SPH most of the simulations are still done with single particle resolution. More attentions should be paid to the development of models with adaptive particle resolutions or particle refinement and derefinement.

A / MPa

B / MPa

C

n

460

320

0.022

0.36

m

ε0 / s−1

D1

D2 - D5

1.0

1.0

0.1

0

2.2.2 Underwater contact explosions by using an SPHFEM coupling technique The shock wave load generated by an underwater explosion is validated in the previous subsection. In this part, the SPH method is coupled with a FEM solver to simulate the contact explosion which is a great threat to the warship or other marine structures. Recently, in Ming et al.[53], a systematic investigation is conducted by using the SPH-FEM technique, in which SPH is devoted to model fluid dynamics and

210

(2) Inflow and outflow boundaries should be extended to curved or other irregular inlet or outlet shapes. At present, almost all the inlet boundary is a straight wall, which in some cases wastes a lot of computational resources. For example in Colagrossi et al.[194], the simulation of a ship bow wave, a straight inflow boundary demands more space in which some particles do not contribute to the bow wave generation or development. In the same paper, the shape of the inflow boundary of Level Set method is a cured one, which shrinks the scale of the fluid domain. (3) Techniques for preventing Lagrangian particle structures (LPS) or tensile instability need improvement. Particle shifting with strict conservations of momentum and energy should be developed to break LPS[149]. For example in the famous benchmark test case called rotating square fluid patch[50], the disordered particle distribution due to LPS can completely make the numerical results go in a wrong direction. In flow past body in high Reynolds number, the particle shifting cannot be sufficient to prevent the tensile instability and therefore more robust techniques should be developed. (4) A particle re-distribution can be developed. In the explosion problems, after the explosive gas is too much exploded, the particle distribution becomes quite diffused and therefore the accuracy of the particle approximation on the interface between explosion products and fluid is restricted. A proper particle re-distribution (similar to re-mesh) which harmonizes the particle volume is hopeful to overcome this problem. (5) SPH can be further extended to model the complex multiphase flow with a free surface considering the development and collapse of the physical cavitation. For example, the supercavitation flow around a high-speed under water vehicle, the water exit of an underwater launched vehicle or the water entry of a high-speed launched missile. (6) The last issue concerns the application of SPH to industry. Proper pre-processing technique needs to be developed in order to meet the request of modeling simulations with complex 3-D boundaries. Post-processing is also important to obtain a data reduction in 3-D simulations and extract flow information from the scattered data, see more in Sun et al.[24]. A user friendly interface for the SPH solver and visualization to demonstrate SPH results are also significant to drive SPH from academic research to industry applications. Some other techniques[20,195,196] evolved from mesh free methods for fluid-structure interactions can also been applied in ocean engineering applications. Coupling SPH with other techniques[18,19,197,198] is also an important way to extend SPH to wider applications. 4. Concluding remarks In the preset work, a review of the SPH method,

from kernel approximations, governing equations, various numerical techniques to different applications in ocean engineering is presented. Formulations of kernel approximations with different orders of accuracy are recalled. Boundary conditions are reviewed. Techniques for achieving a simulation with multi particle resolutions are introduced. Techniques for breaking the Lagrangian particle structure and preventing tensile instability are also reported. We supplied abundant application cases which are representative in fluid-structure interactions in ocean engineering. The applications are classified into two categories: hydrodynamic problems and explosion problems, which are based on different equations of state and have obviously different features. The various cases show the potential of SPH method to be further extended and applied in the fluid-structure interactions in ocean engineering. Acknowledgment This work was supported by the CNR-INSEAN within the Project PANdA: PArticle methods for Naval Applications, protocol number No. 3263, 21 October 2014. References [1] Gingold R. A., Monaghan J. J. Smoothed particle hydrodynamics-theory and application to non-spherical stars [J]. Monthly Notices of the Royal Astronomical Society, 1977, 181: 375-389. [2] Monaghan J., Gingold R. Shock simulation by the particle method SPH [J]. Journal of Computational Physics, 1983, 52(2): 374-389. [3] Monaghan J. J. Smoothed particle hydrodynamics [J]. Annual Review of Astronomy and Astrophysics, 1992, 30: 543-574. [4] Monaghan J. J. Simulating free surface flows with SPH [J]. Journal of Computational Physics, 1994, 110(2): 399-406. [5] Violeau D., Rogers B. D. Smoothed particle hydrodynamics (SPH) for free-surface flows: Past, present and future [J]. Journal of Hydraulic Research, 2016, 54(1): 1-26. [6] Shadloo M., Oger G., Le Touzé D. Smoothed particle hydrodynamics method for fluid flows, towards industrial applications: Motivations, current state, and challenges [J]. Computers and fluids, 2016, 136: 11-34. [7] Liu M. B., Li S. M. On the modeling of viscous incompressible flows with smoothed particle hydro-dynamics [J]. Journal of Hydrodynamics, 2016, 28(5): 731-745. [8] Gotoh H., Khayyer A. Current achievements and future perspectives for projection-based particle methods with applications in ocean engineering [J]. Journal of Ocean Engineering and Marine Energy, 2016, 2(3): 251-278. [9] Tartakovsky A. M., Trask N., Pan K. et al. Smoothed particle hydrodynamics and its applications for multiphase flow and reactive transport in porous media [J]. Computational Geosciences, 2016, 20(4): 807-834. [10] Cao X., Ming F., Zhang A. Sloshing in a rectangular tank based on SPH simulation [J]. Applied Ocean Research, 2014, 47: 241-254.

211

[11] Cercos-Pita J. L. AQUAgpusph, a new free 3D SPH solver accelerated with OpenCL [J]. Computer Physics Communications, 2015, 192: 295-312. [12] Gong K., Liu H., Wang B. L. Water entry of a wedge based on SPH model with an improved boundary treatment [J]. Journal of Hydrodynamics, 2009, 21(6): 750757. [13] Marrone S., Colagrossi A., Antuono M. et al. An accurate SPH modeling of viscous flows around bodies at low and moderate Reynolds numbers [J]. Journal of Computational Physics, 2013, 245: 456-475. [14] Marrone S., Colagrossi A., Antuono M. et al. A 2D+t SPH model to study the breaking wave pattern generated by fast ships [J]. Journal of Fluids and Structures, 2011, 27(8): 1199-1215. [15] Sun P., Ming F., Zhang A. Numerical simulation of interactions between free surface and rigid body using a robust SPH method [J]. Ocean Engineering, 2015, 98: 32-49. [16] Zhang A. M., Yang W. S., Yao X. L. Numerical simulation of underwater contact explosion [J]. Applied Ocean Research, 2012, 34: 10-20. [17] Zhang A. M., Sun P. N., Ming F. R. An SPH modeling of bubble rising and coalescing in three dimensions [J]. Computer Methods in Applied Mechanics and Engineering, 2015, 294: 189-209. [18] Wang S., Khoo B. C., Liu G. et al. Coupling GSM/ALE with ES-FEM-T3 for fluid–deformable structure interactions [J]. Journal of Computational Physics, 2014, 276: 315-340. [19] Wang S., Khoo B., Liu G. et al. An arbitrary Lagrangian– Eulerian gradient smoothing method (GSM/ALE) for interaction of fluid and a moving rigid body [J]. Computers and fluids, 2013, 71: 327-347. [20] Zhang Z. Q., Liu G., Khoo B. C. Immersed smoothed finite element method for two dimensional fluid–structure interaction problems [J]. International Journal for Numerical Methods in Engineering, 2012, 90(10): 1292-1320. [21] Zhang Z. Q., Liu G., Khoo B. C. A three dimensional immersed smoothed finite element method (3D IS-FEM) for fluid–structure interaction problems [J]. Computational Mechanics, 2013, 51(2): 129-150. [22] Colagrossi A., Landrini M. Numerical simulation of interfacial flows by smoothed particle hydrodynamics [J]. Journal of Computational Physics, 2003, 191(2): 448-475. [23] Khayyer A., Gotoh H. Enhancement of performance and stability of MPS mesh-free particle method for multiphase flows characterized by high density ratios [J]. Journal of Computational Physics, 2013, 242: 211-233. [24] Sun P. N., Colagrossi A., Marrone S. et al. Detection of Lagrangian coherent structures in the SPH framework [J]. Computer Methods in Applied Mechanics and Engineering, 2016, 305: 849-868. [25] Colagrossi A., Rossi E., Marrone S. et al. Particle methods for viscous flows: Analogies and differences between the SPH and DVH methods [J]. Communications in Computational Physics, 2016, 20(3): 660-688. [26] Souto-Iglesias A., Macià F., González L. M. et al. On the consistency of MPS [J]. Computer Physics Communications, 2013, 184(3): 732-745. [27] Souto-Iglesias A., Macià F., González L. M. et al. Addendum to “On the consistency of MPS”[Comput. Phys. Comm. 184 (3)(2013) 732–745] [J]. Computer Physics Communications, 2014, 185(2): 595-598. [28] Shao S., Gotoh H. Turbulence particle models for tracking free surfaces [J]. Journal of Hydraulic Research, 2005, 43(3): 276-289.

[29] Zhang A. M., Cao X. Y., Ming F. R. et al. Investigation on a damaged ship model sinking into water based on three dimensional SPH method [J]. Applied Ocean Research, 2013, 42: 24-31. [30] Marrone S., Bouscasse B., Colagrossi A. et al. Study of ship wave breaking patterns using 3D parallel SPH simulations [J]. Computers and fluids, 2012, 69: 54-66. [31] Oger G., Le Touzé D., Guibert D. et al. On distributed memory MPI-based parallelization of SPH codes in massive HPC context [J]. Computer Physics Communications, 2016, 200: 1-14. [32] Domínguez J. M., Crespo A. J., Valdez-Balderas D. et al. New multi-GPU implementation for smoothed particle hydrodynamics on heterogeneous clusters [J]. Computer Physics Communications, 2013, 184(8): 1848-1860. [33] Longshaw S. M., Rogers B. D. Automotive fuel cell sloshing under temporally and spatially varying high acceleration using GPU-based smoothed particle hydrodynamics (SPH) [J]. Advances in Engineering Software, 2015, 83: 31-44. [34] Valdez-Balderas D., Domínguez J. M., Rogers B. D. et al. Towards accelerating smoothed particle hydrodynamics simulations for free-surface flows on multi-GPU clusters [J]. Journal of Parallel and Distributed Computing, 2013, 73(11): 1483-1493. [35] Antuono M., Colagrossi A., Marrone S. Numerical diffusive terms in weakly-compressible SPH schemes [J]. Computer Physics Communications, 2012, 183(12): 25702580. [36] Antuono M., Colagrossi A., Marrone S. et al. Free-surface flows solved by means of SPH schemes with numerical diffusive terms [J]. Computer Physics Communications, 2010, 181(3): 532-549. [37] Bouscasse B., Colagrossi A., Marrone S. et al. Nonlinear water wave interaction with floating bodies in SPH [J]. Journal of Fluids and Structures, 2013, 42: 112-129. [38] Liu X., Lin P., Shao S. An ISPH simulation of coupled structure interaction with free surface flows [J]. Journal of Fluids and Structures, 2014, 48: 46-61. [39] Gui Q., Dong P., Shao S. Numerical study of PPE source term errors in the incompressible SPH models [J]. International Journal for Numerical Methods in Fluids, 2015, 77(6): 358-379. [40] Gui Q., Shao S., Dong P. Wave impact simulations by an improved ISPH model [J]. Journal of Waterway, Port, Coastal, and Ocean Engineering, 2013, 140(3): 04014005. [41] Shao S. Incompressible smoothed particle hydrodynamics simulation of multifluid flows [J]. International Journal for Numerical Methods in Fluids, 2012, 69(11): 17151735. [42] Daly E., Grimaldi S., Bui H. H. Explicit incompressible SPH algorithm for free-surface flow modelling: A comparison with weakly compressible schemes [J]. Advances in Water Resources, 2016, 97: 156-167. [43] Barcarolo D. A. Improvement of the precision and the efficiency of the SPH method: theoretical and numerical study [D]. Doctoral Thesis, Nantes, France: Ecole Centrale de Nantes, 2013. [44] Marrone S., Colagrossi A., Le Touzé D. et al. Fast freesurface detection and level-set function definition in SPH solvers [J]. Journal of Computational Physics, 2010, 229(10): 3652-3663. [45] Xu R., Stansby P., Laurence D. Accuracy and stability in incompressible SPH (ISPH) based on the projection method and a new approach [J]. Journal of Computational Physics, 2009, 228(18): 6703-6725.

212

[46] Liu M., Liu G. Smoothed particle hydrodynamics (SPH): An overview and recent developments [J]. Archives of computational methods in engineering, 2010, 17(1): 2576. [47] Antuono M., Bouscasse B., Colagrossi A. et al. A measure of spatial disorder in particle methods [J]. Computer Physics Communications, 2014, 185(10): 2609-2621. [48] Yao J., Lin T., Liu G. R. et al. An adaptive GSM-CFD solver and its application to shock-wave boundary layer interaction [J]. International Journal of Numerical Methods for Heat and Fluid Flow, 2015, 25(6): 1282-1310. [49] Barcarolo D., Le Touzé D., Oger G. et al. Adaptive particle refinement and derefinement applied to the smoothed particle hydrodynamics method [J]. Journal of Computational Physics, 2014, 273: 640-657. [50] Sun P. N., Colagrossi A., Marrone S. et al. The δ plusSPH model: Simple procedures for a further improvement of the SPH scheme [J]. Computer Methods in Applied Mechanics and Engineering, 2017, 315: 25-49. [51] Chiron L., Oger G., Touze D. L. et al. Improvements on particle refinement method with SPH [C]. Proceeding of the 11th International SPHERIC Workshop. Munich, Germany, 2016. [52] Zhang Z., Wang L., Silberschmidt V. V. et al. SPH-FEM simulation of shaped-charge jet penetration into double hull: A comparison study for steel and SPS [J]. Composite Structures, 2016, 155: 135-144. [53] Ming F. R., Zhang A. M., Xue Y. Z. et al. Damage characteristics of ship structures subjected to shockwaves of underwater contact explosions [J]. Ocean Engineering, 2016, 117: 359-382. [54] Zhang A. M., Yang W. S., Huang C. et al. Numerical simulation of column charge underwater explosion based on SPH and BEM combination [J]. Computers and fluids, 2013, 71: 169-178. [55] Liu M. B., Liu G. R., Lam K. Y. et al. Smoothed particle hydrodynamics for numerical simulation of underwater explosion [J]. Computational Mechanics, 2003, 30(2): 106-118. [56] Liu M. B., Liu G. R., Zong Z. et al. Computer simulation of high explosive explosion using smoothed particle hydrodynamics methodology [J]. Computers and fluids, 2003, 32(3): 305-322. [57] Liu G. R., Liu M. B. Smoothed particle hydrodynamics: A meshfree particle method [M]. Singapore: World Scientific, 2003. [58] Liu M. B., Liu G. R., Lam K. Y. Constructing smoothing functions in smoothed particle hydrodynamics with applications [J]. Journal of Computational and Applied Mathematics, 2003, 155(2): 263-284. [59] Liu M. B., Xie W. P., Liu G. R. Modeling incompressible flows using a finite particle method [J]. Applied Mathematical Modelling, 2005, 29(12): 1252-1270. [60] Jiang T., Ouyang J., Ren J. L. et al. A mixed corrected symmetric SPH (MC-SSPH) method for computational dynamic problems [J]. Computer Physics Communications, 2012, 183(1): 50-62. [61] Ren J., Jiang T., Lu W. et al. An improved parallel SPH approach to solve 3D transient generalized Newtonian free surface flows [J]. Computer Physics Communications, 2016, 205: 87-105. [62] Long T., Hu D., Yang G. et al. A particle-element contact algorithm incorporated into the coupling methods of FEMISPH and FEM-WCSPH for FSI problems [J]. Ocean Engineering, 2016, 123: 154-163. [63] Serván-Camas B., Cercós-Pita J., Colom-Cobb J. et al.

Time domain simulation of coupled sloshing–seakeeping problems by SPH–FEM coupling [J]. Ocean Engineering, 2016, 123: 383-396. [64] Hu D., Long T., Xiao Y. et al. Fluid–structure interaction analysis by coupled FE–SPH model based on a novel searching algorithm [J]. Computer Methods in Applied Mechanics and Engineering, 2014, 276: 266-286. [65] Li Z., Leduc J., Nunez-Ramirez J. et al. A non-intrusive partitioned approach to couple smoothed particle hydrodynamics and finite element methods for transient fluidstructure interaction problems with large interface motion [J]. Computational Mechanics, 2015, 55(4): 697-718. [66] Ming F. R., Zhang A. M., Wang S. P. Smoothed particle hydrodynamics for the linear and nonlinear analyses of elastoplastic damage and fracture of shell [J]. International Journal of Applied Mechanics, 2015, 7(2): 1550032. [67] Ming F. R., Zhang A. M., Cao X. Y. A robust shell element in meshfree SPH method [J]. Acta Mechanica Sinica, 2013, 29(2): 241-255. [68] Hwang S. C., Park J. C., Gotoh H. et al. Numerical simulations of sloshing flows with elastic baffles by using a particle-based fluid–structure interaction analysis method [J]. Ocean Engineering, 2016, 118: 227-241. [69] Hwang S. C., Khayyer A., Gotoh H. et al. Development of a fully lagrangian MPS-based coupled method for simulation of fluid–structure interaction problems [J]. Journal of Fluids and Structures, 2014, 50: 497-511. [70] Yang X., Liu M., Peng S. et al. Numerical modeling of dam-break flow impacting on flexible structures using an improved SPH–EBG method [J]. Coastal Engineering, 2016, 108: 56-64. [71] Yang X., Liu M., Peng S. Smoothed particle hydrodynamics and element bending group modeling of flexible fibers interacting with viscous fluids [J]. Physical Review E, 2014, 90(6): 063011. [72] YANG X., LIU M. B. Bending modes and transition criteria for a flexible fiber in viscous flows [J]. Journal of Hydrodynamics, 2016, 28(6): 1043-1048. [73] Marrone S., Colagrossi A., Di Mascio A. et al. Analysis of free-surface flows through energy considerations: Singlephase versus two-phase modeling [J]. Physical Review E, 2016, 93(5): 053113. [74] Thiagarajan K., Rakshit D., Repalle N. The air–water sloshing problem: Fundamental analysis and parametric studies on excitation and fill levels [J]. Ocean Engineering, 2011, 38(2): 498-508. [75] Gong K., Shao S., Liu H. et al. Two-phase SPH simulation of fluid–structure interactions [J]. Journal of Fluids and Structures, 2016, 65: 155-179. [76] Lugni C., Brocchini M., Faltinsen O. Evolution of the air cavity during a depressurized wave impact. II. The dynamic field [J]. Physics of Fluids, 2010, 22(5): 056102. [77] Lugni C., Miozzi M., Brocchini M. et al. Evolution of the air cavity during a depressurized wave impact. I. The kinematic flow field [J]. Physics of Fluids, 2010, 22(5): 056101. [78] Chen Z., Zong Z., Liu M. B. et al. An SPH model for multiphase flows with complex interfaces and large density differences [J]. Journal of Computational Physics, 2015, 283: 169-188. [79] Luo X. W., Ji B., Tsujimoto Y. A review of cavitation in hydraulic machinery [J]. Journal of Hydrodynamics, 2016, 28(3): 335-358. [80] Sedlar M., Ji B., Kratky T. et al. Numerical and experimental investigation of three-dimensional cavitating flow around the straight NACA2412 hydrofoil [J]. Ocean Engi-

213

neering, 2016, 123: 357-382. [81] Grenier N., Le Touze D., Colagrossi A. et al. Viscous bubbly flows simulation with an interface SPH model [J]. Ocean Engineering, 2013, 69: 88-102. [82] Sun P. N., Li Y. B., Ming F. R. Numerical simulation on the motion characteristics of freely rising bubbles using smoothed particle hydrodynamics method [J]. Acta Physica Sinica, 2015, 64(17): 174701-174701. [83] Hu X., Adams N. A multi-phase SPH method for macroscopic and mesoscopic flows [J]. Journal of Computational Physics, 2006, 213(2): 844-861. [84] Grenier N., Antuono M., Colagrossi A. et al. An Hamiltonian interface SPH formulation for multi-fluid and free surface flows [J]. Journal of Computational Physics, 2009, 228(22): 8380-8393. [85] Ming F. R., Sun P. N., Zhang A. M. Numerical investigation of rising bubbles bursting at a free surface through a multiphase SPH model [J]. Meccanica, 2017, 1-20. [86] Adami S., Hu X., Adams N. A new surface-tension formulation for multi-phase SPH using a reproducing divergence approximation [J]. Journal of Computational Physics, 2010, 229(13): 5011-5021. [87] Colagrossi A., Antuono M., Le Touzé D. Theoretical considerations on the free-surface role in the smoothed-particle-hydrodynamics model [J]. Physical Review E, 2009, 79(5): 056701. [88] Dehnen W., Aly H. Improving convergence in smoothed particle hydrodynamics simulations without pairing instability [J]. Monthly Notices of the Royal Astronomical Society, 2012, 425(2): 1068-1082. [89] Randles P., Libersky L. Smoothed particle hydrodynamics: Some recent improvements and applications [J]. Computer Methods in Applied Mechanics and Engineering, 1996, 139(1): 375-408. [90] Shao J. R., Li H. Q., Liu G. R. et al. An improved SPH method for modeling liquid sloshing dynamics [J]. Computers and Structures, 2012, 100: 18-26. [91] Huang C., Lei J. M., Liu M. B. et al. A kernel gradient free (KGF) SPH method [J]. International Journal for Numerical Methods in Fluids, 2015, 78(11): 691-707. [92] Marrone S., Antuono M., Colagrossi A. et al. δ - SPH model for simulating violent impact flows [J]. Computer Methods in Applied Mechanics and Engineering, 2011, 200(13): 1526-1542. [93] Ren B., He M., Dong P. et al. Nonlinear simulations of wave-induced motions of a freely floating body using WCSPH method [J]. Applied Ocean Research, 2015, 50: 1-12. [94] Zhang G., Batra R. Symmetric smoothed particle hydrodynamics (SSPH) method and its application to elastic problems [J]. Computational Mechanics, 2009, 43(3): 321-340. [95] Xu X. An improved SPH approach for simulating 3D dam-break flows with breaking waves [J]. Computer Methods in Applied Mechanics and Engineering, 2016, 331: 723–742. [96] Xu X., Deng X. L. An improved weakly compressible SPH method for simulating free surface flows of viscous and viscoelastic fluids [J]. Computer Physics Communications, 2016, 201: 43-62. [97] Benz W. Smoothed particle hydrodynamics: A review [J]. The Numerical Modelling of Nonlinear Stellar Pulsations, 1989, 302: 269-288. [98] Adams B. Simulation of ballistic impacts on armored civil vehicles [D]. Doctoral Thesis, Eindhoven, The Netherlands: Eindhoven University of Technology, 2003.

[99] Ming F. R., Sun P. N., Zhang A. M. Investigation on charge parameters of underwater contact explosion based on axisymmetric SPH method [J]. Applied Mathematics and Mechanics, 2014, 35: 453-468. [100] Shibata K., Koshizuka S., Sakai M. et al. Lagrangian simulations of ship-wave interactions in rough seas [J]. Ocean Engineering, 2012, 42: 13-25. [101] Guo K., Sun P. N., Cao X. Y. et al. A 3-D SPH model for simulating water flooding of a damaged floating structure [J]. Journal of Hydrodynamics, 2017(in Press). [102] Kajtar J., Monaghan J. J. SPH simulations of swimming linked bodies [J]. Journal of Computational Physics, 2008, 227(19): 8568-8587. [103] Shao S., Gotoh H. Simulating coupled motion of progressive wave and floating curtain wall by SPH-LES model [J]. Coastal Engineering Journal, 2004, 46(2): 171-202. [104] Cercos-Pita J., Dalrymple R., Herault A. Diffusive terms for the conservation of mass equation in SPH [J]. Applied Mathematical Modelling, 2016, 40(19-20): 8722-8736. [105] Molteni D., Colagrossi A. A simple procedure to improve the pressure evaluation in hydrodynamic context using the SPH [J]. Computer Physics Communications, 2009, 180(6): 861-872. [106] Vila J. On particle weighted methods and smooth particle hydrodynamics [J]. Mathematical models and methods in applied sciences, 1999, 9(02): 161-209. [107] Lanson N., Vila J. P. Renormalized meshfree schemes I: Consistency, stability, and hybrid methods for conservation laws [J]. SIAM Journal on Numerical Analysis, 2008, 46(4): 1912-1934. [108] Ferrari A., Dumbser M., Toro E. F. et al. A new 3D parallel SPH scheme for free surface flows [J]. Computers and fluids, 2009, 38(6): 1203-1217. [109] Khayyer A., Gotoh H. Enhancement of stability and accuracy of the moving particle semi-implicit method [J]. Journal of Computational Physics, 2011, 230(8): 30933118. [110] Khayyer A., Gotoh H. A higher order Laplacian model for enhancement and stabilization of pressure calculation by the MPS method [J]. Applied Ocean Research, 2010, 32(1): 124-131. [111] Colagrossi A., Antuono M., Souto-Iglesias A. et al. Theoretical analysis and numerical verification of the consistency of viscous smoothed-particle-hydrodynamics formulations in simulating free-surface flows [J]. Physical Review E, 2011, 84(2): 026705. [112] Dilts G. A. Moving least-squares particle hydrodynamics II: Conservation and boundaries [J]. International Journal for Numerical Methods in Engineering, 2000, 48(10): 1503-1524. [113] Haque A., Dilts G. A. Three-dimensional boundary detection for particle methods [J]. Journal of Computational Physics, 2007, 226(2): 1710-1730. [114] Zheng X., Duan W. Y., Ma Q. W. A new scheme for identifying free surface particles in improved SPH [J]. Science China Physics, Mechanics and Astronomy, 2012, 55(8): 1454-1463. [115] Tang Z. Y., Zhang Y. L., Wan D. C. Numerical simulation of 3-D free surface flows by overlapping MPS [J]. Journal of Hydrodynamics, 2016, 28(2): 306-312. [116] De Leffe M., Le Touzé D., Alessandrini B. Normal flux method at the boundary for SPH [C]. 4th ERCOFTAC SPHERIC Workshop. Nantes , France, 2009. [117] Ferrand M., Laurence D., Rogers B. et al. Unified semianalytical wall boundary conditions for inviscid, laminar

214

or turbulent flows in the meshless SPH method [J]. International Journal for Numerical Methods in Fluids, 2013, 71(4): 446-472. [118] Monaghan J., Kajtar J. SPH particle boundary forces for arbitrary boundaries [J]. Computer Physics Communications, 2009, 180(10): 1811-1820. [119] Cummins S. J., Rudman M. An SPH projection method [J]. Journal of Computational Physics, 1999, 152(2): 584-607. [120] Adami S., Hu X., Adams N. A generalized wall boundary condition for smoothed particle hydrodynamics [J]. Journal of Computational Physics, 2012, 231(21): 70577075. [121] Liu M. B., Shao J. R., Chang J. Z. On the treatment of solid boundary in smoothed particle hydrodynamics [J]. Science China-Technological Sciences, 2012, 55(1): 244254. [122] Liu M. B., Shao J. R., Li H. Q. An SPH model for free surface flows with moving rigid objects [J]. International Journal for Numerical Methods in Fluids, 2014, 74: 684697. [123] Cercos-Pita J., Antuono M., Colagrossi A. et al. SPH energy conservation for fluid–solid interactions [J]. Computer Methods in Applied Mechanics and Engineering, 2017, 317: 771-791. [124] Maciá F., Antuono M., González L. M. et al. Theoretical analysis of the no-slip boundary condition enforcement in SPH methods [J]. Progress of theoretical physics, 2011, 125(6): 1091-1121. [125] Shao S., Lo E. Y. Incompressible SPH method for simulating Newtonian and non-Newtonian flows with a free surface [J]. Advances in Water Resources, 2003, 26(7): 787-800. [126] Takeda H., Miyama S. M., Sekiya M. Numerical simulation of viscous flow by smoothed particle hydrodynamics [J]. Progress of theoretical physics, 1994, 92(5): 939-960. [127] Federico I., Marrone S., Colagrossi A. et al. Simulating 2D open-channel flows through an SPH model [J]. European Journal of Mechanics-B/Fluids, 2012, 34: 3546. [128] Kazemi E., Nichols A., Tait S. et al. SPH modelling of depth-limited turbulent open channel flows over rough boundaries [J]. International Journal for Numerical Methods in Fluids, 2017, 83(1): 3-27. [129] Ferrand M., Joly A., Kassiotis C. et al. Unsteady open boundaries for SPH using semi-analytical conditions and Riemann solver in 2D [J]. Computer Physics Communications, 2016, 210: 29-44. [130] Leroy A., Violeau D., Ferrand M. et al. A new open boundary formulation for incompressible SPH [J]. Computers and Mathematics with Applications, 2016, 72(9): 2417-2432. [131] Kunz P., Hirschler M., Huber M. et al. Inflow/outflow with dirichlet boundary conditions for pressure in ISPH [J]. Journal of Computational Physics, 2016, 326: 171-187. [132] Tan S. K., Cheng N. S., Xie Y. et al. Incompressible SPH simulation of open channel flow over smooth bed [J]. Journal of Hydro-environment Research, 2015, 9(3): 340-353. [133] Lastiwka M., Basa M., Quinlan N. J. Permeable and nonreflecting boundary conditions in SPH [J]. International Journal for Numerical Methods in Fluids, 2009, 61(7): 709-724. [134] Wen H., Ren B. 3D Numerical wave basin based on

parallelized SPH method [C]. ASME 33rd International Conference on Ocean, Offshore and Arctic Engineering. San Francisco, USA, 2014. [135] Colagrossi A., Bouscasse B., Antuono M. et al. Particle packing algorithm for SPH schemes [J]. Computer Physics Communications, 2012, 183(8): 1641-1653. [136] Antuono M., Colagrossi A., Marrone S. et al. Propagation of gravity waves through an SPH scheme with numerical diffusive terms [J]. Computer Physics Communications, 2011, 182(4): 866-877. [137] Rossi E., Colagrossi A., Durante D. et al. Simulating 2D viscous flow around geometries with vertices through the diffused vortex hydrodynamics method [J]. Computer Methods in Applied Mechanics and Engineering, 2016, 302: 147-169. [138] Marsh A., Oger G., Touzé D. l. et al. Validation of a conservative variable-resolution SPH scheme including ∇h terms [C]. Proceedings of the 6th International SPHERIC Workshop. Hambourg, Germany, 2011. [139] Koukouvinis P. K., Anagnostopoulos J. S., Papantonis D. E. Simulation of 2D wedge impacts on water using the SPH–ALE method [J]. Acta Mechanica, 2013, 224(11): 2559-2575. [140] Tang Z., Wan D., Chen G. et al. Numerical simulation of 3D violent free-surface flows by multi-resolution MPS method [J]. Journal of Ocean Engineering and Marine Energy, 2016, 2(3): 355-364. [141] Feldman J., Bonet J. Dynamic refinement and boundary contact forces in SPH with applications in fluid flow problems [J]. International Journal for Numerical Methods in Engineering, 2007, 72(3): 295-324. [142] Omidvar P., Stansby P. K., Rogers B. D. Wave body interaction in 2D using smoothed particle hydrodynamics (SPH) with variable particle mass [J]. International Journal for Numerical Methods in Fluids, 2012, 68(6): 686705. [143] Omidvar P., Stansby P. K., Rogers B. D. SPH for 3D floating bodies using variable mass particle distribution [J]. International Journal for Numerical Methods in Fluids, 2013, 72(4): 427-452. [144] López Y. R., Roose D., Morfa C. R. Dynamic particle refinement in SPH: Application to free surface flow and non-cohesive soil simulations [J]. Computational Mechanics, 2013, 51(5): 731-741. [145] Vacondio R., Rogers B., Stansby P. et al. Variable resolution for SPH in three dimensions: Towards optimal splitting and coalescing for dynamic adaptivity [J]. Computer Methods in Applied Mechanics and Engineering, 2016, 300: 442-460. [146] Vacondio R., Rogers B., Stansby P. et al. Variable resolution for SPH: A dynamic particle coalescing and splitting scheme [J]. Computer Methods in Applied Mechanics and Engineering, 2013, 256: 132-148. [147] Monaghan J. J. SPH without a tensile instability [J]. Journal of Computational Physics, 2000, 159(2): 290311. [148] Le Touzé D., Colagrossi A., Colicchio G. et al. A critical investigation of smoothed particle hydrodynamics applied to problems with free-surfaces [J]. International Journal for Numerical Methods in Fluids, 2013, 73(7): 660-691. [149] Oger G., Marrone S., Le Touzé D. et al. SPH accuracy improvement through the combination of a quasiLagrangian shifting transport velocity and consistent ALE formalisms [J]. Journal of Computational Physics, 2016, 313: 76-98.

215

[150] Tsuruta N., Khayyer A., Gotoh H. A short note on dynamic stabilization of moving particle semi-implicit method [J]. Computers and fluids, 2013, 82: 158-164. [151] Adami S., Hu X., Adams N. A transport-velocity formulation for smoothed particle hydrodynamics [J]. Journal of Computational Physics, 2013, 241: 292-307. [152] Litvinov S., Hu X., Adams N. Towards consistence and convergence of conservative SPH approximations [J]. Journal of Computational Physics, 2015, 301: 394-401. [153] Lind S. J., Xu R., Stansby P. K. 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]. Journal of Computational Physics, 2012, 231(4): 1499-1523. [154] Liu X., Lin P., Shao S. ISPH wave simulation by using an internal wave maker [J]. Coastal Engineering, 2015, 95: 160-170. [155] Bouscasse B., Antuono M., Colagrossi A. et al. Numerical and experimental investigation of nonlinear shallow water sloshing [J]. International Journal of Nonlinear Sciences and Numerical Simulation, 2013, 14(2): 123138. [156] Delorme L., Colagrossi A., Souto-Iglesias A. et al. A set of canonical problems in sloshing, Part I: Pressure field inforced roll comparison between experimental results and SPH [J]. Ocean Engineering, 2009, 36(2): 168-178. [157] Marrone S., Colagrossi A., Di Mascio A. et al. Prediction of energy losses in water impacts using incompressible and weakly compressible models [J]. Journal of Fluids and Structures, 2015, 54: 802-822. [158] Chen Z., Zong Z., Li H. et al. An investigation into the pressure on solid walls in 2D sloshing using SPH method [J]. Ocean Engineering, 2013, 59: 129-141. [159] Luo M., Koh C., Bai W. A three-dimensional particle method for violent sloshing under regular and irregular excitations [J]. Ocean Engineering, 2016, 120: 52-63. [160] Wei Z., Hu C. Experimental study on water entry of circular cylinders with inclined angles [J]. Journal of Marine Science and Technology, 2015, 20(4): 722-738. [161] Wei Z., Hu C. An experimental study on water entry of horizontal cylinders [J]. Journal of Marine Science and Technology, 2014, 19(3): 338-350. [162] Nguyen V. T., Vu D. T., Park W. G. et al. Navier-Stokes solver for water entry bodies with moving Chimera grid method in 6DOF motions [J]. Computers and fluids, 2016, 140: 19-38. [163] Zhu X., Faltinsen O. M., Hu C. Water entry and exit of a horizontal circular cylinder [J]. Journal of Offshore Mechanics and Arctic Engineering, 2007, 129(4): 253-264. [164] Zhao R., Faltinsen O. Water entry of two-dimensional bodies [J]. Journal of Fluid Mechanics, 1993, 246(1): 593-612. [165] Wu G., Sun H., He Y. Numerical simulation and experimental study of water entry of a wedge in free fall motion [J]. Journal of Fluids and Structures, 2004, 19(3): 277-289. [166] Sun S., Wu G. Oblique water entry of a cone by a fully three-dimensional nonlinear method [J]. Journal of Fluids and Structures, 2013, 42: 313-332. [167] Sun H., Faltinsen O. M. Water impact of horizontal circular cylinders and cylindrical shells [J]. Applied Ocean Research, 2006, 28(5): 299-311. [168] Vandamme J., Zou Q., Reeve D. E. Modeling floating object entry and exit using smoothed particle hydrodynamics [J]. Journal of Waterway, Port, Coastal, and Ocean

Engineering, 2011, 137(5): 213-224. [169] Colicchio G., Greco M., Miozzi M. et al. Experimental and numerical investigation of the water-entry and waterexit of a circular cylinder [C]. Proceedings of the 24th International Workshop on Water Waves and Floating Bodies. Zelenogorsk, Russia, 2009, 19-22. [170] Zhu X. Application of the CIP method to strongly nonlinear wave-body interaction problems [D]. Doctoral Thesis, Trondheim, Norway: Norwegian University of Science and Technology, 2006. [171] Ni B. Y., Zhang A. M., Wu G. X. Simulation of complete water exit of a fully-submerged body [J]. Journal of Fluids and Structures, 2015, 58: 79-98. [172] Wang S., Zhang G., Feng S. et al. Numerical simulation for 2D water-exit problem based on boundary element method [J]. Journal of Dalian Maritime University, 2016, 42(4): 26-32(in Chinese). [173] Gang M. Hydrodynamic forces and dynamic responses of circular cylinders on wave zones [D]. Doctoral Thesis, Trondheim, Norway: Norwegian Institute of Technology, 1989. [174] Pan K., IJzermans R., Jones B. et al. Application of the SPH method to solitary wave impact on an offshore platform [J]. Computational Particle Mechanics, 2016, 3(2): 155-166. [175] Zhao X., Hu C. Numerical and experimental study on a 2-D floating body under extreme wave conditions [J]. Applied Ocean Research, 2012, 35: 1-13. [176] Zhao X., Ye Z., Fu Y. et al. A CIP-based numerical simulation of freak wave impact on a floating body [J]. Ocean Engineering, 2014, 87: 50-63. [177] Zhao X. Z., Hu C. H. Numerical and experimental study on a 2-D floating body under extreme wave conditions [J]. Applied Ocean Research, 2012, 35: 1-13. [178] Idelsohn S., Onate E., Del Pin F. et al. Fluid–structure interaction using the particle finite element method [J]. Computer Methods in Applied Mechanics and Engineering, 2006, 195(17): 2100-2123. [179] Le Touzé D., Marsh A., Oger G. et al. SPH simulation of green water and ship flooding scenarios [J]. Journal of Hydrodynamics, 2010, 22(5): 231-236. [180] Rossi E., Colagrossi A., Bouscasse B. et al. The diffused vortex hydrodynamics method [J]. Communications in Computational Physics, 2015, 18(2): 351-379. [181] Bouscasse B., Colagrossi A., Marrone S. et al. Viscous flow past a circular cylinder below a free surface [C]. ASME 33rd International Conference on Ocean, Offshore and Arctic Engineering. San Francisco, USA, 2014. [182] Bouscasse B., Colagrossi A., Marrone S. et al. High froude number viscous flow past a circular cylinder [C]. ASME 34th International Conference on Ocean, Offshore and Arctic Engineering. St. Johnʼs, Newfoundland, Canada, 2015. [183] Cercos-Pita J. L., Colagrossi A., Souto-Iglesias A. Low Reynolds flow past a circular cylinder close to a freesurface with vertical motion dynamics [C]. ASME 35th International Conference on Ocean, Offshore and Arctic Engineering, Busan, Korea, 2016. [184] Bouscassea B., Colagrossi A., Marronea S. et al. SPH modelling of viscous flows past a circular cylinder interacting with a free surface [J]. Computers and fluids, 2017. [185] Sun P. N., Colagrossi A., Marrone S. et al. SPH formulation of FTLE for the detection of Lagrangian coherent structures [C]. Proceeding of the 11th Int. SPHERIC Workshop. Munich, Germany, 2016.

216

[186] Landrini M., Colagrossi A., Greco M. et al. The fluid mechanics of splashing bow waves on ships: A hybrid BEM–SPH analysis [J]. Ocean Engineering, 2012, 53: 111-127. [187] Morris J. P. Simulating surface tension with smoothed particle hydrodynamics [J]. International Journal for Numerical Methods in Fluids, 2000, 33(3): 333-353. [188] Adelsberger J., Esser P., Griebel M. et al. 3D incompressible two-phase flow benchmark computations for rising droplets [C]. Proceedings of the 11th World Congress on Computational Mechanics (WCCM XI). Barcelona, Spain, 2014. [189] Zhang A., Zhou W., Wang S. et al. Dynamic response of the non-contact underwater explosions on naval equipment [J]. Marine Structures, 2011, 24(4): 396-411. [190] Zhang A. M., Zeng L. Y., Cheng X. D. et al. The evaluation method of total damage to ship in underwater explosion [J]. Applied Ocean Research, 2011, 33(4): 240-251. [191] Zhang A. M., Liu Y. L. Improved three-dimensional bubble dynamics model based on boundary element method [J]. Journal of Computational Physics, 2015, 294: 208-223. [192] Li S., Han R., Zhang A. M. et al. Analysis of pressure field generated by a collapsing bubble [J]. Ocean Engineering, 2016, 117: 22-38.

[193] Li S., Li Y. B., Zhang A. M. Numerical analysis of the bubble jet impact on a rigid wall [J]. Applied Ocean Research, 2015, 50: 227-236. [194] Colagrossi A., Marrone S., Bouscasse B. et al. Numerical simulations of the flow past surface-piercing objects [J]. International Journal of Offshore and Polar Engineering, 2015, 25(1): 13-18. [195] Zhang Z. Q., Yao J., Liu G. R. An immersed smoothed finite element method for fluid–structure interaction problems [J]. International Journal of Computational Methods, 2011, 8(4): 747-757. [196] Yao J., Liu G., Narmoneva D. A. et al. Immersed smoothed finite element method for fluid–structure interaction simulation of aortic valves [J]. Computational Mechanics, 2012, 50(6): 789-804. [197] He Z. C., Liu G. R., Zhong Z. H. et al. A coupled ES-FEM/BEM method for fluid–structure interaction problems [J]. Engineering Analysis with Boundary Elements, 2011, 35(1): 140-147. [198] Marrone S., Di Mascio A., Le Touzé D. Coupling of smoothed particle hydrodynamics with finite volume method for free-surface flows [J]. Journal of Computational Physics, 2016, 310: 161-180.