Finding the best combination of numerical schemes for 2D SPH simulation of wedge water entry for a wide range of deadrise angles
 Author: Farsi Mohammad, Ghadimi Parviz
 Organization: Farsi Mohammad; Ghadimi Parviz
 Publish: International Journal of Naval Architecture and Ocean Engineering Volume 6, Issue3, p638~651, 30 Sep 2014

ABSTRACT
Main aim of this paper is to find the best combination of numerical schemes for 2D SPH simulation of wedge water entry. Diffusion term is considered as laminar, turbulent, and artificial viscosity. Density filter that seriously affects the pressure distribution is investigated by adopting no filter, first order filter, and second order filter. Validation of the results indicates that turbulent model and first order density filter can lead to more reasonable solutions. This simulation was then conducted for wedge water entry with wide range of deadrise angles including 10 degrees, 20 degrees, 30 degrees, 45 degrees, 60 degrees and 81 degrees, with extreme deadrise angles of 10 degrees, 60 degrees and 81 degrees being considered. Comparison of SPH results with BEM solutions has displayed favorable agreement. In two particular cases where experimental data are available, the SPH results are shown to be closer to the experiments than BEM solution. While, accuracy of the obtained results for moderate deadrise angles is desirable, numerical findings for very small or very large deadrise angles are also very reasonable .

KEYWORD
Water entry , SPH , Pressure distribution , Free surface.

INTRODUCTION
Water impact pressure arising during the impact and especially during the initial stages of water entry when the maximum pressure occurs, is extremely important in the design of marine structures. Abrate (2011) has presented an indepth review of the state of the art on the hull slamming or water entry problem. Based on this article and other recent papers, research conducted in the field of water entry is hereby surveyed.
Initial study of water entry problem for determination of water entry impact forces and pressures was done by Von Karman (1929). He used simple principles such as conservation of momentum and the concept of added mass. Wagner (1932) continued the theoretical solution of water impact problem. He analyzed the vertical water entry of a 2D wedge of small deadrise angle. Sedov (1934), on the other hand, extended Wagner’s work to study the impact of 2D wedge of larger deadrise angle.
In order to simulate the water entry problem, many numerical methods have been implemented by different researchers (Greenhow, 1988; Zhao et al., 1997; Yan and Ma, 2007; Zhang et al., 2010; Ghadimi et al., 2012; Ghadimi et al., 2013). Common methods of simulation have been gridbased. These schemes have encountered some computational difficulties when phenomena such as flow surface piercing, separation and large movement are involved in the modeling. As a result of this, it has become difficult to capture the movement of fluid and free surface flow due to a moving body. To overcome this difficulty, Lagrangian computational method such as Smoothed Particle Hydrodynamics (SPH) was adopted which is a mesh free method. SPH method was first developed for simulation of the particle motion in astrophysics by Gringold and Monaghan (1977) and Lucy (1977). Later, SPH was employed for simulating flow through porous media, large deformation of free surface flows, impact problem, Fluid Structure Interaction (FSI) and water wave generation (Shao, 2010; Dalrymple and Herault, 2009; Ferrari, 2010; GómezGesteira et al., 2005; GómezGesteira and Dalrymple, 2004; GómezGesteira et al., 2010; Monaghan, 1994; Rogers et al., 2008).
Application of the SPH method has also been extended to simulation of water impact problem. Oger et al. (2006) modeled symmetric wedge water impact of 30
degrees deadrise angle and asymmetric wedge water impact of 20degrees deadrise angle by SPH method and presented a new scheme to evaluate the pressure on solid boundaries. This work involved a spatially varying resolution based on variable smoothing length technique and a new formulation of the equations was proposed. Oger et al. (2008) also simulated the water impact of a sphere at high impact velocity. They used hundreds of millions of particles to model the sphere water impact. Kai et al. (2009) also studied 2D wedge water impact of 30degrees deadrise angle based on SPH method. They proposed a nonreflecting boundary condition to reduce the reflection of sound waves. A series of cases with different initial velocities was performed and maximum impact force for different conditions was obtained. Vepa et al. (2011) calculated pressure distribution of a 2D rigid cylinder during water entry by SPH method. They presented time history of pressure distribution and maximum pressure. They showed that by increasing the number of particles, the accuracy of solution increases. Veen and Gourlay (2012) developed a 2D SPH scheme and applied it to the ship slamming problem. Vandamme et al. (2011) modeled the entry and exit of a floating body using weakly compressible SPH method. They examined the water entry, fluid motions, hydrodynamic forces and movement of objects. They modeled 2D wedge of 30degrees and 45degrees deadrise angles drop tests and achieved results which showed good agreement with the experimental data.SPH method was also previously used by the current authors in order to simulate hydrodynamic phenomena and a particular WCSPH code was developed. This code has been considered and validated in two stages and for different applications: (1) the produced free surface elevation was validated by a dambreak problem (Ghadimi et al., 2012), and (2) the obtained pressure distribution was validated by a planing flat plate problem (Ghadimi et al., 2012). Based on these studies, it is clear that the WCSPH method is quite appropriate for accurately simulating the water entry problem. In the present study, the developed WCSPH code is used in order to investigate wedge water entry at a wide range of deadrise angles ranging from very small to moderate and to very large angles to find the best numerical schemes within SPH formulation. Accordingly, water impact of 2D wedge of 45
degrees of deadrise angle was first simulated under different numerical configurations. For this purpose, a 2DSPH model in conjunction with artificial viscosity, laminar viscosity and turbulent model were implemented. Subsequently, zero and first order density filter were utilized to remove pressure fluctuations. Furthermore, to obtain more accurate results, the number of particles used in the simulation, was increased. Pressure distribution and free surface elevation were presented. After finding the best numerical scheme, wedge water entry of 10degrees , 20degrees , 30degrees , 60degrees and 81degrees deadrise angles were simulated. To better present who has used SPH method and the extent of their work on wedge water entry, a summary of the SPH studies by different authors as well as the current study is presented in table 1.Based on the comparison presented in Table 1, one can conclude that the current SPH simulations of wedge water entry include a much wider range of deadrise angles (10
degrees , 20degrees , 30degrees , 45degrees , 60degrees and 81degrees , especially at the critical extreme angles of 10degrees , 60degrees and 80degrees at which the physical phenomenon gets very complicated) and a more extensive set of numerical configurations. Utilization of the mesh free weakly compressible SPH method as well as extensive physical and numerical features of the conducted studies can be considered as important contributions of the present paper.SPH FORMULATION
> Integral presentation of a function
SPH interpolation of a quantity
f(x) is based on the integral interpolantwhere the function
W is the kernel anddx is a differential volume element. The interpolant reproducesf exactly, if the kernel is a Dirac delta function, i.e.W =δ .In this study, cubic spline is used as a kernel function which is given by
where
α_{D} is in two dimensions> Particle approximation
To apply SPH interpolation for a fluid flow, computational domain is divided into a set of particles, as shown in Fig. 1. Every element has a mass
m_{j} , densityρ_{j} and positionx_{j} . It may be concluded that the functionf(x) can be written in the following form of discretized particle approximation:where
m_{j} andρ_{i} represent the mass and density of the jth particle, respectively.GOVERNING EQUATIONS
Fundamental physical laws of conservation are the basic governing equations of fluid dynamics which are conservation of mass, momentum and energy. For this problem, conservation of mass and momentum should be used.
In the SPH method, derivative of the density for particle
i must be determined based on the continuity equation (Monaghan, 1994) as inwhere the sum extends over all neighboring particles and
W is the smoothing kernel evaluated at the distance between particlesi andj .Also, velocity can be found from the momentum equation (Monaghan, 1994) as in
NUMERICAL FORMULATION OF SPH METHOD
> Viscosity treatment
To consider the diffusion term in the momentum equation, three different approaches including (1) artificial viscosity, (2) laminar viscosity and (3) laminar viscosity plus subparticle scale turbulence, are investigated.
Artificial viscosity
Monaghan (1992) introduced the artificial viscosity. Based on this concept, the momentum conservation equation can be written as
where is the gravitational acceleration.
The viscosity term, Π_{ij} , can also be given as
where 𝜖 ~ 0.01 is introduced to prevent a singularity when
r_{ij} =0 andv is defined bywhere and
h is the smoothing length.Laminar viscosity
Momentum equation in the form of laminar viscous stresses is given by
Laminar stress term is simplified by Lo and Shao (2002) to
where
v _{0} is the kinetic viscosity of laminar flow. As a result, Eq. (9) can be written asLaminar viscosity and subparticle scale (SPS) turbulence
Gotoh et al. (2004) introduced SubParticle Scale approach for modeling the turbulent effects. Momentum conservation equation can be presented as
where 𝜏 represents the subparticle scale stress tensor.
The eddy viscosity assumption is used to model the subparticle scale stress tensor
where
τ_{ij} is the subparticle stress tensor,v_{t} = [C_{s} Δl ]^{2} S  is the turbulence eddy viscosity,k is the subparticle scale turbulence kinetic energy,C_{I} = 0.0066 ,C_{s} is the Smagorinsky constant, Δl is the particleparticle spacing and , whereS_{ij} is the element of subparticle scale strain tensor (Issa, 2004; Pope, 2000). Dalrymple and Rogers (2006) presented Eq. (12) in SPH notation as> Time stepping scheme
SPH algorithm reduces the original partial differential equations to a set of ordinary differential equations. Then, any time stepping scheme for ordinary differential equations can be used.
In this study, PredictorCorrector algorithm is used. The momentum, density and position equations can be rewritten in the following form:
where
V_{i} represents the velocity contribution from particlei and from neighboring particles.Based on Eq. (14), PredictorCorrector scheme predicts the evolution in time, as in
Also can be calculated according equation of state by GomezGesteira et al. (2010). These values are then modified using forces at the half step, as in
At the end of any time step, the values are calculated as in
The pressure is calculated from density using
> Density filter
In the SPH simulation, pressure oscillations may be observed in the obtained results. To overcome this difficulty, Colagrossi and Landrini (2003) applied a filter over the density of the particles and reassigned a density to each particle. In the current paper, zero order filter and first order filter are considered as compiling options.
The following procedure is known as zero order filter and is applied every 100 time steps:
where the kernel has been corrected using a zerothorder correction:
Dilts (1999) developed the first order filter (moving least squares filter). This filter was applied successfully by Colagrossi and Landrini (2003) and Panizzo (2004). First order filter is a first order correction. Thus, the linear variation of the density field can be exactly reproduced as
The corrected kernel is also evaluated as follows:
> Pressure calculation
In the studies conducted by Monaghan (1994) and Batchelor (1974), the relationship between pressure and density is assumed to be
which is known as Tait's equation of state. Coefficient
B is related to the speed of sound. Using a value corresponding to the real value of the speed of sound in water, a very small time step must be chosen for the numerical modeling. Monaghan showed that the speed of sound could be artificially slowed significantly for fluids, without affecting the fluid motion. However, Monaghan (1994) suggested that the minimum sound speed should be approximately ten times greater than the maximum expected flow speeds. Therefore, the speed of sound is hereby set 15 times greater than the entrance velocity of the wedge for each studied case. Considering the fact that the entrance velocity in all cases is equal to 2m/s , the speed of sound is set to 30m/s for each case.RESULTS AND DISCUSSIONS
> Computational setup
The geometry of the considered problem is illustrated in Fig. 2. Width and height of the domain are equal to 1
m . Depth of water is considered to be 0.6m , while height and breadth of the wedge are set to be 0.3m and 0.15m , respectively. At the initial time, apex touches the undisturbed free surface of water.In the current study, the effects of density filter and viscosity treatment have been examined. Combinations of these parameters for each considered case is shown in Table 2. Focus of this study is on the accuracy of pressure distribution. Pressure distribution was also measured and compared with similarity solution of Zhao et al. (1997). Similarity solution is obtained through an analytical process and has been used by many researchers like Dobrovol’skaya (1969) and Zhao et al. (1997) to study water entry problem. Through this method, pressure distribution on the wedge can be calculated.
> Results
In the present work, pressure distribution of the wedge water entry is studied. Nine tests including entry of constant velocity of wedge of 45
degrees deadrise angle are investigated which are listed in Table 2. Velocity of the wedge was kept as 2m/s . In order to complete the simulation of water entry, both gravity and viscosity are included in these formulations. Simulation time was limited to 0.2sec . Free surface level at time0.05 sec . for each of the considered cases is presented in Fig. 3.It is quite apparent that the momentum treatment can affect the free surface level. Comparison of the obtained results of free surface against other numerical results (Muzaferija et al., 1999) shows that the obtained free surface level by laminar and turbulent viscosity is better than the artificial viscosity. The results in cases 4 and 5 also have good agreement with other results (Muzaferija et al., 1999). Pressure distribution in cases 4 and 5 are compared with similarity solution (Zhao et al., 1997) in Fig. 4. There is no critical difference between the obtained results of cases 4 and 5. Therefore, it is concluded that both laminar viscosity and turbulent viscosity can be used to simulate the water entry problem, because the duration of this phenomenon is very short.
In order to obtain more accurate pressure distribution, cases 4 and 5 with about 500,000 particles were repeated again. These results are presented in Fig. 4. On the horizontal axis, Z/vt is the nondimensional length in which v is the entrance velocity of the wedge and t is the time at which the maximum pressure occurs.
There seems to be an improvement in the results. As expected, by increasing the number of particles, the difference between the obtained results of laminar viscosity and turbulent viscosity became even more smaller.
It is also quite apparent that turbulent modeling and using first order density filter lead to best results. Therefore, simulation of water entry of wedges of 10
degrees , 20degrees , 30degrees , 45degrees , 60degrees and 81degrees deadrise angles along with the mentioned compiling options was conducted and pressure distribution was obtained.Water impact generates sound wave that reaches numerical domain boundaries. Reflection of this wave from boundaries goes back to the wedge and may affect pressure distribution. Thus, at this stage, tank size increases to 4
m . These simulations were performed by 2,000,000 particles that helped improving the results, substantially.Fig. 6 shows the effect of domain size on pressure distribution of the wedge of 45
degrees deadrise angle. It is clear that by increasing the tank size, effect of reflecting sound wave decreases and accuracy of the pressure distribution increases. This strategy has been pointed out by Oger et al. (2006).Simulation of wedge water entry of 10, 60 and 81 deadrise angles by SPH method had not previously been reported by other researchers and was conducted in the current study. Pressure distributions of these wedges are presented in Fig. 7. As evidenced in this figure, there is good agreement between the current SPH results and the boundary element solution (Zhao and Faltinsen, 1993). On the other hand, a comparison is presented between the current results and experimental results of Zhao et al. (1997) for the wedge deadrise angle of 30
degrees . Based on the obtained results, one can conclude that accuracy of the pressure distribution for wedges with very small and very large deadrise angles is less than the accuracy of pressure distribution for wedges with moderate deadrise angle. This may be attributed to more complex nature of this phenomenon at the extreme angles.Differences between SPH results and BEM solutions are shown in Table 3. It is observed that at extreme deadrise angles, i.e. 10
degrees , 60degrees and 81degrees , relative error between pressure peak of the current study and BEM result is large. This has been reported by some researchers to be due to high nonlinearity of this phenomenon at these extreme deadrise angles. It is also quite apparent that, as deadrise angle starts increasing from 10degrees to 81degrees , relative error between pressure peak of the current study and BEM results decreases in oscillating form. This trend is also repeated by the force coefficient curve representing the area under pressure distribution curve.Free surface elevations related to wedges of 10
degrees , 20degrees , 30degrees , 45degrees , 60degrees and 81degrees deadrise angles are presented in Fig. 8. Nondimensional free surface is obtained by capturing the free surface at the time of maximum pressure occurrence.As evidenced in Fig. 8, favorable agreement is observed between SPH results and BEM solutions. Comparison between the results of the current study against experimental results (Tveitnes et al., 2008) for wedge of 45
degrees deadrise angle is presented in Fig. 8. The SPH results are shown to be closer to the experimental data than BEM solution. In all the considered cases, predicted free surfaces by BEM method have a long thin jet flow at the top of the wedge as opposed to the SPH solutions. In fact, speed of the jet flow in BEM solutions is more than the SPH results. This may be attributed to viscosity and gravity. The combined effect of viscosity and gravity reduce the jet velocity and slow down the free surface from rising up. In BEM method, viscosity and gravity not considered. However, these parameters are taken into consideration in the SPH method. As a result of the influence of these parameters, a difference in response between the two methods is quite natural. On the other hand, as the deadrise angle increases, relative error between BEM solutions and SPH results decreases.CONCLUSIONS
In this paper, wedge water entry for a wide range of deadrise angles is simulated by SPH method. In order to find the best numerical configuration, nine different test cases are considered. These cases describe the water entry of wedge of 45
degrees deadrise angle at constant speed. Pressure distribution and free surface level of the considered cases are presented. It is found that first order density filter and turbulent modeling is the best configuration. Pressure distribution is compared against similarity solution and favorable agreement is displayed between them. Also, the effect of number of particles is studied. By increasing the number of particles, results substantially improve and difference between laminar viscosity result and the result of turbulent viscosity is reduced. Furthermore, effect of tank size is studied for wedge of 45degrees deadrise angle. It is found that larger tank with no reflected sound wave leads to better results.In order to assure that the obtained numerical configuration can also lead to good results for wide range of deadrise angles, wedge water entry of 10
degrees , 20degrees , 30degrees , 45degrees , 60degrees and 81degrees deadrise angles is simulated with 2,000,000 particles. Pressure distribution of these wedges displayed good agreement with BEM solutions.While the achieved numerical configuration within SPH formulation has lead to desirably accurate simulation of wedge water entry for a wide range of moderate deadrise angles, it is also quite apparent that numerical findings for particular cases of relatively small and relatively large deadrise angles are also very reasonable.

[Table 1] Summary of wedge water entry studies by SPH method.

[]

[]

[Fig. 1] Particle approximations using particles within the support domain of the smoothing function W for particle i.

[]

[]

[]

[]

[]

[]

[]

[]

[]

[]

[]

[]

[]

[]

[]

[]

[]

[]

[]

[]

[]

[Fig. 2] Geometry of simulation.

[Table 2] Compiling options.

[Fig. 3] Free surface level at 0.05 sec.

[Fig. 4] Pressure distribution of cases 4 and 5 (about 80,000 particles) compared with similarity solution (Zhao et al., 1997).

[Fig. 5] Pressure distribution of cases 4 and 5 (about 500,000 particles) compared with similarity solution.

[Fig. 6] Effect of tank size on pressure distribution for a wedge of 45 degrees deadrise angle.

[Fig. 7] Pressure distribution (about 2,000,000 particles) compared against boundary element solutions (Zhao and Faltinsen, 1993) for different deadrise angles (in the case of 30 degrees deadrise angle, experimental results (Zhao et al., 1997) are also included).

[Table 3] Comparison of SPH results and BEM solutions.

[Fig. 8] Free surface elevation for wedges of deadrise angles 10°, 20°, 30°, 45°, 60° and 81° using about 2,000,000 particles. Nondimensional free surfaces obtained by the current study and BEM method (Zhao et al., 1993) are plotted on the right side. Experimental result of free surface for 45 degrees wedge water entry is related to the work of (Tveitnes et al., 2008).