Abdelraheem M. Aly , Mitsuteru Asai

1. Department of Mathematics, Faculty of Science, King Khalid University, Abha, Saudi Arabia

2. Department of Mathematics, Faculty of Science, South Valley University, Qena, Egypt

3. Department of Civil Engineering, Kyushu University, Fukuoka, Japan

Abstract: In this paper, we simulated the vertical impact of spheres on a water surface using three-dimensional incompressible smoothed particle hydrodynamics (3-D ISPH) method. The sphere motion is taken to be a rigid body motion and it is modeled by ISPH method. The governing equations are discretized and solved numerically using ISPH method. A stabilized incompressible SPH method by relaxing the density invariance condition is adopted. Here, we computed the motions of a rigid body by direct integration of the fluid pressure at the position of each particle on the body surface. The equations of translational and rotational motion were integrated in time domain to update the position of the rigid body at each time step. In this study, we improved the boundary treatment between fluid and fixed solid boundary by using virtual marker technique. In addition, an improved algorithm based on the virtual marker technique for the boundary particles is proposed to treat the moving boundary of the rigid body motion. The force exerted on the moving rigid boundary particles by the surrounding particles, is calculated by the SPH approximation at the virtual marker points. The applicability and efficiency of the current ISPH method are tested by comparison with reference experimental results.

Key words: Incompressible smoothed particle hydrodynamics (ISPH), free surface flow, sphere, rigid body, water entry

Introduction

The water entry problem has many applications such as naval hydrodynamics, ship slamming[1-2],stone skipping[3]and the locomotion of water walking creatures[4]. Numerous experimental, theoretical and numerical studies have been performed to study the water entry problems. Greenhow and Lin[5]conducted a series of experiments to show the considerable differences in the free surface deformation for the entry and exit of a circular cylinder. Zhao et al.[6]used both experiment and potential flow theory to investigate the water entry of a falling wedge.Kleefsman et al.[7], Panahi et al.[8]computed the water entry of a cylinder by solving the Navier-Stokes equation with a volume-of-fluid surface tracking using a finite volume formulation. Lin[9]used the concept of a locally relative stationary in his Reynolds-averaged Navier-Stokes (RANS) modeling to study the water entry of a circular cylinder with prescribed falling velocity. For a review of the water-entry literature, see Seddon and Moatamedi[10], Aristoff and Bush[11].Smoothed particle hydrodynamics (SPH) is a meshfree Lagrangian computational method that has been used for simulating fluid flows Gingold and Monoghan[12]. In this approach, the fluid is discretized into particles, properties of the particle are defined over a spatial distance, and the interaction of the particles is defined using equations of state. The particle-based nature of the definition is advantageous for capturing large deformations as it avoids problems such as mesh distortion associated with Lagrangian mesh-based methods. It is also advantageous compared with Eulerian fixed-mesh methods, as only the material domain is required to be meshed[13]. A comprehensive review of SPH is presented by Liu and Liu[14], which includes detailed descriptions , comparison with other fluid modelling approaches and almost 400 references.

Oger et al.[15]employed the 2-D SPH model with a fluid–solid coupling technique to study the water entry of a wedge with different degrees of freedom.The numerical model used a highly robust spatially varying particle resolution to improve the computational accuracy and efficiency. Liu et al.[16]implemented the two phase SPH model to simulate water entry of a wedge. A 2-D SPH model is implemented to study the water entry problem of a wedge entering the free surface as Gong et al.[17]. Aly et al.[18]adapted stabilized incompressible smoothed particle hydrodynamics (ISPH) method to simulate free falling rigid body into water domain. In the recent work of Koh et al.[19], the consistent particle method was proposed to eliminate pressure fluctuation in solving large-amplitude, free-surface motion. In this method, accompanied by an alternating of the kernel function by the Taylor series expansion-based partial differential operators, a zero-density-variation condition and a velocity-divergence-free condition is also combined with a source term of PPE to enforce fluid incompressibility. Aly and Asai[20]simulated fluid-structure interaction (FSI) on free surface flows using ISPH method. In their study, the rigid body is modeled using ISPH method by two different techniques. In the first technique, the solid particles are treated initially as fluid particles and after corrector step in projection method, the solid constraint is applied to get the rigid body motion. In the second technique, they computed the motions of a rigid body by direct integration of fluid pressure at the position of each particle on the body surface. They reported that, the second technique is more straightforward in a particle approach, but it still requires more improvement in terms of calculating the pressure exerted on the body surface and boundary treatment between moving rigid body and fluid.

The objective of the current study is to perform numerical simulations of water entry at different densities impacting bodies using stabilized ISPH method. We computed the motions of a rigid body by direct integration of fluid pressure at the position of each particle on the body surface. The equations of translational and rotational motions were integrated with respect to time to update the position of the rigid body at each time step. In this study, we improved the boundary treatment between both of fluid and fixed solid boundary and also between the fluid and moving rigid body by using virtual marker technique. The force exerted on the moving rigid boundary particles by the surrounding particles, is calculated by SPH approximation at its mapping virtual marker points.The applicability and efficiency of the current ISPH method are tested by comparison with reference experimental results.

1. SPH

The SPH approach is based on smoothing the hydrodynamic properties of a fluid through a smoothing function/kernel function. The fluid in the solution domain is represented by moving particles,which carry all relevant properties.

A spatial discretization using scattered particles,which is based on the SPH, is summarized. First, a physical scalar functionat a sampling pointican be represented by the following integral form

where W is a weight function called by smoothing kernel function in the SPH literature. In the smoothing kernel function,and h are the distance between neighboring particles and the smoothing length respectively. For the SPH numerical analysis, the integral Eq. (1) is approximated by a summation of contributions from neighboring particles in the support domain.

where the subscripts i and j indicate positions of the labeled particle andmeans representative mass related to particle j. Note that the triangle bracketmeans SPH approximation of a function Φ. The gradient of the scalar function can be assumed by using the above defined SPH approximation as follows

Similarity, the divergence of a vector functioncan be computed by

In this paper, the cubic spline function is utilized as a kernel function.

where η is a parameter to avoid a zero denominator,and its value is usually given by η2=0.0001h2.

2. Governing equations

Here, we consider the fluid is Newtonian, viscous and incompressible. The continuity and momentum equations for the fluid in the 3-D can be written in dimensional form as follows:

where u is the fluid velocity vector, ρ is the fluid density, p is the fluid pressure,is the dynamic viscosity of fluid, g is the gravity acceleration, and t indicates time. The turbulence stress τ is necessary to represent the effects of turbulence with coarse spatial grids, and its application in the particle simulation has been initially developed by Gotoh et al.[22].

3. Projection method

In the projection method[23], the velocity-pressure coupled problem has been solved separately for the velocity and pressure. Here, all the state variables may update from a previous time step to current time step.In the below, superscripts ()n and (+1)n indicate previous and current time step respectively. In the first predictor step, the intermediate velocity field can be evaluated by solving the following equation:

And the corrector step will introduce an effect of the pressure gradient term as follows:

In this study, according to Asai et al.[24], the pressure is obtained by solving the pressure Poisson equation (PPE) as follows

where (01)α≤≤ is the relaxation coefficient,is the temporal velocity. Hence, the corrector step can be implemented by substituting the pressure gradient with the solution of PPE.

4. ISPH formulations

The ISPH algorithm is implemented in a semiimplicit form in order to solve the incompressible viscous flow equations. The ISPH method is based on the calculation of an intermediate velocity from a momentum equation where the pressure gradients are omitted. Then, the pressure is evaluated through solving the PPE. The PPE after SPH interpolation is solved by a preconditioned diagonal scaling conjugate gradient PCG method[25]with a convergence tolerance(=1.0×10-9). Finally, the velocity is corrected using the evaluated pressure.

Here, the gradient of pressure and the divergence of the velocity are approximated as follows:

Laplacian operator for velocity and pressure can be approximated as follows:

5. Boundary treatments

The boundary condition of the rigid bodies has an important role to prevent the penetration and to reduce the error related to the truncation of the kernel function. In the literature, there mainly exists three solid boundary treatment methods, i.e., dynamic boundary[26-27], mirror boundary particles[15,28], and repulsive forces boundary[29-31]. Takeda et al.[32],Morris et al.[21]have introduced a special wall particles which can satisfy imposed boundary conditions.

In this study, a new boundary treatment using a virtual marker as Tanabe et al.[33]was proposed between the fluid and fixed boundary and as an extra studies between the fluid and rigid body. The concept of this treatment is to give a wall particle accurate physical properties, velocity and pressure. The procedure is summarized as follows:

Wall particle is placed on a grid-like structure with equally spaced in a solid boundary. Virtual marker is positioned in a symmetrical line to the wall particle across its solid boundary. Based on the concept of weighted average of neighboring particles,the velocity and pressure on the marker are interpolated from the fundamental equation of SPH method.The virtual marker was used as a computational point for giving the wall particle accurate physical properties and it is not directly related to the SPH approximation. Therefore, the density of the virtual marker does not effect on accuracies in SPH approximation and this will be helpful for robust the boundary condition. Moreover, the computational cost can be reduced compared with the ghost particle method because the virtual marker is created only once at the pre-process for the fixed wall particles case.

In order to satisfy the slip condition, the wall particle needs to be given the velocity, which is mirror-symmetric to the one on the virtual marker.This mirroring processingis given by the following equation.

where M is a second order tensor to implement the mirroring processing, and it is represented by the use of inward normal vector of the walland the kronecker delta δ as follows

In addition, in order to satisfy the non-slip condition, we will give the velocity for the wall particle from the point-symmetrical to the one on the virtual marker. Assuming that, R is a mirror-symmetric tensor and then the velocity is given as the Eq. (18).

Figure 1 presents the examples of velocity vectors for the wall particles to satisfy slip or no-slip conditions.

Fig. 1 (Color online) Virtual markers for the slip and no-slip boundary conditions

Fig. 2 (Color online) Configurations of the exerted pressure at the virtual marker on the body surface

Fig. 3 (Color online) Initial schematic diagram for water entry of sphere with diameter 0.0254 m

Table 1 Densities of the spheres used in this study. Each sphere has diameter 0.0254 m

Table 2 Analysis parameters for the water entry of sphere problem

In order to satisfy the Neumann pressure boundary conditions, we will give the wall boundary particles accurate pressure by mapping the virtual marker. Since the normal component of the velocity on the solid boundary is equal to zero, we should satisfy the following equation.

Fig. 4 (Color online) Time histories of sphere depth for different several spheres with densities 0.86, 1.14, 2.30 and 7.86, respectively

wherewsu is the velocity of the solid boundary. The next non-uniform pressure Neumann condition needs to be satisfied as follows

He toiled11 over a steep rocky shoulder of a hill, and there, just below him, was a stream dashing down a precipitous glen, and, almost beneath his feet, twinkling and flickering12 from the level of the torrent13, was a dim light as of a lamp

For satisfying the non-uniform pressure Neumann condition in SPH method, the pressure distributions at the wall boundary particles are evaluated by the following equation.

wherevp ,vf are the pressure and external force on the virtual markers evaluated by SPH approximation,wd is the distance from a solid boundary to the target wall particle andwρ is the density of the wall particles.

Fig. 6 (Color online) Snapshots of water entry for Nylon sphere with density ratio 1.14

6. Boundary treatment between fluid and moving rigid

Regarding the moving solid boundary of a falling body, difficulty is the accurate estimation of the forces exerted on the moving body by the surrounding water particles. Oger et al.[15], Liu et al.[34]employed the mirror particle method to enforce a free-slip boundary condition on the moving solid surface. The external fluid forces are evaluated by pressure integration of the fluid particles in the vicinity of the solid body boundary. Although this approach produces accurate results, it can be computationally expensive as reported by Shao[35]and may become unwieldy for corners or other geometrically complex solid surfaces[36]. Ren et al.[37]proposed dynamic boundary particles to treat the moving boundary of the floating body.

Fig. 7 (Color online) Snapshots of water entry for Teflon sphere with density ratio 2.30

In this study, the virtual marker is implemented for the moving boundary particles. The virtual marker is positioned in symmetrical to the moving boundary particles and this position is updated at each time step during the whole simulation. Pressure distributions on the marker are interpolated based on the concept of weighted average of neighboring particles, which is the fundamental equation of SPH method. Here, the normal component is calculated numerically from the color function and it updated at each time step depending on the position of moving boundary particles as follows:

Fig. 8 (Color online) Snapshots of water entry for steel sphere with density ratio 7.86

First, the color parameter is set for the surface boundary particles of the rigid body as:

The unit normal vector is computed as follows

The pressure in the moving boundary particles is given by

7. Treatment of moving rigid body

Fig. 9 (Color online) The horizontal velocity, vertical velocity and pressure distributions of water entry for steel sphere with density ratio 7.86

Table 3 Analysis parameters for free falling of several spheres on free surface cavity

Fig. 10 (Color online) Time histories for free falling of several spheres on the free surface cavity filled with water at density ratio 1.5

Here, we computed the motion of a rigid body by direct integration of the fluid pressure at the position of each particle on the body surface as hydrodynamic forces. The equations of translational and rotational motions were integrated in time to update the position of the rigid body surface at each time step. The equations of translation motions are described as follows:

where M is the mass of the body,is the hydrodynamic forces acting on the body surface andis the other external forces. The equation of rotational motions is described as

where θ˙ is the angular velocity,is the hydrodynamic moment andis the external moment.Here, the hydrodynamic forces and the hydrodynamic moment are calculated as:

where NS is the number of the body surface particles andis the area of the body surface at particle iis the position vector of the mass center. In addition, the computation of motion related to the hydrodynamics force acting on body surface has been shown in Fig. 2. In which, the stress vectoris calculated at each time step and it is given by the pressure at the virtual marker multiplied by the i-th component of the normal vector on the body surface.denotes the area considered as a discrete element of the body surface at particle i.

Fig. 11 (Color online) Horizontal velocity distributions for free falling of several spheres on the free surface cavity filled with water at density ratio 1.5

8. Results and discussion

Here, several different simulations of fluid-structure interactions have been introduced and discussed in details. The impact of sphere with different densities on the water tank has been simulated using ISPH method. And free falling of several spheres over free surface cavity has been simulated with two different density ratios between falling spheres and fluid.

Fig. 12 (Color online) Vertical velocity distributions for free falling of several spheres on the free surface cavity filled with water at density ratio 1.5

8.1 Water entry of decelrating spheres with different densities

Fig. 13 (Color online) Pressure distributions for free falling of several spheres on the free surface cavity filled with water at density ratio 1.5

In this section, the impacts of spheres with different densities on the water surface have been simulated using ISPH method. The sphere has diameter 0.0254 m and the initial setting for the current problem is taken as an experimental test for Aristoff et al.[38]. Initial schematic diagram for water entry of a sphere with diameter 0.0254 m has been shown in Fig.3. The water tank has dimensions (30×50×60)×10-6m3.The sphere is released from the rest and falls toward the water, reaching it with approximate speed≈The height of a sphere over water tank is 0.24025 m and then the impact speed is 2.17 m/s. Four one-inch diameter spheres, each made from a different material, were used in the present study. Their densities are reported in Table 1. The numerical model parameters such as initial particle distance and number of particles, etc. are reported in Table 2. The current numerical code was optimized and parallized to achieve simulation up to 107particles as introduced in the current model. Figure 4 shows the time histories of four impacting spheres that differ in density, but they have similar diameter and impact speed. In this work,we are interested only in measuring the depth of impacting spheres with different densities on the water tank. In Fig. 4, as the sphere density decreases the depth of the falling sphere inside water decreases.The current nu- merical results using improved boundary treatment between fluid and moving rigid body in ISPH method showed a good agreement in the cases of Polypropylene, Nylon and Teflon. In the case of a steel sphere with high density ratio 7.86, the current numerical results are still having a small gap compared to the experimental results, which can be improved in terms of boundary condition for the moving rigid body and high density ratio treatment as a future work.

The snapshots of the water entry for the four spheres, Polypropylene, Nylon, Teflon and steel have been shown in Figs. 5-8. Here, times from the sphere center to the free surface (=0)t are shown. Note that, the impacting spheres differ in density, but they have the same radius and impact speed. The evolution of the splash curtain is described by Aristoff and Bush[11]. In these figures, as the sphere descends into the calm water, the air cavity adjoins the sphere near its equator and its radial extents of the order of the sphere radius. Corresponding to the experimental results Aristoff et al.[38], the inertial expansion of the fluid is resisted by hydrostatic pressure, which eventually reverses the direction of the radial flow, thereby initiating cavity collapse. The collapse accelerates until the moment of pinch-off, at which the cavity is divided into two separate cavities. The upper cavity continues collapsing in such a way that a vigorous vertical jet is formed that may ascend well above the initial drop height of the sphere. The lower cavity remains attached to the sphere and may undergo oscillations. A relatively weak downward jet may also be observed to penetrate this lower cavity from above.The most obvious differences between the four impact sequences are the trajectories of the spheres,shown in Fig. 4, and the cavity shapes near pinch-off,as are highlighted in Figs. 5-8. As the sphere density decreases, several trends are readily apparent. First,the depth of pinch-off decreases. Second, the depth of the sphere at pinch-off decreases. Third, the pinch-off depth approaches the sphere depth at pinch-off.Finally, the pinch-off time decreases. The current ISPH simulations show a reasonable agreement with the experimental results as Aristoff et al.[38]. Figure 9 presents the horizontal velocity, vertical velocity and the pressure distributions of the water entry for the steel sphere with density ratio 7.86, respectively. It is found that, there are symmetry for the horizontal velocity of water particles around the impacted sphere.At impact region, the free surface is deformed and the particles over the impacted sphere raised up with high vertical velocity. The pressure distributions are increase at impacted region under the effects of impact sphere. One may notice that the pressure noise is redu-ced well.

Fig. 14 (Color online) Time histories for free falling of several spheres on free surface cavity filled with water with density ratio 0.5

8.2 Free falling of several spheres on free surfacecavity

In this section, we predicted numerically the free falling of several spheres on the free surface cavity filled with calm water with two different densities 1.5 and 0.5, respectively. The cubic cavity has dimensions 1 m and the particle distance is taken as 0.01 m. The analysis parameters for the current test have been shown in Table 3. The current setting of several spheres over free surface cavity has been listed as follows:

(1) Each sphere has diameter 0.1 m.

(2) The first sphere is putted at dimensions ( :X 0.1 m, :Y 1.1 m, :Z 0.15 m)

(3) Each sphere is separated from another by distance 0.1 m, then i.e. the second sphere is putted at dimensions ( :X 0.3 m, :Y 1.1 m, :Z 0.15 m) and the final sphere (number 18) is putted at dimensions( :X 0.7 m, :Y 1.3 m, :Z 0.75 m).

Fig. 15 (Color online) Horizontal velocity distributions for free falling of several spheres on the free surface cavity filled with water at density ratio 0.5

Figure 10 introduces the time histories for free falling of the several spheres on the free surface cavity filled with water at density ratio 1.5. It is clear that,the lower nine spheres impact the water surface generating hall inside calm water, which make an easy way for higher spheres to enter freely inside the water tank. Then, the higher spheres collide with the lower spheres inside water cavity around times 0.3-0.4 s.Here, since, the spheres densities have density ratios around 1.5 compared to the water density, then, the falling spheres are still going down until the bottom of the cavity. Due to their asymmetric mass distribution,the spheres start to rotate after released. Figures 11, 12 and 13 present the horizontal velocity, vertical velocity and the pressure distributions for free falling of several spheres on the free surface cavity filled with calm water at density ratio 1.5. The free surface above the falling spheres is nearly to hit each other and a cavity will form at the tail of the falling top spheres.Since here a single phase ISPH model is applied, the pressure at free surface and inside the cavity is always zero. As time going on, the falling spheres sink with almost unchanged rotation angles. In this study, the pressure distributions are smooth along the water tank and also around falling spheres at the hitting times until their reach to the bottom of cavity.

Fig. 16 (Color online) Vertical velocity distributions for free falling of several spheres on the free surface cavity filled with water at density ratio 0.5

In the second test, we simulated free falling of several spheres into the calm water with density ratio 0.5 as shown in Fig. 14. In this test, nine higher spheres collide with the nine lower spheres at times 0.25-0.50 s. After the collision, the spheres start to rise up to the water surface and as the time goes, the spheres continue with rising up until all the spheres rise up over the free surface. Then the spheres are still floating over the water tank. Figures 15, 16 and 17 present the horizontal velocity, vertical velocity and the pressure distributions for free falling of several spheres on the free surface cavity filled with water at density ratio 0.5. Here, the deformation of free surface with both of horizontal and vertical velocity distributions around the falling spheres have small values compared to the first numerical test (density ratio=1.5).The free surface deformation is almost negligible due to lower masses of the spheres. From Fig. 17, the contour of pressure around falling spheres at hitting times and sinking is acceptable which exhibiting a well sphere motion and floating.

Fig. 17 (Color online) Pressure distributions for free falling of several spheres on the free surface cavity filled with water at density ratio 0.5

Note, the virtual marker position is updated at each time step around the moving solid boundary and the evaluated pressures at virtual marker have been introduced to prevent penetration between the spheres during the collision. In this test, the evaluated pressure at the virtual marker is coming from surrounding fluid particles or from surrounding solid particles during the collision. Due to the limitations of experimental results for free falling of several spheres, we performed numerical prediction only without validations.

9. Conclusions

Improved incompressible smoothed particle hydrodynamics with a new treatment of boundary condition for moving solid body is proposed to simulate water entry of decelerating spheres. The pressure exerted on the moving boundary particles is calculated from the pressure of neighboring fluid particles around virtual marker points.

We computed the motions of a rigid body by direct integration of the fluid pressure at the position of each particle on the body surface. The equations of translational and rotational motions were integrated in time domain to update the position of the rigid body at each time step. The force exerted on the moving rigid boundary particles by the particles surrounding it is calculated by SPH approximation at its mapping virtual marker points. Impact of spheres with different densities over water tank is well simulated using the proposed model. The current numerical code was optimised and parallized to achieve simulation up to 107particles. The current numerical results using improved boundary treatment between fluid and moving rigid body have a good agreement with the experimental results. In addition, free falling of several spheres on the free surface cavity is also simulated using modified ISPH method. Here, the evaluated pressures at a virtual marker in solid have been introduced to prevent penetration between the spheres during their collision.

Acknowledgement

The first author (A. M. Aly) would like to express his gratitude to King Khalid University, Saudi Arabia for providing administrative and technical support.