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–mail: onate@cimne.upc.es

Web page: 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:

1. The difficulty of solving numerically the incompressible fluid dynamic equations which typically include intrinsic nonlinearities except for the simplest and limited potential flow model.
2. The obstacles in solving the constraint equation stating that the fluid particles remain on the free surface boundary which position is in turn unknown.
3. 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.

This paper extends recent work of the authors 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 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 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 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 15. In addition, the advective operator has been modified to include some “upwinding” effects 18. Recent procedures based on Galerkin Least Square (GLS) 26, Characteristic Galerkin 28, Variational Multiscale [30–32] and Residual Free Bubbles [33–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–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 38.

Applications of stabilized GLS FEM to fluid-structure interaction problems, mainly of aerolastic type, have been reported in 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 51. Another successful application of stabilized FEM to free surface wave problems was reported in 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 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 7.

Ilinca et al. 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 1. Idelsohn et al. 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 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 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 ${\textstyle \Omega }$. The balance of flux ${\textstyle q}$ over a domain of finite size belonging to ${\textstyle \Omega }$ can be written as (Figure 1)

$q_{A}-q_{B}=0$

(1)

where ${\textstyle A}$ and ${\textstyle B}$ are the end points of the finite size domain of length ${\textstyle h}$. In eq.(1) ${\textstyle q_{A}}$ and ${\textstyle q_{B}}$ represent the values of the flux ${\textstyle q}$ at points ${\textstyle A}$ and ${\textstyle B}$, respectively.

Equilibrium of fluxes in a finite balance domain

Figure 1: Equilibrium of fluxes in a finite balance domain

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

The flux ${\textstyle q_{A}}$ can be expressed in terms of the values at point ${\textstyle B}$ by the following Taylor series expansion

$q_{A}=q_{B}-h{\partial q \over \partial x}\vert _{B}+{h^{2} \over 2}{d^{2}q \over dx^{2}}\vert _{B}+Oh^{3}$

(2)

Substituting (2) into (1) gives after simplification and neglecting cubic terms in ${\textstyle h}$

${dq \over dx}-{\underline {{h \over 2}{d^{2}q \over dx^{2}}}}{h \over 2}{d^{2}q \over dx^{2}}=0$

(3)

where all terms are evaluated at the arbitrary point ${\textstyle B}$.

Eq. (3) is the finite form of the balance equation over the domain ${\textstyle AB}$. The underlined term in eq.(3) introduces the necessary stabilization for the discrete solution of eq.(3) using any numerical technique. Distance ${\textstyle h}$ 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 6). Note that for ${\textstyle h\to 0}$ the standard infinitesimal form of the balance equation ${\textstyle \left({dq \over dx}=0\right)}$ is recovered.

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

$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$

(4)

where ${\textstyle r_{i}}$ is the standard form of the ${\textstyle i}$th differential equation for the infinitesimal problem, ${\textstyle h_{j}}$ are the dimensions of the domain where balance of fluxes, forces, etc. is enforced, and ${\textstyle j=1,2,3}$ 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 6. Details of the derivation of eq.(4) for steady-state and transient advective-diffusive and fluid flow problems can be found in 6. Applications of the FIC approach to the solution of these problems using Galerkin finite element and meshless procedures are reported in 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 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 2

Momentum

$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$

(5)

Mass balance

$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$

(6)

Free surface

$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$

(7)

where

$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}}$ (8) $r_{d}={\partial u_{i} \over \partial x_{i}}\qquad i=1,2,3$ (9) $r_{\beta }={\partial \beta \over \partial t}+v_{i}{\partial \beta \over \partial x_{i}}-v_{3}\qquad i=1,2$ (10) ${\hbox{and}}v_{i}=u_{i}-u_{i}^{m}$ (11)

Above, ${\textstyle u_{i}}$ is the velocity along the ${\textstyle i}$th global reference axis, ${\textstyle u_{i}^{m}}$ is the velocity of the mesh nodes and ${\textstyle v_{i}}$ is the relative velocity between the moving mesh and the fluid point ${\textstyle i}$, ${\textstyle \rho }$ is the (constant) density of the fluid, ${\textstyle p}$ is the dynamic pressure defined as ${\textstyle p={1 \over \rho }p_{a}-gx_{3}}$ where ${\textstyle p_{a}}$ is the absolute pressure and ${\textstyle x_{3}}$ is the vertical coordinate, ${\textstyle \beta }$ is the wave elevation (measured with respect to a reference flat surface) and ${\textstyle \tau _{ij}}$ are the viscous stresses related to the viscosity ${\textstyle \mu }$ by the standard expression

$\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)$

(12)

where ${\textstyle \delta _{ij}}$ is the Kronecker delta.

The boundary conditions for the stabilized problem are written as

$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}$

(13)

$u_{j}-u_{j}^{p}=0\qquad {\hbox{on }}\Gamma _{u}$

(14)

where ${\textstyle n_{j}}$ are the components of the unit normal vector to the boundary and ${\textstyle t_{i}}$ and ${\textstyle u_{j}^{p}}$ are prescribed tractions and displacements on the boundaries ${\textstyle \Gamma _{t}}$ and ${\textstyle \Gamma _{u}}$, respectively.

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

The characteristic length distances ${\textstyle h_{mj}}$ and ${\textstyle h_{dj}}$ represent the dimensions of the finite domain where balance of momentum and mass is enforced. On the other hand, the characteristic distances ${\textstyle h_{\beta j}}$ 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)–(7) and (13) ensure a positive value of the characteristic length distances. The parameters ${\textstyle \delta }$ and ${\textstyle \gamma }$ in eqs.(5) and (7) have dimensions of time. Details of the derivation of eqs. (5)–(7) can be found in 2. As an example, the stabilized equation for the free surface (eq.(7)) is derived in the Appendix.

Eqs.(5–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 6.

Remark 1

In reference 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

${\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$

(15)

where

${\bar {r}}_{m_{i}}=r_{m_{i}}+{\mu \over 3}{\partial r_{d} \over \partial x_{i}}$

(16)

and ${\textstyle r_{m_{i}}}$ is given by eq.(8).

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

$r_{d}-g_{ii}{\partial ^{2}p \over \partial x_{i}\partial x_{i}}+r_{p}=0$

(17)

with

$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$

(18)

In above

$g_{ij}=\left({4\mu \over 3h_{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$

(19)

Note that for ${\textstyle h_{m_{i}}=h_{d_{i}}=h}$ where ${\textstyle h}$ is a typical grid dimension (i.e. the average element size), the value of ${\textstyle g_{ii}}$ is simply

$g_{ii}=\left({4\mu \over 3h^{2}}+{2\rho u_{i} \over h}\right)^{-1}$

The stabilization parameter ${\textstyle g_{ii}}$ has now the form traditionally used in the GLS formulation for the viscous (Stokes) limit (${\textstyle u_{i}=0}$) and the inviscid (Euler) limit (${\textstyle \mu =0}$) and deduced from ad-hoc extensions of the 1D advective-diffusive problems 18. Note, however, that the general form of the stabilization parameter ${\textstyle g_{ij}}$ 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

$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 ]}$

(20)

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

$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}$ (21) $u_{i}^{n+1}=u_{i}^{*}-{\Delta t \over \rho }{\partial p^{n} \over \partial x_{i}}$ (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

$\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}$

(23)

where

$r_{d}^{*}={\partial u_{i}^{*} \over \partial x_{i}}$

(24)

Standard fractional step procedures neglect the contribution from the terms involving ${\textstyle g_{ii}}$ in eq. (23). These terms have an additional stabilization effect which improves the numerical solution when the values of ${\textstyle \Delta t}$ are small. Also the influence of the stabilization term ${\textstyle g_{ii}}$ 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 59.

Note that the cross-derivative terms ${\textstyle \displaystyle {\partial ^{2}p \over \partial x_{i}\partial x_{j}}}$ have been kept within the term ${\textstyle r_{p}}$ 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

$\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_{\beta }^{n} \over \partial x_{j}}-{\gamma \over 2}{\partial r_{\beta }^{n} \over \partial t}\right]\qquad i,j=1,2$

(25)

A typical solution in time includes the following steps.

Step 1. Solve explicitely for the so called fractional velocities ${\textstyle u_{i}^{*}}$ using eq. (21).

Step 2. Solve for the dynamic pressure field ${\textstyle p^{n}}$ 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 ${\textstyle u_{i}^{n+1}}$ at the updated configuration for each mesh node using eq.(22)

Step 4. Compute the new position of the free surface elevation ${\textstyle {\beta }^{n+1}}$ 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 ${\textstyle p^{n}}$ and the viscous stresses ${\textstyle \tau _{ij}^{n+1}}$.

Step 6. Compute the new position of mesh nodes in the fluid domain at time ${\textstyle n+1}$ 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 ${\textstyle n+1}$ obtained from the stress equilibrium condition (neglecting surface tension effects) as

$p_{a}=\tau _{33}$

(26)

The dynamic pressure at the free surface is computed by

${p}^{n+1}=\left[{\tau _{33} \over \rho }-g\beta \right]^{n+1}$

(27)

where ${\textstyle g}$ 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 15. The stabilized formulation described allows an equal order interpolation of velocities and pressure 10. A linear interpolation over four node tetrahedra for both ${\textstyle u_{i}}$ and ${\textstyle p}$ is chosen in the examples shown in the paper. Similarly, linear triangles are chosen to interpolate ${\textstyle \beta }$ on the free surface mesh. The velocity and pressure fields are interpolated within each element in the standard finite element manner as

$u_{i}=\sum \limits _{j}N_{j}{\bar {(u_{i})}}_{j}$ (28) $p_{i}=\sum \limits _{j}N_{j}({\bar {p}}_{j})$ (29)

where ${\textstyle N_{j}}$ are the linear shape functions interpolating the velocity and pressure fields, respectively, and ${\textstyle {\bar {(\cdot )}}}$ denote nodal values 15.

Similarly the wave height is discretized as

$\beta =\sum \limits _{j}N_{\beta _{j}}({\bar {\beta }}_{j})$

(30)

where ${\textstyle N_{\beta _{j}}}$ 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

${\bar {u}}^{*}={M}^{-1}{f}_{1}^{n}$

(31)

with

${M}_{ij}=\int _{\Omega }N_{i}N_{j}d\Omega$

(32)

${f}^{n}=\left\{{\begin{array}{c}{f}_{1}\\{f}_{2}\\\vdots \\{f}_{n}\end{array}}\right\}^{n}\qquad ,\qquad {f}_{k}^{n}=\left\{{\begin{array}{c}f_{k_{1}}\\f_{k_{2}}\\f_{k_{3}}\end{array}}\right\}^{n}$

(33)

$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 ]}^{n}d\Omega +$ $+\int _{\Omega }{\Delta t \over \rho }{\partial N_{k} \over \partial x_{j}}\tau _{ij}^{n}d\Omega -\int _{\Gamma _{t}}{\Delta t \over \rho }N_{k}t_{i}^{n}d\Gamma \quad ,\quad i=1,2,3$

(34)

The solution of eq.(31) can be speeded up by diagonalizing matrix ${\textstyle {M}}$. 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 ${\textstyle u_{i}^{*}}$ 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

${H}{\bar {p}}^{n}={q}^{n}$

(35)

$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$

(36)

$q_{k}^{n}=\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_{i}^{n}n_{i}d\Gamma$

(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

${\bar {u}}^{n+1}={M}^{-1}{\bar {f}}^{n}$

(38)

where ${\textstyle {M}}$ is given by eq.(35) and

${\bar {f}}_{k_{i}}^{n+1}=\int _{\Omega }N_{k}{\bigg [}u_{i}^{*}-{\Delta t \over \rho }{\partial p^{n} \over \partial x_{i}}{\bigg ]}d\Omega$

(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 ${\textstyle \beta ^{n+1}}$ in the fluid domain is computed as

${\bar {\boldsymbol {\beta }}}^{n+1}={M}_{\beta }^{-1}{s}^{n}$

(40)

with

${M}_{\beta }=\int _{\Gamma _{\beta }}{N}_{\beta }^{T}{N}_{\beta }d\Gamma$

(41)

$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$

(42)

In the derivation of eq.(42) the assumption that ${\textstyle r_{\beta }=0}$ 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 ${\textstyle {h}_{m}}$ has the direction of the velocity field 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–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

${h}_{m}={h}_{d}$

(43)

The problem remains now finding the value of the characteristic length vectors ${\textstyle {h}_{m}}$. Indeed, the components of ${\textstyle {h}_{m}}$ 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

${h}_{m_{i}}=h_{s}{{u} \over {u}}+h_{c}{{\boldsymbol {\nabla }}u \over \vert {\boldsymbol {\nabla }}u\vert }\qquad i=1,2,3$

(44)

where ${\textstyle u=\vert {u}\vert }$ and ${\textstyle h_{s}}$ and ${\textstyle h_{c}}$ are the “streamline” and “cross wind” contributions given by

$h_{s}=\max({l}_{j}^{T}{u})/{u}$ (45) $h_{c}=\max({l}_{j}^{T}{\boldsymbol {\nabla }}u)/\vert {\boldsymbol {\nabla }}u\vert \quad ,\quad j=1,n_{s}$ (46)

where ${\textstyle {l}_{j}}$ are the vectors defining the element sides (${\textstyle n_{s}=6}$ for tetrahedra).

An alternative method for computing vector ${\textstyle {h}_{m}}$ in a more consistent manner is explained in 60.

As for the free surface equation the following value of the characteristic length vector ${\textstyle {h}_{\beta }}$ has been taken

${h}_{\beta }={\bar {h}}_{s}{{u} \over {u}}+{\bar {h}}_{c}{{\boldsymbol {\nabla }}\beta \over \vert {\boldsymbol {\nabla }}\beta \vert }$

(47)

The streamline parameter has been obtained by eq.(45) using the value of the velocity vector ${\textstyle u}$ over the 3 node triangles discretizing the free surface and ${\textstyle n_{s}=3}$.

The cross wind parameter has been computed by

${\bar {h}}_{c}=\max[{l}_{j}^{T}{\boldsymbol {\nabla }}\beta ]{1 \over \vert {\boldsymbol {\nabla }}\beta \vert }\quad ,\quad j=1,2,3$

(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–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 60.

Regarding the time stabilization parameters ${\textstyle \delta }$ and ${\textstyle \gamma }$ in eqs.(5) and (7) the value ${\textstyle \delta =\gamma =\Delta t}$ 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 41.

Chiandussi, Bugeda and Oñate 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 ${\textstyle {\bar {E}}}$ and the Poisson coefficient ${\textstyle \nu }$. Once a discretized finite element problem has been solved using, for instance, standard ${\textstyle C^{o}}$ linear triangles (in 2D) or linear tetraedra (in 3D), the principal stresses ${\textstyle {}^{1}\sigma _{i}}$ at the center of each element are obtained as

${}^{1}\sigma _{i}={\bar {E}}[\varepsilon _{i}-\nu (\varepsilon _{j}+\varepsilon _{k})]\qquad i,j=1,2,3~~{\hbox{for 3D }}$

(49)

where ${\textstyle \varepsilon _{i}}$ are the principal strains.

Let us assume now that a uniform strain field ${\textstyle \varepsilon _{i}={\bar {\varepsilon }}}$ throughout the mesh is sought. The principal stresses are then given by

${}^{2}\sigma _{i}=E{\bar {\varepsilon }}(1-2\nu )\qquad i=1,2,3~~{\hbox{for 3D }}$

(50)

where ${\textstyle E}$ is the unknown Young modulus for the element.

A number of criteria can be now used to find the value of ${\textstyle E}$. The most effective approach found in 14 is to equate the element strain energies in both analysis. Thus

$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})]$ (51) $U_{2}={}^{2}\sigma _{i}\varepsilon _{i}=3E{\bar {\varepsilon }}^{2}(1-2\nu )$ (52)

Equaling eqs.(51) and (52) gives the sought Young modulus ${\textstyle E}$ as

$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})]$

(53)

Note that the element Young modulus is proportional to the element deformation as desired. Also recall that both ${\textstyle {\bar {E}}}$ and ${\textstyle {\bar {\varepsilon }}}$ 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 ${\textstyle {\bar {E}}}$ and ${\textstyle \nu }$. 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 ${\textstyle {\bar {\varepsilon }}}$. Repeat the finite element solution of the linear elastic problem with prescribed boundary displacements using the new values of ${\textstyle E}$ 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 ${\textstyle E}$ can be found in 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).

Changes in the fluid interface in a floating body

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 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 ${\textstyle \alpha =5^{\circ }}$ angle of attack is studied. This configuration was tested experimentally by Duncan 55 for high Reynolds numbers (Re=400000) and modelled numerically using the Euler equations by several authors 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 ${\textstyle Fr={{u} \over {\sqrt {gL}}}=0.5672}$ where ${\textstyle u}$ 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 55.

File:Draft 177 703272400-fig3a.png	Submerged NACA0012 profile. a) Detail of the mesh of 70000 linear tetrahedra chosen. b) Pressure contours. c) Stationary wave profile
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 47.

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

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

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)

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 ${\textstyle y=0.5B(1-4x^{2})(1-z^{2}/D^{2})}$ where ${\textstyle B}$ and ${\textstyle D}$ are the beam and the draft of the ship hull at still water.

The same configuration was tested experimentally in 57 and modelled numerically by several authors 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.

File:Draft 177 703272400-figura6b.png	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: 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 (${\textstyle L_{pp}=6m,F_{n}=0.316,\mu =10^{-3}Kg/m.s}$). 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.

Table. 1 Wigley Hull. Total resistance coefficient

	Experimental	Numerical
Test 1	5.2 ${\textstyle 10^{-3}}$	4.9 ${\textstyle 10^{-3}}$
Test 2	5.2 ${\textstyle 10^{-3}}$	5.3 ${\textstyle 10^{-3}}$
Test 3	4.9 ${\textstyle 10^{-3}}$	5.1 ${\textstyle 10^{-3}}$

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

[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.

[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.

[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.

[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–September 3, Munich, Germany, 1999.

[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.

[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–267, 1998.

[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–1407, 1997.

[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.

[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–221, 1999.

[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–2, 355–370, 2000.

[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.

[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–292, 1988.

[13] E. Oñate, C. Sacco and S. Idelsohn, “A finite point method for incompressible flow problems”, Computing and Visualization in Science, 2, 67–75, 2000.

[14] G. Chiandusi, G. Bugeda and E. Oñate, “A simple method for update of finite element meshes”, Commun, Numer. Meth. Engng., 16, 1–9, 2000.

[15] O.C. Zienkiewicz and R.C. Taylor,The finite element method, 5th Edition, 3 Volumes, Butterworth–Heinemann, 2000.

[16] R. Codina, A finite element model for incompressible flow problems, Ph.D. Thesis, Univ. Politécnica de Catalunya, Barcelona, Spain, June 1992.

[17] M. Fortin and F. Thomasset, “Mixed finite element methods for incompressible flow problems”, J. Comput. Phys., 31, 113–145, 1979.

[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–259, 1982.

[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–328, 1986.

[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–192, 1990.

[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–89, 1986.

[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–233, 1992.

[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–292, 1994.

[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–67, 1997.

[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–255, 1999.

[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–189, 1989.

[27] T.E. Tezduyar, S. Mittal, S.E. Ray and R. Shih, “Incompressible flow computations with stabilized bilinear and linear equal order interpolation velocity–pressure elements”, Comput. Methods Appl. Mech. Engrg., 95, 221–242, 1992.

[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).

[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.

[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–401, 1995.

[31] R. Codina, “A stabilized finite element method for generalized stationary incompressible flows”, Publication PI-148, CIMNE, Barcelona, February 1999.

[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–301, 2000.

[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–129, 1992.

[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.

[35] F. Brezzi, L.P. Franca, T.J.R. Hughes and A. Russo, ``$b=\int g$, Comput. Methods Appl. Mech. Engrg., 145, 329–339, 1997.

[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–336, 1986.

[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–342, 1993.

[38] R. Codina, “Comparison of some finite element methods for solving the diffusion-convection-reaction equation”, Comput. Methods Appl. Mech. Engrg., 156, 185–210, 1998.

[39] R. Codina, “On stabilized finite element methods for linear systems of convection-diffusion-reaction equation”, Comput. Methods Appl. Mech. Engrg., 188, 61–83, 2000.

[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–257, 2000.

[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–351, 1992.

[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–371, 1992.

[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–282, 1994.

[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–953, 1995.

[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–1369, 1997.

[46] N. Maman and C. Farhat, “Matching fluid and structure meshes for aeroelastic computations: A parallel approach”, Computers and Structures, 54, 779–785, 1995.

[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–321, 1997.

[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.

[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–602, 1998.

[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–1012, 1998.

[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–508, 1999.

[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–293, 1999.

[53] O.C. Zienkiewicz, J.Z. Zhu, “The Superconvergent patch recovery (SPR) and adaptive finite element refinement”, Comput. Methods Appl. Mech. Engrg., 101, 207–224, 1992.

[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–242, 1997.

[55] J.H. Duncan, “The breaking and non-breaking wave resistance of a two-dimensional hydrofoil”, J. Fluid Mech., Vol. 126, 1983.

[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.

[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.

[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.

[59] R. Codina, “Pressure stability in fractional step finite element methods for incompressible flows”. Submitted to Comput. Methods Appl. Mech. Engrg., 2000.

[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.