Finite calculus formulations for finite element analysis of incompressible flows. Eulerian, ALE and Lagrangian approaches

==Abstract==

We present a general formulation for incompressible fluid flow analysis using the finite element method (FEM). The standard Eulerian formulation is described first. The necessary stabilization for dealing with convective effects and the incompressibility condition are introduced via the Finite Calculus (FIC) method.  A simple extension of the fluid flow equations to an arbitrary Lagrangian-Eulerian (ALE) frame adequate for treating fluid-structure interaction problems is briefly presented. A fully Lagrangian formulation called the Particle Finite Element Method (PFEM) is also described. The PFEM is particularly attractive for fluid-structure interaction problems involving large motions of the free surface and breaking waves. Examples of application of the Eulerian, the ALE and the fully lagrangian PFEM formulations are presented.

'''Keywords''': Stabilized formulation, incompressible fluid, finite calculus, finite element method, particle finite element method, fluid-structure interaction, lagrangian flow.

==1 INTRODUCTION==

The development of efficient and robust numerical methods for analysis of incompressible flows has been a subject of intensive research in the 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 condition.

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]]] and preconditioning techniques <span id='citeF-7'></span>[[#cite-7|[7]]]. Many FEM schemes for fluid flow analysis 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]]]. A general class of Galerkin FEM has been recently developed where the standard Galerkin variational form of the momentum and mass balance equations is extended with adequate residual-based terms in order to achieve a stabilized numerical scheme. Among the many finite element methods of this kind  we can name the Streamline Upwind Petrov Galerkin (SUPG) method <span id='citeF-1'></span>[[#cite-1|[1]]], the Galerkin Least Square (GLS) method <span id='citeF-19'></span>[[#cite-19|[19]]], the Taylor-Galerkin method <span id='citeF-21'></span>[[#cite-21|[21]]], the Characteristic Galerkin method <span id='citeF-22'></span>[[#cite-22|[22]]] and its variant the Characteristic Based Split (CBS) method <span id='citeF-25'></span>[[#cite-25|[25]]], the pressure gradient operator method <span id='citeF-27'></span>[[#cite-27|[27]]] and the Subgrid Scale (SS) method <span id='citeF-28'></span>[[#cite-28|[28]]]. 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 is 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 momentum and mass balance 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 FEM. The resulting stabilized formulation overcomes the problems due to the convective and incompressibility terms and it allows  the use of low order finite elements (such as linear triangles and tetrahedra) with equal order approximations for the velocity and pressure variables <span id='citeF-33'></span>[[#cite-33|33]]. Application of the FIC approach to the solution of fluid flow problems using a meshless methods are presented in <span id='citeF-42'></span>[[#cite-42|[42]]].

The layout of the paper 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 formulated in an Eulerian frame are presented. The finite element discretization is introduced and the resulting matrix equations are detailed. A fractional step scheme and a predictor-corrector scheme for the transient solution are presented.

The Eulerian formulation is extended to account for free surface wave effects by using an arbitrary Eulerian-Lagrangian (ALE) description and introducing the free surface boundary conditions. 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 to follow the displacement of the mesh nodes accordingly to  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 paper  the fully Lagrangian formulation for fluid flow analysis is presented as a particular case of the ALE form. More specifically, in this paper we describe an innovative Lagrangian formulation  to solve problems involving the interaction between fluids and solids in a unified manner. The procedure, called the ''Particle Finite Element Method'' (PFEM) <span id='citeF-44'></span>[[#cite-44|[44]]], treats the mesh nodes in the fluid and solid domains as dimensionless particles which can freely move an even separate from the main fluid domain representing, for instance, the effect of water drops. A finite element mesh connects the nodes defining the discretized domain where the governing equations are solved in the standard FEM fashion. The PFEM  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.

The main  advantage of the Lagrangian flow formulation is that the convective terms do not enter into the fluid equations. The difficulty is however transferred to the problem of adequately (and efficiently) moving the mesh nodes.  We use a  mesh regeneration procedure blending elements of different shapes based on an extended Delaunay tesselation <span id='citeF-49'></span>[[#cite-49|[49]]].

The examples show the efficiency of the Eulerian, ALE and PFEM formulations to solve classical fluid flow problems, as well as more complex fluid-structure interaction problems 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 1) is written as

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

where <math display="inline">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">u</math> is the velocity and <math display="inline">k</math> is the diffusitivity of the material.

<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_881612373-majesus2.png|316px|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
|}


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 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 after simplification

{| 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}}\frac{h}{2} \frac{d^2q}{dx^2}=0 </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (2)
|}

where <math display="inline">h=d_1-d_2</math> and all the 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>. The new balance equation (2) derived assuming that the balance domain has a finite size, incorporates now the underlined term which introduces the ''characteristic length'' <math display="inline">h</math>. Distance <math display="inline">h</math> can be interpreted as a free parameter depending on the location of point <math display="inline">C</math> (note that <math display="inline">-d\le h \le d</math>). Eq.(2) is the starting point to derive numerical schemes with enhanced stability properties simply by computing the characteristic length parameter from an adequate “optimality” rule.

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

{| 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;" | (3)
|}

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 [1,2,8,32, 36]. In the discretized problem, the characteristic length <math display="inline">h</math> is typically expressed as a function of the cell or element dimensions. The critical and optimal values 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]]].

Equation (3) can be extended to account for source effects. The FIC  equation can be then written in compact form as

{| 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}}\frac{h}{2} \frac{dr}{dx}=0 </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (4)
|}

with

{| 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;" | (5)
|}

where <math display="inline">Q</math> is the external source. A consistent  “finite” form of the Neumann boundary condition  can be  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 in the FIC method is <span id='citeF-32'></span>[[#cite-32|[32]]]

{| 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}\frac{h}{2}r=0 \quad \hbox{at } \, \Gamma _q </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (6)
|}

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.(4) and (6) 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]]].

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

The starting point in the next section are the FIC equations for a viscous incompressible fluid accounting for space stabilization terms only.

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

The FIC governing equations for a viscous incompressible fluid can be written in an Eulerian frame of reference as

'''Momentum'''

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>r_{m_i} - \underline{{1\over 2} h_j{\partial r_{m_i} \over \partial x_j}}{1\over 2} h_j{\partial r_{m_i} \over \partial x_j}=0 \qquad \hbox{in }\Omega  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (7)
|}

'''Mass balance'''

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

where

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>r_{m_i} = \rho \left({\partial u_i \over \partial t}+u_j{\partial u_i \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;" | (9)
|-
| 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;" | (10)
|}
|}

Above <math display="inline">\Omega </math> is the analysis domain, <math display="inline">n_d</math> is the number of space dimensions (<math display="inline">n_d=2</math> for 2D problems), <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

{| 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;" | (11)
|}

where <math display="inline">\delta _{ij}</math> is the Kronecker delta and the strain rates <math display="inline">\dot \varepsilon _{ij}</math> are

{| 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;" | (12)
|}

The FIC boundary conditions  are

{| 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}}{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;" | (13)
|}

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

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

Summation convention for repeated indices in products and derivatives is used unless otherwise specified.

In Eqs.(13) and (14) <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.(13) is positive due to the definition of <math display="inline">r_{m_i}</math> in Eq.(9).

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.(13) these lengths define the domain where equilibrium of boundary tractions is established <span id='citeF-32'></span>[[#cite-32|[32]]]. In the discretized problem the <math display="inline">h_{i,s}</math> become of the order of a typical element dimension as described in Section 9. Note that by making <math display="inline">h_i=0</math> the standard infinitessimal form of the fluid mechanics equations is recovered <span id='citeF-1'></span>[[#cite-1|[1]]].

Eqs.(7)&#8211;(14) are the starting  point for deriving stabilized FEM  for solving the incompressible Navier-Stokes equations. The underlined FIC terms in Eq.(7) are essential to overcome the numerical instabilities due to the convective terms in the momentum equations, whereas the underlined terms in Eq.(8) take care of the instabilities due to the incompressibility constraint. An interesting feature of the FIC formulation is that it allows to use  equal order interpolation for the velocity and pressure variables <span id='citeF-37'></span>[[#cite-37|[37]]].

===3.1 Stabilized integral forms===

From the momentum equations it can be obtained <span id='citeF-37'></span>[[#cite-37|[37]]]

{| 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}\simeq {h_j\over 2a_i} {\partial r_{m_i} \over \partial x_j}\quad ,\quad \hbox{no sum in }i </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (15)
|}

where

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>a_i = {2\mu \over 3} +{\rho u_i h_i\over 2}\quad ,\quad \hbox{no sum in }i </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (16)
|}

Substituting Eq.(15) into Eq.(8) 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 alternative expression for the stabilized mass balance equation

{| 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;" | (17)
|}

with

{| 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;" | (18)
|}

The <math display="inline">\tau _i</math>'s in Eq.(18) when scaled by the density are termed in the stabilization literature ''intrinsic time parameters''. The interest of Eq.(17) is that it introduces the  first space derivatives of the momentum equations (and the corresponding laplacian of pressure terms) into the mass balance equation. These terms have intrinsic good stability properties as explained next.

The weighted residual form of the momentum and mass balance equations (Eqs.(7) and (17)) is written as

{| 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;" | (19)
|}

{| 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;" | (20)
|}

where <math display="inline">\delta u_i</math> and <math display="inline">q</math> are arbitrary weighting functions representing virtual velocities and virtual pressure fields. Integrating by parts the <math display="inline">r_{m_i}</math> terms  leads to

{| 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 + \int _{\Omega } {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;" |  (21a)
|}

{| 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;" |  (21b)
|}

We will neglect hereonwards the third integral in Eq.(21b) 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.(21a) are integrated by parts in the usual manner. The resulting momentum and mass balance equations are

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\begin{array}{r} \displaystyle \int _\Omega \left[\delta u_i\rho \left({\partial u_i \over \partial t}+u_j {\partial {u_i} \over \partial x_j}\right)+ {\partial \delta u_i \over \partial x_j}\left(\mu {\partial u_i \over \partial x_j}-\delta _{ij}p \right)\right]d\Omega -  \int _{\Omega } \delta u_i b_i d\Omega - \int _{\Gamma _t} \delta u_i t_id\Gamma +\qquad \\ \displaystyle + \int _{\Omega } {h_j\over 2}{\partial \delta u_i \over \partial x_j} r_{m_i} d\Omega =0\qquad  \end{array}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" |  (22a)
|}

{| 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;" |  (22b)
|}

In the derivation of the viscous term in Eq.(22a) we have used the following identity (prior to the integration by parts)

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>{\partial s_{ij} \over \partial x_j}=2\mu {\partial \varepsilon _{ij} \over \partial x_j}=\mu {\partial ^2u_i\over \partial x_j \partial x_j} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (23)
|}

Eq.(23) is identically true for the exact incompressible limit <math display="inline">\left({\partial u_i \over \partial x_i}=0\right)</math>.

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

The computation of the residual terms can be simplified if we introduce  the convective and pressure gradient projections <math display="inline">c_i</math> and <math display="inline">\pi _i</math>, respectively defined as

{| 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;" | (24)
|}

We can express <math display="inline">r_{m_i}</math> in  Eqs.(22a) and (22b) 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.(24). This gives the final system of governing equation as:

{| 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)+ {\partial \delta u_i \over \partial x_j}\left(\mu {\partial u_i \over \partial x_j}-\delta _{ij}p \right)\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> + \int _{\Omega } {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\qquad  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (25)
|}

{| 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;" | (26)
|}

{| 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;" | (27)
|}

{| 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;" | (28)
|}

with <math display="inline">i,j,k=1,n_d</math>.  In Eqs.(27) and (28) <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.

We note that accounting for the convective and pressure gradient projections enforces the consistency of the formulation as it ensures that the stabilization terms in Eqs.(25&#8211;28) have a residual form which vanishes for the “exact” solution. Neglecting these terms can reduce the accuracy of the numerical solution and it makes the formulation more sensitive to the value of the stabilization parameters as shown in Example 10.1 and in references <span id='citeF-60'></span>[[#cite-60|[60]]].

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

{| 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 =  N^k \bar u_i^k \quad , \quad p = N^k \bar p^k\\ \displaystyle c_i = N^k \bar c_i^k \quad , \quad \pi _i =  N^k \bar \pi _i^k \end{array} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (29)
|}

where the sum goes over the number of nodes of each element <math display="inline">n</math> (<math display="inline">n=3/4</math> for triangles/tetrahedra), <math display="inline">\bar {(\cdot )}^k</math> denotes nodal variables and <math display="inline">N^k</math> are the linear shape functions <span id='citeF-1'></span>[[#cite-1|[1]]].

Substituting the approximations (29) into Eqs.(25&#8211;28) 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

{| 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 H} \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;" | (30a)
|}

{| 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} + \hat{\boldsymbol L}\bar {\boldsymbol p}+{\boldsymbol Q}\bar {\boldsymbol \pi }={\boldsymbol 0}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (30b)
|}

{| 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;" | (30c)
|}

{| 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;" | (30d)
|}

where

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>{\boldsymbol H}={\boldsymbol A}+{\boldsymbol K}+\hat {\boldsymbol K} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (31)
|}

If we denote the node indexes with superscripts <math display="inline">a,b</math>, the space indices with subscripts <math display="inline">i,j</math>, the element contributions to the components of the arrays involved in these equations are

{| 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}^{ab}= \left(\int _{\Omega ^e} \rho N^a N^b d\Omega \right)\delta _{ij} \quad ,\quad A_{ij}^{ab}= \left(\int _{\Omega ^e} \rho N^a ({\boldsymbol u}^T {\nabla } N^b) d\Omega \right)\delta _{ij}\\ \\ \displaystyle {K}_{ij}^{ab} = \left(\int _{\Omega ^e} \mu {\boldsymbol \nabla }^T N^a{\boldsymbol \nabla }N^b  d\Omega \right)\delta _{ij} \quad ,\quad {\boldsymbol \nabla } = \left[{\partial  \over \partial x_1},{\partial  \over \partial x_2},{\partial  \over \partial x_3}\right]^T\\ \\ \displaystyle \hat{K}_{ij}^{ab} = \left({1\over 2} \int _{\Omega ^e} ({\boldsymbol h}^T {\boldsymbol \nabla } N^a)(\rho {\boldsymbol u}^T {\boldsymbol \nabla } N^b)d\Omega \right)\delta _{ij}\quad ,\quad {G}_{i}^{ab}=\int _{\Omega ^e} {\partial N^a \over \partial x_i}N^b d\Omega \\ \\ \displaystyle {\boldsymbol C}= \left[\begin{matrix}{\boldsymbol C}_1\\ {\boldsymbol C}_2\\ {\boldsymbol C}_3\\\end{matrix}\right]\quad ,\quad {C}_1^{ab}={C}_2^{ab}={C}_3^{ab} = {1\over 2}  \int _{\Omega ^e} [{\boldsymbol h}^T {\boldsymbol \nabla } N^a]N^bd\Omega \\ \\ \displaystyle \hat L^{ab}= \int _{\Omega ^e} ({\boldsymbol \nabla }^T N^a) [\tau ] {\boldsymbol \nabla } N^b d\Omega \quad ,\quad [\tau ]= \left[\begin{matrix}\tau _1 &0 &0 \\ 0 & \tau _2&0\\ 0&0& \tau _3\\\end{matrix}\right]\\ \\ \displaystyle {\boldsymbol Q}= [{\boldsymbol Q}_1,{\boldsymbol Q}_2,{\boldsymbol Q}_3] \quad ,\quad  \displaystyle Q_{i}^{ab} = \int _{\Omega ^e}\tau _i {\partial N^a \over \partial x_i} N^b d\Omega \quad \quad \hbox{no sum in }i \end{array}</math>
|}
|}

{| 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 \hat{\boldsymbol C}= [\hat{\boldsymbol C}_1,\hat{\boldsymbol C}_2,\hat{\boldsymbol C}_3] \quad ,\quad  \displaystyle \hat{C}_1^{ab}=\hat{C}_2^{ab}=\hat{C}_3^{ab} =  \int _{\Omega ^e} \rho ^2 N^a ({\boldsymbol u}^T {\boldsymbol \nabla }N^b)d\Omega \\ \\ \displaystyle  \hat {M}^{ab}_{ij}= \left(\int _{\Omega ^e} \tau _i N^a N^b d\Omega \right)\delta _{ij}\quad ,\quad  \displaystyle {f}_i^a = \int _{\Omega ^e} N^a f_i d\Omega + \int _{\Gamma ^e} N^a t_i d\Gamma  \end{array}</math>
|}
|}

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | 
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (32)
|}

It is understood that all the arrays are matrices (except <math display="inline">{\boldsymbol f}</math>, which is a vector) whose components are obtained by grouping together the left indices in the previous expressions (<math display="inline">a</math> and possibly <math display="inline">i</math>) and the right indices (<math display="inline">b</math> and possibly <math display="inline">j</math>).

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

===4.1 Transient solution scheme===

The solution in time  of the system of Eqs.(30) can be written in general form as

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>{\boldsymbol M} \displaystyle{1\over \Delta t} (\bar{\boldsymbol u}^{n+1}-\bar{\boldsymbol u}^{n}) + {\boldsymbol H}^{n+\theta }\bar{\boldsymbol u} ^{n+\theta } - {\boldsymbol G}\bar{\boldsymbol p}^{n+\theta }+ {\boldsymbol C}^{n+\theta }\bar {\boldsymbol c}^{n+\theta }={\boldsymbol f}^{n+\theta }</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" |  (33a)
|}

{| 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}^{n+\theta }+\hat {\boldsymbol L}^{n+\theta } \bar {\boldsymbol p}^{n+\theta }+{\boldsymbol Q}\bar {\boldsymbol \pi }^{n+\theta }={\boldsymbol 0}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" |  (33b)
|}

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\hat {\boldsymbol C}^{n+\theta } \bar {\boldsymbol u}^{n+\theta }+ {\boldsymbol M} \bar {\boldsymbol c}^{n+\theta }={\boldsymbol 0}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" |  (33c)
|}

{| 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 p}^{n+\theta }+\hat {\boldsymbol M}^{n+\theta }\bar {\boldsymbol \pi }^{n+\theta }={\boldsymbol 0}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" |  (33d)
|}

where <math display="inline">{\boldsymbol H}^{n+\theta }={\boldsymbol H} ({\boldsymbol u}^{n+\theta })</math>, etc and the parameter <math display="inline">\theta \in [0,1]</math>. The direct monolitic solution of Eqs.(33) is possible using an adequate iterative scheme. However, it is usually more convenient to make use of a fractional step method or a predictor-corrector method. Two interesting approaches of this kind implemented by the authors are described next.

===4.2 Fractional step method===

A fractional step scheme is derived by noting that the discretized momentum equation (33a) can be split into the two following equations

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>{\boldsymbol M} \displaystyle{1\over \Delta t} (\tilde{\boldsymbol u}^{n+1}-\bar{\boldsymbol u}^{n}) + {\boldsymbol H}^{n+\theta }\bar{\boldsymbol u} ^{n+\theta } - \alpha {\boldsymbol G}\bar{\boldsymbol p}^{n} + {\boldsymbol C}^{n+\theta }\bar {\boldsymbol c}^{n+\theta }={\boldsymbol f}^{n+\theta }</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" |  (34a)
|}

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>{\boldsymbol M} \displaystyle{1\over \Delta t} (\bar{\boldsymbol u}^{n+1}-\tilde{\boldsymbol u}^{n+1})- {\boldsymbol G}(\bar{\boldsymbol p}^{n+1}-\alpha \bar {\boldsymbol p}^{n})={\boldsymbol 0}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" |  (34b)
|}

In Eqs.(34) <math display="inline">\tilde{\boldsymbol u}^{n+1}</math> is a predicted value of the velocity at time <math display="inline">n+1</math> and <math display="inline">\alpha </math> is a variable whose values of interest are zero and one. For <math display="inline">\alpha =0</math> (first order scheme) the splitting error is of order <math display="inline">0 (\Delta t)</math>, whereas for <math display="inline">\alpha =1</math> (second order scheme) the error is of order <math display="inline">0 (\Delta t^2)</math> <span id='citeF-52'></span>[[#cite-52|[52]]].

Eqs.(34) are completed with the following three equations emanating from Eqs.(33b-d)

{| 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}^{n+1}+\hat {\boldsymbol L}^n \bar{\boldsymbol p}^{n+1}+{\boldsymbol Q}\bar {\boldsymbol \pi }^{n}= {\boldsymbol 0}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" |  (35a)
|}

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

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>{\boldsymbol Q}^T \bar{\boldsymbol p}^{n+1}+ \hat {\boldsymbol M}^{n+1} \bar {\boldsymbol \pi }^{n+1}= {\boldsymbol 0}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" |  (35c)
|}

The value of <math display="inline">\bar{\boldsymbol u}^{n+1}</math> obtained from Eq.(34b) is substituted into Eq.(35a) to give

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>{\boldsymbol G}^T\tilde{\boldsymbol u}^{n+1}+ \Delta t {\boldsymbol G}^T {\boldsymbol M}^{-1} {\boldsymbol G} (\bar{\boldsymbol p}^{n+1} - \alpha \bar{\boldsymbol p}^n)+ \hat {\boldsymbol L}^n {\boldsymbol p}^{n+1}+{\boldsymbol Q}\bar {\boldsymbol \pi }^n={\boldsymbol 0} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (36)
|}

The product <math display="inline">{\boldsymbol G}^T {\boldsymbol M}^{-1} {\boldsymbol G}</math> can be approximated by a laplacian matrix, i.e.

{| 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}={1\over \rho } {\boldsymbol L}\quad \hbox{with }L^{ab}=\int _{ \Omega ^e} {\boldsymbol \nabla }^T N^a {\boldsymbol \nabla } N^b d\Omega  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (37)
|}

where <math display="inline">L^{ab}</math> are the element contributions to <math display="inline">{\boldsymbol L}</math>.

The steps of the fractional step scheme are:

'''Step 1'''

Eq.(34a) is linearized in the following form

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>{\boldsymbol M} {\tilde{\boldsymbol u}^{n+1} - \bar{\boldsymbol u}^{n}\over \Delta t} + \tilde{\boldsymbol H}^{n+\theta }\tilde{\boldsymbol u}^{n+\theta } - \alpha {\boldsymbol G} \bar{\boldsymbol p}^n + \tilde{\boldsymbol C}^{n+\theta }\bar{\boldsymbol c}^{n}= \bar{\boldsymbol f}^{n+\theta }</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (38)
|}

where <math display="inline">\tilde{\boldsymbol u}^{n+\theta }=\theta \tilde{\boldsymbol u}^{n+1} + (1-\theta )\bar{\boldsymbol u}^{n}</math>, <math display="inline"> \tilde{\boldsymbol H}^{n+\theta } = {\boldsymbol H} (\tilde{\boldsymbol u}^{n+\theta })</math>, and <math display="inline">\tilde{\boldsymbol C}^{n+\theta }= {\boldsymbol C} (\tilde{\boldsymbol u}^{n+\theta })</math>. We have chosen in our work <math display="inline">\theta =0</math>. For this value, the fractional nodal velocities <math display="inline">\tilde {\boldsymbol u}^{n+1}</math> can be  explicitely computed from Eq.(38) by

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\tilde{\boldsymbol u}^{n+1} = \bar{\boldsymbol u}^{n} - \Delta t  {\boldsymbol M}^{-1}_d [\tilde{\boldsymbol H}^{n}\bar{\boldsymbol u}^{n} - \alpha {\boldsymbol G} \bar{\boldsymbol p}^n + {\boldsymbol C}^{n}\bar{\boldsymbol c}^{n} - \bar{\boldsymbol f}^{n}]</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (39)
|}

where <math display="inline">{\boldsymbol M}_d</math> is the lumped diagonal form of '''M'''.

<br/>

'''Step 2''' Compute <math display="inline">\bar{\boldsymbol p}^{n+1}</math> from Eq.(36) as

{| 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}= -[\hat{\boldsymbol L}^n + {\Delta t \over \rho } {\boldsymbol L}]^{-1} [{\boldsymbol G}^T\tilde{\boldsymbol u}^{n+1} - \alpha {\Delta t \over \rho }{\boldsymbol L}\bar{\boldsymbol p}^n +{\boldsymbol Q} \bar{\boldsymbol \pi }^{n}] </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (40)
|}

'''Step 3''' Compute <math display="inline"> \bar{\boldsymbol u}^{n+1}</math> explicitly from Eq.(34a) as

{| 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}=\tilde{\boldsymbol u}^{n+1}+ \Delta t  {\boldsymbol M}_d^{-1} {\boldsymbol G} (\bar {\boldsymbol p}^{n+1}- \alpha \bar {\boldsymbol p}^n) </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (41)
|}

<br/>

'''Step 4''' Compute <math display="inline"> \bar{\boldsymbol c}^{n+1}</math> explicitly from Eq.(35b) as

{| 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;" | (42)
|}

'''Step 5''' Compute <math display="inline">  \bar{\boldsymbol \pi }^{n+1}</math> explicitly from Eq.(35c) as

{| 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;" | (43)
|}

Above algorithm has improved stabilization properties versus the standard segregation methods due to the introduction of the laplacian matrix <math display="inline">\hat{\boldsymbol L}</math> in Eq.(40). We note that this matrix emanates from the FIC stabilization terms.

The boundary conditions are applied as follows. No condition is applied in the computation of the fractional velocities <math display="inline">\tilde{\boldsymbol u}^{n+1}</math> in Eq.(39). The prescribed velocities at the boundary are applied when solving for <math display="inline">\bar{\boldsymbol u}^{n+1}</math> in the step 3. The prescribed pressures at the boundary are imposed by making <math display="inline">\bar{\boldsymbol p}^n</math> equal to the prescribed pressure values.

===4.3 Three steps fractional scheme===

Steps 4 and 5 can be elliminated by substituting the expression of <math display="inline">\bar{\boldsymbol c}^{n+1}</math> and <math display="inline">\bar{\boldsymbol \pi }^{n+1}</math> from Eqs.(42) and (41) into (39) and (40), respectively, where <math display="inline">\bar{\boldsymbol c}</math> and <math display="inline">\bar{\boldsymbol \pi }</math> are now sampled at <math display="inline">t=n+1</math>. The resulting three steps scheme has few advantages versus the five steps scheme described above, as the solution for <math display="inline">\bar{\boldsymbol u}^{n+1}</math> in Eq.(39) can not longer be made explicit and it requires the inversion of a non symmetric matrix.

===4.4 Predictor-corrector scheme===

The fractional step method (of Section 4.2) can be taken as the basis for deriving a predictor-multicorrector scheme which converges to the monolithic time discretized problem. Denoting by <math display="inline">i</math> the <math display="inline">i</math>th iteration of the scheme the resulting linearized system is

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>{\boldsymbol M} \displaystyle{1\over \theta \Delta t} (\bar{\boldsymbol u}^{n+\theta ,i+1}-\bar{\boldsymbol u}^{n}) + [{\boldsymbol H}^{n+\theta ,i}]\bar{\boldsymbol u} ^{n+\theta ,i+1} - {\boldsymbol G}\bar{\boldsymbol p}^{n+1,i}+ \left[{\boldsymbol C}^{n+\theta ,i}\right]{\boldsymbol c}^{n+\theta ,i}={\boldsymbol f}^{n+1}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (44a)
|}

{| 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}^{n+\theta ,i+1} + \left[\hat{\boldsymbol L}^{n+1,i+1} +\displaystyle {\Delta t\over \rho }{\boldsymbol L}\right]\bar {\boldsymbol p}^{n+1,i+1}-  \displaystyle{\Delta t\over \rho } {\boldsymbol L} \bar {\boldsymbol p}^{n+1,i}+ {\boldsymbol Q}\bar {\boldsymbol \pi }^{n+1,i}={\boldsymbol 0}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (44b)
|}

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\left[\hat{\boldsymbol C}^{n+\theta ,i+1}\right]\bar{\boldsymbol u}^{n+\theta ,i+1}+{\boldsymbol M}\bar {\boldsymbol c}^{n+\theta ,i+1}={\boldsymbol 0}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (44c)
|}

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>{\boldsymbol Q}^T \bar {\boldsymbol p}^{n+1,i+1}+[\hat {\boldsymbol M}^{n+1,i+1}]\bar {\boldsymbol \pi }^{n+1,i+1}={\boldsymbol 0}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (44d)
|}

In our work we have chosen <math display="inline">\theta =1</math> (backward differencing). Indeed other schemes are possible (i.e. <math display="inline">\theta = 1/2</math>, Crank-Nicolson, etc) <span id='citeF-52'></span>[[#cite-52|[52]]].

Eqs.(44a-d) are solved iteratively in a staggered manner for the values of <math display="inline">\bar{\boldsymbol u}^{n+1}</math>, <math display="inline">\bar{\boldsymbol p}^{n+1}</math>, <math display="inline">\bar{\boldsymbol c}^{n+1}</math> and <math display="inline">\bar {\boldsymbol \pi }^{n+1}</math>, respectively.

The difference with a standard iterative scheme for the monolithic problem comes from the terms involving the laplacian matrix <math display="inline">{\boldsymbol L}</math> in Eq.(44b). These terms emanate from Eq.(40) in the fractional step scheme by making <math display="inline">\alpha =1</math> and <math display="inline">\bar{\boldsymbol p}^{n}\equiv \bar{\boldsymbol p}^{n+1,i}</math>. This idea was originally proposed by Codina <span id='citeF-52'></span>[[#cite-52|[52]]] and it is here extended in the context of the FIC formulation.

===4.5 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.(18)):

{| 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;" | (45)
|}

The resulting discretized system of equations can be written as (see Eqs.(30))

{| 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} + \hat{\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;" | (46)
|}

The algorithms of previous section can now be implemented. We note that convergence of the predictor-corrector scheme is now faster due to the absence of the non linear convective terms in the momentum equation.

The steady-state form of Eqs.(46) can be expressed in matrix form as

{| 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 & - \hat{\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;" | (47)
|}

The system is symmetric and always positive definite and therefore leads to a non singular solution. 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

{| 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 & - (\hat{\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;" | (48)
|}

The reduction process is simplified by using a diagonal form of matrix <math display="inline">\hat {\boldsymbol M}</math>. Applications of this scheme to incompressible solid mechanics problems are reported in <span id='citeF-60'></span>[[#cite-60|[60]]].

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

===5.1 General coupled solution scheme===

The  algorithms of Section 4 can be readily extended for fluid-structure interaction analysis. The solution process in all  cases involves 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;fluidflow&nbsp;forces===

This implies solving the dynamic equations of motion for the structure written as

{| 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;" | (49)
|}

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.(49) 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.(49) 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.Computethenew&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  next.

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

Tezduyar ''et al.'' <span id='citeF-54'></span>[[#cite-54|[54]]] and Chiandussi ''et al.'' <span id='citeF-55'></span>[[#cite-55|[55]]] have proposed simple method for  moving the  mesh nodes based on the iterative solution of a fictitious linear elastic problem on the mesh domain. In the method introduced in <span id='citeF-54'></span>[[#cite-54|[54]]], the mesh deformation is handled selectively based on the element sizes and deformation modes, with the objective to increase stiffening of the smaller elements, which are typically located near solid surfaces. In Chiandusi ''et al.'' <span id='citeF-55'></span>[[#cite-55|[55]]] 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. A simple and effective procedure is to select the Poisson's ratio <math display="inline">\nu =0</math> and compute the “equivalent” Young modulus for each element by

{| 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;" | (50)
|}

where <math display="inline">\varepsilon _i</math> are the principal strains, <math display="inline">\bar E</math> is an arbitrary value of the Young modulus  and <math display="inline">\bar \varepsilon </math> is a prescribed uniform strain field. <math display="inline">\bar E</math> and <math display="inline">\bar  \varepsilon </math> are constant for all the elements in the  mesh.

In our work we have adopted the solution scheme proposed by Chiandusi ''et al.'' <span id='citeF-55'></span>[[#cite-55|[55]]]. This basically involves the following two steps.

''Step 1''. Consider the FE mesh as a linear elastic solid with homogeneous material properties characterized by a prescribed uniform strain field <math display="inline">\bar  \varepsilon </math>, an arbitrary Young modulus <math display="inline">\bar E</math> and <math display="inline">\nu =0</math>. Solve a linear elastic problem with imposed displacements at the mesh boundary defined by the actual movement of the boundary nodes. An approximate solution to this linear elastic problem, such as that given by the first iterations of a conjugate gradient solution scheme, suffices for practical purposes.

''Step 2''. Compute the principal strains  in each element. Repeat the (approximate) FE solution of the linear elastic problem with prescribed boundary displacements using the values of <math display="inline">E</math> of Eq.(50).

The  algorithm is able to treat the movement of the mesh due to changes in position of fully submerged and semi-submerged bodies such as ships. However if the floating body intersects the free-surface, the changes in the analysis domain  can be very important as emersion or immersion of significant parts of the body can occur within a time step. A possible solution to this problem is to remesh the analysis domain. However, for most  problems a mapping of the moving surfaces linked to a mesh updating algorithm described above can avoid remeshing <span id='citeF-38'></span>[[#cite-38|[38]]].

===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;" | (51)
|}

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 approximates 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 <span id='citeF-41'></span>[[#cite-41|[41]]].

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-38'></span>[[#cite-38|[38]]]

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

where

{| 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 ,\quad i=1,2 </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (53)
|}

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.(51). The underlined term in Eq.(52) introduces the necessary stabilization for the solution of the highly convective (and non linear) equation defining the evolution of the wave elevation.

The solution in time of Eq.(52) can be expressed in terms of the nodal velocities computed from the flow solution, as

{| 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} {\partial \beta ^n \over \partial x_i}-v_3^{n+1}-{h_{\beta _i}\over 2} {\partial r_\beta ^n \over \partial x_i}\right]\quad ,\quad  i=1,2 </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (54)
|}

Eq.(54) 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-38'></span>[[#cite-38|[38]]]. This solution step  preceeds the computation of the structure motion in the case of a fluid-structure interaction problem. Typically the general algorithm is as follows:

<ol>

<li>Solve for the nodal velocities <math display="inline">\bar{\boldsymbol u}^{n+1}</math> and the nodal 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.(54). </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.(49)). </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

{| 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;" | (55)
|}

where <math display="inline">g</math> is the  gravity constant. Eq.(55) 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 <span id='citeF-41'></span>[[#cite-41|[41]]].

==7 TURBULENCE MODELLING==

The  discussion of the treatment of turbulent effects in the flow equations in the Eulerian and ALE formulations  falls outside the objective of this paper as many of the existing turbulence models are applicable.

In the examples presented in the paper 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 Reynolds 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>.

For the definition of <math display="inline">\mu _T</math> 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, among others <span id='citeF-56'></span>[[#cite-56|[56]]].

==8 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 is also 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.(51) and the convective terms vanish in the momentum equations, while the rest of the fluid flow equations remain unchanged. The resulting governing equations have an identical form as those of the Stokes flow problem, with the motion of the flow particles being referred now to a Lagrangian coordinate frame.

The FEM algorithms for solving the Lagrangian flow equations are very similar to those for the Eulerian or ALE descriptions presented earlier. For preciseness we present a particular class of Lagrangian formulation  to solve problems involving the interaction between free surface flows  and solids in a unified manner. The approach, called the ''particle finite element method'' (PFEM) treats the mesh nodes in the fluid and solid domains as  dimensionless particles which can freely move and even separate from the main fluid domain representing, for instance, the effect of water drops. A finite element mesh connects the nodes defining the discretized domain where the governing equations are solved in the standard FEM fashion. The PFEM is the natural evolution of recent work of the authors for  the solution of FSI problems using Lagrangian finite element and meshless methods <span id='citeF-44'></span>[[#cite-44|[44]]].

===8.1 The Particle Finite Element Method (PFEM)===

As mentioned in the previous section in the PFEM approach both the fluid and the solid domains are modelled using a Lagrangian formulation. The finite element method (FEM) is used to solve the continuum equations in both domains. Hence a mesh discretizing these domains must be generated in order to solve the governing equations for both the fluid and solid problems in the standard FEM fashion. We note  that the nodes discretizing the fluid and solid domains can be viewed as material points which motion is tracked during the transient solution.

The Lagrangian formulation allows to track the motion of each  single fluid point (a node). This is useful to model the separation of water particles from the main fluid domain and to follow their subsequent motion as individual dimensionless particles with an initial velocity and subject to gravity forces.

The quality of the numerical solution   depends on the discretization chosen as in the standard FEM. Adaptive mesh refinement techniques can be used to improve the solution in zones where large motions of the fluid or the structure occur.

For each time increment the PFEM involves the following steps.

<ol>

<li>Starting with an initial collection of points (nodes), identify the external boundaries for both the fluid and solid domains. This is  an essential step as some boundaries (such as the free surface in fluids) may be severely distorted during the solution process including separation and re-entring of nodes. The Alpha Shape method <span id='citeF-57'></span>[[#cite-57|[57]]] is used for the boundary definition. </li>
<li>Discretize the fluid and solid domains with a finite element mesh. For the mesh generation process we use and extended Delaunay technique   <span id='citeF-50'></span>[[#cite-50|[50]]].  </li>
<li>Solve the coupled Lagrangian equations of motion for  the fluid and the solid domains. Compute the relevant state variables in both domains at each time step: velocities, pressure and viscous stresses in the fluid and displacements, stresses and strains in the solid.  </li>

The solution scheme chosen in this work is a generalization of the fractional step algorithm of Section 4.2. In summary the solution steps are the following.

''Step 3.1'' Compute the predicted velocities (viz, Eq.(39) for <math display="inline">\alpha =1</math> and <math display="inline">{\boldsymbol C}= {\boldsymbol 0}</math>)

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>

\tilde{\boldsymbol u}^{n+1,i+1} = \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;" | (56)
|}

''Step 3.2'' Compute <math display="inline">\bar {\boldsymbol p}^{n+1,i+1}</math> from Eq.(40) for <math display="inline">\alpha =1</math>.

<br/>

''Step 3.3'' Compute explicitely <math display="inline"> \bar{\boldsymbol u}^{n+1,i+1}</math> from Eq.(41) with <math display="inline">\alpha =1</math>.

<br/>

''Step 3.4'' Compute <math display="inline"> \bar{\boldsymbol \pi }^{n+1,i+1}</math> explicitely from Eq.(43).

<br/>

''Step 3.5'' Solve for the motion of the structure by integrating Eq.(51).

<br/>

''Step 3.6'' Estimate a new position of the mesh nodes in terms of the time increment size as

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>

{\boldsymbol x}_j^{n+1,i+1} = {\boldsymbol x}_j^{n}+\bar {\boldsymbol u}_j^{n+1,i+1} \Delta t </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (57)
|}

where index <math display="inline">j</math> denotes the node number.

It is important to note that all matrices in Steps 3.1&#8211;3.5 are evaluated at the last predicted configuration <math display="inline">\Omega ^{n+1,i}</math>.
<li>In steps 3.1&#8211;3.6 superindex <math display="inline">i</math> denotes the iteration within a time increment.  </li>

<br/>

''Step 3.7'' Check the convergence of the velocity and pressure fields in the fluid and the displacements strains and stresses in the structure. If convergence is achieved frozen the final position of the mesh nodes and move to the next time increment, otherwise return to step 3.1 for the next iteration.

<li>Go back to step 1 and repeat the solution process for the next time increment.  </li>

</ol>

Above algorithm can be found to be analogous to the standard ''updated lagrangian'' scheme typically used in non linear solid mechanics problems <span id='citeF-1'></span>[[#cite-1|[1]]].

Despite the motion of the nodes within the iterative process, in general there is no need to regenerate the mesh at each iteration. ''In the examples presented in the paper the mesh in the fluid domain has been regenerated at each time increment''. A cheaper alternative is to generate a new mesh only after a prescribed number of converged time increments, or when  the nodal displacements induce significant geometrical distorsions in some elements.

The boundary conditions are applied as described in Section 4.2.

In the examples presented in the paper the time increment size has been chosen as

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\Delta t =\min (\Delta t_i ) \quad \hbox{with}\quad \Delta t_i ={\vert {\boldsymbol v}\vert \over h_i^{\min }} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (58)
|}

where <math display="inline">h_i^{\min }</math> is the distance between node <math display="inline">i</math> and the closest node in the mesh.

===8.2 Treatment of contact between fluid and solid interfaces===

The  condition of prescribed velocities or pressures at the solid boundaries in the PFEM are  applied in strong form to the boundary nodes. These nodes might belong to fixed external boundaries or to moving boundaries linked to the interacting solids. In some problems it is useful to define a layer of nodes adjacent to the external boundary in the fluid where the condition of prescribed velocity is imposed. These nodes typically remain fixed during the solution process. Contact between water particles and the solid boundaries is accounted for by the incompressibility condition which ''naturally prevents the water nodes to penetrate into the solid boundaries''. This simple way to treat the water-wall contact is another attractive feature of the PFEM.

===8.3 Generation of a new mesh===

One of the key points for the success of the PFEM  is the fast regeneration of a mesh at every time step on the basis of the position of the nodes in the space domain. As mentioned previously, in our work the mesh is generated using the so called extended Delaunay tesselation (EDT) <span id='citeF-50'></span>[[#cite-50|[50]]]. The EDT allows one to generate non standard meshes combining elements of  arbitrary polyhedrical shapes  (triangles, quadrilaterals and other polygons in 2D and tetrahedra, hexahedra and arbitrary polyhedra in 3D) in a computing time of order <math display="inline">n</math>, where <math display="inline">n</math> is the total number of nodes in the mesh. The <math display="inline">C^\circ </math> continuous shape functions of the elements are obtained using the so called ''meshless finite element interpolation'' (MFEM) <span id='citeF-49'></span>[[#cite-49|[49]]].

===8.4 Identification of boundary surfaces===

One of the main tasks  in the PFEM is the correct definition of the boundary domain. Sometimes, boundary nodes are explicitly marked  differently from internal nodes. In other cases, the total set of nodes is the only information available and the algorithm must recognize the boundary nodes.

Considering that the nodes follow a variable <math display="inline">h(x)</math>  distribution, where <math display="inline">h(x)</math> is the minimum distance between two nodes, the following criterion has been used. ''All nodes on an empty sphere with a radius  greater than <math>\alpha h</math>, are considered as boundary nodes''. In practice, <math display="inline">\alpha </math>  is a parameter close to, but greater than one. This criterion coincides with the Alpha Shape concept <span id='citeF-46'></span>[[#cite-46|[46]]].

In this work, the boundary surface is defined by  all the polyhedral surfaces (or polygons in 2D) having all their nodes on the boundary and belonging to just one polyhedron.

The method  also allows to identify isolated fluid particles outside the main fluid domain. These particles are treated as part of the external boundary where the pressure is fixed to the atmospheric value.

Figure 2 shows a schematic example of the process to identify individual particles (or a group of particles) starting from a given collection of nodes.

<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_881612373-Figure5.png|450px|Identification of individual particles (or a group of particles) starting from a given collection of nodes using the Alpha Shape method.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 2:''' Identification of individual particles (or a group of particles) starting from a given collection of nodes using the Alpha Shape method.
|}

==9 COMPUTATION OF THE CHARACTERISTIC LENGTHS==

The evaluation of the stabilization parameters is one of the crucial issues in stabilized methods. Excellent results have been obtained in all problems solved  using linear  elements with the characteristic length vector  defined by

{| 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;" | (59)
|}

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

{| 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;" | (60)
|-
| 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;" | (61)
|}
|}

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 in the ALE formulation the following value of the characteristic length vector <math display="inline">{\boldsymbol h}_\beta </math> in Eq.(52) has been taken

{| 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;" | (62)
|}

The streamline parameter <math display="inline">\bar h_s</math> has been obtained by Eq.(60) 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

{| 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;" | (63)
|}

The cross-wind terms in eqs.(59) and (63) 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-58'></span>[[#cite-58|[58]]].

==10 EXAMPLES==

The examples chosen show the applicability of the Eulerian, ALE and Lagrangian  (PFEM) formulations presented to solve fluid flow problems.  The first two examples are typical flow problems solved with the standard Eulerian formulation using the predictor-corrector scheme of Section 4.4 with <math display="inline">\theta = 1</math>. The first example is the flow of an inviscid fluid past a cylinder.  Then the standard problem of flow past a NACA airfoil is solved with mesh adaptivity.

The next two examples fall within the category of ship hydrodynamics problems solved with the ALE formulation and the fractional step scheme of Section 4.2 with <math display="inline">\alpha =1</math>. The examples are the analysis of  a scale model of a commercial ship and of an American Cup racing sail boat.  In both cases the free surface equation is solved together with the flow equations as described in Section 6. Numerical results are compared with experimental data.

The last series of examples show applications of the Lagrangian PFEM formulation to the simulation of the collapse of a water column, a sloshing problem, a ship hit by a wave, a solid cube falling into a water recipient and a mixing problem.

===10.1 Inviscid flow past a cylinder===

The problem was solved with the mesh of 10096 linear triangles shown in Figure 3. The initial conditions were an horizontal velocity of value one and zero pressure. A value of  the Reynolds number of infinity was taken. The solution was progressed in time until steady state was found. The “correct” solution in this case is a symmetric velocity and pressure field.

The steady state results in Figure 4 show the relevance of accounting for the convective projection terms in the momentum equations. Results for the velocity and pressure contours shown in Figure 5  were obtained with the formulation described in the paper, whereas the results shown in Figure 6 were obtained by neglecting the convective projection terms (i.e. neglecting the terms involving the <math display="inline">\boldsymbol c</math> variables in Eq.(44a) and skipping the solution of Eq.(44c)). It can be clearly seen that neglecting the convective projection terms leads to a deterioration of the results in the region  behind the cylinder where some oscillations for the pressure and velocity fields are found. These oscillations disappear and the correct symmetric solution is found if the convective projection terms are taken into account as described in the paper.

Evidence of the importance of taking into account the pressure gradient projection terms for the accuracy of the incompressible solution has also been found as reported in <span id='citeF-60'></span>[[#cite-60|[60]]].

<div id='img-3'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[File:Draft_Samper_881612373_3747_Figure3.jpg|450px|Inviscid flow past a cylinder. Mesh of 10.096 three node triangles.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 3:''' Inviscid flow past a cylinder. Mesh of 10.096 three node triangles.
|}

<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_881612373-fig4a.png|600px|]]
|[[Image:Draft_Samper_881612373-fig4b.png|600px|Contours of velocity (upper picture and pressure) obtained accounting for the convective projection terms]]
|- style="text-align: center; font-size: 75%;"
| colspan="2" | '''Figure 4:''' Contours of velocity (upper picture and pressure) obtained accounting for the convective projection terms
|}

<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_881612373-fig5a.png|600px|]]
|[[Image:Draft_Samper_881612373-fig5b.png|600px|Contours of velocity (upper picture and pressure) obtained  not taking into account the convective projection terms]]
|- style="text-align: center; font-size: 75%;"
| colspan="2" | '''Figure 5:''' Contours of velocity (upper picture and pressure) obtained  not taking into account the convective projection terms
|}

<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_881612373-fig1.png|600px|Viscous flow past a NACA airfoil. Characteristic length =1m, Re=100. Original mesh]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 6:''' Viscous flow past a NACA airfoil. Characteristic length =1m, <math>Re=100</math>. Original mesh
|}

===10.2 Flow past a NACA airfoil. Adaptative solution===

The problem is the analysis of a viscous flow past a NACA 12 airfoil. This example sows the performance of the viscous FIC formulation with a mesh refinement scheme. The initial  mesh of 2574 nodes and 4784 linear triangular elements is shown in Figure 6. The problem is solved for a value of the Reynolds number of 100.

The mesh refinement algorithm is based on the equidistribution of the global error over the finite element mesh. The error estimation method  is based on the evaluation of the energy rate (the power) of the FE residuals of the momentum and the incompressibility equations. The residuals are computed using recovered values of the derivatives and pressure variables obtained via a  nodal derivative recovery technique. A nodal-based approach is used for computing the residual power integrals <span id='citeF-64'></span>[[#cite-64|[64]]]. The refinement was performed at a time <math display="inline">t=2</math>s during the transient solution process. The resultant mesh after five refinement steps is shown in  Figures 7 and Figure 8. The corresponding velocity field contours are shown in Figures 9  and 10.

<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_881612373-fig2.png|470px|Mesh obtained after five refinement steps using the equidistribution of the global error.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 7:''' Mesh obtained after five refinement steps using the equidistribution of the global error.
|}

<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_881612373-fig3.png|300px|]]
|[[Image:Draft_Samper_881612373-fig3_b.png|300px|Details of the refined mesh of Figure 7 in the vecinity of the airfoil]]
|- style="text-align: center; font-size: 75%;"
| colspan="2" | '''Figure 8:''' Details of the refined mesh of Figure 7 in the vecinity of the airfoil
|}

<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_881612373-naca_line_1.png|530px|Velocity field after five refinement steps]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 9:''' Velocity field after five refinement steps
|}

<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_881612373-naca_line_2.png|530px|Detail of the velocity field of Figure 9]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 10:''' Detail of the velocity field of Figure 9
|}


Table 1 shows the value of the global error computed as the percentage of the energy rate of the FE residuals (<math display="inline">\Vert P\Vert _\Omega </math>)  versus the total energy rate <math display="inline">U</math>, the errors in the drag and lift values and the number of elements and nodes during the mesh refinement process.

Further examples of  mesh adaptivity in incompressible flows using the FIC formulation can be found in <span id='citeF-64'></span>[[#cite-64|[64]]].


{|  class="floating_tableSCP wikitable" style="text-align: center; margin: 1em auto;min-width:50%;"
|+ style="font-size: 75%;" |Table. 1 Flow past a NACA airfoil. Convergence of the global error  with the mesh refinement
|- style="border-top: 2px solid;"
| style="border-left: 2px solid;border-right: 2px solid;" |  Mesh 
| style="border-left: 2px solid;border-right: 2px solid;" | Global error <math display="inline">\left({\Vert P\Vert _\Omega \over U}\right)</math>(%)  
| style="border-left: 2px solid;border-right: 2px solid;" | Nodes 
| style="border-left: 2px solid;border-right: 2px solid;" | Elements
| style="border-left: 2px solid;border-right: 2px solid;" | Drag error %
| style="border-left: 2px solid;border-right: 2px solid;" | Lift error% 
|- style="border-top: 2px solid;"
| style="border-left: 2px solid;border-right: 2px solid;" |  1
| style="border-left: 2px solid;border-right: 2px solid;" | 5.55
| style="border-left: 2px solid;border-right: 2px solid;" | 2574
| style="border-left: 2px solid;border-right: 2px solid;" | 4784
| style="border-left: 2px solid;border-right: 2px solid;" | 29.1810942 
| style="border-left: 2px solid;border-right: 2px solid;" | 2.99902143
|- style="border-top: 2px solid;"
| style="border-left: 2px solid;border-right: 2px solid;" |  2
| style="border-left: 2px solid;border-right: 2px solid;" | 2.93
| style="border-left: 2px solid;border-right: 2px solid;" | 5340
| style="border-left: 2px solid;border-right: 2px solid;" | 10680
| style="border-left: 2px solid;border-right: 2px solid;" | 28.7443059
| style="border-left: 2px solid;border-right: 2px solid;" | 2.5500113
|- style="border-top: 2px solid;"
| style="border-left: 2px solid;border-right: 2px solid;" |  3
| style="border-left: 2px solid;border-right: 2px solid;" | 2.87
| style="border-left: 2px solid;border-right: 2px solid;" | 7942
| style="border-left: 2px solid;border-right: 2px solid;" | 15884
| style="border-left: 2px solid;border-right: 2px solid;" | 22.1078024
| style="border-left: 2px solid;border-right: 2px solid;" | 2.00851784
|- style="border-top: 2px solid;"
| style="border-left: 2px solid;border-right: 2px solid;" |  4
| style="border-left: 2px solid;border-right: 2px solid;" | 2.77
| style="border-left: 2px solid;border-right: 2px solid;" | 11948
| style="border-left: 2px solid;border-right: 2px solid;" | 23896
| style="border-left: 2px solid;border-right: 2px solid;" | 15.2893198
| style="border-left: 2px solid;border-right: 2px solid;" | 1.54898957
|- style="border-top: 2px solid;"
| style="border-left: 2px solid;border-right: 2px solid;" |  5
| style="border-left: 2px solid;border-right: 2px solid;" | 1.81
| style="border-left: 2px solid;border-right: 2px solid;" | 17142
| style="border-left: 2px solid;border-right: 2px solid;" | 34284
| style="border-left: 2px solid;border-right: 2px solid;" | 10.6478917
| style="border-left: 2px solid;border-right: 2px solid;" | 1.17033653
|- style="border-top: 2px solid;"
| style="border-left: 2px solid;border-right: 2px solid;" |  6
| style="border-left: 2px solid;border-right: 2px solid;" | 1.81
| style="border-left: 2px solid;border-right: 2px solid;" | 17255
| style="border-left: 2px solid;border-right: 2px solid;" | 34510
| style="border-left: 2px solid;border-right: 2px solid;" | 8.54878424
| style="border-left: 2px solid;border-right: 2px solid;" | 1.02567233
|- style="border-top: 2px solid;border-bottom: 2px solid;"
| style="border-left: 2px solid;border-right: 2px solid;" |  7
| style="border-left: 2px solid;border-right: 2px solid;" | 1.71
| style="border-left: 2px solid;border-right: 2px solid;" | 17372
| style="border-left: 2px solid;border-right: 2px solid;" | 34744
| style="border-left: 2px solid;border-right: 2px solid;" | 6.93290243
| style="border-left: 2px solid;border-right: 2px solid;" | 0.86514698

|}

===10.3 KVLCC2 hull model===

The example is the analysis of the so called KVLCC2 benchmark model proposed by the Korean Research Institute of Ship and Ocean Engineering (KRISO) <span id='citeF-64'></span>[[#cite-64|[64]]]. Here a partially wetted tramsom stern is expected due to the low Froude number of the test. Figure&nbsp;11 shows the NURBS geometry of the ship provided by KRISO. The obtained results are compared with the experimental data available  <span id='citeF-65'></span>[[#cite-65|[65]]].

<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_881612373-KCSgeo.png|600px|KVLCC2 model. Geometrical definition based on NURBS surfaces.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 11:''' 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 12. A total of 550000 tetrahedra were used to model the virtual towing basin.  The two following test cases were analyzed.

<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_881612373-KCSmesh.png|600px|KVLCC2 model. Surface mesh used in the analysis.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 12:''' KVLCC2 model. Surface mesh used in the analysis.
|}


''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/sec, Froude Number: 0.142.
* Viscosity: <math display="inline">0.00126 Kg/m\cdot sec</math>, Density: <math display="inline">1000 Kg/m^3</math>, Reynolds number: <math display="inline">4.63\times 10^6</math>.

The turbulence model chosen in this case was the <math display="inline">K</math> model. Figures 13 and 14 show the wave profiles on the hull and in  a cut at y/L = 0.082 compared to the experimental data.  The numerical results are quantitatively good close to the hull. A loss 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-13'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Samper_881612373-KCSprof1.png|600px|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 13:''' KVLCC2 model. Wave profile on the hull compared to experimental data. Thick line shows numerical results
|}

<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_881612373-KCSprof2.png|600px|KVLCC2 model. Wave profile on a cut at y/L=0.0964 compared to experimental data <span id='citeF-64'></span>[[#cite-64|64]]. Thick line shows numerical results]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 14''' 
|}


''Test 2.-'' ''Wake analysis at different planes''. Several turbulence models were used  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:

* 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 15&#8211;16 present results corresponding to the test 2.  Figure 15  shows the contours of the axial (X) component of the velocity on a plane at 2.71 m  from the orthogonal aft. Figure 16 shows the maps of the kinetic energy on this plane. Experimental results are shown for comparison in all cases. Further results for this problem and other similar ones can be found in <span id='citeF-38'></span>[[#cite-38|[38]]].

<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_881612373-KCSvelx1.png|600px|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 15:''' 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-16'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Samper_881612373-KCSk.png|600px|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 16:''' 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.4 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 17.  The numerical scheme is the same as that used for the previous example. The results presented in Figures&nbsp;17&#8211;20 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-67'></span>[[#cite-67|[67]]].

<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_881612373-ACmesh.png|600px|Bravo ~Espãna sail racing boat. Mesh used in the analysis.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 17:''' <math>Bravo ~Espa\tilde{n}a</math> sail racing boat. Mesh used in the analysis.
|}

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

<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_881612373-ACres2.png|400px|Bravo ~Espãna. Streamlines.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 19:''' <math>Bravo ~Espa\tilde{n}a</math>. Streamlines.
|}

<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_881612373-ACgraph.png|700px|Bravo~Espãna. Resistance test. Comparison of numerical results with experimental data.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 20:''' <math>Bravo~Espa\tilde{n}a</math>. Resistance test. Comparison of numerical results with experimental data.
|}

===10.5 Collapse of a water column===

The first problem solved to show the potential of the Lagrangian PFEM is the study of the collapse of a water column. This problem was solved by Koshizuka and Oka <span id='citeF-69'></span>[[#cite-69|[69]]] both experimentally and numerically. It has became a classical example to  validate the Lagrangian formulation for fluid flows. The water is initially kept within a rectangular container including a removable vertical board. A double layer of nodes in the solid walls is used in order to prevent water nodes from exiting the analysis domain. The boundary conditions impose zero velocity at the wall nodes and zero (atmospheric) pressure at the free surface. The method allows one to follow the large motion of the water particles including separation of some water drops. The collapse starts at time t = 0, when the board is removed. Viscosity and surface tension are neglected in the analysis. Figure&nbsp;21 shows the point positions at different time steps. The dark points represent the free surface detected with the algorithm described in Section 8. The internal points are shown in a gray colour and the fixed points in black.

<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_881612373-Fig1ArtBathe.png|450px|Water column collapse at different time steps.]]
|-
|[[Image:Draft_Samper_881612373-Fig1contBathe.png|450px|cont.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 21:''' Water column collapse at different time steps.
|}


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

Figure 22 shows the finite element mesh generated at a time step. We recall that this mesh is used to solve the equations of motion of the fluid  particles as described in the previous sections.

The 3D solution of the same case is shown in Figure 23. More information on the PFEM solution of this problem can be found in <span id='citeF-45'></span>[[#cite-45|[45]]].

<div id='img-23'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Samper_881612373-Figure2c.png|420px|Finite element mesh discretizing the fluid domain and the container walls at a certain time step.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 22:''' Finite element mesh discretizing the fluid domain and the container walls at a certain time step.
|}

<div id='img-24'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Samper_881612373-Fig7_1.png|600px|]]
|-
|(a) t = 0 sec. (b) t = 0.2 sec.
|-
|[[Image:Draft_Samper_881612373-Fig7_2.png|600px|]]
|-
|(c) t = 0.4 sec. (d) t = 0.6 sec.
|-
| [[Image:Draft_Samper_881612373-Fig7_3.png|600px|]]
|-
|(e) t = 0.8 sec. (f) t = 1.1 sec.
|- style="text-align: center; font-size: 75%;"
| colspan="2" | '''Figure 23:''' Water column collapse in a 3D domain.
|}

===10.6 Sloshing problems===

The simple problem of the free oscillation of an incompressible liquid in a container is considered next. A numerical and analytical solution for this problem can be found in <span id='citeF-70'></span>[[#cite-70|[70]]]. Figure 24 shows a schematic view of the problem and the point distribution in the initial position. The dark points represent the fixed points on the walls where the velocity is fixed to zero.

<div id='img-24'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Samper_881612373-Fig3.png|300px|]]
|[[Image:Draft_Samper_881612373-Fig4.png|500px|]]
|- style="text-align: center; font-size: 75%;"
| colspan="2" | '''Figure 24:''' Sloshing. Initial point distribution and comparison of the numerical and analytical solutions.
|}

Figure 24 also shows the  time evolution of the amplitude compared with the analytical results for the near inviscid case.  Little numerical viscosity is observed on the phase wave and amplitude in spite of the relative coarse distribution of nodes.

The analytical solution is only acceptable for small wave amplitudes. For larger amplitudes, additional waves are overlapping and, finally, the wave breaks and also some particles separate from the fluid domain due to their large velocity.  Figure 25 shows the numerical results obtained with the PFEM for this case. Breaking waves as well as separation effects can be seen on the free surface. This particular and very complicated effect is  well represented by the PFEM.

In order to further test the potential  of the PFEM the same sloshing problem was solved in 3D.  Figure 26 shows the different point positions at two time steps. Each point position was represented by a sphere and only a half of the fixed recipient is represented on the figure. This  representation is only used in order to improve the visualization of the numerical results.

{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Samper_881612373-Fig5.png|600px|]]
|- style="text-align: center; font-size: 75%;"
| '''Figure 25:''' PFEM results for a large amplitude sloshing problem
|}


{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Samper_881612373-Fig6.png|600px|]]
|- style="text-align: center; font-size: 75%;"
| '''Figure 26:''' 3D sloshing problem 
|}

===10.7 Container ship hit by an incoming wave===

Figure 27 shows the analysis of the motion of the transverse sections of a container ship hit by an incoming wave. The dynamic motion of the ship is induced by the resultant of the pressure and the viscous forces acting on the ship boundaries.  The section of the ship analyzed corresponds  to that of a real container ship. The  ship is  assumed to be rigid and is free to move laterally due to the sea wave forces. The objective of the study was to asses the influence of the stabilizers in the ship roll. The figures show clearly how the PFEM predicts the ship and wave motions in a realistic manner. Further examples of this type solved with the PFEM can be found in <span id='citeF-46'></span>[[#cite-46|[46]]].


<div id='img-27'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Samper_881612373-Figure33.png|600px|Ship with stabilizers hit by a lateral wave ]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 27:''' Ship with stabilizers hit by a lateral wave 
|}

===10.8 Rigid square falling in a recipient with water===

In the next  example a  solid square is initially  free and falls down within a water recipient. The square is modelled a  rigid solid subjected to  pressure and viscous forces acting in its boundaries. The resultant of the fluid forces and the weight of the square are applied to the center of the square. These forces govern the displacement of the square which is computed by solving  the dynamic equations of motion as described in the fractional step algorithm of Section 6, similarly as for the rigid ship of the  previous example. Here again  the moving square contours define a boundary condition for the fluid particles at each time step.

Initially the solid falls down freely due to the gravity forces (Figure 28). Once in contact with the water surface, the Alpha-Shape method recognizes the different boundary contours which are shown with a thick line in the figure.  The pressure  and viscous forces are evaluated in all the domain and in particular on the square contours. The fluid forces introduce a negative acceleration in the vertical motion until, once the square is completely inside the water, the vertical  velocity  becomes zero. Then,  buoyancy  forces bring the square up to the free-surface. It is interesting to observe that there is a rotation of the square. The reason is that the center of the floating forces is higher in the rotated position than in the initial ones.

Figure 29 shows a repetition of the same problem showing now all the finite elements in the mesh discretizing the fluid. We recall that in all the  problems here described the mesh in the fluid domain ''is regenerated at each time step'' combining linear triangles and quadrilateral elements as described in Section 8.3. Note that some fluid particles separate from the fluid domain. These particles are treated as free boundary points with zero pressure and hence fall down due to gravity.

It is interesting to see that the final position of the square is different from that of Figure 28. This is due to the unstable character of the square motion. A small difference in the numerical computations (for instance in the mesh generation process) shifts the movement of the square towards the right or the left. Note  that a final rotated equilibrium position is found in both cases. For further details see <span id='citeF-46'></span>[[#cite-46|[46]]].


<div id='img-26'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Samper_881612373-Fig13.png|600px|Square falling into a recipient with water. The square is modelled as a rigid solid. Motion of the square and free surface positions at different time steps.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 28:''' Square falling into a recipient with water. The square is modelled as a rigid solid. Motion of the square and free surface positions at different time steps.
|}


<div id='img-27'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Samper_881612373-Figure10_con231.png|450px|]]
|-
|[[Image:Draft_Samper_881612373-Figure10cont_con231.png|450px|]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 29:''' Square falling in a water recipient. The square is modeled as a  rigid solid. The finite element meshes generated at the selected instants are shown.
|}

===10.9 Mixing of particles within a fluid===

Figure 30 shows an example of application of the PFEM to the mixing of a collection of particles within a container containing a fluid of a higher density. Initially the particles are thrown into the container and mix with the fluid as shown. As time evolves the particles move up naturally towards the surface of the fluid due to their smaller density. This  example clearly shows the possibilities of the PFEM for analysis of fluid mixing situations.


<div id='img-30'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Samper_881612373-Fig_mixing.png|458px|Mixing of particles in a fluid. Evolution of the particles during the mixing process]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 30:''' Mixing of particles in a fluid. Evolution of the particles during the mixing process
|}

==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. Fractional step and predictor-corrector algorithms with intrinsic stabilization properties can be readily derived from the FIC governing equations. Free surface wave effects and fluid-structure interaction situations  can be accounted for in a straightforward manner within the general flow solution algorithm. The ALE formulation is particularly adequate for analysis of problems involving free surface waves of moderate amplitude, typical of ship hydrodynamics situations. The Lagrangian PFEM formulation is very effective for fluid flow problems involving large motions of the free surface, splashing of waves,  complex fluid-structure interactions and mixing problems.

==ACKNOWLEDGEMENTS==

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

Examples 10.1&#8211;10.4 were analyzed with the finite element code Tdyn based on the FEM formulation here presented [70].

Thanks are also given to Profs. R.L. Taylor and R. Löhner, Dr. R. Flores and Mr. R. Aubry for many useful discussions.

==BIBLIOGRAPHY==

<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>
'''[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;99, 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>
'''[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>
'''[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”,  ''Int. J. Num. Meth. Engng.'', '''60'''(1), 255&#8211;281, 2004.

<div id="cite-40"></div>
'''[40]''' J. García and E. Oñate, “An unstructured finite element solver for ship hydrodynamic problems”, ''Journal of Applied Mechanics'', '''70''', 18&#8211;26, January 2003.

<div id="cite-41"></div>
'''[[#citeF-41|[41]]]'''  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-42"></div>
'''[[#citeF-42|[42]]]''' 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, 1998.

<div id="cite-43"></div>
'''[43]'''  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-44"></div>
'''[[#citeF-44|[44]]]'''  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”, ''Fifth World Congress on Computational Mechanics'', Mang HA, Rammerstorfer FG and Eberhardsteiner J. (eds.), July 7&#8211;12, Viena, Austria, 2002.

<div id="cite-45"></div>
'''[[#citeF-45|[45]]]''' S.R. Idelsohn, E. Oñate and F. Del Pin, “A lagrangian meshless finite element method applied to fluid-structure interaction problems”, ''Computer and Structures'', '''81''', 655&#8211;671, 2003.

<div id="cite-46"></div>
'''[[#citeF-46|[46]]]'''  E. Oñate, S.R. Idelsohn, F. Del Pin and R. Aubry,  “The particle finite element method. An overview”, Submitted to ''Int. J. Comput. Meth.'', '''1''' ('''2'''), 2004.

<div id="cite-47"></div>
'''[47]'''  R. Aubry, S.R. Idelsohn  and E. Oñate, “Particle finite element method in fluid mechanics including thermal convection-diffusion”. ''Computer & Structures'', submitted, 2004.

<div id="cite-48"></div>
'''[48]'''  S.R. Idelsohn, E. Oñate and F. Del Pin, “The Particle Finite Element Method:  a powerful tool to solve incompressible flows with free-surfaces and breaking waves”, ''Int. J. Num. Meth. Engng.'', (Submitted), 2004.

<div id="cite-49"></div>
'''[[#citeF-49|[49]]]''' S.R. Idelsohn, E. Oñate, N. Calvo and F. del Pin, “The meshless finite element method”,  ''Int. J. Num. Meth. Engng.'', Vol. '''58''', pp. 893- 912, 2003.

<div id="cite-50"></div>
'''[[#citeF-50|[50]]]'''  N. Calvo, S.R. Idelsohn and E. Oñate, “The extended Delaunay tesselation”, ''Engineering Computations'', '''20''' (5/6), 583&#8211;600, 2003.

<div id="cite-51"></div>
'''[51]''' S.R. Idelsohn, N. Calvo and E. Oñate, “Polyhedrization of an arbitrary 3D point set”, ''Comput. Meth. Appl. Mech. Engng.'', '''192''', 2649&#8211;2667, 2003.

<div id="cite-52"></div>
'''[[#citeF-52|[52]]]'''  R. Codina, “Pressure stability in fractional step finite element method for incompressible flows”, ''J. Comput. Physics'', '''170''', 112&#8211;140, 2001.

<div id="cite-53"></div>
'''[53]'''  R. Codina and O. Soto, “Approximation of the incompressible Navier-Stokes equations using orthogonal-subscale stabilization and pressure segregation on anisotropic finite element meshes”, ''Comput. Meth. Appl. Mech. Eng.'', to appear.

<div id="cite-54"></div>
'''[[#citeF-54|[54]]]'''  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-55"></div>
'''[[#citeF-55|[55]]]''' 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-56"></div>
'''[[#citeF-56|[56]]]''' D.C. Wilcox, ''Turbulence modeling for CFD'', DCW Industries Inc., 1994.

<div id="cite-57"></div>
'''[[#citeF-57|[57]]]'''  H. Edelsbrunner  and E.P. Mucke, E.P., “Three dimensional alpha shapes”, ''ACM Trans. Graphics'', '''13''', 43&#8211;72, 1999.

<div id="cite-58"></div>
'''[[#citeF-58|[58]]]'''  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-59"></div>
'''[59]'''  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-60"></div>
'''[[#citeF-60|[60]]]'''   E. Oñate, J. Rojek, R.L. Taylor and O.C. Zienkiewicz, “Finite calculus formulation for incompressible solids using linear triangles and tetrahedra”, ''Int. J. Num. Meth. Engng.'',  '''59''', 1473&#8211;1500, 2004.

<div id="cite-61"></div>
'''[61]'''   E. Oñate, R.L. Taylor, O.C. Zienkiewicz and J. Rojek, “A residual correction method based on finite calculus”, ''Engineering Computations'', '''20''', (5/6), 629&#8211;658, 2003.

<div id="cite-62"></div>
'''[62]'''  J.C. Cante, J. Oliver and C. González, “Powder forming simulation with the particle finite element method”, ''Computational Mechanics'', Z.H. Yao, M.W. Yuan and W.X. Zhong (Eds.), CD-Rom Proceedings of the WCCM VI, Beijing, September 5&#8211;10, 2004.

<div id="cite-63"></div>
'''[63]'''  E. Oñate, J. Arteaga, J, García and R. Flores, “Error estimation and mesh adaptivity in incompressible viscous flows using a residual power approach”, Submitted to ''Comput. Meth. Appl. Mech. Engrg.'', 2004.

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

<div id="cite-65"></div>
'''[[#citeF-65|[65]]]'''  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-66"></div>
'''[66]''' 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 Yacht Design Conference'', Auckland, 4&#8211;6, December, 2002.

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

<div id="cite-68"></div>
'''[68]'''  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-69"></div>
'''[[#citeF-69|[69]]]'''  R. Radovitzki  and M. Ortiz,  “Lagrangian finite element analysis of a Newtonian flows”, ''Int.  J. Num. Meth. Engng.'', '''43''', 607&#8211;619, 1998.

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