Geofluids

Volume 2018 (2018), Article ID 5980386, 18 pages

https://doi.org/10.1155/2018/5980386

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

Correspondence should be addressed to Mohammad Reza Atrechian

Received 19 August 2017; Revised 17 October 2017; Accepted 3 January 2018; Published 6 February 2018

Academic Editor: Stefano Lo Russo

Copyright © 2018 Mohammad Hossein Taghizadeh Valdi et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

#### 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.

#### 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 with a variable velocity and finding for , 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 where is density, is the stress, is strain, is shear stress, is the volume of the domain, and is the surface, which is divided into and , which is where traction and displacement boundary conditions are imposed, respectively. Equation (1) can be rewritten as where represents internal force and represents external force, which is expressed as

In the finite element method, the coordinates, velocity, acceleration, and virtual displacement are interpolated from their values at the nodes using shape functions, , where is the node number. They are expressed as

Substituting the relations from (2) into the inertial term in the weak form gives

The external force vector is also derived by a direct substitution of the appropriate relations from (6):

The internal force vector requires the evaluation of the virtual strain, *δ*. Again, using the relations from (9)

These relations may be put in matrix form, which defines the matrix, , asand the internal force may be written aswhere is the section of the matrix associated with node . To simplify, the mass matrix is defined as where 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

Since the virtual displacements are arbitrary, the value is in accordance with

##### 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 :where the velocity vector is with respect to Cartesian coordinate , the pressure is , the dynamic viscosity is *, *and the Kronecker symbol is . With this the conservation laws for mass, momentum, and energy take the form ofwhere is the fluid density, is the specific internal energy, and and are external forces and heat sources, respectively. In energy balance (19) for the heat flux vector Fourier’s law with thermal conductivity is used; that is, a constant specific heat capacity is assumed and the work done by pressure and friction forces is neglected:

System (17)–(19) has to be completed by two equations of state of the form:which 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 function at :where is the time step.

Using (15), acceleration in time can be written as

Now using (22) for velocity-acceleration, velocity at time is obtained as

Then using (22) for velocity-displacement, displacement at time is obtained as

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 towhere is the length of the shortest element and 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, , 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 where is the pressure stress defined as positive in compression, is the pressure stress due to the bulk viscosity, is the deviator stress tensor, is the deviator part of strain rate, and is the heat rate per unit mass. According to Abaqus documentation, flow modeling of compressible fluid can be achieved using the linear - 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, , according to which defines all the equilibrium states that can exist in a material. The internal energy can be eliminated from (27), to obtain a pressure versus volume relationship or, equivalently, a versus 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 - states achievable behind a shock. The Hugoniot pressure, , 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 where and 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].