MODELLING WATER ENTRY OF A WEDGE BY MULTIPHASE SPH METHOD

The hydrodynamic problem of two-dimensional wedge entering water is studied based on SPH model. A nonreflection boundary treatment for SPH method is proposed to reduce the size of computational domain. The details of water entry and enclosing are simulated using two phase SPH model. Both the water flow around the cavity and wedge and air flow inside the cavity are predicted simultaneity. Good agreement is obtained comparing experimental data in terms of both the hydrodynamics force exerting on the wedge and geometry of the cavity and jet.


INTRODUCTIONS
Water entry is part of the general fluid-structure impact problem in the field of naval architecture.SPH method is attractive on simulating the violent wave impact problems, e.g.Colagrossi & Landrini (2003) and Oger et al (2006).Applying SPH method on water entry, dambreaking, sloshing and breaking of solitary waves for both 2D, axisymmetry 2D and 3D problems are carried out by Gong et al (2008).For water entry problems, the air cavity enclosed by the water may significantly affect the local free surface profile and flow field, and then the hydrodynamics loads.Correctly simulating the multiphase flows may improve the force prediction, especially after cavity enclosing.

NUMERICAL MODELING OF MULTIPHASE SPH
Assuming ideal fluid flow, Euler equations are used for multiphase flows.Using general SPH formulations, e.g. , Monaghan and Kocharyan(1995), the governing equations can be written as where  is the fluid particle density, t is time, v  is the particle velocity, p is the pressure and g  is the gravitational acceleration.
ab W is cubic spline kernel function.ab  is the artificial viscosity term for enhancing the stability properties of numerical algorithms and takes the same form as Colagrossi & Landrini (2003).Comparing different artificial terms, e.g.Monaghan & Gingold(1983), Balsara (1991) and Colagrossi & Landrini (2003), preliminary simulation results show that the method used by Colagrossi & Landrini (2003) is the best one on the problem of dam-breaking.As a weak compressible model, the pressure is calculated from the density using the state equation.To prevent tensile instability, a background pressure  is added for all the particles, i.e.
prevent mutual-penetration of particles and to regularize the weakly compressible treatment of fluids, displacements of particles are calculated using XSPH variant given by Monaghan(1989).The value of  takes 0.5 in all the simulations presented in this paper.Following Colagrossi & Landrini (2003), a first order accurate interpolation scheme to re-initialize the density field is used ensuring the mass conservation is strictly satisfied.This procedure is applied only for the particles of the same phase.

SPONGER LAYER FOR SOUND WAVE
For weak compressible model, the speed of sound wave limits not only the size of time step but also the size of computational domain.For time step, the Courant number, based on the character value of particle distance and summation of sound speed and flow velocity, should not exceed an order of 0.3 in general scheme.In order to simulating incompressible flows using artificial compressible method, the Mach number should be controlled less than 0.1 in computational domain, as discussed by many SPH research, e.g.Colagrossi & Landrini (2003).Thus the sound speed should be large enough for consideration of limitation of Mach number and time step should be small enough for the limitation of Courant number.This is the first challenge for the computational expanse.Moreover, the disturbance of sound wave will travels a long distance when a finite time span is needed for flow simulating.Hence the computational domain should be large enough to avoid the reflected sound wave reaching the flow region of interested.This is another challenge for the computational expanse.
To save the computational cost, a sponger layer along the boundary of computation region is designed using the damping function to reduce the reflection of sound wave, as illustrated in Fig. the width of the sponger layer is set to be length of d .The density of particles inside the sponger layers is calculated using  Flow fields with and without this kind relaxation technique are compared.Two snapshots of pressure distribution during water entry are illustrated in Fig 2 .Without the sponger layers, the front of sound will reach the bottom wall at t=0.028s, and reach the wedge at t=0.031s as shown in the left panel.
Once the pressure disturbance hit the wedge, the pressure will deviated from the real physics value.Therefore the numerical results after t=0.031s are not physically correct.The remedy is that we have to enlarge the computation domain to delay the reflection of sound wave reaching the moving body.With the sponger layer implemented, the sound wave could be removed from the computation domain.In the comparison, nearly no reflection sound wave is observed at t=0.0248.At time instant of t=0.031s, the reflected pressure disturbance in the right panel is much more weak than the left panel in Fig 2, indicating the effectiveness of the proposed sponger layers.Furthermore, the wave front of the pressure disturbance is lower in the right panel as well, providing the chance of advancing computations.This kind technique can be applied for other flow variables such as velocity components in specified directions in a straightforward ways.

Rising bubbles
To verify the two phase flow solvers using SPH method, simulating of rising air bubble (phase Y) in a water column (phase X) is carried out, which is similar as the test case conducted by Colagrossi & Lanrini (2003).A circular bubble of air is free to rise in quiescent water.The density ratio is 1:1000 in present work.The radius of the air bubble is R , rising form a distant of 2R away from the bottom.
The state parameters are / 114  x From the comparison with Sussman et al (1994), it is seen that good agreements between two phase SPH model and data are obtained, especially for non-dimensional time less than 4 / t R g  .After that moment, the main part of the air bubble region match the data still rather well, as shown in Fig. 4.Only the detached small bubbles are smaller than the data.With SPH method, the size of the detached small bubbles is smaller than that predicted by level set method.Similar results were shown in Colagrossi & Lanrini (2003).At time instant of 6 / t R g  , small pressure fluctuations are observed at two particular spots near the small detached bubbles in Fig. 5, which is thought to be caused by the resolution is not enough for such small flow structures.Therefore, the numerical simulation is terminated at

Initial stage of water entry
To demonstrate the SPH code is correctly programmed, verification is conducted to verify the weak compressible SPH code and various boundary implementations.A symmetric wedge vertical impacting the free surface is simulated.To save the computational cost, a non-reflection boundary is designed using the damping function to reduce the reflection of sound wave, as discussed in previous section.The boundary pressure is obtained using an improved coupling boundary treatment approach.The initial particle distance and time step are . Total particles number reduces from 4.5million to 2.0 million through decreasing the computational domain size when sponger layer for sound wave is employed.Computational results were compared with the experimental and analytical results given by Zhao & Faltinsen (1993), see the bottom panel of Fig. 6.The pressure distribution along the wedge edge fits well the analytical results of boundary integration solution of potential theory.It should be pointed out that both these two kinds numerical solutions are deviate experimental data, which has been analyzed and explained as three dimensional effects in Zhao & Faltinsen (1993).By integrating the pressure along the wet surface of the wedge, the hydrodynamics force in vertical direction can be obtained, as shown in Fig. 7.As a consequence of pressure under prediction shown in Fig. 6, the drag force is under-predicted using SPH method in this study.

SIMULATION OF CAVITY ENCLOSING
In some circumstances, the single phase model could not predict the physical phenomenon after enclosing, wherein the pressure is zero even in the enclosed cavity.However, with two phase SPH model, the flows of entrapped air could be well simulated after enclosing of the cavity.In a water tank of width 0.9 meter, a wedge of 10cm is released at latitude of 0.26m (measured form the top of the wedge).The water entry process is recorded by a high speed video with 600 fps.The experiments are terminated after the cavity is closed.In numerical simulation, the particle distance is 1.25mm and time step is 1  s.There are 345,600 water particles and 214,250 air particles in the flow domain.Total flow time is 0.3s.Typical simulation of 300,000 time steps takes around 230 hours in a DELL T5400 workstation using 4 cores in parallel computation.To save computational efforts, absorbing boundary is implemented to remove the sound disturbance from the computational domain.By this approach, the computational time can be extended without the limitation of sound wave's reflection from solid boundaries.The position of the wedge can be fitted from the experimental data using the formulae, which is used as input in the numerical simulations.
The interface between air and water is digitized from the recording video and shown as hollow squares in the right panel of Fig. 8.With the resolution of 1.25mm, even the concave cavity shape around the corners of the wedge can be captured.Comparison between numerical results and physical experiments show a good agreement, as the particles distribution shown in Fig. 8.Only water particles are shown in the figure.
The flow field and pressure is shown in Fig. 9.The pressure of vented cavity equals the atmosphere pressure, shown in relative value.After enclosing, the pressure in the sealed cavity increases rapidly and equals to the surrounding water's pressure.Note the re-entrant jets are formed upward and downward after the deep enclosure of the cavity.In Fig. 9, the pressure inside the cavity varies a little bit after enclosing.The density varies less than O(1%) for the air in the cavity.By summation over all the particles inside the cavity, the total volume varies less than O(2%).With weak compressible model, the error magnitude of a few percent is acceptable.

CONCLUSIONS
To simulate the enclosing process after water entry, a two phase SPH method is applied.With fine particle distribution, details of water entry, including surface profile, pressure distribution and total force etc, could be well predicted.The enclosing of water entry is successfully simulated with proposed SPH model, providing a powerful Lagrangian approach for violent free surface flow calculation.
a and b take the value of 9 and 50 in present work.the ratio of the distance of the i-th particle and the neighbor solid boundary and width of sponger layer.

Figure 1
Figure 1 Setup of the sponger layers along the computational domain.
90,339 water particles and 5,022 air particles in the computational domain and 3864 boundary particles.Smoothing length number is 1time steps rising process takes around 36 hours in a DELL T5400 workstation using 4 cores in parallel computation.

Figure 2 Figure 3 Figure 4
Figure 2 Pressure distribution in the whole flow region without (left) and with (right) sponger layers.Top panel, t=0.0248s, bottom panel, t=0.031s.

Figure 5
Figure 5 Snapshots of the pressure and velocity field during the rising process.

Figure 6 Figure 7
Figure 6 Pressure field in the computational domain and comparison of the pressure distribution on the wedge.

Figure 8 Figure 9
Figure 8 Comparisons between experimental results and computed results of water phase at t=0.235s and t=0.267s.