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 
A stabilized semiimplicit fractional step finite element method for solving coupled fluidstructure 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 semiimplicit algorithm for the analysis of fluidstructure interaction problems in totally or partially submerged bodies is presented.
keywords Fluidstructure interaction, surface waves, finite element method.
Accurate prediction of the fluidstructure 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 fluidstructure interaction problem in this case are mainly due to the following reasons:
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 advectivediffusive transport and fluid flow problems 1.
The stabilized governing equations are written in an arbitrary LagrangianEulerian (ALE) form to account for the effect of relative movement between the mesh and the fluid points. These equations are solved in spacetime using a semiimplicit 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 semiimplicit 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 fluidstructure interaction problems with free surface waves.
The finite element solution of the incompressible NavierStokes equations with the classical Galerkin method may suffer from numerical instabilities from two main sources. The first is due to the advectivediffusive 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 velocitypressure interpolations satisfying the infsup 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 advectivediffusivereactive problem 38.
Applications of stabilized GLS FEM to fluidstructure 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 fluidstructure 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 advectivediffusive 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 NavierStokes 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 semiimplicit 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 shipwave problems, including coupled fluidstructure interaction situations, are presented.
For the sake of preciseness, the basic ideas of the FIC method are given next.
Let us consider a sourceless transport problem over a one dimensional domain . The balance of flux over a domain of finite size belonging to can be written as (Figure 1)

(1) 
where and are the end points of the finite size domain of length . In eq.(1) and represent the values of the flux at points and , respectively.
Equilibrium of fluxes in a finite balance domain 
Figure 1: Equilibrium of fluxes in a finite balance domain 
For instance, in an 1D advectivediffusive problem the flux , where is the transported variable (i.e. the temperature in a thermal problem), is the advective velocity and and are the advective and diffusive material parameters, respectively.
The flux can be expressed in terms of the values at point by the following Taylor series expansion

(2) 
Substituting (2) into (1) gives after simplification and neglecting cubic terms in

(3) 
where all terms are evaluated at the arbitrary point .
Eq. (3) is the finite form of the balance equation over the domain . The underlined term in eq.(3) introduces the necessary stabilization for the discrete solution of eq.(3) using any numerical technique. Distance 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 the standard infinitesimal form of the balance equation is recovered.
The above process can be extended to derive the stabilized balance differencial equations for any problem in fluid or solid mechanics as

(4) 
where is the standard form of the th differential equation for the infinitesimal problem, are the dimensions of the domain where balance of fluxes, forces, etc. is enforced, and 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 steadystate and transient advectivediffusive 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.
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 LagrangianEulerian (ALE) form as 2
Momentum

(5) 
Mass balance

(6) 
Free surface

(7) 
where

Above, is the velocity along the th global reference axis, is the velocity of the mesh nodes and is the relative velocity between the moving mesh and the fluid point , is the (constant) density of the fluid, is the dynamic pressure defined as where is the absolute pressure and is the vertical coordinate, is the wave elevation (measured with respect to a reference flat surface) and are the viscous stresses related to the viscosity by the standard expression

(12) 
where is the Kronecker delta.
The boundary conditions for the stabilized problem are written as

(13) 

(14) 
where are the components of the unit normal vector to the boundary and and are prescribed tractions and displacements on the boundaries and , respectively.
The underlined terms in eqs.(5)–(7) introduce the necessary stabilization for the approximated numerical solution.
The characteristic length distances and represent the dimensions of the finite domain where balance of momentum and mass is enforced. On the other hand, the characteristic distances 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 and 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 NavierStokes 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.
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.
Starting from eq.(5) we can obtain after algegraic rearrangement

(15) 
where

(16) 
and is given by eq.(8).
Inserting eq.(15) into eq.(6) gives

(17) 
with

(18) 
In above

(19) 
Note that for where is a typical grid dimension (i.e. the average element size), the value of is simply

The stabilization parameter has now the form traditionally used in the GLS formulation for the viscous (Stokes) limit () and the inviscid (Euler) limit () and deduced from adhoc extensions of the 1D advectivediffusive problems 18. Note, however, that the general form of the stabilization parameter 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.
The momentum equations (5) are first discretized in time using the following scheme

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

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

(23) 
where

(24) 
Standard fractional step procedures neglect the contribution from the terms involving in eq. (23). These terms have an additional stabilization effect which improves the numerical solution when the values of are small. Also the influence of the stabilization term 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 crossderivative terms have been kept within the term 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

(25) 
A typical solution in time includes the following steps.
Step 1. Solve explicitely for the so called fractional velocities using eq. (21).
Step 2. Solve for the dynamic pressure field 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 at the updated configuration for each mesh node using eq.(22)
Step 4. Compute the new position of the free surface elevation 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 and the viscous stresses .
Step 6. Compute the new position of mesh nodes in the fluid domain at time 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 obtained from the stress equilibrium condition (neglecting surface tension effects) as

(26) 
The dynamic pressure at the free surface is computed by

(27) 
where 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.
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 and is chosen in the examples shown in the paper. Similarly, linear triangles are chosen to interpolate on the free surface mesh. The velocity and pressure fields are interpolated within each element in the standard finite element manner as

where are the linear shape functions interpolating the velocity and pressure fields, respectively, and denote nodal values 15.
Similarly the wave height is discretized as

(30) 
where 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.

(31) 
with

(32) 

(33) 

(34) 
The solution of eq.(31) can be speeded up by diagonalizing matrix . 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 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.

(35) 

(36) 

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

(38) 
where is given by eq.(35) and

(39) 
The kinematic boundary conditions on the nodal velocities (eq.(14)) are imposed when solving eq.(38).
The new free surface elevation in the fluid domain is computed as

(40) 
with

(41) 

(42) 
In the derivation of eq.(42) the assumption that at the boundary line of the free surface domain has been made.
Steps 5 and 6 follow the process described in the previous section.
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 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

(43) 
The problem remains now finding the value of the characteristic length vectors . Indeed, the components of 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

(44) 
where and and are the “streamline” and “cross wind” contributions given by

where are the vectors defining the element sides ( for tetrahedra).
An alternative method for computing vector in a more consistent manner is explained in 60.
As for the free surface equation the following value of the characteristic length vector has been taken

(47) 
The streamline parameter has been obtained by eq.(45) using the value of the velocity vector over the 3 node triangles discretizing the free surface and .
The cross wind parameter has been computed by

(48) 
Note that the crosswind 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 “shockcapturing” 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 and in eqs.(5) and (7) the value has been taken for solution of the problems presented in a next section.
Different techniques have been proposed for dealing with mesh updating in fluidstructure 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 and the Poisson coefficient . Once a discretized finite element problem has been solved using, for instance, standard linear triangles (in 2D) or linear tetraedra (in 3D), the principal stresses at the center of each element are obtained as

(49) 
where are the principal strains.
Let us assume now that a uniform strain field throughout the mesh is sought. The principal stresses are then given by

(50) 
where is the unknown Young modulus for the element.
A number of criteria can be now used to find the value of . The most effective approach found in 14 is to equate the element strain energies in both analysis. Thus

Equaling eqs.(51) and (52) gives the sought Young modulus as

(53) 
Note that the element Young modulus is proportional to the element deformation as desired. Also recall that both and 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 and . 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 . Repeat the finite element solution of the linear elastic problem with prescribed boundary displacements using the new values of 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 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 semisubmerged 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 (inplane) 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.
All the examples shown next have been solved in a standard PC Pentium II 450 Mhz with a memory of 128 Mb.
A 2D submerged NACA0012 profile at 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 where is the incoming flow velocity at infinity.
The stationary free surface and the pressure distribution in the domain are shown in Figure 3. The nondimensional wave heights compare well with the experimental results of 55.
File:Draft 177 844678375fig3a.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 
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) 
The last case considered here is the well known Wigley Hull, given by the analytical formula where and 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 844678375figura6b.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 (). 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.
Experimental  Numerical  
Test 1  5.2  4.9 
Test 2  5.2  5.3 
Test 3  4.9  5.1 
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.
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 semiimplicit fractional step finite element method provides a straightforward and stable algorithm for fluidstructure interaction analysis.
The meshmoving scheme presented ensures minimum mesh distortion for large mesh displacements. The stabilized finite element method developed is adequate for solving large scale fluidstructure 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.
Support for this work was provided by the European Community through projects BriteEuram BR 9674342 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.
[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 fluidstructure interaction problems involving surface waves, in Computational Methods for FluidStructure Interaction, T. Kvamsdal et al., (Eds.), TAPIR Publishers, Norway, 1999.
[4] E. Oñate and J. García, “A methodology for analysis of fluidstructure 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 advectivediffusive transport and fluid flow problems”, Comput. Methods Appl. Mech. Engrg., Vol. 151, 12, 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 advectivediffusive problems Int. J. Num. Meth. Fluids, Vol. 25, pp. 1385–1407, 1997.
[8] E. Oñate, J. García and S. Idelsohn, An alphaadaptive approach for stabilized finite element solution of advectivediffusive 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 spacetime finite element methods for advectivediffusive 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 advectivediffusive 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/PetrovGalerkin formulation for convection dominated flows with particular emphasis on the incompressible NavierStokes 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 advectivediffusive systems”, Comput Methods Appl. Mech. Engrg., 58, pp. 305–328, 1986.
[20] P. Hansbo and a. Szepessy, “A velocitypressure streamline diffusion finite element method for the incompressible NavierStokes 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 BabuskaBrezzi condition: A stable PetrovGalerkin 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 NavierStokes 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/leastsquares method for advectivediffusive 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, 86985, (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. 89, 886913, 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 PI148, CIMNE, Barcelona, February 1999.
[32] R. Codina and J. Blasco, “Stabilized finite element method for the transient NavierStokes 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 residualfree 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, ``, 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 advectivediffusive system, Comput. Methods Appl. Mech. Engrg., 58, 329–336, 1986.
[37] R. Codina, “A discontinuitycapturing crosswind dissipation for the finite element solution of the convectiondiffusion equation”, Comput. Methods Appl. Mech. Engrg., 110, 325–342, 1993.
[38] R. Codina, “Comparison of some finite element methods for solving the diffusionconvectionreaction equation”, Comput. Methods Appl. Mech. Engrg., 156, 185–210, 1998.
[39] R. Codina, “On stabilized finite element methods for linear systems of convectiondiffusionreaction 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 advectivediffusive 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 deformingspatialdomain/spacetime 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 deformingspatialdomain/spacetime procedure: II. Computation of freesurface 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 fluidbody 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 Ramair 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 fluidparticle 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 nonlocal (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 shipwave problem”, Int. J. Num. Meth. Engng., 45, 503–508, 1999.
[52] R. Löhner, C. Yang, E. Oñate and S. Idelsohn, “An unstructured gridbased 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 nonbreaking wave resistance of a twodimensional 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. 173194 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, 11751182, 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 fluidstructure interaction with surface waves using a finite calculus formulation”. To appear in Comput. Methods Appl. Mech. Engrg.