Bow Flare Water Entry Impact Prediction and Simulation Based on Moving Particle Semi-Implicit Turbulence Method

Prandtl’s mixing lengthmethod and the k-epsilonmethod are introduced into theMoving Particle Semi-Implicit (MPS)method for the purpose of modeling turbulence effects associated with water entries of two-dimensional (2D) bow flare section.The presented numerical method is validated by comparing its numerical prediction with experimental data and other numerical results obtained from the Boundary Element Method (BEM). The time histories of the pressure and the vertical slamming force acting on the dropping ship section subjected to various conditions with different dropping velocity and inclined angles are analyzed.The results show that both the pressure and the vertical slamming force are in good agreement with the experimental data.


Introduction
For container ships the bow flare slamming can easily happen in rough seas; the slamming impact may lead to a series of instantaneous vibration and structure damage.
The slamming calculation theory is firstly proposed by Von Karman [1] and Wagner [2], which is mainly used to solve the 2D water entry problem.Wagner's theory is focused on flat plate water entry problem; based on the ideal incompressible fluid potential flow theory the linearized Bernoulli equation is used to calculate the pressure.Based on 2D vee-wedge theory, Bisplinghoff and Doherty [3] derived the unsteady state force using the flow pattern about an expanding prism.Zhao and Faltinsen [4] improved Wagner's theory to solve symmetric wedge water entry problem not considering the fluid separation; the body surface and nonlinear free surface conditions are simulated by numerical analysis; by this way the wedges of dead rise angle from 4 ∘ to 81 ∘ water entry problems are computed.In order to solve arbitrary shape water entry problem, Zhao et al. [5] proposed Boundary Element Method (BEM) to realize the practical numerical analysis; the slamming pressure is calculated by solving the nonlinear Bernoulli equation.This method can satisfy the moving boundary condition accurately when calculating the velocity potential.
When the dead rise angle is large, the original Wagner method overestimates the slamming pressure; in order to solve this problem, Korobkin [6] modified Logvinovich model (MLM), and the body boundary surface condition is solved by the Taylor expansion term of the velocity potential of the original Wagner model.Further, Korobkin [6] and Tassin et al. [7] expanded MLM to 3D calculation method to solve slamming problem.
On aspect of experiments, Marintek Sintef Group used a 30 ∘ dead rise angle wedge to accomplish the water entry impact experiment; the slamming impact is recorded in time domain; the results show that the pressure reached the maximum value before the flow separation, and after the flow separation the pressure rapidly attenuated at the separation point.Chuang [8] studied the air cushion effect from 1 ∘ to 15 ∘ wedge; when the maximum pressure appeared at the oblique elevation angle equal to 3 ∘ , the air between the wedge and the surface of the water more easily escapes than the smaller inclined angle.Aarsnes (1996) designed a bow flare section water entry experiment; four pressure sensors and two force sensors were installed on the model; the section dropped at different heights and inclined to simulate the real slamming phenomenon; the pressures and the vertical slamming forces were recorded in time domain.In this paper the numerical results will be compared with this experiment data.

Shock and Vibration
In this paper the main numerical method is Moving Particle Semi-Implicit method, which is a meshless method assuming the fluid as particles.Meshless method is a new numerical solution; the representative method is Smoothed Particle Hydrodynamic Method, which is firstly proposed by Lucy [9].Nowadays the typical meshless method includes many different realization solutions.Diffuse Element Method (DEM) is proposed by Nayroles et al. [10], which is a method for function approximation, by which the partial differential equations and fluid dynamic problems can be solved.Belytschko et al. (1994) used an element-free Galerkin method to solve elasticity and heat conduction problems, in which moving least-squares interpolants are used to construct the trial and test functions for the variation principle, and the moving least-squares interpolants and the choices of the weight functions are also discussed in detail.Batina [11] applied clouds of points in CFD algorithms, which are not required to be connected.This method is named Hpclouds gridless CFD method, which offers great potential for solving viscous flows about complex fluid calculations.Oñate et al. [12] proposed a procedure named "Finite Point Method" to solve convection-diffusion and fluid flow problem.Based on reproducing kernel and wavelet analysis a multiple scale method was developed by Liu and Chen [13], which can be used to simulate the wave numbers corresponding to spatial variables or the frequencies corresponding to the temporal variables.The Radial Basis Function (RBF) collocation method is a numerical approach for the solution of partial differential equations (PDE), in which no mesh is required.To transform the Finite Element Method to obtain a method which satisfies the meshless definitions, Idelsohn et al. [14] proposed Meshless Finite Element Method (MFEM).Hao et al. [15] introduced a shape function and derivatives to reproduce any order of monomial based on particle method, the approximation is constructed at individual discrete points where an approximation is desired, and this method is Moving Particle Finite Element Method (MPFEM).Zhu et al. [16] combined boundary integral equation method and moving least-squares approximation to form a direct meshless boundary integral equation method.
Moving Particle Semi-Implicit Method is firstly presented by Koshizuka [17] to calculate incompressible fluid flow problem, which is a meshless numerical method based on Lagrange particles.Further, Koshizuka and Oka [17,18] applied this method to simulate the dam breaking flow, and the results were compared with the experiments; in their research the kernel functions and the particle distance were discussed through comparison.Yoon [19] modified MPS method to combine the Lagrange particles with Euler particles; the typical advantage is that the arrangement is not necessary to be uniform; this is a meshless advection using flow-direction local-grid method (MAFL).Nowadays, MPS method is widely applied in different fields.Gotoh and Sakai [20] applied MPS method to simulate wave propagation problem; the breaking waves after flowing across the vertical wall were simulated.Chikazawa et al. [21] applied particle interaction models for differential operators to simulate the interaction between fluid and elastic plastic structures based on MPS method.Khayyer and Gotoh [22] proposed a corrected MPS method to track water surface in breaking waves; the formulations were revisited from the view of momentum conservation; modifications and corrections were made to ensure the momentum conservation in a particle-based calculation of viscous incompressible free surface flows.Zhang and Wan [23] investigated the roll motion of a two-dimensional floating body in low-amplitude regular waves through applying an improved MPS method; the efficiency and the accuracy in solving the interaction between the floating body and regular waves are improved.K. S. Kim and M.-H.Kim [24] extended MPS method to a multiphase system with multiple interfaces to solve the multifluid multi-interface problems; the results are compared with corresponding experiment results with three liquids.Xiang and Chen [25] proposed an interparticle force stabilization and consistency restoring MPS method to overcome the compressive instability that occurred under compressive stress states.A new and simple free surface detection criterion is proposed by Wang et al. [26] to enhance the free surface recognition in the MPS method.Sun et al. [27] used a weak coupling between MPS and BEM method to improve the computational efficiency for wave structure interaction simulation; better accuracy is obtained compared with published results.A multiphase MPS method is presented by Nabian and Farhadi [28] to develop a straightforward, robust, stable, and accurate meshfree numerical technique for modeling the dynamic behavior of free surface, incompressible, multiphase granular flows.Wada et al. [29] proposed a distributed point source method (DPSM) and the least square Moving Particle Semi-Implicit (LS-MPS) method to handle the phenomenon including fluid dynamics within the droplet, and the rotation of the droplet is successfully reproduced numerically and its acceleration is discussed and compared with those in the literature.Nath et al. [30] solved the time-dependent Navier-Stokes equations in stream function-vorticity form for lid-driven cavity flows using a meshless method based on fundamental and particular solution (MFS-MPS).Xu and Jin [31] in their study presented a local rheology model to calculate the effective viscosity and shear stress in granular flows in the numerical method; the final wave front in the simulations showed good agreement with the relationship obtained from experiments.
MPS method can be applied in water entry problems; Shibata et al. [32] calculated the acceleration of a free-fall lifeboat at water entry in three dimensions, in which the large deformation of the water surface, the splash of water at water entry, and the rigid-fluid interactions are predicted numerically.Alam et al. [33] used MPS method to solve the unsteady Navier-Stokes equation for incompressible fluid flows with and without the surface tension effect; equilateral prism-shaped object models were simulated to fall onto the free surface of the water; the hydrodynamic behaviors of water splash with and without the surface tension effect were presented to show the comparison and differences.Zhang et al. [34] presented an in-house solver based on MPS method to predict the process of free-falling wedge impacting on water; the numerical pressures, free surface elevations, and velocities of wedge show agreement with experimental data.
Turbulence models can be applied in meshless calculation method, based on eddy viscosity assumption and Langevin process.Violeau et al. [35] modified SPH method to prove the effective application in reproducing a large variety of flows.Shao and Gotoh [36] combined Large Eddy Simulation (LES) model with SPH method to simulate the coupled motion between progressive wave and floating curtain-wall type breakwater.Large Eddy Simulation (LES) filtering procedure was modified as a sort of LES Lagrangian filtering by Mascio et al. [37]; the closure formulas were derived for the additional terms to show the features to be used in SPH method.Dalrymple and Rogers [38] applied a two-dimensional LES model to wave propagation and interaction with coastal defense.Shao [39,40] applied two-equation - model to deal with turbulence and vortices during wave breaking and overtopping and it is coupled with the incompressible SPH numerical scheme.In order to reproduce cnoidal wave breaking on a slope under two different breaking conditions, spilling and plunging, Shao [39,40] coupled - model with the incompressible SPH method to realize the calculation; the results are in good agreement with the experiment data.Based on SPH method a subparticle-scale model is applied to account for the effect of turbulence, which is proposed by Kazemi et al. [41] to simulate depth-limited turbulent open channel flows over hydraulically rough beds.
The meshless method is now widely applied in computational fluid dynamics (CFD) to promise the complex free surface flow.However, the increasing flow complexity requires appropriate approaches for taking account of turbulent effects.The above recently developed turbulence models adapted to the SPH method are presented; each proposal is tested and validated on the basis of laboratory data.Theoretical or numerical solutions are available in the general field of turbulence free surface incompressible flows; they gave satisfactory results, but some progress should still be made in different methods in the future.
Compared with the other methods, in MPS method the fluid is separated into moving particles and diffusion terms, which are calculated during the motions.Therefore, the numerical diffusion problem in Eulerian methods does not take place.In this method, incompressibility is presented by keeping the density of particles constant during the computational time.The equations of continuity and momentum are converted into equations of interaction between particles in which all interactions are limited to a specific distance and weighing of interaction between two particles with distance of  is determined based on kernel function.The movement of each particle due to the interaction of neighboring particles is calculated using kernel weighing function approximation.Consequently, Laplacian gradient and divergence operators can be changed to consider the effect of moving particles.
In this study, MPS method is expanded to investigate the bow flare section of container ship water entry prediction.The effect of turbulence is considered using two turbulence models of Prandtl's mixing length method and - method.The modified method is based on the equations updated by Khayyer et al. [42] to provide more stability and accuracy; the results are compared with those obtained by BEM method and experiment data to verify the feasibility of this proposed method.

Realization of MPS Method
2.1.Governing Equations.As presented by Kenyon and Seeger [43], for a Newtonian incompressible viscous fluid, the continuity and momentum equations can be presented as follows: where u is velocity vector,  is time,  is fluid density,  is pressure, ] is kinematic viscosity, and F is body force vector.
For MPS method, as proposed by Koshizuka [17], the simulation process consists of two steps: the first step is explicit location modification of particles, and after the first step the number density is changed.The second step is implicit modification; the particle number density is modified to the initial value, so as to maintain the fluid incompressible.According to the derivative principle (2) can be written as where u +1 and u  are the velocity vectors of  + 1 and  steps after modification, u * is the intermediate velocity between two time steps, Δu * and Δu  are defined as After the velocity modification the particle number density is changed; in order to maintain the density invariable, (1) can be written as where  +1 and   are the particle number densities of  + 1 and  steps after modification,  * is the intermediate particle number density between two time steps, and Δ * and Δ  are defined as Shock and Vibration Substituting ( 5) and ( 3) into ( 1) and ( 2), in the first explicit modification, the fluid is considered as compressible.
Based on (7), the pressure can be obtained through the following equation: where  * is the particle number density after the implicit modification;  0 is the initial particle number density.The pressure obtained from ( 8) satisfies the incompressible condition, according to pressure value to modify the velocities and locations implicitly, the particle number density resumes the initial value, and the incompressible condition is satisfied.
According to velocity increment the location and velocity can be modified as where  is the location of particles, when the particle number density  * is obtained, Poisson's equation ( 8) can be solved, and the velocity is calculated implicitly; at last the velocity and location are modified again as follows: Through implicit pressure Poisson equation, the pressure for each particle can be presented as where  0 is the initial particle number density and  is the dimension number, and  + According to Prandtl's mixing length theory, the eddy viscosity is written as For MPS method the above equation can be presented as The - equation is deployed in the model, for each particle the kinetic energy can be defined as   , and the energy dissipation rate is defined as   ; the eddy viscosity is written as where   is empirical constant.
The turbulence energy can be expressed based on Lagrangian transport equation In scalar advection-diffusion equation ( 16),   is kinetic energy of particle , which is an energy source term, and   is energy dissipation of particle , which is a sink term.The Schmidt number   is equal to the ratio between the viscosity and energy diffusivity.
In continuous formulation the pressure and velocity can be written as Reynolds-averaged form; the Navier-Stokes equation is written as Based on Boussinesq assumption the Reynolds stress tensor R takes the following form: where   is eddy viscosity,  =       /2,   is the turbulent fluctuation velocity, and where letters  and  here denote spatial coordinates, and the production of kinetic energy  = −R : S, so  can be written as In this equation  is the scalar mean rate of strain, in water entry problem the jet flow is difficult to simulate, in this situation  is overestimated in case of large , the turbulence anisotropy is bounded by √  , and in MPS method for particle  the production of kinetic energy is written as In MPS method the rate of strain is presented as The energy dissipation rate is written as MPS form: The constant coefficients are valued as follows, given by Shakibaeinia and Jin [44]: After the implicit calculation velocities are obtained, the location of the particles should be modified as follows: The MPS realization considering the turbulence effect is presented in Figure 1.

Free Surface Condition.
The definition of the free surface particles is based on the particle number density normally; the free surface particles are defined as where  *  is the particle number density after the first velocity and location modification;  is determination parameter below 1.
In MPS method the boundary pressure particles are involved in the solution of Poisson equation, satisfying (25); in the simulation process these particles are defined as free surface particles incorrectly.So in the determination procedure a mixed method is applied; in the simulation the fluid is possible to be broken and the jet flow phenomenon happens; in these situations the particle number density is less than the others, so these particles can be defined as free surface particles directly; the determination equation can be expressed as The neighbor search method is proposed by Pan et al. [45], the initial searching is not necessary, the neighbor particles of the target particle between 1-2.1 particle distance should be searched, and this region is divided into eight parts, as shown in Figure 2.
When there is no particle in two adjacent parts of the divided region, the target particle is defined as the free surface particle.The particles can be predefined based on (26); if the particle number density is less than , this particle is considered to be located in the broken fluid or jet flow area.When the particle number density is greater than , the particle must be located in the internal fluid; here the coefficient  is valued between 0.97 and 0.99, which is proposed and testified by Koshizuka and Oka [17,18].

Fluid-Solid Boundary Condition.
In order to avoid the particle penetration and erroneous judgment of free surface particles, two kinds of particles are arranged on the fluid-solid boundary: the outside particles are wall particles which are not involved in the pressure Poisson calculation; the other kind of particles is wall pressure particles involved in the Poisson calculation.One layer of wall particles is arranged outside and one layer of wall pressure particles is arranged inside, as shown in Figure 3.

Kernel Function.
Kernel function is very important in the MPS method; the force between particles is defined based on kernel function, which is a weight function in number density model, Laplacian model, and gradient model.Different kernel functions have great influence on the MPS numerical simulation results; based on the conclusion obtained by Ataie-Ashtiani and Farhadi [46], (27) can give the best stability for MPS method, which is proposed by Shao and Lo [47].

Calculation Condition and Numerical Model
In order to compare the calculated results by MPS method with the experiment results, the models used in the calculation are the same as the models applied by other scholars.Aarsnes (1996) applied one section shape of container ship bow flare to complete a series of experiments; the section shapes are shown in Figure 4.The length of the section applied by Aarsnes is 1 m, the area measured in the test is 0.1 m at the middle of the section, the total weight of the section and equipment is 261 kg, and  1 - 4 are the pressure measured points.For the numerical calculation the cases are presented in Table 1.The size of the tank for numerical simulation is length 2.0 m × height 2.0 m; the water depth is 1 m; the initial particle distance is 0.01 m.

Numerical Results and Comparison
The initial particle arrangement and the boundary conditions are presented in Figure 5; the pressure is expressed in different color.The treatment of the moving solid and liquid is based on the research proposed by Koshizuka et al. [48].
There is no big difference between the simulation snapshots obtained by Prandtl mixing model and - model, so just the simulations obtained by - model of cases 1-4 are shown in Figures 6, 8, 10, and 12.It can be seen from Figure 6 that after the duration of 0.1 s the bottom of the section emerges in the water; the pressure on the bottom is  Here ℎ is the drop height,  0 is the initial water entry speeds given by Aarsnes, and the roll angle is the section entry angle to the vertical direction.
increasing sharply.After the duration of 0.13 s the flare part emerges in the water, the water resistance and the velocity of the particles are increasing, and the pressure continues to increase; the particles are moving along the tangential orientation.At the moment of the 0.16 s the jet flow is separating from the flare section; the particles are splashed outside.
In Figure 7 the BEM results are obtained by Sun [49]; the experiment data is given by Aarsnes (1996).As shown in Figure 7, compared with experiment the maximum values obtained by different methods appeared at different times; the peak value of the experiment may be prolonged because of the existence of the frictional resistance from the equipment. 1 is located at the lowest point of the bottom; as  1 touched the water the pressure increased to the peak value sharply, and then the pressure decreased with the slamming effect weakening.As seen in Figures 7(b) and 7(c) the pressures decreased when the flare section emerged in the water and became more stable, and after 0.14 s the section entry velocity reduced.Through the comparison it can be seen that the results obtained by MPS is close to the results calculated by    BEM but may get a better agreement with the experiment results.
According to the experiment presented by Aarsnes (1996), the section entry is with different angles and velocities, as shown in Figures 8 and 9, in which the entry angle is 9.8 ∘ and the velocity is 0.61 m/s and 2.43 m/s, the pressure reached the peak value when the jet flow touched flare part, on two sides of the section jet flow existed, and on the right side the jet flow is more serious.On the left side of the section the existence of air cavity leads to the second slamming on the section; although the pressure caused by air and water mixture cannot be predicted by this method, the main effect of the slamming is on the right side, which is calculated accurately.Compared with the results by BEM, MPS method  gives a better numerical prediction of the pressures, which is close to the experiment results.
As seen in Figures 10 and 12, the entry angle is larger, the right side experienced serious slamming; compared with case 2 and case 3 the pressure is larger.As seen in Figures 11 and  13, in the left side the pressure is lower because the dead rise angle for the right side is smaller, and for point  1 the pressure may be minus for a large entry angle, which is proved to be correct by the experiments.
The vertical slamming force is an overall effect; as shown in Figure 14, the slamming force is obtained by MPS methods which are close to the experiment data and BEM results.In case 2 the peak value of the slamming force happens when the jet flow contacts the flare part, and when the jet flow  at the left side arrives at the section surface the slamming force reaches another peak value; the two peak values are not equal which is caused by the asymmetry between the left and the right sides.As shown in case 2 and case 3, it is clear that the slamming force is more serious when the entry velocity is higher.After the bottom slamming the flow is separated to the right side and the left side; the second slamming leads to the two peak values on the vertical slamming force curve.Through comparison with the experiment and BEM method, the results calculated by MPS method considering the turbulence models are reasonable and more accurate.In order to analyze the influence of entry velocity on pressure, the water entry velocities 0.61 m/s, 0.9 m/s, 1.2 m/s, 1.8/s, 2.1 m/s, and 2.43 m/s are computed in vertical water entry and oblique water entry with an angle of 9.8 ∘ , and the results are shown in Figures 15 and 16.
In vertical water entry condition, as shown in Figure 15, when the velocity increases the maximum pressure increases; through comparison between the four pressure test points the pressure on  1 is obviously the largest, but when the velocity is increasing, the pressure on the flare is more close to the bottom pressure.In oblique water entry condition, when the velocity is small, the maximum slamming pressure is on the flare part on the entry side, and when the velocity increases the maximum pressure occurs through the flare part to the bottom.It is different from the vertical water entry; when the velocity is 0.61 m/s, the pressure on  1 is the smallest between the four points, but when the velocity is more than 2 m/s the pressure on  1 exceeds the pressures on the other points.When the velocity is equal to 1.2 m/s the pressure on  1 is smallest; the reason is that the asymmetry of the section shape leads to the variation of the dead rise angle, which is also testified by Chuang [8] in the experiment.
To compare with the results not considering turbulence effect, the above cases are calculated by original MPS method again.As shown in Figures 7, 9, 11, 13, and 14, the pressure values obtained by original MPS method are larger than the results considering the turbulence effect, and the peak values appeared a little earlier.The results also show that the pressure oscillation can be alleviated through coupling with turbulence models.In Figures 17 and 18, the vertical slamming forces obtained by MPS not considering turbulence are obviously larger than those obtained by MPS method coupled with Prandtl and - models.Compared with BEM method and experiment data, results calculated by the modified MPS method proposed in this paper are more close to the experiment results.

Conclusions
In this study, the turbulence models are considered in MPS method, which is applied in the bow flare water entry numerical simulation; the verification of the accuracy of this method is evaluated through comparison with the experiment data and the results obtained by BEM method.A typical bow flare Initial k and  Source term of k- equation Turbulence kinetic energy k Coefficient of energy dissipation  i

Figure 1 :
Figure 1: MPS realization process of MPS method considering turbulence effect.

Figure 4 :
Figure 4: Bow flare section and pressure sensor locations used in MPS simulation applied by Aarsnes (1996).

Figure 5 :
Figure 5: Initial particle arrangement and boundary condition.

Figure 6 :Figure 7 :
Figure 6: Particles colored by pressure after the impact with entry angle 0 with velocity 0.61 m/s.

Figure 8 :
Figure 8: Particles colored by pressure after the impact with entry angle 9.8 ∘ with velocity 0.61 m/s.

Figure 10 :
Figure 10: Particles colored by pressure after the impact with entry angle 9.8 ∘ with velocity 2.43 m/s.

Figure 12 :
Figure 12: Particles colored by pressure after the impact with entry angle 20.3 ∘ with velocity 0.75 m/s.
Figure 17: Vertical slamming force of vertical water entry.section is used in the calculation; the vertical and oblique water entry with different velocities are simulated by this method.Through computation of model in case of inclusion of turbulence, the comparison indicates that inclusion of turbulence in modeling free surfaces of water entry problem is feasible and the pressure numerical predictions tend to be more accurate.Through the comparison between the pressure results obtained by Prandtl mixing length model and - model, there is no obvious difference, but the computational cost of Prandtl model is lower than the - model.The accuracy of MPS method considered turbulence closures of Prandtl's mixing length theory and - model for calculating eddy viscosity that can be increased.Moreover, the changing regularity of the pressures on different points and the vertical slamming force in different water entry velocities are presented to analyze the effect of the velocity on slamming action; the results show that with the increasing of water entry velocity the slamming is more violent and because of the asymmetry of the section shape the maximum pressure position is different in oblique water entry condition.Consequently, the MPS method in case of turbulence closures Figure 18: Vertical slamming force of water entry with angle 9.8 ∘ . is feasible and effective in numerical simulation of water entry problem.: The intermediate particle number density between two time steps : The acceleration of gravity [9.81 m/s 2 ]   : Particle of the kinetic energy   : Energy dissipation rate   : Kinetic energy of particle    : The empirical constant *