A finite element method for fluid-structure interaction with surface waves

==1 Introduction==

Accurate prediction of the fluid-structure interaction effects for a totally or partially submerged body in a flowing liquid including a free surface is a problem if great relevance in offshore engineering and naval architecture among many other fields.

The difficulties in accurately solving the coupled fluid-structure interaction problem in this case are mainly due to the following reasons:

'''1.'''  The difficulty of solving numerically the incompressible fluid dynamic equations which typically include intrinsic non linearities except for the simplest and limited potential flow model.

'''2.''' The obstacles in solving the constraint equation stating that at the free surface boundary the fluid particles remain on that surface which position is in turn unknown.

'''3.'''  The difficulties in solving the problem of motion of the submerged body due to the interaction forces while minimizing the distorsion of the finite elements discretizing the fluid domain thus reducing the need of remeshing.

This paper presents a stabilized finite element method which allows to overcome above three obstacles. The starting point are the modified governing differential equations for the incompressible viscous flow and the free surface condition incorporating the necessary stabilization terms via a ''finite increment calculus'' (FIC) procedure developed by the authors [1&#8211;4]. The FIC approach has been successfully applied to the finite element and meshless solution of a range of advective-diffusive transport and fluid flow problems [1&#8211;5].

The modified governing equations are solved in space-time using a semi-implicit fractional step approach and the finite element method (FEM). Free surface wave effects are accounted for via the introduction of a prescribed pressure at the free surface computed from the wave height.

The movement of the submerged body within the fluid due to the interaction forces is treated by solving a structural dynamic problem using the fluid forces as input loads. A method to update the mesh for the fluid domain following the movement of the submerged body which minimizes element distorsion is presented. The mesh update procedure is based on the iterative finite element solution of a linear elastic problem on the mesh domain where fictitions elastic properties are assigned so that elements suffering higher movements are stiffer [6].

The content of the paper is structured as follows. First details of the stabilized semi-implicit fractional step approach using the FEM is described. Next the mesh updating procedure is presented. Finally some examples of a coupled fluid-interaction problem are given.

==2 Stabilized finite element formulation for the fluid flow equations == 
       
We consider the motion around a body of a viscous incompressible fluid including a free surface.

The stabilized form of the governing differential equations for the three dimensional  (3D) problem can be written as

''Momentum''

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>r_{m_i} - \underline{{1\over 2} h_{mj} {\partial r_{m_i}\over \partial x_ j}}{1\over 2} h_{mj} {\partial r_{m_i}\over \partial x_ j} =0 \qquad \hbox{on } \Omega \quad i,j=1,2,3 </math>
|}
| style="width: 5px;text-align: right;" |  (1a)
|}

''Mass balance''

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>r_d - \underline{{1\over 2} h_{dj}{\partial r_d \over \partial x_j}}{1\over 2} h_{dj}{\partial r_d \over \partial x_j} =0 \qquad \hbox{on } \Omega \quad j=1,2,3</math>
|}
| style="width: 5px;text-align: right;" |  (1b)
|}

''Free surface''

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>r_\beta - \underline{{1\over 2} h_{\beta _j}{\partial r_\beta \over \partial x_j}}{1\over 2} h_{\beta _j}{\partial r_\beta \over \partial x_j} =0 \qquad \hbox{on } \Gamma _\beta \quad j=1,2</math>
|}
| style="width: 5px;text-align: right;" |  (1c)
|}

where

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math> r_{m_i} =</math>
| style="text-align: right;" | <math>\rho \left[{\partial u_i\over \partial t} + {\partial \over \partial x_j} (u_iu_j)\right]+ {\partial p\over \partial x_i} - {\partial \tau _{ij}\over \partial x_j} -b_i</math>
| <math> (2a)</math>
|-
| style="text-align: center;" | <math> r_d=</math>
| style="text-align: right;" | <math>\rho {\partial u_i\over \partial x_i} \qquad i=1,2,3 </math>
| <math>(2b)</math>
|-
| style="text-align: center;" | <math> r_\beta =</math>
| style="text-align: right;" | <math>{\partial \beta \over \partial t} + u_i {\partial \beta \over \partial x_i}-u_3 \qquad i=1,2 </math>
| <math>(2c)</math>
|-
| style="text-align: center;" | 
|}
|}

In above <math display="inline">u_i</math> is the velocity along the ith global reference axis, <math display="inline">\rho </math> the (constant) density of the fluid, <math display="inline">p</math> the pressure, <math display="inline">\beta </math> the wave elevation, <math display="inline">b_i</math> the body forces acting in the fluid and <math display="inline">\tau _{ij}</math> the viscous stresses related to the viscosity <math display="inline">\mu </math> by the standard expression

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>\tau _{ij} = \mu \left({\partial u_i\over \partial x_j}+{\partial u_j\over \partial x_i}-\delta _{ij} {2\over 3} {\partial u_k\over \partial x_k}\right)</math>
|}
| style="width: 5px;text-align: right;" |  (3)
|}

The underlined terms in eqs.(1) introduce the necessary stabilization for the approximated numerical solution.

The distances <math display="inline">h_{mj},h_{dj}</math> and <math display="inline">h_{\beta _j}</math> are termed ''characteristic lengths''  and represent the dimensions of the finite domain where balance of momentum, mass and transport of fluid particles is enforced. Details of the derivation of eqs.(2) can be found in [1].

A more convenient form of equation (1b) can be written by assuming <math display="inline">h_{dj}=-2\tau _du_j</math> where <math display="inline">\tau _d</math> is an  intrinsic time parameter. Under this assumption and using eq.(1a) the stabilized form of the mass balance equation can be written as (neglecting high order terms) [4]

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>r_d -\tau _d  {\partial \bar r_{m_i}\over \partial x_i}=0 </math>
|}
| style="width: 5px;text-align: right;" |  (4)
|}

where

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>\bar r_{m_i} =\rho \left({\partial u_i\over \partial t} + u_j{\partial u_i\over \partial x_j}\right)+ {\partial p\over \partial x_i} - {\partial \tau _{ij}\over \partial x_j}    </math>
|}
| style="width: 5px;text-align: right;" |  (5)
|}

The boundary conditions for the stabilized problem are written as

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>n_j \tau _{ij} + t_i + \underline{{1\over 2} h_{mj} n_j r_{m_i}}{1\over 2} h_{mj} n_j r_{m_i} =0 \qquad \hbox{on }\Gamma _t </math>
|}
| style="width: 5px;text-align: right;" |  (6a)
|}

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>u_j-u_j^p=0 \qquad \hbox{on }\Gamma _u </math>
|}
| style="width: 5px;text-align: right;" |  (6b)
|}

where <math display="inline">n_j</math> are the components of the unit normal vector to the boundary and <math display="inline">t_i</math> and <math display="inline">u_j^p</math> are prescribed tractions and displacements on the boundaries <math display="inline">\Gamma _t </math> and <math display="inline">\Gamma _u</math>, respectively. The underlined stabilized tersm appearing in  the Neumann boundary condition (6a) are obtained via the FIC approach [1].

Eqs.(1&#8211;6) are the starting point for deriving a variety of stabilized numerical methods for solving the incompressible Navier-Stokes equations. It can be shown that a number of standard stabilized finite element methods allowing equal order interpolations for the velocity and pressure fields can be recovered from the modified form of the momentum and mass balance equations given above [4]. A semi-implicit fractional step finite element procedure for solution of eqs.(1a),(1c),(4) and (6) is presented in next section.

===2.1 Stabilized fractional step method=== 

Let us discretize in time the stabilized momentum equation (1a) as

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>\rho \left[{u_i^{n+1}-u_i^n\over \Delta t} + {\partial \over \partial x_j}(u_iu_j)^n\right]+ {\partial p^{n+1} \over \partial x_i}-{\partial \tau ^n_{ij}\over \partial x_j}-b_i^n - {1\over 2} h_{mj} {\partial r_{m_i}^n \over \partial x_j}=0</math>
|}
| style="width: 5px;text-align: right;" |  (7)
|}

A fractional step method (also termed “segregation” or “splitting” procedure) can be simply derived by splitting eq.(7) as follows

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math> u_i^* =</math>
| style="text-align: right;" | <math> u_i^n -\Delta t \left[{\partial \over \partial x_j}(u_iu_j)- {1\over \rho } {\partial \tau _{ij}\over \partial x_j}-{1\over \rho }b_i - {1\over 2\rho } h_{mj} {\partial r_{m_i}^n \over \partial x_j}\right]^n </math>
| <math>(8a)</math>
|-
| style="text-align: center;" | <math> u_i^{n+1}=</math>
| style="text-align: right;" | <math> u_i^* - {\Delta t\over \rho } {\partial p^{n+1}\over \partial x_i}</math>
| <math>(8b)</math>
|-
| style="text-align: center;" | 
|}
|}

Note that addition of eqs.(8a) and (8b) gives the original stabilized momentum equation (7).

Substitution of eq.(8b) into eq.(4) gives after some algebra [4]

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>(\Delta t +\tau _d) \Delta p^{n+1}= {\partial u_i^*\over \partial x_i}- \tau _d {\partial g_i^n\over \partial x_i}=0</math>
|}
| style="width: 5px;text-align: right;" |  (9)
|}

where

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>g_i=\rho \left({\partial u_i\over \partial t}+u_j {\partial u_i\over \partial x_j}\right)- {\partial \tau _{ij}\over \partial x_j}-b_i</math>
|}
| style="width: 5px;text-align: right;" |  (10)
|}

Standard fractional step procedures neglect the contribution from the terms involving <math display="inline">\tau _d</math> in eq.(9). It can be shown that these terms have an additional stabilization effect which improves the numerical solution when the values of <math display="inline">\Delta t</math> are small.

The free surface wave equation (1c) can be also discretized in time to give

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>\beta ^{n+1}=\beta ^n - \Delta t \left[u_i^{n+1} {\partial \beta ^n \over \partial x_i}-u_3^{n+1} - {1\over 2}h_{\beta j} {\partial r^n_\beta \over \partial x_j}\right]\qquad i,j=1,2 </math>
|}
| style="width: 5px;text-align: right;" |  (11)
|}

A typical solution in time includes the following steps.

''Step 1''. Solve explicitely for the so called fractional velocities <math display="inline">u^*_i</math> using eq.(8a).

''Step 2''. Solve for the pressure field <math display="inline">p^{n+1}</math> solving the Laplacian equation (9). The pressures at the free surface computed from step 4 below in the previous time step are used as boundary conditions for solution of eq.(9).

''Step 3''. Compute the movement of the submerged body by solving the dynamic equations of motion in the body subjected to the pressure field <math display="inline">p^{n+1}</math> and the viscous stresses <math display="inline">\tau ^n_{ij}</math>.

''Step 4''. Compute the new position of mesh nodes <math display="inline">{x}_j^{n+1}</math> in the fluid domain by using the mesh update algorithm described in next section.

''Step 5''. Compute the fractional velocity and pressure fields at the new position of the nodes. This can be simply done by the following expression

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math> (\hat u_i^*)_j = </math>
| style="text-align: right;" | <math>(u_i^*)_j - [{d}_j^{n+1}]^T ({              \nabla \hbox{ }}{u}_i^*)_j \quad i=1,2,3 </math>
| <math>(12a)</math>
|-
| style="text-align: center;" | <math>  \hat p_j^{n+1}=</math>
| style="text-align: right;" | <math>p_j^{n+1}-[{d}_j^{n+1}]^T ({              \nabla \hbox{ }}{p}_i^n)_j </math>
| <math>(12b)</math>
|-
| style="text-align: center;" | 
|}
|}

In eq.(12) <math display="inline">j</math> is the node number, <math display="inline">\hat {(\cdot )}_j</math> are values of <math display="inline">u^*_i</math> in the updated nodes after mesh movement whereas the r.h.s. includes values of <math display="inline">u^*_i</math> and <math display="inline">p</math> in the original position, <math display="inline">{d}_j=({x}_j^{n+1}-{x}_j^n)</math> and <math display="inline">{              \nabla \hbox{ }}</math> is the gradient vector.

An alternative to eq.(12b) is to compute the nodal pressures in the updated mesh configuration by solving once more eq.(9) using the values of <math display="inline">\hat u_i^*</math> obtained from eq.(12a). This will ensure better satisfaction of the incompressibility condition in the updated mesh.

''Step 6''. Compute the velocity field <math display="inline">u^{n+1}_i</math> at the updated configuration for each mesh node

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>u^{n+1}_i=\hat u_i^* - \Delta t {\partial \hat p^{n+1}\over \partial x_i}</math>
|}
| style="width: 5px;text-align: right;" |  (13)
|}

''Step 7''. Solve for the updated value of the free surface elevation <math display="inline">\beta ^{n+1}</math> using eq.(11). Compute the pressure in the free surface from Benouilli equation as

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>p^{n+1}=p^o +\rho g (\beta ^{n+1}-\beta ^o)</math>
|}
| style="width: 5px;text-align: right;" |  (14)
|}

where <math display="inline">\beta ^o</math> and <math display="inline">p^o</math> are reference values of the free surface elevation and the pressure respectively and <math display="inline">g</math> is the gravity constant.

As already mentioned the effect of changes in the free surface elevation is introduced in the  step 2 of the flow solution  as a prescribed pressure acting on the free surface.

The accuracy of above transient solution process depends on the time step size which should satisfy stability criteria for the coupled solution. Indeed larger time steps can be used if the values at time <math display="inline">n</math> in above equations are computed at <math display="inline">n+1/2</math>. The solution process becomes now implicit and an iteration loop within each time step  is then required.

===2.2 Finite element discretization===

Space discretization is carried out using the finite element method [7]. The velocity and pressure fields are interpolated within each element in the standard finite element manner as

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math> u_i=</math>
| style="text-align: right;" | <math>\sum \limits _j N_{u_j} \bar{(u_i)}_j</math>
|-
| style="text-align: center;" | <math> p_i=</math>
| style="text-align: right;" | <math>\sum \limits _j N_{p_j}\bar p_j</math>
|-
| style="text-align: center;" | 
|}
| style="width: 5px;text-align: right;" |  (15)
|}

where <math display="inline">N_{u_j}</math> and <math display="inline">N_{p_j}</math> are the shape functions interpolating the velocity and pressure fields, respectively and <math display="inline">\bar{(\cdot )}</math> denote nodal values.

Similarly the wave height is discretized as

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>\beta =\sum \limits _j N_{\beta _j}\bar \beta _j </math>
|}
| style="width: 5px;text-align: right;" |  (16)
|}

where <math display="inline">N_{\beta _j}</math> are shape functions defined over the nodes discretizing the free surface.

It is worth noting that the stabilized formulation described allows an equal order interpolation of velocities and pressure [4]. A linear interpolation over triangles (2D) and tetrahedra (3D) for both <math display="inline">u_i</math> and <math display="inline">p</math> is chosen in the examples shown in the paper. Similarly linear elements are chosen to interpolate <math display="inline">\beta </math> on the free surface mesh.

The discretized integral form is obtained by applying the standard Galerkin procedure to eqs.(8a),(8b),(9),(11),(12) and the boundary conditions (6a). The resulting expressions  follow the pattern given in [4].

The computation of the stabilization parameters <math display="inline">h_{m_j},h_{\beta _j}</math> and <math display="inline">\tau _d</math> can be based in the diminishing residual technique explained in [1&#8211;5]. In the example presented in Section&nbsp;4 the simpler option <math display="inline">h_{m_j}=h_{\beta _j}=h^{(e)}</math> and <math display="inline">\tau _d={h^{(e)}\over 2\vert{u}\vert }</math> with <math display="inline">h^{(e)}=[V^{(e)}]^{1/3}</math> where <math display="inline">V^{(e)}</math> is the element volume has been taken.


==3 A simple algorithm for stable updating of mesh nodes==                                 

Finite element solution of fluid-structure interaction problems usually requires the update of the analysis mesh as described in previous section. A typical example is the study of movement of an object within a flowing liquid where the fluid mesh needs to be continuously updated accordingly to the changes in position of the object due to the interaction forces.

Chiandussi, Bugeda and Oñate [6] have recently proposed a simple method for movement of mesh nodes ensuring minimum element distorsion. The method is based on the iterative solution of a fictitions linear elastic problem on the mesh domain. In order to minimize mesh deformation the “elastic” properties of each mesh element are appropiately selected so that elements suffering greater movements are stiffer. The basis of the method is given below.

Let us consider an elastic domain with homogeneous isotropic elastic properties characterized by the Young modulus <math display="inline">\bar E</math> and the Poisson coefficient <math display="inline">\nu </math>. Once a discretized finite element problem has been solved using, for instance, standard <math display="inline">C_o</math> linear triangles (in 2D) or linear tetraedra (in 3D), the principal stresses <math display="inline">{}^1\sigma _i</math> at the center of each element are obtained as

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>{}^1\sigma _i =\bar E [\varepsilon _i - \nu (\varepsilon _j + \varepsilon _k)] \qquad i,j=1,2,3 ~~\hbox{for 3D }</math>
|}
| style="width: 5px;text-align: right;" |  (17)
|}

where <math display="inline">\varepsilon _i</math> are the principal strains.

Let us assume now that a uniform strain field <math display="inline">\varepsilon _i =\bar \varepsilon </math> throughout the mesh is sougth. The principal stresses are then given by

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>{}^2\sigma _i = E \bar \varepsilon (1-2\nu ) \qquad i=1,2,3 ~~\hbox{for 3D }</math>
|}
| style="width: 5px;text-align: right;" |  (18)
|}

where <math display="inline">E</math> is the unknown Young modulus for the element.

A number of criteria can be now used to find the value of <math display="inline">E</math>. The most effective approach found in [6] is to equal the element strain energies in both analysis. Thus

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math> U_1=</math>
| style="text-align: right;" | <math> {}^1\sigma _i \varepsilon _i = \bar E [(\varepsilon _1^2+\varepsilon _2^2 + \varepsilon _3^2)-2\nu (\varepsilon _1 \varepsilon _2+ \varepsilon _2 \varepsilon _3+ \varepsilon _1 \varepsilon _3)] </math>
| <math>(19)</math>
|-
| style="text-align: center;" | <math> U_2=</math>
| style="text-align: right;" | <math> {}^2\sigma _i \varepsilon _i = 3E \bar  \varepsilon ^2  (1-2\nu )</math>
| <math>(20)</math>
|-
| style="text-align: center;" | 
|}
|}

Equaling eqs.(19) and (20) gives the sought Young modulus <math display="inline">E</math> as

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>E = {\bar E \over 3 \bar  \varepsilon ^2 (1-2\nu )}[(\varepsilon _1^2+\varepsilon _2^2 + \varepsilon _3^2)-2\nu (\varepsilon _1 \varepsilon _2+ \varepsilon _2 \varepsilon _3+ \varepsilon _1 \varepsilon _3)]</math>
|}
| style="width: 5px;text-align: right;" |  (21)
|}

Note that the element Young modulus is proportional to the element deformation as desired. Also recall that both <math display="inline">\bar E </math> and <math display="inline">\bar  \varepsilon </math> are constant for all elements in the mesh.

The solution process includes the following two steps.

''Step 1''. Consider the finite element mesh as a linear elastic solid with homogeneous material properties characterized by <math display="inline">\bar E</math> and <math display="inline">\nu </math>. Solve the corresponding elastic problem with imposed displacements at the mesh boundary. These displacements can be due to a prescribed motion of a body within a fluid, to changes in the shape of the domain in an optimum design problem, etc.

''Step 2''. Compute the principal strains  and the values of the new Young modulus in each element using eq.(21) for a given value of <math display="inline">\bar \varepsilon </math>. Repeat the finite element solution of the linear elastic problem with prescribed boundary displacements using the new values of <math display="inline">E</math> for each element.

The movement of the mesh nodes obtained in the second step ensures a quasi uniform mesh distorsion. Further details on this method including other alternatives for evaluating the Young modulus <math display="inline">E</math> can be found in [6].

The previous algorithm for movement of mesh nodes is able to treat the movement of the mesh due to changes in position of fully submerged and semi-submerged bodies. Note however that if the floating body intersects the free surface, the changes in the analysis domain geometry can be very important. From one time step to other emersion or inmersion of significant parts of the body can occur.

A posible solution to this problem is to remesh the analysis domain. However for m      ost  problems, a mapping of the moving surfaces linked to mesh updating algorithm described above can avoid remeshing (Figure 1).

The surface mapping technique used in this work is based on transforming 3D curved surfaces into reference planes. This allows to compute within each plane the local (in-plane) coordinates of the nodes for the final surface mesh accordingly to the changes in the floating line. The final step is to transform back the local coordinates of the surface mesh in the reference plane to the final curved configuration which incorporates the new floating line.

[[File:Example.jpg]]

<span style="text-align: center; font-size: 75%;">                                  =                                                                                                                                                                                                                                              >                                                                                                          =                                                                                                                                                                                                                                                  Figure 1.&nbsp;  Changes in the fluid interface in a floating body.                                 </span>           
         
==4 Example. Movement of a submerged sphere in an open channel==                                  

Figure 2 shows the geometry of the channel and the position of the sphere of 2m diameter with a weight of 1000 N and a rotational inertia of 1000 kgm<math display="inline">^2</math>. A mesh of 19870 linear tetrahedra with 4973 nodes has been used for the analysis.

The problem has been analyzed for values of Reynolds number = 200 and Froude number = 0.71 corresponding to velocity of 1m/s at the inlet.

It is assumed that the sphere can only move vertically and rotate a fluid around the global <math display="inline">y</math> axes due to the forces induced by the fluid. The vertical displacement is constrained by a spring linking the sphere to the ground. An initial vertical velocity of 1m/s for the sphere has been taken.

Figure 3 shows a plot of the time evolution of the vertical displacement of the sphere. The position of the sphere at different time intervals is shown in Figure&nbsp;4. A plot of velocity vectors in the fluid displayed on Figure&nbsp;5. Pressure contours in the fluid domain are shown in Figure&nbsp;6.

<span style="text-align: center; font-size: 75%;">                                  =                                                                                                                                                                                                                                              >                                                                                                          =                                                                                                                                                                                                                                                  Figure 2.&nbsp;  Geometry of the chanel with submerged sphere.                                 </span>                   

<span style="text-align: center; font-size: 75%;">                                  =                                                                                                                                                                                                                                              >                                                                                                          =                                                                                                                                                                                                                                                  Figure 3.&nbsp;  Time evolution of vertical displacement of shpere.                                 </span>                   

The streamlines in the vecinity of the sphere at a certain time are shown in Figure 7. Figure 8 shows contours of the surface wave elevation at two particular times.

<span style="text-align: center; font-size: 75%;">                                  =                                                                                                                                                                                                                                              >                                                                                                          =                                                                                                                                                                                                                                                  Figure 4.&nbsp;  Position of sphere and mesh at different times. a) t=0.47s, b) t=0.94s, c) t=1.22s, d) t=1.83s, e) t=2.57s, f) t=3.16s.                                 </span>                                                                                                                                                                                                                       >                                                                           >          <br/>                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                     

<span style="text-align: center; font-size: 75%;">                                  =                                                                                                                                                                                                                                              >                                                                                                          =                                                                                                                                                                                                                                                  Figure 5.&nbsp;  Velocity vectors at different times. a) t=0.47s, b) t=0.94s, c) t=1.22s, d) t=1.83s, e) t=2.57s, f) t=3.16s.                                 </span>                                                                                                                                                                                                                       >                                                                           >          <br/>                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                     

<span style="text-align: center; font-size: 75%;">                                  =                                                                                                                                                                                                                                              >                                                                                                          =                                                                                                                                                                                                                                                  Figure 6.&nbsp;  Pressure distribution in fluid domain at different times. a) t=0.47s, b) t=0.94s, c) t=1.22s, d) t=1.83s, e) t=2.32s, f) t=2.85s.                                 </span>                                               <span style="text-align: center; font-size: 75%;">                                  =                                                                                                                                                                                                                                              >                                                                                                          =                                                                                                                                                                                                                                                  Figure 7.&nbsp;  Streamline tracks at times: t=1.83s                                 </span>                                                                                                                                                                                                                       >                                                                           >          <br/>                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                     

<span style="text-align: center; font-size: 75%;">                                  =                                                                                                                                                                                                                                              >                                                                                                          =                                                                                                                                                                                                                                                  Figure 8.&nbsp;  Free surface elevation at different times: a) t=0.47s, b) t=3.16s.                                 </span>                                                                                                                                                                                                                       >                                                                           >          <br/>                                                                                                                                                                                                                                                                                                                                                                 ==5 Conclusions==                                      

A stabilized semi-implicit fractional step finite element method for analysis of coupled fluid interaction problems involving free surface waves has been presented. A procedure for automatic movement of mesh nodes during the coupled solution process has been developed. The method is adequate for solving large scale fluid-structure interaction situations  in naval architecture and offshare engineering problems.

==6 Acknowledgements==                             

This work was partially supported by Empresa Nacional BAZAN de Construcciones Navales y Militares S.A. This support is gratefully acknowledged.

==7. References==      

'''1.'''  E. Oñate - Derivation of stabilized equations for  advective-diffusive transport and fluid flow problems.  ''Comput. Meth. Appl. Mech. Engng.'',  Vol. 151, 1-2, pp. 233&#8211;267, 1998.

'''2.''' E. Oñate, J. Garcia and S. Idelsohn - Computation  of the stabilization  parameter for the finite  element solution of advective-diffusive problems.  ''Int. J. Num. Meth. Fluids'', Vol. 25, pp. 1385&#8211;1407, 1997.

'''3.'''  E. Oñate, J. Garcia and S. Idelsohn - An alpha-adaptive  approach for stabilized finite element solution of advective-diffusive  problems with sharp gradients. ''New Adv. in Adaptive  Comp. Met. in Mech.'', P. Ladeveze and J.T. Oden (Eds.),  Elsevier, 1998.

'''4.'''  E. Oñate - A finite element method for incompressible viscous flows using a finite increment calculus formulation. ''Research Report'' N. 150, CIMNE, Barcelona, January 1999.

'''5.'''  E. Oñate and S. Idelsohn - A mesh free finite point method for advective-diffusive transport and fluid flow problems. ''Computational Mechanics'', 21, 283&#8211;292, 1988.

'''6.'''  Chiandusi, G., Bugeda, G. and Oñate, E. - A simple method for update of finite element meshes. Research Report 147, CIMNE, Barcelona, January 1999.

'''7.'''   Zienkiewicz, O.C. and Taylor, R.C. - ''The finite element method'', 4th Edition, Vol. 1, McGraw Hill, 1989.