## Abstract

The free surface capturing method based on Cartesian cut cell mesh is extended to the water-entry and -exit fields with body–fluid interaction. The governing equations are the incompressible Euler equations for a variable density fluid system with a free surface, which is treated as a contact discontinuity. The solver is based on the artificial compressibility method with a dual time-stepping technique for time advancing and the finite-volume method with a high-resolution Godunov-type upwind scheme on spatial discretization. For solving the numerical problem caused by the extension, the acceleration term of body is introduced for the new pressure condition of moving solid boundary, and the exact Riemann solution is used to calculate the flux of the solid boundary. In addition, a new solution of gradients based on the least-squares idea is presented for simply calculating the gradients. Finally, test cases show that the present method is available and can be successfully applied in various water-entry and -exit phenomena.

## 1. Introduction

The slamming (water impact) phenomenon is of concern in many marine applications, including different kinds of slamming on ships, sloshing, green-water impact on deck, ship or platform launch and offshore installation. During impact between a body and water, impulse loads with high-pressure peaks occur, which can result in local damage and even global trajectory change for some marine structures. Owing to its practical importance in the marine engineering field, the water impact problem has been widely studied by many scientists and engineers using various methods. Therein, the two-dimensional method with strip theory for many marine structures is not only simple and practical but also a theoretical basis for three-dimensional study, and so it is still necessary to study the two-dimensional water impact models.

The important pioneering works for the two-dimensional water impact by Von Karman (1929) and Wagner (1932) laid a theoretical foundation for further studies. With the development of power and capacity of computers, many computational fluid dynamics methods get substantial improvement. Therein, the boundary-element method (Zhao & Faltinsen 1993; Zhao *et al*. 1997) is a representative method, by which the load on the body and free surface shape can be calculated. But because the boundary-element method is based on potential, assume the rotational motion near free surface cannot be well simulated. In the present study, in order to predict the force and motion of the body and achieve the fluid information near the free surface, the free surface capturing method based on Cartesian cut cell mesh is improved and used to simulate water impact problems.

Compared with traditional methods to treat the free surface such as the free surface fitting method (Glimm *et al*. 1988; Lemos 1992; Farm *et al*. 1994) and the free surface tracking method (Hirt & Nichols 1981; Youngs 1982; Kleefsman *et al*. 2005), the free surface capturing method as a novel method was proposed initially by Kelecy & Pletcher (1997) and developed by Pan & Chang (2000) and Qian *et al*. (2006). In this approach, both liquid and gas fields are solved simultaneously, and the location of the free surface is captured as a contact discontinuity in the density field. The enforcement of a conservation law eliminates the need for complex surface tracking and reconstruction procedures.

On the other hand, because the Cartesian cut cell mesh (Coirier & Powell 1993; Yang *et al*. 1997*a*,*b*; Causon *et al*. 2001; Qian *et al*. 2006) has the advantage to treat the moving bodies by updating a few cut cells locally rather than remeshing the whole-flow domain, here this mesh system is taken as an alternative to structured and unstructured grid methods for the problem with moving boundaries.

Based on the free surface capturing method and Cartesian cut cell mesh, this paper employs the two-dimensional incompressible Euler equations with a variable density fluid as governing equations in §2. Then the finite-volume method is applied in spatial discretization and Roe's approximate Riemann solver (Qian *et al*. 2006) is used to calculate the numerical flux in §3*b*(i). The artificial compressibility method and dual time-stepping technique (Soh & Doodrich 1988; Beddhu *et al*. 1994) are adopted for time advancing in §3*c*.

The above schemes originated from Qian *et al*. (2006). When the original method is extended to the water-entry and -exit problems with body–fluid interaction in §4 the numerical problem appears, which is described in §5*e*. In order to solve the problem, as the improvement of the original method, here the acceleration term is introduced for the new pressure condition of solid boundary in §2 and the exact Riemann solution is computed to solve the boundary flux in §3*b*(iii). In addition, a new simplified solution for gradients based on least-squares is presented in §3*b*(ii).

Finally, five test cases of water entry and exit are designed in §5. The first two cases of water entry with a constant velocity are the fundamental numerical tests, in cases (iii) and (iv), the present method is successfully applied in the water-entry and -exit problems with body–fluid interaction, and the last test directly shows the advantages of the present method compared with the original method (Qian *et al*. 2006).

## 2. Mathematical model

For the two-dimensional incompressible, unsteady, inviscid fluid system with variable density field, the governing equations can be written in conservation form as follows.

The mass conservation equation(2.1)The *x*-direction momentum equation(2.2)The *y*-direction momentum equation(2.3)The incompressibility constraint equation(2.4)where *ρ* is the fluid density; *u*, *v* represent the *x*- and *y*-direction fluid velocity, respectively; *p* is the fluid pressure; and *B*_{x}, *B*_{y} are the *x*- and *y*-directional acceleration of body force.

For better use of the artificial compressibility method, a pressure term is introduced into the incompressibility constraint equation. And the equations (2.1)–(2.4) are modified and rewritten as the integral conservation form on Cartesian cut cell mesh system(2.5)where *I*_{ta}=diag[1, 1, 1, 0]; *Ω* is the domain of control volume; *S* is the boundary surrounding *Ω*; ** n** is the unit vector in outward normal direction of

*S*;

*Q=*[

*ρ*,

*ρu*,

*ρv*,

*p*/

*β*]

^{T}is the fluid vector of conserved variables; is the inviscid flux function through

*S*, where

*f*

^{I}=[

*ρ*(

*u*−

*u*

_{b}),

*ρu*(

*u*−

*u*

_{b})+

*p*,

*ρv*(

*u*−

*u*

_{b}),

*u*]

^{T};

*h*

^{I}=[

*ρ*(

*v*−

*v*

_{b}),

*ρu*(

*v*−

*v*

_{b}),

*ρv*(

*v*−

*v*

_{b})+

*p*,

*v*]

^{T};

*B*=[0, 0, −

*ρg*, 0]

^{T}is the source term by assuming that the only body force is gravity;

*β*is the coefficient of artificial compressibility method;

*g*is the gravitational acceleration;

*u*

_{b}and

*v*

_{b}are

*x*- and

*y*-direction velocity of the boundary

*S*, which are 0 when the boundary is stationary.

In the present study, various kinds of boundary conditions can be classified.

*Outlet or open boundary*. A zero gradient condition is applied in the velocity and density, and the pressure at this boundary is fixed, which allows fluid to enter or leave the computational domain freely according to the local flow velocity and direction.*Solid body boundary*. The no-penetration condition can be applied to the velocity, and the density is assumed to have a zero normal gradient. For the pressure boundary condition, here the momentum equations (2.2) and (2.3) are projected on the normal direction of moving solid boundary and can be rewritten as(2.6)where is the fluid velocity at the solid boundary; is the normal unit vector of solid edge; and*u*_{fn}=*un*_{x}+*vn*_{y}is the normal projection of fluid velocity.At the solid boundary, the normal projection

*u*_{sn}of the body velocity is equal to*u*_{fn}, so equation (2.6) can be described as(2.7)Assume that the rigid body moves with one degree (up or down) with an acceleration and*v*_{S}is the vertical velocity of body with positive direction upwards. Thus, the second term on the left-hand side in equation (2.7) is 0 and the pressure boundary condition is(2.8)

## 3. Numerical method

### (a) The Cartesian cut cell mesh

To deal with the moving boundaries, here the Cartesian cut cell mesh system is used, which can be generated by cutting solid regions out of a background Cartesian mesh. Thus, three kinds of cells are formed, namely the cut cell, the fluid cell and the solid cell.

In order to be convenient for deleting and inserting the grid cell, here the adjacent cells are connected in the form of a linked list. Furthermore, the data in each cell are stored together as a whole, including geometry and fluid information. This kind of data structure in the cell is called the structural body. Owing to space limitation, the steps to generate the Cartesian cut cell are not discussed in detail here. For more information, refer to Coirier & Powell (1993), Yang *et al*. (1997*a*,*b*), Causon *et al*. (2001) and Qian *et al*. (2006).

The code of the two-dimensional Cartesian cut cell is programmed by the authors and should be verified. Here we take the circle, wedge and NACA2415 as examples (as shown in figure 1).

Figure 1 is the amplified sketch of the Cartesian cut cell close to the body surface. As the figure shows, boundary-fitted meshes around the body boundary can be produced by using the authors' code; furthermore, the data storage form and the scheme of generating the Cartesian cut cell in the program can be considered to be right.

### (b) Spatial discretization

In the present study, the finite-volume method (Causon *et al*. 2001) has been applied in numerical spatial discretization for equation (2.5), and the discretization equation can be written as (a case study of *cell* (*i*, *j*))(3.1)where *Q*_{i,j} is the average conservation variable that is stored at the centre of *cell* (*i*, *j*); ∂*Ω*_{i,j}, *A*_{i,j} are the boundary and area of *cell* (*i*, *j*); and *R*(*Q*_{i,j}) denotes the spatial term.

On the right side of equation (3.1), the first term is called a convection term and should be discretized as , where *m* is the number of edges in *cell* (*i*, *j*); Δ*l*_{k} is the length of the edge *k*; and is the inviscid flux through outward normal direction of edge *k*.

#### (i) Flux evaluation

The formula (3.1) represents a set of nonlinear hyperbolic equations, which may produce a discontinuous solution. Here, the Roe's approximate Riemann solver as a Godunov-type scheme is used for the reconstruction of the flux , because it has the ability to achieve both discontinuous and smooth solutions without generating spurious oscillations.

Assume a one-dimensional Riemann problem in the direction normal to the cell edge, and Roe's approximate flux is adopted to replace the real inviscid flux (3.2)where *Q*_{R} and *Q*_{L} are the reconstructed conservation variables on the right or left side of edge *k* with anticlockwise direction; denotes the Jacobian matrix of the flux ; is the approximate Jacobian matrix evaluated by Roe's average state on edge *k*; and are the right eigenvectors and the eigenvalue matrix of , respectively; and is the diagonal matrix constructed by the absolute value of . The Jacobian matrix *A* can be deduced according to the definition formula (3.3)where *u*_{b} and *v*_{b} are the *x*- and *y*-direction velocity of cell edges and is the unit vector in outward normal direction of cell boundaries.

For the Roe's approximate Riemann solver, the appropriate average on the cell edge is obtained by the following equation (3.4) for the two-fluid system:(3.4)then is introduced into equation (3.3) and the Jacobian matrix can be achieved.

#### (ii) Calculation of gradients and reconstruction technique of *Q*_{R}, *Q*_{L}

Before the flux is calculated, the variables *Q*_{L} and *Q*_{R} on both sides of edge should be determined. To achieve second-order spatial accuracy, a piecewise linear representation of the cell variables must be reconstructed according to the data stored in the cell centre. Meanwhile, to minimize numerical diffusion and improve resolution at the interface, the primitive variable vector *U*=[*ρ*, *u*, *v*, *p*]^{T} is applied in the linear reconstruction to replace the conservation variable *Q* (Qian *et al*. 2006), so the linear representation can be written as(3.5)where Δ** r** is the position vector from the cell centre to any point (

*x*,

*y*) within cell and

*U*

_{C}, ∇

*U*

_{C}are the primitive variable and gradient stored at the cell centre.

In past studies (Yang *et al*. 1997*a*,*b*; Causon *et al*. 2001; Qian *et al*. 2006), the gradient ∇*U*_{C} in the cut cell is calculated based on the geometry shape of the cut cell. This method has been validated, but there are so many classifications of the cut cell that it is a little tedious to programme the computer code.

Because the physical concept and mathematical model of the least-squares method are quite definite, simple and convenient for programming, here a new method based on the least-squares idea is adopted to calculate the gradients. Take the case in figure 2 as an example to describe the details of the new method.

For the given cell (

*i*,*j*), the beginning and critical work is creating a node stencil which is used to calculate the gradients. The selection standard of the stencil is described as follows: a constant*d*_{0}is prescribed, and then the distance from neighbour cell centre to the given cell centre is calculated. If the distance is less than*d*_{0}, the neighbour cell should be chosen into the node stencil of the given cell. In figure 2, based on the*d*_{0}as shown, the neighbour cell (*i*,*j*−1), (*i*−1,*j*), (*i*,*j*+1) and (*i*+1,*j*) are selected into the node stencil of the given cell (*i*,*j*).For the cut cell, the influence of the moving boundary should be taken into account and the concept of the fictional cell is introduced into the node stencil. The central position of fictional cell

*R*is determined as the mirror symmetric of the given cut cell centre about the solid boundary and the variables in the fictional cell*R*can be obtained by (figure 2)(3.6)(3.7)(3.8)where*ρ*_{i,j},*U*_{i,j},*p*_{i,j}are the fluid variables in the given cut cell; is the body velocity, here assume*u*_{S}=0 and*v*_{S}is the vertical velocity of body; and is the unit vector in the outward normal direction of body boundaries.After the suitable node stencil is determined, the information of each node in the stencil including the coordinates and fluid variables is introduced into equation (3.5). Then a linear least-squares system is formed. By using the iterative refinement algorithm that is described by Golub & van Loan (1983), the least-squares system is solved and the unknown gradients

*U*_{x}and*U*_{y}are calculated.In order to maintain monotony of the numerical scheme and control spurious oscillations, a parameter called the ‘limiter’ is introduced into equation (3.5)(3.9)The chosen limiter

*Φ*=min(*Φ*_{j}(*r*_{j})) is a constant in the given cell;*j*denotes each cell edge; and*Φ*_{j}(*r*_{j}) is the limiter function. Here*Φ*_{j}(*r*_{j})=max[min(2*r*_{j}, 1), min(*r*_{j}, 2)] called the Roe's Superbee limiter is taken, and*r*_{j}can be calculated by(3.10)where , ;*U*_{j}is the unlimited value on the cell edge; and*U*_{0}is the value at cell centre.

Finally, a case to test the new method is designed in §5*a*, for comparison with the solution of Qian *et al*. (2006) for gradients.

#### (iii) The Riemann solutions on a moving boundary

Based on the reconstructed fluid variables at the solid boundary and the body velocity, how to calculate the Riemann variables and flux on a moving boundary is described in this section. Here, the Riemann solutions on the edge are defined with the superscript ‘*’, which are normal and tangential velocity , , density *ρ*^{*}, pressure *p*^{*} and normal flux .

Take a cut cell with a moving boundary as shown in figure 3 as an example; at the solid boundary, *U*_{f}, *ρ*_{f}, *p*_{f} are the reconstructed fluid velocity, density and pressure, respectively, and *U*_{S} is the velocity of solid boundary. Now project both *U*_{f} and *U*_{S} into the normal direction ** n** and tangential direction

**of the solid boundary and generate**

*t**u*

_{fn},

*u*

_{ft},

*u*

_{sn}and

*u*

_{st}.

In the tangential direction, any difference between *u*_{ft} and *u*_{st} can be treated as a shear wave superimposed on a contact discontinuity.

In the normal direction, the Riemann solution for a moving piston is used and can be determined,(3.11)

The type of moving wave near the solid boundary is related to the relative magnitude of *u*_{fn} and *u*_{sn}, so before *p*^{*} is calculated, the two velocities *u*_{fn} and *u*_{sn} should be compared.

If *u*_{sn}≤*u*_{fn}, both left and right moving waves are shocks. Thus, *p*^{*} can be obtained from the normal shock relations deduced in detail in appendix A(3.12)

If *u*_{sn}>*u*_{fn}, both left and right moving waves are rarefactions. Then *p*^{*} can be achieved from the rarefaction relations, which can be referred to in appendix A(3.13)where is the artificial speed of sound; *β* is the artificial compressibility coefficient that has also been defined in equation (2.5) in §2.

Therefore, after *p*^{*} is calculated by using equations (3.12) or (3.13), the inviscid flux through the outward normal direction of the solid boundary in the cut cell is obtained(3.14)where ** n**, Δ

*l*are the unit vector through the outward normal direction and the length of solid edge in the cut cell, respectively.

### (c) The dual time-stepping technique

To achieve a time-accurate solution for unsteady flow problems, based on the artificial compressibility method, here the implicit dual-time iteration technique is used. A fictitious time derivative of variables is introduced into equations and the solution at each real time step is obtained by solving a steady-state problem in the pseudotime domain.

The first-order Euler implicit difference scheme and a pseudotime derivative of variables are employed for equation (3.1)(3.15)where *t* is the physical time; *τ* is the pseudo time; and the spatial term *R*(*Q*^{n+1,m+1}) is linearized according to Kelecy & Pletcher (1997).

If the edge *k* is the fluid boundary, the linearization of should be(3.16)where [*A*^{+}] and [*A*^{−}] are defined as(3.17)

(3.18)

If the edge *k* is the solid edge, the linearization of should be(3.19)and(3.20)where *ρ*^{*}, *u*^{*}, *v*^{*}and *p*^{*} are the Riemann variables on the solid edge as defined in §3*b*(iii).

For the fluid or solid edge, the linearization of body force can be represented as(3.21)where(3.22)

For a given cell, combine equation (3.16) or (3.19) with (3.21) and then construct the equation system that can be written in matrix form as(3.23)where *D* is a block diagonal matrix; *L* is a block lower triangular matrix; *U* is a block upper triangular matrix; and *RHS* is a vector that is independent of time. In matrix *D*, *L* and *U*, the data storage form for each cell is a 4×4 matrix.

An approximate LU factorization scheme (Pan & Lomax 1988) is adopted for equation (3.23), and the inverse of the equation can be written as(3.24)Equation (3.24) can be solved by introducing an intermediate variable *δQ*^{*}, and at each pseudo time step, *δQ*^{S} should be calculated by forward and backward sweeps(3.25)(3.26)For each physical time step, according to equation (3.15), when *δQ*^{S} is iterated to 0, i.e. the ((*QA*)^{n+1,m+1}−(*QA*)^{n+1,m})/Δ*τ* tends to 0, the density, momentum and incompressibility constraint equations are satisfied. Thus, here the norm of *δQ*^{S} vector is selected as the convergence standard of subiteration. Subiteration is not terminated until *L*_{2} is less than a specified limit *ϵ* that is equal to 10^{−4} here.

## 4. Interactive body–fluid motion

The present research on the water-entry and -exit problems include not only that the solid body moves with a prescribed velocity but also that the motion of the solid body is calculated from the body–fluid interaction.

In order to ensure enough numerical stability, here a fully body–fluid coupling method for the body velocity proposed by Kleefsman *et al*. (2005) is taken as follows.

The new body velocity is calculated by(4.1)where Δ*t* is the physical time step; *m* is the mass of body; and is the vertical body velocity at the physical time *n*. *ω* is a relaxation parameter to control the stability of subiteration, which is in the range from 0 to 1. The larger the *ω* is, the faster the convergence rate of subiteration is.

The equation (4.1) describes a subiterative process of body–fluid coupling at each physical moment. At the physical time step *n*, the iterative initial value is taken as . For the subiterative *k* step, based on the vertical body velocity , a new pressure field is computed and then is obtained by direct integration of the new pressure over the body surface. When is less than 10^{−3}, the subiteration is convergence and is taken as the new velocity boundary condition of the solid body for the next physical time step *n*+1.

The new vertical displacement of the solid body is calculated by(4.2)where is the new displacement condition of the solid body for the next physical time step *n*+1.

## 5. Results of numerical tests and discussions

In this section, various test cases of water entry and exit have been simulated to validate the numerical accuracy of the current method.

Water entry of a two-dimensional rigid wedge with a constant velocity.

Water entry of a two-dimensional rigid circular cylinder with a constant velocity.

Water entry of a two-dimensional free-falling rigid wedge.

Water exit of a two-dimensional free-rising rigid circular cylinder.

Comparison of numerical results corresponding to different pressure solutions in the cut cell between this paper and Qian

*et al*. (2006).

Here the density ratio of water to air is taken as 1000 : 1 and viscous effects have been ignored. The gravitational acceleration is taken as 9.81 m s^{−2}. The location of free surface was determined as the density contour with the average value of the two fluid densities *ρ*_{0}=(*ρ*_{1}+*ρ*_{2})/2, and all free surface plots were created by using 10 evenly spaced contours between 0.9*ρ*_{0} and 1.1*ρ*_{0}. All the calculations presented here were performed on a computer with Intel Pentium (R) D CPU 3.00 GHz and 0.99 GB physical memory. In the numerical simulations, it was found that the artificial compressibility coefficient *β* value between 200 and 2000 can produce good results for the following test cases.

### (a) Water entry of a two-dimensional rigid wedge with a constant velocity

In this test case, the basic computation domain was a square of 3×3 m occupied half by water and half by air. The water-entry cases can be designed as

a rigid wedge with 45° deadrise angle (the angle between the wedge slope and the horizontal or undisturbed water surface) 0.4 m wide and 0.2 m high,

a rigid wedge with 30°, 0.6928 m wide and 0.2 m high,

a rigid wedge with 60°, 0.2309 m wide and 0.2 m high, and

a rigid wedge with 81°, 0.0634 m wide and 0.2 m high.

The zero time (*t*=0) of water entry was chosen as the moment that the apex of wedge contacts the initial calm free surface, and the prescribed velocity of wedge *V* was 1 m s^{−1}. The physical and pseudo time steps Δ*t*, Δ*τ* were 0.0001 and 0.01 s, and *β* was taken as 500.

In order to investigate the influence of grid refinement, case (i) is simulated and the results at *t*=0.09 s are plotted in figures 4 and 5.

From figure 4 it can be found that whether the mesh is coarse or refined it has an influence on the profile of the free surface. The finer the grid, the better the jets around wedge are resolved. In figure 5, *y* is the vertical coordinate of wedge edges, *C*_{p}=(*p*−*p*_{a})/0.5*ρV*^{2} is the pressure coefficient and *p*_{a} is the atmospheric pressure. When the grid size reduces to a certain extent, the pressure distributions for different meshes are approximately considered to be the same. Thus, based on calculation precision and speed, the grid discretization scheme 300×300 (Δ*x*=Δ*y*=0.01 m) should be applied in this numerical test. In addition, approximately 9 hours of CPU time was spent on advancing 2000 time steps by using the 300×300 grid scheme.

The results of pressure distribution on the edge of wedge with 45° at *t*=0.08, 0.10 and 0.12 s are plotted in figure 6. In the figure it is found that the numerical solutions at different water-entry times are almost identical, which shows that during the whole-water entry, the self-similar characteristic of the wedge can be well simulated by using the present method. So in the following calculations of this test, it is feasible to select a classic moment of water entry to illustrate the behaviour of the wedge.

To test the method for gradients based on the least-squares idea in §3*b*(ii), case (i) was simulated. As the calculation advances 1950 time steps, the geometry information and the calculated pressure field near *cell* (*i*, *j*) are shown in figure 7 and in table 1.

In table 1, from *cell* (*i*, *j*) to *cell* (*i*, *j*+1), the ‘*pressure*’ values have been calculated and are used as known conditions for pressure reconstruction. The *pressure* value of *R* is obtained by the equation (3.8). For the *point* 1–4, the ‘*pressure*1’ values are reconstructed by the method in §3*b*(ii), the ‘*pressure*2’ values are computed by the solution of Qian *et al*. (2006) and the ‘*error*’ values are calculated by *error*=|*pressure*1−*pressure*2|/*pressure*2. As shown in table 1, the error is very small and the results calculated by the two methods from this paper and from Qian *et al*. (2006) are a little different from each other, which shows that, based on the calculation precision requirement, the new solution for the reconstruction of fluid variables can replace the method of Qian *et al*. (2006).

In figure 8, *y* is the vertical coordinate of the wedge edge and *C*_{p}=(*p*−*p*_{a})/0.5*ρV*^{2} is the pressure coefficient. The pressure distributions on the wedge edge with different angles at *t*=0.10 s are plotted by using the numerical method of this paper, along with the results of the analytic solution (Mei *et al*. 1999), the similarity method (Dobrovol'skaya 1969) and the nonlinear numerical simulation (Zhao & Faltinsen 1993). In the four methods, the jet flow by the analytic solution (Mei *et al*. 1999) cannot be modelled because in this method a single-valued height function is used to locate the free surface. But by the authors' method, the jet can be well simulated, and the location and magnitude of the maximum impact pressure agree better with the solution of Zhao & Faltinsen (1993) than with the similarity solution (Dobrovol'skaya 1969).

### (b) Water entry of a two-dimensional rigid circular cylinder with a constant velocity

This test model was based on a square domain of 2×2 m occupied half by water and half by air, and then the grid scheme Δ*x*=Δ*y*=0.01 m was used. A cylinder with 0.1 m radius entered water with the constant velocity 1 m s^{−1}. The zero time (*t*=0) of this case was selected as the moment that the bottom of the cylinder is 0.1 m up from the initial free surface. The physical and pseudo time steps Δ*t* and Δ*τ* were 0.0001 and 0.01 s, *β* was 500 and the computational time was approximately 18 hours using the 200×200 grid scheme and 10 000 time steps.

In figure 9, *C*_{S}=*F*/(*ρRV*^{2}) is the slamming coefficient; *F* is the vertical hydrodynamic; *V* denotes the velocity of the body; and *R* is the radius of the circular cylinder. In the present method, the hydrodynamic *F* on the cylinder can be obtained by direct integration of the pressure over the body surface and is plotted in the figure, along with the experiment data (Campbell & Weynberg 1980) and classic theories (Von Karman 1929; Wagner 1932). Compared with the results from Von Karman and Wagner, the authors' solutions are closer to the experiment data (Campbell & Weynberg 1980). But it can be found that by the present method, the moment corresponding to the maximum slamming coefficient *C*_{S} lags a little behind *t*=0.10 s. The reason may be that before the cylinder as a blunt body contacts the initial calm free surface, due to the effect of airflow in the air-gap, the deformation of the free surface has been produced, which will delay the moment of water entry.

Here, during water entry, the free surface profile and velocity vectors at various moments are plotted in figure 10 to describe the water-entry phenomenon of the two-dimensional cylinder.

At *t*=0.00 s, the initial positions of the cylinder and the calm free surface are shown in the first graph. Then the cylinder begins to move downwards in air, and the air near the body surface is in motion together with the cylinder. When the cylinder is close to the free surface (*t*=0.09 s), the air-gap between the body surface and free surface is relatively small. Thus, the influence of the air in the air-gap on the free surface becomes very important, which causes the deformation of the free surface. That is why the cylinder is also in air at *t*=0.10 s. Meanwhile, the air escapes from the air-gap with 5–6 m s^{−1} high speed. As the cylinder pierces the surface (*t*=0.11, 0.15 s), the water rises up along the each side of the body with 1–3 m s^{−1} high velocity and the displacement waves are formed. As the body travels further into the water (*t*=0.28, 0.32 s), the two intersections between each edge of the body and the free surface continue to move up to the top of the cylinder and then towards each other. After the free surface is reconnected (*t*=0.38, 0.43, 0.53 s), the water particles near the solid boundary are forced to move together with the body and flow through both sides to the rear point of the cylinder. In addition, a jet flow is created upwards above the cylinder. Then due to the effect of gravity the jet flow begins to move downwards (*t*=0.63 s) and almost vanishes (*t*=0.93 s).

### (c) Water entry of a two-dimensional free-falling rigid wedge

The test model was created on a square domain of 3×3 m occupied half by water and half by air, and the wedge with 30° was 0.5 m wide, 0.29 m high and 241 kg weight. The zero time (*t*=0) of this case was selected as the moment that the apex of the wedge contacts the initial calm free surface. Here for comparison with experimental data of Zhao *et al*. (1997), the initial velocity of the wedge was prescribed as −6.15 m s^{−1}, and then the wedge began to fall freely. The physical and pseudo time steps Δ*t* and Δ*τ* were 0.0001 and 0.01 s, the grid scheme Δ*x*=Δ*y*=0.01 m was taken, *β* was 500, *ω* was taken as 0.8, and this numerical test advancing 4000 time steps by using 300×300 grid scheme took up approximately 34 hours CPU time.

In this test case, the motion of wedge need to be calculated from the body–fluid interaction as described in §4. As shown in figure 11 compared with the numerical solutions of Kleefsman *et al*. (2005), the results of vertical force and velocity by the present method are in better agreement with the experimental data of Zhao *et al*. (1997). But it is a pity that by this method the absolute value of vertical velocity is still a bit underestimated, which may be caused by the effect of different dimensions between numerical simulation of this paper (two dimensions) and experiment (three dimensions).

By the present method during the whole process of water entry, the time history of force and velocity can be simulated and the results are plotted in figure 12.

### (d) Water exit of a two-dimensional free-rising rigid circular cylinder

A square domain of 10×10 m occupied half by water and half by air was selected as the basic test model. The radius of the circular cylinder was 1.0 m and the weight was approximately 2820.9 kg. At the moment *t*=0, under the calm free surface, the cylinder was placed 5.0 m distance from the centre of the cylinder to the water surface. Because the buoyancy of the body was larger than the gravity, the cylinder rose up freely from the initial resting position. The physical and pseudo time steps Δ*t* and Δ*τ* were 0.0005 and 0.05 s, *β* was taken as 500, the grid scheme Δ*x*=Δ*y*=0.1 m was adopted, *ω* was taken as 0.8, and the computational time was approximately 52 hours by using the 100×100 grid scheme and advancing 48 000 time steps.

First the dimensionless parameters in figures 13 and 14 should be discussed. *R* denotes the radius of cylinder; *V* is the vertical velocity of cylinder; *a* is the vertical acceleration of cylinder; *t* is the physical time; the dimensionless time *T*=*t*(*g*/*R*)^{1/2}; the Froude number *F*_{r}=*V*/(*gR*)^{1/2}; and the dimensionless ratio of acceleration *b*=*a*/*g*.

The curves of velocity and acceleration are basically consistent between the present method and the nonlinear calculation of Moyo & Greenhow (2000) from figure 13.

But there are still two differences for the two methods. One is that near the free surface, the acceleration and the moment versus the maximum velocity are different, and the other is whether the calculation will break down or not. For the first difference, as shown in figure 13, the acceleration by the present method is larger than that in nonlinear method near the free surface and the time versus the maximum velocity delays. That is because at the same time by the present method the distance between the cylinder and free surface is larger, which can be shown in figure 15. For the second difference, in the nonlinear calculation of Moyo & Greenhow (2000), the Rayleigh–Taylor instability can lead to termination of the numerical calculation when the cylinder is almost totally above the initial free-surface position. But in the present method, the calculation will not break down as shown in figure 14, and it is found that the cylinder may vibrate up and down.

When the cylinder was close to the free surface, the deformation of the free surface was produced. As shown in figure 15, the thickness between the free surface and cylinder by the present method is larger than that by the nonlinear method at the same time.

### (e) Comparison of numerical results corresponding to different pressure solutions in the cut cell between this paper and Qian *et al*. (2006)

Because the authors' method is an improvement of the method by Qian *et al*. (2006), here the two methods should be compared with each other using two cases, i.e. the case (ii) in §5*b* and the case (iii) in §5*c*.

First, the difference of the two methods should be described. In the method of Qian *et al*. (2006), the pressure boundary condition is and the flux function of body boundary is directly calculated by the reconstructed fluid variables on solid boundary. But in this paper, the pressure boundary condition is in §2, and the flux function of body boundary is calculated by using the Riemann solution on the solid edge which is described in §3*b*(iii). In figures 16–18, the results by the method of Qian *et al*. (2006) are expressed by ‘B’. On the other hand ‘A’ is the solution of this paper.

As shown in figure 16, the results from the two methods agree well with each other, which shows that if the solid body moves with a constant velocity, the influence of body acceleration and moving waves can be neglected. However, for the case of a free-falling wedge, the results of the two methods are different from each other. In figure 17 compared with B, the results of force and velocity by A are in better agreement with the experiment data. Furthermore, as shown in figure 18 by using method B, when the physical time advances to 0.0760 s, the force curve begins to oscillate with high frequency. Then until 0.0855 s the calculation is terminated and the cell with negative density appears, which shows that the calculation has been divergent. But by the present method A, the calculation will not break down.

The above results show that for the case of free-falling water entry, the influence of the body acceleration and moving wave becomes important, and neglecting this influence will cause error. If the error accumulates to a certain extent with time stepping, the numerical calculations may be divergent and stop.

## 6. Conclusions

The knowledge about the projection method of momentum equations and the Riemann solution for a moving piston is introduced into the two-dimensional incompressible Euler equations of a two-fluid system with a variable density field, and then the new pressure solution on moving solid boundary including the effect of the moving waves and body acceleration is deduced, and successfully applied in various water impact models.

In addition, based on the least-squares theory, this paper presents a new calculation method for gradients. Owing to eliminating the consideration for the detail of different cut cells, it is convenient for programming to adopt this method.

Finally, the results in test cases (i)–(v) (§5*a*–*e*) show that the present method has a good ability to simulate two-dimensional water-entry and -exit problems.

Furthermore, the method can be developed further on the following aspects:

introducing the concept of adaptation based on the quadtree refinement algorithm into the Cartesian cut cell mesh system to improve computing capability,

being closer to the practical problems by solving the N–S equations for validating the viscous effect on the calculation precision and time,

extending the present method to three dimensions for more extensive application about ocean engineering problems, and

taking the structural elastic deformation of a solid body into account and combining this method with structural mechanics for the hydroelastic problem.

## Acknowledgments

The authors are grateful to Dr Qian of the Manchester Metropolitan University in the UK for providing many very useful theoretical help and suggestions about the numerical method. The results of the analytical and numerical methods are kindly provided by Mei of MIT in USA, Dobrovol'skaya of CCASU in Moscow, Moyo of Brunel University in the UK, Kleefsman of Groningen University in the Netherlands and Zhao of MARINTEK in Norway. The authors also would like to thank Faltinsen of NUST in Norway, Campbell of Southampton University in the UK and Greenhow of Brunel University in the UK for experimental data. This work was financially supported by the National 863 Plan Foundation under grant no. 2006AA 09A 104 (China).

## Appendix A

Here with reference to Yang *et al*. (1997*a*,*b*) and Causon *et al*. (2001), take the one-dimensional governing equations to demonstrate how the Riemann pressure on a moving solid boundary is deduced.

The form of one-dimensional governing equations is shown as(A1)

(A2)

(A3)

If *u*_{sn}≤*u*_{fn}, both left and right moving waves are shocks, and *p*^{*} can be obtained from the normal shock discontinuous relations (Ranking–Hugoniot relations)(A4)(A5)(A6)where *D* is the shock velocity.

According to (A 4) and (A 5), equation (A 7) can be achieved(A7)Integrate (A 4) and (A 6) into one equation as(A8)Combine (A 7) with (A 8), then eliminate *ρ*^{*}(A9)(A10)where is the artificial speed of sound.

If *u*_{sn}>*u*_{fn}, both left and right moving waves are rarefactions. Then *p*^{*} can be achieved from the rarefaction relation.

According to equations (A 1)–(A 3), the governing equations with primitive variables can be rewritten in the matrix form as(A11)where(A12)and the column vector *V*=(*ρ*, *u*, *p*)^{T}.

The eigenvalues of *A* are *λ*_{1}=*u*, and . The eigenvector versus *λ*_{1} is (1, 0, 0)^{T}, the eigenvector corresponding to *λ*_{2} is (−*ρ*^{2}*λ*_{2}/*β*, 1, −*ρλ*_{3})^{T} and for *λ*_{3} the eigenvector is (−*ρ*^{2}*λ*_{3}/*β*, 1, −*ρλ*_{2})^{T}. Thus, the matrix *A* can be diagonalized as *A*=*TΛT*^{−1}, where *Λ*=diag(*λ*_{1}, *λ*_{2}, *λ*_{3}) and matrix *T* is formed according to the three eigenvectors. Then the rarefaction relations can be obtained from the total differentiation for the above mentioned diagonalizable matrix and *p*^{*}can be approximately achieved as(A13)where is the artificial speed of sound.

## Footnotes

- Received October 15, 2008.
- Accepted February 20, 2009.

- © 2009 The Royal Society