# Hydrodynamic impact: Numerical and experimental investigations

Bernard Peseux, Laurent Gornet, Bundy Donguy

► To cite this version:

Bernard Peseux, Laurent Gornet, Bundy Donguy. Hydrodynamic impact: Numerical and experimental investigations. Journal of Fluids and Structures, 2005, 21 (3), pp.277-303. 10.1016/j.jfluidstructs.2005.04.011. hal-01006935

HAL Id: hal-01006935

https://hal.science/hal-01006935

Submitted on 3 Feb 2017

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Hydrodynamic impact: Numerical and experimental investigations

B. Peseux, L. Gornet, B. Donguy

Institut de Recherches en Génie Civil et Mécanique, Equipe Structure et Simulation, Ecole Centrale de Nantes, 1, rue de la Noë,
44321 Nantes, Cedex 3, France

This paper deals with the slamming phenomenon experienced by ships during impact between the bow and the water free surface. Slamming loads on ships may be sufficiently important to induce plastic deformation of the hull external structure. For extreme loading cases, they have been identified as the cause of ship loss. The problem to be solved is transient and highly nonlinear due to the complex water flow conditions. In the present paper, the three-dimensional Wagner problem is solved numerically using the finite element method. A numerical analysis is performed for both rigid and deformable structures. After this numerical analysis, an original experimental investigation is presented. It consists of a series of free fall drop-tests of rigid and deformable cone-shaped samples with different deadrise angles and thickness. Distribution and evolution of pressure are analyzed. Finally, our numerical results are successfully compared with experimental data.

## 1. Introduction

Rough sea conditions can induce high-amplitude ship motions. Different factors, such as boat speed, heading and sea state, may cause the emergence of the ship bow from the water. In this context, slamming refers to the hydrodynamic impact between the water free surface and the hull as the bow re-enters the incoming wave. In this situation, the hull undergoes hydrodynamic loads that may induce localized plastic deformation. For extreme loading cases, it has been recognized to be responsible for ship loss. The resulting stress levels induced by swell waves can significantly reduce the life duration of the hull. Nowadays, interest in the subject has been revived by the advent of new types of high-speed vessels, such as hydrofoils, small waterplane area twin hull (SWATH), sea effect surface (SES) boats and hovercraft. The speed of these boats reaches 35-40 knots and results in particular ship motion behaviour.

An overview of these different types of fast boats exhibits a great variety of aspects. For boat design, two options can be considered: on the one hand, time-tested experimental techniques, whose advantages and disadvantages are well-known, can be used, and on the other hand, simulation approaches are more attractive, even though all their aspects are not completely controlled. Then, numerical tools must be developed not only to calculate ship dimensions but also to provide comparative data in order to choose the best technological solutions. Thus, boat design can be improved by considering slamming loads. For these types of problems, two approaches can be considered. First, the Navier-Stokes
equations can be solved in terms of velocity and pressure formulation, and second a velocity potential can be introduced under certain assumptions. Impact loads are commonly estimated within the framework of potential flow theory. This approach is one of the simplest ways to describe the violent motion of the free surface in the vicinity of a moving surface-piercing body. Despite simplifications, the resulting problem remains highly nonlinear, mainly due to the fact that free-surface boundary conditions must be applied on a boundary whose position is unknown.

Von Karman (1929) was one of the first who investigated the hydrodynamic impact problem by studying seaplane floats during landing. He introduced the concept of added apparent mass due to the water motion induced by the body. This work was followed by Wagner (1932) who considered piled-up water effects. The body was idealized as a flat plate, whose width is determined by the intersection between the elevation of the free surface and the body position. The resulting force is obtained by integrating the pressure distribution. It is noted that this distribution admits a singularity at the plate edges. To correct this singularity and describe the jet flow, a complementary problem should be formulated. Wagner's work was truly remarkable because it took place before the mathematical development of matched asymptotic expansion techniques. These methods are nowadays commonly used to approximate the solution of differential equations.

First-order asymptotic solutions were found for the case of a wedge with small deadrise angles (Wagner, 1932). Wanatabe (1986) corrected the singularity of Wagner's solution by considering an inner problem in the region of the spray source. Solutions for cylinders (Cointe and Armand, 1987), and more generally for arbitrary two-dimensional blunt body shapes (Cointe, 1989; Howison et al., 1991) were recently established. These two-dimensional asymptotic results were successfully validated by experiments in the cases of wedges (Chuang, 1967; Fontaine and Cointe, 1997), cylinders (Cointe and Armand, 1987), and realistic ship cross-sections (Zhao et al., 1996; Magee and Fontaine, 1998).

For three-dimensional problems, first-order analytical solutions have only been established for particular solid shapes (Korobkin, 1985; Scolan and Korobkin, 2001).

## 2. Hydrodynamic impact

In this section, assumptions and equations which lead to the general formulation of the hydrodynamic impact problem are recalled.

The structure which occupies the domain $\Omega_S$ with boundary $\partial\Omega_S$ comes into contact with a fluid domain $\Omega_F$, bounded or not. The fluid free-surface is denoted by $\Gamma_L$. In the case of an unbounded fluid, the problem is approximated by introducing an artificial boundary $\Gamma_\infty$ located sufficiently far from the structure to consider that its influence on both fluid and structure responses in the vicinity of the impact zone can be neglected. Under these assumptions, infinity boundary conditions can be written on this boundary.

During hydrodynamic impact, the fluid-structure interface (the wetted body surface), denoted as $\Gamma_{FS}$, belongs to both the structure boundary $\partial\Omega_S$ and the initial fluid free-surface $\Gamma_L$. The principal difficulty of the problem is that the interface $\Gamma_{FS}$ and the free-surface $\Gamma_L$ evolve during impact and consequently are unknown variables. The fluid and the structure are described in the same coordinate system $(O, \vec{i}, \vec{j}, \vec{k})$, the origin $O$ belonging to the initial free-surface of the fluid and the $\vec{k}$-axis pointing upward, as shown in Fig. 1.

### 2.1. Fluid problem formulation

For the hydrodynamic impact problem, the fluid description is carried out by adopting classical assumptions (Batchelor, 1967; Darrozes and François, 1982). The fundamental hypotheses of the mathematical description of the fluid are recalled: (a) the fluid is initially at rest, (b) the fluid is inviscid (the viscosity $\mu$ is equal to zero), (c) the fluid is incompressible (the mass density $\rho_f$ is constant), (d) the fluid flow is irrotational, thus the velocity $\vec{v}$ is governed by $\text{rot} \vec{v} = \vec{0}$.

In the following, $p$ and $\vec{f}_f$ stand, respectively, for the pressure and the external body forces for the fluid. For the applications under consideration, only gravitational forces are considered: $\vec{f}_f = -\rho_f g \vec{k}$, where $g$ is the acceleration of gravity.

With these assumptions and notations, the conservation equations of the mass and momentum are written according to the velocity of a flowing fluid element. For an incompressible fluid, the mass density is constant and the two basic hydrodynamic equations, namely continuity equations in the fluid domain $\Omega_F$ are:

$$ \mathrm{div} \, \vec{v} = 0, \qquad (1) $$

$$ \rho_f \left( \frac{\partial \vec{v}}{\partial t} + \vec{v} \cdot \nabla \vec{v} \right) = - \nabla p + \vec{f}_f. \qquad (2) $$
<figure><img src="image_2.png" /><figcaption>Fig. 1. Geometric definition of the problem.</figcaption></figure>

### 2.1.1. Potential flow theory

The preceding Eqs. (1) and (2) are rewritten by introducing an auxiliary mathematical function $\Phi = \Phi(M, t)$, the velocity potential defined at the fluid current point $M(x,y,z)$. Assuming that the fluid is initially at rest and irrotational, then it will remain irrotational, and then:

$$ \vec{v} = \nabla\Phi = \text{grad } \Phi. \qquad (3) $$

By using the mass conservation equation, the velocity potential satisfies the Laplace equation,

$$ \Delta\Phi = 0. \qquad (4) $$

By inserting Eq. (3) in the conservation equation of momentum, Eq. (2), the nonstationary version of the Bernoulli equation is established. This equation furnishes the pressure $p$ in terms of the velocity potential $\Phi$:

$$ \frac{\partial\Phi}{\partial t} + \frac{1}{2}(\nabla\Phi)^2 + \frac{p}{\rho_f} + gz = 0. \qquad (5) $$

### 2.1.2. Boundary conditions

In this subsection, the boundary conditions to be applied to the various border parts of the fluid are recalled, as follows.

#### 2.1.2.1. Free surface boundary condition. On the free-surface $\Gamma_L$, the pressure is null:

$$ p = 0. \qquad (6) $$

If gravity effects are neglected with regard to both the potential rate and quadratic terms, Eq. (5) on $\Gamma_L$ can be written as follows:

$$ \frac{\partial\Phi}{\partial t} + \frac{1}{2}(\nabla\Phi)^2 = 0. \qquad (7) $$

This condition is retained for the hydrodynamic impact formulation.

Note that if gravity effects are taken into account and assuming a small amplitude motion of the free surface and low velocities of fluid particles, the free-surface boundary conditions on $\Gamma_L$ are written as (Morand and Ohayon, 1992):

$$ \frac{\partial^2\Phi}{\partial t^2} + g\frac{\partial\Phi}{\partial z} = 0. \qquad (8) $$

The free surface is a material surface. Denoting by $z = h(x, y, t)$ its elevation relative to its initial position, then the material surface condition reduces to:

$$ \frac{d}{dt}(z - h(x, y, t)) = 0. \qquad (9) $$
2.1.2.2. Fluid-structure interface conditions. Normal velocities of fluid and solid particles are equal on the fluid-structure interface $\Gamma_{FS}$. Denoting by $\vec{V}_s$ the solid particle velocity, the boundary conditions are:

$$ \vec{v} \cdot \vec{n} = \vec{V}_s \cdot \vec{n}, \qquad (10) $$

and, replacing the fluid particle velocity by its expression in terms of the velocity potential, this condition can be expressed with the normal derivative of the velocity potential:

$$ \frac{\partial \Phi}{\partial n} = \vec{V}_s \cdot \vec{n}. \qquad (11) $$

2.1.2.3. Infinity conditions. In the case of an unbounded fluid domain, the fluid is physically at rest at a broad distance from the impact zone, and then it can be considered that

$$ |\nabla \Phi(x, y, z, t)| \to 0 \quad \text{as} \quad (x^2 + y^2 + z^2)^{1/2} \to \infty. \qquad (12) $$

As mentioned above, this condition can be replaced by a boundary condition expressed on the artificial boundary $\Gamma_{\infty}$.

2.1.3. Initial conditions

Initial conditions must be specified at impact time $t = 0$, and the far-field behaviour must be described. The potential at the initial time can be set equal to zero:

$$ \Phi(x, y, z, t = 0) = 0. \qquad (13) $$

Then, the free-surface elevation is also null initially:

$$ h(x, y, t = 0) = 0. \qquad (14) $$

2.2. Equilibrium equations of the structure problem

The deformable structure is studied under the classical small-strain isotropic elasticity assumption. Thus, no distinction between the deformed and undeformed configurations is required. Also, if $\vec{u}_s(x, y, z, t) = u_1(x, y, z, t)\vec{i} + u_2(x, y, z, t)\vec{j} + u_3(x, y, z, t)\vec{k}$ represents the displacement of a current point of the solid, then the equations of the structure motion are:

(i) The equilibrium equations in $\Omega_s$:

$$ \overrightarrow{\mathrm{div}}_{\mathrm{d}} \vec{\sigma}_{s}+\vec{f}_{s}=\rho_{s} \vec{u}_{s}, \qquad (15) $$

where $\vec{\sigma}_s$ denotes the Cauchy stress tensor, $\rho_s$ is the solid density, $\vec{u}$ stands for the acceleration vector, and $\vec{f}_s$ corresponds to body forces. Even if gravity effects do not contribute to local deformations, they influence the global motion of the body before impact.

(ii) The constitutive equation:

Hooke's law is used to express the relationship between Cauchy stresses and linearized strain tensor $\vec{\varepsilon}$ given by

$$ \varepsilon_{i,j} = \frac{1}{2}(u_{i,j} + u_{j,i}). \qquad (16) $$

Hooke's law may be expressed as

$$ \vec{\sigma}_s = \overline{\overline{D}} : \vec{\varepsilon}, \qquad (17) $$

where $\overline{\overline{D}}$ is the fourth-order elasticity tensor.

2.3. Coupled problem

Interface conditions on $\Gamma_{FS}$ are the continuity of normal velocities (Eq. (10) with $\vec{V}_s = \vec{u}_s$):

$$ \frac{\partial \Phi}{\partial n} = \vec{u}_s \cdot \vec{n} \qquad (18) $$
and the continuity of normal stress vectors, defined with the solid upward normal vector:

$$ \bar{\bar{\sigma}}_f \cdot \vec{n} = \bar{\bar{\sigma}}_s \cdot \vec{n}. \qquad (19) $$

Considering the previously described assumptions regarding the fluid, Eq. (19) becomes

$$ \bar{\bar{\sigma}}_s \cdot \vec{n} = -p\vec{n}. \qquad (20) $$

### 3. Analytical method/asymptotic method

As mentioned in Section 1, the resolution of the hydrodynamic impact problem was investigated in many works. Here, the emphasis is laid on the steps which lead to the most classical analytical approach. Analytical solutions corresponding to different domains are connected and completed to be valid in the whole fluid domain. They will be used as a reference to validate our numerical developments.

In the following, only analytical solutions related to the two-dimensional and three-dimensional rigid solids of simple geometries are given.

Counte (1989) and Wilson (1989) formalized the asymptotic analysis for the impact problem of a blunt and a rigid body. The small parameter used in the asymptotic expansion is the ratio between the immersion and the characteristic length scale of the body wetted surface. The technique of matched asymptotic expansions leads to the definition of three different zones (shown in Fig. 2), in which three asymptotic expansions are performed successively:

(i) the far-field domain, where the flow is similar to the one obtained around a flat plate of unknown width immersed in an unbounded fluid; this problem is called the outer problem;

(ii) the spray source domain near the contact line where the flow overturns to create a jet; the corresponding problem is referred to the internal problem;

(iii) the jet domain; the contribution of this jet domain to the impact problem is not significant and so it is not developed here.

#### 3.1. The outer problem

This problem reduces to a Taylor development at the first-order of the conservation equations in the neighbourhood of zero. Boundary conditions are projected onto the plane $z = 0$. The system to be solved is deduced from Eqs. (4), (7), (9) and (11):

$$ \Delta\Phi = 0 \quad \text{in } \Omega_f, \qquad (21) $$

<figure><img src="image_3.png" /><figcaption>Fig. 2. The different zones for the study of fluid-structure impact.</figcaption></figure>
$$ \frac{\partial \Phi}{\partial z} = -V \quad \text{on } \Gamma_{FS}, \qquad (22) $$

$$ \frac{\partial \Phi}{\partial z} = \frac{\partial h}{\partial t} \quad \text{on } \Gamma_L, \qquad (23) $$

$$ \frac{\partial \Phi}{\partial t} = 0 \quad \text{on } \Gamma_L, \qquad (24) $$

where $V$ is the vertical velocity of the undeformed structure. The fluid is initially at rest, so initial conditions at $t = 0$ reduce to

$$ \Phi = 0, \qquad (25) $$

$$ h = 0. \qquad (26) $$

The wetted surface is unknown in this problem, but it can be calculated as follows. First, the free-surface condition is written as

$$ \Phi = 0, \qquad (27) $$

then the pressure is obtained with the linearized Lagrange relation:

$$ p = -\rho_f \frac{\partial \Phi}{\partial t}. \qquad (28) $$

The problem so formulated is a Neumann–Dirichlet problem for the velocity potential $\Phi$. Fig. 3 shows the domain including the different equations of the problem.

<figure><img src="image_4.png" /><figcaption>Fig. 3. Outer problem for a symmetrical impact in 2D.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Fig. 4. Composite pressure field in 2D. $p_f$—inner pressure, $p_c$—outer pressure, $p_{c,\text{ composite}}$—composite pressure.</figcaption></figure>
The problem can be solved analytically by introducing a complex function (Lamb, 1932). It leads to the determination of the free-surface elevation and of the velocity potential:

$$h(x, t) = -Vt + V \int_{0}^{t} \frac{x}{x^2 - d^2(\tau)} d\tau \qquad (29)$$

and

$$\Phi = -V \sqrt{d^2(\tau) - x^2(t)}, \qquad (30)$$

where $d(t)$ represents the projection of the wetted line on the undeformed water surface $z = 0$. Thus, the pressure is given by

$$p_e(x, t) = \rho_f \frac{Vd(t)\dot{d}(t)}{\sqrt{d^2(t) - x^2}}. \qquad (31)$$

The pressure field $p_e$ for the wetted surface is shown in Fig. 4.

It should be noted that the pressure field presents a singularity at the intersection between wetted and free surfaces. At the corresponding points, the pressure tends to infinity.

## 3.2. Determination of the wetted surface

As mentioned previously, the wetted surface value $d$ is unknown in this problem. Following Wagner's (1932) proposal, $d$ is calculated by imposing that the free-surface elevation at the point $x = d$ should be equal to the body position for the solid point of the same abscissa. This can be summarized by

$$h(d(t), t) = f(d(t)) - Vt. \qquad (32)$$

Substituting Eq. (32) into Eq. (29) gives a singular integral equation. This solution was adopted by some authors such as Wagner (1932), Fabula (1957), Wanatabe (1986) and Greenhow and Yanbao (1987). Finally Wilson (1989) expressed $t$ as a function of $d$:

$$t(d) = \frac{2}{\pi} \int_{0}^{d} \frac{f(\xi)}{(d^2 - \xi^2)} d\xi, \qquad (33)$$

where $f$ is a function depending on the geometry of the body.

## 3.3. The inner problem

The inner solution of Wilson (1989) is obtained with a Taylor expansion of the variable at point $x = d$. Thus, the jet thickness is equal to

$$\delta(t) = \frac{\pi V^2 d(t)}{8 \dot{d}^2(t)}, \qquad (34)$$

its velocity is given by

$$V_{\text{jet}} = 2\dot{d}(t) \qquad (35)$$

the inner pressure can be expressed as function of a variable $u$:

$$p_i(x, y, t) = \rho_f \frac{\dot{d}(t)}{2} \left( 1 - \left( \frac{1-u}{1+u} \right)^2 \right) \qquad (36)$$

with

$$y - d(t) = \frac{\delta(t)}{\pi} \left( 2 \ln u - \frac{4}{u} - \frac{1}{u^2} + 5 \right) \quad \text{for } y < d(t), u \in [0; 1], \qquad (37)$$

$$y - d(t) = \frac{\delta(t)}{\pi} - (2 \ln u - 4u - u^2 + 5) \quad \text{for } y > d(t), u \in [1; +\infty]. \qquad (38)$$

The inner distribution $p_i$ is shown in Fig. 4.
### 3.4. Composite solution

On the one hand, the outer solution was obtained using a Taylor development in the vicinity of the origin, and so it is only valid near this point. On the other hand, the inner solution is valid at the intersection of the wetted surface and the free surface. Thus, these two solutions have to be matched. In order to obtain the solution in the complete fluid domain, using an analogous approach to that of Van Dyke (1975) or of Zhao and Faltinsen (1993) for example, the composite solution is derived.

The sum of the outer and inner solutions is corrected by subtracting their common parts. We consider that the two solutions are equivalent far from their zone of validity, and so their limits at $-\infty$ must be equal.

Thus, as $u$ tends toward zero

$$ (p_i)_e \sim 2\rho_f d^2 u \qquad (39) $$

for

$$ y - d(t) \sim -\frac{\delta(t)}{\pi} u^{-2}. \qquad (40) $$

Replacing $u$ and $\delta(t)$, the common pressure reduces to

$$ (p_i)_e = \rho_f \frac{Vd(t)\dot{d}(t)}{\sqrt{2d(t)(d(t) - y)}}. \qquad (41) $$

The final pressure solution $p_e$ can be expressed by the sum of the outer and inner pressures and the subtraction of their common parts:

$$ p_e = p_e + p_i - (p_i)_e. \qquad (42) $$

Fig. 4 shows how the pressure field is deduced from the inner and outer pressures.

Considering that the velocity does not remain constant over time, this should be considered in the wetted surface calculation (Wilson, 1989). The outer pressure is given by

$$ p_e(x, t) = \rho_f \frac{V(t)d(t)\dot{d}(t)}{\sqrt{d^2(t) - x^2}} + \rho_f \dot{V}(t)\sqrt{(d^2(t) - x^2)}. \qquad (43) $$

This matching technique does not modify either the jet thickness or the inner pressure expression.

## 4. Numerical approach

The finite element method is a popular technique that leads to computer solution of complex problems in engineering. Here, it is used to solve the full set of fluid and structural equations simultaneously. In the following, we present the discrete set of equations.

### 4.1. General discretization of the coupled problem

Utilizing the notations of Zienkiewicz and Taylor (1991), the finite element equations of the fluid structure problem are briefly recalled.

#### 4.1.1. Variational formulation of the fluid problem

The bounded fluid domain $\Omega_F$ is divided into elementary sub-domains $\Omega_e$ in which the potential velocity $\Phi$ is interpolated according to nodal potential variables $\{\Phi\}_e$:

$$ \Phi = \langle N_f \rangle \{\Phi\}_e = \langle N_1, \dots, N_i, \dots \rangle \{\Phi\}_e, \qquad (44) $$

where the interpolation functions $\langle N_f \rangle$ depend on the type of element, and $\{\Phi\}_e$ denotes the nodal potential vector of the finite element (e). The potential $\Phi$ is determined by using the weak formulation of the fluid problem.
For this purpose, let $\Psi$ denote an arbitrary function; then, Eq. (45) is the basis of the finite element discretization process:

$$ \int_{\Omega_f} \nabla \Psi \cdot \nabla \Phi \, dV = \int_{\partial \Omega_2} \Psi \left( \frac{\partial \Phi}{\partial n} \right)_d \, d\Gamma. \qquad (45) $$

The boundary $\partial\Omega_f$ consisting of two complementary parts, part $\Gamma_\infty$ on which the potential function $\Phi$ is given (Dirichlet boundary condition) and part $\partial\Omega_2 = \Gamma_L \cup \Gamma_{SF}$ on which the normal derivative of the potential $\partial\Phi/\partial n$ is enforced (Neumann boundary condition).

On the fluid structure interface $\Gamma_{SF}$, Eq. (11) becomes

$$ \nabla \Phi \cdot \vec{n} = \left( \frac{\partial \Phi}{\partial n} \right)_d = -\vec{V}_S \cdot \vec{n}, \qquad (46) $$

where $\vec{n}$ is the normal direction to $\Gamma_{SF}$, pointing upward of the structure.

The Galerkin method consists of developing the function $\Psi$ on the same shape functions as those used for the interpolation of velocity potential $\Phi$. Hence, Eq. (45) leads to the general matrix system according to

$$ [M_f]\{\ddot{\Phi}\} + [H_f]\{\Phi\} = \{G\}, \qquad (47) $$

where matrices are determined by connecting the elemental matrices calculated on the subdomains $\Omega_e$ and on the elemental surface $\Gamma_e$, the general terms of these elemental matrices being given by

$$ \begin{aligned} (h_{ij})_e &= \int_{\Omega_e} \nabla N_i \cdot \nabla N_j \, dV, \\ (g_i)_e &= \int_{\Gamma_e} N_i(-\vec{V}_s \cdot \vec{n}) \, d\Gamma, \quad \text{with } \Gamma_e \in \Gamma_{SF}. \end{aligned} \qquad (48) $$

For applications involving gravity waves, the boundary condition on the free surface (Eq. (8)) reduces to

$$ \frac{\partial \Phi}{\partial n} = - \frac{1}{g} \frac{\partial^2 \Phi}{\partial t^2} $$

and the matrix $[M_f]$ is determined by connecting elemental matrices calculated on the elemental free surface $\Gamma_{Le}$:

$$ (m_{ij})_e = \int_{\Gamma_{Le}} N_i \frac{1}{g} N_j \, d\Gamma. \qquad (49) $$

For the impact problem, the free-surface condition deduced from Eq. (7) leads to $[M_f] = [0]$, and the condition $p = 0$ permits to calculate the free surface elevation.

### 4.1.2. Variational formulation of the structural problem

Similarly, in order to obtain an approximation of the displacement field, the body $\Omega_s$ is divided into finite elements. The displacement is interpolated in each element:

$$ \{u\} = [N_s]\{a\}_e, \qquad (50) $$

where $\{u\}$ represents the components of the displacement vector $\vec{u}$ in the body $\Omega_s$, $\{a\}_e$ is the nodal displacement vector of the nodes of the finite element $\Omega_e$, and $[N_s]$ is the matrix of interpolation functions.

Then Eq. (16), the relation between nodal strain and displacements, gives

$$ \{e\} = [B]\{a\}_e, \qquad (51) $$

where $[B]_e$ is the strain-displacement matrix. The discrete form of Hooke's law is easily established by replacing $\{e\}$ in Eq. (17):

$$ \{\sigma\} = [D][B]\{a\}_e, \qquad (52) $$

where $[D]$ is the usual set of elastic constants. Differentiating the strain energy expression on each finite element domain leads to the stiffness matrix:

$$ [k]_e = \int_{\Omega_e} [B]^T [D] [B] \, dV. \qquad (53) $$
Similarly, the elemental mass matrix is obtained from the differentiation of the kinetic energy in $\Omega_e$:

$$[m]_e = \int_{\Omega_e} [N_s]^T \rho_s [N_s] dV. \qquad (54)$$

The virtual work of external forces leads to the definition of nodal forces. For the fluid-structure impact applications under consideration, there are external distributed loads. The nodal vector forces is

$$\{\phi_s\}_e = \int_{\Omega_e} [N_s]^T \{q_s\} dV. \qquad (55)$$

The matrix form of the whole structure is obtained by the complete set of each finite element contribution, and by taking the variations of the resulting quantities:

$$[M_s]\{\ddot{a}\} + [K_s]\{a\} = \{F_{\text{ext}}\} + \{\phi_s\}_{\text{ass}}, \qquad (56)$$

where $\{F_{\text{ext}}\}$ is the external nodal vector forces directly applied to the nodes of the structure.

### 4.1.3. General matrix system of the coupled problem

The distributed loads $\{q_s\}$ are the dynamic pressures ($-p\bar{n}$) on the $\Gamma_{SF}$ interface due to the fluid loading on the structure.

The associated nodal forces ($\{\phi_h\}_e = \{\phi_s\}_e$) expression is

$$\{\phi_h\}_e = \int_{\Gamma_{SFe}} [N_s]^T \{n\} (-p) dS. \qquad (57)$$

In the same way, on each fluid finite element in contact with the interface $\Gamma_{SF}$, the interpolation of the potential velocity according to the elemental nodal values Eq. (44) leads to the determination of the dynamic pressure in any point of a fluid finite element. The pressure is then defined by

$$p = -\rho_f \langle N_f \rangle \langle \dot{\Phi} \rangle_e \qquad (58)$$

and on the elementary fluid-structure interface $\Gamma_{SFe}$, nodal forces corresponding to the hydrodynamic pressure are represented by the vector

$$\{\phi_h\}_e = \rho_f \int_{\Gamma_{SFe}} [N_s]^T \{n\} \langle N_f \rangle dS \dot{\{\Phi\}}_e, \qquad (59)$$

which can be written as

$$\{\phi_h\}_e = \rho_f [s f]_e \{\dot{\Phi}\}_e. \qquad (60)$$

The $[s f]_e$ matrix represents the elemental fluid-structure coupling matrix associated to the finite element approach:

$$[s f]_e = \int_{\Gamma_{SFe}} [N_s]^T \{n\} \langle N_f \rangle dS. \qquad (61)$$

Proceeding in a similar way for the fluid discretization analysis, one can express normal velocity in an unspecified point of the interface $\Gamma_{SFe}$ belonging to a finite element (e) by

$$V_n = \vec{V}_s \cdot \bar{n} = \{n\}^T [N_s]_{\Gamma_{SFe}} \{\dot{a}\}_e, \qquad (62)$$

thus, the $(g_i)_e$ terms calculated on $\Gamma_{SFe} \in \Gamma_{SF}$ are written as

$$\{g\}_e = - \int_{\Gamma_{SFe}} \langle N_f \rangle^T \{n\}^T [N_s] dS \{\dot{a}\}_e, \qquad (63)$$

or

$$\{g\}_e = -[s f]_e \{\dot{a}\}_e = -[s f]_e^T \{\dot{a}\}_e. \qquad (64)$$

If the fluid and the structure interface discretization do not coincide, the calculation of the preceding quantities is performed considering integrals on the elementary fluid structure interfaces. With elementary contributions coming from the weak forms, we obtain the following two linear sets of coupled systems equations:

$$\begin{aligned} [M_s]\{\ddot{a}\} + [K_s]\{a\} &= \rho_f [SF]\{\dot{\Phi}\}, \\ [M_f]\{\dot{\Phi}\} + [H_f]\{\Phi\} &= -[FS]\{\dot{a}\} \end{aligned} \qquad (65)$$
with $[FS]$ or $[SF] = [FS]^T$ assembled fluid-structure matrices. In the general case this system is solved by gathering the nodal unknown factors $\{a\}$ and $\{\Phi\}$ (Morand and Ohayon, 1992) in the same nodal unknown vector $\{W\}({\{W\}}^T = \{\{a\}^T, \{\Phi\}^T\})$ such that:

$$[M]\{\ddot{W}\} + [C]\{\dot{W}\} + [H]\{W\} = \{Q\} \qquad (66)$$

with

$$[M] = \begin{bmatrix} [M_s] \\ [M_f] \end{bmatrix}, \quad [C] = \begin{bmatrix} -\rho_f [SF] \end{bmatrix}, \quad [H] = \begin{bmatrix} [K_s] \\ [H_f] \end{bmatrix} \qquad (67)$$

and here for the impact problems under consideration, $\{Q\} = \{0\}$. In the following, the validation part, we shall present several examples solved by the coupled fluid-structure finite-element approach presented above.

## 4.2. Determination of the wetted surface

As noted above, one of the major difficulties of the problem is due to the determination of the position of the moving boundaries. Hence, the first step of the resolution consists in the determination of the fluid-structure interface motion during the impact.

To evaluate the right-hand side of Eq. (48), it is necessary to know the wetted surface of the body, or, similarly, the position of the contact line between the body and the fluid. Due to the deformation of the free surface, the wetted surface is an unknown variable of the problem. An additional equation is therefore needed to solve the problem. Practically, the contact line is determined by forcing the solution to satisfy the principle of volume conservation (Wilson, 1989; Fontaine and Cointe, 1992) or, similarly, by enforcing the existence of an intersection point between the free surface elevation and the body, as intuitively done by Wagner (1932). The volume conservation principle should be automatically satisfied since the Laplace equation was used. Nevertheless, the problem is singular at the initial time when the contact line reduced to a single point, and thus an extra closure equation is needed.

As suggested by Korobkin (1982), it is easier to introduce the displacement potential given by

$$\Psi(x, y, z, t) = \int_{0}^{t} \Phi(x, y, z, s) \, ds, \qquad (68)$$

than of the velocity potential. The boundary value problem satisfied by $\Psi$ is similar to the one satisfied by $\Phi$, except for the body boundary conditions. The main advantage of introducing this transformation is that the value of $\Psi$ in $z = 0$ is defined in terms of the free-surface elevation motion. The determination of $\Psi$ does not require the knowledge of the temporal evolution of the free surface elevation. From a numerical point of view, the discretization scheme used for the displacement potential is the same at the one applied for the velocity potential Eq. (45) and leads to

$$[H_f]\{\Psi\} = \{G_\psi\}, \qquad (69)$$

where the boundary condition on the interface $\Gamma_{FS}$ can be expressed as

$$\{g_\psi\}^e = \int_{\Gamma_{Sfe}} \{N_f\}(-\vec{u}_s \cdot \vec{n}) \, ds. \qquad (70)$$

The wetted surface is determined using an iterative procedure. First, the wetted surface is expressed as a function of a previous wetted surface. Then, the elevation of the free surface is obtained and, a new wetted surface defined by the intersection between the free surface and the structure is calculated. The resolution is performed until the free surface elevation is equal to the position of the body at the boundary $\Gamma(t)$:

$$h(x, y, t) = \vec{u}_s \cdot \vec{z}. \qquad (71)$$

Once convergence is ensured, the position of the contact line is known. The linear system for the velocity potential is then solved. It is noted that the matrix $[H_f]$ only depends on the previously constructed geometry and thus does not have to be computed once again to solve the velocity potential problem. From one time step to the other, the boundary mesh does not change, and therefore no substantive derivative has to be introduced. As a result, a simple first-order backward finite-difference scheme provides an accurate estimation of the pressure, even close to singularity, as will be demonstrated in the next section.

Fig. 5 shows the fluid domain finite element discretization for an example of a wedge impacting water. The mesh must be refined near the body and near the free surface. A numerical convergence study (Donguy, 2002) showed that the size
<figure><img src="image_6.png" /><figcaption>Fig. 5. Wedge impact: example of finite element discretization in the fluid domain.</figcaption></figure>

of the fluid-structure interface elements has to be equal to approximately 1/6 of the increase in wet length between two time increments. The wetted interface is composed of approximately 200 nodes. The number of nodes on the free surface does not affect the final results.

## 4.3. Validation

### 4.3.1. Validation of the numerical procedure: rigid structure

In this section, calculations are carried out using the finite element code Cast3M developed by the French Atomic Agency (Verpeaux et al., 1988).

#### 4.3.1.1. Two-dimensional case.
In the following, we will consider the example of a two-dimensional parabolic hull defined by $z = ay^2$, which penetrates a free surface initially at rest with a constant vertical velocity $V$. The numerical method previously described is applied to this case. Comparisons between numerical and analytical results are presented in Figs. 6(a) (free surface elevation) and (b) (pressure distribution). Triangular and quadrilateral quadratic interpolation elements are used for the numerical computation.

Both meshes usually coincide with the interface, but in the more general case of distinct meshes, it is necessary to project physical fields. Several authors have proposed ways addressing this issue. For example,

(i) Nayroles et al. (1992) have proposed to use diffuse approximation, one of the numerous meshless methods;

(ii) Moyroud (1998) and Moyroud et al. (2000) have proposed a grid method for both the projection of fluid and structure field;

(iii) finally, Tanizawa (1998) proposed a specific interpolation method for the interface problem.

For impact problems, as presented here, no specific procedures are developed. Elementary coupling matrices (Eqs. (61) and (63)) are determined as presented in Section 4.1.3, by using both fluid and structure interpolation fields.
<figure><img src="image_7.png" /><figcaption>Fig. 6. Parabolic hull impact. (a) Comparison between numerical and analytical solutions for the free surface elevation. (b) Comparison between numerical and analytical solutions for the outer pressure distribution. The error distribution is plotted on the right.</figcaption></figure>

The analytical expressions for the wetted width, the free surface elevation and the outer expansion of the pressure are given by Cointe (1989):

$$d^2(t) = 2 \frac{Vt}{a}, \qquad (72)$$

$$h(y, t) = -Vt - \left[ \sqrt{ay}(ay^2 - 2Vt)^{1/2} - ay^2 \right], \qquad (73)$$

$$p_e(y, t) = \rho_f V^2 (a(2Vt - ay^2))^{1/2}. \qquad (74)$$

Although the solution is singular at the intersection, the numerical scheme gives a relatively good estimation of the pressure distribution in this region. More precisely, the increase in the error near the singularity is not consequential, since the numerical solution agrees with the nonlinear spray root solution.

#### 4.3.1.2. Axisymmetric rigid impact case.
This problem solved by Donguy et al. (2001b) relates to the impact case of a rigid cone ($\beta = 10°$ angle) with an impact velocity of 6 m/s. Fig. 7 shows the mesh adopted for the fluid. Fig. 8(a) presents a comparison of numerical and experimental solutions for the free surface elevation at time $t = 10$ ms and Fig. 8(b) compares numerical and experimental outer pressures and shows the error between these two solutions. Hexahedral elements are used for the numerical simulation. Good agreement between numerical and analytical results is revealed. All parameters being comparable to the previous two-dimensional computation, the relative error remains of the same order of magnitude, although the mesh is three dimensional this time. In Fig. 9, the pressure distribution at times $t_1, t_2, t_3, t_4$, for a two-dimensional wedge and a cone immersed equivalently are roughly compared to estimate three-dimensional effects. The figure also shows that the two-dimensional computations is conservative.

To conclude this section, it can be said that analytical and numerical results are very close to each other for both 2-D and 3-D cases. Nevertheless, the evaluation of the free-surface elevation remains sensitive to the boundary between the structure and the fluid. To represent the singular behaviour of the solution in the vicinity of this singularity part of the boundary, the finite-element mesh has to be refined in this zone. The regularity of the mesh is also an important parameter, which controls the accuracy of the numerical solution.

The outer pressure is corrected by the analytical solution of the inner pressure and by considering that the outer pressure is accurate close to the origin and the inner solution is accurate close to the intersection of the wetted surface and the free surface.
<figure><img src="image_8.png" /><figcaption>Fig. 7. Three-dimensional mesh used in finite element method around a cone.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Fig. 8. Cone shaped hull impact. (a) Comparison between numerical and analytical solutions for the free surface elevation. (b) Comparison between numerical and analytical solutions for the outer pressure distribution.</figcaption></figure>

4.3.1.3. A three-dimensional case for a paraboloid ellipse. The calculation (Tourbier et al., 2002b) is performed for a paraboloid elliptic case for which the analytical solution is known (Scolan and Korobkin, 2001). Fig. 10(a) gives a glimpse of the paraboloid ellipse, and Fig. 10(b) the motions of the free and wetted surfaces. Fig 11(a) compares numerical and analytical results for the elevation of the free surface and Fig. 11(b) the outer pressure distribution. The results are close to each other.
<figure><img src="image_10.png" /><figcaption>Fig. 9. Comparison between the pressure distributions along a cone and a wedge penetrating the free surface with same velocity. At each time, the immersion is the same for the two bodies.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Fig. 10. (a) Representation of a paraboloid elliptic. (b) Mesh of the wetted and free surface.</figcaption></figure>

### 4.3.2. Validation of the numerical procedure: deformable body

In this section, the previous analysis is extended to solve the case of deformable body. First, the equations governing the structural deformations are recalled. A coupled method to solve the fluid-structure interaction problem is then proposed.

#### 4.3.2.1. Validation of the coupling matrix [FS] in 2-D.
The numerical resolution of the general matrix system Eq. (66) taking into account the gravity effect is performed using the commercial Finite Element software Cast3M (Verpeaux et al., 1988). An external procedure is implemented for the evaluation of the coupling matrix [FS]. The correct implementation of the approach is checked for the test-case problem of sloshing in a rectangular tank including an
<figure><img src="image_12.png" /><figcaption>Fig. 11. (a) Comparison between analytical and numerical results for the elevation of the free surface. (b) Comparison between analytical and numerical results for the outer pressure.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Fig. 12. Elastic beam in a heavy fluid problem.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>Fig. 13. Elastic beam in a heavy fluid problem: first eigenmodes shapes for the coupled problem.</figcaption></figure>

elastic plate in its center, as shown in Fig. 12 for a two-dimensional case. This problem was experimentally studied by Chai (1996). The problem dimensions, material and fluid properties are: $L = 0.266 \text{ m}$, $e = 6 \times 10^{-4} \text{ m}$, $l = 0.231 \text{ m}$, $h_f = 0.210 \text{ m}$, $E = 2.1 \times 10^{11} \text{ Pa}$, $\rho_s = 7800 \text{ kg m}^{-3}$, $\rho_f = 1000 \text{ kg m}^{-3}$.
<table><caption>Table 1
Eigenfrequencies of immersed beam (in Hertz)</caption><thead><tr><td>Frequency (Hz)</td><td>Reference case (Peseux et al., 1999)</td><td>C1</td></tr></thead><tbody><tr><td>f<sub>1</sub></td><td>1.49</td><td>1.60</td></tr><tr><td>f<sub>2</sub></td><td>1.86</td><td>1.84</td></tr><tr><td>f<sub>3</sub></td><td>2.39</td><td>2.46</td></tr><tr><td>f<sub>4</sub></td><td>2.66</td><td>2.63</td></tr><tr><td>f<sub>5</sub></td><td>2.97</td><td>3.06</td></tr><tr><td>f<sub>6</sub></td><td>3.33</td><td>3.29</td></tr><tr><td>f<sub>7</sub></td><td>3.46</td><td>3.51</td></tr></tbody></table>

<figure><img src="image_15.png" /><figcaption>Fig. 14. (a) Elastic plate in a heavy fluid problem. (b) Mesh of the elastic plate problem.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>Fig. 15. Elastic plate in a heavy fluid problem: first eigenmodes shapes for the coupled problem.</figcaption></figure>

Mechanical loading conditions consist in a harmonic force applied to the beam extremity. Frequency scanning is performed to evaluate the eigenfrequencies, which correspond to peaks of beam displacement. Both beam modal shape and free surface elevation for the first four modes are shown in Fig. 13. Eigenfrequencies are presented in Table 1, column C1. They are compared with the results of Peseux et al. (1999). It can be seen that the frequencies obtained by
<table><caption>Table 2
Eigenfrequencies of immersed plate (in Hertz)</caption><thead><tr><th>Frequency (Hz)</th><th>Reference case (Peseux et al., 1999)</th><th>C1</th><th>Experiment</th></tr></thead><tbody><tr><td>f<sub>1</sub></td><td>171.1</td><td>170.1</td><td>170.9</td></tr><tr><td>f<sub>2</sub></td><td>266.4</td><td>258.9</td><td>265.3</td></tr><tr><td>f<sub>3</sub></td><td>500.1</td><td>525.6</td><td>484.4</td></tr><tr><td>f<sub>4</sub></td><td>610.5</td><td>617.9</td><td>597.5</td></tr><tr><td>f<sub>6</sub></td><td>1037</td><td>1083</td><td>941</td></tr><tr><td>f<sub>7</sub></td><td>1052</td><td>1089</td><td>998</td></tr><tr><td>f<sub>8</sub></td><td>1341</td><td>1265</td><td>1228</td></tr></tbody></table>

the different methods are relatively close to each other, especially the results in the two last columns. The computations by Peseux et al. were performed using another FE code, which may explain the small differences between results. The agreement between the three methods is very satisfactory.

### 4.3.2.2. Validation of the coupling matrix in 3-D.
To validate the coupling matrix in a three-dimensional case, we consider the case of sloshing of an elastic plate in a rectangular tank, as shown in Fig. 14(a). The dimensions, material and fluid properties are: $L = 0.40 \text{ m}$, $e = 5.7 \times 10^{-3} \text{ m}$, $h = 0.10 \text{ m}$, $E = 2.1 \times 10^{11} \text{ Pa}$, $\rho_s = 7800 \text{ kg m}^{-3}$, $\rho_f = 1000 \text{ kg m}^{-3}$.

Fig. 14(b) shows the finite elements mesh used, with hexahedral elements for the fluid and quadrilateral elements for the structure. As above, the mechanical loading conditions consist of a harmonic force applied to the plate extremity. A frequency scan is performed to evaluate the eigenfrequencies which correspond to the peaks of the plate displacement. The free-surface elevation and modal shapes of the plate for the first four modes of the coupled problem are shown in Fig. 15. The numerical values of the eigenfrequencies for the two dimensional case are reported in Table 2 column C1. They are compared with the results of Peseux (1989). Finally, experimental results for plate vibrations (Peseux, 1989) are given in column “Experiment”. Comparison of frequencies shows good agreement, which lends confidence to the external numerical procedure developed in the present study.

### 4.3.2.3. Temporal evolution.
The fluid-structure interaction problem, Eq. (66), is then integrated starting from the initial conditions:

$$W(t=0) = W_0, \tag{69'}$$

$$\dot{W}(t=0) = \dot{W}_0. \tag{70'}$$

A numerical temporal implicit scheme is used. Nevertheless, when studying rigid body slamming, it appears that time step $\delta t$ for the fluid problem must be less than one fifth of the time step used to calculate accurately the pressure distribution. Choosing the same time step for the coupled problem simulation is too time consuming.

Consequently, the fluid-structure interaction problem is solved using an iterative procedure (Tombier et al., 2002a,b). The fluid problem is first calculated with a step time $\delta t$. After determining the force due to pressure, the structural resolution can be done using a Newmark temporal integration scheme with the step time of $\Delta t$ with $\Delta t = 50 \delta t$.

Fig. 16 presents the time scheme used for this simulation and Fig. 17 shows the structure deformation obtained during simulation.

## 5. Experiments

In this part, experimental tests due to Donguy (2002) are presented. These tests were conducted to provide data to validate the previous numerical approach. We first describe the experimental procedure. Then, results on rigid and deformable structures are presented, and finally comparisons with the numerical results are carried out.
<figure><img src="image_17.png" /><figcaption>Fig. 16. Iterative schema used to solve the coupled problem.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>Fig. 17. Shape deformation of the structure with the coupled resolution.</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>Fig. 18. Impact tower: experimental set-up.</figcaption></figure>
<figure><img src="image_20.png" /><figcaption>Fig. 19. Detailed views of the rigid cone shaped models.</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>Fig. 20. Detailed views of the deformable cone shaped models.</figcaption></figure>

## 5.1. *Experimental set-up*

The investigation tests consist of free-fall drop tests of cone-shaped models (Donguy et al., 2001a). Drop tests with impact velocities varying from 2 to 8 m/s are carried out from a free falling rig system of 3.5 m height (see Fig. 18). The cone falls in a water tank 1 m deep and 1.20 m in diameter. Foam lines the tank bottom to cushion the cone's impact.

Cones-shaped models are rigidly attached to a cylindrical support of 322 mm diameter (see Fig. 19). Three different cone-shaped models with deadrise angles of 6°, 10°, and 14° are successively tested. To study the deformations, both rigid bodies and elastic structures are tested. The effective thickness of the steel models is sufficient (from 25 to 50 mm), to ensure that they do not deform during impact. For the elastic structures, three thicknesses of 0.5, 1 and 1.5 mm (see Fig. 20) are considered.

A relatively large experimental facility allows a typical size of 0.32 m for the tested models and therefore reduces the possible influence of surface tension effects. The impact velocity is set to increase by realistic values from approximately 2–6 m/s.

Drop velocity and pressure measurements are performed during impact. Two quartz ICP compensated pressure sensors are placed at 40 and 90 mm from the cone symmetry axis (Fig. 20). These sensors are well suited for impact measurements: their sampling frequency is up to 400 kHz, and the measurement range is 0–69 bars. They also allow for the use of relatively long wires (20 m) without altering the electric signal, so that the data-processing system does not need to be too close to the water. The drop velocity is measured by means of a cable sensor connected to the cylindrical support. Finally, pressure and velocity signals are plotted and recorded by a numerical memory oscilloscope. The cones
<table><caption>Table 3
Comparison between the peak pressure values for impact velocity equal to 5.2 m/s for rigid cone for three experimental tests</caption><thead><tr><th>Cone</th><th>Test</th><th>V<sub>impact</sub> (m/s)</th><th>Pressure P<sub>1</sub> (bar)</th><th>Pressure P<sub>2</sub> (bar)</th></tr></thead><tbody><tr><td rowspan="3">6°</td><td>E<sub>1</sub></td><td>5.3</td><td>10.8</td><td>11.2</td></tr><tr><td>E<sub>2</sub></td><td>5.4</td><td>10.9</td><td>9.8</td></tr><tr><td>E<sub>3</sub></td><td>5.3</td><td>10.4</td><td>12.2</td></tr><tr><td rowspan="3">10°</td><td>E<sub>1</sub></td><td>5.1</td><td>4.1</td><td>5.6</td></tr><tr><td>E<sub>2</sub></td><td>5.3</td><td>3.9</td><td>5.4</td></tr><tr><td>E<sub>3</sub></td><td>5.2</td><td>3.8</td><td>5.4</td></tr><tr><td rowspan="3">14°</td><td>E<sub>1</sub></td><td>5.1</td><td>2.6</td><td>2.7</td></tr><tr><td>E<sub>2</sub></td><td>5.2</td><td>2.6</td><td>2.7</td></tr><tr><td>E<sub>3</sub></td><td>5.1</td><td>2.6</td><td>2.6</td></tr></tbody></table>

<table><caption>Table 4
Comparison between the peak pressure values for impact velocity equal to 2.6 m/s for elastic cone for three experimental tests</caption><thead><tr><th rowspan="2">Cone</th><th rowspan="2">Test</th><th colspan="2">Pressure P<sub>1</sub> (bar)</th><th colspan="2">Pressure P<sub>2</sub> (bar)</th></tr><tr><th>First peak</th><th>Second peak</th><th>First peak</th><th>Second peak</th></tr></thead><tbody><tr><td rowspan="3">6° (t = 1.5 mm)</td><td>E<sub>1</sub></td><td>1.3</td><td>0.8</td><td>0.2</td><td>1.5</td></tr><tr><td>E<sub>2</sub></td><td>1.5</td><td>0.5</td><td>0.1</td><td>0.9</td></tr><tr><td>E<sub>3</sub></td><td>1.5</td><td>0.5</td><td>0.3</td><td>0.6</td></tr><tr><td rowspan="3">10° (t = 1.0 mm)</td><td>E<sub>1</sub></td><td>0.7</td><td>1.1</td><td>0.7</td><td>1.1</td></tr><tr><td>E<sub>2</sub></td><td>1.1</td><td>1.8</td><td>1.0</td><td>2.2</td></tr><tr><td>E<sub>3</sub></td><td>1.2</td><td>1.8</td><td>1.0</td><td>2.4</td></tr><tr><td rowspan="3">14° (t = 0.5 mm)</td><td>E<sub>1</sub></td><td>0.6</td><td>1.0</td><td>0.7</td><td>1.0</td></tr><tr><td>E<sub>2</sub></td><td>0.6</td><td>0.9</td><td>0.9</td><td>1.0</td></tr><tr><td>E<sub>3</sub></td><td>0.6</td><td>0.8</td><td>1.0</td><td>1.1</td></tr></tbody></table>

are dropped into calm water. At least three drops were performed for each test condition to ensure repeatability of the measurements. Repeatability of tests is illustrated by the experimental values presented in Table 3 which correspond to the drop tests of rigid cones with an average impact velocity of 5.6 m/s. Table 4 gives results obtained for the deformable cones.

The results presented below are representative of these drop-test series.

## 5.2. Signal on the rigid body

Fig. 21 presents the evolution of pressure recorded by the two sensors $p_1$, $p_2$ and the velocity during the impact of a 10° deadrise angle cone with an impact velocity of 2.5 m/s. For other velocities and cone characteristics, the resulting curves are similar. Fig. 22 shows the position of the free surface during impact and renders the pressure evolution shown in Fig. 21 easier to understand. Firstly, the global velocity impact is controlled, then it remains constant during water entry. The cable sensor used to measure velocity causes some small signal disturbance. A similar evolution of the measured pressure is observed with both sensors. First, we observe a slight depression, just before a pressure peak. This pressure depression was observed by some other authors and is due to the jet flow. Next, pressure peaks are reached when the sensors are in the jet formation region and more precisely at the stagnation point. The increase of the pressure level depends on the velocity at impact and the deadrise angle of the cone. Its values range from 10 µs for the fastest and 200 µs for the slowest. After the peak pressure, we note a slow pressure decrease while the cone is progressing into the water. This is in accordance with the fact that the sensors are in the water. Finally, a sudden decrease in pressure occurs
<figure><img src="image_22.png" /><figcaption>Fig. 21. Pressure and velocity evolution during impact of the 10° rigid cone, $V_{impact} = 2.5$ m/s.</figcaption></figure>

<figure><img src="image_23.png" /><figcaption>Fig. 22. Correspondence between pressure history and cone location.</figcaption></figure>

<figure><img src="image_24.png" /><figcaption>Fig. 23. Correspondence between the sudden pressure decrease and the cone position in the water tank.</figcaption></figure>
<figure><img src="image_25.png" /><figcaption>Fig. 24. Pressure evolution during impact of the 6° deformable cone: (a) $e = 0.5$ mm, (b) $e = 1$ mm, (c) $e = 1.5$ mm, $V_{\text{impact}} = 2.5$ m/s.</figcaption></figure>

<figure><img src="image_26.png" /><figcaption>Fig. 25. Pressure evolution during impact of the 10° deformable cone: (a) $e = 0.5$ mm, (b) $e = 1$ mm, (c) $e = 1.5$ mm, $V_{\text{impact}} = 2.5$ m/s.</figcaption></figure>

<figure><img src="image_27.png" /><figcaption>Fig. 26. Pressure evolution during impact of the 14° deformable cone: (a) $e = 0.5$ mm, (b) $e = 1$ mm, (c) $e = 1.5$ mm, $V_{\text{impact}} = 2.5$ m/s.</figcaption></figure>

simultaneously in both sensors. At this instant, the jet zone is on the cylindrical support and the cone is totally immersed. The jet separates when airborne, so that the pressure on the cone suddenly decreases. This phenomenon is shown in Fig. 23.

### 5.3. Signals on the elastic structure

Figs. 24–26 show the evolution of pressure with time. The measured signals are representative of whole measurements made for all cone-shaped models. The tests performed on the elastic structure show more disparity than the rigid ones.
<table><caption>Table 5
Comparison between predicted and experimental peak pressure values for impact velocity equal to 5.2 m/s for rigid cone</caption><thead><tr><th rowspan="2">V<sub>impact</sub> = 5.2 m/s</th><th colspan="2">Experiment (bar)</th><th rowspan="2">Numerical (bar)</th></tr><tr><th>P<sub>1</sub></th><th>P<sub>2</sub></th></tr></thead><tbody><tr><td>Cone 6°</td><td>10.3</td><td>10.7</td><td>19.8</td></tr><tr><td>Cone 10°</td><td>3.9</td><td>5.5</td><td>7</td></tr><tr><td>Cone 14°</td><td>2.6</td><td>2.7</td><td>3.5</td></tr></tbody></table>

<table><caption>Table 6
Comparison between predicted and experimental peak pressure values for impact velocity equal to 2.5 m/s for deformable cones with a thickness equal to 1 mm</caption><thead><tr><th rowspan="2"></th><th colspan="2">Sensor P1</th><th colspan="2">Sensor P2</th></tr><tr><th>Experiment (bar)</th><th>Numerical (bar)</th><th>Experiment (bar)</th><th>Numerical (bar)</th></tr></thead><tbody><tr><td>Cone 6°</td><td>1.4</td><td>2.2</td><td>1.3</td><td>2.8</td></tr><tr><td>Cone 10°</td><td>1.2</td><td>1.3</td><td>1</td><td>1.5</td></tr><tr><td>Cone 14°</td><td>0.6</td><td>0.6</td><td>0.8</td><td>0.7</td></tr></tbody></table>

In particular, some barely perceptible phenomena in the case of rigid bodies become noticeable and highlight the influence of structural deformation.

New phenomena occur on all elastic body shapes; this corresponds to the occurrence of a secondary pressure peak. Figs. 24–26 show the monitoring of the pressure for the 6°, 10° and 14° deadrise angles for specimens with several thicknesses. For each deadrise, three thicknesses are considered ((a) $e = 0.5$ mm, (b) $e = 1$ mm, (c) $e = 1.5$ mm). These secondary peaks of pressure or depression may have amplitude greater than the amplitude of the first pressure peaks. These secondary peaks appear when the cylindrical support reaches the free surface. The interpretation of the experimental results is difficult. New tests with larger dimensions will permit to reach a conclusion this secondary peak phenomena.

## 5.4. Comparisons between theory and experiments

Tables 5 and 6 show the experimental and numerical maximum values of the pressure on sensors $p_1$ and $p_2$ for the three deadrise angles ($\beta = 6°$, $10°$ and $14°$). Numerical and experimental results for the pressure levels are found to be in reasonable agreement although the pressure is always slightly overpredicted, indicating that the numerical model is conservative.

For the first sensor, the experimental value is far overpredicted for all deadrise angles. The difference is less pronounced for experimental values from the second sensor. A certain reduction in the pressure can be explained by the physical process of air entrapment. This phenomenon was also observed by Chuang and Milne (1971) on cones with deadrise angles smaller than 1°, and by Hagiwara and Yuhara (1974) for wedges with deadrise angles lower than 3°.

Examples of the temporal pressure evolution on the two probes are presented in Fig. 27 with the different deadrise angles ((a) $\beta = 6°$, (b) $\beta = 10°$ (c) $\beta = 14°$) for rigid bodies on the left and for deformable bodies with a thickness of 1 mm on the right. Globally, the experimental distribution corresponds well to numerical results. Both signals exhibit a pressure peak traveling along the body, followed by a constant value.

The location of the peak is clearly quite well described by this simplified theory. The wetted length seems to be correctly estimated. This local quantity is derived from the global volume conservation property, as mentioned previously. The pressure rise immediately after impact is also well-reproduced. The simplified model is indeed asymptotically valid for small times. The level of the peak, i.e., the maximum pressure, is reasonably estimated for the biggest deadrise angle, but the error increases with a decrease of deadrise angle. According to the asymptotic solution,
<figure><img src="image_28.png" /><figcaption>Fig. 27. Pressure histories—comparison between numerical (straight line) and experimental (dotted line) results for rigid cones (on the left) and for deformable structures—Thickness = 1 mm (on the right): (a) Cone $\beta = 6^{\circ}$; (b) Cone $\beta = 10^{\circ}$; (c) Cone $\beta = 14^{\circ}$.</figcaption></figure>

the level of the pressure peak should not vary along the body. Experiments clearly indicate that the peak intensity increases between the two probes. This phenomenon is not reproduced in the numerical simulation where the impact velocity remains constant.
Comparisons between numerical and experimental results as well as descriptions of the pressure distribution on elastic structures are a challenge. On the one hand, the depression zone cannot be reproduced numerically because the jet phenomenon has not been taken into account. On the other hand, because the second peak zone occurs after the period of time we are interested in, we have not reproduced it. Thus, the major difference between numerical results and an experimental test comes from the velocity of the peak progression. Indeed, we observed that the instant at which the peak pressure occurred is different in the experimental and numerical cases. In spite of this dissimilitude, the global numerical pressure distribution fits well. Moreover, we remarked (Tables 5 and 6) that the comparison is closer for the cases of large deadrise angle cones than for the smaller ones where water perturbations were more significant.

## 6. Conclusion

In this paper, we presented a numerical approach for solving two-dimensional and three-dimensional slamming problems. Firstly, we presented the fluid problem in terms of a velocity potential formulation. The asymptotic method, used to solve these equations, was described. A variational formulation together with a finite element method was adopted to solve the so-called Wagner problem. The wetting correction was obtained through an iterative procedure. The numerical resolution was validated with simple problems, such as the water impact problem of a wedge with small deadrise in two-dimensional or with a cone for the axisymmetric case. Good agreement between numerical and analytical results was achieved. Next, the structural governing equations were discussed and solved with the finite element method. The fully coupled problem was then resolved with the same method. This formulation introduced a coupling matrix. The correct evaluation of its terms was checked for the sloshing problem of a tank with an elastic beam for the two-dimensional case and with an elastic plate for the three-dimensional case.

The second part of this paper dealt with experiments. It described the series of drop-test experiments. These drop-tests consist in impact of rigid or deformable cones with a free surface initially at rest. Experimental equipment and results were described in detail. Experimental results, were well-predicted by the numerical pressure simulations, although some differences were observed. In particular, dissimilarity occurred between the peak pressure values (the first peak being weaker than the second), as a result of an increase in the velocity. The delay effect, resulting from impact loads, appeared later in the experiments. In conclusion, numerical results were well confirmed by the experimental ones.

## References

* Batchelor, G.K., 1967. Introduction to Fluid Dynamics. Cambridge University Press, Cambridge, UK.

* Chai, X.J., 1996. Influence de la gravité sur les interactions fluide-structure pour un fluide dans un domaine borné à surface libre. Ph.D. Thesis, Institut National Polytechnique de Lorraine, Nancy, France.

* Chuang, S.L., 1967. Experiments on slamming of wedge-shaped bodies. Journal of Ship Research 11, 190-198.

* Chuang, S.L., Milne, D.T., 1971. Drop tests of cones to investigate the three-dimensional effects of slamming. Navy Naval Ship Research and Development Center, Report 3543.

* Cointe, R., 1989. Two dimensional water solid impact. ASME Journal of Offshore Mechanics and Arctic Engineering 111, 109-114.

* Cointe, R., Armand, J.L., 1987. L'impact hydrodynamique en deux dimensions. 1erés Journées de l'Hydrodynamique, Nantes, France.

* Darrozes, J.S., François, C., 1982. Mécanique des fluids incompressibles. Lecture Notes in Physics, vol. 163. Springer, Berlin.

* Donguy, B., 2002. Etude de l'interaction fluide structure lors de l'impact hydrodynamique. Ph.D. Thesis, Ecole Centrale de Nantes, France.

* Donguy, B., Peseux, B., Fontaine, E., 2001a. Numerical and experimental investigation of a rigid cone striking the free surface of an incompressible fluid. Proceedings of ASME Pressure Vessels and Piping Conference, Atlanta, USA.

* Donguy, B., Peseux, B., Fontaine, E., Gornet, L., 2001b. Three-dimensional hydroelastic water entry: preliminary results. 11th International Offshore and Polar Engineering Conference, Stavanger, Norway.

* Fabula, A., 1957. Ellipse-fitting approximation of impact rigid bodies on water. Proceedings of the Fifth Midwestern Conference on Fluid Mechanics, pp. 299-315.

* Fontaine, E., Cointe, R., 1992. A second-order solution for the wedge entry with small deadrise angle. Seventh International Workshop on Water Waves and Floating Bodies, Val de Reuil, France.

* Fontaine, E., Cointe, R., 1997. Asymptotic theory of water entry. NATO Conference on High Speed Body Motion in Water, Kiev.

* Greenhow, M., Yanbao, L., 1987. Added masses for circular cylinders near or penetrating fluid boundaries-review, extension and application to water-entry, exit and slamming. Ocean Engineering 14, 325-348.

* Hagiwara, K., Yuhara, T., 1974. Fundamental study of wave impact loads on ship bow. Journal of the Society of Naval Architects of Japan 14, 73-85.
Howison, S.D., Ochenson, J.R., Wilson, S.K., 1991. Incompressible water entry problems at small deadrise angles. Journal of Fluid Mechanics 222, 215-230.

Korobkin, A.A., 1982. Formulation of penetration problem as a variational inequality. Dinamika Sploshnoi Sredy 58, 73-79.

Korobkin, A.A., 1985. Initial asymptotics solution of three-dimensional problem on a blunt body penetration in ideal liquid. Doklady Akademii nauk SSSR 283, 838-842.

Lamb, H., 1932. Hydrodynamics, sixth ed. Dover Publications, New York.

Magge, A., Fontaine, E., 1998. A coupled approach for the evaluation of slamming loads on ships. Proceedings of the Seventh International Symposium on Practical Design of Ships and Mobile Units, The Hague, Netherlands.

Morand, H.J.-P., Ohayon, R., 1992. Interactions fluides-structures. Recherche en Mathématiques Appliquées, Masson, Paris, France.

Moyroud, F., 1998. Fluid-structure integrated computational methods for turbomachinery blade flutter and forced response predictions. Ph.D. Thesis, Institut National des Sciences Appliquées de Lyon, France, Royal Institute of Technology of Stockholm, Sweden.

Moyroud, F., Cosme, N., Jocker, M., Franson, T.H., Lornage, D., Jacquet-Richardet, G., 2000. A fluid-structure interfacing technique for computational aeroelastic simulations. In: Ferrand, P., Aubert, S. (Eds.), Proceedings of the Ninth International Symposium on Unsteady Aerodynamics, Aeroacoustics and Aeroelasticity of Turbomachines. Lyon, France, pp. 721-738.

Nayroles, B., Touzot, G., Villon, P., 1992. Generalizing the finite element method: diffuse approximation and diffuse elements. Computational Mechanics 10, 307-318.

Peseux, B., 1989. Contribution à l'étude de structures partiellement immergées en matériau homogène ou en composite. Ph.D. Thesis, Ecole Nationale Supérieure de Mécanique de Nantes, France.

Peseux, B., Cartraud, P., Argouar'ch, Y., 1999. Dynamic response of partially filled tanks—explicite time integration method. In: Jerome, D.M. (Ed.), Proceedings of ASME Pressure Vessels and Piping Conference, PVP-vol. 394, pp. 191-199.

Scolan, Y.M., Korobkin, A.A., 2001. Three-dimensional theory of water impact. Part 1: Inverse Wagner problem. Journal of Fluid Mechanics 440, 293-326.

Tanizawa, K., 1998. A time-domain simulation method for hydroelastic impact problem. In: Proceedings of the Second International Conference on Hydroelasticity in Marine Technology, Fukuoka, Japan, pp. 119-127.

Tourbier, M., Peseux, B., Donguy, B., Gornet, L., 2002a. Modelling and simulation of the three-dimensional hydrodynamic problem. In: Moody, F.J. (Ed.), Proceedings of ASME Pressure Vessels and Piping Conference, PVP-vol. 435, pp. 173-179.

Tourbier, M., Peseux, B., Donguy, B., Gornet, L., 2002b. Numerical modelling of the three dimensional slamming problem. In: Païdoussis, M.P., et al. (Eds.), Proceedings of IMECE2002, ASME International Mechanical Engineering Congress, Exposition Fifth International Symposium on FSI, AE & FIV+N, IMECE2002-39043, pp. 1-9.

Van Dyke, M., 1975. Perturbation Method in Fluid Mechanics. Parabolic Press, Stanford, CA.

Verpeaux, P., Charras, T., Millard, A., 1988. Castem 2000: une approche moderne du calcul des structures. In: Fouet, J.M., Ladevèze, P., Ohayon, R. (Eds.), Calcul des Structures et Intelligence Artificielle, vol. 2. Pluralis, Paris, pp. 227-261.

Von Karman, T., 1929. The impact on a seaplane floats during landing. Technical Notes for National Advisory Committee for Aeronautics, N.A.C.A. TN321, Washington.

Wagner, H., 1932. Über stoss end an der oberfache von Flussigkeiten. Zeitschrift für Angewandte Mathematik und Mechanik 12, 193-215.

Wanatabe, I., 1986. On hydrodynamic impact pressure acting upon flat bottomed ships. Journal of the Society of Naval Architects of Japan 159, 81-92.

Wilson, S., 1989. The mathematics of ship slamming. Ph.D. Thesis, University of Oxford, UK.

Zhao, R., Faltinsen, O.M., 1993. Water entry of two-dimensional bodies. Journal of Fluid Mechanics 246, 593-612.

Zhao, R., Faltinsen, O.M., Aarnes, J.V., 1996. Water entry of arbitrary two-dimensional sections with and without flow separation. In: Proceedings of the 21st Symposium on Naval Hydrodynamics, Trondheim, Norway, pp. 408-423.

Zienkiewicz, O.C., Taylor, R.L., 1991. The Finite Element Method, fourth ed., vol. 2: Solid and Fluid Mechanics, Dynamics and Non Linearity. McGraw-Hill, New York.