# Numerical Investigation of Water Entry Problem of Pounders with Different Geometric Shapes and Drop Heights for Dynamic Compaction of the Seabed

Authors: Mohammad Hossein Taghizadeh Valdi; Mohammad Reza Atrechian; Ata Jafary Shalkoohy; Elham Chavoshi
Journal: Geofluids
 (2018)

---

## Abstract

The water entry problem of three-dimensional pounders with different geometric shapes of cube, cylinder, sphere, pyramid, and cone was numerically simulated by the commercial software Abaqus, and the effects of pounder shape and drop height from the free surface of water on deepwater displacement and velocity as well as pinch-off time and depth were investigated. An explicit dynamic analysis method was employed to model fluid-structure interactions using a Coupled Eulerian-Lagrangian (CEL) formulation. The simulation results are verified by showing the computed shape of the air cavity, displacement of sphere, pinch-off time, and depth which all agreed with the experimental results. The results reveal that the drag force of water has the highest and lowest effect on cubical and conical pounders, respectively. Increasing the pounder drop height up to the critical height leads to increased pounder velocity while impacting the model bed and more than the critical drop height has a reverse effect on pounder impact velocity. Pinch-off time is a very weak function of pounder impact velocity; but pinch-off depth increases linearly with increased impact velocity.


---

## Body

## 1. Introduction

Dynamic compaction is regarded as an appropriate method for soil improvement after its introduction in 1970 by Louis Menard. This method is based on the drop of heavy pounders of 8 to 25 tons from a height of 10 to 20 m on loose soil surfaces, although lighter or heavier weights and other drop heights are occasionally used [1]. The pounders impact on the soil surface and the resulting vibration on the surface and underlying layers leads to the compaction of the soil grains. As a result, the relative density of soil and its stiffness and penetration resistance increase. This leads to increased soil shear strength and bearing capacity and decreased liquefaction. The pounders used in this method have different shapes and weights due to the operations implementation on the coast and in the sea. Onshore dynamic compaction operations mostly make use of cubical or cylindrical pounders; but in offshore dynamic compaction operations, a large amount of pounder velocity depreciates due to its impact on the water’s surface and deepwater movement, meaning that the pounder does not have the energy necessary for soil vibration and compaction of the surface and its underlying layers. Therefore, the geometric shape of the pounder must be selected so that the least velocity depreciation occurs during deepwater movement to improve a proper radius and depth of soil after seabed impact and soil layer vibration.


### 1.1. Experimental Studies

The water entry problem of solid projectiles has been an interesting topic in fluid mechanics for more than 80 years. The study of the underwater motion of projectiles has wide applications in military fields, especially for naval industries since World War II. Furthermore, accurate prediction of water impact forces has particular importance in the design of marine structures and projectiles employed in dynamic compaction of seabed. In the last three decades, a lot of researchers addressed the water entry problem using various methods. Watanabe [2, 3] investigated the impact of cones upon water. Cones of different weights were dropped into water from different heights, while being equipped with a piezoelectric sensor connected to a vibrometer, to register force changes resulting from water impact. Such experimental results were very helpful in approving analytical results and relations. May and Woodhull [4–6] published three important papers to find the drag coefficient and the scaling relationship in water entry of sphere projectiles by considering the Reynolds and Froude numbers. They used an experimental method to study the air cavity formed during projectile impact on water using high-speed photography. May’s work indicated that Froude scaling is a good first approximation to use when describing cavity behaviours such as deep closure. New et al. [7] investigated the impact and hydrodynamic characteristics of water entry associated with prismatic bodies with different fore section shapes and entry angles. Transient loading on the fore section of the prismatic body was measured with a triaxial accelerometer. Splashes, cavity formation, and closure procedures were observed with high-speed photography and video recordings. Anghileri and Spizzica [8] studied the normal impact of a rigid sphere with a water surface using finite element methods and certain experiments to gain variations in sphere acceleration during impact in order to verify the applied numerical method. Their results were in good agreement with the experimental data. Engle and Lewis [9] compared the results of the numerical method for predicting hydrodynamic impacts with experimental results of symmetrical wedges with different initial impact velocities. Their studies catered for the accuracy and validity of the various methods. Aristoff et al. [10] performed an experimental and theoretical study by impacting spheres with different densities on a water surface. A theoretical model was developed that yielded simple expressions for pinch-off time and depth, as well as the volume of air entrained by the sphere. Their study demonstrated that the dynamics of the air cavity formed in the wake of a sphere falling through the water surface can be altered considerably by the density of the sphere.


### 1.2. Analytical Solution

Von Karman [11] developed analytical solutions for the water entry problem of two-dimensional wedge bodies. He theoretically presented formulas using asymptotic theory. Wagner [12] presented a more precise analytical model which took the free surface elevation of the water into consideration. He developed an asymptotic solution for water entry by two-dimensional projectiles with small local deadrise angles. Vorus [13] proposed a novel analytical model for the water entry problem of bodies. In his model, nonlinearity in the boundary conditions was considered but the geometrical nonlinearity of the impact problem was neglected. Miloh [14] obtained analytical solutions based on several mathematical and physical assumptions for the small-time slamming coefficient and wetting factor of a rigid spherical shape in a vertical and oblique water entry. Aristoff and Bush [15] presented the results of a combined experimental and theoretical investigation of the normal impact of hydrophobic spheres on a water surface by determining the shape of the air cavity behind the sphere, pinch-off time, and depth. Ghadimi et al. [16] employed an analytical method and a numerical solution to study the water entry problem of a circular section. They demonstrated a procedure for the step by step derivation of the analytical formulas and presented a comparison between the linear and nonlinear solutions. Tassin et al. [17] investigated the water impact problem of bodies by CFD code OpenFOAM. They studied the water entry and exit of a two-dimensional body with a time-varying shape by using a modified Logvinovich model during the entry stage and the Von Karman model during the exit stage. Yao et al. [18] developed a theoretical model to describe the evolution of the cavity shape, including the time evolution of the cavity on fixed locations and its location evolution at fixed times.


### 1.3. Numerical Simulations

Park et al. [19] presented a numerical method for calculating the impact and buckling forces of bodies after a quick water entry. By neglecting fluid viscosity, potential flow was assumed to be nonviscous. This assumption greatly reduced computational time while solution accuracy was preserved by considering nonlinear free surface conditions. They solved this issue using a source panel method, in which an object is divided into panels or elements and each panel can be fully floating, partially floating, or fully out of the water. The results were compared with laboratory data and valuable results were obtained. Battistin and Iafrati [20] conducted a numerical study on normal entry into water of a two-dimensional symmetrical or asymmetrical object with a desirable form. This study considered potential flow of an incompressible fluid by ignoring the effects of gravity and adhesion. In this work, simulation was done using the boundary element method, and the boundary conditions of the free surface were fully nonlinear. The focus of this investigation was on evaluating pressure distribution and overall hydrodynamic loads applied to the colliding body. They verified the results for cylinder, sphere, and cone forms. Korobkin and Ohkusu [21] performed an interesting study on the hydroelastic coupling of finite element models using Wagner’s theory for water impacts. The aim of the study was to measure direct coupling possibilities using finite element methods for the structure and Wagner’s theory for hydrodynamic loads of impact of an elastic object with the water’s surface. Kleefsman et al. [22] numerically investigated the dambreak problem and water entry problems. The numerical method was based on Navier-Stokes equations that describe the flow of an incompressible viscous fluid. The equations were discretized onto a fixed Cartesian grid using the finite volume method. Kim et al. [23] analyzed the water entry problem of symmetrical bodies using the particle hydrodynamic method. Yang and Qiu [24] investigated the two-dimensional water entry problems of symmetric and asymmetric wedges with various deadrise angles using the Constrained Interpolation Profile (CIP) method. The effect of air compressibility for small deadrise angles was also discussed in their work. Fairlie-Clarke and Tveitnes [25] conducted a study on impact of wedge-shaped parts with a water surface at a constant velocity. In their study, Computational Fluid Dynamic (CFD) analysis was used to determine the momentum of added mass, momentum of flow, and gravity’s effect during the entry of wedge-shaped objects at a constant velocity and a dead angle of 5 to 45 degrees into water. A numerical solution based on finite volume method was implemented using Fluent software to solve equations of mass conservation and momentum to obtain flow fields. Hydrodynamic forces were not significantly changed by the effect of gravity on the flow fields. Yang and Qiu [26] computed slamming forces on two- and three-dimensional bodies entering calm water with a constant velocity. The nonlinear water entry problem was governed by Navier-Stokes equations and solved using a CIP-based finite difference method on a fixed Cartesian grid. Wu [27] developed a high order boundary element method for the complex velocity potential problem. He focused on application of the water entry problem for a wedge with varying speeds. His method ensured not only the continuity of the potential at the nodes of each element but also the velocity, allowing it to be applied to a variety of velocity potential problems. The continuity of the velocity achieved herein was particularly important for this kind of nonlinear free surface flow problem, because when the time stepping method is used, the free surface is updated using the velocity obtained at each node and the accuracy of the velocity is therefore crucial. Ahmadzadeh et al. [28] developed a three-dimensional numerical model to simulate a transient free-falling sphere impacting a free water surface using the Coupled Eulerian-Lagrangian formulations included in the commercial software Abaqus with consideration for important parameters such as viscosity, compressibility inertia, and variable downward velocity. Their study was successful in determining the shape of the air cavity, displacement of sphere versus time, pinch-off time, and depth compared to the experimental results found by Aristoff et al. [10]. Nguyen et al. [29] studied the water impact of various three-dimensional geometries, namely, a hemisphere, two cones, and a free-falling wedge, with an implicit algorithm based on a dual-time pseudo-compressibility method. Flow fields of incompressible viscous fluids were solved using unsteady RANS equations. A second-order VOF interface tracking algorithm was developed in a generalized curvilinear coordinate system to track the interface between the two phases in the computational domain. Sensitivity analysis with respect to the spatial grid resolution and the time step was performed to assess the accuracy of the results. Free surface deformation, pressure coefficients, impact velocities, and vertical acceleration during impact were compared with available experimental data and asymptotic theories, showing good agreement at a reasonable computational cost. S. Kim and N. Kim [30] performed integrated dynamic modeling of supercavitating vehicles and established 6-DOF equations by defining hydrodynamic forces and moments. These simulations also demonstrated the instability of supercavitating vehicles caused by the asymmetric immersion depth of coupled horizontal or vertical fins. Nguyen et al. [31] used a moving Chimera grid method to predict the real-time motion of water entry bodies by combining the 6-DOF rigid body motion model and the numerical calculation of the multiphase flow field, by which the coupled effect of supercavity and moving body was obtained. Mirzaei et al. [32] held the idea that existing planning force models weaken nonlinear interactions among the solid, liquid, and gaseous phases and are often too simple and therefore incorrect. They established the 6-DOF motion model for supercavitating vehicles by redefining the spatial planning force and developed a hybrid control scheme based on online planning force identification, which increased the stabilization of supercavitating vehicles.In this study, the water entry problem of three-dimensional pounders with different geometric shapes including cube, cylinder, sphere, pyramid, and cone was investigated numerically using commercial finite element code Abaqus 6.14-2 and an explicit dynamic analysis method was employed to model fluid-structure interactions using a Coupled Eulerian-Lagrangian (CEL) formulation [33]. This investigation focused on the problem of rigid pounders touching the free surface of a stationary viscous fluid at time t=0 with a variable velocity and finding for t>0, the shape of cavity, pinch-off time, and depth of the rigid pounders. It should be noted that the fluid was assumed to be viscous and compressible. Numerical simulations of solid-fluid interactions were performed using the Coupled Eulerian-Lagrangian (CEL) formulation and the remeshing techniques available in Abaqus 6.14-2 software [34]. The pounder was considered as a rigid solid and its mesh was created in Lagrangian form. The water was regarded as a compressible fluid and its mesh was created in Eulerian form. After validating the numerical results with recent experimental data, the effect of the pounder shape on its deepwater displacement and velocity was investigated. Then the effect of increasing the pounder drop height from the water’s surface on displacement, submersion velocity, pinch-off time, and depth was investigated. The critical drop height of the pounder, which is indicative of the maximum effect of the pounders drop height from the water’s surface on deepwater movement, was finally determined.


## 1.1. Experimental Studies

The water entry problem of solid projectiles has been an interesting topic in fluid mechanics for more than 80 years. The study of the underwater motion of projectiles has wide applications in military fields, especially for naval industries since World War II. Furthermore, accurate prediction of water impact forces has particular importance in the design of marine structures and projectiles employed in dynamic compaction of seabed. In the last three decades, a lot of researchers addressed the water entry problem using various methods. Watanabe [2, 3] investigated the impact of cones upon water. Cones of different weights were dropped into water from different heights, while being equipped with a piezoelectric sensor connected to a vibrometer, to register force changes resulting from water impact. Such experimental results were very helpful in approving analytical results and relations. May and Woodhull [4–6] published three important papers to find the drag coefficient and the scaling relationship in water entry of sphere projectiles by considering the Reynolds and Froude numbers. They used an experimental method to study the air cavity formed during projectile impact on water using high-speed photography. May’s work indicated that Froude scaling is a good first approximation to use when describing cavity behaviours such as deep closure. New et al. [7] investigated the impact and hydrodynamic characteristics of water entry associated with prismatic bodies with different fore section shapes and entry angles. Transient loading on the fore section of the prismatic body was measured with a triaxial accelerometer. Splashes, cavity formation, and closure procedures were observed with high-speed photography and video recordings. Anghileri and Spizzica [8] studied the normal impact of a rigid sphere with a water surface using finite element methods and certain experiments to gain variations in sphere acceleration during impact in order to verify the applied numerical method. Their results were in good agreement with the experimental data. Engle and Lewis [9] compared the results of the numerical method for predicting hydrodynamic impacts with experimental results of symmetrical wedges with different initial impact velocities. Their studies catered for the accuracy and validity of the various methods. Aristoff et al. [10] performed an experimental and theoretical study by impacting spheres with different densities on a water surface. A theoretical model was developed that yielded simple expressions for pinch-off time and depth, as well as the volume of air entrained by the sphere. Their study demonstrated that the dynamics of the air cavity formed in the wake of a sphere falling through the water surface can be altered considerably by the density of the sphere.


## 1.2. Analytical Solution

Von Karman [11] developed analytical solutions for the water entry problem of two-dimensional wedge bodies. He theoretically presented formulas using asymptotic theory. Wagner [12] presented a more precise analytical model which took the free surface elevation of the water into consideration. He developed an asymptotic solution for water entry by two-dimensional projectiles with small local deadrise angles. Vorus [13] proposed a novel analytical model for the water entry problem of bodies. In his model, nonlinearity in the boundary conditions was considered but the geometrical nonlinearity of the impact problem was neglected. Miloh [14] obtained analytical solutions based on several mathematical and physical assumptions for the small-time slamming coefficient and wetting factor of a rigid spherical shape in a vertical and oblique water entry. Aristoff and Bush [15] presented the results of a combined experimental and theoretical investigation of the normal impact of hydrophobic spheres on a water surface by determining the shape of the air cavity behind the sphere, pinch-off time, and depth. Ghadimi et al. [16] employed an analytical method and a numerical solution to study the water entry problem of a circular section. They demonstrated a procedure for the step by step derivation of the analytical formulas and presented a comparison between the linear and nonlinear solutions. Tassin et al. [17] investigated the water impact problem of bodies by CFD code OpenFOAM. They studied the water entry and exit of a two-dimensional body with a time-varying shape by using a modified Logvinovich model during the entry stage and the Von Karman model during the exit stage. Yao et al. [18] developed a theoretical model to describe the evolution of the cavity shape, including the time evolution of the cavity on fixed locations and its location evolution at fixed times.


## 1.3. Numerical Simulations

Park et al. [19] presented a numerical method for calculating the impact and buckling forces of bodies after a quick water entry. By neglecting fluid viscosity, potential flow was assumed to be nonviscous. This assumption greatly reduced computational time while solution accuracy was preserved by considering nonlinear free surface conditions. They solved this issue using a source panel method, in which an object is divided into panels or elements and each panel can be fully floating, partially floating, or fully out of the water. The results were compared with laboratory data and valuable results were obtained. Battistin and Iafrati [20] conducted a numerical study on normal entry into water of a two-dimensional symmetrical or asymmetrical object with a desirable form. This study considered potential flow of an incompressible fluid by ignoring the effects of gravity and adhesion. In this work, simulation was done using the boundary element method, and the boundary conditions of the free surface were fully nonlinear. The focus of this investigation was on evaluating pressure distribution and overall hydrodynamic loads applied to the colliding body. They verified the results for cylinder, sphere, and cone forms. Korobkin and Ohkusu [21] performed an interesting study on the hydroelastic coupling of finite element models using Wagner’s theory for water impacts. The aim of the study was to measure direct coupling possibilities using finite element methods for the structure and Wagner’s theory for hydrodynamic loads of impact of an elastic object with the water’s surface. Kleefsman et al. [22] numerically investigated the dambreak problem and water entry problems. The numerical method was based on Navier-Stokes equations that describe the flow of an incompressible viscous fluid. The equations were discretized onto a fixed Cartesian grid using the finite volume method. Kim et al. [23] analyzed the water entry problem of symmetrical bodies using the particle hydrodynamic method. Yang and Qiu [24] investigated the two-dimensional water entry problems of symmetric and asymmetric wedges with various deadrise angles using the Constrained Interpolation Profile (CIP) method. The effect of air compressibility for small deadrise angles was also discussed in their work. Fairlie-Clarke and Tveitnes [25] conducted a study on impact of wedge-shaped parts with a water surface at a constant velocity. In their study, Computational Fluid Dynamic (CFD) analysis was used to determine the momentum of added mass, momentum of flow, and gravity’s effect during the entry of wedge-shaped objects at a constant velocity and a dead angle of 5 to 45 degrees into water. A numerical solution based on finite volume method was implemented using Fluent software to solve equations of mass conservation and momentum to obtain flow fields. Hydrodynamic forces were not significantly changed by the effect of gravity on the flow fields. Yang and Qiu [26] computed slamming forces on two- and three-dimensional bodies entering calm water with a constant velocity. The nonlinear water entry problem was governed by Navier-Stokes equations and solved using a CIP-based finite difference method on a fixed Cartesian grid. Wu [27] developed a high order boundary element method for the complex velocity potential problem. He focused on application of the water entry problem for a wedge with varying speeds. His method ensured not only the continuity of the potential at the nodes of each element but also the velocity, allowing it to be applied to a variety of velocity potential problems. The continuity of the velocity achieved herein was particularly important for this kind of nonlinear free surface flow problem, because when the time stepping method is used, the free surface is updated using the velocity obtained at each node and the accuracy of the velocity is therefore crucial. Ahmadzadeh et al. [28] developed a three-dimensional numerical model to simulate a transient free-falling sphere impacting a free water surface using the Coupled Eulerian-Lagrangian formulations included in the commercial software Abaqus with consideration for important parameters such as viscosity, compressibility inertia, and variable downward velocity. Their study was successful in determining the shape of the air cavity, displacement of sphere versus time, pinch-off time, and depth compared to the experimental results found by Aristoff et al. [10]. Nguyen et al. [29] studied the water impact of various three-dimensional geometries, namely, a hemisphere, two cones, and a free-falling wedge, with an implicit algorithm based on a dual-time pseudo-compressibility method. Flow fields of incompressible viscous fluids were solved using unsteady RANS equations. A second-order VOF interface tracking algorithm was developed in a generalized curvilinear coordinate system to track the interface between the two phases in the computational domain. Sensitivity analysis with respect to the spatial grid resolution and the time step was performed to assess the accuracy of the results. Free surface deformation, pressure coefficients, impact velocities, and vertical acceleration during impact were compared with available experimental data and asymptotic theories, showing good agreement at a reasonable computational cost. S. Kim and N. Kim [30] performed integrated dynamic modeling of supercavitating vehicles and established 6-DOF equations by defining hydrodynamic forces and moments. These simulations also demonstrated the instability of supercavitating vehicles caused by the asymmetric immersion depth of coupled horizontal or vertical fins. Nguyen et al. [31] used a moving Chimera grid method to predict the real-time motion of water entry bodies by combining the 6-DOF rigid body motion model and the numerical calculation of the multiphase flow field, by which the coupled effect of supercavity and moving body was obtained. Mirzaei et al. [32] held the idea that existing planning force models weaken nonlinear interactions among the solid, liquid, and gaseous phases and are often too simple and therefore incorrect. They established the 6-DOF motion model for supercavitating vehicles by redefining the spatial planning force and developed a hybrid control scheme based on online planning force identification, which increased the stabilization of supercavitating vehicles.In this study, the water entry problem of three-dimensional pounders with different geometric shapes including cube, cylinder, sphere, pyramid, and cone was investigated numerically using commercial finite element code Abaqus 6.14-2 and an explicit dynamic analysis method was employed to model fluid-structure interactions using a Coupled Eulerian-Lagrangian (CEL) formulation [33]. This investigation focused on the problem of rigid pounders touching the free surface of a stationary viscous fluid at time t=0 with a variable velocity and finding for t>0, the shape of cavity, pinch-off time, and depth of the rigid pounders. It should be noted that the fluid was assumed to be viscous and compressible. Numerical simulations of solid-fluid interactions were performed using the Coupled Eulerian-Lagrangian (CEL) formulation and the remeshing techniques available in Abaqus 6.14-2 software [34]. The pounder was considered as a rigid solid and its mesh was created in Lagrangian form. The water was regarded as a compressible fluid and its mesh was created in Eulerian form. After validating the numerical results with recent experimental data, the effect of the pounder shape on its deepwater displacement and velocity was investigated. Then the effect of increasing the pounder drop height from the water’s surface on displacement, submersion velocity, pinch-off time, and depth was investigated. The critical drop height of the pounder, which is indicative of the maximum effect of the pounders drop height from the water’s surface on deepwater movement, was finally determined.


## 2. Coupled Eulerian-Lagrangian (CEL) Method

Eulerian materials and Lagrangian elements can interact using the Eulerian-Lagrangian contact. Simulations that include this type of contact are called Coupled Eulerian-Lagrangian (CEL) analyses [28]. The CEL method in Abaqus provides engineers and scientists with the ability to simulate a class of problems where the interaction between structures and fluids is important. This capability does not rely on the coupling of multiple software products but instead solves the fluid-structure interaction (FSI) simultaneously within Abaqus. This powerful, easy to use feature of Abaqus/Explicit general contact enables fully coupled multiphysics simulations such as fluid-structure interactions. The attractiveness of this method is that the fluid is modeled in an Eulerian frame while the structure can be modeled in a Lagrangian frame as is typical for solid mechanics applications. Specifically, in CEL, the governing equation of the solid uses the conservation of momentum, while the governing equation for the fluid uses a general Navier-Stokes form [36]. The Eulerian mesh is typically a simple rectangular grid of elements extending well beyond the Eulerian material boundaries, giving the material space in which to move and deform. If any Eulerian material moves outside the Eulerian mesh, it is lost from the simulation. The Eulerian material boundary must, therefore, be computed during each time increment and it generally does not correspond to an element boundary. The interaction between the structure and the water is handled remarkably well by the solver using the general contact algorithm. In traditional Lagrangian analysis, nodes are fixed within the material, and elements deform as the material deforms. Lagrangian elements are always 100% full of a single material, so the material boundary coincides with the element boundary. By contrast, in Eulerian analysis, nodes are fixed in space, and material flows through elements that do not deform. Eulerian elements may not always be 100% full of material and many may be partially or completely void.The Eulerian implementation in Abaqus/Explicit is based on the Volume Of Fluid (VOF) method. In numerical analyses using this CEL method, the Eulerian material is tracked as it flows through the mesh by computing its Eulerian Volume Fraction (EVF). Each Eulerian element is designated a percentage, which represents the portion of that element filled with a material. If an Eulerian element is completely filled with a material, its EVF is 1 and if there is no material in the element, its EVF is 0. Eulerian elements may simultaneously contain more than one material. If the sum of all material volume fractions in an element is less than one, the remainder of the element is automatically filled with a void material that has neither mass nor strength. Contact between Eulerian materials and Lagrangian materials is enforced using a general contact that is based on a penalty contact method. The Lagrangian elements can move through the Eulerian mesh without resistance until they encounter an Eulerian element filled with material (EVF ≠ 0). A CEL method that attempts to capture the advantages of both the Lagrangian and the Eulerian method is implemented in Abaqus [33].


### 2.1. Time Integration Scheme

The CEL method implemented in Abaqus/Explicit uses an explicit time integration scheme. The central difference rule is employed for the solution of a nonlinear system of differential equations. The unknown solution for the next time step can be found directly from the solution of the previous time step, so no iteration is needed. Another advantage of choosing an explicit time integration scheme for the present problem is robustness regarding difficult contact conditions [33].


### 2.2. Formulation of Eulerian-Lagrangian Contact

By using the CEL method the contact between the Eulerian domain and Lagrangian domain is discretized using general contact, which is based on the penalty contact method. This method is less stringent compared to the kinematic contact method. The penalty contact method is used in Abaqus to handle contact in CEL analysis [33].


### 2.3. Governing Equations in Lagrangian Method

The three fundamental conservation equations are the conservation of mass, the conservation of momentum, and the conservation of energy. The Lagrangian formulation trivially enforces the conservation of mass by defining the density as the ratio of the element mass to the current volume. The conservation of energy, in the absence of heat conduction, is enforced locally at selected points within each element or control volume. Different finite elements and finite difference methods are therefore defined in terms of how they handle the conservation of momentum.Finite element methods use the weak form of the equation, which is also known as the principle of virtual work or the principle of virtual power within the solid mechanics community. In this method, the principle of virtual work is presented as(1)∫Vρx¨δx+σ:δεdV=∫VρbδxdV+∫SττδxdS,where ρ is density, σ is the stress, ε is strain, τ is shear stress, V is the volume of the domain, and S is the surface, which is divided into Sτ and Su, which is where traction and displacement boundary conditions are imposed, respectively. Equation (1) can be rewritten as (2)∫Vρx¨δxdV=Fext-Fint,where Fint represents internal force and Fext represents external force, which is expressed as (3)Fint=∫Vσ:δεdV,Fext=∫VρbδxdV+∫SττδxdS.In the finite element method, the coordinates, velocity, acceleration, and virtual displacement are interpolated from their values at the nodes using shape functions,NA, where A is the node number. They are expressed as (4)x=∑ANAxA,x˙=∑ANAx˙A,x¨=∑ANAx¨A,δx=∑ANAδxA.Substituting the relations from (2) into the inertial term in the weak form gives(5)∫Vρx¨δxdV=∫Vρ∑ANAx¨A∑BNBδxBdV=∑A∑B∫VρNANBdVx¨AδxB.The external force vector is also derived by a direct substitution of the appropriate relations from (6):(6)Fext=∫VρbδxdV+∫SττδxdS,(7)Fext=∫Vρb∑BNBδxBdV+∫Sττ∑BNBδxBdS,(8)Fext=∑B∫VρbNBdV+∫SττNBdSδxB.The internal force vector requires the evaluation of the virtual strain,δ. Again, using the relations from (9)(9)δεij=12∂xi∂xj+∂xj∂xi,(10)δεij=∑B12∂NB∂xjδxBi+∂NB∂xiδxBj.These relations may be put in matrix form, which defines theB matrix, B, as(11)εδ=Bδxand the internal force may be written as(12)Fint=∫Vσ:δεdV,Fint=∑B∫VBBσdVδxB,where BB is the section of the B matrix associated with node B. To simplify, the mass matrix is defined as (13)MABc=∫VρNANBdV,where Mc is called the consistent mass matrix because its derivation is consistent with the finite element method. All the terms in the weak form have now been written in terms of unknown coordinates and virtual displacements. Collecting the terms gives(14)∑Bx¨B-1MB∫VρbNBdV+∫SττNBdS-∫VBBσdVδxB=0.Since the virtual displacements are arbitrary, thex¨B value is in accordance with (15)x¨B=1MB∫VρbNBdV+∫SττNBdS-∫VBBσdV.


### 2.4. Governing Equations in Eulerian Method

For the description of fluid flows usually Eulerian formulation is employed, because we are usually interested in the properties of the flow at certain locations in the flow domain. In Eulerian formulation, an independent equation is used to express conservation of mass. Inertia forces also include displacement sentences that when present in numerical solution lead to a matrix with nonsymmetric coefficients that depends on calculated velocities. Eulerian formulation allows for the use of simple strain and stress rates and velocity calculations instead of displacement values (in the ultrasmall displacement analysis) [36, 37]. Linear viscous isotropic fluids are known as Newtonian fluids, which are by far the most important fluids for practical applications. Newtonian fluids are characterized by the following material law for the Cauchy stress tensor T:(16)Tij=μ∂vi∂xj+∂vj∂xi-23∂vk∂xkδij-pδij,where the velocity vector is vi with respect to Cartesian coordinate xi, the pressure is p, the dynamic viscosity is μ,and the Kronecker symbol is δij. With this the conservation laws for mass, momentum, and energy take the form of(17)∂p∂t+∂pvi∂xi=0,(18)∂pvi∂t+∂pvivj∂xj=∂∂xjμ∂vi∂xj+∂vj∂xi-23∂vk∂xkδij-∂p∂xj+ρfi,(19)∂pe∂t+∂pvie∂xi=μ∂vi∂xj∂vi∂xj+∂vj∂xi-23∂vi∂xi2-p∂vi∂xi+∂∂xik∂T∂xi+ρq,where ρ is the fluid density, e is the specific internal energy, andfi and q are external forces and heat sources, respectively. In energy balance (19) for the heat flux vector hi Fourier’s law with thermal conductivity κf is used; that is, a constant specific heat capacity is assumed and the work done by pressure and friction forces is neglected:(20)hi=-κf∂T∂xi.System (17)–(19) has to be completed by two equations of state of the form:(21)p=pρ,T,e=eρ,Twhich define the thermodynamic properties of the fluid, that is, the thermal and caloric equation of state, respectively. When defining the viscosity for the material, it is treated as a Newtonian fluid, meaning that the viscosity depends only on temperature. This is the default setting in Abaqus. Since the temperature is constant throughout the analysis, the viscosity also is constant.


### 2.5. Explicit Time-Dependent Solution

There are implicit and explicit algorithms for dynamic analysis. In the implicit method the Newmark method is used, the matrices need to be reversed, and more implementation time is required. However, this method enjoys high accuracy and its convergence to linear problems is guaranteed. In the explicit method, the matrices do not need be reversed and the analysis speed is much higher. Convergence in the explicit method is achieved if the time step is less than the smallest element divided by the speed of sound in the material [37].The solution was advanced in time using the central difference method, which is also called the Verlet or leap frog method. It is based on the central difference approximation for evaluating the derivative of a functionf(t) at t:(22)dftdt=ft+1/2Δt-ft-1/2ΔtΔt,where Δt is the time step.Using (15), acceleration in time t can be written as(23)x¨t=FtMB.Now using (22) for velocity-acceleration, velocity at time t+1/2Δt is obtained as (24)x˙t+1/2Δt=x˙t-1/2Δt+x¨tΔt.Then using (22) for velocity-displacement, displacement at time t+Δt is obtained as(25)xt+Δt=xt+x˙t+ΔtΔt.The central difference method has its step size limited by stability considerations like all explicit integration methods, and the stable time step size was limited to(26)Δt≤lc,where l is the length of the shortest element and c is the speed of sound.


## 2.1. Time Integration Scheme

The CEL method implemented in Abaqus/Explicit uses an explicit time integration scheme. The central difference rule is employed for the solution of a nonlinear system of differential equations. The unknown solution for the next time step can be found directly from the solution of the previous time step, so no iteration is needed. Another advantage of choosing an explicit time integration scheme for the present problem is robustness regarding difficult contact conditions [33].


## 2.2. Formulation of Eulerian-Lagrangian Contact

By using the CEL method the contact between the Eulerian domain and Lagrangian domain is discretized using general contact, which is based on the penalty contact method. This method is less stringent compared to the kinematic contact method. The penalty contact method is used in Abaqus to handle contact in CEL analysis [33].


## 2.3. Governing Equations in Lagrangian Method

The three fundamental conservation equations are the conservation of mass, the conservation of momentum, and the conservation of energy. The Lagrangian formulation trivially enforces the conservation of mass by defining the density as the ratio of the element mass to the current volume. The conservation of energy, in the absence of heat conduction, is enforced locally at selected points within each element or control volume. Different finite elements and finite difference methods are therefore defined in terms of how they handle the conservation of momentum.Finite element methods use the weak form of the equation, which is also known as the principle of virtual work or the principle of virtual power within the solid mechanics community. In this method, the principle of virtual work is presented as(1)∫Vρx¨δx+σ:δεdV=∫VρbδxdV+∫SττδxdS,where ρ is density, σ is the stress, ε is strain, τ is shear stress, V is the volume of the domain, and S is the surface, which is divided into Sτ and Su, which is where traction and displacement boundary conditions are imposed, respectively. Equation (1) can be rewritten as (2)∫Vρx¨δxdV=Fext-Fint,where Fint represents internal force and Fext represents external force, which is expressed as (3)Fint=∫Vσ:δεdV,Fext=∫VρbδxdV+∫SττδxdS.In the finite element method, the coordinates, velocity, acceleration, and virtual displacement are interpolated from their values at the nodes using shape functions,NA, where A is the node number. They are expressed as (4)x=∑ANAxA,x˙=∑ANAx˙A,x¨=∑ANAx¨A,δx=∑ANAδxA.Substituting the relations from (2) into the inertial term in the weak form gives(5)∫Vρx¨δxdV=∫Vρ∑ANAx¨A∑BNBδxBdV=∑A∑B∫VρNANBdVx¨AδxB.The external force vector is also derived by a direct substitution of the appropriate relations from (6):(6)Fext=∫VρbδxdV+∫SττδxdS,(7)Fext=∫Vρb∑BNBδxBdV+∫Sττ∑BNBδxBdS,(8)Fext=∑B∫VρbNBdV+∫SττNBdSδxB.The internal force vector requires the evaluation of the virtual strain,δ. Again, using the relations from (9)(9)δεij=12∂xi∂xj+∂xj∂xi,(10)δεij=∑B12∂NB∂xjδxBi+∂NB∂xiδxBj.These relations may be put in matrix form, which defines theB matrix, B, as(11)εδ=Bδxand the internal force may be written as(12)Fint=∫Vσ:δεdV,Fint=∑B∫VBBσdVδxB,where BB is the section of the B matrix associated with node B. To simplify, the mass matrix is defined as (13)MABc=∫VρNANBdV,where Mc is called the consistent mass matrix because its derivation is consistent with the finite element method. All the terms in the weak form have now been written in terms of unknown coordinates and virtual displacements. Collecting the terms gives(14)∑Bx¨B-1MB∫VρbNBdV+∫SττNBdS-∫VBBσdVδxB=0.Since the virtual displacements are arbitrary, thex¨B value is in accordance with (15)x¨B=1MB∫VρbNBdV+∫SττNBdS-∫VBBσdV.


## 2.4. Governing Equations in Eulerian Method

For the description of fluid flows usually Eulerian formulation is employed, because we are usually interested in the properties of the flow at certain locations in the flow domain. In Eulerian formulation, an independent equation is used to express conservation of mass. Inertia forces also include displacement sentences that when present in numerical solution lead to a matrix with nonsymmetric coefficients that depends on calculated velocities. Eulerian formulation allows for the use of simple strain and stress rates and velocity calculations instead of displacement values (in the ultrasmall displacement analysis) [36, 37]. Linear viscous isotropic fluids are known as Newtonian fluids, which are by far the most important fluids for practical applications. Newtonian fluids are characterized by the following material law for the Cauchy stress tensor T:(16)Tij=μ∂vi∂xj+∂vj∂xi-23∂vk∂xkδij-pδij,where the velocity vector is vi with respect to Cartesian coordinate xi, the pressure is p, the dynamic viscosity is μ,and the Kronecker symbol is δij. With this the conservation laws for mass, momentum, and energy take the form of(17)∂p∂t+∂pvi∂xi=0,(18)∂pvi∂t+∂pvivj∂xj=∂∂xjμ∂vi∂xj+∂vj∂xi-23∂vk∂xkδij-∂p∂xj+ρfi,(19)∂pe∂t+∂pvie∂xi=μ∂vi∂xj∂vi∂xj+∂vj∂xi-23∂vi∂xi2-p∂vi∂xi+∂∂xik∂T∂xi+ρq,where ρ is the fluid density, e is the specific internal energy, andfi and q are external forces and heat sources, respectively. In energy balance (19) for the heat flux vector hi Fourier’s law with thermal conductivity κf is used; that is, a constant specific heat capacity is assumed and the work done by pressure and friction forces is neglected:(20)hi=-κf∂T∂xi.System (17)–(19) has to be completed by two equations of state of the form:(21)p=pρ,T,e=eρ,Twhich define the thermodynamic properties of the fluid, that is, the thermal and caloric equation of state, respectively. When defining the viscosity for the material, it is treated as a Newtonian fluid, meaning that the viscosity depends only on temperature. This is the default setting in Abaqus. Since the temperature is constant throughout the analysis, the viscosity also is constant.


## 2.5. Explicit Time-Dependent Solution

There are implicit and explicit algorithms for dynamic analysis. In the implicit method the Newmark method is used, the matrices need to be reversed, and more implementation time is required. However, this method enjoys high accuracy and its convergence to linear problems is guaranteed. In the explicit method, the matrices do not need be reversed and the analysis speed is much higher. Convergence in the explicit method is achieved if the time step is less than the smallest element divided by the speed of sound in the material [37].The solution was advanced in time using the central difference method, which is also called the Verlet or leap frog method. It is based on the central difference approximation for evaluating the derivative of a functionf(t) at t:(22)dftdt=ft+1/2Δt-ft-1/2ΔtΔt,where Δt is the time step.Using (15), acceleration in time t can be written as(23)x¨t=FtMB.Now using (22) for velocity-acceleration, velocity at time t+1/2Δt is obtained as (24)x˙t+1/2Δt=x˙t-1/2Δt+x¨tΔt.Then using (22) for velocity-displacement, displacement at time t+Δt is obtained as(25)xt+Δt=xt+x˙t+ΔtΔt.The central difference method has its step size limited by stability considerations like all explicit integration methods, and the stable time step size was limited to(26)Δt≤lc,where l is the length of the shortest element and c is the speed of sound.


## 3. Finite Element Modeling

### 3.1. Equation of Energy and Hugoniot Curve

The equation for conservation of energy equates the increase in internal energy per unit mass,Em, to the rate at which work is being done by stresses and the rate at which heat is being added. In the absence of heat conduction, the energy equation can be written as (27)ρ∂Em∂t=p-pbν1ρ∂ρ∂t+s:e˙+ρQ˙,where p is the pressure stress defined as positive in compression, pbν is the pressure stress due to the bulk viscosity, s is the deviator stress tensor, e˙ is the deviator part of strain rate, and Q˙ is the heat rate per unit mass. According to Abaqus documentation, flow modeling of compressible fluid can be achieved using the linear Us-Up form of the Mie-Gruneisen equation of state. The equation of state assumes pressure as a function of the current density, ρ,and the internal energy per unit mass, Em, according to (28)p=fρ,Emwhich defines all the equilibrium states that can exist in a material. The internal energy can be eliminated from (27), to obtain a pressure (p) versus volume (V) relationship or, equivalently, a p versus 1/ρ relationship that is unique to the material described by the equation of state model. This equation is uniquely dependent on the material defined by the equation of state. This unique relationship is called the Hugoniot curve and it is the locus of p-V states achievable behind a shock. The Hugoniot pressure, pH, is only a function of density and can be defined, in general, from fitting experimental data. An equation of state is said to be linear in energy when it can be written as (29)p=f+gEm,where f(ρ) and g(ρ) are only functions of density and depend on the particular equation of state model. Figure 1 illustrates a schematic expression of a Hugoniot curve [35].Figure 1
Schematic representation of a Hugoniot curve [35].


### 3.2. Mie-Gruneisen Equation of State

As mentioned previously, an equation of state is used to express the behaviour of Eulerian materials. For a Mie-Gruneisen equation of state for linear energy, the most common form is written as [38](30)p-pH=ΓρEm-EH,where pH and EH are the Hugoniot pressure and specific energy (per unit mass) and are only functions of density, and Γ is the Gruneisen ratio defined as (31)Γ=Γ0ρ0ρ,where Γ0 is a material constant and ρ0 is the reference density. Hugoniot energy, EH, is related to the Hugoniot pressure, pH, by (32)EH=pHη2ρ0,where η is the nominal volumetric compressive strain, written as (33)η=1-ρ0ρ.Elimination ofΓ and EH, from the above equations, yields (34)p=pH1-Γ0η2+Γ0ρ0Em.The equation of state and the energy equation represent coupled equations for pressure and internal energy. Abaqus solves these equations simultaneously at each material point using an explicit method.


### 3.3. LinearUs-Up Hugoniot Form

Normally, theUs-Up formulation of equation of state (EOS) is used to simulate shocks in solid materials. In this study it is used to define fluid materials. A common fit to the Hugoniot data is given by (35)pH=ρ0C02η1-sη2,where c0 and s define the linear relationship between the linear shock velocity, Us, and the particle velocity, Up, using (36)Us=c0+sUp.The Mie-Gruneisen EOS is normally used for processes with high pressures and high sound velocities. The relationship between density and pressure is governed by an EOS model. A linearUs-Up Hugoniot form of the Mie-Gruneisen equation of state is used to evaluate pressure in the water in the simulation. Using this implementation, the water is modeled as a compressible viscous fluid. With the above assumptions the linear Us-Up Hugoniot form is written as (37)p=ρ0C02η1-sη21-Γ0η2+Γ0ρ0Em,where ρ0c02 is equivalent to the elastic bulk modulus at small nominal strains.When modeling water, Abaqus prescribess=0 and Γ0=0. Filling this simplifies the relation between pressure p and density ρ to an inversely proportional equation and a simplification of the EOS. This simplification suffices for the modeling of water [39]:(38)ps=0∧Γ=0=ρ0C02η1-sη21-Γ0η2+Γ0ρ0Ems=0∧Γ=0=ρ0c02η.Recalling thatη=1-ρ0/ρ we get(39)p=ρ0C021-ρ0ρ.When defining a viscosity for the material, it is treated as a Newtonian fluid, meaning that viscosity only depends on the temperature. This is the default setting in Abaqus. Since the temperature is constant throughout the analysis, the viscosity is also constant.There is a limiting compression given by the denominator of this form of the equation of state. The denominator of (34) should not be equal to zero; therefore, an extreme value is defined for ρ and η as (40)ηlim=1s,ρlim=sρ0s-1.At this limit there is a tensile minimum; thereafter, negative sound speeds are calculated for the material.The linear Hugoniot formUs-Up equation of state can be employed to model laminar viscous flows governed by Navier-Stokes equations. The volumetric response is governed by the equations of state, where the bulk modulus acts as a penalty parameter for the constraint [38]. The Eulerian part (water) is simulated using the linear Us-Up Hugoniot form of the Mie-Gruneisen equation of state.The solving method in Eulerian equations is based on the division of the governing equations. In the Lagrangian method, the nodes are considered to be fixed relative to the material at any time interval and elements are deformed with the material, while in the Eulerian method element deformation is calculated at any time interval, but deformation is applied as a remeshing. The material is also displaced between each element and its neighboring elements at any time interval. Eulerian-Lagrangian contact allows the Eulerian materials to be combined with traditional nonlinear Lagrangian analyses. In the Eulerian-Lagrangian contact, the Lagrangian structure occupies void regions inside the Eulerian mesh. The contact algorithm automatically computes and tracks the interface between the Lagrangian structure and the Eulerian materials. A great benefit of this method is that there is no need to generate a conforming mesh for the Eulerian domain. In fact, a simple regular grid of Eulerian elements often yields the best accuracy. If the Lagrangian body is initially positioned inside the Eulerian mesh, it must be ensured that the underlying Eulerian elements contain voids after material initialization. During analysis the Lagrangian body pushes material out of the Eulerian elements that it passes through, and they become filled with void. Similarly, Eulerian material flowing towards the Lagrangian body is prevented from entering the underlying Eulerian elements. This formulation ensures that two materials never occupy the same physical space.In this study, CEL analysis allowed a direct coupling to solve this problem. The fact that CEL implementation uses the full Navier-Stokes equations implies that it was able to handle flow appropriately. It is important to note that CEL is not structured in the same way as a CFD code and was not intended to fill this role. While it does include viscous effects with extension to laminar flows, it does not include any turbulence effects, but in the context of an impact, which has a short duration and is highly transient (with turbulence effects being of less importance to the structural loading), CEL is well suited [28].


## 3.1. Equation of Energy and Hugoniot Curve

The equation for conservation of energy equates the increase in internal energy per unit mass,Em, to the rate at which work is being done by stresses and the rate at which heat is being added. In the absence of heat conduction, the energy equation can be written as (27)ρ∂Em∂t=p-pbν1ρ∂ρ∂t+s:e˙+ρQ˙,where p is the pressure stress defined as positive in compression, pbν is the pressure stress due to the bulk viscosity, s is the deviator stress tensor, e˙ is the deviator part of strain rate, and Q˙ is the heat rate per unit mass. According to Abaqus documentation, flow modeling of compressible fluid can be achieved using the linear Us-Up form of the Mie-Gruneisen equation of state. The equation of state assumes pressure as a function of the current density, ρ,and the internal energy per unit mass, Em, according to (28)p=fρ,Emwhich defines all the equilibrium states that can exist in a material. The internal energy can be eliminated from (27), to obtain a pressure (p) versus volume (V) relationship or, equivalently, a p versus 1/ρ relationship that is unique to the material described by the equation of state model. This equation is uniquely dependent on the material defined by the equation of state. This unique relationship is called the Hugoniot curve and it is the locus of p-V states achievable behind a shock. The Hugoniot pressure, pH, is only a function of density and can be defined, in general, from fitting experimental data. An equation of state is said to be linear in energy when it can be written as (29)p=f+gEm,where f(ρ) and g(ρ) are only functions of density and depend on the particular equation of state model. Figure 1 illustrates a schematic expression of a Hugoniot curve [35].Figure 1
Schematic representation of a Hugoniot curve [35].


## 3.2. Mie-Gruneisen Equation of State

As mentioned previously, an equation of state is used to express the behaviour of Eulerian materials. For a Mie-Gruneisen equation of state for linear energy, the most common form is written as [38](30)p-pH=ΓρEm-EH,where pH and EH are the Hugoniot pressure and specific energy (per unit mass) and are only functions of density, and Γ is the Gruneisen ratio defined as (31)Γ=Γ0ρ0ρ,where Γ0 is a material constant and ρ0 is the reference density. Hugoniot energy, EH, is related to the Hugoniot pressure, pH, by (32)EH=pHη2ρ0,where η is the nominal volumetric compressive strain, written as (33)η=1-ρ0ρ.Elimination ofΓ and EH, from the above equations, yields (34)p=pH1-Γ0η2+Γ0ρ0Em.The equation of state and the energy equation represent coupled equations for pressure and internal energy. Abaqus solves these equations simultaneously at each material point using an explicit method.


## 3.3. LinearUs-Up Hugoniot Form

Normally, theUs-Up formulation of equation of state (EOS) is used to simulate shocks in solid materials. In this study it is used to define fluid materials. A common fit to the Hugoniot data is given by (35)pH=ρ0C02η1-sη2,where c0 and s define the linear relationship between the linear shock velocity, Us, and the particle velocity, Up, using (36)Us=c0+sUp.The Mie-Gruneisen EOS is normally used for processes with high pressures and high sound velocities. The relationship between density and pressure is governed by an EOS model. A linearUs-Up Hugoniot form of the Mie-Gruneisen equation of state is used to evaluate pressure in the water in the simulation. Using this implementation, the water is modeled as a compressible viscous fluid. With the above assumptions the linear Us-Up Hugoniot form is written as (37)p=ρ0C02η1-sη21-Γ0η2+Γ0ρ0Em,where ρ0c02 is equivalent to the elastic bulk modulus at small nominal strains.When modeling water, Abaqus prescribess=0 and Γ0=0. Filling this simplifies the relation between pressure p and density ρ to an inversely proportional equation and a simplification of the EOS. This simplification suffices for the modeling of water [39]:(38)ps=0∧Γ=0=ρ0C02η1-sη21-Γ0η2+Γ0ρ0Ems=0∧Γ=0=ρ0c02η.Recalling thatη=1-ρ0/ρ we get(39)p=ρ0C021-ρ0ρ.When defining a viscosity for the material, it is treated as a Newtonian fluid, meaning that viscosity only depends on the temperature. This is the default setting in Abaqus. Since the temperature is constant throughout the analysis, the viscosity is also constant.There is a limiting compression given by the denominator of this form of the equation of state. The denominator of (34) should not be equal to zero; therefore, an extreme value is defined for ρ and η as (40)ηlim=1s,ρlim=sρ0s-1.At this limit there is a tensile minimum; thereafter, negative sound speeds are calculated for the material.The linear Hugoniot formUs-Up equation of state can be employed to model laminar viscous flows governed by Navier-Stokes equations. The volumetric response is governed by the equations of state, where the bulk modulus acts as a penalty parameter for the constraint [38]. The Eulerian part (water) is simulated using the linear Us-Up Hugoniot form of the Mie-Gruneisen equation of state.The solving method in Eulerian equations is based on the division of the governing equations. In the Lagrangian method, the nodes are considered to be fixed relative to the material at any time interval and elements are deformed with the material, while in the Eulerian method element deformation is calculated at any time interval, but deformation is applied as a remeshing. The material is also displaced between each element and its neighboring elements at any time interval. Eulerian-Lagrangian contact allows the Eulerian materials to be combined with traditional nonlinear Lagrangian analyses. In the Eulerian-Lagrangian contact, the Lagrangian structure occupies void regions inside the Eulerian mesh. The contact algorithm automatically computes and tracks the interface between the Lagrangian structure and the Eulerian materials. A great benefit of this method is that there is no need to generate a conforming mesh for the Eulerian domain. In fact, a simple regular grid of Eulerian elements often yields the best accuracy. If the Lagrangian body is initially positioned inside the Eulerian mesh, it must be ensured that the underlying Eulerian elements contain voids after material initialization. During analysis the Lagrangian body pushes material out of the Eulerian elements that it passes through, and they become filled with void. Similarly, Eulerian material flowing towards the Lagrangian body is prevented from entering the underlying Eulerian elements. This formulation ensures that two materials never occupy the same physical space.In this study, CEL analysis allowed a direct coupling to solve this problem. The fact that CEL implementation uses the full Navier-Stokes equations implies that it was able to handle flow appropriately. It is important to note that CEL is not structured in the same way as a CFD code and was not intended to fill this role. While it does include viscous effects with extension to laminar flows, it does not include any turbulence effects, but in the context of an impact, which has a short duration and is highly transient (with turbulence effects being of less importance to the structural loading), CEL is well suited [28].


## 4. Implementation of Coupled Eulerian-Lagrangian (CEL) Method in Abaqus

In order to simulate the water entry problem of the pounder and implement the Coupled Eulerian-Lagrangian (CEL) method in Abaqus, a geometric model for the pounder and water was first created. A partition was created to define water and void regions at the water surface. Since the fluid (water) undergoes a large deformation, an Eulerian grid is suitable. The rigid body (pounder) has small deformations relative to the fluid, making an Lagrangian grid suitable. Therefore, the water was modeled as Eulerian and the pounder was modeled as solid geometry with a rigid body constraint applied to it. Then, the material properties for the pounder and water were defined. Properties include density(ρ), Young’s modulus (E), and Poisson’s ratio (ν) for the pounder and density (ρw) and dynamic viscosity (η) for the water. In the Coupled Eulerian-Lagrangian method, when the pounder is considered as rigid, Young’s modulus and Poisson ratios values should be defined for the Abaqus, but they have no effect on the results. In the Coupled Eulerian-Lagrangian method, Us-Up equation of state (EOS) was selected to define the water properties in Abaqus. The Us-Up equation of state is described in detail in Abaqus 6.14 documentation, and it is used to simulate Navier-Stokes flow when turbulence flow is negligible. Because the major force imposed on the rigid body in water entry problems is pressure force, this equation of state can be used as an appropriate approximation. Parameters such as C0 (speed of sound in water), s (constant coefficient in shock velocity equation), and Γ0 (material constant in Gruneisen equation of state) were defined for the Us-Up equation of state. In the next step, the pounder and water were entered into the assembly section and the pounder was transferred to the waters center tangent to its surface. In order to simulate a free-falling rigid body into water, the pounder was first considered in a position tangent to the water surface; then the secondary velocity of pounder at the moment of impact was assigned. In the next step, explicit dynamic analysis was selected and the time period for pounder movement between entering the water and reaching the model bed was considered as 0.5 seconds. Nlgeom was also activated to control the nonlinear effects of large displacements. Then, by defining a reference point and assigning it to the pounder, required outputs such as displacement and velocity values for the vertical direction were selected. Based on its defined algorithm, Abaqus considers the maximum hardness between the Lagrangian elements (pounder) and Eulerian material (water). The Eulerian solution method uses the meshes defined for pounder movement between water meshes, and Abaqus detects the boundary between the Lagrangian and Eulerian materials. Considering the maximum hardness between the Lagrangian and Eulerian boundaries, the contact between these two materials is established using the penalty method. Concentrated mass was assigned to the pounder and its moment of inertia was defined for X, Y, and Z. In the next step, the pounder was defined as a rigid body that can move freely in all directions. Gravitational acceleration (Gravity) was defined as −9.8 m/s2 in the Y direction, and a nonreflecting Eulerian boundary was assigned to the sides and bottom surfaces of the computational water domain. This boundary condition is specific to the Eulerian method and prevents wave propagation due to rigid body impact on the water’s surface from returning to the computational domain after reaching the boundaries. With these boundary conditions the water is regarded as infinite. Assuming that the pounder is dropped from a height H above a free water surface with an initial velocity (V0) of 0 m/s, its secondary velocity (V) at the moment of impact was calculated according to the free fall equation of bodies. This velocity was vertically applied to the pounder. Eulerian and void regions were specified using material assignment and the mesh size was defined for the water and pounder. An Eulerian mesh was selected as the water mesh and an explicit mesh was selected as the pounder mesh using structured techniques. Finally, the problem was solved and the results were analyzed.


### 4.1. Model Geometry

Before applying the model to pounder water entry, a three-dimensional numerical model was developed to simulate a transient free-falling wedge impacting the free surface of water, in order to check and calibrate the adopted computational model for simpler problems. Then, the accuracy of the numerical model was examined by comparing the numerical results with the experimental data of Aristoff et al. [10] for the water entry of a sphere. Figure 2 illustrates a schematic for the computational domain, boundary conditions, and grid distribution for the sphere water entry problem. The dimensions for the computational domain are considered to be large enough, so that with further increases no considerable changes are observed in the numerical results. The sphere was modeled as solid geometry and a rigid body constraint was subsequently applied to it. The fluid domain was modeled as a cube with dimensions of 30 × 50 × 60 cm3. Mesh refinement was performed after assigning it to the computational domain; therefore the Eulerian mesh grid size was gradually decreased until no considerable changes were observed in the numerical results. The Eulerian domain (water) was modeled using quarter symmetry and contains 746172 elements. The sphere was also modeled as a Lagrangian solid using 15000 elements. The Eulerian domain was divided into two parts including the upper and lower regions. These regions were defined as void and stationary water, respectively. A solid sphere with a diameter of 2.54 cm and an initial velocity of 2.17 m/s in the normal direction was placed on the free water surface. The sphere can move freely in all directions with six degrees of freedom. The physical properties of water used in the numerical simulations are given in Table 1.Table 1
Equation of state parameters for water.
Density(ρw) (kg/m3)
Dynamic viscosity(η) (kgs/m2)
Speed of sound(c0) (m/sec)
Γ
0
s
1000
0.001
1450
0
0Figure 2
The computational domain, boundary conditions, and grid distribution for the water entry problem of sphere.Γ
0 and s are material constants in the Gruneisen equation of state and the constant coefficient in shock velocity equation, respectively.


## 4.1. Model Geometry

Before applying the model to pounder water entry, a three-dimensional numerical model was developed to simulate a transient free-falling wedge impacting the free surface of water, in order to check and calibrate the adopted computational model for simpler problems. Then, the accuracy of the numerical model was examined by comparing the numerical results with the experimental data of Aristoff et al. [10] for the water entry of a sphere. Figure 2 illustrates a schematic for the computational domain, boundary conditions, and grid distribution for the sphere water entry problem. The dimensions for the computational domain are considered to be large enough, so that with further increases no considerable changes are observed in the numerical results. The sphere was modeled as solid geometry and a rigid body constraint was subsequently applied to it. The fluid domain was modeled as a cube with dimensions of 30 × 50 × 60 cm3. Mesh refinement was performed after assigning it to the computational domain; therefore the Eulerian mesh grid size was gradually decreased until no considerable changes were observed in the numerical results. The Eulerian domain (water) was modeled using quarter symmetry and contains 746172 elements. The sphere was also modeled as a Lagrangian solid using 15000 elements. The Eulerian domain was divided into two parts including the upper and lower regions. These regions were defined as void and stationary water, respectively. A solid sphere with a diameter of 2.54 cm and an initial velocity of 2.17 m/s in the normal direction was placed on the free water surface. The sphere can move freely in all directions with six degrees of freedom. The physical properties of water used in the numerical simulations are given in Table 1.Table 1
Equation of state parameters for water.
Density(ρw) (kg/m3)
Dynamic viscosity(η) (kgs/m2)
Speed of sound(c0) (m/sec)
Γ
0
s
1000
0.001
1450
0
0Figure 2
The computational domain, boundary conditions, and grid distribution for the water entry problem of sphere.Γ
0 and s are material constants in the Gruneisen equation of state and the constant coefficient in shock velocity equation, respectively.


## 5. Validation

A schematic of the experimental apparatus of Aristoff et al. [10] is presented in Figure 3. A sphere was held by a camera like aperture at height H above a water tank. The tank had dimensions of 30 × 50 × 60 cm3 and was illuminated by a bank of twenty 32 W fluorescent bulbs. A diffuser was used to provide uniform lighting, and care was taken to keep the water surface free of dust. Tank dimensions were selected so that the walls had negligible effects on the results. The sphere was released from rest and fell towards the water, reaching it with an approximate speed V≈2gh. The impact sequence was recorded at 2000 frames per second using a high-speed camera. The resolution was set to 524 × 1280 pixels (px) with a field of view of 11.28 × 27.55 cm2, yielding a 46.46 px/cm magnification. The trajectory of the sphere and its impact speed were determined with subpixel accuracy through a cross correlation and Gaussian peak-fitting method yielding position estimates accurate to 0.025 px (0.0005 cm) and impact speeds accurate to ±4%.Figure 3
Schematic of the experimental apparatus of Aristoff et al. [10].Figure4 illustrates a comparison of cavity shapes between numerical simulation and the experimental photographs taken by Aristoff et al. [10] for a density ratio of 7.86 and an similar period of time. The impact velocity of the sphere was 2.17 m/s and t=0 refers to when the spheres center passes through the free surface. This figure shows the cavity shapes for both experimental and numerical simulation for similar time periods, except that cavity shape for the pinch-off time of the numerical simulation was added to the pictures.Figure 4
Comparison of numerical results with experimental photographs of Aristoff et al. [10] for the sphere water entry problem.An symmetric air cavity formed behind the sphere after entering the water. Air cavity formation consists of several steps including cavity creation, cavity expansion, cavity contraction behind the sphere, and cavity collapse. As the sphere moves deeper through the water, it exerts a force on surrounding fluid in radial directions and transfers its momentum to the fluid. However, fluid expansion is confronted with fluid hydrostatic pressure resistance and as a result the radial flow direction is reversed, leading to cavity contraction and collapse. Cavity collapse occurs when the cavity wall moves inward, until the moment of pinch-off (t = 70.8 ms), and the cavity is divided into two separated parts [34]. While the lower cavity sticks to the sphere and moves along with it, the upper cavity continues its contraction and moves towards the water’s surface. As shown in Figure 4, the numerical results and experimental observations are in good agreement.Abaqus numerical simulation methods efficiently model water splashes and jet flow. Pictures of these phenomena are depicted in Figure4. It was concluded from this figure that pinch-off in the numerical simulation happens with a delay time compared to the experimental method. Increased numerical errors during solution time are one of the causes for differences between experimental and numerical results.In Figure5, the numerical results of this study are compared with the theoretical and experimental data of Aristoff et al. as a displacement-time graph of a sphere in different water depths [10]. As can be seen, the numerical results are in good agreement with those of the analytics and experiments. However, the presently computed sphere center depth is lower than that found experimentally. It seems that, due to lack of turbulent modeling ability in Abaqus Eulerian formulation and evidently increasing turbulence in flow by time, there is just less difference between experiment and numerical results in the early time of simulation.Figure 5
Comparison of displacement-time graph of sphere in water depth with the theoretical and experimental data of Aristoff et al. [10].


## 6. Materials and Methods

### 6.1. Eulerian Model Definition

Figure6 illustrates a schematic of the computational domain, boundary conditions, and grid distribution for the water entry problem of a spherical pounder. The dimensions of the computational domain are considered to be large enough, so that with further increases no considerable changes are observed in the numerical results. The Eulerian domain was divided into the upper and lower regions. These regions are defined as void and stationary water, respectively. In water entry problems for rigid bodies, fluid domain dimensions must be at least eight times the dimensions of the rigid body. Hence, in this study, the fluid domain was modeled as a cube with dimensions of 2 × 2 × 1 m3. The pounders were placed exactly tangent to the water surface and their velocity at the moment of water surface impact was determined according to the free fall equation of bodies. The pounders can move freely in all directions with six degrees of freedom. Us-Up equation of state was used to define water material in the Coupled Eulerian-Lagrangian method. This equation of state is useful for Navier-Stokes flow simulations when turbulence flow is negligible. In cases where the impact of a rigid body on a water surface is investigated, the major force imposed on water body is the pressure force. Therefore, with appropriate approximations, this equation of state can be used to investigate the impact of rigid bodies on water surfaces.Figure 6
The computational domain, boundary conditions, and grid distribution for water entry problem of spherical pounder.


### 6.2. Eulerian Mesh Size

As shown in Figure7, with decreases in the Eulerian domain mesh size to less than 1.5 cm, the displacement-time graph of the spherical pounder shows negligible changes. Therefore, to improve analysis time, the Eulerian domain was modeled using quarter symmetry with a mesh size of 1.5 cm and contains 1538943 elements.Figure 7
Comparison of the effect of the Eulerian mesh size on displacement-time graph of spherical pounder.


### 6.3. Lagrangian Model Definition

Figure8 illustrates pounder dimensions for different geometric shapes. These pounders are identical in terms of material, mass density, Young’s modulus, Poisson’s ratio, volume, mass, and drop height from the water surface. The only difference is their geometric shape. These parameters were regarded as identical in all pounders to investigate the effect of pounder shape on deepwater displacement and velocity, air cavity shape, pinch-off time, and depth and the effect of water drag force on velocity depreciation when impacting the model bed.Figure 8
The dimensions of pounders with different geometric shapes.Pounder dimensions were selected so that they had a volume of 523 cm3 and equal mass due to their identical material and mass density. Table 2 illustrates the moment of inertia for pounders with different geometric shapes.Table 2
Moment of inertia of pounders with different geometric shapes.
Geometric shape of the pounder
Cube
Cylinder
Sphere
Pyramid
Cone
Moment of inertia (cm4)
68.5
51.37
41.1
30.82
41.1Table3 illustrates the material properties of the pounders. These parameters are regarded as identical for all pounders with the mentioned geometric shapes.Table 3
Material properties of pounders.
Material type
Mass (kg)
Densityρ (kg/m3)
Young’s modulus(E) (kg/m2)
Poisson’s ratioυ
Steel
4.11
7850
2.1 × 1010
0.3


### 6.4. Lagrangian Mesh Size

Due to the negligible effect of Lagrangian mesh size on the results and to improve analysis time, the Lagrangian mesh size for all pounders was set to 0.5 cm.


### 6.5. Boundary Condition

Wall boundary conditions were applied to the sides and bottom of the Eulerian domain to simulate a water tank. Pounder displacement was not constrained in any direction and the pounders could freely move in all directions. The sides and bottom boundaries of the three-dimensional Eulerian domain were set as nonreflecting boundaries, so that the waves released due to pounder impact on water surface did not return to the computational domain after reaching these boundaries and being absorbed by them. This boundary condition is written as [37](41)dp-pcdu=0,dxdt=-c,where ρ is density, c is the speed of sound in the fluid, p is pressure, u is velocity perpendicular to the wave, and x is the direction perpendicular to the boundary.


### 6.6. Drop Height and Impact Velocity of Pounders

After dropping the pounders from a height(h) above the water surface, the velocity of the pounders at the moment of impact was determined using (42)V2-V02=2gh,where V0 is the pounder initial velocity before being dropped from height h above the water’s surface, V is pounder secondary velocity when impacting the water surface, and g is the gravitational acceleration of the Earth.From (42), which expresses a time independent equation for the free fall of objects, the initial velocity and the terminal velocity of the object, drop height, and the gravitational acceleration of the Earth are interdependent in the absence of time. In this numerical simulation, the pounders were dropped from a height of 1 m above the water’s surface with an initial velocity (V0) equal to 0 m/s. Therefore, the secondary velocity of the pounders when impacting the water surface was obtained using (43)V=√2gh.Instead of simulating a full dropping event from the initial position, the pounders were placed exactly tangent to water’s surface and their secondary velocity at the moment of impact was determined according to (43). A gravitational acceleration (gravity) of −9.8 m/s2 was applied in the normal direction of the pounders. Therefore, a secondary velocity of 4.43 m/s in the normal direction of the water’s surface corresponds to the speed that would be obtained for the pounders.


## 6.1. Eulerian Model Definition

Figure6 illustrates a schematic of the computational domain, boundary conditions, and grid distribution for the water entry problem of a spherical pounder. The dimensions of the computational domain are considered to be large enough, so that with further increases no considerable changes are observed in the numerical results. The Eulerian domain was divided into the upper and lower regions. These regions are defined as void and stationary water, respectively. In water entry problems for rigid bodies, fluid domain dimensions must be at least eight times the dimensions of the rigid body. Hence, in this study, the fluid domain was modeled as a cube with dimensions of 2 × 2 × 1 m3. The pounders were placed exactly tangent to the water surface and their velocity at the moment of water surface impact was determined according to the free fall equation of bodies. The pounders can move freely in all directions with six degrees of freedom. Us-Up equation of state was used to define water material in the Coupled Eulerian-Lagrangian method. This equation of state is useful for Navier-Stokes flow simulations when turbulence flow is negligible. In cases where the impact of a rigid body on a water surface is investigated, the major force imposed on water body is the pressure force. Therefore, with appropriate approximations, this equation of state can be used to investigate the impact of rigid bodies on water surfaces.Figure 6
The computational domain, boundary conditions, and grid distribution for water entry problem of spherical pounder.


## 6.2. Eulerian Mesh Size

As shown in Figure7, with decreases in the Eulerian domain mesh size to less than 1.5 cm, the displacement-time graph of the spherical pounder shows negligible changes. Therefore, to improve analysis time, the Eulerian domain was modeled using quarter symmetry with a mesh size of 1.5 cm and contains 1538943 elements.Figure 7
Comparison of the effect of the Eulerian mesh size on displacement-time graph of spherical pounder.


## 6.3. Lagrangian Model Definition

Figure8 illustrates pounder dimensions for different geometric shapes. These pounders are identical in terms of material, mass density, Young’s modulus, Poisson’s ratio, volume, mass, and drop height from the water surface. The only difference is their geometric shape. These parameters were regarded as identical in all pounders to investigate the effect of pounder shape on deepwater displacement and velocity, air cavity shape, pinch-off time, and depth and the effect of water drag force on velocity depreciation when impacting the model bed.Figure 8
The dimensions of pounders with different geometric shapes.Pounder dimensions were selected so that they had a volume of 523 cm3 and equal mass due to their identical material and mass density. Table 2 illustrates the moment of inertia for pounders with different geometric shapes.Table 2
Moment of inertia of pounders with different geometric shapes.
Geometric shape of the pounder
Cube
Cylinder
Sphere
Pyramid
Cone
Moment of inertia (cm4)
68.5
51.37
41.1
30.82
41.1Table3 illustrates the material properties of the pounders. These parameters are regarded as identical for all pounders with the mentioned geometric shapes.Table 3
Material properties of pounders.
Material type
Mass (kg)
Densityρ (kg/m3)
Young’s modulus(E) (kg/m2)
Poisson’s ratioυ
Steel
4.11
7850
2.1 × 1010
0.3


## 6.4. Lagrangian Mesh Size

Due to the negligible effect of Lagrangian mesh size on the results and to improve analysis time, the Lagrangian mesh size for all pounders was set to 0.5 cm.


## 6.5. Boundary Condition

Wall boundary conditions were applied to the sides and bottom of the Eulerian domain to simulate a water tank. Pounder displacement was not constrained in any direction and the pounders could freely move in all directions. The sides and bottom boundaries of the three-dimensional Eulerian domain were set as nonreflecting boundaries, so that the waves released due to pounder impact on water surface did not return to the computational domain after reaching these boundaries and being absorbed by them. This boundary condition is written as [37](41)dp-pcdu=0,dxdt=-c,where ρ is density, c is the speed of sound in the fluid, p is pressure, u is velocity perpendicular to the wave, and x is the direction perpendicular to the boundary.


## 6.6. Drop Height and Impact Velocity of Pounders

After dropping the pounders from a height(h) above the water surface, the velocity of the pounders at the moment of impact was determined using (42)V2-V02=2gh,where V0 is the pounder initial velocity before being dropped from height h above the water’s surface, V is pounder secondary velocity when impacting the water surface, and g is the gravitational acceleration of the Earth.From (42), which expresses a time independent equation for the free fall of objects, the initial velocity and the terminal velocity of the object, drop height, and the gravitational acceleration of the Earth are interdependent in the absence of time. In this numerical simulation, the pounders were dropped from a height of 1 m above the water’s surface with an initial velocity (V0) equal to 0 m/s. Therefore, the secondary velocity of the pounders when impacting the water surface was obtained using (43)V=√2gh.Instead of simulating a full dropping event from the initial position, the pounders were placed exactly tangent to water’s surface and their secondary velocity at the moment of impact was determined according to (43). A gravitational acceleration (gravity) of −9.8 m/s2 was applied in the normal direction of the pounders. Therefore, a secondary velocity of 4.43 m/s in the normal direction of the water’s surface corresponds to the speed that would be obtained for the pounders.


## 7. Results and Discussion

### 7.1. Effect of the Geometric Shape of the Pounder on Deepwater Displacement

Figure9 illustrates the displacement-time graph of pounders with different geometric shapes. The graph continues for the cubical pounder from the moment of water surface impact at time zero up to the depth of 1 m at 415 ms. The displacement-time graph of the cylindrical pounder indicates that after entering the water it reaches the model bed in a shorter amount of time than the cubical pounder. The graph for this pounder continues from the moment of water surface impact at time zero up to the depth of 1 m at 355 ms. The linear displacement-time graph for the spherical and pyramidal pounders indicates the faster deepwater movement of these pounders compared to the cubical and cylindrical pounders. The deepwater movement time of the spherical and pyramidal pounders after entering the water to the model bed is 285 ms and 265 ms, respectively. The linear displacement-time graph of the conical pounder indicates that it has a faster deepwater movement compared to other pounders as it takes 235 ms to reach the model bed after entering the water.Figure 9
Displacement-time graph of pounders with different geometric shapes.


### 7.2. Mechanism of Air Cavity Formation and Pinch-Off due to Deepwater Movement of Pounders with Different Geometric Shapes

Figure10 illustrates the movement trajectory of pounders with different geometric shapes, from the moment of water impact to when the pounders reached the model bed. After entering the water, a symmetric air cavity is formed behind the pounders. Air cavity formation involves several steps including air cavity creation and expansion behind the pounder, air cavity contraction, and air cavity collapse [34]. As the pounder moves downward in water depth, it imposes a force onto surrounding fluid in its radial directions and transfers its momentum to the fluid. This extension is faced with the hydrostatic pressure resistance of fluid. The direction of the radial flow is reversed and eventually leads to air cavity contraction and collapse. In the other words, the air cavity contraction is accelerated until pinch-off occurs and the air cavity is divided into two distinct parts. The upper part of the air cavity continues its contraction and moves towards the water’s surface, whereas the lower part of the air cavity clings to the pounder and moves along with it [34]. After impacting the water surface, the cubical and cylindrical pounders due to their flat contact surfaces lose some of their velocity. This velocity depreciation slows down and continues until pinch-off. The spherical pounder with its curved contact surface has less velocity depreciation compared to the cubical and cylindrical pounders. However, some pounder velocity is always depreciated due to impacting the water surface. The velocity depreciation process for the spherical pounder during deepwater movement continues until pinch-off. After impacting the water’s surface, the pyramidal pounder with its pointy contact surface has a slight velocity depreciation but continues to move deeper into the water with approximately the same secondary velocity until pinch-off. By imposing the drag force of water onto the pounders surface, its velocity is gradually decreased until it reaches the model bed. After impacting the water surface, the conical pounder with its pointy contact surface continues to move deeper into the water with an increasing velocity. The increasing velocity of the conical pounder during deepwater movement continues even after pinch-off. When the pounder moves towards the model bed, the waters drag force is imposed onto pounder surface, causing velocity depreciation. The pinch-off time for the pounders with the above-mentioned shapes is between 225 and 250 ms.Figure 10
Air cavity formation due to deepwater movement of pounders with different geometric shapes and pinch-off time and depth.


### 7.3. Effect of the Pounder Geometric Shape on Deepwater Velocity

Figure11 illustrates the velocity-time graph for pounders with different geometric shapes. The cubical and cylindrical pounders lose some of their secondary velocity at the early moments of water entry. Hence, the velocity-time graph of these pounders illustrates great velocity depreciation up to 15 ms. The velocity reduction of the cubical and cylindrical pounders gradually continues until pinch-off. The velocity of the spherical pounder decreases up to 45 ms after water entry, but its velocity depreciation is less than the cubical and cylindrical pounders. This velocity reduction continues at a slower rate until pinch-off. In the initial moments of water entry, the pyramidal pounder has a slight reduction in secondary velocity up to 55 ms. Then velocity depreciation continues at a faster rate due to the drag force imposed onto pounders surface, and the velocity-time graph of the pyramidal pounder shows a gradual velocity reduction until pinch-off. After impacting the water surface, the conical pounder continues its deepwater movement with an increasing velocity rate. Due to its very low slope, the velocity-time graph indicates an increase in velocity for the conical pounder; but due to the waters drag force its velocity begins to decrease at 145 ms. The velocity-time graph indicates a gradual decrease in velocity for conical pounder, and velocity depreciation continues until pinch-off. After pinch-off, pounder velocities gradually decrease until they reach the model bed. Eventually, with a velocity of 3.88 m/s and 1.57 m/s at the moment of model bed impact, the conical and cubical pounders have the highest and lowest velocity, respectively. Water drag force is imposed in the opposite direction of pounder movement, and it has the most and least influence on the cubical and conical pounders, respectively.Figure 11
Velocity-time graph of pounders with different geometric shapes.


## 7.1. Effect of the Geometric Shape of the Pounder on Deepwater Displacement

Figure9 illustrates the displacement-time graph of pounders with different geometric shapes. The graph continues for the cubical pounder from the moment of water surface impact at time zero up to the depth of 1 m at 415 ms. The displacement-time graph of the cylindrical pounder indicates that after entering the water it reaches the model bed in a shorter amount of time than the cubical pounder. The graph for this pounder continues from the moment of water surface impact at time zero up to the depth of 1 m at 355 ms. The linear displacement-time graph for the spherical and pyramidal pounders indicates the faster deepwater movement of these pounders compared to the cubical and cylindrical pounders. The deepwater movement time of the spherical and pyramidal pounders after entering the water to the model bed is 285 ms and 265 ms, respectively. The linear displacement-time graph of the conical pounder indicates that it has a faster deepwater movement compared to other pounders as it takes 235 ms to reach the model bed after entering the water.Figure 9
Displacement-time graph of pounders with different geometric shapes.


## 7.2. Mechanism of Air Cavity Formation and Pinch-Off due to Deepwater Movement of Pounders with Different Geometric Shapes

Figure10 illustrates the movement trajectory of pounders with different geometric shapes, from the moment of water impact to when the pounders reached the model bed. After entering the water, a symmetric air cavity is formed behind the pounders. Air cavity formation involves several steps including air cavity creation and expansion behind the pounder, air cavity contraction, and air cavity collapse [34]. As the pounder moves downward in water depth, it imposes a force onto surrounding fluid in its radial directions and transfers its momentum to the fluid. This extension is faced with the hydrostatic pressure resistance of fluid. The direction of the radial flow is reversed and eventually leads to air cavity contraction and collapse. In the other words, the air cavity contraction is accelerated until pinch-off occurs and the air cavity is divided into two distinct parts. The upper part of the air cavity continues its contraction and moves towards the water’s surface, whereas the lower part of the air cavity clings to the pounder and moves along with it [34]. After impacting the water surface, the cubical and cylindrical pounders due to their flat contact surfaces lose some of their velocity. This velocity depreciation slows down and continues until pinch-off. The spherical pounder with its curved contact surface has less velocity depreciation compared to the cubical and cylindrical pounders. However, some pounder velocity is always depreciated due to impacting the water surface. The velocity depreciation process for the spherical pounder during deepwater movement continues until pinch-off. After impacting the water’s surface, the pyramidal pounder with its pointy contact surface has a slight velocity depreciation but continues to move deeper into the water with approximately the same secondary velocity until pinch-off. By imposing the drag force of water onto the pounders surface, its velocity is gradually decreased until it reaches the model bed. After impacting the water surface, the conical pounder with its pointy contact surface continues to move deeper into the water with an increasing velocity. The increasing velocity of the conical pounder during deepwater movement continues even after pinch-off. When the pounder moves towards the model bed, the waters drag force is imposed onto pounder surface, causing velocity depreciation. The pinch-off time for the pounders with the above-mentioned shapes is between 225 and 250 ms.Figure 10
Air cavity formation due to deepwater movement of pounders with different geometric shapes and pinch-off time and depth.


## 7.3. Effect of the Pounder Geometric Shape on Deepwater Velocity

Figure11 illustrates the velocity-time graph for pounders with different geometric shapes. The cubical and cylindrical pounders lose some of their secondary velocity at the early moments of water entry. Hence, the velocity-time graph of these pounders illustrates great velocity depreciation up to 15 ms. The velocity reduction of the cubical and cylindrical pounders gradually continues until pinch-off. The velocity of the spherical pounder decreases up to 45 ms after water entry, but its velocity depreciation is less than the cubical and cylindrical pounders. This velocity reduction continues at a slower rate until pinch-off. In the initial moments of water entry, the pyramidal pounder has a slight reduction in secondary velocity up to 55 ms. Then velocity depreciation continues at a faster rate due to the drag force imposed onto pounders surface, and the velocity-time graph of the pyramidal pounder shows a gradual velocity reduction until pinch-off. After impacting the water surface, the conical pounder continues its deepwater movement with an increasing velocity rate. Due to its very low slope, the velocity-time graph indicates an increase in velocity for the conical pounder; but due to the waters drag force its velocity begins to decrease at 145 ms. The velocity-time graph indicates a gradual decrease in velocity for conical pounder, and velocity depreciation continues until pinch-off. After pinch-off, pounder velocities gradually decrease until they reach the model bed. Eventually, with a velocity of 3.88 m/s and 1.57 m/s at the moment of model bed impact, the conical and cubical pounders have the highest and lowest velocity, respectively. Water drag force is imposed in the opposite direction of pounder movement, and it has the most and least influence on the cubical and conical pounders, respectively.Figure 11
Velocity-time graph of pounders with different geometric shapes.


## 8. Critical Drop Height of the Pounder

In the dynamic compaction of the seabed, determining the critical drop height of the pounder from the free water surface is an effective parameter for the pounder to reach maximum deepwater velocity. After dropping the pounder from heightH above the free water surface with an initial velocity (V0) of 0 m/s, its velocity gradually increases until it reaches the secondary velocity (V) at the moment of water impact. After entering the water, pounder velocity is decreased due to water drag force. According to the free fall equation for objects, increasing the drop height from the water’s surface leads to increased pounder velocity at the moment of water surface impact, but it does not necessarily lead to an increase in pounder velocity during seabed impact because the critical drop height of the pounder is a determining factor of its maximum deepwater velocity. According to Figure 12, a conical pounder with a mass of 4.11 kg was dropped from different heights (H), and its velocity when impacting the waters was measured according to the time independent equation for the free fall of objects. Table 4 illustrates the velocity of the conical pounder with different drop heights when impacting the water’s surface.Table 4
Velocity of the conical pounder with different drop heights at the moment of impacting the water surface.
Drop height of the pounder from water surface (m)
Pounder velocity at the moment of impacting water surface (m/s)
0.25
2.21
0.5
3.13
0.75
3.83
1
4.43
1.25
4.95
1.5
5.42
1.75
5.85
2
6.26
2.25
6.64
2.5
7Figure 12
Schematic of drop height of the conical pounder from water surface.


### 8.1. Effect of the Drop Height of the Pounder from Water Surface on Deepwater Displacement

Figure13 illustrates the displacement-time graph of the conical pounder with different drop heights from the water’s surface. As can be seen, increasing the drop height of the pounder leads to displacement-time graphs that are closer to each other. In fact, an increase in the drop height of the pounder from the water’s surface leads to an increase in its secondary velocity when impacting the water’s surface. Hence, the pounder enters the water with a higher velocity and reaches the model bed in less time.Figure 13
Displacement-time graph of the conical pounder with different drop heights.


### 8.2. Effect of the Drop Height of the Pounder from Water Surface on Its Deepwater Velocity

Increasing the drop height of the pounder from the water’s surface leads to an increase in its velocity before impacting the water; but the pounder velocity gradually decreases after entering the water due to the its drag force. Figure14 illustrates the velocity-time graph for the conical pounder with different drop heights from the water’s surface. When the pounder is dropped from a low height, due to its low velocity, it moves deepwater due to its weight. Moreover, because of the low velocity of the pounder when entering the water, pinch-off occurs in less time and at lower depths, and its velocity increases after pinch-off until it reaches terminal velocity. Increases in pounder drop heights increase water entry velocity. Hence, it will have greater velocity depreciation in the first moments of water entry.Figure 14
Velocity-time graph of the conical pounder with different drop heights.Figure15 illustrates the terminal velocity-time graph for the conical pounder with different drop heights when impacting the model bed. As can be seen, increasing the drop height of the pounder from the water’s surface leads to an increase in its terminal velocity when impacting the model bed. These conditions are dominant as long as the pounder reaches a critical drop height from the water surface, but increasing the pounder drop height higher than the critical drop height leads to decreased terminal velocity when impacting the model bed. In this numerical simulation, after dropping the pounder from a height of 2 m above the water’s surface, its secondary velocity was 6.26 m/s when impacting the water surface. After entering the water, the secondary velocity of the pounder decreased until it impacts the model bed with a terminal velocity of 4.54 m/s. Therefore, a drop height of 2 m from above the water’s surface has the highest terminal velocity when the pounder impacts the model bed. Increasing the drop height by more than 2 m leads to decreased pounder terminal velocity when impacting the model bed. The height which has the greatest effect on pounder terminal velocity when impacting the model bed is called the critical drop height. Increasing the drop height of the pounder to greater than its critical drop height leads to decreased terminal velocity when impacting the model bed.Figure 15
Terminal velocity of conical pounder with different drop heights versus time.Figure16 illustrates the secondary and terminal velocity of the conical pounder with different drop heights. As can be observed, increasing the drop height of the pounder from water’s surface leads to increased differences between the secondary and terminal velocities of the pounder. Increasing the drop height to more than the critical drop height leads to increased pounder secondary velocity when impacting the water’s surface and decreased terminal velocity when impacting the model bed. Hence, the difference between the secondary and terminal velocities of the pounder is not reasonable, and a significant amount of the pounders secondary velocity is depreciated when entering the water. The great difference between the secondary and terminal velocities of the pounder reveals the necessity of determining the critical drop height of the pounder before dynamic seabed compaction operations. Increasing the pounder drop height to be greater than the critical drop height not only has a reverse result on pounder velocity when impacting the seabed, but it also requires special equipment such as lengthy cranes and huge barges leading to increased implementation costs.Figure 16
Secondary and terminal velocities of conical pounder versus drop heights.Figure17 illustrates pinch-off time based on the secondary velocity of the pounder at the moment of water surface impact. As can be seen, pinch-off time is a weak function of pounder impact velocity and increasing impact velocity leads to a slight reduction in pinch-off time. This behaviour was also observed in the analytical and numerical results found by Lee et al. [40].Figure 17
Pinch-off time versus impact velocity of conical pounder.Figure18 illustrates pinch-off depth based on the secondary velocity of pounder at the moment of water surface impact. As can be seen, the pinch-off depth is a linear function of pounder impact velocity. These results are in line with the analytical equations of Lee et al. [40].Figure 18
Pinch-off depth versus impact velocity of conical pounder.


## 8.1. Effect of the Drop Height of the Pounder from Water Surface on Deepwater Displacement

Figure13 illustrates the displacement-time graph of the conical pounder with different drop heights from the water’s surface. As can be seen, increasing the drop height of the pounder leads to displacement-time graphs that are closer to each other. In fact, an increase in the drop height of the pounder from the water’s surface leads to an increase in its secondary velocity when impacting the water’s surface. Hence, the pounder enters the water with a higher velocity and reaches the model bed in less time.Figure 13
Displacement-time graph of the conical pounder with different drop heights.


## 8.2. Effect of the Drop Height of the Pounder from Water Surface on Its Deepwater Velocity

Increasing the drop height of the pounder from the water’s surface leads to an increase in its velocity before impacting the water; but the pounder velocity gradually decreases after entering the water due to the its drag force. Figure14 illustrates the velocity-time graph for the conical pounder with different drop heights from the water’s surface. When the pounder is dropped from a low height, due to its low velocity, it moves deepwater due to its weight. Moreover, because of the low velocity of the pounder when entering the water, pinch-off occurs in less time and at lower depths, and its velocity increases after pinch-off until it reaches terminal velocity. Increases in pounder drop heights increase water entry velocity. Hence, it will have greater velocity depreciation in the first moments of water entry.Figure 14
Velocity-time graph of the conical pounder with different drop heights.Figure15 illustrates the terminal velocity-time graph for the conical pounder with different drop heights when impacting the model bed. As can be seen, increasing the drop height of the pounder from the water’s surface leads to an increase in its terminal velocity when impacting the model bed. These conditions are dominant as long as the pounder reaches a critical drop height from the water surface, but increasing the pounder drop height higher than the critical drop height leads to decreased terminal velocity when impacting the model bed. In this numerical simulation, after dropping the pounder from a height of 2 m above the water’s surface, its secondary velocity was 6.26 m/s when impacting the water surface. After entering the water, the secondary velocity of the pounder decreased until it impacts the model bed with a terminal velocity of 4.54 m/s. Therefore, a drop height of 2 m from above the water’s surface has the highest terminal velocity when the pounder impacts the model bed. Increasing the drop height by more than 2 m leads to decreased pounder terminal velocity when impacting the model bed. The height which has the greatest effect on pounder terminal velocity when impacting the model bed is called the critical drop height. Increasing the drop height of the pounder to greater than its critical drop height leads to decreased terminal velocity when impacting the model bed.Figure 15
Terminal velocity of conical pounder with different drop heights versus time.Figure16 illustrates the secondary and terminal velocity of the conical pounder with different drop heights. As can be observed, increasing the drop height of the pounder from water’s surface leads to increased differences between the secondary and terminal velocities of the pounder. Increasing the drop height to more than the critical drop height leads to increased pounder secondary velocity when impacting the water’s surface and decreased terminal velocity when impacting the model bed. Hence, the difference between the secondary and terminal velocities of the pounder is not reasonable, and a significant amount of the pounders secondary velocity is depreciated when entering the water. The great difference between the secondary and terminal velocities of the pounder reveals the necessity of determining the critical drop height of the pounder before dynamic seabed compaction operations. Increasing the pounder drop height to be greater than the critical drop height not only has a reverse result on pounder velocity when impacting the seabed, but it also requires special equipment such as lengthy cranes and huge barges leading to increased implementation costs.Figure 16
Secondary and terminal velocities of conical pounder versus drop heights.Figure17 illustrates pinch-off time based on the secondary velocity of the pounder at the moment of water surface impact. As can be seen, pinch-off time is a weak function of pounder impact velocity and increasing impact velocity leads to a slight reduction in pinch-off time. This behaviour was also observed in the analytical and numerical results found by Lee et al. [40].Figure 17
Pinch-off time versus impact velocity of conical pounder.Figure18 illustrates pinch-off depth based on the secondary velocity of pounder at the moment of water surface impact. As can be seen, the pinch-off depth is a linear function of pounder impact velocity. These results are in line with the analytical equations of Lee et al. [40].Figure 18
Pinch-off depth versus impact velocity of conical pounder.


## 9. Conclusions

In this study, the water entry problem of three-dimensional pounders with different geometric shapes of cube, cylinder, sphere, pyramid, and cone was numerically simulated by commercial finite element code Abaqus 6.14-2. An explicit dynamic analysis method was employed to model the fluid-structure interaction using a Coupled Eulerian-Lagrangian (CEL) method. Before applying the model to water entry of the pounders, a three-dimensional numerical model was developed to simulate a transient free-falling wedge impacting the free surface of the water in order to calibrate the adopted computational model. The accuracy of the numerical model was examined by comparing numerical results with the experimental data of Aristoff et al. The good agreement between numerical simulation results with the theoretical and experimental data of Aristoff et al. demonstrates the capability of the employed numerical methodology to simulate the water entry problem for three-dimensional geometries with acceptable accuracy. After water entry by a pounder dropped from a height of 1 m, a symmetric air cavity is formed behind the pounder. As the pounder moves downward in water depth, it imposes a force onto surrounding fluid in the radial directions and transfers its momentum to the fluid. This extension is faced with the hydrostatic pressure resistance of the fluid, causing the direction of the radial flow to be reversed, leading to air cavity contraction and collapse. The results show that the water drag force had the highest and lowest effect on the cubical and conical pounders, respectively, and increasing the drop height to greater than the critical drop height leads to increased pounder velocity when impacting the water’s surface and decreased velocity when impacting the model bed. Pinch-off time is a very weak function of pounder velocity when impacting the water’s surface; but the pinch-off depth increases linearly along with increased pounder impact velocities. It is recommended that the critical drop height for the pounder is determined before dynamic seabed compaction operations and the adjustment of crane heights to this height. In addition, it is also recommended that perforated pounders be used in order to reduce velocity depreciation when impacting the water’s surface and to reduce the effect of water drag force.



---

*Source: 5980386-2018-02-06.xml*