# Further Application of Surface Capturing Method and Cartesian Cut Cell Mesh on Hydroelastic Water-Entry Problems of Free-Falling Elastic Wedge

Authors: Wen-hua Wang; Yi Huang; Yan-ying Wang
Journal: Mathematical Problems in Engineering
 (2014)

---

## Abstract

In order to study the interactions between fluid and elastic structure (such as marine lifeboat falling down and ship), this paper presents a new CFD method on hydroelastic water-entry problem of free-falling elastic wedge, which can more conveniently handle moving solid boundaries. In the CFD solver, a surface capturing method and the Cartesian cut cell mesh are employed to deal with the moving free surface and solid boundaries, respectively. On the other hand, in structural analysis, the finite element method and lath-beam structural model are introduced to calculate the elastic response. Furthermore, based on the current CFD and structural solver, a particular data transfer method and coupling strategy are presented for the fluid-structure interaction. Finally, by comparing numerical results with experimental data, the present method is validated to be available and feasible for hydroelastic water-entry problem and further successfully adopted to analyze the motion characteristics of free-falling elastic wedge.


---

## Body

## 1. Introduction

Interactions between fluid and elastic structure, especially water entry problem of elastic body, are of concern and practical importance in many ocean and marine engineering applications, such as seaplane landing on sea, ship slamming, sloshing, green-water impacting on deck, and marine lifeboat falling down. For the water entry problem, it is complicated and difficult to numerically study moving solid bodies and free surface and the challenge is further higher if global motion and elastic deformation of solid body need to be calculated from the fluid-structure interaction.The water entry problem has been widely studied by many scientists and engineers using various methods since the initial study of von Karman [1], such as Wagner [2], Yamamoto et al. [3], Worthington [4], Bottomley [5], Mayo [6], Stoffmacher [7], Chuang [8], Huges [9], Ochi and Motter [10], Belyschko and Mullen [11], Korobkin and Pukhnachov [12–15], Greenhow and Yanbao [16], Cointe and Armand [17], Wilson [18], Miloh [19], Zhao and Faltinsen [20, 21], Takagi [22], Fraenkel and McLeod [23], Wu et al. [24–26], Greenhow and Moyo [27], Sames et al. [28], Mei et al. [29], Lu et al. [30], Scolan and Korobkin [31], Aquelet and Souli [32–34], Battistin and Iafrati [35], Kleefsman et al. [36], Sun [37], Faltinsen and Semenov [38], Tveitnes et al. [39], Truscott and Techet [40], and Xu et al. [41].Therein boundary element method is based on the theory assumption of potential flow, so the rotational motion near free surface cannot be well simulated. Thus in order to better predict the force and motion of body and achieve fluid information near free surface simultaneously, a CFD method should be applied in water entry problem. Furthermore, compared with the traditional CFD method ([36, 42–44] and so on), the surface capturing method [45–48] based on Cartesian cut cell mesh [47–52] as a novel method can not only well simulate complicated cases of wave breaking and trapping, but also has higher computational efficiency to treat moving boundaries by updating a few cut cells locally rather than remeshing the whole flow domain. On the other hand, so far there are only few documents which consider the interaction among fluid field, global motion, and local deformation of elastic body simultaneously for the water entry problem. Therefore, a new CFD method of surface capturing method and Cartesian cut cell mesh should be taken and developed to handle hydroelastic water-entry problem of free-falling wedge.This paper presents surface capturing method and Cartesian cut cell mesh to treat moving free surface and solid boundaries. In the CFD solver, incompressible Euler equations are presented as governing equations for a variable density fluid and the location of free surface can be captured as a contact discontinuity in the density field. Furthermore, finite volume method is applied in numerical spatial discretization. Therein Roe’s approximate Riemann solver is adopted to evaluate numerical flow flux and dual-time stepping technique with artificial compressibility method is used for time advancing. On the other hand, in the structural solver finite element method and lath-beam structural model are used to analyze structural response of elastic wedge. Particularly based on the geometry of Cartesian cut cell mesh and lath-beam element here a particular data transfer method on solid boundary is deduced and presented. Furthermore, the flow field, structural response, and global motion are calculated in a coupled manner, which not only computes hydrodynamic force on solid boundary, but also allows the feedback of elastic deformation and global motion into the CFD solver. Finally, some test cases of water entry for various kinds of free-falling elastic wedge are numerically simulated and the calculated results show the feasibility and availability of the present method. Furthermore, by this method the global motion and local response of elastic wedge are studied and analyzed during water-entry phase.


## 2. Description of CFD Solver

### 2.1. Governing Equation and Boundary Condition

For the 2D, incompressible, unsteady, inviscid fluid system with a variable density field, thex- and y-coordinate axis are along and vertical to free surface, respectively; thus the governing equations (mass conservation, x, y-directional momentum, and incompressibility constraint equations) for surface capturing method can be written in conservation form, which have been verified and described by Kelecy and Pletcher [45], Pan and Chang [46], Qian et al. [47], and W. Wang and Y. Wang [48]:
(1)∂ρ∂t+∂(ρvx)∂x+∂(ρvy)∂y=0(2)∂(ρvx)∂t+∂(ρvxvx+p)∂x+∂(ρvxvy)∂y=ρBx(3)∂(ρvy)∂t+∂(ρvyvx)∂x+∂(ρvyvy+p)∂y=ρBy(4)∂vx∂x+∂vy∂y=0,
where, ρ is the fluid density, vx and vy represent the x and y-directional fluid velocities, p is the fluid pressure, Bx and By are the x and y-directional accelerations of body force, and, if assume the only body force is gravity, the source term Bx=0 and By=-ρg, respectively.In the present study, various kinds of boundary conditions can be classified.(1)
Outlet or open boundary: a zero gradient condition is applied in the velocity and density and the pressure at this boundary is fixed to be static pressure, which allows fluid to enter or leave the computational domain freely according to the local flow velocity and direction.(2)
Solid body boundary: the no-penetration condition can be applied in velocity and the density is assumed to have a zero normal gradient. For pressure boundary condition, here the momentum equations (2) and (3) are projected on the normal direction of solid boundary and can be rewritten as
(5)ρ∂(vsn)∂t+ρvf·∇vsn=-∂p∂n-ρgny,
where vf=vxi+vyj is fluid velocity at the solid boundary, n=nxi+nyj is normal unit vector of solid edge, and vsn is normal projection of average velocity of solid edge. Furthermore, a detailed description about how to achieve vsn and ∇vsn based on the fluid-structure interaction on solid boundary should be demonstrated in Section 4.1.2 of data transfer method.


### 2.2. Numerical Solver

In the present study, Cartesian cut cell mesh is employed for spatial discretization, which can be generated by cutting solid bodies out of a background Cartesian mesh. Thus fluid cell, solid cell and cut cell are accordingly created, which include various storage data (such as identification of cell type, geometry, and fluid variables). For more information, refer to Coirier and Powell [49], Yang et al. [50, 51], Causon et al. [52], and Qian et al. [47].Based on the Cartesian cut cell mesh, here a cell-central finite volume method is applied in numerical discretization. In this numerical scheme, Roe’s approximate Riemann solver is adopted to calculate the numerical flux on each edge of fluid cell where fluid variables are reconstructed by using a piecewise linear upwind scheme. In the linear representation, a least-square method [48] is taken to achieve variable gradients and Superbee limiter [47] is used to maintain monotony of numerical scheme and control spurious oscillations. Then the flux on solid boundary can be achieved by means of exact Riemann discontinuous solution [48]. Furthermore, dual-time stepping technique with artificial compressibility method is used for time advancing. Finally, an approximate LU factorization (ALU) scheme [53] is applied to solve linear equations. The details about numerical scheme can be further found from the study of W. Wang and Y. Wang [48].


## 2.1. Governing Equation and Boundary Condition

For the 2D, incompressible, unsteady, inviscid fluid system with a variable density field, thex- and y-coordinate axis are along and vertical to free surface, respectively; thus the governing equations (mass conservation, x, y-directional momentum, and incompressibility constraint equations) for surface capturing method can be written in conservation form, which have been verified and described by Kelecy and Pletcher [45], Pan and Chang [46], Qian et al. [47], and W. Wang and Y. Wang [48]:
(1)∂ρ∂t+∂(ρvx)∂x+∂(ρvy)∂y=0(2)∂(ρvx)∂t+∂(ρvxvx+p)∂x+∂(ρvxvy)∂y=ρBx(3)∂(ρvy)∂t+∂(ρvyvx)∂x+∂(ρvyvy+p)∂y=ρBy(4)∂vx∂x+∂vy∂y=0,
where, ρ is the fluid density, vx and vy represent the x and y-directional fluid velocities, p is the fluid pressure, Bx and By are the x and y-directional accelerations of body force, and, if assume the only body force is gravity, the source term Bx=0 and By=-ρg, respectively.In the present study, various kinds of boundary conditions can be classified.(1)
Outlet or open boundary: a zero gradient condition is applied in the velocity and density and the pressure at this boundary is fixed to be static pressure, which allows fluid to enter or leave the computational domain freely according to the local flow velocity and direction.(2)
Solid body boundary: the no-penetration condition can be applied in velocity and the density is assumed to have a zero normal gradient. For pressure boundary condition, here the momentum equations (2) and (3) are projected on the normal direction of solid boundary and can be rewritten as
(5)ρ∂(vsn)∂t+ρvf·∇vsn=-∂p∂n-ρgny,
where vf=vxi+vyj is fluid velocity at the solid boundary, n=nxi+nyj is normal unit vector of solid edge, and vsn is normal projection of average velocity of solid edge. Furthermore, a detailed description about how to achieve vsn and ∇vsn based on the fluid-structure interaction on solid boundary should be demonstrated in Section 4.1.2 of data transfer method.


## 2.2. Numerical Solver

In the present study, Cartesian cut cell mesh is employed for spatial discretization, which can be generated by cutting solid bodies out of a background Cartesian mesh. Thus fluid cell, solid cell and cut cell are accordingly created, which include various storage data (such as identification of cell type, geometry, and fluid variables). For more information, refer to Coirier and Powell [49], Yang et al. [50, 51], Causon et al. [52], and Qian et al. [47].Based on the Cartesian cut cell mesh, here a cell-central finite volume method is applied in numerical discretization. In this numerical scheme, Roe’s approximate Riemann solver is adopted to calculate the numerical flux on each edge of fluid cell where fluid variables are reconstructed by using a piecewise linear upwind scheme. In the linear representation, a least-square method [48] is taken to achieve variable gradients and Superbee limiter [47] is used to maintain monotony of numerical scheme and control spurious oscillations. Then the flux on solid boundary can be achieved by means of exact Riemann discontinuous solution [48]. Furthermore, dual-time stepping technique with artificial compressibility method is used for time advancing. Finally, an approximate LU factorization (ALU) scheme [53] is applied to solve linear equations. The details about numerical scheme can be further found from the study of W. Wang and Y. Wang [48].


## 3. Introduction of Lath-Beam Structural Model

For the bottom plate with cylindrical bending of elastic wedge structure, a lath-beam model can be used to analyze structural deformation. Compared with Euler-Bernoulli beam, due to bilateral constraints, the cross-section of lath-beam model cannot freely deform and keep the original shape. On this basis, only if by introducingE1=E/(1-μ2) (E is the material elastic modulus and μ is Poisson ratio) as new elastic modulus into Euler-Bernoulli beam theory can the elastic deformation of lath-beam model be studied. For more information, please refer to Shu and Tan [54]. Thus, in structural solver, here finite element method and lath-beam model based on Euler-Bernoulli beam theory should be applied to calculate structural response on solid boundaries of free-falling elastic wedge.Initially, the global (fixed) and local (moving with body) coordinate systems are created. Then, lath-beam structure is divided into many beam elements with equal length, which can support bending and tension-compression effects simultaneously. Here in the local coordinate system, the nodal displacement vector{δ}k of element k can be demonstrated as
(6){δ}k=[u1w1θ1,u2w2θ2]k,
where u1, w1, u2, and w2 are translational displacements, and θ1 and θ2 are rotation angles on two nodes.Next, in local coordinate system the corresponding mass and stiffness matrix of lath-beam elements are generated and then load vectors of element nodes are calculated from the fluid pressure on solid boundary. Subsequently, the mass matrix, stiffness matrix, and load vector are transformed into global coordinate system. For the water-entry model of elastic wedge, the elastic deformations basically take place on sloping sides; therefore here on the corner of sloping and straight sides the rigid connection is selected as boundary constraint of lath-beam structure. Then the global dynamical equation without structural damping can be written as(7)[M]{δ¨}+[K]{δ}={F},
where {F} is the nodal load vector and the detail calculation method can be referred to in Section 4.1.1. In matrix [M] and [K], for the structural element k, the [M]k and [K]k can be expressed as
(8)[M]k=ρshl420[140007000015622l054-13l022l4l2013l-3l270001400005413l0156-22l0-13l-3l20-22l4l2],[K]k=E1l[A00-A00012Il26Il0-12Il26Il06Il4I0-6Il2I-A00A000-12Il2-6Il012Il2-6Il06Il2I0-6Il4I],
where, for lath-beam structural model, ρs is the structural density, h is the element thickness, l is the element length, A=1.0*h is the cross-section area, I=1.0*h3/12 is the moment of inertia, E1=E/(1-μ2), μ is the Poisson ratio, and E is the material elastic modulus.Finally, by means of numerical Newmark method, the nodal vector of displacement, velocity, and acceleration are computed and fed back to CFD solver and meanwhile local structural response on elastic boundary can be obtained.


## 4. Analysis of Fluid-Structure Interaction

For the water-entry model of free-falling elastic wedge, the global motion and local deformation of elastic body should be calculated from the fluid-structure interaction, so here the data transfer method on elastic boundary and fluid-structure coupling strategy are discussed, respectively.


### 4.1. Data Transfer Method on Elastic Boundary

In order to study and analyze the fluid-structure interaction, it is necessary to transfer the fluid pressure to the element node and then feed the global motion and structural deformation back to CFD solver.In this paper, Cartesian cut cell mesh in CFD solver and lath-beam element in structural solver are different from each other. So here a particular data transfer method on solid boundary is deduced and demonstrated, which includes transfer of distributional fluid pressure to structural nodes and feedback of structural response on CFD boundary condition.


#### 4.1.1. Transfer of Distributional Fluid Pressure to Structural Nodes

Based on respective geometries of Cartesian cut cell mesh and lath-beam element, there should be many kinds of position relations between the solid boundaryAB in cut cell and structural element L1L2 to be considered (three classic cases as shown in Figure 1), which are full inclusion of AB in L1L2 (Figure 1(a)) or L1L2 in AB (Figure 1(b)), partial intersection between L1L2 with AB (Figure 1(c)), and some other situations (such as that AB coincides with L1L2). Nevertheless, the principle of transferring the pressure distribution to the FEM nodes is the same in all cases. Thus, the important issue is how to estimate the pressure distribution along the beam elements. Here take the structural element L1L2 in Figure 1(a), for example, to demonstrate how to calculate the nodal forces F1 and F2 of the nodes L1 and L2 from fluid pressure distributed on AB.Different position relations of cut cell and structural elements; (a) full inclusion ofAB in L1L2; (b) full inclusion of L1L2 in AB; (c) partial intersection between L1L2 and AB.
(a)
(b)
(c)Firstly, local coordinate systemx′o′y′ should be created on structural L1L2. Therein coordinateorigin o′ coincides with node L1, coordinates axis o′x′ is clockwise and, along with element L1L2, axis o′y′ points from fluid to structure. Thus unit vectors of o′x′ and o′y′ can be written:
(9)o′x′=x2-x1(x2-x1)2+(y2-y1)2i+y2-y1(x2-x1)2+(y2-y1)2j(10)o′y′=-y2-y1(x2-x1)2+(y2-y1)2i+x2-x1(x2-x1)2+(y2-y1)2j,
where (x1,y1) and (x2,y2) are coordinates of nodes L1 and L2, respectively.Next, by using the coordinates(xo,yo) of centre O in cut cell, position vector OL1 can be achieved:
(11)OL1=(x1-xo)i+(y1-yo)j.Then by combining (9) with (11), position vector Δr (from centre O to arbitrary point C on AB) can be written:
(12)Δr=OL1+L1C=OL1+x′·o′x′=[(x1-xo)+(x2-x1)·x′(x2-x1)2+(y2-y1)2]i+[(y1-yo)+(y2-y1)·x′(x2-x1)2+(y2-y1)2]j,
where x′=|L1C| is the coordinate value along coordinates axis o′x′.Finally, according to the equivalence principle about total force and moment, the distributional fluid pressure onAB can be transferred to the equivalent concentrated forces F1 and F2 on nodes L1 and L2 of structural lath-beam element L1L2, which can be described as follows:
(13)F1=∫L1AL1B[(po+Δp·Δr)·(L1L2-x′)]dx′L1L2F2=∫L1AL1B[(po+Δp·Δr)·x′]dx′L1L2,
where Po and Δp=pxi+pyj are the fluid pressure and its gradient stored in centre O of cut cell, which can be obtained by CFD solver. Similarly, the equivalent concentrated forces on element nodes can be calculated for the other situations of position relations.


#### 4.1.2. Feedback of Structural Response on CFD Boundary Condition

Here among many kinds of position relations of Cartesian cut cell mesh and lath-beam element as mentioned in Section4.1.1 (Figure 1), take one (as shown in Figure 1(a)), for example, to describe how to adopt structural response to update boundary condition of CFD solver.In Figure1(a), structural element L1L2 fully includes the solid boundary AB. Here the velocities of nodes L1(u˙1,w˙1,θ˙1) and L2(u˙2,w˙2,θ˙2) in local coordinate system x′o′y′ have been calculated by finite element analysis.Firstly based on the shape function of structural elementL1L2, the velocities of point A in local coordinate system can be obtained:
(14)u˙A=N5u˙1+N6u˙2w˙A=N1w˙1+N2θ˙1+N3w˙2+N4θ˙2,
where N1~N4 and N5~N6 are the shape function versus bending and tension-compression effects, respectively. For the lath-beam element, N1(ξ1)=1-3ξ12+2ξ13,  N2(ξ1)=(ξ1-2ξ12+ξ13)L1L2,  N3(ξ1)=3ξ12-2ξ13,  N4(ξ1)=(ξ13-ξ12)L1L2,  N5(ξ2)=(1-ξ2)/2, and N6(ξ2)=(1+ξ2)/2, ξ1=L1A/L1L2 and ξ2=2(L1A-L1L2/2)/L1L2 are natural coordinates of structural element L1L2. Then by the same method, the velocities u˙B and w˙B of point B also can be calculated.Here it should be noted thatvsn and ∇vsn terms in (5) are the sum velocity, which includes local elastic deformation and global rigid translation. Subsequently, the average normal velocity vsn and tangential velocity vst of solid boundary AB can be calculated:
(15)vsn=w˙A+w˙B2+vr·n,(16)vst=u˙A+u˙B2+vr·t,
where n and t are normal and tangential unit vectors of solid boundary and vr=vrxi+vryj is the global rigid velocity of solid boundary, which can be calculated by a subiteration coupling approach in Section 4.2.Meanwhile, gradient of normal velocity∇vsn=(∂vsn/∂x)i+(∂vsn/∂y)j on solid boundary AB can be achieved:
(17)∂vsn∂x=w˙B-w˙AxB-xA,(18)∂vsn∂y=w˙B-w˙AyB-yA.Finally, the tangential velocity, average normal velocity and its gradient can be obtained from (15)~(18), and then should be applied in the new boundary condition (5) of cut cell.


### 4.2. Fluid-Structure Coupled Solution Strategy

In the present study, the fluid-structure interaction of water-entry model for free-falling elastic wedge includes the coupling calculation of fluid with global body motion and local structural deformation.On the one hand, in order to ensure enough numerical stability, a subiteration coupling approach for the global velocity of wedge proposed by Kleefsman et al. [36] is taken. Here assume that the elastic wedge moves with one degree (up and down) with an acceleration; the x-direction velocity is zero and the y-directional global rigid velocity of solid boundary can be computed by
(19)(vryn+1)k+1=ω(vryn+ΔtFbodykm)+(1-ω)(vryn+1)k,
where Δt is the physical time step, m is the wedge mass, and vryn is the y-directional global rigid velocity of solid edge 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. Equation (19) demonstrates a subiterative process of fluid-global body motion coupling at each physical moment. At physical time step n, the iterative initial value (vryn+1)0 is taken as vryn. For the subiterative k step, based on the vertical body velocity (vryn+1)k, a new pressure field is calculated and then Fbodyk is obtained by direct integration of the new pressure over body surface. When dV=(vryn+1)k+1-(vryn+1)k<10-3, the subiteration is convergence and (vryn+1)k+1 is taken as the new velocity boundary condition vryn+1 of the solid body for the next physical time step n+1.On the other hand, compared with military missile at high speed, ship and ocean engineering structures usually enter water with medium-low velocity and therefore the deformation of elastic boundary is relatively small and its reaction to the fluid field can be neglected. Thus in order to increase the computational efficiency, here a two-way “weak” coupling method as shown in Figure2 is used for the coupled solution of fluid and local structural deformation.Figure 2
Flow chart of iteration procedure for the fluid-structure coupled solution strategy.


## 4.1. Data Transfer Method on Elastic Boundary

In order to study and analyze the fluid-structure interaction, it is necessary to transfer the fluid pressure to the element node and then feed the global motion and structural deformation back to CFD solver.In this paper, Cartesian cut cell mesh in CFD solver and lath-beam element in structural solver are different from each other. So here a particular data transfer method on solid boundary is deduced and demonstrated, which includes transfer of distributional fluid pressure to structural nodes and feedback of structural response on CFD boundary condition.


### 4.1.1. Transfer of Distributional Fluid Pressure to Structural Nodes

Based on respective geometries of Cartesian cut cell mesh and lath-beam element, there should be many kinds of position relations between the solid boundaryAB in cut cell and structural element L1L2 to be considered (three classic cases as shown in Figure 1), which are full inclusion of AB in L1L2 (Figure 1(a)) or L1L2 in AB (Figure 1(b)), partial intersection between L1L2 with AB (Figure 1(c)), and some other situations (such as that AB coincides with L1L2). Nevertheless, the principle of transferring the pressure distribution to the FEM nodes is the same in all cases. Thus, the important issue is how to estimate the pressure distribution along the beam elements. Here take the structural element L1L2 in Figure 1(a), for example, to demonstrate how to calculate the nodal forces F1 and F2 of the nodes L1 and L2 from fluid pressure distributed on AB.Different position relations of cut cell and structural elements; (a) full inclusion ofAB in L1L2; (b) full inclusion of L1L2 in AB; (c) partial intersection between L1L2 and AB.
(a)
(b)
(c)Firstly, local coordinate systemx′o′y′ should be created on structural L1L2. Therein coordinateorigin o′ coincides with node L1, coordinates axis o′x′ is clockwise and, along with element L1L2, axis o′y′ points from fluid to structure. Thus unit vectors of o′x′ and o′y′ can be written:
(9)o′x′=x2-x1(x2-x1)2+(y2-y1)2i+y2-y1(x2-x1)2+(y2-y1)2j(10)o′y′=-y2-y1(x2-x1)2+(y2-y1)2i+x2-x1(x2-x1)2+(y2-y1)2j,
where (x1,y1) and (x2,y2) are coordinates of nodes L1 and L2, respectively.Next, by using the coordinates(xo,yo) of centre O in cut cell, position vector OL1 can be achieved:
(11)OL1=(x1-xo)i+(y1-yo)j.Then by combining (9) with (11), position vector Δr (from centre O to arbitrary point C on AB) can be written:
(12)Δr=OL1+L1C=OL1+x′·o′x′=[(x1-xo)+(x2-x1)·x′(x2-x1)2+(y2-y1)2]i+[(y1-yo)+(y2-y1)·x′(x2-x1)2+(y2-y1)2]j,
where x′=|L1C| is the coordinate value along coordinates axis o′x′.Finally, according to the equivalence principle about total force and moment, the distributional fluid pressure onAB can be transferred to the equivalent concentrated forces F1 and F2 on nodes L1 and L2 of structural lath-beam element L1L2, which can be described as follows:
(13)F1=∫L1AL1B[(po+Δp·Δr)·(L1L2-x′)]dx′L1L2F2=∫L1AL1B[(po+Δp·Δr)·x′]dx′L1L2,
where Po and Δp=pxi+pyj are the fluid pressure and its gradient stored in centre O of cut cell, which can be obtained by CFD solver. Similarly, the equivalent concentrated forces on element nodes can be calculated for the other situations of position relations.


### 4.1.2. Feedback of Structural Response on CFD Boundary Condition

Here among many kinds of position relations of Cartesian cut cell mesh and lath-beam element as mentioned in Section4.1.1 (Figure 1), take one (as shown in Figure 1(a)), for example, to describe how to adopt structural response to update boundary condition of CFD solver.In Figure1(a), structural element L1L2 fully includes the solid boundary AB. Here the velocities of nodes L1(u˙1,w˙1,θ˙1) and L2(u˙2,w˙2,θ˙2) in local coordinate system x′o′y′ have been calculated by finite element analysis.Firstly based on the shape function of structural elementL1L2, the velocities of point A in local coordinate system can be obtained:
(14)u˙A=N5u˙1+N6u˙2w˙A=N1w˙1+N2θ˙1+N3w˙2+N4θ˙2,
where N1~N4 and N5~N6 are the shape function versus bending and tension-compression effects, respectively. For the lath-beam element, N1(ξ1)=1-3ξ12+2ξ13,  N2(ξ1)=(ξ1-2ξ12+ξ13)L1L2,  N3(ξ1)=3ξ12-2ξ13,  N4(ξ1)=(ξ13-ξ12)L1L2,  N5(ξ2)=(1-ξ2)/2, and N6(ξ2)=(1+ξ2)/2, ξ1=L1A/L1L2 and ξ2=2(L1A-L1L2/2)/L1L2 are natural coordinates of structural element L1L2. Then by the same method, the velocities u˙B and w˙B of point B also can be calculated.Here it should be noted thatvsn and ∇vsn terms in (5) are the sum velocity, which includes local elastic deformation and global rigid translation. Subsequently, the average normal velocity vsn and tangential velocity vst of solid boundary AB can be calculated:
(15)vsn=w˙A+w˙B2+vr·n,(16)vst=u˙A+u˙B2+vr·t,
where n and t are normal and tangential unit vectors of solid boundary and vr=vrxi+vryj is the global rigid velocity of solid boundary, which can be calculated by a subiteration coupling approach in Section 4.2.Meanwhile, gradient of normal velocity∇vsn=(∂vsn/∂x)i+(∂vsn/∂y)j on solid boundary AB can be achieved:
(17)∂vsn∂x=w˙B-w˙AxB-xA,(18)∂vsn∂y=w˙B-w˙AyB-yA.Finally, the tangential velocity, average normal velocity and its gradient can be obtained from (15)~(18), and then should be applied in the new boundary condition (5) of cut cell.


## 4.1.1. Transfer of Distributional Fluid Pressure to Structural Nodes

Based on respective geometries of Cartesian cut cell mesh and lath-beam element, there should be many kinds of position relations between the solid boundaryAB in cut cell and structural element L1L2 to be considered (three classic cases as shown in Figure 1), which are full inclusion of AB in L1L2 (Figure 1(a)) or L1L2 in AB (Figure 1(b)), partial intersection between L1L2 with AB (Figure 1(c)), and some other situations (such as that AB coincides with L1L2). Nevertheless, the principle of transferring the pressure distribution to the FEM nodes is the same in all cases. Thus, the important issue is how to estimate the pressure distribution along the beam elements. Here take the structural element L1L2 in Figure 1(a), for example, to demonstrate how to calculate the nodal forces F1 and F2 of the nodes L1 and L2 from fluid pressure distributed on AB.Different position relations of cut cell and structural elements; (a) full inclusion ofAB in L1L2; (b) full inclusion of L1L2 in AB; (c) partial intersection between L1L2 and AB.
(a)
(b)
(c)Firstly, local coordinate systemx′o′y′ should be created on structural L1L2. Therein coordinateorigin o′ coincides with node L1, coordinates axis o′x′ is clockwise and, along with element L1L2, axis o′y′ points from fluid to structure. Thus unit vectors of o′x′ and o′y′ can be written:
(9)o′x′=x2-x1(x2-x1)2+(y2-y1)2i+y2-y1(x2-x1)2+(y2-y1)2j(10)o′y′=-y2-y1(x2-x1)2+(y2-y1)2i+x2-x1(x2-x1)2+(y2-y1)2j,
where (x1,y1) and (x2,y2) are coordinates of nodes L1 and L2, respectively.Next, by using the coordinates(xo,yo) of centre O in cut cell, position vector OL1 can be achieved:
(11)OL1=(x1-xo)i+(y1-yo)j.Then by combining (9) with (11), position vector Δr (from centre O to arbitrary point C on AB) can be written:
(12)Δr=OL1+L1C=OL1+x′·o′x′=[(x1-xo)+(x2-x1)·x′(x2-x1)2+(y2-y1)2]i+[(y1-yo)+(y2-y1)·x′(x2-x1)2+(y2-y1)2]j,
where x′=|L1C| is the coordinate value along coordinates axis o′x′.Finally, according to the equivalence principle about total force and moment, the distributional fluid pressure onAB can be transferred to the equivalent concentrated forces F1 and F2 on nodes L1 and L2 of structural lath-beam element L1L2, which can be described as follows:
(13)F1=∫L1AL1B[(po+Δp·Δr)·(L1L2-x′)]dx′L1L2F2=∫L1AL1B[(po+Δp·Δr)·x′]dx′L1L2,
where Po and Δp=pxi+pyj are the fluid pressure and its gradient stored in centre O of cut cell, which can be obtained by CFD solver. Similarly, the equivalent concentrated forces on element nodes can be calculated for the other situations of position relations.


## 4.1.2. Feedback of Structural Response on CFD Boundary Condition

Here among many kinds of position relations of Cartesian cut cell mesh and lath-beam element as mentioned in Section4.1.1 (Figure 1), take one (as shown in Figure 1(a)), for example, to describe how to adopt structural response to update boundary condition of CFD solver.In Figure1(a), structural element L1L2 fully includes the solid boundary AB. Here the velocities of nodes L1(u˙1,w˙1,θ˙1) and L2(u˙2,w˙2,θ˙2) in local coordinate system x′o′y′ have been calculated by finite element analysis.Firstly based on the shape function of structural elementL1L2, the velocities of point A in local coordinate system can be obtained:
(14)u˙A=N5u˙1+N6u˙2w˙A=N1w˙1+N2θ˙1+N3w˙2+N4θ˙2,
where N1~N4 and N5~N6 are the shape function versus bending and tension-compression effects, respectively. For the lath-beam element, N1(ξ1)=1-3ξ12+2ξ13,  N2(ξ1)=(ξ1-2ξ12+ξ13)L1L2,  N3(ξ1)=3ξ12-2ξ13,  N4(ξ1)=(ξ13-ξ12)L1L2,  N5(ξ2)=(1-ξ2)/2, and N6(ξ2)=(1+ξ2)/2, ξ1=L1A/L1L2 and ξ2=2(L1A-L1L2/2)/L1L2 are natural coordinates of structural element L1L2. Then by the same method, the velocities u˙B and w˙B of point B also can be calculated.Here it should be noted thatvsn and ∇vsn terms in (5) are the sum velocity, which includes local elastic deformation and global rigid translation. Subsequently, the average normal velocity vsn and tangential velocity vst of solid boundary AB can be calculated:
(15)vsn=w˙A+w˙B2+vr·n,(16)vst=u˙A+u˙B2+vr·t,
where n and t are normal and tangential unit vectors of solid boundary and vr=vrxi+vryj is the global rigid velocity of solid boundary, which can be calculated by a subiteration coupling approach in Section 4.2.Meanwhile, gradient of normal velocity∇vsn=(∂vsn/∂x)i+(∂vsn/∂y)j on solid boundary AB can be achieved:
(17)∂vsn∂x=w˙B-w˙AxB-xA,(18)∂vsn∂y=w˙B-w˙AyB-yA.Finally, the tangential velocity, average normal velocity and its gradient can be obtained from (15)~(18), and then should be applied in the new boundary condition (5) of cut cell.


## 4.2. Fluid-Structure Coupled Solution Strategy

In the present study, the fluid-structure interaction of water-entry model for free-falling elastic wedge includes the coupling calculation of fluid with global body motion and local structural deformation.On the one hand, in order to ensure enough numerical stability, a subiteration coupling approach for the global velocity of wedge proposed by Kleefsman et al. [36] is taken. Here assume that the elastic wedge moves with one degree (up and down) with an acceleration; the x-direction velocity is zero and the y-directional global rigid velocity of solid boundary can be computed by
(19)(vryn+1)k+1=ω(vryn+ΔtFbodykm)+(1-ω)(vryn+1)k,
where Δt is the physical time step, m is the wedge mass, and vryn is the y-directional global rigid velocity of solid edge 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. Equation (19) demonstrates a subiterative process of fluid-global body motion coupling at each physical moment. At physical time step n, the iterative initial value (vryn+1)0 is taken as vryn. For the subiterative k step, based on the vertical body velocity (vryn+1)k, a new pressure field is calculated and then Fbodyk is obtained by direct integration of the new pressure over body surface. When dV=(vryn+1)k+1-(vryn+1)k<10-3, the subiteration is convergence and (vryn+1)k+1 is taken as the new velocity boundary condition vryn+1 of the solid body for the next physical time step n+1.On the other hand, compared with military missile at high speed, ship and ocean engineering structures usually enter water with medium-low velocity and therefore the deformation of elastic boundary is relatively small and its reaction to the fluid field can be neglected. Thus in order to increase the computational efficiency, here a two-way “weak” coupling method as shown in Figure2 is used for the coupled solution of fluid and local structural deformation.Figure 2
Flow chart of iteration procedure for the fluid-structure coupled solution strategy.


## 5. Results of Test Cases and Discussions

### 5.1. Validation

In this section, the current numerical method for hydroelastic analysis on water-entry problems of free-falling wedge should be tested by comparing calculated results with experimental data of Sun et al. [55] about three different water-entry test cases.Based on the parameters of water tank in Sun’s experiment, 2D CFD computational domain was a square of 0.8 m × 1.45 m and the water depth was 1.1 m. The material properties of elastic wedges were 1.059 × 103 kg/m3 density, 2.67 × 109 kg/(ms2) elastic module, and 0.357 passion ratio. Furthermore, geometric parameters of three different wedges were in Table 1 and the diagrammatic sketch is shown in Figure 3, where L represents length of elastic wedge, B is width, and α is deadrise angle.Table 1
Geometric parameters of three different elastic wedges.
Number
Parameter
Unit
Value
1
Length
m
0.6
Width
m
0.2
Thickness
mm
1.5
Mass
kg
7.7870
Deadrise angle
45°
2
Length
m
0.6
Width
m
0.2
Thickness
mm
1.5
Mass
kg
7.4780
Deadrise angle
20°
3
Length
m
0.6
Width
m
0.2
Thickness
mm
3.0
Mass
kg
7.7890
Deadrise angle
20°Figure 3
Sketch of elastic model and description of geometric parameters.In the numerical calculation, the zero physical time (t=0) was selected as the moment in which the bottom of the wedge is 0.1 m up from the initial free surface and then the elastic wedge began to fall freely. The physical time Δt was 0.0001 s, fictitious time Δτ was 0.01 s, the artificial compressibility coefficient β was 500, the gravity acceleration was 9.81 m/s2, and the length of structural element on solid boundary was 0.01 m. Furthermore, the inhomogeneous meshes were generated in fluid field. Near solid boundary the local meshes with equal spacing Δx=Δy=0.01 m were refined and regular and then the mesh size gradually enlarged away from boundary. By taking the model 1, for example, an inhomogeneous mesh system 66 × 90 is shown in Figure 4.Global mesh (a) and local grid (b) near solid boundary.
(a)
(b)In Figure5, the time histories of global acceleration and local strain of observation point on the sloping edge versus numerical solution and experiment data for three different elastic wedges are compared. Therein Figures 5(a), 5(b), and 5(c) represent various wedge models (number 1~3 in Table 1), respectively. From the figure, it can be found that the numerical solutions by this paper and experimental data are in good agreement with each other, except acceleration peak and negative strain of elastic wedge. On the one hand, for the sake of resistance of experimental equipment, wedge accelerations in experiment before water entry are less than the numerical simulation, which should cause a little difference about acceleration peak. On the other hand, by means of numerical method, the absolute value of negative strain on solid boundary is a bit underpredicted. The reason may be that the two-way “weak” coupling strategy between fluid and local structure artificially increases the numerical stiffness of water-entry model.Time history of global acceleration and local strain of observation point on the sloping edge for three different elastic wedges, (a) 1.5 mm thickness, 7.787 kg mass, and 45° deadrise angle; (b) 1.5 mm thickness, 7.478 kg mass, and 20° deadrise angle; (c) 3.0 mm thickness, 7.789 kg mass, and 20° deadrise angle.
(a)
(b)
(c)


### 5.2. Description of Free Surface Profile and Velocity Vectors

By using the present method to study water-entry model 1 of Table1, the developing process of free surface shape and velocity vectors in flow field can be successfully numerically simulated and the results at various water-entry moments (t=0.00 s, 0.12 s, 0,18 s, 0.24 s, 0.30 s, 0.42 s, 0.48 s, and 0.60 s) are shown in Figure 6 to describe the water-entry phenomenon of the free-falling elastic wedge.Free surface profile and velocity vectors for a free-falling elastic wedge with 1.5 mm thickness, 7.787 kg mass, and 45° deadrise angle; (a)t=0.00 s, (b) 0.12 s, (c) 0.18 s, (d) 0.24 s, (e) 0.30 s, (f) 0.42 s, (g) 0.48 s, and (h) 0.60 s.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)Att=0.00 s the first graph shows the initial location of elastic wedge and calm free surface. Then the elastic wedge freely falls down (t=0.12 s) and the air near solid surface moves with the wedge. As the body pieces the free surface (t=0.18 s), the water rises up along each sloping side and the displacement waves are formed originating at the vertex of wedge. Meanwhile, the air between the solid boundary and free surface is forced around wedge edge and separates at the corner to form a pair of symmetric vortices behind the body in air. As the body travels further into the water (t=0.24 s and 0.30 s), the two intersections between wedge and free surface move up along both straight sides and subsequently overturn towards each other (t=0.42 s). After the free surface is reconnected (t=0.48 s), a jet flow is created upwards above the wedge (t=0.60 s).


### 5.3. Analysis on Hydroelasticity of Free-Falling Wedge

In order to study the global dynamic behaviors and local structural responses of elastic wedge, here, based on the geometric parameters of model 1 in Table1, two kinds of water-entry models versus rigid and elastic wedge are created simultaneously. By comparing with the numerical results of different models, the hydroelasticity of free-falling wedge during the water-entry phase is studied and analyzed.


#### 5.3.1. Global Dynamic Behaviors of Elastic Wedge

Firstly, by means of the current numerical method with or without elastic deformation, the global hydrodynamic acceleration and velocity versus elastic and rigid wedge model were calculated and compared as shown in Figure7.Figure 7
Time history of global hydrodynamic acceleration and velocity versus elastic and rigid wedge model with 7.787 kg mass and 45° deadrise angle.In Figure7, for the sake of structural deformation, some high-frequency vibrations exist in the variation curves of global hydrodynamic and acceleration for the elastic model. In order to directly and clearly study the dynamic behaviors of elastic wedge, here a Savitzky-Golay method with 5 points and 4 degrees was used to filter high-frequency oscillation and smooth variable curves. For the rigid and elastic model, the time history of hydrodynamic and acceleration are in substantial agreement except two little differences. One is the peak value of hydrodynamic and acceleration curves and the results of rigid model are a little larger than the elastic model. The reason may be that the elastic deformation can provide buffer action and weaken fluid hydrodynamic during the water-entry phase. The other is that there are many high-frequency vibrations in hydrodynamic and acceleration curves, which can be caused by local deformation of elastic structure. However, the amplitudes of these vibrations are too small to influence the global motion of wedge. Thus, the velocity curves versus two models are basically in agreement with each other.On this basis, in order to further study the long global motion of free-falling wedge during the water entry phase, here the water entry model was numerically simulated for a long physical time and the results of wedge velocity were shown in the last picture of Figure7. From the picture, the motion process of free-falling wedge during water entry phase can be described as follows: firstly the wedge in air free falls straight down and then slows down due to large hydrodynamic at the initial moment of water-entry. With the continuous hydrodynamic acting on the wedge, the velocity gradually decreases until zero. Next, because the buoyancy of wedge is larger than its gravity, along with the closure of free surface, the wedge begins to move upwards. Subsequently, the jet successively grows upwards, falls downwards, and finally disappears on free surface, which will influence the global motion of wedge. Furthermore, by combining the effect of fluid force and gravity, wedge should vibrate up and down with less and less amplitude.Figure8 shows the pressure contour in fluid field at classic water-entry moments (t=0.15 s, 0.18 s, 0.21 s, and 0.24 s). In the figure, at the initial stage of water-entry (t=0.15 s), wedge has just begun to pierce free surface and the pressure peak appears near the vertex of wedge. Furthermore, the air escapes with high velocity from the gap of free surface and solid boundary to both sides and hence forms a strong air-flow to press on the wedge edges. With the wedge continuously entering water (t=0.18 s), free surface rises up along each sloping side and the amplitude and influence sphere of fluid pressure on wedge edge gradually enlarge. Furthermore, it can be found that the variation rate of pressure distribution on the wedge edge is related to the location of free surface. Therein the pressure slowly changes in water and rapidly reduces in air. As the body travels further into the water (t=0.18 s and 0.24 s), the two intersections between wedge edge and free surface move up through corners and along both straight sides. Because of the continuous effect of fluid hydrodynamic, the wedge velocity becomes smaller and smaller, which results in the reduction of pressure distribution on wedge edge. Furthermore, it should be noted that, because near the unsmooth corners the fluid velocity become very large, the negative pressure on solid boundary should take place over there.Figure 8
Pressure contour in fluid field for elastic wedge with 7.787 kg mass and 45° deadrise angle;t=0.15 s, 0.18 s, 0.21 s, and 0.24 s (from left to right and from up to bottom).


#### 5.3.2. Local Structural Response on Elastic Edge

The time history and frequency analysis of local pressure and strain on the midpoint of wedge edge are shown in Figure9. In order to better discuss the characteristic of local structural deformation; here fast Fourier transform method is applied in time-frequency analysis and the local pressure and strain in frequency domain are achieved and shown. From the frequency analysis, it can be found that the variations of local pressure and stain on the midpoint mainly include two characteristics.Figure 9
Time history and frequency analysis of local pressure and strain on the midpoint of wedge edge versus elastic and rigid model with 7.787 kg mass and 45° deadrise angle.One is the pressure peak with high amplitude and narrow range in low frequency domain. Based on the effect of max pressure (frequency tends to zero) in Figure9 and global hydrodynamic in Figure 7, elastic boundaries should be deformed with 1st mode of vibration. Furthermore, the corresponding frequency of max strain is 1.3021 Hz and approximately equal to the first order natural frequency (1.2241 Hz) of the elastic model which is calculated by ANSYS software.The other is the pressure oscillation with approximately 72~85 Hz range and low amplitude in high frequency domain which includes two crests and one trough. From the modal analysis, the pressure oscillation possibly causes the midpoint to locally vibrate with 11th and 12th natural frequencies 73.13 and 86.822 Hz which are calculated by ANSYS software. Furthermore, because the two natural frequencies are very close, the vibration phenomenon of “beats” should appear on the midpoint of wedge edge. Here according to the properties of “beats” phenomenon, the theoretical frequencies versus amplitude profile(ω1-ω2)/4π and high-frequency vibration (ω1+ω2)/4π are 6.846 and 79.976 Hz, respectively. Correspondingly, from the variation of local strains in frequency and time domain as shown in Figure 9, it can be found that the calculated frequencies (7.33 and 75.52 Hz) are close to the theoretical values, which theoretically and numerically verifies the “beats” phenomenon of elastic structural vibration.Thus the local strain of midpoint on wedge edge can be divided into main deformation versus 1st principal mode and local deformation of “beat” phenomenon versus 11th and 12th modes.In order to evaluate the hydroelastic effect of free-falling wedge, here case 1 is numerically simulated without feedback loop of elastic deformation, which means that the local structural response is handled as just a “pure” postprocessing step. Then the results are shown and compared with two-way “weak” coupling method in Figures10 and 11.Pressure distribution on elastic boundaries of a free-falling elastic wedge with 1.5 mm thickness, 7.787 kg mass, and 45° deadrise angle at initial water entry moments: (a)t=0.15 s, (b) 0.16 s, (c) 0.17 s, (d) 0.18 s, (e) 0.19 s, and (f) 0.20 s.
(a)
(b)
(c)
(d)
(e)
(f)Comparison of local pressure (a) and strain (b) on the midpoint of wedge edge with 45° deadrise angle versus two handling method on local elastic deformation (two-way “weak” coupling and “pure” postprocessing method).
(a)
(b)At initial water-entry moments (fromt=0.15 s to 0.20 s), for the pressure distributed on elastic boundaries (as shown in Figure 10) and local pressure on the midpoint (from Figure 11(a)), the results of “pure” postprocessing treatment without feedback loop are mostly larger than the two-way “weak” coupling method, which consequently causes the strain of “pure” postprocessing treatment to irrationally and impractically enlarge during water-entry phase as shown in Figure 11(b). Hence, from Figures 10 and 11, it can be found that the feedback of structural response on CFD boundary condition is very important to ensure the calculation precision. If the feedback loop is neglected, the local structural deformation on elastic boundaries should be overvalued.In Figure12, the time history of local structural deformation and global y-displacement of the wedge vertex are calculated and drawn. From the figure, during the water entry phase, the x- and y-structural deformation and rotation-angle are under 0.08 μm, 10 μm, and 1.5 μrad, respectively. By comparing with the wedge thickness 0.0015 m, the deformations of wedge vertex are very small and can be neglected. Furthermore, from the last picture, it also can be found that the two time histories of global y-displacement of the wedge vertex versus rigid and elastic model are almost the same. Thus, a conclusion can be made that, for some water-entry cases, the vertex of elastic wedge can be handled as rigid-point and so the simplified half-wedge model can be efficiently applied in the numerical simulation to enhance computational efficiency.Figure 12
Time history of local structural deformation and global displacement on the vertex of wedge with 1.5 mm thickness, 7.787 kg mass and 45° deadrise angle.


## 5.1. Validation

In this section, the current numerical method for hydroelastic analysis on water-entry problems of free-falling wedge should be tested by comparing calculated results with experimental data of Sun et al. [55] about three different water-entry test cases.Based on the parameters of water tank in Sun’s experiment, 2D CFD computational domain was a square of 0.8 m × 1.45 m and the water depth was 1.1 m. The material properties of elastic wedges were 1.059 × 103 kg/m3 density, 2.67 × 109 kg/(ms2) elastic module, and 0.357 passion ratio. Furthermore, geometric parameters of three different wedges were in Table 1 and the diagrammatic sketch is shown in Figure 3, where L represents length of elastic wedge, B is width, and α is deadrise angle.Table 1
Geometric parameters of three different elastic wedges.
Number
Parameter
Unit
Value
1
Length
m
0.6
Width
m
0.2
Thickness
mm
1.5
Mass
kg
7.7870
Deadrise angle
45°
2
Length
m
0.6
Width
m
0.2
Thickness
mm
1.5
Mass
kg
7.4780
Deadrise angle
20°
3
Length
m
0.6
Width
m
0.2
Thickness
mm
3.0
Mass
kg
7.7890
Deadrise angle
20°Figure 3
Sketch of elastic model and description of geometric parameters.In the numerical calculation, the zero physical time (t=0) was selected as the moment in which the bottom of the wedge is 0.1 m up from the initial free surface and then the elastic wedge began to fall freely. The physical time Δt was 0.0001 s, fictitious time Δτ was 0.01 s, the artificial compressibility coefficient β was 500, the gravity acceleration was 9.81 m/s2, and the length of structural element on solid boundary was 0.01 m. Furthermore, the inhomogeneous meshes were generated in fluid field. Near solid boundary the local meshes with equal spacing Δx=Δy=0.01 m were refined and regular and then the mesh size gradually enlarged away from boundary. By taking the model 1, for example, an inhomogeneous mesh system 66 × 90 is shown in Figure 4.Global mesh (a) and local grid (b) near solid boundary.
(a)
(b)In Figure5, the time histories of global acceleration and local strain of observation point on the sloping edge versus numerical solution and experiment data for three different elastic wedges are compared. Therein Figures 5(a), 5(b), and 5(c) represent various wedge models (number 1~3 in Table 1), respectively. From the figure, it can be found that the numerical solutions by this paper and experimental data are in good agreement with each other, except acceleration peak and negative strain of elastic wedge. On the one hand, for the sake of resistance of experimental equipment, wedge accelerations in experiment before water entry are less than the numerical simulation, which should cause a little difference about acceleration peak. On the other hand, by means of numerical method, the absolute value of negative strain on solid boundary is a bit underpredicted. The reason may be that the two-way “weak” coupling strategy between fluid and local structure artificially increases the numerical stiffness of water-entry model.Time history of global acceleration and local strain of observation point on the sloping edge for three different elastic wedges, (a) 1.5 mm thickness, 7.787 kg mass, and 45° deadrise angle; (b) 1.5 mm thickness, 7.478 kg mass, and 20° deadrise angle; (c) 3.0 mm thickness, 7.789 kg mass, and 20° deadrise angle.
(a)
(b)
(c)


## 5.2. Description of Free Surface Profile and Velocity Vectors

By using the present method to study water-entry model 1 of Table1, the developing process of free surface shape and velocity vectors in flow field can be successfully numerically simulated and the results at various water-entry moments (t=0.00 s, 0.12 s, 0,18 s, 0.24 s, 0.30 s, 0.42 s, 0.48 s, and 0.60 s) are shown in Figure 6 to describe the water-entry phenomenon of the free-falling elastic wedge.Free surface profile and velocity vectors for a free-falling elastic wedge with 1.5 mm thickness, 7.787 kg mass, and 45° deadrise angle; (a)t=0.00 s, (b) 0.12 s, (c) 0.18 s, (d) 0.24 s, (e) 0.30 s, (f) 0.42 s, (g) 0.48 s, and (h) 0.60 s.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)Att=0.00 s the first graph shows the initial location of elastic wedge and calm free surface. Then the elastic wedge freely falls down (t=0.12 s) and the air near solid surface moves with the wedge. As the body pieces the free surface (t=0.18 s), the water rises up along each sloping side and the displacement waves are formed originating at the vertex of wedge. Meanwhile, the air between the solid boundary and free surface is forced around wedge edge and separates at the corner to form a pair of symmetric vortices behind the body in air. As the body travels further into the water (t=0.24 s and 0.30 s), the two intersections between wedge and free surface move up along both straight sides and subsequently overturn towards each other (t=0.42 s). After the free surface is reconnected (t=0.48 s), a jet flow is created upwards above the wedge (t=0.60 s).


## 5.3. Analysis on Hydroelasticity of Free-Falling Wedge

In order to study the global dynamic behaviors and local structural responses of elastic wedge, here, based on the geometric parameters of model 1 in Table1, two kinds of water-entry models versus rigid and elastic wedge are created simultaneously. By comparing with the numerical results of different models, the hydroelasticity of free-falling wedge during the water-entry phase is studied and analyzed.


### 5.3.1. Global Dynamic Behaviors of Elastic Wedge

Firstly, by means of the current numerical method with or without elastic deformation, the global hydrodynamic acceleration and velocity versus elastic and rigid wedge model were calculated and compared as shown in Figure7.Figure 7
Time history of global hydrodynamic acceleration and velocity versus elastic and rigid wedge model with 7.787 kg mass and 45° deadrise angle.In Figure7, for the sake of structural deformation, some high-frequency vibrations exist in the variation curves of global hydrodynamic and acceleration for the elastic model. In order to directly and clearly study the dynamic behaviors of elastic wedge, here a Savitzky-Golay method with 5 points and 4 degrees was used to filter high-frequency oscillation and smooth variable curves. For the rigid and elastic model, the time history of hydrodynamic and acceleration are in substantial agreement except two little differences. One is the peak value of hydrodynamic and acceleration curves and the results of rigid model are a little larger than the elastic model. The reason may be that the elastic deformation can provide buffer action and weaken fluid hydrodynamic during the water-entry phase. The other is that there are many high-frequency vibrations in hydrodynamic and acceleration curves, which can be caused by local deformation of elastic structure. However, the amplitudes of these vibrations are too small to influence the global motion of wedge. Thus, the velocity curves versus two models are basically in agreement with each other.On this basis, in order to further study the long global motion of free-falling wedge during the water entry phase, here the water entry model was numerically simulated for a long physical time and the results of wedge velocity were shown in the last picture of Figure7. From the picture, the motion process of free-falling wedge during water entry phase can be described as follows: firstly the wedge in air free falls straight down and then slows down due to large hydrodynamic at the initial moment of water-entry. With the continuous hydrodynamic acting on the wedge, the velocity gradually decreases until zero. Next, because the buoyancy of wedge is larger than its gravity, along with the closure of free surface, the wedge begins to move upwards. Subsequently, the jet successively grows upwards, falls downwards, and finally disappears on free surface, which will influence the global motion of wedge. Furthermore, by combining the effect of fluid force and gravity, wedge should vibrate up and down with less and less amplitude.Figure8 shows the pressure contour in fluid field at classic water-entry moments (t=0.15 s, 0.18 s, 0.21 s, and 0.24 s). In the figure, at the initial stage of water-entry (t=0.15 s), wedge has just begun to pierce free surface and the pressure peak appears near the vertex of wedge. Furthermore, the air escapes with high velocity from the gap of free surface and solid boundary to both sides and hence forms a strong air-flow to press on the wedge edges. With the wedge continuously entering water (t=0.18 s), free surface rises up along each sloping side and the amplitude and influence sphere of fluid pressure on wedge edge gradually enlarge. Furthermore, it can be found that the variation rate of pressure distribution on the wedge edge is related to the location of free surface. Therein the pressure slowly changes in water and rapidly reduces in air. As the body travels further into the water (t=0.18 s and 0.24 s), the two intersections between wedge edge and free surface move up through corners and along both straight sides. Because of the continuous effect of fluid hydrodynamic, the wedge velocity becomes smaller and smaller, which results in the reduction of pressure distribution on wedge edge. Furthermore, it should be noted that, because near the unsmooth corners the fluid velocity become very large, the negative pressure on solid boundary should take place over there.Figure 8
Pressure contour in fluid field for elastic wedge with 7.787 kg mass and 45° deadrise angle;t=0.15 s, 0.18 s, 0.21 s, and 0.24 s (from left to right and from up to bottom).


### 5.3.2. Local Structural Response on Elastic Edge

The time history and frequency analysis of local pressure and strain on the midpoint of wedge edge are shown in Figure9. In order to better discuss the characteristic of local structural deformation; here fast Fourier transform method is applied in time-frequency analysis and the local pressure and strain in frequency domain are achieved and shown. From the frequency analysis, it can be found that the variations of local pressure and stain on the midpoint mainly include two characteristics.Figure 9
Time history and frequency analysis of local pressure and strain on the midpoint of wedge edge versus elastic and rigid model with 7.787 kg mass and 45° deadrise angle.One is the pressure peak with high amplitude and narrow range in low frequency domain. Based on the effect of max pressure (frequency tends to zero) in Figure9 and global hydrodynamic in Figure 7, elastic boundaries should be deformed with 1st mode of vibration. Furthermore, the corresponding frequency of max strain is 1.3021 Hz and approximately equal to the first order natural frequency (1.2241 Hz) of the elastic model which is calculated by ANSYS software.The other is the pressure oscillation with approximately 72~85 Hz range and low amplitude in high frequency domain which includes two crests and one trough. From the modal analysis, the pressure oscillation possibly causes the midpoint to locally vibrate with 11th and 12th natural frequencies 73.13 and 86.822 Hz which are calculated by ANSYS software. Furthermore, because the two natural frequencies are very close, the vibration phenomenon of “beats” should appear on the midpoint of wedge edge. Here according to the properties of “beats” phenomenon, the theoretical frequencies versus amplitude profile(ω1-ω2)/4π and high-frequency vibration (ω1+ω2)/4π are 6.846 and 79.976 Hz, respectively. Correspondingly, from the variation of local strains in frequency and time domain as shown in Figure 9, it can be found that the calculated frequencies (7.33 and 75.52 Hz) are close to the theoretical values, which theoretically and numerically verifies the “beats” phenomenon of elastic structural vibration.Thus the local strain of midpoint on wedge edge can be divided into main deformation versus 1st principal mode and local deformation of “beat” phenomenon versus 11th and 12th modes.In order to evaluate the hydroelastic effect of free-falling wedge, here case 1 is numerically simulated without feedback loop of elastic deformation, which means that the local structural response is handled as just a “pure” postprocessing step. Then the results are shown and compared with two-way “weak” coupling method in Figures10 and 11.Pressure distribution on elastic boundaries of a free-falling elastic wedge with 1.5 mm thickness, 7.787 kg mass, and 45° deadrise angle at initial water entry moments: (a)t=0.15 s, (b) 0.16 s, (c) 0.17 s, (d) 0.18 s, (e) 0.19 s, and (f) 0.20 s.
(a)
(b)
(c)
(d)
(e)
(f)Comparison of local pressure (a) and strain (b) on the midpoint of wedge edge with 45° deadrise angle versus two handling method on local elastic deformation (two-way “weak” coupling and “pure” postprocessing method).
(a)
(b)At initial water-entry moments (fromt=0.15 s to 0.20 s), for the pressure distributed on elastic boundaries (as shown in Figure 10) and local pressure on the midpoint (from Figure 11(a)), the results of “pure” postprocessing treatment without feedback loop are mostly larger than the two-way “weak” coupling method, which consequently causes the strain of “pure” postprocessing treatment to irrationally and impractically enlarge during water-entry phase as shown in Figure 11(b). Hence, from Figures 10 and 11, it can be found that the feedback of structural response on CFD boundary condition is very important to ensure the calculation precision. If the feedback loop is neglected, the local structural deformation on elastic boundaries should be overvalued.In Figure12, the time history of local structural deformation and global y-displacement of the wedge vertex are calculated and drawn. From the figure, during the water entry phase, the x- and y-structural deformation and rotation-angle are under 0.08 μm, 10 μm, and 1.5 μrad, respectively. By comparing with the wedge thickness 0.0015 m, the deformations of wedge vertex are very small and can be neglected. Furthermore, from the last picture, it also can be found that the two time histories of global y-displacement of the wedge vertex versus rigid and elastic model are almost the same. Thus, a conclusion can be made that, for some water-entry cases, the vertex of elastic wedge can be handled as rigid-point and so the simplified half-wedge model can be efficiently applied in the numerical simulation to enhance computational efficiency.Figure 12
Time history of local structural deformation and global displacement on the vertex of wedge with 1.5 mm thickness, 7.787 kg mass and 45° deadrise angle.


## 5.3.1. Global Dynamic Behaviors of Elastic Wedge

Firstly, by means of the current numerical method with or without elastic deformation, the global hydrodynamic acceleration and velocity versus elastic and rigid wedge model were calculated and compared as shown in Figure7.Figure 7
Time history of global hydrodynamic acceleration and velocity versus elastic and rigid wedge model with 7.787 kg mass and 45° deadrise angle.In Figure7, for the sake of structural deformation, some high-frequency vibrations exist in the variation curves of global hydrodynamic and acceleration for the elastic model. In order to directly and clearly study the dynamic behaviors of elastic wedge, here a Savitzky-Golay method with 5 points and 4 degrees was used to filter high-frequency oscillation and smooth variable curves. For the rigid and elastic model, the time history of hydrodynamic and acceleration are in substantial agreement except two little differences. One is the peak value of hydrodynamic and acceleration curves and the results of rigid model are a little larger than the elastic model. The reason may be that the elastic deformation can provide buffer action and weaken fluid hydrodynamic during the water-entry phase. The other is that there are many high-frequency vibrations in hydrodynamic and acceleration curves, which can be caused by local deformation of elastic structure. However, the amplitudes of these vibrations are too small to influence the global motion of wedge. Thus, the velocity curves versus two models are basically in agreement with each other.On this basis, in order to further study the long global motion of free-falling wedge during the water entry phase, here the water entry model was numerically simulated for a long physical time and the results of wedge velocity were shown in the last picture of Figure7. From the picture, the motion process of free-falling wedge during water entry phase can be described as follows: firstly the wedge in air free falls straight down and then slows down due to large hydrodynamic at the initial moment of water-entry. With the continuous hydrodynamic acting on the wedge, the velocity gradually decreases until zero. Next, because the buoyancy of wedge is larger than its gravity, along with the closure of free surface, the wedge begins to move upwards. Subsequently, the jet successively grows upwards, falls downwards, and finally disappears on free surface, which will influence the global motion of wedge. Furthermore, by combining the effect of fluid force and gravity, wedge should vibrate up and down with less and less amplitude.Figure8 shows the pressure contour in fluid field at classic water-entry moments (t=0.15 s, 0.18 s, 0.21 s, and 0.24 s). In the figure, at the initial stage of water-entry (t=0.15 s), wedge has just begun to pierce free surface and the pressure peak appears near the vertex of wedge. Furthermore, the air escapes with high velocity from the gap of free surface and solid boundary to both sides and hence forms a strong air-flow to press on the wedge edges. With the wedge continuously entering water (t=0.18 s), free surface rises up along each sloping side and the amplitude and influence sphere of fluid pressure on wedge edge gradually enlarge. Furthermore, it can be found that the variation rate of pressure distribution on the wedge edge is related to the location of free surface. Therein the pressure slowly changes in water and rapidly reduces in air. As the body travels further into the water (t=0.18 s and 0.24 s), the two intersections between wedge edge and free surface move up through corners and along both straight sides. Because of the continuous effect of fluid hydrodynamic, the wedge velocity becomes smaller and smaller, which results in the reduction of pressure distribution on wedge edge. Furthermore, it should be noted that, because near the unsmooth corners the fluid velocity become very large, the negative pressure on solid boundary should take place over there.Figure 8
Pressure contour in fluid field for elastic wedge with 7.787 kg mass and 45° deadrise angle;t=0.15 s, 0.18 s, 0.21 s, and 0.24 s (from left to right and from up to bottom).


## 5.3.2. Local Structural Response on Elastic Edge

The time history and frequency analysis of local pressure and strain on the midpoint of wedge edge are shown in Figure9. In order to better discuss the characteristic of local structural deformation; here fast Fourier transform method is applied in time-frequency analysis and the local pressure and strain in frequency domain are achieved and shown. From the frequency analysis, it can be found that the variations of local pressure and stain on the midpoint mainly include two characteristics.Figure 9
Time history and frequency analysis of local pressure and strain on the midpoint of wedge edge versus elastic and rigid model with 7.787 kg mass and 45° deadrise angle.One is the pressure peak with high amplitude and narrow range in low frequency domain. Based on the effect of max pressure (frequency tends to zero) in Figure9 and global hydrodynamic in Figure 7, elastic boundaries should be deformed with 1st mode of vibration. Furthermore, the corresponding frequency of max strain is 1.3021 Hz and approximately equal to the first order natural frequency (1.2241 Hz) of the elastic model which is calculated by ANSYS software.The other is the pressure oscillation with approximately 72~85 Hz range and low amplitude in high frequency domain which includes two crests and one trough. From the modal analysis, the pressure oscillation possibly causes the midpoint to locally vibrate with 11th and 12th natural frequencies 73.13 and 86.822 Hz which are calculated by ANSYS software. Furthermore, because the two natural frequencies are very close, the vibration phenomenon of “beats” should appear on the midpoint of wedge edge. Here according to the properties of “beats” phenomenon, the theoretical frequencies versus amplitude profile(ω1-ω2)/4π and high-frequency vibration (ω1+ω2)/4π are 6.846 and 79.976 Hz, respectively. Correspondingly, from the variation of local strains in frequency and time domain as shown in Figure 9, it can be found that the calculated frequencies (7.33 and 75.52 Hz) are close to the theoretical values, which theoretically and numerically verifies the “beats” phenomenon of elastic structural vibration.Thus the local strain of midpoint on wedge edge can be divided into main deformation versus 1st principal mode and local deformation of “beat” phenomenon versus 11th and 12th modes.In order to evaluate the hydroelastic effect of free-falling wedge, here case 1 is numerically simulated without feedback loop of elastic deformation, which means that the local structural response is handled as just a “pure” postprocessing step. Then the results are shown and compared with two-way “weak” coupling method in Figures10 and 11.Pressure distribution on elastic boundaries of a free-falling elastic wedge with 1.5 mm thickness, 7.787 kg mass, and 45° deadrise angle at initial water entry moments: (a)t=0.15 s, (b) 0.16 s, (c) 0.17 s, (d) 0.18 s, (e) 0.19 s, and (f) 0.20 s.
(a)
(b)
(c)
(d)
(e)
(f)Comparison of local pressure (a) and strain (b) on the midpoint of wedge edge with 45° deadrise angle versus two handling method on local elastic deformation (two-way “weak” coupling and “pure” postprocessing method).
(a)
(b)At initial water-entry moments (fromt=0.15 s to 0.20 s), for the pressure distributed on elastic boundaries (as shown in Figure 10) and local pressure on the midpoint (from Figure 11(a)), the results of “pure” postprocessing treatment without feedback loop are mostly larger than the two-way “weak” coupling method, which consequently causes the strain of “pure” postprocessing treatment to irrationally and impractically enlarge during water-entry phase as shown in Figure 11(b). Hence, from Figures 10 and 11, it can be found that the feedback of structural response on CFD boundary condition is very important to ensure the calculation precision. If the feedback loop is neglected, the local structural deformation on elastic boundaries should be overvalued.In Figure12, the time history of local structural deformation and global y-displacement of the wedge vertex are calculated and drawn. From the figure, during the water entry phase, the x- and y-structural deformation and rotation-angle are under 0.08 μm, 10 μm, and 1.5 μrad, respectively. By comparing with the wedge thickness 0.0015 m, the deformations of wedge vertex are very small and can be neglected. Furthermore, from the last picture, it also can be found that the two time histories of global y-displacement of the wedge vertex versus rigid and elastic model are almost the same. Thus, a conclusion can be made that, for some water-entry cases, the vertex of elastic wedge can be handled as rigid-point and so the simplified half-wedge model can be efficiently applied in the numerical simulation to enhance computational efficiency.Figure 12
Time history of local structural deformation and global displacement on the vertex of wedge with 1.5 mm thickness, 7.787 kg mass and 45° deadrise angle.


## 6. Conclusion

(1)
In this paper, surface capturing method and Cartesian cut cell mesh are successfully developed to handle the interaction among fluid field, global motion, and local deformation of elastic body simultaneously for water entry model of free-falling wedge.
(2)
In finite element analysis, based on the different characteristic of fluid and structural mesh, a particular data transfer method on solid boundary is successfully deduced and presented. Furthermore, an effective fluid-structure coupled solution strategy is taken to solve the interaction among fluid field, global motion, and local deformation of elastic body.
(3)
Based on the present method, some test cases of water entry for free-falling elastic wedge are numerically simulated. By comparing with experimental data, the results show that the present method has a good ability to study hydroelastic water-entry problems of free-falling elastic wedge.
(4)
Through further analysis, the global hydrodynamic performance and physical phenomenon of water entry problem (such as free surface rising, pressure distribution, jet-flow, and negative pressure) for elastic wedge are discussed in detail.
(5)
For the local structural response on elastic boundaries, especial “beats” phenomenon of elastic vibration versus certain model without structural damping are discussed. Furthermore, based on the results of vertex for elastic wedge, a conclusion has been made that, for some test cases of elastic wedge, simplified half model can be applied in the numerical simulation to enhance computational efficiency.Furthermore, the method should be developed further on the following aspects.(1)
Although the present method is available and accurate for the water-entry model of wedge, it has trouble in water-entry problems of solid body with complex shape because of the invalidation of lath-beam structural model and a more finite element model (shell or plate) should be applied. Furthermore, the current fluid-structure coupling strategy should be improved and gradually transformed into other advanced coupling approaches (such as two-way “tight” coupling method).(2)
Because here, the structural finite element analysis is based on linear elastic theory, the present method is not suitable for water entry problems of elastic body with large deformation. Hence in future the geometric nonlinearity should be taken into account to extend the application fields.



---

*Source: 545642-2014-07-23.xml*