A general stabilized formulation for incompressible fluid flow using finite calculus and the finite element method

<!-- metadata commented in wiki content
==A GENERAL STABILIZED FORMULATION FOR INCOMPRESSIBLE FLUID FLOW USING FINITE CALCULUS AND THE FINITE ELEMENT METHOD==

'''E. Oñate, J. García, G. Bugeda '''and''' S.R. Idelsohn'''

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

==Abstract==

We present a general formulation for incompressible fluid flow analysis using the finite element method (FEM). The necessary stabilization for dealing with convective effects and the incompressibility condition are introduced via the so called finite calculus (FIC) method. The extension of the standard eulerian form of the equations to an arbitrary lagrangian-eulerian (ALE) frame adequate for treating fluid-structure interaction problems is presented. The fully lagrangian form is also discussed. Details of an effective mesh updating procedure are presented together with a method for dealing with free surface effects of importance for ship hydrodynamic analysis and many other fluid flow problems. Examples of application of the eulerian, the ALE and the fully lagrangian flow descriptions are presented.

'''Keywords:''' Stabilized formulation, incompressible fluid flow, finite calculus, finite element method, ship hydrodynamics.

==1 INTRODUCTION==

The development of efficient and robust numerical methods for analysis of incompressible flows has been a subject of intensive research in last decades. Much effort has been spent in developing the so called stabilized numerical methods overcoming the two main sources of instability in incompressible flow analysis, namely those originated by the high values of the convective terms and those induced by the difficulty in satisfying the incompressibility conditions.

The solution of above problems in the context of the finite element method (FEM) has been attempted in a number of ways <span id='citeF-1'></span>[[#cite-1|[1]]]. The underdiffusive character of the Galerkin FEM for high convection flows (which incidentaly also occurs for the central finite difference (FD) and finite volume (FV) methods <span id='citeF-2'></span>[[#cite-2|[2]]]) has been corrected by adding some kind of artificial viscosity terms to the standard Galerkin equations.

A popular way to overcome the problems with the incompressibility constraint is by introducing a pseudo-compressibility in the flow and using implicit and explicit algorithms developed for this kind of problems such as artificial compressibility schemes <span id='citeF-4'></span>[[#cite-4|[4]]-<span id='citeF-6'></span>[[#cite-6|6]]] and preconditioning techniques <span id='citeF-7'></span>[[#cite-7|[7]]]. Other FEM schemes with good stabilization properties for the convective and incompressibility terms  are based in Petrov-Galerkin (PG) techniques. The background of PG methods are the non-centred (upwind) schemes for computing the first derivatives of the convective operator in FD and FV methods <span id='citeF-2'></span>[[#cite-2|[2]],<span id='citeF-8'></span>[[#cite-8|8]]]. More recently a general class of Galerkin FEM has been developed where the standard Galerkin variational form is extended with adequate residual-based terms in order to achieve a stabilized numerical scheme. Among the many methods of this kind in the finite element universe we can name the Streamline Upwind Petrov Galerkin (SUPG) method <span id='citeF-1'></span>[[#cite-1|[1]],<span id='citeF-29'></span>[[#cite-9|9]]-<span id='citeF-18'></span>[[#cite-18|18]]] the Galerkin Least Square (GLS) method <span id='citeF-19'></span>[[#cite-19|[19]],<span id='citeF-20'></span>[[#cite-20|20]]], the Taylor-Galerkin method <span id='citeF-21'></span>[[#cite-21|[21]]], the Characteristic Galerkin method <span id='citeF-22'></span>[[#cite-22|[22]]-<span id='citeF-24'></span>[[#cite-24|24]]] and its variant the Characteristic Based Split (CBS) method <span id='citeF-25'></span>[[#cite-25|[25]],<span id='citeF-26'></span>[[#cite-26|26]]], pressure gradient operator methods <span id='citeF-27'></span>[[#cite-27|[27]]] and the Subgrid Scale (SS) method <span id='citeF-28'></span>[[#cite-28|[28]]-<span id='citeF-30'></span>[[#cite-30|30]]]. A good review of these methods can be found in <span id='citeF-31'></span>[[#cite-31|[31]]].

In this paper a stabilized finite element formulation for incompressible flows is derived in a different manner. The starting point are the modified governing differential equations of the fluid flow problem formulated via a finite calculus (FIC) approach <span id='citeF-32'></span>[[#cite-32|[32]]]. The FIC method is based in invoking the balance of fluxes in a fluid domain of finite size. This introduces naturally additional terms in the classical differential equations of infinitesimal fluid mechanics which are a function of the balance domain dimensions. The new terms in the modified governing equations provide the necessary stabilization to the discrete equations obtained via the standard Galerkin finite element method <span id='citeF-33'></span>[[#cite-33|[33]]-<span id='citeF-39'></span>[[#cite-39|39]]].

The layout of the chapter is the following. In the next section, the main concepts of the FIC approach are introduced via a simple 1D convection-diffusion model problem. Then the basic FIC equations for incompressible flow problems are presented. The finite element discretization is introduced and the resulting matrix formulation is detailed. Both monolithic and fractional step schemes for the transient solution are presented.

The basic formulation is extended to account for free surface wave effects by using an arbitrary eulerian-lagrangian (ALE) frame and introducing the free surface boundary conditions. Here the numerical treatment of the free surface equation using the FIC method is presented. The analysis of fluid-structure interaction problems involving the movement of floating or submerged solids in a fluid is also discussed. These problems require the displacement of the mesh nodes in accordance with the motion of the structure or the free surface and here a simple and effective algorithm for updating the mesh nodes is described. In the last part of the chapter the fully lagrangian formulation for fluid flow analysis is presented as a particular case of the ALE form. The lagrangian description has many advantages for tracking the displacement of fluid particles in flows where large motions of the fluid surface occur such in the case of breaking waves, splashing of water, filling of moulds, etc.  A positive feature of the lagrangian formulation is that the convective terms dissapear in the governing equations of the fluid flow; in return the updating of the mesh at almost every time step is now a necessity and efficient algorithms for mesh generation must be used.

The examples show the efficiency of the eulerian, ALE and fully lagrangian formulations to solve classical fluid flow problems, as well as fluid-structure interaction situations involving contact with moving solids, waves around ships and large motions of the free surface, among others.

==2 THE FINITE CALCULUS METHOD==

We will consider a convection-diffusion problem in a 1D domain <math display="inline">\Omega </math> of length <math display="inline">L</math>. The equation of balance of fluxes in a subdomain of size <math display="inline">d</math> belonging to <math display="inline">\Omega </math> (Figure [[#img-1|1]]) is written as

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

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

where <math display="inline">q_A</math> and <math display="inline">q_B</math> are the incoming and outgoing fluxes at points <math display="inline">A</math> and <math display="inline">B</math>, respectively. The flux <math display="inline">q</math> includes both convective and diffusive terms; i.e. <math display="inline">q=u\phi - k{d\phi \over dx}</math>, where <math display="inline">\phi </math> is the transported variable, <math display="inline">v</math> is the velocity and <math display="inline">k</math> is the diffusitivity of the material.

Let us express now the fluxes <math display="inline">q_A</math> and <math display="inline">q_B</math> in terms of the flux at an arbitrary point <math display="inline">C</math> within the balance domain (Figure [[#img-1|1]]). Expanding <math display="inline">q_A</math> and <math display="inline">q_B</math> in Taylor series around point <math display="inline">C</math> up to second order terms gives

<div id='eq-2'></div>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>q_A= q_C - d_1 \frac{dq}{d x}\vert _C+ \frac{d^2_1}{2}\frac{d^2q}{dx^2}\vert _C + O(d^3_1)\quad ,\quad  q_B= q_C + d_2 \frac{dq}{d x}\vert _C+\frac{d^2_2}{2}\frac{d^2q}{dx^2}\vert _C + O(d^3_2) </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (2)
|}

Substituting eq.([[#eq-2|2]]) into eq.([[#eq-1|1]]) gives after simplification

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

where <math display="inline">h=d_1-d_2</math> and all derivatives are computed at point <math display="inline">C</math>.

Standard calculus theory assumes  that the domain <math display="inline">d</math> is of infinitesimal size and the resulting balance equation is simply <math display="inline">{dq\over dx}=0</math>. We will relax this assumption and allow the balance domain to have a ''finite size''. The new balance equation ([[#eq-3|3]]) incorporates now the underlined term which introduces the ''characteristic length'' <math display="inline">h</math>. Obviously, accounting for higher order terms in eq.([[#eq-2|2]]) would lead to new terms in eq.([[#eq-3|3]]) involving higher powers of <math display="inline">h</math>.

Distance <math display="inline">h</math> in eq.([[#eq-3|3]]) can be interpreted as a free parameter depending, of course, on the location of point <math display="inline">C</math> (note that <math display="inline">-d\le h \le d</math>). However, the fact that eq.([[#eq-3|3]]) is the exact balance equation (up to second order terms) for any 1D domain of finite size and that the position of point <math display="inline">C</math> is arbitrary, can be used to derive numerical schemes with enhanced properties simply by computing the characteristic length parameter from an adequate “optimality” rule.

Consider, for instance, the modified equation ([[#eq-3|3]]) applied to the convection-diffusion problem. Neglecting third order derivatives of <math display="inline">\phi </math>, eq.([[#eq-3|3]]) can be written in an explicit form as

<div id='eq-4'></div>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>-u \frac{d\phi }{dx}+\left(k+\frac{u h}{2}\right)\frac{d^2\phi }{dx^2}=0 </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (4)
|}

We see  that the modified equation via the FIC method introduces ''naturally'' an additional diffusion term into the standard convection-diffusion equation. This is the basis of the popular “artificial diffusion” procedure <span id='citeF-1'></span>[[#cite-1|[1]],<span id='citeF-2'></span>[[#cite-2|2]],<span id='citeF-8'></span>[[#cite-8|8]],<span id='citeF-36'></span>[[#cite-36|36]]]. The characteristic length <math display="inline">h</math> is typically expressed as a function of the cell or element dimensions. The optimal or critical value of <math display="inline">h</math> for each cell or element can be computed from numerical stability conditions such as obtaining a physically meaningful solution, or even obtaining “exact” nodal values <span id='citeF-32'></span>[[#cite-32|[32]]-<span id='citeF-39'></span>[[#cite-39|39]]].

Equation ([[#eq-3|3]]) can be extended to account for source effects. The full stabilized equation can be then written in compact form as

<div id='eq-5'></div>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>r- \underline{\frac{h}{2} \frac{dr}{dx}}=0 </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (5)
|}

with

<div id='eq-6'></div>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>r = -u \frac{d\phi }{dx}+ \frac{d}{dx}\left(k \frac{d\phi }{dx}\right)+ Q </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (6)
|}

where <math display="inline">Q</math> is the external source. For consistency a “finite” form of the Neumann boundary condition should be used. This can be readily obtained by invoking balance of fluxes in a domain of finite size next to the boundary <math display="inline">\Gamma _q</math> where the external (diffusive) flux is prescribed to a value <math display="inline">q_p</math>. The modified Neumann boundary condition can be written as <span id='citeF-32'></span>[[#cite-32|[32]]]

<div id='eq-7'></div>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>k \frac{d\phi }{dx}+ q_p- \underline{\frac{h}{2}r}=0 \quad \hbox{at } \, \Gamma _q </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (7)
|}

The definition of the problem is completed with the standard Dirichlet condition prescribing the value of <math display="inline">\phi </math> at the boundary <math display="inline">\Gamma _\phi </math>.

The underlined terms in Eqs.([[#eq-5|5]]) and ([[#eq-7|7]]) introduce the necessary stabilization in the discrete solution of the problem using whatever numerical scheme. For details see <span id='citeF-32'></span>[[#cite-32|[32]]-<span id='citeF-41'></span>[[#cite-41|41]]].

The time dimension can be  introduced in the FIC method by considering the balance equation in a space-time slab domain <span id='citeF-32'></span>[[#cite-32|[32]],<span id='citeF-35'></span>[[#cite-35|35]],<span id='citeF-36'></span>[[#cite-36|36]]]. Quite generally the FIC equation can be written for any problem in mechanics as <span id='citeF-32'></span>[[#cite-32|[32]]].

<div id='eq-8'></div>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>r_i - \underline{\frac{h_j}{2}{\partial r_i \over \partial x_j}}-\underline{\frac{\delta }{2}{\partial r_i \over \partial t}}=0 \quad , \begin{array}{l}i=1,n_b\\ j=1,n_d \end{array} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (8)
|}

where <math display="inline">r_i</math> is the ith standard differential equation of the infinitesimal theory, <math display="inline">h_j</math> are characteristic length parameters, <math display="inline">\delta </math> is a time stabilization parameter and <math display="inline">t</math> the time; <math display="inline">n_b</math> and <math display="inline">n_d</math> are respectively the number of balance equations and the number of dimension of the problem along which balance of fluxes or forces is invoked (i.e., <math display="inline">n_d =2</math> for 2D problems, etc.).

For example, in the case of the convection-diffusion problem <math display="inline">n_b=1</math>, Eq.([[#eq-8|8]]) is particularized as

<div id='eq-9'></div>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>r- \underline{\frac{h_j}{2}{\partial r \over \partial x_j}}-\underline{\frac{\delta }{2}{\partial r \over \partial t}}=0 \qquad ,\quad j=1,n_d </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (9)
|}

with <math display="inline">r:=-\left(\displaystyle {\partial \phi  \over \partial t}+u_i {\partial \phi  \over \partial x_i}\right)+ \displaystyle \frac{d}{dx_k}\left(k \frac{d\phi }{dx_k}\right)+ Q</math>.  

The modified Neumann boundary conditions in the FIC formulation can be expressed in the general case as

<div id='eq-10'></div>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>q_{ij}n_j - t_i - \underline{h_jn_j r_i}=0\quad \hbox{on } \Gamma _q \quad i=1,n_d </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (10)
|}

where <math display="inline">q_{ij}</math> are the generalized “fluxes” (such as the heat fluxes in a heat transfer problem or the stresses in solid or fluid mechanics), <math display="inline">t_i</math> are the prescribed values of the boundary fluxes and <math display="inline">n_j</math> are the components of the outward normal to the Neumann boundary <math display="inline">\Gamma _q</math>. For the transient case the initial boundary condition  should also be specified <span id='citeF-32'></span>[[#cite-32|[32]],<span id='citeF-35'></span>[[#cite-35|35]],<span id='citeF-36'></span>[[#cite-36|36]]].

In Eqs.([[#eq-8|8]])-([[#eq-10|10]]) we have underlined once more the terms introduced by the FIC approach which are essential for deriving stabilized numerical formulations.

The starting point in the next section are the FIC equation for a viscous incompressible fluid. A simplified version of the equation will be chosen neglecting the time stabilization term as this is not relevant for the purposes of this work.

==3 GENERAL FIC EQUATIONS FOR VISCOUS INCOMPRESSIBLE FLOW==

The FIC governing equations for a viscous incompressible fluid can be written as

'''Momentum'''
<span id="eq-11"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>r_{m_i} - \underline{{1\over 2} h_j{\partial r_{m_i} \over \partial x_j}}=0 </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (11)
|}

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

where
<span id="eq-13"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>r_{m_i} = \rho \left({\partial u_i \over \partial t}+{\partial (u_iu_j) \over \partial x_j}\right)+ {\partial p \over \partial x_i}- {\partial s_{ij} \over \partial x_j}-b_i</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (13)
|}
<span id="eq-14"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math> r_d = {\partial u_i \over \partial x_i}\qquad i,j = 1, n_d </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (14)
|}

Above <math display="inline">u_i</math> is the velocity along the ith global axis, <math display="inline">\rho </math> is the (constant) density of the fluid, <math display="inline">p</math> is the absolute pressure (defined positive in compression), <math display="inline">b_i</math> are the body forces and <math display="inline">s_{ij}</math> are the viscous deviatoric stresses related to the viscosity <math display="inline">\mu </math> by the standard expression
<span id="eq-15"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>s_{ij}=2\mu \left(\dot \varepsilon _{ij} - \delta _{ij} {1\over 3} {\partial u_k \over \partial x_k}\right) </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (15)
|}

where <math display="inline">\delta _{ij}</math> is the Kronecker delta and the strain rates <math display="inline">\dot \varepsilon _{ij}</math> are
<span id="eq-16"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\dot \varepsilon _{ij}={1\over 2} \left({\partial u_i \over \partial x_j}+{\partial u_j \over \partial x_i}\right) </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (16)
|}

The FIC boundary conditions  are
<span id="eq-17"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>n_j \sigma _{ij} -t_i + \underline{{1\over 2} h_j n_j r_{m_i}}=0 \quad \hbox{on }\Gamma _t </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (17)
|}
<span id="eq-18"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>u_j - u_j^p =0 \quad \hbox{on }\Gamma _u </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (18)
|}

and the initial condition is <math display="inline">u_j =u_j^0</math> for <math display="inline">t=t_0</math>.

In Eqs.([[#eq-17|17]]) and ([[#eq-18|18]]) <math display="inline">t_i</math> and <math display="inline">u_j^p</math> are surface tractions and prescribed displacements on the boundaries <math display="inline">\Gamma _t</math> and <math display="inline">\Gamma _u</math>, respectively, <math display="inline">n_j</math> are the components of the unit normal vector to the boundary and <math display="inline">\sigma _{ij}</math> are the total stresses given by <math display="inline">\sigma _{ij}=s_{ij}-\delta _{ij}p</math>.  The sign in front the stabilization term in Eq.([[#eq-17|17]]) is positive due to the definition of <math display="inline">r_{m_i}</math> in Eq.([[#eq-13|13]]).

The <math display="inline">h_i's</math> in above equations are characteristic lengths of the domain where balance of momentum and mass is enforced. In Eq.([[#eq-17|17]]) these lengths define the domain where equilibrium of boundary tractions is established <span id='citeF-32'></span>[[#cite-32|[32]]].

Eqs.([[#eq-11|11]])-([[#eq-18|18]]) are the starting  point for deriving stabilized finite element methods for solving the incompressible Navier-Stokes equations using equal order interpolation for the velocity and pressure variables <span id='citeF-37'></span>[[#cite-37|[37]]-<span id='citeF-39'></span>[[#cite-39|39]]]. Application of the FIC formulations to meshless analysis of fluid flow problems using the finite point method can be found in <span id='citeF-40'></span>[[#cite-40|[40]],<span id='citeF-41'></span>[[#cite-41|41]]].

We note that the “conservative” form of the convective terms in Eq.([[#eq-13|13]]) and the presence of the volumetric strain rate in the constitutive equation ([[#eq-15|15]]) do not take advantage of the incompressibility condition. These forms are useful for obtaining the relationship between the derivative of the volumetric strain rate and the  momentum equations as shown in the next section. However, the standard  form of the governing equations for incompressible flows will be used for the final derivation of the discretized FEM equations.

===3.1 Stabilized integral forms===

From Eq.([[#eq-11|11]]) it can be obtained (taking into account Eq.([[#eq-15|15]]))

<span id="eq-19"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>{\partial r_d \over \partial x_i}=-{1\over a_i} \left[\bar r_{m_i}-{h_j\over 2} {\partial r_{m_i} \over \partial x_j}\right]-  {\rho u_i h_k\over 2a_i} {\partial r_d \over \partial x_k}\quad i,j=1,n_d\quad ,\quad  k\not =i </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (19)
|}

where
<span id="eq-20a"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>a_i = {2\mu \over 3} +{\rho u_i h_i\over 2}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (20a)
|}

and
<span id="eq-20b"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\bar r_{m_i}=\rho \left({\partial u_i \over \partial t}+u_j {\partial u_i \over \partial x_j}\right)+ {\partial p \over \partial x_i}-{\partial  \over \partial x_j} (2\mu \dot \varepsilon _{ij}) - b_i</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (20b)
|}

Substituting Eq.([[#eq-19|19]]) into Eq.([[#eq-12|12]]) and retaining the terms involving the derivatives of <math display="inline">r_{m_i}</math> with respect to <math display="inline">x_i</math> only, leads to the following expression for the stabilized mass balance equation
<span id="eq-21"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>r_d - \sum \limits _{i=1}^{n_d} \tau _i {\partial r_{m_i} \over \partial x_i}=0 </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (21)
|}

with
<span id="eq-22"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\tau _i = \left({8\mu \over 3h_i^2}+{2\rho u_i\over h_i}\right)^{-1} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (22)
|}

The <math display="inline">\tau _i</math>'s in Eq.([[#eq-21|21]]) are termed in the stabilization literature ''intrinsic time parameters''. It is interesting to note that  these parameters take here the  values of <math display="inline">\tau _i =\displaystyle{3h_i^2\over 8\mu }</math> and <math display="inline">\tau _i =\displaystyle{h_i\over 2\rho u_i}</math> for the viscous  limit (Stokes flow) and the inviscid limit (Euler flow), respectively. Similar values for <math display="inline">\tau _i</math> (usually <math display="inline">\tau _i =\tau </math> is taken) are used in other works from ad-hoc extensions of the 1D advective-diffusive problem <span id='citeF-12'></span>[[#cite-12|[12]]-<span id='citeF-31'></span>[[#cite-31|31]]]. It is remarkable that the intrinsic time parameters have been deduced here from the general FIC formulation and this shows the possibilities of the method.

The weighted residual form of the momentum and mass balance equations (Eqs.([[#eq-11|11]]) and ([[#eq-21|21]])) is written as

'''Momentum'''
<span id="eq-23"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\int _\Omega \delta u_i \left[r_{m_i} - {h_j\over 2} {\partial r_{m_i} \over \partial x_j}\right]+ \int _{\Gamma _t} \delta u_i (\sigma _{ij} n_j - t_i + {h_j \over 2} n_j r_{m_i}) d\Gamma =0 </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (23)
|}

'''Mass balance'''
<span id="eq-24"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\int _\Omega q \left[r_d - \sum \limits _{i=1}^{n_d} \tau _i {\partial r_{m_i} \over \partial x_i}\right]d\Omega =0 </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (24)
|}

where <math display="inline">\delta u_i</math> and <math display="inline">q</math> are arbitrary weighting functions representing virtual velocity and virtual pressure fields. Integration by parts of the <math display="inline">r_{m_i}</math> terms  leads to
<span id="eq-25a"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\int _\Omega \delta u_i r_{m_i} + \int _{\Gamma _t} \delta u_i (\sigma _{ij} n_j - t_i)d\Gamma + \sum \limits _e \int _{\Omega ^e} {h_j\over 2}{\partial \delta u_i \over \partial x_j} r_{m_i} d\Omega =0 </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (25a)
|}
<span id="eq-25b"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\int _\Omega q r_d d\Omega + \int _\Omega \left[\sum \limits _{i=1}^{n_d}\tau _i {\partial q \over \partial x_i}r_{m_i} \right]d\Omega - \int _\Gamma \left[\sum \limits _{i=1}^{n_d} q \tau _i n_i r_{m_i}\right]d\Gamma =0</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (25b)
|}

The third integral in Eq.([[#eq-25a|25a]]) is expressed as a sum of the element contributions to allow for discontinuities in the derivatives of <math display="inline">r_{m_i}</math> along the element interfaces.

Also in Eq.([[#eq-25b|25b]]) we will neglect hereonwards the third integral by assuming that <math display="inline">r_{m_i}</math> is negligible on the boundaries. The deviatoric stresses and the pressure terms in the first integral of Eq.([[#eq-25a|25a]]) are integrated by parts in the usual manner. The resulting equations are

'''Momentum'''
<span id="eq-26"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\int _\Omega \left[\delta u_i\rho \left({\partial u_i \over \partial t}+u_j {\partial {u_i} \over \partial x_j}\right)+ \delta \dot \varepsilon _{ij}(s_{ij}- \delta _{ij}p )\right]d\Omega -  \int _{\Omega } \delta u_i b_i d\Omega - \int _{\Gamma _t} \delta u_i t_id\Gamma +</math>
|-
| style="text-align: center;" | <math> + \sum \limits _e \int _{\Omega ^e} {h_j\over 2}{\partial \delta u_i \over \partial x_j} r_{m_i} d\Omega =0 </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (26)
|}

'''Mass balance'''
<span id="eq-27"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\int _\Omega q {\partial u_i \over \partial x_i} d\Omega + \int _\Omega \left[\sum \limits _{i=1}^{n_d} \tau _i {\partial q \over \partial x_i} r_{m_i}\right]d\Omega =0 </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (27)
|}

We note that in Eq.([[#eq-26|26]]) we use now the standard form of the convective operator for incompressible flows (i.e. neglecting the contribution from the volumetric strain rate <math display="inline">\displaystyle {\partial u_i \over \partial x_i}</math>). Also in Eq.([[#eq-26|26]])

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\delta \dot \varepsilon _{ij} ={1\over 2} \left({\partial \delta u_i \over \partial x_j}+{\partial \delta u_j \over \partial x_i}\right)</math>
|}
|}

===3.2 Convective and pressure gradient projections===

The computation of the residual terms can be simplified if we introduce now the convective and pressure gradient projections <math display="inline">c_i</math> and <math display="inline">\pi _i</math>, respectively defined as
<span id="eq-28"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\begin{array}{l}\displaystyle c_i = r_{m_i} -\rho u_j {\partial u_i \over \partial x_j}\\ \displaystyle \pi _i = r_{m_i} - {\partial p \over \partial x_i} \end{array} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (28)
|}

We can express <math display="inline">r_{m_i}</math> in  Eqs.([[#eq.26|26]]) and ([[#eq-27|27]]) in terms of <math display="inline">c_i</math> and <math display="inline">\pi _i</math>, respectively which then become additional variables. The system of integral equations is now augmented in the necessary number of additional equations by imposing that the residual <math display="inline">r_{m_i}</math> vanishes (in average sense) for both forms given by Eqs.([[#eq-28|28]]). This gives the final system of governing equation as:
<span id="eq-29"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\int _\Omega \left[\delta u_i\rho \left({\partial u_i \over \partial t}+u_j {\partial {u_i} \over \partial x_j}\right)+ \delta \dot \varepsilon _{ij}(s_{ij}- \delta _{ij}p )\right]d\Omega -  \int _{\Omega } \delta u_i b_i d\Omega - \int _{\Gamma _t} \delta u_i t_id\Gamma +</math>
|-
| style="text-align: center;" | <math> + \sum \limits _e \int _{\Omega ^e} {h_k\over 2}{\partial (\delta u_i) \over \partial x_k} \left(\rho u_j {\partial {u_i} \over \partial x_j} + c_i\right)d\Omega =0 </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (29)
|}
<span id="eq-30"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\int _\Omega q {\partial u_i \over \partial x_i} d\Omega + \int _\Omega \sum \limits _{i=1}^{n_d} \tau _i {\partial q \over \partial x_i} \left({\partial p \over \partial x_i}+\pi _i\right)d\Omega =0 </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (30)
|}
<span id="eq-31"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\int _\Omega \delta c_i \rho \left(\rho u_j {\partial {u_i} \over \partial x_j} + c_i\right)d\Omega =0 \qquad \hbox{no sum in }i </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (31)
|}
<span id="eq-32"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\int _\Omega \delta \pi _i \tau _i \left({\partial p \over \partial x_i}+\pi _i\right)d\Omega =0\qquad \hbox{no sum in }i </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (32)
|}

with <math display="inline">i,j,k=1,n_d</math>.  In Eqs.([[#eq-31|31]]) and ([[#eq-32|32]]) <math display="inline">\delta c_i</math> and <math display="inline">\delta \pi _i</math> are appropriate weighting functions and the <math display="inline">\rho </math> and <math display="inline">\tau _i</math> weights are introduced for convenience.

==4 FINITE ELEMENT DISCRETIZATION==

We choose <math display="inline">C^\circ </math> continuous linear interpolations of the velocities, the pressure, the convection projections <math display="inline">c_i</math> and the pressure gradient projections <math display="inline">\pi _i</math> over three node triangles (2D) and four node tetrahedra (3D). The linear interpolations are written as
<span id="eq-33"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\begin{array}{l}\displaystyle u_i = \sum \limits _{j=1}^n N_j \bar u_i^j \quad , \quad p = \sum \limits _{j=1}^n N_j \bar p^j\\ \displaystyle c_i = \sum \limits _{j=1}^n N_j \bar c_i^j \quad , \quad \pi _i = \sum \limits _{j=1}^n N_j \bar \pi _i^j \end{array} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (33)
|}

where <math display="inline">n=3</math> (4) for triangles (tetrahedra), <math display="inline">\bar {(\cdot )}^j</math> denotes nodal variables and <math display="inline">N_j</math> are the linear shape functions <span id='citeF-1'></span>[[#cite-1|1]].

Substituting the approximations ([[#eq-33|33]]) into Eqs.([[#eq-29|29]]-[[#eq-32|32]]) and choosing the Galerking form with <math display="inline">\delta u_i =q=\delta c_i=\delta \pi _i =N_i</math> leads to following system of discretized equations
<span id="eq-34"></span>
<span id="eq-34a"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\displaystyle {\boldsymbol M}\dot{\bar {\boldsymbol u}} + ({\boldsymbol A}+{\boldsymbol K}+\hat{\boldsymbol K}) \bar {\boldsymbol u} - {\boldsymbol G}\bar {\boldsymbol p}+{\boldsymbol C}\bar {\boldsymbol c}={\boldsymbol f}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (34a)
|}
<span id="eq-34b"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\displaystyle {\boldsymbol G}^T \bar {\boldsymbol u} + {\boldsymbol L}\bar {\boldsymbol p}+{\boldsymbol Q}\bar {\boldsymbol \pi }={\boldsymbol 0}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (34b)
|}
<span id="eq-34c"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\displaystyle \hat {\boldsymbol C}\bar{\boldsymbol u}+ {\boldsymbol M}\bar {\boldsymbol c}={\boldsymbol 0}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (34c)
|}
<span id="eq-34d"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\displaystyle {\boldsymbol Q}^T \bar {\boldsymbol p} + \hat {\boldsymbol M}\bar {\boldsymbol \pi }={\boldsymbol 0}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (34d)
|}

where the element contributions are given by (for 2D problems)
<span id="eq-35"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\begin{array}{l}\displaystyle M_{ij}= \int _{\Omega ^e} \rho N_iN_j d\Omega \quad ,\quad A_{ij}= \int _{\Omega ^e} N_i \rho {\boldsymbol u}^T {\boldsymbol \nabla } N_j d\Omega \\ \\ \displaystyle {\boldsymbol K}_{ij}= \int _{\Omega ^e}{\boldsymbol B}_i^T {\boldsymbol D}{\boldsymbol B}_j d\Omega \quad ,\quad \hat {\boldsymbol K}_{ij}= \int _{\Omega ^e} \rho {h_ku_l\over 2}\displaystyle {\partial N_i \over \partial x_k} \displaystyle {\partial N_j \over \partial x_l}d\Omega \\ \\ \displaystyle {\boldsymbol G}_{ij}= \int _{\Omega ^e}({\boldsymbol \nabla } N_i)N_j d\Omega \quad ,\quad {\boldsymbol \nabla }= \left[{\partial  \over \partial x_1},{\partial  \over \partial x_2}\right]^T \quad ,\quad  {C}_{ij} = \int _{\Omega ^e}{h_k \over 2} \displaystyle {\partial N_i \over \partial x_k}N_j d\Omega \, d\Omega \\ \\ \displaystyle L_{ij}=\int _{\Omega ^e} {\boldsymbol \nabla }^T N_i [\tau ] {\boldsymbol \nabla } N_j d\Omega \quad ,\quad [\tau ]= \left[\begin{matrix}\tau _1 &0\\ 0 & \tau _2\\\end{matrix}\right]\quad ,\quad  \hat {C}_{ij}=\int _{\Omega ^e} N_i \rho ^2 u_k \displaystyle {\partial N_j \over \partial x_k}d\Omega \\ \\ \displaystyle {\boldsymbol Q}= [{\boldsymbol Q}^1,{\boldsymbol Q}^2] \quad ,\quad  \displaystyle Q_{ij}^k = \int _{\Omega ^e}\tau _k {\partial N_i \over \partial x_k} N_jd\Omega \\ \\ \displaystyle \hat {\boldsymbol M}= \left[\begin{matrix}\hat {\boldsymbol M}^1 &{\boldsymbol 0}\\  {\boldsymbol 0} & \hat {\boldsymbol M}^2\\\end{matrix}\right]\quad ,\quad  \hat {\boldsymbol M}_{ij} = \int _{\Omega ^e} \tau _k N_i N_j d\Omega \\ \\ \displaystyle {\boldsymbol f}_i = \int _{\Omega ^e} N_i {\boldsymbol b}d\Omega + \int _{\Gamma ^e}N_i {\boldsymbol t} d\Gamma \quad ,\quad  {\boldsymbol b}=[b_1,b_2]^T \quad ,\quad {\boldsymbol t}=[t_1,t_2]^T \end{array} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (35)
|}

with <math display="inline">i,j=1,n</math> and <math display="inline">k,l=1,n_d</math>.  

In above <math display="inline">{\boldsymbol B}_i</math> is the standard strain rate matrix and <math display="inline">{\boldsymbol D}</math> the deviatoric constitutive matrix (assuming <math display="inline">\displaystyle {\partial u_i \over \partial x_i}=0</math>). For 2D problems
<span id="eq-36"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>{\boldsymbol B}_i =\left[\begin{matrix}\displaystyle {\partial N_i \over \partial x_1}&0\\ 0 & \displaystyle {\partial N_i \over \partial x_2}\\ \displaystyle {\partial N_i \over \partial x_2} & \displaystyle {\partial N_i \over \partial x_1}\\\end{matrix}\right]\quad ,\quad {\boldsymbol D} =\mu \left[\begin{matrix}2&0&0\\ 0&2&0\\ 0&0&1\\\end{matrix}\right] </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (36)
|}

Note that the stabilization matrix <math display="inline">\hat{\boldsymbol K}</math> adds an additional orthotropic diffusivity of value <math display="inline">\rho \displaystyle{h_ku_l\over 2}</math>.

===4.1 Quasi-implicit scheme===

It can be seen that matrices <math display="inline">{\boldsymbol A}, \hat {\boldsymbol K}</math> and <math display="inline">\hat {\boldsymbol C}</math> are dependent on the velocity field. The solution process can be advanced in time in a (quasi-nearly) implicit iterative manner using the following scheme.

<span id="step-1"></span>
'''Step 1'''
<span id="eq-37"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\bar{\boldsymbol u}^{n+1,i} = \bar{\boldsymbol u}^n - \Delta t {\boldsymbol M}^{-1} [({\boldsymbol A}^{n+\theta _1,i-1}+  {\boldsymbol K} + \hat {\boldsymbol K}^{n+\theta _1,i-1})\bar{\boldsymbol u}^{n+\theta _1,i-1}-</math>
|-
| style="text-align: center;" | <math> - {\boldsymbol G}\bar{\boldsymbol p}^{n+\theta _2,i-1} + {\boldsymbol C} \bar {\boldsymbol c}^{n+\theta _3,i-1}-{\boldsymbol f}^{n+1}] </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (37)
|}
<span id="step-2"></span>
'''Step 2'''
<span id="eq-38"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\bar{\boldsymbol p}^{n+1,i}=-{\boldsymbol L}^{-1} [{\boldsymbol G}^T \bar {\boldsymbol u}^{n+1,i}+{\boldsymbol Q}\bar {\boldsymbol \pi }^{n+\theta _4,i-1}] </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (38)
|}
<span id="step-3"></span>
'''Step 3'''
<span id="eq-39"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\bar {\boldsymbol c}^{n+1,i}= - {\boldsymbol M}^{-1} \hat {\boldsymbol C}^{n+1,i}\bar {\boldsymbol u}^{n+1,i} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (39)
|}
<span id="step-4"></span>
'''Step 4'''
<span id="eq-40"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\bar {\boldsymbol \pi }^{n+1,i}= - \hat {\boldsymbol M}^{-1}{\boldsymbol Q}^T \bar {\boldsymbol p}^{n+1,i} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (40)
|}

where <math display="inline">0\le \theta _i\le 1</math>.

In above <math display="inline">\bar{(\cdot )}^{n,i}</math> denote nodal values at the <math display="inline">n</math>th time step and the ith iteration. Note that <math display="inline">{\boldsymbol A}^{n+\theta _1,i-1}\equiv {\boldsymbol A}(\bar {\boldsymbol u}^{n+\theta _1,i-1})</math> etc. Also <math display="inline">(\cdot )^{n+\theta _i,0}\equiv (\cdot )^n</math> for the computations in step 1 at the onset of the iterations.

Steps [[#step-1|1]], [[#step-3|3]] and [[#step-4|4]] can be solved explicitely by choosing a ''lumped (diagonal) form'' of matrices <math display="inline">{\boldsymbol M}</math> and <math display="inline">\hat {\boldsymbol M}</math>. In this manner the main computational cost is the solution of step [[#step-2|2]] involving the inverse of a Laplacian matrix. This can be solved very effectively using an iterative method such as the conjugate gradient method or similar.

For <math display="inline">\theta _i\not =0</math> the iterative proces is unavoidable. The iterations follow until convergence is reached in an adequate error norm in terms of the velocity and pressure variables, or the residuals <math display="inline">r_{m_i}</math> and <math display="inline">r_d</math>. Indeed some ot the <math display="inline">\theta _i</math>'s in Eqs.([[#eq-37|37]])-([[#eq-40|40]]) can be made equal to zero. Note that for <math display="inline">\theta _2 =0</math> the algorithm is unconditionally unstable. A particularity interesting and simple semi-implicit form is obtained by making <math display="inline">\theta _1 = \theta _3=\theta _4=0</math>. Now all steps can be solved explicitely with exception of Step [[#step-2|2]] for the pressure, which still requires the solution of a simultaneous system of equations.

Convergence of this solution scheme is however difficult for some problems. An enhanced version of the algorithm can be obtained by simply adding the term <math display="inline">\hat {\boldsymbol  L} (\bar {\boldsymbol p}^{n+1,i} - \bar {\boldsymbol p}^{n+1,i-1})</math> where <math display="inline">\hat L_{ij} =\Delta t \int _{\Omega ^e} {\boldsymbol \nabla }^T N_i {\boldsymbol \nabla } N_j d\Omega </math> to the equation for the computation of the pressure in the second step. The new term acts as a preconditioner of the pressure equation given now by
<span id="eq-41"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\bar {\boldsymbol p}^{n+1,i}=-[{\boldsymbol L} + \hat {\boldsymbol L}]^{-1}[{\boldsymbol G}^T \bar {\boldsymbol u}^{n+1,i}+\hat {\boldsymbol L}\bar {\boldsymbol p}^{n+1,i-1} +{\boldsymbol Q}\bar {\boldsymbol \pi }^{n+\theta _4,i-1}] </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (41)
|}

Note that the added term vanishes for the converged solution (i.e. when <math display="inline">\bar {\boldsymbol p}^{n+1,i}= \bar {\boldsymbol p}^{n+1,i-1}</math>).

An alternative to above algorithm is to use the fractional step method described in the nex section.

===4.2 Fractional step method===

An alternative algorithm can be obtained by splitting the pressure from the momentum equations as follow
<span id="eq-42a"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\bar {\boldsymbol u}^* = \bar {\boldsymbol u}^n -\Delta t {\boldsymbol M}^{-1}[({\boldsymbol A}^{n+\theta _1}+  {\boldsymbol K} + \hat {\boldsymbol K}^{n+\theta _1})\bar{\boldsymbol u}^{n+\theta _1}-\alpha {\boldsymbol G}\bar{\boldsymbol p}^{n} + {\boldsymbol C} \bar {\boldsymbol c}^{n+\theta _3}-{\boldsymbol f}^{n+1}]</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (42a)
|}
<span id="eq-42b"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\bar{\boldsymbol u}^{n+1}= \bar{\boldsymbol u}^* + \Delta t {\boldsymbol M}^{-1}{\boldsymbol G}\delta \bar{\boldsymbol p}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (42b)
|}

In Eq.([[#eq-42a|42a]]) <math display="inline">\alpha </math> is a variable taking values equal to zero or one. For <math display="inline">\alpha =0</math>, <math display="inline">\delta p \equiv p^{n+1}</math> and for <math display="inline">\alpha =1</math>, <math display="inline">\delta p =\Delta p</math>. Note that in both cases the sum of Eqs.([[#eq-42a|42a]]) and ([[#eq-42b|42b]]) gives the time discretization of the momentum equations with the pressures computed at <math display="inline">t^{n+1}</math>. The value of <math display="inline">\bar {\boldsymbol u}^{n+1}</math> from Eq.([[#eq-42b|42b]]) is substituted now into Eq.([[#eq-34b|34b]]) to give
<span id="eq-43a"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>{\boldsymbol G}^T\bar {\boldsymbol u}^* + \Delta t {\boldsymbol G}^T {\boldsymbol M}^{-1}{\boldsymbol G}\delta \bar{\boldsymbol p} + {\boldsymbol L} \bar{\boldsymbol p}^{n+1}+ {\boldsymbol Q}\bar {\boldsymbol \pi }^{n+\theta _4}=0</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (43a)
|}

The product <math display="inline">{\boldsymbol G}^T {\boldsymbol M}^{-1}{\boldsymbol G}</math> can be approximated by a laplacian matrix, i.e.
<span id="eq-43b"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>{\boldsymbol G}^T {\boldsymbol M}^{-1}{\boldsymbol G}=\hat {\boldsymbol L} \quad \hbox{with } \hat{\boldsymbol L} \simeq \int _{\Omega ^e} {1\over \rho } {\boldsymbol \nabla }^T N_i {\boldsymbol \nabla } N_j \,d\Omega  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (43b)
|}
<br />
A semi-implicit algorithm can thus be derived as follows.<br/>
<span id="step-1a"></span>
'''Step 1''' Compute <math display="inline">{\boldsymbol u}^*</math> explicitely from Eq.([[#eq-42a|42a]]) with <math display="inline">{\boldsymbol M}={\boldsymbol M}_d</math> where subscript <math display="inline">d</math> denotes hereonwards a diagonal matrix.

<br/>
<span id="step-2a"></span>
'''Step 2''' Compute <math display="inline">\delta \bar {\boldsymbol p}</math> from Eq.([[#eq-43a|43a]]) as
<span id="eq-44"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\delta \bar {\boldsymbol p} =-({\boldsymbol L}+\Delta t \hat  {\boldsymbol L})^{-1} [{\boldsymbol G}^T\bar{\boldsymbol u}^* +{\boldsymbol Q}\bar {\boldsymbol \pi }^{n+\theta _4}+ \alpha {\boldsymbol L} \bar {\boldsymbol p}^n] </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (44)
|}
<span id="step-3a"></span>
'''Step 3''' Compute <math display="inline"> \bar{\boldsymbol u}^{n+1}</math> explicitely from Eq.([[#eq-42b|42b]]) with <math display="inline">{\boldsymbol M}={\boldsymbol M}_d</math>

<br/>
<span id="step-4a"></span>
'''Step 4''' Compute <math display="inline"> \bar{\boldsymbol c}^{n+1}</math> explicitely from Eq.([[#eq-39|39]]) with
<span id="eq-45"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\bar{\boldsymbol c}^{n+1}=-{\boldsymbol M}_d^{-1}\hat {\boldsymbol C}^{n+1}\bar{\boldsymbol u}^{n+1} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (45)
|}
<span id="step-5a"></span>
'''Step 5''' Compute <math display="inline">  \bar{\boldsymbol \pi }^{n+1}</math> explicitely from Eq.([[#eq-40|40]]) as
<span id="eq-46"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\bar{\boldsymbol \pi }^{n+1}=- \hat {\boldsymbol M}_d^{-1} {\boldsymbol Q}^T \bar {\boldsymbol p}^{n+1} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (46)
|}

This algorithm has an additional step than  the iterative scheme of Section [[#4.1 Quasi-implicit scheme|4.1.]] The advantage is that now Steps [[#step-1a|1]] and [[#step-2a|2]] can be fully linearized by choosing <math display="inline">\theta _1 = \theta _3=\theta _4=0</math>. Also the equation for the pressure variables in Step [[#step-2a|2]] has improved stabilization properties due to the additional laplacian matrix <math display="inline">\hat {\boldsymbol L}</math>.

The boundary conditions are applied as follow. No condition is applied in the computation of the fractional velocities <math display="inline">{\boldsymbol u}^*</math> in Eq.([[#eq-44|44]]). The prescribed velocities at the boundary are applied when solving for <math display="inline">\bar{\boldsymbol u}^{n+1}</math> in the step [[#step-3a|3]]. The prescribed pressures at the boundary are imposed by making zero the pressure increments at the relevant boundary nodes and making <math display="inline">\bar{\boldsymbol p}^n</math> equal to the prescribed pressure values.

===4.3 Stokes flow===

The formulation for a Stokes flow can be readily obtained simply by neglecting the convective terms in the general Navier-Stokes formulation. This also implies neglecting the convective stabilization terms in the momentum equations and, consequently, the convective projection variables are not larger necessary. Also the intrinsic time parameters <math display="inline">\tau _i</math> take now the simpler form (see Eq.([[#eq-22|22]])):
<span id="eq-47"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\tau _i={3h_i^2\over 8\mu } </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (47)
|}

The resulting discretized system of equations can be written as (see Eqs.([[#eq-34|34]]))
<span id="eq-48"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\begin{array}{l}\displaystyle {\boldsymbol M}\dot{\bar {\boldsymbol u}} + {\boldsymbol K}\bar{\boldsymbol u} - {\boldsymbol G}\bar {\boldsymbol p}={\boldsymbol f}\\  \displaystyle {\boldsymbol G}^T \bar {\boldsymbol u} + {\boldsymbol L}\bar {\boldsymbol p}+{\boldsymbol Q}\bar {\boldsymbol \pi }={\boldsymbol 0}\\ \displaystyle {\boldsymbol Q}^T \bar {\boldsymbol p} + \hat {\boldsymbol M}\bar {\boldsymbol \pi }={\boldsymbol 0} \end{array} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (48)
|}

The iterative algorithm of Section [[#Quasi-implicit scheme|4.1]] can now be implemented. Convergence is now faster due to the absence of the non linear convective terms in the momentum equation.

The steady-state form of Eqs.([[#eq-48|48]]) can be expressed in matrix form as
<span id="eq-49"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\left[\begin{matrix}{\boldsymbol K}&-{\boldsymbol G}&{\boldsymbol 0}\\ -{\boldsymbol G}^T & - {\boldsymbol L} &-{\boldsymbol Q}\\ {\boldsymbol 0}& -{\boldsymbol Q}^T&-\hat {\boldsymbol M}\\\end{matrix}\right]\left\{\begin{matrix}\bar {\boldsymbol u}\\ \bar {\boldsymbol p}\\ \bar {\boldsymbol \pi }\\\end{matrix}\right\}= \left\{\begin{matrix}{\boldsymbol f}\\ {\boldsymbol 0}\\{\boldsymbol 0}\\ \end{matrix}\right\} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (49)
|}

The system is symmetric and always positive definite and therefore leads to a non singular solution. We note that this property holds for ''any interpolation function'' chosen for <math display="inline">\bar {\boldsymbol u},\bar {\boldsymbol p}</math> and <math display="inline">\bar {\boldsymbol \pi }</math>, therefore overcoming the Babuŝka-Brezzi (BB) restrictions <span id='citeF-1'></span>[[#cite-1|[1]]].

A reduced velocity-pressure formulation can be obtained by eliminating the pressure gradient projection variables <math display="inline">\bar {\boldsymbol \pi }</math> from the last equation to give
<span id="eq-50"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\left[\begin{matrix}{\boldsymbol K}&-{\boldsymbol G}\\ -{\boldsymbol G}^T & - ({\boldsymbol L} -{\boldsymbol Q}\hat {\boldsymbol M}^{-1}{\boldsymbol Q}^T)\\\end{matrix}\right]\left\{\begin{matrix}\bar {\boldsymbol u}\\ \bar {\boldsymbol p}\\\end{matrix}\right\}= \left\{\begin{matrix}{\boldsymbol f}\\ {\boldsymbol 0}\\\end{matrix}\right\} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (50)
|}

The reduction process is simplified by using a diagonal form of matrix <math display="inline">\hat {\boldsymbol M}</math>. Obviously above reduction is also applicable to the transient case.

==5 FLUID-STRUCTURE INTERACTION. MESH UPDATING. ALE FORMULATION==

===5.1 General coupled solution scheme===

The two algorithms of previous section can be readily extended for fluid-structure interaction analysis. The solution process in both cases includes the two additional steps.

'''Step&nbsp; A1.&nbsp; Solve&nbsp; for&nbsp; the&nbsp; movement&nbsp; of&nbsp; the&nbsp; structure &nbsp;due &nbsp;to&nbsp; the&nbsp; fluid&nbsp; flow &nbsp;forces'''

This implies solving the dynamic equations of motion for the structure written as
<span id="eq-51"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>{\boldsymbol M}_s \ddot {\boldsymbol d}+ {\boldsymbol K}_s {\boldsymbol d}={\boldsymbol f}_{ext} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (51)
|}

where <math display="inline">{\boldsymbol d}</math> and <math display="inline">\ddot {\boldsymbol d}</math> are respectively the displacement and acceleration vectors of the nodes discretizing the structure, <math display="inline">{\boldsymbol M}_s</math> and <math display="inline">{\boldsymbol K}_s</math> are the mass and stiffness matrices of the structure and <math display="inline">{\boldsymbol f}_{ext}</math> is the vector of external nodal forces accounting for the fluid flow forces induced by the pressure and the viscous stresses. Clearly the main driving forces for the motion of the structure is the fluid pressure which acts in the form of a surface traction on the structure. Indeed Eq.([[#eq-51|51]]) can be augmented with an appropriate damping term. The form of all the relevant matrices and vectors can be found in standard books on FEM for structural analysis <span id='citeF-1'></span>[[#cite-1|[1]]].

Solution of Eq.([[#eq-51|51]]) in time can be performed using implicit or fully explicit time integration algorithms. In both cases the values of the nodal displacement, velocities and accelerations at <math display="inline">t^{n+1}</math> are found.

'''Step&nbsp;A2. Compute the new&nbsp;position&nbsp;of&nbsp;the&nbsp;mesh&nbsp;nodes'''

Movement of a structure in a fluid originates a distorsion in the mesh defining the control volume where the fluid equations are solved. Clearly a new mesh can be regenerated at each time step and this option is discussed in a later section dealing with lagrangian flows. A cheaper alternative is to update the position of the mesh nodes once the iterative process for the fluid and solid variables has converged. A simple algorithm for updating the mesh nodes is described in the next section.

===5.2 A simple algorithm for updating the mesh nodes===

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

Chiandussi, Bugeda and Oñate <span id='citeF-45'></span>[[#cite-45|[45]]] have 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 the 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 arbitrary homogeneous isotropic elastic properties characterized by the Young modulus <math display="inline">\bar E</math> and the Poisson coefficient <math display="inline">\nu </math>. Once a discretized finite element problem has been solved using, for instance, standard <math display="inline">C^o</math> linear triangles (in 2D) or linear tetraedra (in 3D), the principal stresses <math display="inline">{}^1\sigma _i</math> at the center of each element can be obtained. For 3D problems
<span id="eq-52"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>{}^1\sigma _i ={\bar E (1-\nu )\over ( 1+\nu )(1-2\nu )}\left[\varepsilon _i + {\nu \over  1-\nu } (\varepsilon _j + \varepsilon _k)\right]\qquad i,j,k=1,2,3 ~~i\not =j\not = k </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (52)
|}

where <math display="inline">\varepsilon _i</math> are the principal strains. As the values of <math display="inline">\bar E</math> and <math display="inline">\nu </math> are arbitrary it is useful to select <math display="inline">\nu =0</math>. Eq.([[#eq-32|32]]) simplifies in this case to <math display="inline">\sigma _i=\bar E \varepsilon _i</math>.

Let us assume now that a uniform strain field <math display="inline">\varepsilon _i =\bar \varepsilon </math> throughout the mesh is sought. The principal stresses are then given by
<span id="eq-53"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>{}^2\sigma _i =  E \bar \varepsilon \qquad i=1,2,3 </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (53)
|}

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

A number of criteria can be now used to find the value of <math display="inline">E</math>. An effective approach found in <span id='citeF-45'></span>[[#cite-45|[45]]] is to make equal the element strain energy densities in both analysis. Thus (for <math display="inline">\nu =0</math>)
<span id="eq-54"></span>
<span id="eq-54a"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>U_1 =  {{}^1\sigma _i \varepsilon _i\over 2} = {\bar E [(\varepsilon _1^2+\varepsilon _2^2 + \varepsilon _3^2)]\over 2} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (54a)
|}
<span id="eq-54b"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>U_2 =  {{}^2\sigma _i \varepsilon _i \over 2}= {3E\bar \varepsilon ^2\over 2}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (54b)
|}

Equaling eqs.(54a) and (54b) gives the sought Young modulus <math display="inline">E</math> as
<span id="eq-55"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>E = {\bar E \over 3 \bar  \varepsilon ^2 }(\varepsilon _1^2+\varepsilon _2^2 + \varepsilon _3^2) </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (55)
|}

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

The solution process includes the following two steps.

'''Step 1'''. Consider the finite element mesh as a linear elastic solid with homogeneous material properties characterized by <math display="inline">\bar E</math> and arbitrary Young modulus <math display="inline">\bar E</math> and the Poisson ratio <math display="inline">\nu =0</math>. Solve the corresponding elastic problem with imposed displacements at the mesh boundary.

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

The movement of the mesh nodes obtained in the second step ensures a quasi uniform mesh distortion. Further details on this method including other alternatives for evaluating the Young modulus <math display="inline">E</math> can be found in <span id='citeF-45'></span>[[#cite-45|[45]]].

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 solution to this problem is to remesh the analysis domain. However, for most  problems, a mapping of the moving surfaces linked to the mesh updating algorithm described above can avoid remeshing <span id='citeF-38'></span>[[#cite-38|[38]],<span id='citeF-46'></span>[[#cite-46|46]]-<span id='citeF-49'></span>[[#cite-49|49]]].

===5.3 ALE formulation===

The movement of the mesh defining the fluid domain requires accounting for the relative motion of the fluid particles with respect to the moving mesh. This can be dealt with by an arbitrary lagrangian-eulerian (ALE) formulation. This basically implies redefining the convective transport term in the momentum equation as

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>v_j {\partial u_i \over \partial x_j}\quad \hbox{with } v_j =u_j - u_j^m </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (56)
|}

where <math display="inline">v_j</math> is the relative velocity between the moving mesh and the fluid point and <math display="inline">u_j^m</math> is the velocity of the mesh nodes. This velocity can be simply computed dividing by <math display="inline">\Delta t</math> the displacement vector of the nodes in the mesh obtained from the mesh updating algorithm previously described.

==6 FREE SURFACE WAVE EFFECTS==

Many problems of practical importance involve a free surface in the fluid. In general the position of such a free surface is unknown and has to be determined. Typical problems of this kind are water flow around ships, flow under and over water control structures, mould filling processes, etc.

On the free surface <math display="inline">\Gamma _\beta </math> we must ensure al all times that (1) the pressure (which approximate the normal traction) equals the atmospheric pressure <math display="inline">p_a</math> and the tangential tractions are zero (unless specific otherwise) and (2) that the material particles of the fluid belong to the free surface.

Condition (1) is simply fulfilled by imposing <math display="inline">p=p_a</math> on <math display="inline">\Gamma _\beta </math> during the solution for the nodal pressures.

The free surface condition (2) can be written in the FIC formulation  (neglecting time stabilization effects) as <span id='citeF-46'></span>[[#cite-46|[46]]-<span id='citeF-49'></span>[[#cite-49|49]]]
<span id="eq-57"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>r_\beta - \underline{{1\over 2} h_{\beta _j} {\partial r_\beta  \over \partial x_j}}=0\quad j=1,2 </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (57)
|}

where
<span id="eq-58"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>r_\beta := {\partial \beta  \over \partial t}+v_i {\partial \beta  \over \partial x_i}-v_3 \quad i=1,2 </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (58)
|}

where <math display="inline">\beta </math> is the wave elevation (measured with respect to a reference surface of height <math display="inline">\beta _{ref}</math>) and <math display="inline">v_i</math> is the relative velocity defined in Eq.([[#eq-56|56]]). The underlined term in Eq.([[#eq-57|57]]) introduces the necessary stabilization for the solution of the highly convective (and non linear) equation defining the evolution of the wave elevation. Note that neglecting the stabilization term,  the steady state form of Eq.([[#eq-57|57]]) <math display="inline">\left(\displaystyle {\partial \beta  \over \partial t}= 0\right)</math> simply states that the fluid particles move in the tangential direction to the free surface (in 2D: <math display="inline">\displaystyle {\partial \beta  \over \partial x}=\displaystyle{u_2\over u_1} =\hbox{tg } \delta </math> where <math display="inline">\delta </math> is the angle which the velocity vector forms with the horizontal axis).

The solution in time of Eq.([[#eq-57|57]]) can be expressed in terms of the nodal velocities computed from the flow solution, as
<span id="eq-59"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\beta ^{n+1} = \beta ^n -\Delta t \left[v_i^{n+1,i} {\partial \beta ^n \over \partial x_i}-v_3^{n+1,i}-{h_{\beta _j}\over 2} {\partial r_\beta ^n \over \partial x_j}\right] </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (59)
|}

Eq.([[#eq-59|59]]) can now be discretized in space using the standard Galerkin method and solved ''explicitely'' to give the nodal wave heights at <math display="inline">t^{n+1}</math> <span id='citeF-46'></span>[[#cite-46|[46]]-<span id='citeF-49'></span>[[#cite-49|49]]]. This solution step should preceed the computation of the structure motion in the case of a fluid-structure interaction problem. Typically the general algorithm will be as follows:

<ol>

<li>Solve for the nodal velocities <math display="inline">\bar{\boldsymbol u}^{n+1}</math> and the pressures <math display="inline">\bar {\boldsymbol p}^{n+1}</math> in the fluid domain using any of the algorithms of Section 4. When solving for the pressure variables impose <math display="inline">p^{n+1}=p_a</math> at the free surface <math display="inline">\Gamma _\beta </math>. </li>
<li>Solve for the free surface elevation <math display="inline">\beta ^{n+1}</math> via Eq.([[#eq-59|59]]). </li>
<li>Compute the movement of the fully or semi-submerged or floating structure by solving the dynamic equations of motion of the structure (Eq.([[#eq-51|51]])). </li>
<li>Compute the new position of the mesh nodes in the fluid domain at time <math display="inline">t^{n+1}</math>. Alternatively, regenerate a new mesh. </li>

</ol>

The mesh updating proces can also include the free surface nodes, although this is not strictly necessary. An ''hydrostatic adjustement'' can be implemented once the new free surface elevation is computed by simple imposing the pressure at the nodes on the reference surface as
<span id="eq-60"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>p^{n+1}=p_a + \rho \vert g\vert \Delta \beta \quad \hbox{with } \Delta \beta = \beta ^{n+1}- \beta _{ref} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (60)
|}

where <math display="inline">g</math> is the  gravity constant. Eq.([[#eq-60|60]]) allows to take into account the changes in the free surface without the need of updating the reference surface nodes. A higher accuracy in the solution of the flow problem can be obtained by updating the reference surface nodes after a number of time steps.

==7 LAGRANGIAN FLOW FORMULATION==

The Lagrangian formulation is an effective (and relatively simple) procedure for modelling the flow of fluid particles undergoing severe distorsions such as water jets, high amplitude waves, water splashing, breaking waves, filling of cavities, etc. Indeed the lagrangian formulation seems to be an excellent procedure for treating fluid-structure interaction problems where the structure has large displacements. An obvious “a priori” advantage of the lagrangian formulation is that both the structure and the fluid motions are defined in the same frame of reference.

The lagrangian fluid flow equations can be simply obtained by noting that the velocity of the mesh nodes and that of the fluid particles are the same. Hence the relative velocity <math display="inline">v_i</math> is zero in Eq.([[#eq-56|56]]) and the convective terms vanish in the momentum equations, while the rest of the fluid flow equations defined in Section [[#3 GENERAL FIC EQUATIONS FOR VISCOUS INCOMPRESSIBLE FLOW|3]] remain unchanged.

The FEM algorithms for solving the lagrangian flow equations is very similar to those for the eulerian or ALE description presented earlier and only the main differences will be given here. For preciseness we will focus in the semi-implicit fractional step algorithm of Section [[#4.2 Fractional step method|4.2]] (for <math display="inline">\theta _1 = \theta _4=0</math> and <math display="inline">\alpha =1</math>) accounting also for fluid-structure interaction effects.

<br/>

'''Step 1''' Compute explicitely a predicted value of the velocities <math display="inline">{\boldsymbol u}^*</math> as
<span id="eq-61"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\bar{\boldsymbol u}^{*} = \bar{\boldsymbol u}^n - \Delta t {\boldsymbol M}^{-1}_d [{\boldsymbol K} \bar{\boldsymbol u}^{n} - {\boldsymbol G}\bar {\boldsymbol p}^{n} -{\boldsymbol f}^{n+1}] </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (61)
|}

'''Step 2''' Compute <math display="inline">\delta \bar {\boldsymbol p}</math> from Eq.([[#eq-44|44]]).<br/>

'''Step 3''' Compute explicitely <math display="inline"> \bar{\boldsymbol u}^{n+1}</math> from Eq.([[#eq-42b|42b]]) with <math display="inline">{\boldsymbol M}={\boldsymbol M}_d</math>.

<br/>

'''Step 4''' Compute <math display="inline"> \bar{\boldsymbol \pi }^{n+1}</math> explicitely from Eq.([[#eq-40|40]]).

<br/>

'''Step 5''' Solve for the motion of the structure by integrating Eq.([[#eq-51|51]]).

<br/>

'''Step 6''' Update the mesh nodes in a lagrangian manner as
<span id="eq-62"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>{\boldsymbol x}_i^{n+1} = {\boldsymbol x}_i^{n}+\bar {\boldsymbol u}_i^{n+1} \Delta t </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (62)
|}

'''Step 7''' Generate a new mesh. This can be effectively performed using the extended Delaunay Tesselation described in <span id='citeF-50'></span>[[#cite-50|[50]]]. Indeed the mesh regeneration can take place after a prescribed number of time steps or when  the nodal displacements induce significant distorsions in some element shapes.

Details of the treatment of the boundary conditions in the lagrangian flow formulation can be found in <span id='citeF-51'></span>[[#cite-51|[51]],<span id='citeF-52'></span>[[#cite-52|52]]].

==8 TURBULENCE MODELLING==

The detailed discussion on the treatment of turbulent effects in the flow equation falls outside the objective of this chapter as any of the existing turbulence model is applicable.

In the examples presented next we have chosen a turbulence model based on the Reynolds averaged Navier-Stokes equations where the deviatoric stresses are computed as sum of the standard viscous contributions and the so called Reynold stresses. Here we have chosen the Boussinesq assumption leading to a modification of the viscosity in the standard Navier-Stokes equations as sum of the “physical” viscosity <math display="inline">\mu </math> and a turbulent viscosity <math display="inline">\mu _T</math>.

One of the simplest and more effective choices for <math display="inline">\mu _T</math> is the Smagorinski LES model giving

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\mu _T=C_l h^e (2\varepsilon _{ij}\varepsilon _{ij})^{1/2} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (63)
|}

where <math display="inline">h^e</math> is the element size and <math display="inline">C</math> is a constant (<math display="inline">C\simeq 0.01</math>).

Indeed other many options are possible such as the one and two equations turbulence models (i.e. the <math display="inline">k</math> model and the <math display="inline">k - \varepsilon </math> and <math display="inline">k - w</math> models) and the algebraic stress models and the reader is refered to specialized books on this matter <span id='citeF-53'></span>[[#cite-53|[53]]].

==9 COMPUTATION OF THE CHARACTERISTIC LENGTHS==

The evaluation of the stabilization parameters is one of the crucial issues in stabilized methods. Most of existing methods use expressions which are direct extensions of the values obtained for the simplest 1D case. It is also usual to accept the so called SUPG assumption, i.e. to admit that vector <math display="inline">{\boldsymbol h}</math> has the direction of the velocity field <span id='citeF-32'></span>[[#cite-32|[32]],<span id='citeF-37'></span>[[#cite-37|37]]]. This unnecessary restriction leads to instabilities when sharp layers transversal to the velocity direction are present. This deficiency is usually corrected by adding a shock capturing or crosswind stabilization term <span id='citeF-54'></span>[[#cite-54|[54]],<span id='citeF-55'></span>[[#cite-55|55]]]. Indeed, in the FIC formulation the components of <math display="inline">{\boldsymbol h}</math> 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  defined by
<span id="eq-X"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>{\boldsymbol h}=h_s {{\boldsymbol u}\over {u}}+h_{c} {{\boldsymbol \nabla } u\over \vert{\boldsymbol \nabla }u\vert } </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (64)
|}

where <math display="inline">u=\vert {\boldsymbol u}\vert </math> and <math display="inline">h_s</math> and <math display="inline">h_{c}</math> are the “streamline” and “cross wind” contributions given by
<span id="eq-65"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>h_s=\max ({\boldsymbol l}^T_j {\boldsymbol u})/{u} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (65)
|}
<span id="eq-66"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math> h_{c}=\max ({\boldsymbol l}^T_j {\boldsymbol \nabla }u)/ \vert {\boldsymbol \nabla }u\vert \quad , \quad  j=1,n_s </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (66)
|}

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

As for the free surface equation the following value of the characteristic length vector <math display="inline">{\boldsymbol h}_\beta </math> has been taken
<span id="eq-67"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>{\boldsymbol h}_\beta =\bar h_s {{\boldsymbol u}\over {u}}+\bar h_c {{\boldsymbol \nabla }\beta \over \vert {\boldsymbol \nabla }\beta \vert } </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (67)
|}

The streamline parameter <math display="inline">\bar h_s</math> has been obtained by Eq.([[#eq-65|65]]) using the value of the velocity vector <math display="inline">\boldsymbol u</math> over the 3 node triangles discretizing the free surface and <math display="inline">n_s=3</math>.

The cross wind parameter <math display="inline">\bar h_c</math> has been computed by
<span id="eq-68"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\bar h_c = \max [{\boldsymbol l}_j^T {\boldsymbol \nabla }\beta ] {1\over \vert {\boldsymbol \nabla }\beta \vert } \quad ,\quad j=1,2,3 </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (68)
|}

The cross-wind terms in eqs.([[#eq-64|64]]) and ([[#eq-67|67]]) 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 <span id='citeF-54'></span>[[#cite-54|[54]],<span id='citeF-55'></span>[[#cite-55|55]]].

A more consistent evaluation of <math display="inline">{\boldsymbol h}</math> based on a diminishing residual technique can be found in <span id='citeF-37'></span>[[#cite-37|[37]]].

==10 EXAMPLES==

The examples chosen show the applicability of the Eulerian, ALE and lagrangian formulations presented to solve fluid flow problems. Linear tetrahedra and triangles have been used in the 3D and 2D analysis shown. The fractional step algorithm of Section [[#4.2 Fractional step method|4.2]] for <math display="inline">\theta _1 = \theta _3=\theta _4=0</math> and <math display="inline">\alpha =1</math> has been used in all cases. The first example is the standard square cavity problem solved in 3D using an Euler formulation (<math display="inline">\mu =0</math>). The second example is the flow past a submerged NACA 0012 profile. Here the free surface equation is solved together with the flow equations.

The next four examples fall within the category of fluid-structure interaction problems. The first is the analysis of a sphere falling in a tube filled with liquid where the mesh updating procedure of Section 5.2 has been used. Then three of ship hydrodynamics problems are solved including the analysis of a Wigley hull, a scale model of a commercial ship and an American Cup racing sail boat. Numerical results are compared with experimental data in all cases.

The last series of examples show applications of the Lagrangian formulation to the simulation of the collapse of a water column, a semi-submerged rotating water mill and a solid cube falling into a water recipient.

===10.1 Square cavity problem===

The purpose of this example is to test the stabilized formulation presented in the solution of a standard benchmark problem <span id='citeF-17'></span>[[#cite-17|[17]]]. Figure [[#img-2|2]] shows the definition of the problem solved with an unstructured 3D mesh of 7395 linear tetrahedra for a Reynolds number value of 1.

<div id='img-2'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Samper_187067199-Fig3acon183.png|251px|]]
|-
|[[Image:Draft_Samper_187067199-Fig3bcon183.png|251px|]]
|-
| colspan="2"|[[Image:Draft_Samper_187067199-Fig3ccon183.png|251px|Square cavity problem. a) Problem definition. b) Unstructured mesh of 7395  linear tetrahedra. c) velocity field for Re = 1.]]
|- style="text-align: center; font-size: 75%;"
| colspan="2" | '''Figure 2:''' Square cavity problem. a) Problem definition. b) Unstructured mesh of 7395  linear tetrahedra. c) velocity field for <math>Re = 1</math>.
|}

Results in Figure [[#img-3|3a,b]] are tabulated for the horizontal velocity along the vertical centerline of the mid-section and for vertical velocity and pressure along the horizontal centerline of the same section. Numerical results are fully stable and agree well with similar solutions reported <span id='citeF-17'></span>[[#cite-17|[17]]]. The effect of the <math display="inline">\tau _i</math> stabilization term  in the pressure equation (see eq.([[#eq-30|30]])) is seen clearly in Figure [[#img-3|3c]]. The curves in this figure show the convergence towards steady state of the <math display="inline">L\infty </math> norm of the nodal pressures with time. The curve listed as “standard” is obtained neglecting the <math display="inline">\tau _i</math> stabilization term, whereas the second curve shows the convergence when this term is taken into account. The difference between the two curves is noticeable as the error obtained with the fully stabilized solution is several orders of magnitude smaller than that obtained neglecting the <math display="inline">\tau _i</math> term.

<div id='img-3'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Samper_187067199-Fig4acon183.png|351px|]]
|-
|[[Image:Draft_Samper_187067199-Fig4bcon183.png|351px|]]
|-
|[[Image:Draft_Samper_187067199-Fig4ccon183.png|375px|Square cavity. a) Distribution of pressure along horizontal centerline of mid-section. b) Distribution of velocity along horizontal centerline of mid-section. c) Convergence histories of the nodal pressure norm (L_∞) for the stabilised (accounting for τ<sub>i</sub>) and the standard (τ<sub>i</sub>=0) schemes.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 3:''' Square cavity. a) Distribution of pressure along horizontal centerline of mid-section. b) Distribution of velocity along horizontal centerline of mid-section. c) Convergence histories of the nodal pressure norm (<math>L_\infty </math>) for the stabilised (accounting for <math>\tau _i</math>) and the standard (<math>\tau _i=0</math>) schemes.
|}

===10.2 Submerged NACA 0012 profile===

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

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

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

===10.3 Sphere falling in a tube filled with water===

The movement of a sphere falling by gravity in a cylindrical tube filled with water 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 [[#img-5|5]]).

Figures [[#img-5|5]] and [[#img-6|6]] 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 6c. 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 greater number of spheres has been solved by Johnson and Tezduyar <span id='citeF-60'></span>[[#cite-60|[60]]].

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

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

===10.4 Wigley hull===

The next problem case considered here is the study of the hydrodynamics of the well known  Wigley Hull.

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

<div id='img-7'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Samper_187067199-figura6a.png|450px|]]
|-
|[[Image:Draft_Samper_187067199-figura6b.png|400px|]]
|-
| colspan="1"|[[Image:Draft_Samper_187067199-figura6c.png|400px|Wigley hull. a) Pressure distribution and mesh deformation of the wigley hull (free model). b) Numerical and experimental body wave profiles. c) Free surface contours for the truly free ship motion.]]
|- style="text-align: center; font-size: 75%;"
| colspan="2" | '''Figure 7:''' 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 [[#img-7|7]] also shows the results of the viscous  analysis of the Wigley model in three different cases for (<math display="inline">L=6m, F=0.316, \mu=10^{-3} Kg/m.s</math>). In the first case the volume mesh was considered fixed, not allowing free surface nor ship movements. Secondly, the volume mesh was updated due to free surface movement, considering the model fixed. The third case corresponds to the analysis of a real free model including the mesh updating due to free surface evaluation and ship movement (sinkage and trim). A Smagorinsky turbulence model was used in  the three cases.

Table [[#table-1|1]] shows the obtained total resistance coefficient in the three cases studied compared with the experimental data.

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

|}

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

Figure 7a 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 7b.

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

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

===10.5 KVLCC2 hull model===

The example is the analysis of the KVLCC2 benchmark model. Here a partially wetted tramsom stern is expected due to the low Froude number of the test. Figure 8 shows the NURBS geometry used obtained from the Hydrodynamic Performance Research team of Korea (KRISO). The obtained results are compared with the experimental data available in the KRISO database <span id='citeF-63'></span>[[#cite-63|63]].

<div id='img-8'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Samper_187067199-KCSgeo.png|100px|KVLCC2 model. Geometrical definition based on NURBS surfaces.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 8:''' KVLCC2 model. Geometrical definition based on NURBS surfaces.
|}

The smallest element size used was 0.001 m and the largest 0.50 m. The surface mesh chosen  is shown in Figure 9. A total of 550.000 tetrahedra were used in the analysis.  The tramsom stern flow model presented  in the previous section was used.

''Test 1.-'' ''Wave pattern calculation''. The main characteristics of the analysis are listed below:

* Length: 5.52 m, Beam (at water plane): 0.82 m, Draught: 0.18 m, Wetted Surface: <math display="inline">8.08 m^2</math>.
* Velocity: 1.05 m/seg, Froude Number: 0.142.
* Viscosity: <math display="inline">0.00126 Kg/m\cdot seg</math>, Density: <math display="inline">1000 Kg/m^3</math>, Reynolds number: <math display="inline">4.63\cdot 10^6</math>.

<div id='img-9'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Samper_187067199-KCSmesh.png|100px|KVLCC2 model. Surface mesh used in the analysis.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 9:''' KVLCC2 model. Surface mesh used in the analysis.
|}

The turbulence model used in this case was the <math display="inline">K</math> model. Figures 10 and 11 show the wave profiles on the hull and in  a cut at y/L = 0.082 obtained in Test 1, compared to the experimental data. The obtained results are quantitatively good close to the hull. A lost of accuracy is observed in the profiles away from the hull. This is probably due to the fact that the element sizes are not small enough in this area.

<div id='img-10'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Samper_187067199-KCSprof1.png|100px|KVLCC2 model. Wave profile on the hull compared to experimental data. Thick line shows numerical results]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 10:''' KVLCC2 model. Wave profile on the hull compared to experimental data. Thick line shows numerical results
|}

<div id='img-11'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Samper_187067199-KCSprof2.png|100px|KVLCC2 model. Wave profile on a cut at y/L=0.0964 compared to experimental data <span id='citeF-26'></span>[[#cite-26|26]]. Thick line shows numerical results]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 11''' 
|}

''Test 2.-'' ''Wake analysis at different planes''. Several turbulence models were used (''Smagorinsky'', <math display="inline">K</math> and <math display="inline">K-\epsilon </math> model) in order to verify the quality of the results. Here, only the results from the <math display="inline">K-\epsilon </math> model are shown. We note that the velocity maps obtained even for the simplest <math display="inline">Smagorinsky</math> model were qualitatively good, showing the accuracy of the fluid solver scheme used. The main characteristics of this analysis are listed below:

* Length: 2.76 m, Beam (at water plane): 0.41 m, Draught: 0.09 m, Wetted Surface: <math display="inline">2.02 m^2</math>.
* Velocity: 25 m/seg.
* Viscosity: <math display="inline">3.05\cdot 10^{-5} Kg/m\cdot seg</math>, Density: <math display="inline">1.01 Kg/m^3</math>, Reynolds number: <math display="inline">4.63\cdot 10^6</math>.

Figures 12&#8211;13 present results corresponding to the test 2. Figure 12  shows the contours of the axial (X) component of the velocity on a plane at 2.71 m  from the orthogonal aft. Figure 13 shows the maps of the kinetic energy on this plane. Experimental results are shown for comparison in all cases. Further results for this problem can be found in <span id='citeF-48'></span>[[#cite-48|48]].

<div id='img-12'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Samper_187067199-KCSvelx1.png|100px|KVLCC2 model. Map of the X component of the velocity on a plane at 2.71 m from the orthogonal aft. Experimental results shown in the right figure.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 12:''' KVLCC2 model. Map of the X component of the velocity on a plane at 2.71 m from the orthogonal aft. Experimental results shown in the right figure.
|}

<div id='img-13'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Samper_187067199-KCSk.png|100px|KVLCC2 model. Map of the eddy kinetic energy (K) on a plane at 2.71 m from the orthogonal aft. Experimental data shown in the right figure.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 13:''' KVLCC2 model. Map of the eddy kinetic energy (<math>K</math>) on a plane at 2.71 m from the orthogonal aft. Experimental data shown in the right figure.
|}

===10.6 American Cup BRAVO ESPAA Model===

The next example is the analysis of the Spanish American Cup racing sail boat ''Bravo España''. The finite element mesh used is shown in Figure 14. The results presented in Figures&nbsp;14&#8211;17 correspond to the analysis of a non symmetrical case including appendages. Good comparison between the experimental data and the numerical results was again obtained.

Other results of the hydrodynamic analysis of American Cup racing boats carried out with the FEM formulation presented in the paper can be seen in <span id='citeF-64'></span>[[#cite-64|64]].

<div id='img-14'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Samper_187067199-ACmesh.png|100px|Bravo ~Espãna sail racing boat. Mesh used in the analysis.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 14:''' <math>Bravo ~Espa\tilde{n}a</math> sail racing boat. Mesh used in the analysis.
|}

<div id='img-15'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Samper_187067199-ACres1.png|100px|Bravo ~Espãna. Velocity contours.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 15:''' <math>Bravo ~Espa\tilde{n}a</math>. Velocity contours.
|}

<div id='img-16'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Samper_187067199-ACres2.png|100px|Bravo ~Espãna. Streamlines.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 16:''' <math>Bravo ~Espa\tilde{n}a</math>. Streamlines.
|}

<div id='img-17'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Samper_187067199-ACgraph.png|100px|Bravo~Espãna. Resistance test. Comparison of numerical results with experimental data.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 17:''' <math>Bravo~Espa\tilde{n}a</math>. Resistance test. Comparison of numerical results with experimental data.
|}

===10.7 Lagrangian flow===

The first problem solved with the lagrangian formulation is the study of the collapse of a water column. This problem was solved by Koshizu and Oka <span id='citeF-65'></span>[[#cite-65|65]] both experimentally and numerically. It has became a classical example to test the validation of the lagrangian formulation in fluid flows. The water is initially located on the left supported by a removable board. The collapse starts at time t = 0, when the removable board is removed. Viscosity and surface tension are neglected. Figure&nbsp;18 shows the point positions at different time steps. The dark points represent the free-surface detected with an alpha-shape algorithm <span id='citeF-51'></span>[[#cite-51|51]]. The internal points are gray and the fixed points are black.

The water is running on the bottom wall until, near 0.3 sec, it impinges on the right vertical wall. Breaking waves appear at 0.6 sec. Around 1 sec. the water reaches the left wall. Agreement with the experimental results of <span id='citeF-65'></span>[[#cite-65|65]] both in the shape of the free surface as well as in the time evolution are excellent.

<div id='img-18'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Samper_187067199-Fig1ArtBathe.png|351px|Water column collapse at different time steps.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 18:''' Water column collapse at different time steps.
|}

<div id='img-19'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Samper_187067199-Fig1contBathe.png|351px|cont.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 19:''' cont.
|}

<div id='img-20'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Samper_187067199-Fig4ArtBathe.png|351px|Rotating water mill.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 20:''' Rotating water mill.
|}

<div id='img-21'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Samper_187067199-Fig5ArtBathe.png|351px|Solid cube falling into a recipient with water.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 21:''' Solid cube falling into a recipient with water.
|}

The second example is the analysis of a rotating water mill semi submerged in water. A schematic representation of a water mill is presented in Figure 19. The blades of the mill have an imposed rotating velocity, while the water is initially in a stationary and flat position. Fluid structure interactions with free-surfaces and water fragmentation are well reproduced in this example.

The last example represents a free cube falling down into a recipient full of water. The solid cube was modeled by introducing a high viscosity parameter in the elements in the following way: all the polyhedral elements formed by nodes contained in the cube have a high viscosity value. The other elements are inviscid. The results of Figure 20 represent correctly the contact problem when the cube hits the water and also the  speed during the sinking  process.

More examples showing applications of the Lagrangian formulation previously described can be found in <span id='citeF-51'></span>[[#cite-51|51]].

==11 CONCLUSIONS==

The finite calculus form of the fluid mechanics equations  is a good starting point for deriving stabilized finite algorithms for solving a variety of fluid flow problems using Euler, ALE and fully lagrangian descriptions. Both monolithic and fractional step algorithms with intrinsic stabilization properties can be readily derived as shown here. Free surface wave effects and fluid-structure interaction situations  can be accounted for in a straight forward manner within the general flow solution schemes. The ALE formulation is particularly adequate for analysis of problems involving free surface waves of moderate amplitude typical of ship hydrodynamics situations. The lagrangian formulation allows to solve in an effective manner fluid flow problems involving large motions of the free surface and complex fluid-structure interactions.

==ACKNOWLEDGEMENTS==

The authors are grateful to Copa America Desafio Español SA for providing the geometry and experimental data of the racing boat analyzed.

Examples 10.1&#8211;10.7 were analyzed with the finite element code Tdyn based on the FEM formulation here presented <span id='citeF-66'></span>[[#cite-66|66]].

Thanks are also given to Dr. Roberto Flores for many useful discussions.

==References==

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

<div id="cite-2"></div>
[[#citeF-2|[2]]] C. Hirsch, ''Numerical computation of internal and external flow'', J. Wiley, Vol. 1 1988, Vol. 2, 1990.

<div id="cite-3"></div>
[3] A.J. Chorin,  “A numerical solution for solving incompressible viscous flow problems” ''J. Comp. Phys.'', 2, 12&#8211;26, 1967.

<div id="cite-4"></div>
[[#citeF-4|[4]]] W.R. Briley, S.S. Neerarambam and D.L. Whitfield, “Multigrid algorithm for three-dimensional incompressible high-Reynolds number turbulent flows”, ''AIAA Journal'',  33 (1), 2073&#8211;2079, 1995.

<div id="cite-5"></div>
[5] J. Peraire, K. Morgan  and J. Peiro, “The simulation of 3D incompressible flows using unstructured grids”, In ''Frontiers of Computational Fluid Dynamics'', Caughey DA and Hafez MM. (eds), Chapter 16, J. Wiley, 1994.

<div id="cite-6"></div>
[6]  C. Sheng,  L.K. Taylor and D.L. Whitfield,  “Implicit lower-upper/approximate-factorization schemes for incompressible flows” ''Journal of Computational Physics'', 128 (1), 32&#8211;42, 1996.

<div id="cite-7"></div>
[[#citeF-7|[7]]] M. Storti, N. Nigro  and S.R. Idelsohn, “Steady state incompressible flows using explicit schemes with an optimal local preconditioning”, ''Computer Methods in Applied Mechanics and Engineering'', 124, 231&#8211;252, 1995.

<div id="cite-8"></div>
[8] J.C. Heinrich, P.S. Hayakorn and O.C. Zienkiewicz, “An upwind finite element scheme for two dimensional convective transport equations”, ''Int. J. Num. Meth. Engng.'', 11, 131&#8211;143, 1977.

<div id="cite-9"></div>
[9] S.R. Idelsohn, M. Storti and N. Nigro, “Stability analysis of mixed finite element formulation with special mention to stabilized equal-order interpolations” ''Int. J. for Num. Meth. in Fluids'', 20, 1003-1022, 1995.

<div id="cite-10"></div>
[10] S.R. Idelsohn, N. Nigro, M. Storti and G. Buscaglia, “A Petrov-Galerkin formulation for advection-reaction-diffusion”, ''Comput. Meth. Appl. Mech. Engrg.'', 136, 27&#8211;46, 1996.

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

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

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

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

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

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

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

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

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

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

<div id="cite-21"></div>
[[#citeF-21|[21]]] J. Donea, “A Taylor-Galerkin method for convective transport problems”, ''Int. J. Num. Meth. Engng.'', 20, 101&#8211;119, 1984.

<div id="cite-22"></div>
[[#citeF-22|[22]]] J. Douglas, T.F. Russell, “Numerical methods for convection dominated diffusion problems based on combining the method of characteristics with finite element or finite difference procedures”, ''SIAM J. Numer. Anal.'', 19, 871, 1982.

<div id="cite-23"></div>
[23] O. Pironneau, “On the transport-diffusion algorithm and its applications to the Navier-Stokes equations”, ''Numer. Math.'', 38, 309, 1982.

<div id="cite-24"></div>
[24] R. Löhner, K. Morgan, O.C. Zienkiewicz, “The solution of non-linear hyperbolic equation systems by the finite element method”, ''Int. J. Num. Meth. in Fluids'', 4, 1043, 1984.

<div id="cite-25"></div>
[[#citeF-25|[25]]]  R. Codina, M. Vazquez and O.C. Zienkiewicz,  “A general algorithm for compressible and incompressible flow - Part III. The semi-implicit form” ''Int. J. Num. Meth. in Fluids'', 27, 13&#8211;32, 1998.

<div id="cite-26"></div>
[[#citeF-26|[26]]]  R. Codina and O.C. Zienkiewicz, “CBS versus GLS stabilization of the incompressible Navier-Stokes equations and the role of the time step as stabilization parameter”, ''Communications in Numerical Methods in Engineering'', 18 (2), 99&#8211;112, 2002.

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

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

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

<div id="cite-30"></div>
[30]  R. Codina, “Stabilized finite element approximation of transient incompressible flows using orthogonal subscales”, ''Comput. Methods Appl. Mech. Engrg.'', 191, 4295&#8211;4321, 2002.

<div id="cite-31"></div>
[[#citeF-31|[31]]]  J. Donea and A. Huerta, “Finite element method for flow problems”, J. Wiley, 2003.

<div id="cite-32"></div>
[[#citeF-32|[32]]] E. Oñate, Derivation of stabilized equations for advective-diffusive transport and fluid flow problems, ''Comput. Meth. Appl. Mech. Engng.'', Vol. 151, pp. 233&#8211;267, (1998).

<div id="cite-33"></div>
[[#citeF-33|[33]]]  E. Oñate, J. García and S. Idelsohn, ''Computation of the stabilization parameter for the finite element solution of advective-diffusive problems'', Int. J. Num. Meth. Fluids, Vol. 25, pp. 1385&#8211;1407, (1997).

<div id="cite-34"></div>
[34]  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 Comp. Met. in Mech., P. Ladeveze and J.T. Oden (Eds.), Elsevier, (1998).

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

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

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

<div id="cite-38"></div>
[[#citeF-38|[38]]]  E. Oñate and J. García, “A finite element method for  fluid-structure interaction with surface waves using a finite calculus formulation”, ''Comput. MethodsAppl. Mech. Engrg.'', 191, 635&#8211;660, (2001).

<div id="cite-39"></div>
[39]  E. Oñate, “Possibilities of finite calculus in computational mechanics”, Submitted to ''Int. J. Num. Meth. Engng.'', 2002.

<div id="cite-40"></div>
[[#citeF-40|[40]]] E. Oñate and S. Idelsohn, A mesh free finite point method for advective-diffusive transport and fluid flow problems,, ''Computational Mechanics'', 21, 283&#8211;292, 1988.

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

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

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

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

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

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

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

<div id="cite-48"></div>
[[#citeF-48|[48]]] J. García and E. Oñate, “An unstructured finite element solver for ship hydrodynamic problems”, Accepted for publication in ''J. Appl. Mech.'', 2002.

<div id="cite-49"></div>
[49]  E. Oñate, J. García and S.R. Idelsohn, “Ship hydrodynamics”, in ''Encyclopedia of Computational Mechanics'', E. Stein, R. de Borst and T.J.R. Hughes (Eds), J. Wiley, 2004.

<div id="cite-50"></div>
[[#citeF-50|[50]]] S.R. Idelsohn, E. Oñate, N. Calvo and F. del Pin, “The meshless finite element method”,  Accepted ofr publication in ''Int. J. Num. Meth. Engng.'', 2002.

<div id="cite-51"></div>
[[#citeF-51|[51]]] S.R. Idelsohn, E. Oñate, F. Del Pin and N. Calvo, “Lagrangian formulation: the only way to solve some free-surface fluid mechanics problems”, ''Fith World Congress on Computational Mechanics'', Mang HA, Rammerstorfer FG and Eberhardsteiner J. (eds), July 7&#8211;12, Viena, Austria, 2002, Web...

<div id="cite-52"></div>
[52] S.R. Idelsohn, E. Oñate and F. Del Pin, “A lagrangian meshless finite element method applied to fluid-structure interaction problems”, Submitted to ''Computer and Structures'', 2002.

<div id="cite-53"></div>
[[#citeF-53|[53]]] D.C. Wilcox, ''Turbulence modeling for CFD'', DCW Industries Inc., 1994.

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

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

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

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

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

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

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

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

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

<div id="cite-63"></div>
[[#citeF-63|[63]]]  Korea Research Institute of Ships and Ocean Engineering (KRISO).http://www.iihr.uiowa.edu/gothenburg2000/KVLCC/tanker.html.

<div id="cite-64"></div>
[[#citeF-64|[64]]] J. García, R. Luco-Salman, M. Salas, M. López-Rodríguez and E. Oñate E, ``An advanced finite element method for fluid-dynamic analysis of America's Cup boat'', ''High Performance Yatch Design Conference'', Auckland, 4&#8211;6, December, 2002.

<div id="cite-65"></div>
[[#citeF-65|[65]]]  S. Koshizuka and Y. Oka, “Moving particle semi-implicit method for fragmentation of incompressible fluid”, ''Nuclear Engineering Science'', 123, 421-434, 1996.

<div id="cite-66"></div>
[[#citeF-66|[66]]]  ''Tdyn. A finite element code for fluid-dynamic analysis'', COMPASS Ingeniería y Sistemas SA,  <u>www.compassis.com</u>, (2002).