Numerical simulation of water entry problems considering air effect using a multiphase Riemann-SPH model

The water entry is a classic fluid-structure interaction problem in ocean engineering. The prediction of impact loads on structure during the water entry is critical to some engineering applications. In this paper, a multiphase Riemann-SPH model is developed to investigate water entry problems. In this model, a special treatment, a cut-off value for the particle density, is arranged to avoid the occurrence of negative pressure. A remarkable advantage of the present multiphase SPH model is that the real speed of sound in air can be allowed when simulating water-air flows. In the present work, considering the air effect, several typical water entry problems are studied, and the evolution of multiphase interface, the motion characteristic of structure and complex fluid-structure interactions during the water entry are analyzed. Compared with the experimental data, the present multiphase SPH model can obtain satisfactory results, and it can be considered as a reliable tool in reproducing some fluid-structure interaction problems.


Introduction
The study on water entry problems has been a long-standing subject in ocean engineering. The prediction of the flow fields and the structure dynamic responses is of great use for many applications, such as the ship launching and the ditching of aircrafts [1]. The process of water entry involves some complicated physical problems including the evolution of multiphase interface and the fluid-structure interaction [2]. These problems bring some difficulties to traditional mesh-based methods. For instance, to capture the multiphase interface and to track the moving structure, some numerical treatments need to be added [3]. However, meshfree methods are advantageous in handling these difficulties. As a popular variant of meshfree methods, Smoothed Particle Hydrodynamics (SPH) method [4] has been applied to model numerous hydraulic flows (see e.g. [5][6][7][8][9][10]).
For water entry problems, some investigations have been carried out by SPH practitioners. Oger et al. [11] first applied a weakly compressible SPH (WCSPH) method to simulate wedge water entries. In their work, a new technic based on a particle sampling method is developed to evaluate the pressure on solid boundaries, and a variable smoothing length technique is proposed to improve the efficiency of the simulations. Skillen et al. [12] used an incompressible SPH method (ISPH) with modified shifting technology to study the water entry of a wedge and a cylinder. The authors developed a diffusion-based shifting to model the water entry problems for stability. To reduce temporal noises, the effect of particles near the free surface was introduced smoothly into the pressure Poisson equation. Sun et al. [13,14] conducted the 2D and 3D simulations of water entry problems using SPH method. The authors proposed a new SPH version named δplus-SPH model where a shifting technology suitable for WCSPH was developed and an adaptive particle refinement method was applied to save the computational overhead. Marrone et al. [15] adopted a Riemann-ALE-SPH method to simulate high-speed water entry of plates with different ditching configurations. They considered three main parameters, the horizontal velocity, the pitch angle and the approach angle during these water entries. Recently, Cheng et al. [16] analyzed the slamming loads in the process of ship hulls dropping into the water with the help of SPH method. The change characteristics of slamming loads and the evolution of free surface in the water entry process of a bow-flare ship section under different roll angles were analyzed.
However, most of studies are performed by the single-phase SPH model which neglects the effect of air on results. Due to the existence of the large density ratio between water and air, the multiphase simulations on water entries are hard to implement. Therefore, there are few related published literatures considering air effect on the water entries. Gong et al. [17] performed a water entry experiment of a wedge and employed a twophase SPH model to reproduce the entire wedge entry process. The wedge cavity and the air flow inside obtained by the SPH method show good agreements with the experimental results. Marrone et al. [18] conducted a comprehensive numerical study on predicting slamming loads of a LNG tank insulation panel with a deadrise angle of 4 • using the Riemann-ALE-SPH method. For this case, 3D effects and air-cushioning play relevant roles because of its small deadrise angle. Recently, Yang et al. [3] applied a new multiphase SPH model based on an improved Riemann solver to simulate some water entries. They proposed a switch-function-based Riemann solver dissipation to improve the instability of multiphase interfaces, and one-sided Riemann problem is considered to deal with the fluid-structure interactions. Although good numerical results can be obtained, for the sake of numerical stability, the SPH practitioners mentioned above did not consider the real speed of sound in air when simulating water entries. This paper is dedicated to providing a study on water entries with air flow based on a multiphase Riemann-SPH model proposed by Meng et al. [19], which is our previous work. In this model, the standard SPH governing equations are recast with the introduction of a Riemann problem. Therefore, the interaction between particle pair is determined by the solution to a Riemann problem. In the present work, this problem is solved by Roe's approximate Riemann solver with a designed dissipation limiter. To better simulate these water entries, the multiphase SPH model is further developed. In order to avoid the occurrence of negative pressure, a cut-off value for the particle density is set here. Due to the good robustness of this multiphase model, the real speed of sound in air can be considered. In the present work, the water entry processes of two cylinders, a wedge and a vessel section are investigated. The evolution of multiphase interface, the motion characteristic of structure and complex fluid-structure interactions are focused. Numerical results show that this multiphase model proves to be accurate and robust for predicting the motion of structure and impact loads during the water entry. Another motivation of this work is to test the performance of the present multiphase SPH model in simulating the whole process of a rigid body into static water.
This paper is arranged as follows. In Section 2, the multiphase Riemann-SPH model is briefly reviewed and developed for water entry problems. Section 3 discusses some numerical results of water entries. Section 4 draws some conclusions.

Governing equations
The multiphase Riemann-SPH model proposed by Meng et al. [19] is constructed combining SPH concepts with a Riemann solver. In this model, a one-dimensional Riemann problem is embedded into the standard SPH formulation. Specifically, every particle pair b and a carrying the left and right states can construct a one-dimensional Riemann problem where the discontinuity is defined at the middle point between this particle pair. Through the work of Meng et al. [19], the governing equations of this model are reformulated as where ρ a and v a denote the density and the velocity associated with the a−th particle, respectively. g is the gravity acceleration and it is set to be 9.81 m/s 2 pointing downward throughout this work. V b is the volume of particle b. ∇ a W ab is the gradient of the kernel function W (x ab , h) where x ab = x a − x b . In the present work, the improved Gaussian kernel function [20] with h x = 1.2 is used. h denotes the smoothing length and x the initial particle spacing.
For the water entries in the present work, the water-air interface is a special Riemann problem i.e., the contact discontinuity where the pressure and velocity are continuous while the density is discontinuous. Therefore, by means of the idea of the Riemann-SPH method, the multiphase interface problem is solved straightforwardly in theory. In the present work, for Riemann variables, the intermediate velocity is Here, the intermediate variable u * is obtained by Roe's approximate Riemann solver, and the solution is where u a = v a ·x ab and u b = v b ·x ab . C ab is the Roe-averaged Lagrangian speed of sound which can be obtained by c a and c b denote the speed of sound of particles a and b, respectively. In this model, the intermediate pressure p * is determined by where p a is the pressure of particle a. φ denotes a dissipation limiter. This limiter is designed to control the numerical dissipations introduced by this Riemann solver (More details can be found in [19]). φ has the following form whereρ stands for the harmonic average of ρ a and ρ b , i.e.ρ = 2ρ a ρ b . λ builds a relation with the artificial viscosity α, In this paper, the artificial viscosity α is set to be 0.1 if not specifically mentioned.
To solve the governing equations, the following equation of state [21] is adopted, where ρ 0 and c 0 denote the initial density and the artificial speed of sound, respectively. γ = 7 for water and γ = 1.4 for air. p bg is the background pressure which is used to avoid the negative pressure and to obtain uniform particle distribution [22]. It is set to be 100 Pa if not specially mentioned in this paper. For water-air flows, under the weakly compressible assumption, the artificial speed of sound in water c w is determined by v 0 and H are the initial velocity of water entry and the initial water depth, respectively. For most of SPH practitioners (see e.g. [23]), the artificial speed of sound in air c g is determined by where the subscripts w and g stand for the water and air, respectively. This constraint is proposed by Colagrossi and Landrini [23] for numerical stability. It can be found that, under this constraint, for water-air flows (i.e., ρ w /ρ g = 1000), the actual speed of sound in air can not be allowed when c w > 25 m/s. Hence, the actual compressibility of air can not be considered for some flows, especially for water entries. In our previous work (see [19]), the multiphase Riemann-SPH model adopted here has been demonstrated that it can allow the actual speed of sound in air when modelling multiphase flows. In this work, the artificial speed of sound in air is 340 m/s if not specifically noted. Because of the use of the real speed of sound in air, the time step of the proposed method is larger than that of other multiphase SPH models for some water entry cases when c w > 25 m/s. Therefore, the computational overhead can be saved.

Treatment on the particle density
Negative pressure is usually encountered in SPH for modelling multiphase flows, especially for long-term simulations. The existence of negative pressure does not accord with real physics, and it may lead to the occurrence of the tensile instability, even the termination of the simulation. Therefore, similar to Chen et al. [24], the particle density is limited to a cut-off value in this multiphase model. Based on the equation of state, to avoid the negative pressure, the minimum density of particle is set to ρ min = ρ 0 . Numerical experiments show that, if not using this cut-off value, some long-term simulations will be stopped. On the contrary, the simulation will be finished successfully.

Boundary conditions and time-stepping
The solid boundary conditions are crucial for predicting the impact pressure and force in the simulations of fluid-structure interactions. In the present work, the dummy particle method [25] is applied. In this method, four layers of particles are placed to model solid boundary, and the pressure of the boundary particles can be interpolated from the neighbouring fluid particles, as follows where the subscripts s and f denote the solid boundary and fluid particles, respectively. a s is the acceleration of the solid boundary particle. To implement the free-slip condition, the viscous force between the solid particle and fluid particle is set to zero, which can be performed by setting φ = 0 in Eq. (4). For the moving rigid body, the motion can be obtained by where M, a, I, , x c represent the mass, the acceleration, the moment of interia, the angular acceleration, the position of the centroid, respectively. More details can be found in [26]. For the time integration scheme, the predictor-corrector scheme [27] is used here, and the time step satisfies the following CFL condition, where c max denotes the maximum of the speed of sound.

Water entry of a half-buoyant cylinder
In this section, a cylinder dropped into calm water is simulated and the results are compared with the experiment by Greenhow and Lin [28]. As shown in Fig. 1, the diameter of the cylinder is D = 0.11 m and the cylinder drops from a height of 0.445 m away from the water surface resulting in the water entry velocity v 0 = 2.95 m/s. The density of water is ρ w = 1000kg/m 3 and the density of air is ρ g = 1kg/m 3 . According to the experiment in [28], when the cylinder is completely submerged, its weight is half of the buoyancy. Therefore, this cylinder is a half-buoyant structure and the mass is 4.75 kg in 2D numerical model. The sketch of this problem is shown in Fig. 1. The width of the water tank is W = 1.0 m. The depth of water is H w = 0.65 m, and the height of air is H g = 0.65 m. In the present work, the speed of sound in water is c w = 113.3 m/s, and the particle spacing x is 1.25 mm.  Figure 2 presents snapshots of the water entry obtained by the SPH method at four time instants. To better distinguish water and air, the pressure field for water and the velocity field for air are displayed in Fig. 2. At the time t g/D =0.22, the cylinder begins to fall freely, forming a region of velocity disturbance around the structure and it can be observed that the pressure field of water is fairly smooth. As time approaching t g/D = 0.87, the velocity field of air extends further and the pressure field of water remains unchanged because the structure is far away from the water surface. As the cylinder drops, at t g/D =2.18, we can find that the velocity field is elongated, which indicates that the air is dragged down by this cylinder. Since the disturbance has not yet spread to the water surface, the pressure field still keeps the initial state at this time. At t g/D =2.83, this cylinder approaches the water surface and it starts to affect the pressure field of water. The following processes of this cylinder entering into the water are shown in Fig. 3.
In Fig. 3, a set of snapshots of the water entry obtained by the SPH method together with the experimental results [28] and the BEM results [29] are provided. The BEM results [29] are shown by solid lines. At the time t g/D =2.88, the cylinder touches the water surface producing a high pressure region near the water surface. As time goes on t g/D =3.02, 3.12, the water surface is expanded outward and the air escapes from the water at a high velocity. At t g/D =3.64, a V-shape open cavity near the water surface is formed. It can be found that the free surface profiles obtained by the SPH method show good agreement with the experimental results [28] and the BEM results [29].  Figure 4 gives time histories of the vertical displacement of this cylinder obtained by the present multiphase SPH model, where y 0 and y denote the initial coordinate and the current coordinate of the centroid of the cylinder, respectively. Note that the time t = 0 in Fig. 3 is the moment when the cylinder touches the water surface. The SPH results are compared against the experimental results [28] and the BEM results [29]. From the results, it can be observed that SPH results agree well with the experimental results [28] and the BEM results [29], which verifies the accuracy of the present SPH method.

Water entry of a fully buoyant cylinder
To further investigate the water entry, another cylinder water entry problem is conducted here. Compared with the previous subsection, when the cylinder is completely submerged, its weight is equal to the buoyancy. Therefore, the mass of this cylinder is 9.5 kg. To save the computational cost, this cylinder begins to fall near the water surface, keeping the initial velocity v 0 = 2.95 m/s. This treatment has no effect on the following water entry process. The sketch of the computational domain is same as that of the 2D half-buoyant cylinder water entry in the last subsection.  3 Numerical results for the 2D half-buoyant cylinder water entry problem. Left: experimental results [28] and BEM results (green line) [29]. Right: SPH results. The typical time instants from the top to the bottom are t g/D =2.88, 3.02, 3.12, 3.64, respectively Figure 5 displays snapshots of the water entry obtained by the SPH method in comparison with experimental data [28] and BEM results [29]. At the initial stage of the water entry t g/D =0.132, 0.84, 1.028, 1.878, similar water surface evolutions are obtained compared with those in the last subsection. During this period, under the impact of the structure, the water surface is torn into a large V-shape forming an open cavity, and these violent flow phenomena are reproduced well by SPH as compared to the BEM [29] results and the experimental data [28]. Figure 6 shows time histories of the vertical displacement of this cylinder obtained by SPH with different resolutions, where y 0 and y denote the initial coordinate and the current coordinate of the centroid of the cylinder, respectively. The SPH results are compared with the experimental results [28] and the BEM results [29]. One can observe that SPH Fig. 4 Time histories of the vertical displacement of the 2D half-buoyant cylinder obtained by the present multiphase SPH model, the BEM method [29] and the experiment [28] result obtained by the finest resolution is similar to that of BEM [29] but shows better agreement with the experimental data [28].

Water entry of a wedge with a deadrise angle of 25 •
In the field of ocean engineering, the process of a ship sailing on the sea is a complicated fluid-structure interaction problem. A sailing ship is subjected to the pressure from the water and the air. When a ship sails with a high speed, a considerable pressure may cause damage to the ship. Therefore, it is necessary for researchers to develop a method to predict these pressure loads. Because the ship bottom can be considered as a 2D wedge [14], in this section, the present multiphase SPH model is applied to predict the pressure loads in the process of a 2D wedge entering still water, which can provide some reference for predicting slamming loads of ship in the future.
The sketch of the wedge and the water tank is depicted in Fig. 7. The cross-sectional length of this wedge is L = 1.2 m and the deadrise angle is θ = 25 • . The mass of this wedge is 94 kg, and the initial water entry velocity is v 0 = 5.05 m/s. For this problem, the density of water is ρ w = 1000kg/m 3 and the density of air is ρ g = 1kg/m 3 . The width of the water tank is W = 1.5 m. The depth of water is H w = 1.1 m and the height of air is H g = 1.0 m. Four pressure measuring points, namely P1, P2, P3, P4 are arranged at the bottom of the wedge. The distance D = 50 mm between these measuring points is the same, as shown in Fig. 7. In the present work, the speed of sound in water is c w = 113.3 m/s, and the particle spacing x is 2.5 mm. As shown in Fig. 8, snapshots of the water entry obtained by the SPH method are provided. The pressure field for water and the velocity field for air are displayed in Fig. 8. At the time t g/L =0.014, the wedge comes into contact with the water, producing a Numerical results for the 2D neutrally buoyant cylinder water entry problem. Left: experimental results [28] and BEM results (green line) [29]. Right: SPH results. The typical time instants from the top to the bottom are t g/D =0.132, 0.84, 1.028, 1.878, respectively high pressure region near the interface and resulting in two air vortices up the vertex of wedge. As time advances to t g/L =0.034, the wedge further penetrates the water, and the pressure field in the water keeps expanding. Meanwhile, the air vortices gradually separate from the wedge. At this time, two jets adjacent to the interface between the wedge and the water are formed, and they make neighboring air escape at a large velocity. At t g/L =0.114, high pressure still exists between the wedge bottom and the water surface, and it should be concerned to ship designers.
The predicted pressure loads of four measuring points by SPH compared with the experimental data [30] are given in Fig. 9. The results obtained by single-phase SPH model in which the air is not considered are also presented. As stated in subsection 2.1, most of SPH practitioners adopted the constraint Eq. (9) to simulate water-air multiphase flows. This water entry was also simulated using this constraint. In this case, the speed  [29] and the experiment [28] of sound in air is 714 m/s higher than the real one. It can be observed that the variation tendency of the pressure obtained by the SPH method agrees well with the experimental data [30]. Specifically, the pressure peaks obtained by the multiphase SPH method with the real speed of sound in air are in better agreement with the experimental data [30] than those obtained by the other SPH methods, especially for the measuring point P1 and P2. Considering air effect, global pressure values by the present multiphase SPH model accord well with the experimental data [30] except for the peak values of P3 and P4 which are underestimated. The undershoot of the peak values was also mentioned in Koukouvinis et al. [31] in which this phenomenon is attributed to the rapid dynamics of the water sheet formation. From these results, it can be observed that the proposed mul-

Water entry of a vessel section
In comparison with previous subsections, a more complex structure, similar to a vessel section, is considered here, which can provide some reference for complex applications in the future. The sketch of the vessel section and the water tank is depicted in Fig. 10. The cross-sectional length of this vessel is L = 0.56 m and the height of side is D = 0.08  [32] m. The deadrise angle is θ = 25 • and the mass of this object is 72.25 kg/m (the mass of unit thickness is 72.25 kg for 2D simulation). The object falls from a height of h = 0.21 m that is measured from the upper edge of the object to the initial water surface, as shown in Fig. 10. The initial water entry velocity is v 0 = 1.22 m/s, which is consistent with the experiment [32]. The density of water is ρ w = 1000kg/m 3 and the density of air is ρ g = 1kg/m 3 . The width of the water tank is W = 1.8 m. The depth of water is H w = 1.0 m and the height of air is H g = 0.8 m. In the present work, the speed of sound in water is c w = 113.3 m/s, and the particle spacing x is 2.5 mm. Snapshots of the water entry obtained by the present multiphase SPH method are provided in Fig. 11. The pressure fields for water and the velocity fields for air are displayed to make the flow field clear. At the time t g/L =0.029, the object starts to drop, producing a small velocity disturbance around itself and resulting in four small air vortices near the vertex of this structure. The object comes into contact with the water causing a small slamming load for the water at t g/L =0.29. Meanwhile, the air vortices evolve into a curved strip shape near the side of the structure. With the entry into the water at t g/L =0.38, two jets adjacent to the interface between the object and the water are formed, and they make neighboring air escape at a large velocity. The pressure field of water is fairly smooth. At t g/L =0.44, 0.49, 0.59, the object continues to enter the water resulting in a violent change for the air domain. The previous air vortices disappear but two air flows with a large velocity occur due to the impact of the water jets. However, in contrast to the rapid change of the air field, the pressure of the whole water domain changes slightly.
The vertical displacement and the vertical acceleration of the 2D vessel section obtained by the single-phase SPH model and the multiphase SPH model along with the experimental data [32] are plotted in Fig. 12. It can be found that the results of the present multiphase SPH model are similar to those of the single-phase SPH model. The variation tendency of these curves obtained by both SPH models accords well with the experimental data [32]. However, there are some differences in terms of the amplitudes of these variables. The main reason may be that this structure is complex and the influence of 3D effect may affect the evolution of this problem.

Conclusions
In this paper, a multiphase Riemann-SPH model is developed to study the water entry problems. Different from most existing multiphase models in SPH, this model takes the real speed of sound in air into account. To avoid the occurrence of the negative pressure which brings some numerical difficulties to multiphase flow simulations, a cut-off value for the particle density is set. Thanks to the good performance of this multiphase model, the real speed of sound in air is considered for these water entries.
Several examples of water entry are well reproduced by the present multiphase model. We first simulate two cylinder water entries where the evolution of the water surface and the motion of the structure are investigated. Moreover, the results obtained by the present multiphase model are compared against the relevant experimental data which demonstrates it has good accuracy. Afterwards, a wedge water entry is modelled, and the impact pressure loads during the process of this water entry are focused. In contrast to the singlephase SPH model, the present multiphase model provides a more satisfactory prediction for the slamming loads. On this basis, the water entry of a more complex structure, a vessel section, is considered. As compared to the experimental data, general agreements are Fig. 12 Time histories of the vertical displacement and the acceleration of the 2D vessel section obtained in terms of the vertical displacement and acceleration of this vessel. To sum up, the present multiphase model can be a reliable tool for simulating violent water entries and engineering applications.
In the future work, the three-dimensional simulation of the water entry will be conducted, which may provide more flow details and satisfactory results. To reduce the computational cost, the efficiency of the SPH model needs to be improved. Besides, the process of high speed water entry in which there are some complex flow phenomena will be investigated.