Two-phase SPH model based on an improved Riemann solver for water entry problems

Ocean Engineering 199 (2020) 107039

Contents lists available at ScienceDirect

Ocean Engineering journal homepage: www.elsevier.com/locate/oceaneng

Two-phase SPH model based on an improved Riemann solver for water entry problems Qiuzu Yang a, b, Fei Xu a, b, *, Yang Yang a, b, **, Jingyu Wang c a

School of Aeronautics, Northwestern Polytechnical University, 710072, Xi’an, PR China Institute for Computational Mechanics and Its Application, Northwestern Polytechnical University, 710072, Xi’an, PR China c Department of Mechanical Engineering, Technical University of Munich, 85748, Garching, Germany b

A R T I C L E I N F O

A B S T R A C T

Keywords: Two-phase flows Fluid-structure coupling SPH (smoothed particle hydrodynamics) Riemann solver Water entry

There are challenges in simulating two-phase flows with large density ratios using traditional SPH. Taking full advantages of Riemann solver in dealing with contact discontinuity problems and combining two-phase flow SPH method with Riemann solver, a new model of two-phase flow coupling with structure based on an improved Riemann solver is described for the water entry problems with the effect of air. The flow momentum equation of the Riemann form is improved to decrease the Riemann dissipation and simulate the two-phase flows with different viscosity ratios. One-sided Riemann problem considering the free-slip and no-slip conditions at a wall is used to deal with the fluid-body interaction. A switch-function-based Riemann solver dissipation is used to improve the instability of interface owing to the strong impact. In order to test the accuracy of the model in fluidstructure coupling problems, five cases, including dam break, sinking of a rectangle body, water entry of a wedge, flat plate slamming on water and water entry of a catamaran, are numerically simulated. The agreements between the present results and other reference results demonstrate that the proposed method is robust and stable to simulate the water entry problems with the air.

1. Introduction In the ocean engineering field, fluid-structure coupling problems, such as diving of sportsman, ship slamming, waves slapping and landing of aircraft on water, are important subject matter. The accurate pre­ diction of the structure dynamics response can provide guide and reference for engineering design (e.g. the ship design in Naval Archi­ tecture and Ocean Engineering). The structure entering water is a typical fluid-structure coupling problem which involves complex hydrody­ namics and kinetics problems (Wang et al., 2015). When the structure enters the water, due to the strong impact, structure will produce rolling, pitching and vertical movement (which causes internal device mal­ function or even local damage of structure), and the water surface will produce splashing, breaking and rolling phenomena. Therefore, it is worthy to make further studies. Numerical study of water entry problems using the traditional Eulerian grid-based numerical methods have been challenging. Considering the influence of the air on the water entry structure, the difficulties will increase. To capture the large free interface deformation

and track the moving objects, special algorithms need to be added, which will complicate the algorithms (Lin, 2007; Wick, 2011). Mesh-free numerical methods, such as Smoothed Particle Hydrody­ namics (SPH) (Monaghan, 1985) and Moving Particle Semi-implicit (MPS) (Koshizuka and Oka, 1996), can overcome above difficulties and have the advantage of modeling the water entry problems. Oger et al. (2006) was the first to apply the weakly compressible SPH (WSPH) to simulate the water entry problems, and a technic was introduced to evaluate fluid pressure on solid boundaries. Then Maruzewski et al. (2010) further developed the SPH method for the water entry problems. Afterwards, Skillen et al. (2013) adopted incompressible SPH (ISPH) to simulate the water entry of the wedge and cylinder. Sun et al. (2018) studied 2D and 3D water-entry cylinder using the developed δþ-SPH, in which the adaptive particle refinement (APR) was applied to reduce the computational cost. To accurately predict the high and localized impact pressure, Marrone et al. (2018) used a Riemann-ALE smoothed particle hydrodynamics model to simulate the high speed water entry of plates. In recent literature, Liu and Zhang (2019) reviewed applications of the SPH method in various water entry problems.

* Corresponding author. School of Aeronautics, Northwestern Polytechnical University, 710072, Xi’an, PR China. ** Corresponding author. School of Aeronautics, Northwestern Polytechnical University, 710072, Xi’an, PR China. E-mail addresses: [email protected] (F. Xu), [email protected] (Y. Yang). https://doi.org/10.1016/j.oceaneng.2020.107039 Received 9 October 2019; Received in revised form 29 December 2019; Accepted 29 January 2020 Available online 7 February 2020 0029-8018/© 2020 Elsevier Ltd. All rights reserved.

Q. Yang et al.

Ocean Engineering 199 (2020) 107039

However, above simulations of water entry problems generally are limited in the scope of a single-phase flow coupling with structure, none of these works investigate the effect of air. Due to the discontinuity of property across the interface, the simulation of two-phase flow problems is difficult with the traditional particle methods, especially in the multiphase flow problems with large density ratios. In view of this, many scholars have conducted studies on the simulation of the multi-phase flow and multi-phase flow coupling with structure by particle methods. Monaghan and Kocharyan (1995) were the first to put forward the concept of multi-phase SPH model for the simulation of two-phase flows with small density ratios. Subsequently, to steady the interface of multi-phase flows with large density ratios, Colagrossi and Landrini (2003), Hu and Adams (2006, 2007) and Grenier et al. (2009, 2013) presented various modifications including density reinitialization, arti­ ficial viscosity, cohesive force, XSPH technique and so on. Khayyer and Gotoh (2013) combined an extended version of SPH density smooth­ ening scheme with MPS method to enhance the numerical stability of multi-phase flows with high density ratios. Later, Khayyer et al. (2019) improved previously schemes and proposed a novel projection-based particle method for simulation of multi-phase flows which was called the improved MPS þ OPS (Optimized Particle Shifting). For the studies of multi-phase flow coupling with structure in Lind et al. (2015), it proved that the air cushion had an important effect on the slamming load subjected to entry structure with a low dead-rise angle like plate. Yan et al. (2015) studied the water entry catamaran hull, where the air was trapped beneath the hull and the effect of air was pronounced again. Khayyer and Gotoh (2016) proposed to adopt a modified Error-Compensating Source term of Poisson Pressure Equations in the framework of multi-phase MPS to model water entry of a rigid plate, and the air cushioning effect was validated. Based on Riemann-ALE smoothed particle hydrodynamics model, Marrone et al. (2017) considered the influence of air cushioning to study the water entry of the 3D complex geometries. In the above simulation of the two-phase flow coupling with structure problems, pressure oscillation is always a tricky problem to deal with. In order to reduce pressure oscillation of multi-phase flow field produced by the impacting of water entry, most researches introduced the artificial viscosity term proposed by Monaghan (1994) into the momentum equation. But the artificial viscosity is not the real viscosity of fluid, which easily affects the real motion characteristics of the fluid and energy conservation. By reconstructing the SPH convolution in­ tegrals, a Godunov SPH (GSPH) algorithm with no artificial viscosity was proposed by Inutsuka (1994) where the Riemann problem was applied to construct an imaginary interface between two particles. Later, various kinds of Godunov SPH schemes had been proposed in the literature (Cha and Whitworth, 2003; Sirotkin and Yoh, 2013; Puri and Ramachandan, 2014; Yang et al., 2020). Puri and Ramachandan (2014) established the equivalence of the dissipative terms in GSPH to the signal-based artificial viscosity in SPH and pointed that the GSPH dissipation helped to prevent particle penetration in some cases, such as subsonic and supersonic flow collision tests (Cha and Whitworth, 2003). Teschner et al. (2019) utilized the numerical dissipation introduced by the Riemann solver to preserve pressure stability of steady state flows. Taking full advantages of Riemann solver in dealing with contact discontinuity problems, Rezavand et al. (2019) proposed a multi-phase SPH method to simulate highly violent flows problems on the base of the Zhang and Hu (2017) study. But this multi-phase SPH model cannot simulate the two-phase flow problems with various viscosity ratios and has limitations in application. Therefore, this paper develops a set of novel gas-liquid two-phase SPH model based on Riemann solver which can simulate the problems of the viscous two-phase flow coupling with structure. In the present work, the momentum equation of the Riemann form is improved to make up for the defects of Rezavand et al. (2019) model and reduce the Riemann dissipation. For the sake of establishing a consistent SPH model based on Riemann solver, one-sided Riemann problem proposed by Zhang and Hu

(2017) is expanded and used to impose the interaction between two-phase flow and structure, which is different from the previous literature (Zhang and Hu, 2017) and can deal with both the no-slip boundary and free-slip boundary. And a switch-function-based Rie­ mann solver dissipation, as a new lower numerical diffusion term, is used to improve the instability of flow field where the strong impact occurs. This imposed dissipation is different from that of the global application of the artificial diffusion used in the previous studies (Yan et al., 2015; Cao et al., 2018). This paper is organized as follows. In section 2, governing equations of two-phase flow and the motion equations of moving object are pre­ sented. The detailed derivations of the proposed two-phase SPH model based on an improved Riemann solver are reported in section 3, including some strategies to ensure the stability and accuracy. In section 4, five representative cases, including dam-break, sinking of a rectangle body, water-entry of a wedge, flat plate slamming on water and water entry of a catamaran, are numerically simulated by the present twophase SPH model. Several conclusions are drawn in section 5. 2. Fluid-structure interact formulation The governing equations of the two-phase fluid and the motion equations of structure are introduced first, the following SPH formula­ tions are descripted, and the two-phase SPH model based on an improved Riemann solver for simulating fluid-structure problems is presented in the section 3. 2.1. Governing equations for two-phase fluid In the present work, the mass and momentum conservation equa­ tions of two-phase fluid based on the compressible Navier-Stokes equations are: 8 dρ > > < dt ¼ ρr⋅v (1) > dv 1 > : ¼ rp þ g dt ρ where ρ, v, p and g refer to the density, the velocity, the pressure and gravitational acceleration, respectively. In order to fulfill the incom­ pressible limit, the two-phase fluid is assumed to be weaklycompressible (Monaghan and Gingold, 1983; Morris et al., 1997). The Eq. (1) is non-closed, the Tait equation of state is used to solve the two-phase fluid evolution. Moreover, the form of the Tait equation of state is applied: � �� �γ ρ c2 ρ p¼ 0 1 þ p0 (2) ρ0 γ where ρ0 is the reference density. For water, the initial density is ρw0 ¼ 1000 kg/m3; while for air, the initial density is ρa0 ¼ 1.29 kg/m3. For the heavier phase, γ ¼ 7 is used and γ ¼ 1.4 is used in the lighter phase as suggested in literature (Zhang et al., 2015; Colagrossi and Landrini, 2003). p0 is the background pressure which can keep the particle uni­ form distribution and prevent the tension instability induced by nega­ tive pressure. The background pressure is only used in two-phase flow problems, and its value should be carefully chosen (Gong et al., 2016). To ensure the initial pressure is smooth at the interface, the different phases would have the same coefficient ρ0c2/γ. c0 is the artificial sound speed which is chosen to keep the variation of density within 1%. So, the sound speed c0w of the heavier phase meets conditions as follows (Marrone et al., 2017): � pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi� c0w > 10max vmax ; pmax =ρ0 (3) where vmax and pmax are the maximum expected velocity and pressure, respectively. And the sound speed c0h of the lighter phase is 2

Q. Yang et al.

Ocean Engineering 199 (2020) 107039

(a)

(b)

Fig. 1. The Riemann problem between a pair of particles i and j in 2D space. (a) The Riemann problem along the line joining particle i to j; (b) wave solution of Riemann solver.

c0h¼(ρ0wc20wγ h/ρ0hγ w)0.5.

chosen to be 1.25Δx, and Δx is the initial particle spacing. W(ri-rj,h) and riW(ri-rj,h) are the kernel function and its gradient with respect to the particle i respectively. In this paper, the kernel function proposed by Wendland (1995) is applied: � � αð2 qÞ4 ð2q þ 1Þ; 0�q<2 W ri rj ; h ¼ (7) 0; q�2

2.2. The motion equations for moving object We consider the rotation of the object during the motion. Based on the Newton’s law of motion, the motion equations of the object’s center of mass can be written as follows: 8 dv0 F > > ¼ þg < M dt (4) > > : dw0 ¼ J I dt

where, q ¼ |ri-rj|/h, one-dimensional problems: α ¼ 3/64h, twodimensional problems: α ¼ 7/64h2, three-dimensional problems: α ¼ 21/256h3. Through the transformation of SPH method in Eq. (6), gov­ erning equations for two-phase fluid in Eq. (1) can be expressed as: � 8 X mj X mj � � � � d ρi > > vi vj ⋅ri Wij ri rj ;h ¼2ρi vi vij ⋅ri Wij ri rj ;h > dt ¼ ρi > ρj ρj < j j

where v0 and w0 denote the velocity and the angular velocity of the object’s center of mass, respectively. F, M, J and I are the summation of forces without gravity, the mass of the object, the total moment of force on the object about the center of mass and the moment of inertia, respectively. The velocity of a point i on the object is:

> � � � dvi 1 X mj 2 X mj > > > p þp r W r r ;h þg¼ p r W r r ;h þg : dt ¼ ρ ρ i j i ij i j ρ ρ ij i ij i j i

i

j

j

where vij ¼ ðvi þvj Þ=2 and pij ¼ ðpi þpj Þ=2 are the average velocity and pressure of particle i and particle j, respectively.

where rio is the vector rio ¼ ri- ro.

3.2. The improved SPH model based on Riemann solver

3. The present two-phase SPH model

The new SPH model, unlike the conventional SPH method describing the interaction between particles only by the weighting summation of the neighboring particles in the support domain, is based on Riemann solver where each pair of particles is regarded as a Riemann problem. This Riemann problem is constructed along the vector eij ¼ -(ri-rj)/| rirj|, and this interaction point is located at the midway point between particle i and particle j (Puri and Ramachandan, 2014; Zhang and Hu, 2017). As shown in Fig. 1(a), the left and right states of a constructed inter-particle Riemann problem are shown as � � ðρL ; uL ; PL Þ ¼ ρi ; vi ⋅eij ; pi � (9) ðρR ; uR ; PR Þ ¼ ρj ; vj ⋅eij ; pj

3.1. Basic SPH Smoothed particle hydrodynamics (SPH) is a Lagrangian meshless method in which the continuous matter field is discretized into particles possessing physical properties (such as density, velocity, acceleration, pressure and viscosity) (Liu and Liu, 2003). The governing equations of two-phase fluid and moving object will be solved by SPH method. An arbitrary field function f(ri) and its derivative rf(ri) at particle i position can be approximated by kernel interpolation in its supported domain as (Monaghan, 1985): Z X mj � � 8 f ðri Þ ¼ f ðr’ ÞWðri r’ ; hÞd r’ � f r W ri rj ; h > > ρj j > > j > Ω > > > Z > < X mj � � rf ðri Þ ¼ f ðr’ Þri Wðri r’ ; hÞdr’ � f rj ri W ri rj ; h (6) ρ > j > j > Ω > > > > X mj � �� � > > : f ðri Þ þ f rj ri W ri rj ; h ¼ j

j

(8)

(5)

vi ¼ v0 þ w0 � rio

j

The contact discontinuity of Riemann problem is located at r ¼ ðri þ rj Þ=2. The first-order Riemann solver is used to approximate the Rie­ mann problem (four states are divided by three waves). As shown in Fig. 1(b). Assume that the intermediate states on both sides of the middle wave are equivalent. Similarly to Zhang and Hu (2017), a linearized Riemann solver is used to approximate the intermediate states: 8 1 PL PR * * * > > > uL ¼ uR ¼ uij ¼ u þ 2 ρ c < ij ij (10) > > > : P* ¼ P* ¼ P* ¼ p þ 1ρ CðuL uR Þ L R ij 2 ij

ρj

where Ω is the support domain of the smooth function. The subscript j represents the j-th particle within the supported domain of the particle i. mj and ρj are the mass and density of particle j respectively. h refers to the smoothing length of the area affected by the kernel function which is 3

Q. Yang et al.

Ocean Engineering 199 (2020) 107039

where u; ρ; c and P are the inter-particle velocity, density, sound speed and pressure, respectively. C is the inter-particular sound speed which is related to the Riemann dissipation in SPH model, see more in Puri and Ramachandan (2014). In the two-phase flows, the abrupt discontinuities of density and viscosity across the interface contribute to the numerical instability, especially in the water entry problems with the effect of air. At the moment of impact, the pressure of two-phase flow suddenly increases near the slamming area, water splash and the large interface deforma­ tion can be observed. So, the reasonable way to assume values of the inter-particle physical quantity, such as velocity, density, sound speed and pressure, has a crucial influence on the stability and precision of numerical simulation. Rezavand et al. (2019) adopted the density-weighted inter-particle velocity and pressure to represent the inter-particle velocity and pressure, respectively. These formulas are reasonable for two-phase flows with small density ratios. However, for the two-phase flows with large density ratios, the introduction of the (density-weighted) average pressure or velocity will reduce the denser-fluid influence on the lighter fluid, which will obviously lead to great error (e.g. when ρa<<ρb, P¼(ρaPbþρbPa)/(ρaþρb)�Pa). Detailed discussion can be found in Chen et al. (2015). In the reference Grenier et al. (2013), it showed the considerable discrepancy between the results of the arithmetic mean viscosity and the harmonic mean viscosity in two-phase flow with large viscosity ratio (e.g. when μa<<μb, the vis­ cosity based on the arithmetic mean is μab ¼(μaþμb)/2�μb/2 and the viscosity based on the harmonic mean is μab ¼(2μaμb)/(μaþμb)�2μa). Therefore, the inter-particle properties based on the harmonic mean are used in our two-phase SPH model: 8 2uL uR > > > u ¼ ðu þ u Þ > > L R > > > > > 2ρi ci ρj cj > > � > ρij cij ¼ > < ρi ci þ ρj cj (11) > 2pi pj > > � p ¼ > > > pi þ pj > > > > > > 2ρi ρj > > � : ρij ¼ ρi þ ρj

P*ij

v

uR Þ;

ðuL

uR Þ > 0 ðuL

(13) uR Þ � 0

where,λ ¼ 5(2 þ d)/4 is the correction coefficient of viscous term, which is similar to the coefficient proposed by Monaghan and Gingold (1983). According to the numerical examples,and d is the dimension; μij ¼ 2μi μj =ðμi þμj Þ is the inter-particle dynamic viscosity,and μi is the dynamic viscosity of particle i. Substituting Eq. (13) in Eq. (12), the improved SPH model based on Riemann solver can be obtained. 3.3. A switch-function-based Riemann solver dissipation The water entry problems involve complicated flow phenomena, namely, overturning, breaking and merging of free surface, and vortex shedding. The traditional grid-based methods used in modeling these phenomenon have difficulties. When the effect of the air is considered, it is even more numerically difficult. Though the particle methods, such as SPH method, can make up the above defects, the numerical instability still remains. To deal with this instability caused by complex interface deformation and the strong impact, we propose to add a switchfunction-based Riemann solver dissipation into the intermediate pres­ sure in Eq. (13). The new intermediate pressure can be written as: 0 1 P*ij ¼ P*ij þ ηρij uðuL 2

uR Þ

(14)

where η is the switch function. If the dissipation term 12ρij uðuL uR Þ in Eq. (14) is used throughout the whole two-phase flow field, the numerical instability induced by the impact can be solved while the excessive dissipation is introduced into the model. The best way is to timely adjust the dissipation, the effec­ tiveness of this dissipation term would be opened in the strong impact domain and be closed in the weak shock parts. The Morris and Mon­ aghan (1997) introduced the switch function to adjust the artificial viscosity, the parameter of artificial viscosity varied with time and location, but it was not fully “switched off” at free position. Referring to the work of Morris and Monaghan (1997) and considering the consis­ tency of Riemann solver in our SPH model, we propose a new switch-function-based Riemann solver dissipation for water entry problems in the two-phase flow. The switch function η is written as: � 1; ϕ>χ η¼ (15) 0; ϕ�χ

After the Riemann solver is obtained in Eq. (10), P*ij and v*ij are used to replace the term (piþ pj)/2 and the term(viþ vj)/2 in Eq. (8) respectively. Through this transformation, the continuity and momentum equation in Eq. (8) can be written as � 8 dρ X� i > vi v*ij ⋅ri Wij Vj ¼ 2 ρi > > > < dt j (12) > X > dvi 2 > > P*ij ri Wij Vj þ g ¼ : dt ρi j where v*ij ¼ u*ij eij þ

8 < p þ λ μij ðu L h ¼ : p;

where ϕ ¼ rP⋅v as a color function can be interpreted as the pressure convection rate. χ is the threshold parameter which is determined ac­ cording to the impact strength of FSI problems. The formula of χ ¼ ρgvmax is used in this paper. If not specially noted, vmax is the maximum particle velocity in the calculation. Finally, combining with the improved SPH model based on Riemann solver and the switch-function-based Riemann solver dissipation, the new intermediate pressure in Eq. (14) is used to replace the pressure term in Eq. (12) and the novel two-phase SPH model based on an improved Riemann solver can be obtained.

� ueij , the harmonic mean continues to be

adopted in the average velocity between particle i and j, namely vij ¼ 2vi vj =ðvi þ vj Þ (it is different from the literature Rezavand et al. (2019)). If the model in Eq. (12) is used directly for simulating the two-phase flow, the numerical dissipation introduced by Riemann solver will in­ fluence the real motion characteristics of the two-phase flows (the dissipation term is 12ρij CðuL uR Þ). For the sake of reducing the numer­ ical dissipation, Zhang and Hu (2017) applied a limiter to replace the inter-particular sound speed and simulate single-phase flows. Later, Rezavand et al. (2019) used it to simulate the multi-phase flows. But when the fluid properties change, such as viscosity, the use of the limiter is restricted. In order to reduce the Riemann dissipation in Eq. (12) and simulate the fluid-structure coupling problems in two-phase flows with different viscosity ratios, the intermediate pressure is improved in the present model and shown as:

3.4. Algorithm for the fluid-body interaction The dummy boundary is adopted (Adami et al., 2012) in this paper. The boundary particles are treated as the fluid particles participating in the calculation of N–S governing equations for two-phase fluid to compensate the lack of particles in the support domain of fluid particles near the boundary. To prevent the fluid particles from penetrating the solid boundary, the reference density of the boundary particles will be the initial density of heavier fluid. The interaction between a fluid particle and a boundary particle is regarded as a one-sided Riemann 4

Q. Yang et al.

Ocean Engineering 199 (2020) 107039

problem, and the direction is along the boundary normal direction n (Wang et al., 2019). The initial states of the fluid particle i and the boundary particle j are written as:

8 > r ¼ r0 þ T � ðrwi0 r0 Þ > > wi > < vwi ¼ v0 þ θ � ðrwi r0 Þ > > dw > > : awi ¼ � ðrwi r0 Þ þ w � ½w � ðrwi dt

(16)

Ui ¼ ðρL ; uL ; PL Þ ¼ ðρi ; vi ⋅n; Pi Þ Uj ¼ ðρR ; uR ; PR Þ ¼ ρj ; vj ⋅n; PL þ ρi ðaw

gÞ⋅rij

�

(17)

cos θ þ e2x ð1

6 T¼6 4 ex ey ð1 ex ez ð1

3 cos θÞ

cos θÞ þ ez sin θ cos θÞ þ ey sin θ

ex ey ð1

cos θÞ þ ez sin θ

cos θ þ e2y ð1 ey ez ð1

cos θÞ

ex ez ð1

cos θÞ þ ey sin θ

7 cos θÞ þ ex sin θ 7 5 cos θ þ e2z ð1 cos θÞ

cos θÞ þ ex sin θ

where eθ ¼ [ex, ey, ez] and θ are the unit vector and modulus of the angular vector θ, respectively. After adding this new technical means, our model can simulate the multi-degree motion of body in water entry problems. In order to eliminate the non-physical pressure reflection from the solid wall, the sponge layer with a certain thickness is arranged inside the solid wall. The sponge layer is initially introduced by Gong et al. (2009),which is applied to the fluid particles. The forms are shown as � � ��dρ � � d ρi 50ξ i ¼ 1 100 0:9 (24) dt s dt

0

no-slip boundary, the viscous force is not ignored in the momentum equation, the intermediate pressure is obtained from the Eq. (13) and the right state velocity uR in Eq. (17) should be written as follows: (18)

uR ¼ vi ⋅n þ 2vj ⋅n

(23)

ey ez ð1

the fluid and boundary and the intermediate pressure is P*ij ¼ p. For the

The final pressure of the boundary particles is obtained through the Shepard interpolation from the neighboring fluid particles, the pressure at each time step is calculated as: � P P� PL þ ρj ðaw gÞ⋅rij Wij PR Wij j2fluid j2fluid P Pw ¼ P ¼ (19) Wij Wij j2fluid

r0 Þ �

where T is the transformation matrix. Note that, the particle velocity vwi is not used to update the particle position. The transformation matrix is given by:

where vj is the boundary velocity and aw is the boundary acceleration. For the free-slip condition, we ignore the viscous interaction between

2

(22)

� � � dvi ¼ 1 dt s

j2fluid

where w represents the boundary and f is the fluid. PL is the pressure of fluid particles around the boundary particle. The density of the bound­ ary particle can be obtained by solving the equation of state (Eq. (2)). The free-slip boundaries are implemented in the cases of this paper unless otherwise stated. In the water entry problems, the object is regarded as the moving boundary. Under the fluid inertia force and the gravity force, the motion of the object has two types: translation and rotation. Through the transformation of Eq. (6), the Eq. (4) can be expressed as: ! � 8 dv 1 X X 2mw * ’ 0 > þg P r W V ¼ > w wj j > M w2object j2fluid ρw wj > dt < (20) ! ! > > > 1 X X 2mw X * ’ > dw0 : rw r0 � P r W V ¼ I w2object j2fluid dt ρw j wj w wj j

100

0:950ξ

��dv � i dt

(25)

where ξ ¼ l/s, l is the perpendicular distance from particle i to solid wall, and s is the sponge layer thickness. 3.5. Time integration scheme To integrate the system of the governing equations of two-phase fluid and object, a velocity-Verlet scheme which is second-order accurate and reversible is adopted (Adami et al., 2012). To ensure the stability simulation, the time step must satisfy the CFL conditions: 8 � � h > > > Δt � 0:25min > > c þ jvmax j > > > > � 2� > < ρh Δt � 0:125min (26) μ > > > ffi ffi ffi ffi ffi s > > > > h > > > : Δt � 0:25min jgj

where r0 is the centre of the object’s mass M. The position of the object is updated by the velocity and the angular velocity of the object: 8 dr0 > > ¼ v0 < dt (21) > > : dθ ¼ w0 dt

where, vmax is the maximum particle velocity in the calculation. The minimum value of the above conditions is the numerical time step. 4. Numerical examples In this section, five cases of the fluid-structure coupling problems are simulated to demonstrate the stability and accuracy of the present twophase SPH model for two-dimensional simulation. Firstly, the two-phase dam break is presented, and the ability of our model to deal with the impact between fluids is tested by the switch-function-based Riemann solver dissipation. Then the sinking of a rectangle body is simulated to verify the accuracy for predicting multi-degree motion of sinking body in the viscous two-phase flow. Then the water entry of a wedge, flat plate slamming on water and water entry of a catamaran are simulated. With

where θ ¼ [θx, θy, θz] denotes the angular vector of rotations at new time step which is the sum of the object rotation angle with respect to the initial position. So, based on the rotation formula of coordinate axis, the position rwi, velocity vwi and the acceleration awi of each object particle can be written as (Leroy et al., 2014):

5

Q. Yang et al.

Ocean Engineering 199 (2020) 107039

Fig. 2. The initial setup of two-phase dam-break model.

Fig. 5. The time histories of surge-wave.

in two-phase flow. It was also modeled by the SPH method (Gong et al., 2016; Chen et al., 2015; Adami et al., 2012; Rezavand et al., 2019) and experimented in literature (Buchner, 2002). Although the air has no effect on the motion of the water, which has been proved by Colagrossi and Landrini (2003), the two-phase flow simulation of dam break can capture the formation and collapse of the closed air cavity which is an exact reproduction of the natural phenomena of dam-break. The initial setup of the two-phase dam break is shown in Fig. 2, which is the same as the experiment in Buchner (2002). The length of the dam break is lw ¼ 1.2 m and the height of the dam break is hw ¼ 0.6 m where g ¼ 9.81 m/s2. The length of the computational domain is Lt ¼ 5.366hw, the height is Ht ¼ 3hw, and the solid wall is simulated with free-slip condition at each wall. The water column of dam break is surrounded by the air. In the SPH simulation, the reference density of water is ρ0w ¼ 1000 kg/m3 and the viscosity is μw ¼ 0.001 Pa s, and the reference density of air is ρ0a ¼ 1 kg/m3 and the viscosity is μa ¼ 1.0 � 10 5 Pa s. The initial particle spacing Δx ¼ 0.005 m is used to discrete the computational domain, the artificial sound speed for water is cw ¼ 10 (ghw)0.5, and the time step is 1.0 � 10 6s. The physical time is non­ dimensionalized by the gravity and dam-break height. The background pressure is P0 ¼ 1000pa. The initial pressure distribution of dam break is also plotted as Fig. 2, which is derived by Greco (2001) and given in Eq. (27). The initial setup is to simulate the process of water released from a gate in the experiment.

Fig. 3. The evolution process of dam break, pressure distribution for water and velocity distribution for air.

the change of the structure shape, the influence of air on fluid-structure coupling problems becomes more and more apparent. It can test the ability of the present two-phase SPH model to deal with the complicated water-entry problems. These cases we select have been studied by many researchers adopted by experiment and numerical simulation. 4.1. The dam-break The dam-break case is a typical violent free surface flow, which in­ volves the phenomena of the flow impacting on the solid wall and overturning, breaking and merging of free surface. This case can validate the accuracy of the proposed two-phase model in predicting impact load, capturing two-phase flow interface and keeping the numerical stability

Pw ðx; yÞ ¼ ρ0w gðhw ∞ � 8ρ ghw X 0w 2

π

(a)

n

yÞ

� �� 2n þ 1 ð2nþ1Þπ x=ð2hw Þ e cos π y 2hw ð2n þ 1Þ2 1

(27)

(b)

Fig. 4. Two-phase flow distribution of dam break at t(g/hw)0.5 ¼ 6.17. (a) Without the switch-function-based Riemann solver dissipation, (b) with the switchfunction-based Riemann solver dissipation. 6

Q. Yang et al.

Ocean Engineering 199 (2020) 107039

wave front of dam break with four different particle resolution, namely Δx ¼ 0.04 m, Δx ¼ 0.02 m, Δx ¼ 0.01 m and Δx ¼ 0.005 m. The results from present model are compared with the experimental data (Buchner, 2002) and an analytical solution (Ritter, 1982). From the Fig. 5, we can see that the speed of front wave is slower in the experiment than that in the present simulation. The main reasons for these differences are attributed to the uncertain factors existing in experiment which include measuring accuracy of the experiment, the roughness of the water tank and so on (the free-slip condition is used in the present simulation). After t(g/hw)0.5 ¼ 1.7, the present results agree well with the results of Shallow water theory (Ritter, 1982). It can also be seen that the results from the present are convergence with the increasing particle spacing. Fig. 6 shows the time histories of a pressure probe located at 0.19hw from the bottom of the right wall. At t(g/hw)0.5 ¼ 2.41, the surge-wave front impacts the right wall in the Fig. 3, the pressure at the probe suddenly reaches the first peak. Then, the water moves along the rigid wall and the pressure remains roughly invariant. At t(g/hw)0.5 ¼ 6.15, the air cavity closes, the pressure reaches the maximum and a water jet is formed. Later, due to the closed air cavity, the pressure curve oscillates. Compared with the experiment, second pressure peak of the present results is slightly larger and the occurrence time is later than the experiment. It is also necessary to see that the present results in general agree well with the experimental results and the present two-phase SPH model is reliable in predicting the impact loads. For the study of the energy conservation of the present two-phase SPH model, the time histories of the total energy variation are shown in Fig. 7. The total energy variation ΔE is calculated by the method of the literature Cao et al. (2018). By contrasting the results of the present model without the switch-function-based Riemann solver dissipation and the literature (Cao et al., 2018), it can be observed that the curves of the energy loss show two significant turning point at the time of the water-wall impact and the water-water impact. For the same resolution (Δx ¼ 0.005), the present two-phase SPH model produces a better conservation than the model in literature (Cao et al., 2018). Without the switch-function-based Riemann solver dissipation, it produces a best conservation, while the simulation may not continue due to the nu­ merical instability caused by strong impact (it is shown in Fig. 4(a)).

Fig. 6. The time histories of impact pressure on the right wall.

Fig. 3 presents the pressure field distribution for water and the ve­ locity field distribution for air at six instants of dam-break evolution. It can predict the main flow characteristics of break-dam evolution. At t(g/ hw)0.5 ¼ 2.41, the water front reaches the right wall and impacts the vertical wall. Then the water climbs along the right wall and the back­ ward plunging jet forms at t(g/hw)0.5 ¼ 5.92. Subsequently, it shows that the air cavity region is enclosed under the action of gravity, then the water jet forms at t(g/hw)0.5 ¼ 6.15. The two-phase interface is captured clearly and the stable pressure field in water is simulated in the process of dam-break evolution. To visualize the effect of the new switch-function-based Riemann solver dissipation in such two-phase flow, we show the results of phase distribution of dam-break in Fig. 4. Without this strategy, the instable phenomena of particle penetration, particle clumping and unnatural void can be observed when the water-water impact occurs. With the strategy, the interface between two phases is captured clearly and smoothly, and the numerical calculation will continue steadily. Through the comparison, the effectiveness of the switch-function-based Riemann solver dissipation in suppressing the numerical instability is validated. In Fig. 5, for convergence study, we show the time histories of surge-

4.2. Sinking of a rectangle body To further validate the effectiveness of the new two-phase SPH method in predicting the translation and rotation of rigid, the sinking of

Fig. 7. The time histories of the total energy variation of the dam-break. 7

Q. Yang et al.

Ocean Engineering 199 (2020) 107039

P0 ¼ 1Pa, and the numerical sound speed of the heavier phase is c0 ¼ 10 (3gh)0.5. The free-slip boundary conditions are applied in the rectangle body and solid wall, and the sponge layer is not used in this case. To reveal the flow field structure of the sinking process of rectangle body, the pressure field distribution for water and the velocity field distribution for air at eight instants are shown in Fig. 9. Under the gravity and buoyancy, the rectangle body with the asymmetric distri­ bution of mass starts to sink and do rotational motion at t ¼ 0.75s. The heavier phase slams the side wall of the rectangle body and the lighter phase cavity forms at t ¼ 1.91s. The slamming load stops the further rotation of the body. Then, the lighter phase in the cavity is seen escaping from the upper right corner of the body, and a lighter phase cavity above the body forms at t ¼ 2.65s. With the continuous sinking, the lighter phase continuously escapes from the upper of the body. Finally, the body continues to sink almost without rotational motions (t > 3.68s). There is no pressure noise, and two-phase flow interface is smooth and clear in the sinking process. In Fig. 10, the time histories of the horizontal displacement and the vertical displacement of the rectangle body are presented together with the results of Finite Volume Particle (FVPM) without lighter phase from Barcaolo (2013). Four examples with different particle spacing, Δx ¼ 0.04 m, Δx ¼ 0.02 m, Δx ¼ 0.01 m, and Δx ¼ 0.005 m, are simulated. As the resolution increases from Δx ¼ 0.04 m to Δx ¼ 0.005 m, the results of the horizontal displacement curve and the vertical displacement curve are convergent and in good coincidence with the FVPM results. Comparing the results of two methods in the horizontal displacement, it can be found that the lighter phase has little effect on the horizontal motion of the rectangle body due to the influence of enclosed cavity of lighter phase. And, the present two-phase flow SPH model can well simulate such fluid-structure coupling problems and can be a new technical means for the study on multi-degree motion of body in two-phase flow.

Fig. 8. The initial configuration of the sinking process of a rectangle body.

a rectangle body with the asymmetric distribution of mass in the twophase flow is simulated (the setup is the same as that of Finite Volume Particle Method (FVPM) in Barcaolo (2013)). Early Sun et al., 2015 used this case to validate their single-fluid SPH model against the FVPM re­ sults (Barcaolo, 2013). The initial configuration of the system is shown in Fig. 8, the rectangle body has the length of l ¼ 1.0 m and the height of h ¼ 0.5l, which is placed in a 4l � 5l container. The depth of water is hw ¼ 3l. The mass center of the rectangle body is G¼(2.25l, 3l), the mass is M ¼ 1.0 kg/m (the mass of unit thickness is 1.0 kg in 2D numerical model), and the inertia moment is I ¼ 0.083 kg m2. The container contains two phases, the density of the heavier phase is set to ρw ¼ 1 kg/m3 and the kinematic viscosity is υ ¼ 2 � 10 3 m2/s, and the density of the other phase is ρa ¼ 1 � 10 3kg/m3 and the kinematic viscosity is υ ¼ 1 � 10 3 m2/s. The case is simulated in two-dimensional space and the rectangle body freely sinks into the heavier phase under the gravity in which the rectangle body has three degrees motion. The initial particle spacing is Δx ¼ 0.01 m, the time step is Δt ¼ 5 � 10 6s, the background pressure is

4.3. Water entry of a wedge Water entry problems are very significant in ocean engineering fields. In this section, the proposed two-phase SPH model is used to simulate the wedge entry into the water, the experiment of which is conducted by Zhao et al. (1996). This case not only validates the accu­ racy of the numerical calculation model in predicting the fluid force subjected on the structure and the motion characteristics of the struc­ ture, but also tests the model’s ability in capturing the flow character­ istics of the lighter phase and in suppressing instability due to the strong impact in water entry problems. As shown in Fig. 11, the total mass of the free-falling wedge with D ¼

Fig. 9. Snapshots of the sinking process of the rectangle body (pressure distribution for water and velocity distribution for air). 8

Q. Yang et al.

Ocean Engineering 199 (2020) 107039

(a)

(b)

Fig. 10. Comparison of the present two-phase SPH predictions against FVPM results (Barcaolo, 2013): (a) the horizontal displacement, and (b) the vertical displacement of body’s center of mass.

velocity field of air. It can be seen that the velocity field of air keeps smooth and stable in the free-falling process of wedge, and air vortex field appears in the rear of wedge due to the pressure reducing. When the wedge is approaching the water surface, the air between wedge and water is seen escaping at high velocity. At t ¼ 0.63s, the wedge reaches the water, which brings high pressure region in the water. Soon the water jet forms at t ¼ 0.64s. To further demonstrate the stability of our SPH model in long-time numerical simulation of water entry problems, it is simulated that the wedge further moves down after slamming water surface. It can be observed that the air vortex generates outside the splashing jets at t ¼ 0.9s. A clear interface between water, air and solid wall, and a stable pressure field are observed during the whole process of wedge entry. In order to further illuminate the effect of the new switchfunction-based Riemann solver dissipation on such strong impact problems, Fig. 13 shows the snapshots of wedge entry at t ¼ 0.67s adding this treatment technology or not. Without this treatment tech­ nology, the unstable interface, particle clumping and unnatural void regions can be observed as shown in Fig. 13(a). The time histories of the fluid force acting on the wedge and the falling velocity are shown in the Fig. 14. Comparing with the experi­ mental data of Zhao et al. (1996), a good agreement can be observed. Minor deviations in the falling velocity can be seen, the falling velocity of the present SPH results goes down more quickly than the experiment, which is likely due to the three-dimensional effect in experiment. It is also noted that there is almost no difference in the contact force and falling velocity between SPH results with air and no-air. The effect of the air on the water entry of wedge can be ignored. The present two-phase SPH model can predict the impacting loads and motion characteristics of structure very well in simulating such fluid-structure coupling problems. In order to study the convergence of the present method in water entry problems, two cases of water entry of a wedge are carried out (Δx ¼ 0.01 m and Δx ¼ 0.02 m). Fig. 15 shows the comparisons of the simu­ lation results. The present method achieves convergence with increasing particle resolution.

Fig. 11. The initial setup of water-entry of wedge.

0.5 m width and 1.0 m length is 241 kg in the experiment (the mass of unit thickness is 241 kg in 2D numerical model), the deadrise angle of the wedge is 30� and the entry-water velocity is 6.15 m/s. The numerical water tank is b ¼ 7D wide and h ¼ 3D deep. In order to improve the computational accuracy and capture the detailed flowing behaviors, the particle spacing Δx ¼ 0.005 m is used (about 0.65 million particles are used in the numerical study) and the time step is Δt ¼ 2 � 10 6s. The density and viscosity of the water are ρw ¼ 1000 kg/m3 and μw ¼ 0.001 Pa s, and that of air are taken as ρa ¼ 1 kg/m3 and μa ¼ 1.0 � 10 5 Pa s respectively. The background pressure is 1500Pa. Besides, the sponge layer is adopted in the solid boundary wall of the water tank to reduce the effect of pressure reflection on the flow field pressure and impact load on wedge. Fig. 12 presents the snapshots of the wedge-entry process. In order to better observe the stability of the present model in simulating multi­ phase flows, each snapshot presents the pressure field of water and the

4.4. Flat plate slamming on water In the present section, the influence of the air on water entry of the flat plate is investigated. Past research (Lind et al., 2015) has demon­ strated that the air has little effect on low-speed water entry problems in 9

Q. Yang et al.

Ocean Engineering 199 (2020) 107039

Fig. 12. The snapshots of the water entry wedge (pressure distribution for water and velocity distribution for air).

(a)

(b)

Fig. 13. The particle distribution at t ¼ 0.68s. (a) Without the switch-function-based Riemann solver dissipation, (b) with the switch-function-based Riemann solver dissipation.

(a)

(b)

Fig. 14. Time histories of the contact force and falling velocity of wedge during the impacting process.

the water-entry early stages, especially the structures with the deep dead-rise angle. For the water-entry plate with a low dead-rise angle, the air trapped underneath has a major effect on the structure subjected to water impact load. The plate we simulate is the same as that in the experimental study is

carried out by Ma et al. (2016). In the experiment, the mass of plate is 32 kg, the size is 0.25 m � 0.25 m and the impact velocity of the free-falling plate is 5.5 m/s. The width of water entry plate is 0.25 m in 2D simu­ lation by our present two-phase SPH model. In order to reduce the large amount of computation, the initial distance from the free-falling plate to 10

Q. Yang et al.

Ocean Engineering 199 (2020) 107039

Before t ¼ 0.12s, the falling plate has no influence on the water surface. The air velocity behind plate increases as the plate moves down. Air vortex zones forms at the both sides of plate, where the velocity is much less than neighboring region. At t ¼ 0.12s, the plate does not touch the water surface, the air between the plate and water is seen escaping quickly from the both sides of plate, and the air is compressed and high pressure zone occurs in the water. When t ¼ 0.13s, the plate slams the water surface and water pressure increases dramatically. There is still a small amount of air that could not have time to escape and is entrapped beneath the plate. At t ¼ 0.132s, the high pressure region of water en­ larges and water jets form because of plate slamming on water. At t ¼ 0.14s, the high pressure region propagates far away and the jets continue to develop. Due to the sponge layer implemented inside the solid wall, non-physical pressure reflection phenomenon could not be seen. Considering the air effect, Fig. 18 shows the time histories of the contact force on the plate from the present single-phase SPH model, twophase SPH model and experiment during the water entry process. The time t ¼ 0 is the time of contact-force peak. It can be seen that the force peak and the changing gradient of the contact-force curve in present single-phase model are much larger than the others. The main reason is that air is not considered in the single-phase model and air cushion effect that can reduce slamming load cannot be shown. In the respect of force peak, the changing gradient of contact force and the duration of impact, the two-phase SPH results are in general agreement with the experi­ mental data (Ma et al., 2016), and the minor deviation of the contact-force peak can be seen due to the influence of three-dimensional effect. It also illustrates that the air can cushion the impact load of structure with a low dead-rise angle in the water-entry process.

Fig. 15. The time histories of the contact force for the different parti­ cle resolution.

water surface is 0.6 m and the initial velocity is 4.3 m/s. This treatment has no effect on the impact phenomenon. The sketch of the simulation is shown in Fig. 16, the computational domain is 1.2 m � 1.5 m and the water depth is h ¼ 0.8 m. The physical properties of water and air are the same as the section 4.3. The particle spacing is Δx ¼ 0.002 m, the nu­ merical sound speed of water is cw ¼ 60 m/s, the background pressure is P0 ¼ 2000Pa and the time step is 1.5 � 10 5s. To investigate the dynamics characteristic of the two-phase flow in the water-entry plate, the Fig. 17 shows several snapshots of the computed process about water pressure field and air velocity field.

Fig. 16. The sketch of the water entry plate. 11

Q. Yang et al.

Ocean Engineering 199 (2020) 107039

Fig. 17. The snapshots of the water entry plate (pressure distribution for water and velocity distribution for air).

Fig. 18. Time histories of the contact force acting on the wedge.

4.5. Water entry of a catamaran

catamaran cannot escape easily in our 2D simulation (but the air can escape from the fore and aft of the 3D model in the Whelan (2004) experiment), the cabin effect will be more obvious when the catamaran touches the water surface. The 2D catamaran model is shown in Fig. 19. In the Fig. 19(b), the air model 1 is the same as that used by Shahraki et al. (2011), the width of this mode is 0.56 m, the height of side wall is 0.19 m, depth of the in­ ward valley is 0.13 m, the dead-rise angle of the intermediate wedge is 25� , and the mass is 72.25 kg/m (the mass of unit thickness is 72.25 kg in

One of the most challenging water entry problems in the two-phase flow is the water entry of a catamaran, which involves many dynamic problems such as the two-phase flow coupling with the complex struc­ ture, the complex air flow, the distribution of flow field structure, the impact loads, the dynamic response of the structure and so on. In this section, the water entry catamaran with a 25� dead-rise angle is simulated. When the catamaran impacts the water, the air beneath the 12

Q. Yang et al.

Ocean Engineering 199 (2020) 107039

(a)

(b)

(c)

Fig. 19. The model of water-entry catamaran: (a) The initial setup of the model; (b) The catamaran dimensions of air model 1; (c) The catamaran dimensions of air model 2.

m, the time step is 1.0 � 10 5s, the sound speed of water is 15 m/s, and the background pressure is 1000 Pa. The time evolutions of the catamaran acceleration are shown in Fig. 20 where the results of no-air and air SPH models are compared to the experiment data (Whelan, 2004). It is clear that the results of the present SPH solution without air have two peaks: the first peak at t ¼ 61.5 ms is due to the impact on the intermediate wedge; the second peak at t ¼ 88.5 is due to the water impact on the corner of the center wedge meeting the side panels (the snapshots are shown in Fig. 21). However, the SPH acceleration of two peaks is higher than the acceleration recorded in the experiment and the experimental time histories are much smoother and gentler. It can be concluded that the air cushion effects cannot be ignored in the simulation of a water-entry catamaran. As in the air model 1, large differences appear for the time histories of acceleration and the peak time relative to the experimental measure­ ments. The reason for these discrepancies is that the experiment builds the three degree model of water-entry catamaran, the air beneath the catamaran could escape easily from the fore and aft of the 3D model. However, in the air model 1, most air is trapped under the catamaran and cannot escape away. At the peak time t ¼ 62.5 ms, the snapshot is reported in Fig. 22, the center wedge does not touch the water surface due to the entrapped air and the acceleration peak occurs earlier than that in the experiment (the same conclusion can be seen in Yan et al. (2015)). The smooth and uniform pressure field across the interface can be observed in Fig. 22(a). In order to reduce the excessive influence of the air cushion effects in 2D numerical model, the air model 2 with shorter side wall is simulated. At t ¼ 82.5 ms, the snapshot of the peak acceleration of air model 2 is shown in Fig. 23, it can be seen that there is only a small amount of air trapped in the catamaran. From Fig. 20 it is possible to observe that the changing gradient and peak time of the acceleration-time histories are much closer to the experimental data than that in the air model 1.

Fig. 20. Time histories of the catamaran acceleration.

2D numerical model). For the study of the effect of air beneath the catamaran, the air model 2 with short side wall in Fig. 19(c) is added and simulated (in addition to the height of side wall, other parameters are the same as the air model 1). The computational domain is 1.8 m � 1.8 m, and the depth of water is 1 m which is the same as the experiment conducted by Whelan (2004). In order to study the influence of air on the water entry catamaran, the impact velocity of the catamaran in all models is 1.03 m/s. For the air model, the initial velocity and initial distance between the inter­ mediate wedge bottom and water surface are 0.85 m/s and 0.082 m respectively, the values of the no-air model are 0.97 m/s and 0.03 m. The same as that in Yan et al. (2015). We adopt the present SPH model to simulate this water entry problem. The initial particle spacing is 0.003

(a )

(b)

Fig. 21. The snapshots of the peak acceleration without air model: (a) The first peak acceleration at t ¼ 61.5 ms; (b) the second peak acceleration at t ¼ 88.5 ms. 13

Q. Yang et al.

Ocean Engineering 199 (2020) 107039

(a)

(b)

Fig. 22. The snapshots of the peak acceleration with air model at t ¼ 62.5 ms: (a) Pressure distribution for water and air; (b) pressure distribution for water and velocity distribution for air.

reproduced by the no-air model. In the air model 1 and 2, the air beneath the catamaran cannot escape in time and is trapped under the catamaran and a smooth and gentle rise is on the pressure curve. The P4 peak pressure is not caused by the direct impact of water, but by the trapped air, which occurs earlier than those in the no-air model and experiment. Comparison between the results of air model 1 and 2, reducing the height of side wall can decrease the influence of air cushion effect of 2D model and delay the appearance of the peak pressure, and the results of 2D air model 2 are closer to the 3D experimental data than the air model 1. Through the numerical simulation of the water entry of a 2D cata­ maran, it is proved that the present two-phase SPH model is suitable for the study of fluid-structure coupling problems considering air cushion effect in engineering applications. Because of the limit of computational scale, this paper only analyzed the 2D model of water-entry catamaran. By extending the 2D two-phase SPH model presented in this paper to 3D, the whole dynamic process of catamaran entry can be reproduced.

Fig. 23. The snapshot of the peak acceleration of air model 2 at t ¼ 82.5 ms.

5. Conclusion In the fluid-structure coupling problems, the effect of the air on the dynamic response of the water-entry rigid structure is significant. In this paper, a new two-phase SPH model based on an improved Riemann solver is presented for simulating the water entry problems of the twophase flow with large density ratios coupling with the structure. Four improved measures are made in this paper. An intermediate pressure of a Riemann form is developed and implemented in the dis­ cretized momentum equation, and this term can decrease the Riemann dissipation and simulate the two-phase flows with different viscosity ratios. The inter-particle properties based on the harmonic mean are adopted to weaken the influence of the interfacial instability and in­ crease precision of numerical computation. Considering the free-slip or no-slip condition at the wall, one-sided Riemann problem is used to impose the fluid-body interaction in two-phase flow, which could guarantee the stability of the intersection region between the solid wall and the two-phase flow interface. A switch-function-based Riemann solver dissipation is introduced in the momentum equation. Such strategy could improve the instability of interface owing to the strong impact. Five fluid-structure coupling examples are simulated by this new two-phase SPH model. Through rigorous comparisons with the existing numerical results and experimental data, it is demonstrated that the present model can give accurate dynamic-response solutions of structure and stable two-phase flow field distribution. It is also found that the effect of air on the water entry structure with a deep dead-rise angle is great, and the cushion effect can reduce the impact loads and protect the

Fig. 24. Time histories of the catamaran pressure at P4.

To further study the dynamic response of catamaran entry, the time histories of the pressure at P4 are shown in Fig. 24 (the location of point P4 is shown in Fig. 19), the results are compared with the experimental data (Whelan, 2004). By comparing the Fig. 20 and the Fig. 24, it is possible to observe that the P4 pressure peak and the acceleration peak appear almost at the same time in our 2D model. In the experiment, the air can escape from the fore and aft of the 3D catamaran model and the peak pressure is caused by the direct impact of water which is similar to the no-air model. The results of the experiment show that the pressure increases slowly in early stage due to the influence of the air cushion effect, and a sharp rise with an inflexion is on the pressure curve and a maximum value occurs. Although the pressure-time curve of our no-air model agrees well with the experimental results (Whelan, 2004), it can be seen from the impact time that the cushioning effect of air cannot be 14

Ocean Engineering 199 (2020) 107039

Q. Yang et al.

water entry structure from failure. Further work will be dedicated to studying the air cushion effect on the high-speed impact problems by the two-phase SPH model described in this paper.

Lind, S.J., Stansby, P.K., Rogers, B.D., Lloyd, P.M., 2015. Numerical predictions of waterair wave slam using incompressible-compressible smoothed particle hydrodynamics. Appl. Ocean Res. 49, 57–71. Lin, P.Z., 2007. A fixed-grid model for simulation of a moving body in free surface flows. Comput. Fluids 36 (3), 549–561. Liu, G.R., Liu, M.B., 2003. Smoothed Particle Hydrodynamics: a Meshfree Particle Method. Word Scientific, Singapore. Liu, M.B., Zhang, Z.L., 2019. Smoothed particle hydrodynamics (SPH) for modeling fluidstructure interactions. Sci. China (Phys. Mech. Astron.) 62 (8), 5–42. Ma, Z.H., Causon, D.M., Qian, L., Mai, T., 2016. Pure and aerated water entry of a flat plate. Phys. Fluids 28(1) (16104). Marrone, S., Colagrossi, A., Chiron, L., De Leffe, M., Le Touz� e, D., 2018. High-speed water impacts of flat plates in different ditching configuration through a RiemannALE SPH model. J. Hydrodyn. 30 (1), 38–48. Marrone, S., Colagrossi, A., Park, J.S., Campana, E.F., 2017. Challenges on the numerical prediction of slamming loads on LNG tank insulation panels. Ocean. Eng. 141, 512–530. Maruzewski, P., Le Touz� e, D., Oger, G., Avellan, F., 2010. SPH high-performance computing simulations of rigid solids impacting the free-surface of water. J. Hydraul. Res. 48 (Suppl. 1), 126–134. Monaghan, J.J., 1985. Particle methods for hydrodynamics. Comput. Phys. Rep. 3, 71–124. Monaghan, J.J., 1994. Simulating free surface flows with SPH. J. Comput. Phys. 110, 399–406. Monaghan, J.J., Gingold, R., 1983. Shock simulation by the particle method SPH. J. Comput. Phys. 52, 374–389. Monaghan, J.J., Kocharyan, A., 1995. SPH simulation of multi-phase flow. Comput. Phys. Commun. 87, 225–235. Morris, J.P., Fox, P.J., Zhu, Y., 1997. Modeling low Reynolds number incompressible flows using SPH. J. Comput. Phys. 136, 214–226. Morris, J.P., Monaghan, J.J., 1997. A switch to reduce SPH viscosity. J. Comput. Phys. 136, 41–50. Oger, G., Doring, M., Alessandrini, B., Ferrant, P., 2006. Two-dimensional SPH simulations of wedge water entries. J. Comput. Phys. 213, 803–822. Puri, K., Ramachandan, P., 2014. Approximate Riemann solvers for the Godunov SPH (GSPH). J. Comput. Phys. 270 (8), 432–458. Rezavand, M., Zhang, C., Hu, X.Y., 2019. A weakly compressible SPH method for violent multi-phase flows with high density ratio. J. Comput. Phys. 109092. Ritter, A., 1982. Die fortpflanzung de Wasserwellen. Zeitschrift Verein Deutscher Ingenieure 36 (33), 947–954. Shahraki, J.R., Penesis, I., Thomas, G., Davis, M.R., Whelan, J.R., 2011. Prediction of slamming behaviour of monohull and multihull forms using smoothed particle hydrodymamics. In: HSMV 2011 Conference Proceedings, May 25–27, 2011, Naples, Italy, pp. 1–10. Sirotkin, F.V., Yoh, J.J., 2013. A Smoothed Particle Hydrodynamics method with approximate Riemann solvers for simulation of strong explosions. Comput. Fluids 88, 418–429. Skillen, A., Lind, S., Stansby, P.K., Rogers, B.D., 2013. Incompressible smoothed particle hydrodynamics (SPH) with reduced temporal noise and generalised Fickian smoothing applied to body-water slam and efficient wave-body interaction. Comput. Methods Appl. Mech. Eng. 265 (9), 163–173. Sun, P.N., Ming, F.R., Zhang, A.M., 2015. Numerical simulation of interactions between free surface and rigid body using a robust SPH method. Ocean. Eng. 98, 32–49. Sun, P.N., Zhang, A.M., Marrone, S., Ming, F.R., 2018. An accurate and efficient SPH modeling of the water entry of circular cylinders. Appl. Ocean Res. 72, 60–75. Teschner, T.R., K€ on€ ozsy, L., Jenkins, K.W., 2019. A generalised and low-dissipative multi-directional characteristics-based scheme with inclusion of the local Riemann problem investigating incompressible flows without free-surfaces. Comput. Phys. Commun. 239, 283–310. Wang, J., Lugni, C., Faltinsen, O.M., 2015. Experimental and numerical investigation of a freefall wedge vertically entering the water surface. Appl. Ocean Res. 51, 181–203. Wang, L., Xu, F., Yang, Y., 2019. SPH scheme for simulating the water entry of an elastomer. Ocean. Eng. 178, 233–245. Wendland, H., 1995. Piecewise polynomial, positive definite and compactly supported radial functions of minimal degree. Adv. Comput. Math. 4 (1), 389–396. Whelan, J.R., 2004. Wetdeck Slamming of High-Speed Catamarans with a Centrebow. PhD Thesis. University of Tasmania. Wick, T., 2011. Fluid-structure interactions using different mesh motion techniques. Comput. Struct. 89 (13–14), 1456–1467. Yan, R., Monaghan, J.J., Valizadeh, A., et al., 2015. The effect of air on solid body impact with water in two dimensions. J. Fluid Struct. 59, 146–164. Zhang, A.M., Sun, P.N., Ming, F.R., 2015. An SPH modeling of bubble rising and coalescing in three dimensions. Comput. Methods Appl. Mech. Eng. 294, 189–209. Yang, Q.Z., Xu, F., Yang, Y., Wang, L., 2020. A multi-phase SPH model based on Riemann solvers for simulation of jet breakup. Eng. Anal. Bound. Elem. 111, 134–147. Zhang, C., Hu, X.Y., 2017. A weakly compressible SPH method based on a lowdissipation Riemann solver. J. Comput. Phys. 335, 605–620. Zhao, Ran, Faltinsen, O., Aarsnes, J., 1996. Water entry of arbitrary two-dimensional sections with and without flow separation. In: Proceedings of the 21st Symposium on Naval Hydrodynamics. National Academy Press, Washington, DC, USA, pp. 408–423. Trondheim, Norway.

Declaration of competing interest The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. CRediT authorship contribution statement Qiuzu Yang: Investigation, Data curation, Writing - original draft, Visualization. Fei Xu: Formal analysis, Resources, Supervision, Funding acquisition. Yang Yang: Software, Writing - review & editing. Jingyu Wang: Methodology, Visualization, Investigation. Acknowledgements This work was supported by the National Natural Science Foundation of China [Grant no. 11972309], the National Natural Science Founda­ tion for Young Scientists of China [Grant no. 11702220], the Funda­ mental Research Funds for the Central Universities [Grant no. 310201901A012],and Overseas Expertise Introduction Project for Discipline Innovation (the 111 Project) [Grant no. BP0719007). References Adami, S., Hu, X.Y., Adams, N.A., 2012. A generalized wall boundary condition for smoothed particle hydrodynamics. J. Comput. Phys. 231 (21), 7057–7075. Barcaolo, D.A., 2013. Improvement of the precision and the efficiency of the SPH method: theoretical and numerical study. PhD Thesis, Ecole Centrale de Nantes. Buchner, B., 2002. Green Water on Ship-type Offshore Structures. Ph.D. thesis. Delft University of Technology. Cao, X.Y., Ming, F.R., Zhang, A.M., Tao, T., 2018. Multi-phase SPH modelling of air effect on the dynamic flooding of a damaged cabin. Comput. Fluids 163, 7–19. Cha, S.H., Whitworth, A.P., 2003. Implementations and tests of Godunov-type particle hydrodynamics. Mon. Not. Roy. Astron. Soc. 340 (1), 73–90. Chen, Z., Zong, Z., Liu, M.B., Zou, L., Li, H.T., Shu, C., 2015. An SPH model for multiphase flows with complex interfaces and large density differences. J. Comput. Phys. 283, 169–188. Colagrossi, A., Landrini, M., 2003. Numerical simulation of interfacial flows by smoothed particle hydrodynamics. J. Comput. Phys. 191, 448–475. Gong, K., Liu, H., Wang, B.L., 2009. Water entry of a wedge based on SPH model with an improved boundary treatment. J. Hydrodyn. Ser. B 21 (6), 750–757. Gong, K., Shao, S.D., Liu, H., Wang, B.L., Tan, S.K., 2016. Two-phase SPH simulation of fluid-structure interactions. J. Fluid Struct. 65, 155–179. Greco, M., 2001. A Two-Dimensional Study of Green-Water Loading. Ph.D. thesis. Norwegian University of Science and Technology, Trondheim, Norway. Grenier, N., Antuono, M., Colagrossi, A., Le Touz�e, D., Alessandrinia, B., 2009. An Hamiltonian interface SPH formulation for multi-fluid and free surface flows. J. Comput. Phys. 228 (22), 8380–8393. Grenier, N., Le Touz�e, D., Colagrossi, A., Antuono, & M., Colicchio, G., 2013. Viscous bubbly flows simulation with an interface SPH model. Ocean. Eng. 69, 88–102. Hu, X.Y., Adams, N.A., 2006. A multiphase SPH method for macroscopic and mesoscopic flows. J. Comput. Phys. 213, 844–861. Hu, X.Y., Adams, N.A., 2007. An incompressible multi-phase SPH method. J. Comput. Phys. 227, 264–278. Inutsuka, S., 1994. Godunov-type SPH. Memor. Soc. Astronom. Ital. 65, 1027–1031. Khayyer, A., Gotoh, H., 2013. Enhancement of performance and stability of MPS meshfree particle method for multiphase flows characterized by high density ratios. J. Comput. Phys. 242, 211–233. Khayyer, A., Gotoh, H., 2016. A multiphase compressible–incompressible particle method for water slamming. Int. J. Offshore Polar Eng. 26 (1), 20–25. Khayyer, A., Gotoh, H., Shimizu, Y., 2019. A projection-based particle method with optimized particle shifting for multiphase flows with large density ratios and discontinuous density fields. Comput. Fluids 179, 356–371. Koshizuka, S., Oka, Y., 1996. Moving particle semi-implicit method for fragmentation of incompressible fluid. Nucl. Sci. Eng. 123, 421–434. Leroy, A., Violeau, D., Ferrand, M., Kassiotis, C., 2014. Unified semi-analytical wall boundary conditions applied to 2-D incompressible SPH. J. Comput. Phys. 261, 106–129.

15