Congress abstract test 1

==A FINITE ELEMENT METHOD FOR FLUID-STRUCTURE INTERACTION WITH SURFACE WAVES ==

'''E. Oñate '''and''' J. García'''

{|
|-
|International Centre for Numerical Methods in Engineering
|-
| Universidad Politécnica de Cataluña 
|-
| Gran Capitán s/n, 08034 Barcelona, Spain 
|-
| E&#8211;mail: [mailto:onate@cimne.upc.es onate@cimne.upc.es]
|-
| Web page: [http://www.cimne.upc.es www.cimne.upc.es]
|}

==Abstract==

A stabilized semi-implicit fractional step finite element method for solving coupled fluid-structure interaction problems involving free surface waves is presented. The stabilized governing equations for the viscous incompressible fluid and the free surface are derived at a differential level via a finite  calculus procedure. A mesh updating technique based on solving a fictitious elastic problem on the moving mesh is described. Examples of the efficiency of the stabilized semi-implicit algorithm for the analysis of fluid-structure interaction problems in totally or partially submerged bodies is presented.

'''keywords''' Fluid-structure interaction, surface waves, finite element method.

==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 of great relevance in civil and 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:

<ol>

<li> The difficulty of solving numerically the incompressible fluid dynamic equations which typically include intrinsic nonlinearities except for the simplest and limited potential flow model.    </li>

<li>The obstacles in solving the constraint equation stating that the fluid particles remain on the free surface boundary which position is in turn unknown.   </li>

<li>The difficulties in predicting the motion of the submerged body due to the interaction forces while minimizing the distortion of the finite elements discretizing the fluid domain, thus reducing the need for remeshing.   </li>

</ol>

This paper extends recent work of the authors <span id='citeF-1'></span>[[#cite-1|1]] to derive a stabilized finite element method which overcomes the 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 calculus'' (FIC) procedure developed by the authors <span id='citeF-6'></span>[[#cite-6|6]]. 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 <span id='citeF-1'></span>[[#cite-1|1]].

The stabilized governing equations are written in an arbitrary Lagrangian-Eulerian (ALE) form to account for the effect of relative movement between the mesh and the fluid points. These equations  are solved in space-time using a semi-implicit fractional step approach and the finite element method (FEM). Free surface wave boundary effects are accounted for in the flow solution either by moving the free surface nodes in a Lagrangian manner, or else via the introduction of a prescribed pressure at the free surface computed from the wave height.

The movement of a fully or partially 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 in the fluid domain following the movement of the submerged body with minimum element distortion is presented. The mesh update procedure is based on the finite element solution of a linear elastic problem on the mesh domain, where fictitious elastic properties are assigned so that elements suffering a larger straining are stiffer <span id='citeF-14'></span>[[#cite-14|14]].

The content of the paper is structured as follows. First, details of the stabilized form of the governing equations for a viscous flow and the free surface using a finite  calculus procedure are given. The semi-implicit fractional step approach using the FEM is then described. Details of the computation of the stabilization parameters are also given. Next the mesh updating procedure is presented. Finally, the efficiency of the method proposed is shown in the 3D analysis of the standard square cavity problem and of several fluid-structure interaction problems with free surface waves.

==2 FINITE CALCULUS FORMULATION OF FLUID-FLOW ANDFREE SURFACE EQUATIONS==

The finite element solution of the incompressible Navier-Stokes equations with the classical Galerkin method may suffer from numerical instabilities from two main sources. The first is due to the advective-diffusive character of the equations which induces oscillations for high values of the velocity. The second source has to do with the mixed character of the equations which limits the choice of finite element interpolations for the velocity and pressure fields.

Solutions for these two problems have been extensively sought in the last years. Compatible velocity-pressure interpolations satisfying the inf-sup condition emanating from the second problem above mentioned have been used <span id='citeF-15'></span>[[#cite-15|15]]. In addition, the advective operator has been modified to include some “upwinding” effects <span id='citeF-18'></span>[[#cite-18|18]]. Recent procedures based on Galerkin Least Square (GLS) <span id='citeF-26'></span>[[#cite-26|26]], Characteristic Galerkin <span id='citeF-28'></span>[[#cite-28|28]], Variational Multiscale [30&#8211;32] and Residual Free Bubbles [33&#8211;35] techniques allow equal order interpolation for velocities and pressure by introducing a Laplacian of pressure term in the mass balance equation, while preserving the upwinding stabilization of the momentum equations. Most of these methods lack enough stability in the presence of sharp layers transversal to the velocity. This deficiency is usually corrected by adding new “shock capturing” stabilization terms to the already stabilized equations [36&#8211;37]. The computation of the stabilization parameters in all these methods is typically based on “ad hoc” generalizations of the parameters for the 1D linear advective-diffusive-reactive problem <span id='citeF-38'></span>[[#cite-38|38]].

Applications of stabilized GLS FEM to fluid-structure interaction problems, mainly of aerolastic type, have been reported in <span id='citeF-41'></span>[[#cite-41|41]].

Introduction of the free surface boundary condition in the flow equations increases considerably the difficulty of solving fluid-structure interation problems using FEM. A review of these difficulties and some solution procedures can be found in <span id='citeF-51'></span>[[#cite-51|51]]. Another successful application of stabilized FEM to free surface wave problems was reported in <span id='citeF-52'></span>[[#cite-52|52]].

This work presents a different approach for deriving stabilized finite element methods for incompressible flow problems with a free surface. The starting point is the stabilized form of the governing differential equations derived via a ''finite calculus'' (FIC) procedure. This technique first presented in <span id='citeF-6'></span>[[#cite-6|6]] is based on writting the different balance equations over a domain of finite size and retaining higher order terms. These terms incorporate the ingredients for the necessary stabilization of any transient and steady state numerical solution ''already at the  differential equations level''. Application of the standard Galerkin formulation to the consistently modified differential equations for the fluid flow problem leads to a stabilized system of discretized equations which overcomes ''the two problems'' above mentioned (i.e. the advective type instability  and that due to lack of compatibility between the velocity and pressure fields). Application of the FIC method to the free surface wave problem leads to a new stabilized governing equation for the free surface which again can be solved numerically by standard Galerkin FEM. In addition, the modified differential equations can be used to derive a numerical scheme for iteratively computing the stabilization parameters <span id='citeF-7'></span>[[#cite-7|7]].

Ilinca ''et al.'' <span id='citeF-40'></span>[[#cite-40|40]] have recently proposed a stabilized FEM for incompressible advective-diffusive transport and fluid flow problems based on applying the standard Galerkin technique to the modified governing differential equations obtained by expanding the residuals around a known finite element solution using Taylor series. The set of modified equations ressembles those obtained by the FIC method using a conceptually different procedure.

Initial applications of the FIC method to solve free surface ship wave problems were reported in <span id='citeF-1'></span>[[#cite-1|1]]. Idelsohn ''et al.'' <span id='citeF-51'></span>[[#cite-51|51]] have shown that starting from the stabilized FIC form of the free surface equation allows the identification of a number of stabilized upwinding finite difference schemes traditionally used for solving free surface problems in naval architecture.

The FIC formulation presented in this paper for incompressible flows with a free surface can be considered an extension of that recently developed in <span id='citeF-10'></span>[[#cite-10|10]] for finite element analysis of incompressible Navier-Stokes flows. A new formulation of the stabilized governing differential equations via the FIC method is here presented which holds for the viscous (Stokes) and zero viscosity (Euler) limit cases. The stabilized fluid flow equations are completed with  the FIC form of the free surface wave equation following the ideas first presented in <span id='citeF-2'></span>[[#cite-2|2]]. The set of stabilized governing equations is first discretized in time and then solved in space using a Galerkin finite element method.

A semi-implicit fractional step procedure is used for the momentum and mass balance equations allowing for equal order linear interpolations of the velocity and pressure variables over tetrahedral elements. Examples of application of the new stabilized finite element formulation to the standard square cavity flow problem and to a number of free surface ship-wave problems, including coupled fluid-structure interaction situations, are presented.

For the sake of preciseness, the basic ideas of the FIC method are given next.

===2.1 Basic concepts of the finite calculus (FIC) method===

Let us consider a sourceless transport problem over a one dimensional domain <math display="inline">\Omega </math>. The balance of flux <math display="inline">q</math> over a domain of finite size belonging to <math display="inline">\Omega </math> can be written as (Figure 1)

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>q_A - q_B=0  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (1)
|}

where <math display="inline">A</math> and <math display="inline">B</math> are the end points of the finite size domain of length <math display="inline">h</math>. In eq.(1) <math display="inline">q_A</math> and <math display="inline">q_B</math> represent the values of the flux <math display="inline">q</math> at points <math display="inline">A</math> and <math display="inline">B</math>, respectively.

<div id='img-1'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:draft_177_844678375-figura1.png|600px|Equilibrium of fluxes in a finite balance domain]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 1:''' Equilibrium of fluxes in a finite balance domain
|}

For instance, in an 1D advective-diffusive problem the flux <math display="inline">q=-cu\phi + k{d\phi \over dx}</math>, where <math display="inline">\phi </math> is the transported variable (i.e. the temperature in a thermal problem), <math display="inline">u</math> is the advective velocity and <math display="inline">c</math> and <math display="inline">k</math> are the advective and diffusive material parameters, respectively.

The flux <math display="inline">q_A</math> can be expressed in terms of the values at point <math display="inline">B</math> by the following Taylor series expansion

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>q_A=q_B-h {\partial q\over \partial x}\vert _B+{h^2\over 2}{d^2q\over dx^2}\vert _B + Oh^3  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (2)
|}

Substituting (2) into (1) gives after simplification and neglecting cubic terms in <math display="inline">h</math>

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>{dq\over dx}-\underline{{h\over 2} {d^2q\over dx^2}}{h\over 2} {d^2q\over dx^2}=0  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (3)
|}

where all terms are evaluated at the arbitrary point <math display="inline">B</math>.

Eq. (3) is the ''finite'' form of the balance equation over the domain <math display="inline">AB</math>. The underlined term in eq.(3) introduces the necessary stabilization for the discrete solution of eq.(3) using ''any'' numerical technique. Distance <math display="inline">h</math> is the characteristic length of the discrete problem and its value depends on the parameters of the discretization method chosen (such as the grid size <span id='citeF-6'></span>[[#cite-6|6]]). Note that for <math display="inline">h\to 0</math> the standard infinitesimal form of the balance equation <math display="inline">\left({dq\over dx}=0\right)</math> is recovered.

The above process can be extended to derive the stabilized balance differencial equations for any problem in fluid or solid mechanics as

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>r_d - \underline{{h_j\over 2} {\partial r_i\over \partial x_j}}{h_j\over 2} {\partial r_i\over \partial x_j}=0  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (4)
|}

where <math display="inline">r_i</math> is the standard form of the <math display="inline">i</math>th differential equation for the infinitesimal problem, <math display="inline">h_j</math> are the dimensions of the domain where balance of fluxes, forces, etc. is enforced, and <math display="inline">j=1,2,3</math> for 3D problems. It is important to note that the numerical solution of eq.(4) (together with the appropriate stabilized  boundary conditions) using Galerkin FE or central finite difference schemes leads to stable results <span id='citeF-6'></span>[[#cite-6|6]]. Details of the derivation of eq.(4) for steady-state and transient advective-diffusive and fluid flow problems can be found in <span id='citeF-6'></span>[[#cite-6|6]]. Applications of the FIC approach to the solution of these problems using Galerkin finite element and meshless procedures are reported in <span id='citeF-1'></span>[[#cite-1|1]].

The underlined stabilization terms in eqs.(3) and (4) are a consequence of accepting that the infinitesimal form of the balance equations is an unreachable limit within the framework of a discrete numerical solution. Indeed eqs.(3) or (4) are ''not longer valid'' for obtaining an analytical solution following traditional integration methods from infinitesimal calculus theory. The meaning of the new stabilized equations makes sense only in the context of a discrete numerical method  yielding approximate values of the solution at a finite set of points within the analysis domain.  Convergence to the ''exact '' analytical value at the points will occur only for the limit case of zero grid size (except for some simple 1D problems <span id='citeF-6'></span>[[#cite-6|6]]) which also implies naturally a zero value of the characteristic length parameters.

===2.2 FIC formulation of viscous flow and free surface equations===

We consider the motion around a body of a viscous incompressible fluid including a free surface.

The stabilized FIC form of the governing differential equations for the three dimensional  (3D) problem can be written in Arbitrary Lagrangian-Eulerian (ALE) form as <span id='citeF-2'></span>[[#cite-2|2]]

''Momentum''

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| 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} - \underline{{1\over 2} \delta {\partial r_{m_i}\over \partial t}}{1\over 2} \delta {\partial r_{m_i}\over \partial t} =0 \qquad \hbox{on } \Omega \quad i,j=1,2,3</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (5)
|}

''Mass balance''

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| 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;white-space: nowrap;" | (6)
|}

''Free surface''

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| 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}-\underline{{1\over 2}\gamma  {\partial r_\beta \over \partial t}}{1\over 2}\gamma  {\partial r_\beta \over \partial t} =0 \qquad \hbox{on } \Gamma _\beta \quad j=1,2</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (7)
|}

where

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>r_{m_i} =\rho \left[{\partial u_i\over \partial t} + v_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;white-space: nowrap;" | (8)
|-
| style="text-align: center;" | <math> r_d={\partial u_i\over \partial x_i} \qquad i=1,2,3 </math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (9)
|-
| style="text-align: center;" | <math> r_\beta ={\partial \beta \over \partial t} + v_i {\partial \beta \over \partial x_i}-v_3 \qquad i=1,2 </math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (10)
|-
| style="text-align: center;" | <math> \hbox{and}  v_i=u_i-u_i^m</math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (11)
|}
|}

Above, <math display="inline">u_i</math> is the velocity along the <math display="inline">i</math>th global reference axis, <math display="inline">u_i^m</math> is the velocity of the mesh nodes and <math display="inline">v_i</math> is the relative velocity between the moving mesh and the fluid point <math display="inline">i</math>, <math display="inline">\rho </math>  is the (constant) density of the fluid, <math display="inline">p</math> is the dynamic pressure defined as <math display="inline">p={1\over \rho } p_a - gx_3</math> where <math display="inline">p_a</math> is the absolute pressure and <math display="inline">x_3</math> is the vertical coordinate, <math display="inline">\beta </math> is the wave elevation (measured with respect to a reference flat surface) and <math display="inline">\tau _{ij}</math> are 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;width: 100%;" 
|-
| 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} {1\over 3} {\partial u_k\over \partial x_k}\right) </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (12)
|}

where <math display="inline">\delta _{ij}</math> is the Kronecker delta.

The boundary conditions for the stabilized problem are written as

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| 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;white-space: nowrap;" | (13)
|}

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

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 terms in eqs.(5)&#8211;(7)  introduce the necessary stabilization for the approximated numerical solution.

The ''characteristic length'' distances <math display="inline">h_{mj}</math> and <math display="inline">h_{dj}</math> represent the dimensions of the finite domain where balance of momentum and mass  is enforced. On the other hand, the characteristic distances <math display="inline">h_{\beta j}</math> in eq.(7) represent the dimensions of a finite domain surrounding a point where the velocity is constrained to be tangent to the free surface. The signs before the stabilization terms in eqs.(5)&#8211;(7) and (13) ensure a positive value of the characteristic length distances. The parameters <math display="inline">\delta </math> and <math display="inline">\gamma </math> in eqs.(5) and (7) have dimensions of time. Details of the derivation of eqs. (5)&#8211;(7) can be found in <span id='citeF-2'></span>[[#cite-2|2]]. As an example, the stabilized equation for the free surface (eq.(7)) is derived in the Appendix.

Eqs.(5&#8211;14) are the starting point for deriving a variety of stabilized numerical methods for solving the incompressible Navier-Stokes equations with a free surface. 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 <span id='citeF-6'></span>[[#cite-6|6]].

====Remark 1====

In reference <span id='citeF-10'></span>[[#cite-10|10]] a modified version of the Dirichlet condition (14) is used including an additional stabilization term. This term is not strictly necessary for the subsequent derivation and will be neglected here.

===2.3 Alternative form of the mass balance equation===

Starting from eq.(5) we can obtain after algegraic rearrangement

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>{\partial r_d\over \partial x_i}=\left({\mu \over 3}+{\rho u_i h_{m_i}\over 2}\right)^{-1} \bigg[\bar r_{m_i} -{h_{m_k}\over 2} {\partial r_{m_i}\over \partial x_k}+{\rho u_i h_{m_i}\over 2}{\partial r_d\over \partial x_i}-{\delta \over 2} {\partial r_{m_i}\over \partial t}\bigg] \quad \hbox{no sum in} \,\,i  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (15)
|}

where

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\bar r_{m_i}=r_{m_i}+ {\mu \over 3}{\partial r_d \over \partial x_i}  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (16)
|}

and <math display="inline">r_{m_i}</math> is given by eq.(8).

Inserting eq.(15) into eq.(6) gives

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>r_d - g_{ii} {\partial ^2p\over \partial x_i \partial x_i}+r_p=0  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (17)
|}

with

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>r_p=c_i \bar r_{m_i}-g_{ij}{\partial \over \partial x_j} \left(r_{m_i}-\delta _{ij} {\partial p\over \partial x_i}\right)+ {\rho u_i h_{m_i}\over 2}{\partial r_d\over \partial x_i}-{\delta \over 2} {\partial r_{m_i}\over \partial t}  \quad \hbox{no sum in} \,\,i  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (18)
|}

In above

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>g_{ij}=\left({4\mu \over 3 h_{d_i}h_{m_j}}+  {2\rho u_i h_{m_i}\over h_{d_i}h_{m_j}}\right)^{-1}\quad \hbox{no sum in}\,\, i  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (19)
|}

Note that for <math display="inline">h_{m_i}=h_{d_i}=h</math> where <math display="inline">h</math> is a typical grid dimension (i.e. the average element size), the value of <math display="inline">g_{ii}</math> is simply

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>g_{ii} =\left({4\mu \over 3h^2} + {2\rho u_i\over h}\right)^{-1}</math>
|}
|}

The stabilization parameter <math display="inline">g_{ii}</math> has now the form traditionally used in the GLS formulation for the viscous (Stokes) limit (<math display="inline">u_i=0</math>) and the inviscid (Euler) limit (<math display="inline">\mu =0</math>) and deduced from ad-hoc extensions of the 1D advective-diffusive problems <span id='citeF-18'></span>[[#cite-18|18]]. Note, however, that the general form of the stabilization parameter <math display="inline">g_{ij}</math> is deduced here from the general FIC formulation without further extrinsic assumptions.

Indeed, the precise computation of the characteristic length values is crucial for the practical application of above stabilized expressions. This problem is dealt with in a later section.

==3 FRACTIONAL STEP APPROACH==

The momentum equations (5) are first discretized in time using the following scheme

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>u_i^{n+1}=u_i^n-{\Delta t\over \rho } \bigg[\rho v_j^n{\partial u_i^n\over \partial x_j}+{\partial p^n\over \partial x_i}- {\partial \tau _{ij}^n \over \partial x_j}  - {h_{m_k}^n\over 2}{\partial r_{m_i}^n \over \partial x_k} - {\delta ^n\over 2} {\partial r_{m_i}^n \over \partial t}\bigg] </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (20)
|}

Eq.(20) is now split into the two following equations

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>u^*_i= u_i^n -{\Delta t\over \rho } \bigg[\rho v_j{\partial u_i\over \partial x_j} -  {\partial \tau _{ij}\over \partial x_j} -{h_{m_k}\over 2} {\partial r_{m_i}\over \partial x_k}-{\delta \over 2}{\partial r_{m_i}\over \partial t}\bigg]^n</math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (21)
|-
| style="text-align: center;" | <math> u_i^{n+1}= u^*_i - {\Delta t \over \rho } {\partial p^n\over \partial x_i}  </math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (22)
|}
|}

Note that the sum of eqs.(21) and (22) gives the original form of eq.(20).

Substituting eq.(22) into the stabilized mass balance equation (20) gives the standard Laplacian of pressure form

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\left({\Delta t \over \rho } + g_{ii}^n\right){\partial ^2 p^n\over \partial x_i \partial x_i}=r_d^*+r_p^n  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (23)
|}

where

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>r_d^*={\partial u_i^*\over \partial x_i}  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (24)
|}

Standard fractional step procedures neglect the contribution from the terms involving <math display="inline">g_{ii}</math> in eq. (23). These terms have an additional stabilization effect which improves the numerical solution when the values of <math display="inline">\Delta t</math> are small. Also the influence of the stabilization term <math display="inline">g_{ii}</math> has proven to be essential for obtaining a fully converged solution in steady state problems (see the square cavity example in a next section). Indeed accounting for this additional stabilization term has lead to improved numerical solutions in all problems solved. Similar conclusions have been reached in a recent work by Codina <span id='citeF-59'></span>[[#cite-59|59]].

Note that the cross-derivative terms <math display="inline">\displaystyle{\partial ^2 p \over \partial x_i \partial x_j}</math> have been kept within the term <math display="inline">r_p</math> in the r.h.s. of eq.(23). The influence of these terms should be studied in more detail in the future.

The stabilized free surface wave equation (7) is discretized in time to give

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

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. (21).

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

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

''Step 4''. Compute the new position of the free surface elevation <math display="inline">{\beta }^{n+1}</math> in the fluid domain by using  eq.(25).

''Step 5''. 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</math> and the viscous stresses <math display="inline">\tau ^{n+1}_{ij}</math>.

''Step 6''. Compute the new position of mesh nodes  in the fluid domain at time <math display="inline">n+1</math> by using the mesh update algorithm described in the next section. The updating process can also include the free surface nodes, although this is not strictly necessary.

Assuming air is at rest, the absolute pressure at the free surface at time <math display="inline">n+1</math> obtained from the stress equilibrium condition  (neglecting surface tension effects) as

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>p_a=\tau _{33}  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (26)
|}

The dynamic pressure at the free surface  is computed by

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>{p}^{n+1}=\left[{\tau _{33}\over \rho } - g \beta \right]^{n+1}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (27)
|}

where <math display="inline">g</math> is the gravity constant.

As already mentioned, the effect of changes in the free surface elevation are introduced in step 2 of the flow solution  as a prescribed dynamic pressure acting on the free surface. Note that eq.(27) allows to take into account the changes in the free surface without the need of updating the free surface nodes. A higher accuracy in the solution of the flow problem can however be obtained if the free surface nodes are updated after a number of time steps.

==4 FINITE ELEMENT DISCRETIZATION==

Spatial discretization is carried out using the finite element method <span id='citeF-15'></span>[[#cite-15|15]]. The stabilized formulation described allows an equal order interpolation of velocities and pressure <span id='citeF-10'></span>[[#cite-10|10]]. A linear interpolation over four node tetrahedra  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 triangles are chosen to interpolate <math display="inline">\beta </math> on the free surface mesh.  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;width: 100%;" 
|-
| style="text-align: center;" | <math>u_i =\sum \limits _j N_{j} \bar{(u_i)}_j</math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (28)
|-
| style="text-align: center;" | <math> p_i=\sum \limits _j N_{j}(\bar p_j)  </math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (29)
|}
|}

where <math display="inline">N_{j}</math> are the linear shape functions interpolating the velocity and pressure fields, respectively, and <math display="inline">\bar{(\cdot )}</math> denote nodal values <span id='citeF-15'></span>[[#cite-15|15]].

Similarly the wave height is discretized as

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

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

The discretized integral form in space is obtained by applying the standard Galerkin procedure to eqs.(21),(22),(23) and (25) and the boundary conditions (13). Solution of the discretized problem follows the pattern given below.

====Step 1. Solve for the nodal fractional velocities====

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\bar {u}^*={M}^{-1}{f}^n_1  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (31)
|}

with

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>{M}_{ij}=\int _\Omega N_iN_j d\Omega   </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (32)
|}

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>{f}^n = \left\{ \begin{array}{c}{f}_1\\ {f}_2\\ \vdots \\ {f}_n \end{array}\right\}^n\qquad , \qquad {f}^n_k = \left\{\begin{array}{c}f_{k_1}\\ f_{k_2}\\  f_{k_3} \end{array}\right\}^n  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (33)
|}

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>f_{k_i}^n=\int _\Omega N_{k} \bigg[u_i -{\Delta t\over \rho } \left(\rho v_j{\partial u_i\over \partial x_j} -{h_{m_k}\over 2} {\partial r_{m_i}\over \partial x_k}-{\delta \over 2}{\partial r_{m_i}\over \partial t}\right)\bigg]^nd\Omega +</math>
|-
| style="text-align: center;" | <math> +\int _\Omega {\Delta t\over \rho } {\partial N_{k}\over \partial x_j} \tau ^n_{ij} d\Omega -\int _{\Gamma _t} {\Delta t\over \rho } N_{k}t_i^n d\Gamma \quad ,\quad  i=1,2,3  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (34)
|}

The solution of eq.(31) can be speeded up by diagonalizing matrix <math display="inline">{M}</math>. Alternatively a simple Jacobi iteration procedure can be used and this has proved to converge in very few iterations.

No boundary condition is applied when solving for the fractional velocities <math display="inline">u_i^*</math> in eq.(31) as these velocities can be interpreted as a predicted value of the actual velocities. The kinematic boundary conditions (14) are applied in step 3 as shown below.

====Step 2. Solve for the nodal pressures at time n====

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>{H} \bar {p}^n={q}^n  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (35)
|}

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>H_{kl}=\int _\Omega {\partial N_{k} \over \partial x_i} \left({\Delta t\over \rho } + g_{ii}^n\right){\partial N_{l} \over \partial x_i}d\Omega </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (36)
|}

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>q^n_k= \int _\Omega {\partial N_{k} \over \partial x_i} u^*_i d\Omega - \int _\Omega N_{k} r_p^n d\Omega + \int _\Gamma N_{k} {\rho \over \Delta t} \left({\Delta t\over \rho } + g_{ii}^n\right)u^n_i n_i d\Gamma   </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (37)
|}

The last integral in eq.(37) can be neglected in solid walls and stationary free surfaces where the normal velocity is zero.

Recall that the dynamic pressures computed from step 6 are used as a boundary condition for solution of eq.(35).

====Step 3. Solve for the nodal velocities at time n+1====

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\bar {u}^{n+1}={M}^{-1}\bar {f}^n  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (38)
|}

where <math display="inline">{M}</math> is given by eq.(35) and

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\bar  f^{n+1}_{k_i}=\int _\Omega N_{k}\bigg[u_i^* - {\Delta t\over \rho }{\partial p^{n}\over \partial x_i}\bigg]d\Omega   </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (39)
|}

The kinematic boundary conditions on the nodal velocities (eq.(14)) are imposed when solving eq.(38).

====Step 4. Solve for the new free surface height at the time n+1====

The new free surface elevation <math display="inline">\beta ^{n+1}</math> in the fluid domain is computed as

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\bar{\boldsymbol \beta }^{n+1}={M}_\beta ^{-1} {s}^{n}  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (40)
|}

with

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>{M}_\beta =\int _{\Gamma _\beta } {N}_\beta ^T  {N}_\beta d\Gamma  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (41)
|}

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>s_i^{n+1}=\int _{\Gamma _\beta } N_{\beta _i} \bigg[\beta ^n - \Delta t \left(v_k^{n} {\partial \beta ^n\over \partial x_k}-v_3^{n}-{\gamma \over 2} {\partial r_\beta ^n\over \partial t}\right)\bigg]d\Gamma + \int _{\Gamma _\beta } {h_{\beta _j}\over 2}{\partial N_{\beta _i} \over \partial x_j}r_\beta d\Gamma  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (42)
|}

In the derivation of eq.(42) the assumption that <math display="inline">r_\beta =0</math>  at the boundary line of the free surface domain has been made.

Steps 5 and 6 follow the process described in the previous section.

==5 COMPUTATION OF THE STABILIZATION PARAMETERS==

Accurate evaluation of the stabilization parameters is one of the crucial issues in stabilized methods. Most of existing methods use expressions which are direct extensions of the values obtained for the simplest 1D case. It is also usual to accept the so called SUPG assumption, i.e. to admit that vector <math display="inline">{h}_m</math> has the direction of the velocity field <span id='citeF-6'></span>[[#cite-6|6]]. This unnecessary restriction leads to instabilities when sharp layers transversal to the velocity direction are present. This additional deficiency is usually corrected by adding a shock capturing or crosswind stabilization term [36&#8211;38].

Let us first assume for simplicity that the stabilization parameters for the mass balance equations are the same as those for the momentum equations. This implies

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>{h}_{m}={h}_{d}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (43)
|}

The problem remains now finding the value of the characteristic length vectors <math display="inline">{h}_{m}</math>. Indeed, the components of <math display="inline">{h}_m</math> can introduce the necessary stabilization along both the streamline and transversal directions to the flow.

Excellent results have been obtained in all problems solved  using linear  tetrahedra with the same value of the characteristic length vector for the three momentum equations defined by

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>{h}_{m_i}=h_s {{u}\over {u}}+h_{c} {{\boldsymbol \nabla } u\over \vert{\boldsymbol \nabla }u\vert }\qquad i=1,2,3   </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (44)
|}

where <math display="inline">u=\vert {u}\vert </math> and <math display="inline">h_s</math> and <math display="inline">h_{c}</math> are the “streamline” and “cross wind” contributions given by

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>h_s=\max ({l}^T_j {u})/{u} </math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (45)
|-
| style="text-align: center;" | <math> h_{c}=\max ({l}^T_j {\boldsymbol \nabla }u)/ \vert {\boldsymbol \nabla }u\vert \quad , \quad  j=1,n_s  </math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (46)
|}
|}

where <math display="inline">{l}_j</math> are the vectors defining the element sides (<math display="inline">n_s=6</math> for tetrahedra).

An alternative method for computing vector <math display="inline">{h}_m</math> in a more consistent manner is explained in <span id='citeF-60'></span>[[#cite-60|60]].

As for the free surface equation the following value of the characteristic length vector <math display="inline">{h}_\beta </math> has been taken

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>{h}_\beta =\bar h_s {{u}\over {u}}+\bar h_c {{\boldsymbol \nabla }\beta \over \vert {\boldsymbol \nabla }\beta \vert }   </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (47)
|}

The streamline parameter has been obtained by eq.(45) using the value of the velocity vector <math display="inline">u</math> over the 3 node triangles discretizing the free surface and <math display="inline">n_s=3</math>.

The cross wind parameter has been computed by

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\bar h_c = \max [{l}_j^T {\boldsymbol \nabla }\beta ] {1\over \vert {\boldsymbol \nabla }\beta \vert } \quad ,\quad j=1,2,3  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (48)
|}

Note that the cross-wind terms in eqs.(44) and (47) account for the effect of the gradient of the solution in the stabilization parameters. This is a standard assumption in most “shock-capturing” stabilization procedures [36&#8211;39].

Oñate and García have proposed a more consistent approach for computing the characteristic length vectors in terms of the residuals of the finite element solution. Details of this approach can be found in <span id='citeF-60'></span>[[#cite-60|60]].

Regarding the time stabilization parameters <math display="inline">\delta </math> and <math display="inline">\gamma </math> in eqs.(5) and (7) the value <math display="inline">\delta =\gamma =\Delta t</math> has been taken for solution of the problems presented in a next section.

==6 A SIMPLE ALGORITHM FOR UPDATING THE MESH NODES==

Different techniques have been proposed for dealing with mesh updating in fluid-structure interaction problems. The general aim of all methods is to prevent element distortion during mesh deformation <span id='citeF-41'></span>[[#cite-41|41]].

Chiandussi, Bugeda and Oñate <span id='citeF-14'></span>[[#cite-14|14]] have recently proposed a simple method for  the movement of mesh nodes ensuring minimum element distortion. The method is based on the iterative solution of a fictious 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;width: 100%;" 
|-
| 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;white-space: nowrap;" | (49)
|}

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 sought. The principal stresses are then given by

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| 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;white-space: nowrap;" | (50)
|}

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 <span id='citeF-14'></span>[[#cite-14|14]] is to equate the element strain energies in both analysis. Thus

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>U_1= {}^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>
| style="width: 5px;text-align: right;white-space: nowrap;" | (51)
|-
| style="text-align: center;" | <math> U_2= {}^2\sigma _i \varepsilon _i = 3E \bar  \varepsilon ^2  (1-2\nu )  </math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (52)
|}
|}

Equaling eqs.(51) and (52) 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;width: 100%;" 
|-
| 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;white-space: nowrap;" | (53)
|}

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.

''Step 2''. Compute the principal strains  and the values of the new Young modulus in each element using eq.(53) 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 distortion. Further details on this method including other alternatives for evaluating the Young modulus <math display="inline">E</math> can be found in <span id='citeF-14'></span>[[#cite-14|14]].

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 possible solution to this problem is to remesh the analysis domain. However, for most  problems, a mapping of the moving surfaces linked to mesh updating algorithm described above can avoid remeshing (Figure 2).

<div id='img-2'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Void_image.svg|600px|Changes in the fluid interface in a floating body]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 2:''' Changes in the fluid interface in a floating body
|}

The surface mapping technique used in this work is based on transforming the 3D curved surfaces into reference planes. This makes it possible 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 <span id='citeF-5'></span>[[#cite-5|5]].

==7 EXAMPLES==

All the examples shown next have been solved in a standard PC Pentium II 450 Mhz with a memory of 128 Mb.

===7.1 Example 1. Submerged NACA 0012 profile===

A 2D submerged NACA0012 profile at <math display="inline">\alpha=5^\circ </math> angle   of attack is studied. This configuration was tested experimentally by Duncan <span id='citeF-55'></span>[[#cite-55|55]]  for high Reynolds numbers (Re=400000) and modelled numerically using the   Euler equations   by several authors <span id='citeF-50'></span>[[#cite-50|50]]. The submerged depth of the airfoil is equal to  the chord and this was used as the length (L) for normalizing the problem.  The Froude number for all the cases tested was set to <math display="inline"> Fr={{u}\over \sqrt{{g L}}}= 0.5672 </math>  where <math display="inline">u</math> is the incoming flow velocity at infinity.

The stationary free   surface and the pressure distribution in the domain are shown in Figure 3.  The non-dimensional wave heights compare well with the experimental  results of <span id='citeF-55'></span>[[#cite-55|55]].

<div id='img-3'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:draft_177_844678375-fig3a.png|570px|]]
|[[Image:draft_177_844678375-fig3b.png|400px|Submerged NACA0012 profile. a) Detail of the mesh of 70000 linear tetrahedra chosen. b) Pressure contours. c) Stationary wave profile]]
|- style="text-align: center; font-size: 75%;"
| colspan="2" | '''Figure 3:''' Submerged NACA0012 profile. a) Detail of the mesh of 70000 linear tetrahedra chosen. b) Pressure contours. c) Stationary wave profile
|}

===7.2 Example 2. Sphere falling in a tube filled with liquid===

The movement of a sphere falling by gravity in a cylindrical tube filled with liquid is studied. The relationship between the diameters of the sphere and the tube is 1:4. The Reynolds number for the stationary speed is 100. The mesh has 85765 elements with 13946 nodes (Figure 4).

Figures 4 and 5 show the mesh deformation and contours of the mesh deformation and of the velocity in the domain for different times, respectively. The evolution of the falling speed is shown in Figure 5c. Note the good agreement with the so called Stokes velocity computed by equaling the weight of the sphere with the resistance  to the movement of the sphere expressed in terms of the velocity. Obviously, this value is slightly greater than the actual one as frictional effects are neglected.

A similar problem for a much greater number of spheres has been solved by Johnson and Tezduyar <span id='citeF-47'></span>[[#cite-47|47]].

<div id='img-4'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:draft_177_844678375-figura4.png|540px|Sphere falling in a tube filled with liquid. a) Geometry definition and detail of the mesh of 85765 linear tetrahedra chosen. b) Mesh deformation during the falling of the sphere]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 4:'''  Sphere falling in a tube filled with liquid. a) Geometry definition and detail of the mesh of 85765 linear tetrahedra chosen. b) Mesh deformation during the falling of the sphere
|}

<div id='img-5'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Void_image.svg|600px|Sphere falling in a tube filled with liquid a) Evolution of contours of the mesh deformation. b) Evolution of contours of velocity module. c) Evolution of falling speed. Straight line indicates the theoretical Stokes speed (1.195 m/s)]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 5:''' Sphere falling in a tube filled with liquid a) Evolution of contours of the mesh deformation. b) Evolution of contours of velocity module. c) Evolution of falling speed. Straight line indicates the theoretical Stokes speed (1.195 m/s)
|}

===7.3 Example 3. Wigley hull===

The last case considered here is the well known  Wigley Hull, given by   the analytical formula <math display="inline">y=0.5B(1 - 4x^2)(1 - z^2/D^2)</math> where <math display="inline">B</math> and <math display="inline">D</math>   are the beam and the draft of the ship hull at still water.

The same configuration was tested experimentally in <span id='citeF-57'></span>[[#cite-57|57]] and modelled   numerically by several authors <span id='citeF-50'></span>[[#cite-50|50]]. We use here an unstructured 3D   finite  element mesh of 65434 linear tetrahedra,   with a reference surface of  7800 triangles, partially   represented in Figure 6.

<div id='img-6'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:draft_177_844678375-figura6b.png|600px|]]
|[[Image:draft_177_844678375-figura6c.png|400px|Wigley hull. a) Pressure distribution and mesh deformation of the wigley hull (free model). b) Numerical and experimental body wave profiles.   c) Free surface contours for the truly free ship motion  ]]
|- style="text-align: center; font-size: 75%;"
| colspan="2" | '''Figure 6:''' Wigley hull. a) Pressure distribution and mesh deformation of the wigley hull (free model). b) Numerical and experimental body wave profiles.   c) Free surface contours for the truly free ship motion  
|}

Figure 6 also shows the results of the viscous  analysis of the Wigley model in three different cases (<math display="inline">L_{pp}=6m, F_n=0.316, \mu=10^{-3} Kg/m.s</math>). In the first case the volume mesh was considered fixed, not allowing free surface nor ship movements. Secondly, the volume mesh was updated due to free surface movement, considering the model fixed. The third case corresponds to the analysis of a real free model including the mesh updating due to free surface evaluation and ship movement (sinkage and trim). A Smagorinsky turbulence model was used in all the cases.

Table 1 shows the obtained total resistance coefficient in the three cases studied compared with the experimental data.

<br/>


{|  class="floating_tableSCP wikitable" style="text-align: center; margin: 1em auto;min-width:50%;"
|+ style="font-size: 75%;" |Table. 1 Wigley Hull. Total resistance coefficient
|- style="border-top: 2px solid;"
| style="text-align: left;border-left: 2px solid;border-right: 2px solid;" |  
| style="border-left: 2px solid;border-right: 2px solid;" | '''Experimental'''
| style="border-left: 2px solid;border-right: 2px solid;" | '''Numerical'''
|- style="border-top: 2px solid;"
| style="text-align: left;border-left: 2px solid;border-right: 2px solid;" |  Test 1 
| style="border-left: 2px solid;border-right: 2px solid;" | 5.2 <math display="inline">10^{-3}</math> 
| style="border-left: 2px solid;border-right: 2px solid;" | 4.9 <math display="inline">10^{-3}</math>
|- style="border-top: 2px solid;"
| style="text-align: left;border-left: 2px solid;border-right: 2px solid;" |   Test 2 
| style="border-left: 2px solid;border-right: 2px solid;" | 5.2 <math display="inline">10^{-3}</math> 
| style="border-left: 2px solid;border-right: 2px solid;" | 5.3 <math display="inline">10^{-3}</math>
|- style="border-top: 2px solid;border-bottom: 2px solid;"
| style="text-align: left;border-left: 2px solid;border-right: 2px solid;" |  Test 3 
| style="border-left: 2px solid;border-right: 2px solid;" | 4.9 <math display="inline">10^{-3}</math> 
| style="border-left: 2px solid;border-right: 2px solid;" | 5.1 <math display="inline">10^{-3}</math>

|}

In the study of the free model the numerical values of sinkage and trim were -0.1% and 0.035, respectively, while experiment gave -0.15% and 0.04.

Figure 6a shows the pressure distribution obtained near the Wigley hull for the free model. A number of streamlines have also been plotted in the figure. The obtained mesh deformation in this case is also presented in Figure 6b.

Comparisons of the obtained body wave profile with the experimental data for the free and fixed models are shown in Figure 6b. Significant differences are found close to stern in the case of the fixed model.

The free surface contours for the truly free ship motion are shown in Figure 6c.

==8 CONCLUSIONS==

The finite calculus method makes it possible to derive stabilized forms of the governing differential equations for a viscous fluid with a free surface. Solution of the new stabilized equations written in ALE form with a semi-implicit fractional step finite element method provides a straight-forward and stable algorithm for fluid-structure interaction analysis.

The mesh-moving scheme presented ensures minimum mesh distortion for large mesh displacements. The stabilized finite element method developed is adequate for solving large scale fluid-structure interaction problems in naval architecture and offshore engineering among others.

An academic version of the software developed using the formulation presented can be freely downloaded from www.cimne.upc.es/shyne.

==9 ACKNOWLEDGEMENTS==

Support for this work was provided by the European Community through projects Brite-Euram BR 967-4342 SHEAKS  and Esprit 24903 FLASH.

The  support of IZAR S.A.  is also  gratefully acknowledged.

Thanks are  given to Mr. J. Royo for his help in computing some of the examples presented.

The authors also thank Prof. S. Idelsohn, Prof. R. Lohner and Dr. C. Sacco for many useful discussions.

===BIBLIOGRAPHY===

<span style="text-align: center; font-size: 75%;">   <div id="cite-1"></div>
'''[[#citeF-1|[1]]]'''  J.  García, E. Oñate, H. Sierra, C. Sacco and S. Idelsohn, “A stabilized numerical method for analysis of ship hydrodynamics”, ECCOMAS CFD98, K. Papaliou et al. (Eds.), J. Wiley, 1998.   </span>

<div id="cite-2"></div>
'''[[#citeF-2|[2]]]'''  E. Oñate, S. Idelsohn, C. Sacco and J.  García, “Stabilization of the numerical solution for the free surface wave equation in fluid dynamics”, ECCOMAS CFD98, K. Papaliou ''et al.'' (Eds.), J. Wiley, 1998.

<div id="cite-3"></div>
'''[3]''' E. Oñate and J. García,A stabilized finite element method for analysis of fluid-structure interaction problems involving  surface waves, in ''Computational Methods for Fluid-Structure Interaction'', T. Kvamsdal ''et al.'', (Eds.), TAPIR Publishers, Norway, 1999.

<div id="cite-4"></div>
'''[4]'''  E. Oñate and J.  García, “A methodology for analysis of fluid-structure interaction accounting for free surface waves”, ''European Conference on Computational Mechanics (ECCM99)'', W. Wunderlich ''et al.'' (Eds.), August 31&#8211;September 3, Munich, Germany, 1999.

<div id="cite-5"></div>
'''[[#citeF-5|[5]]]'''  J.  García, ''A finite element method for analysis of naval structures'' (in Spanish), ''Ph.D. Thesis'', Univ. Politecnica de Catalunya, Barcelona, Spain, December 1999.

<div id="cite-6"></div>
'''[[#citeF-6|[6]]]''' E. Oñate,“Derivation of stabilized equations for   advective-diffusive transport and fluid flow problems”,  ''Comput. Methods Appl. Mech. Engrg.'',   Vol. 151, 1-2, pp. 233&#8211;267, 1998.

<div id="cite-7"></div>
'''[[#citeF-7|[7]]]''' E. Oñate, J. García 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.

<div id="cite-8"></div>
'''[8]''' E. Oñate, J.  García and S. Idelsohn,  An alpha-adaptive   approach for stabilized finite element solution of advective-diffusive   problems with sharp gradients, ''New Adv. in Adaptive   Comput. Methods in Mech.'', P. Ladeveze and J.T. Oden (Eds.),   Elsevier, 1998.

<div id="cite-9"></div>
'''[9]''' E. Oñate and M. Manzan, A general procedure for deriving stabilized space-time finite element methods for advective-diffusive problems, ''Int. J. Num. Meth. Fluids'', '''31''', 203&#8211;221, 1999.

<div id="cite-10"></div>
'''[[#citeF-10|[10]]]'''  E. Oñate, “A stabilized finite element method for incompressible viscous flows using a finite increment calculus formulation”, ''Comput. Methods Appl. Mech. Engrg.'', '''182''', 1&#8211;2, 355&#8211;370, 2000.

<div id="cite-11"></div>
'''[11]''' E. Oñate and M. Manzan, “Stabilization techniques for finite element analysis of convection diffusion problems”, in ''Computational Analysis of Heat Transfer'', G. Comini and B. Sunden (Eds.), WIT Press, Southampton, 2000.

<div id="cite-12"></div>
'''[12]''' 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.

<div id="cite-13"></div>
'''[13]'''  E. Oñate, C. Sacco and  S. Idelsohn, “A finite point method for incompressible flow problems”, ''Computing and Visualization in Science'', '''2''', 67&#8211;75, 2000.

<div id="cite-14"></div>
'''[[#citeF-14|[14]]]''' G. Chiandusi, G. Bugeda  and E. Oñate, “A simple method for update of finite element meshes”, ''Commun, Numer. Meth. Engng.'', '''16''', 1&#8211;9, 2000.

<div id="cite-15"></div>
'''[[#citeF-15|[15]]]''' O.C. Zienkiewicz  and  R.C. Taylor,''The finite element method'', 5th Edition, 3 Volumes, Butterworth&#8211;Heinemann, 2000.

<div id="cite-16"></div>
'''[16]''' R. Codina, ''A finite element model for incompressible flow   problems'', Ph.D. Thesis, Univ. Politécnica de Catalunya, Barcelona, Spain,   June 1992.

<div id="cite-17"></div>
'''[17]'''  M. Fortin and F. Thomasset, “Mixed finite element methods for incompressible flow problems”, ''J. Comput. Phys.'', '''31''', 113&#8211;145, 1979.

<div id="cite-18"></div>
'''[[#citeF-18|[18]]]''' A. Brooks and T.J.R. Hughes, “Streamline upwind/Petrov-Galerkin   formulation for convection dominated flows with particular emphasis   on the incompressible Navier-Stokes equations”, ''Comput. Methods Appl.   Mech. Engrg'', '''32''',  199&#8211;259, 1982.

<div id="cite-19"></div>
'''[19]''' T.J.R. Hughes and M. Mallet, “A new finite element formulations   for computational fluid dynamics: III. The generalized streamline operator   for multidimensional advective-diffusive systems”, ''Comput Methods Appl.   Mech. Engrg.'', '''58''', pp. 305&#8211;328,  1986.

<div id="cite-20"></div>
'''[20]''' P. Hansbo and a. Szepessy, “A velocity-pressure streamline diffusion finite element method for the incompressible Navier-Stokes equations”, ''Comput. Methods Appl. Mech. Engrg.'', '''84''', 175&#8211;192, 1990.

<div id="cite-21"></div>
'''[21]''' T.J.R. Hughes, L.P. Franca and M. Balestra, “A new finite element formulation for computational fluid dynamics. V Circumventing the Babuska-Brezzi condition: A stable Petrov-Galerkin formulation of the Stokes problem accomodating equal order interpolations”, ''Comput. Methods Appl. Mech. Engrg.'', '''59''', 85&#8211;89, 1986.

<div id="cite-22"></div>
'''[22]''' L.P. Franca and S.L. Frey, “Stabilized finite element methods:   II. The incompressible Navier-Stokes equations”, ''Comput. Method Appl.   Mech. Engrg.'', Vol. '''99''', pp. 209&#8211;233, 1992.

<div id="cite-23"></div>
'''[23]''' T.J.R. Hughes, G. Hauke and K. Jansen, “Stabilized finite   element methods in fluids: Inspirations, origins, status and recent   developments”, in: ''Recent Developments in Finite Element Analysis''.  A Book Dedicated to Robert L. Taylor,  T.J.R. Hughes, E. Oñate and   O.C. Zienkiewicz (Eds.), (International Center for Numerical   Methods in Engineering, Barcelona, Spain, pp. 272&#8211;292, 1994.

<div id="cite-24"></div>
'''[24]'''  M.A. Cruchaga and E. Oñate, “A finite element formulation for incompressible flow problems using a generalized streamline operator”, ''Comput. Methods in Appl. Mech. Engrg.'', '''143''', 49&#8211;67, 1997.

<div id="cite-25"></div>
'''[25]'''   M.A. Cruchaga and E. Oñate, “A generalized streamline finite element approach for the analysis of incompressible flow problems including moving surfaces”, ''Comput. Methods in Appl. Mech. Engrg.'', '''173''', 241&#8211;255, 1999.

<div id="cite-26"></div>
'''[[#citeF-26|[26]]]''' T.J.R. Hughes, L.P. Franca and G.M. Hulbert, “A new finite   element formulation for computational fluid dynamics: VIII. The   Galerkin/least-squares method for advective-diffusive equations”, ''  Comput. Methods Appl. Mech. Engrg.'', '''73''',  pp. 173&#8211;189, 1989.

<div id="cite-27"></div>
'''[27]''' T.E. Tezduyar, S. Mittal, S.E. Ray and R. Shih, “Incompressible flow computations with stabilized bilinear and linear equal order interpolation velocity&#8211;pressure elements”, ''Comput. Methods Appl. Mech. Engrg.'', '''95''', 221&#8211;242, 1992.

<div id="cite-28"></div>
'''[[#citeF-28|[28]]]''' O.C. Zienkiewicz  and  R. Codina, “A general algorithm for   compressible and incompressible flow. Part I: The split characteristic   based scheme”, ''Int. J. Num. Meth. in Fluids'', '''20''', 869-85, (1995).

<div id="cite-29"></div>
'''[29]''' O.C. Zienkiewicz, K. Morgan, B.V.K. Satya Sai, R. Codina and M.   Vázquez, “A general algorithm for compressible and incompressible flow. Part II: Tests on the explicit form”, ''Int. J. Num. Meth. in Fluids'', '''  20''', No. 8-9, 886-913, 1995.

<div id="cite-30"></div>
'''[30]''' T.J.R. Hughes, “Multiscale phenomena: Green functions,   subgrid scale models, bubbles and the origins of stabilized methods”,    ''Comput. Methods Appl. Mech. Engrg'', Vol. '''127''', pp. 387&#8211;401, 1995.

<div id="cite-31"></div>
'''[31]'''  R. Codina, “A stabilized finite element method for generalized stationary incompressible flows”, Publication PI-148, CIMNE, Barcelona, February 1999.

<div id="cite-32"></div>
'''[32]'''  R. Codina and J. Blasco, “Stabilized finite element method for the transient Navier-Stokes equations based on a pressure gradient operator”,  ''Comput. Methods in Appl. Mech. Engrg.'', '''182''', 277&#8211;301, 2000.

<div id="cite-33"></div>
'''[33]''' F. Brezzi, M.O. Bristeau, L.P. Franca, M. Mallet and G. Rogé, “A   relationship between stabilized finite element methods and the Galerkin   method with bubble functions”, ''Comput. Methods Appl. Mech. Engrg.'', Vol.   '''96''', pp. 117&#8211;129, 1992.

<div id="cite-34"></div>
'''[34]''' F. Brezzi, D. Marini and A. Russo, “Pseudo residual-free bubbles   and stabilized methods”,  ''Computational  Methods   in Applied  Sciences '96'',  J. Periaux ''et. al.'' (Eds.), J. Wiley, 1996.

<div id="cite-35"></div>
'''[35]''' F. Brezzi, L.P. Franca, T.J.R. Hughes and A. Russo, ``<math>b=\int   g</math>'', ''Comput. Methods Appl. Mech. Engrg.'', '''145''', 329&#8211;339, 1997.

<div id="cite-36"></div>
'''[36]'''  T.J.R. Hughes and M. Mallet, ``A new finite element formulations for computational fluid dynamics: IV. A discontinuity capturing operator for multidimensional advective-diffusive system, ''Comput. Methods Appl. Mech. Engrg.'', '''58''', 329&#8211;336, 1986.

<div id="cite-37"></div>
'''[37]'''  R. Codina, “A discontinuity-capturing crosswind dissipation for the finite element solution of the convection-diffusion equation”, ''Comput.  Methods Appl. Mech. Engrg.'', '''110''', 325&#8211;342, 1993.

<div id="cite-38"></div>
'''[[#citeF-38|[38]]]'''  R. Codina, “Comparison of some finite element methods for solving the diffusion-convection-reaction equation”, ''Comput. Methods Appl. Mech. Engrg.'', '''156''', 185&#8211;210, 1998.

<div id="cite-39"></div>
'''[39]'''  R. Codina, “On stabilized finite element methods for linear systems of convection-diffusion-reaction equation”, ''Comput. Methods Appl. Mech. Engrg.'', '''188''', 61&#8211;83, 2000.

<div id="cite-40"></div>
'''[[#citeF-40|[40]]]'''  F. Ilinca, J.F. Hétu and D. Pelletier, “On stabilized finite element formulation for incompressible advective-diffusive transport and fluid flow problems”, ''Comp. Methods Appl. Mech. Engrg.'', '''188''', 235&#8211;257, 2000.

<div id="cite-41"></div>
'''[[#citeF-41|[41]]]'''  T.E. Tezduyar, M. Behr and J. Liou, “A new strategy for finite element computations involving moving boundaries and interfaces - the deforming-spatial-domain/space-time procedure: I. The concept and the preliminary tests”, ''Comput. Methods in Appl. Mech. and Engrg.'', '''94''', 339&#8211;351, 1992.

<div id="cite-42"></div>
'''[42]'''  T.E. Tezduyar, M. Behr, S. Mittal and J. Liou, “A new strategy for finite element computations involving moving boundaries and interfaces - the deforming-spatial-domain/space-time procedure: II. Computation of free-surface flows, two liquid flows and flows with drifting cylinders”, ''Comput. Methods Appl. Mech. and Engrg.'', '''94''', 353&#8211;371, 1992.

<div id="cite-43"></div>
'''[43]'''  S. Mittal and T.E. Tezduyar, “Massive parallel finite element computation of incompressible flows involving fluid-body interaction”, ''Comput. Methods Appl. Mech. Engrg.'', '''112''', 253&#8211;282, 1994.

<div id="cite-44"></div>
'''[44]'''  S. Mittal and T.E. Tezduyar, “Parallel finite element simulation of 3D incompressible flows - fluid structure interactions”, ''Int. J. Num. Meth. Fluids'', '''21''', 933&#8211;953, 1995.

<div id="cite-45"></div>
'''[45]'''  V. Kalro, S. Aliabadi, W. Garrard, T. Tezduyar, S. Mittal and K. Stein, “Parallel finite element simulation of large Ram-air parachutes”, ''Int. J. Num. Meth. Fluids'', '''24''', 1353&#8211;1369, 1997.

<div id="cite-46"></div>
'''[46]'''  N. Maman and C. Farhat, “Matching fluid and structure meshes for aeroelastic computations: A parallel approach”, ''Computers and Structures'', '''54''', 779&#8211;785, 1995.

<div id="cite-47"></div>
'''[[#citeF-47|[47]]]'''  A.A. Johnson and T.E. Tezduyar, “3D simulation of fluid-particle interaction with the number of particles reaching 100”, ''Comput. Methods Appl. Mech. Engrg.'', '''145''', 301&#8211;321, 1997.

<div id="cite-48"></div>
'''[48]'''  N.A. Wall, M. Bischoff and E. Ramm, “Stabilization techniques for fluid and structural finite elements”, in ''Computational Mechanics. New Trends and Applications'', S.R. Idelsohn, E. Oñate and E.N. Dvorkin (Eds.), CIMNE/IACM, 1998.

<div id="cite-49"></div>
'''[49]'''  M. Storti, J. d'Elia and S.R. Idelsohn, “Algebraic discrete non-local (DNL) absorbing boundary condition for the ship wave resistance problem”, ''J. Computational Physics'', Vol. '''146''', n.2, 570&#8211;602, 1998.

<div id="cite-50"></div>
'''[[#citeF-50|[50]]]'''  M. Storti, J. d'Elia and S.R. Idelsohn, “Computing ship wave resistance form wave amplitude with the DNL absorbing boundary condition”, ''Comun. in Numer. Meth. in Engng'', Vol. '''14''', 997&#8211;1012, 1998.

<div id="cite-51"></div>
'''[[#citeF-51|[51]]]'''  S. Idelsohn, E. Oñate and C. Sacco, “Finite element solution of free surface ship-wave problem”, ''Int. J. Num. Meth. Engng.'', '''45''', 503&#8211;508, 1999.

<div id="cite-52"></div>
'''[[#citeF-52|[52]]]'''  R. Löhner, C. Yang, E. Oñate and S. Idelsohn, “An unstructured grid-based parallel free surface solver”, ''Appl. Num. Math.'', '''31''', 271&#8211;293, 1999.

<div id="cite-53"></div>
'''[53]'''  O.C. Zienkiewicz, J.Z. Zhu, “The Superconvergent patch recovery (SPR) and adaptive finite element refinement”, ''Comput. Methods Appl. Mech. Engrg.'', '''101''', 207&#8211;224, 1992.

<div id="cite-54"></div>
'''[54]'''  N.W. Wiberg, F. Abdulwahab, X.D. Li, “Error estimation  and adaptive procedures based on superconvergent patch recovery”, ''Archives Comput. Meth. Engng.'', '''4''' (3), 203&#8211;242, 1997.

<div id="cite-55"></div>
'''[[#citeF-55|[55]]]'''  J.H. Duncan, “The breaking and non-breaking wave resistance of a two-dimensional hydrofoil”, ''J. Fluid Mech.'', Vol. '''126''', 1983.

<div id="cite-56"></div>
'''[56]'''  T. Hino, L. Martinelli and A. Jameson, “A finite volument method with unstructured grid for free surface flow”, pp. 173-194 in ''Proc. of the 6th Int. Conf. Num. Ship Hydrodynamics'', Iowa City, Iowa, 1993.

<div id="cite-57"></div>
'''[[#citeF-57|[57]]]'''  Procedings 2nd DTNSRDC Workshop on Ship Wave Resistance Computations. David Taylor Naval Ship Research and Development Center. Noblese, F. and McCarthy J.H. (Eds.), Maryland, USA, 1983.

<div id="cite-58"></div>
'''[58]'''  J.R. Farmer, L. Martinelli and A. Jameson, “A fast multigrid method for solving incompressible hydrodynamic problems with free surfaces”, ''AIAA J.'', '''32''', 6, 1175-1182, 1993.

<div id="cite-59"></div>
'''[[#citeF-59|[59]]]'''  R. Codina, “Pressure stability in fractional step finite element methods for incompressible flows”. Submitted to ''Comput. Methods Appl. Mech. Engrg.'', 2000.

<div id="cite-60"></div>
'''[[#citeF-60|[60]]]'''  E. Oñate and J. García, “A finite element method for fluid-structure interaction with surface waves using a finite calculus formulation”. To appear in ''Comput. Methods Appl. Mech. Engrg.''