Langrangian formulation for incompressible fluids using finite calculus and the finite element method

Latest revision as of 17:07, 10 December 2019

Abstract

We present a general formulation for incompressible fluid flow analysis using the finite element method (FEM) and a lagrangian description. The necessary stabilization for dealing with the incompressibility condition is introduced via the so called finite calculus (FIC) method. Both a quasi-implicit algorithm and a fractional step scheme are described. Examples of application of the lagrangian flow description are presented.

Keywords: Lagrangian formulation, incompressible fluid flow, finite calculus, finite element method.

1 INTRODUCTION

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

The problems originated by the convective terms dissapear if a lagrangian description is used. In the lagrangian formulation the motion of each individual flow particle is followed as in solid mechanics problems. The difficulty is tranferred now to the problem if adequately moving the mesh nodes. Indeed for large mesh motions remeshing is necessary after flow time steps.

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 pressure segregation schemes [3], artificial compressibility methods [4-6] and preconditioning techniques [7]. Other FEM schemes with good stabilization properties for the 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 [2,8]. More recently a general class of Galerkin FEM has been developed where the standard Galerkin variational form is extended with adequate residual-based terms in order to achieve a stabilized numerical scheme. Among the many methods of this kind we can name the Streamline Upwind Petrov Galerkin (SUPG) method [1,9-18] the Galerkin Least Square (GLS) method [19,20], the Taylor-Galerkin method [21], the Characteristic Galerkin method [22-24] and its variant the Characteristic Based Split (CBS) method [25,26] the pressure gradient operator technique [27] and the Subgrid Scale (SS) method [28-30]. A good review of these methods can be found in [31].

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

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

The lagrangian description has many advantages for tracking the displacement of fluid particles in flows where large motions of the fluid surface occur such in the case of breaking waves, splashing of water, filling of moulds, etc. As above mentioned, the positive feature of the lagrangian formulation is that the convective terms dissapear in the governing equations of the fluid flow; in return the updating of the mesh at almost every time step is now a necessity and an efficient algorithm for generation of a finite element mesh must be used.

The examples show the efficiency of the lagrangian formulation to solve fluid flow problems involving contact with moving solids, surface 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 ${\textstyle \Omega }$ of length ${\textstyle L}$ . The equation of balance of fluxes in a subdomain of size ${\textstyle d}$ belonging to ${\textstyle \Omega }$ (Figure 1) is written as

(1)

where ${\textstyle q_{A}}$ and ${\textstyle q_{B}}$ are the incoming and outgoing fluxes at points ${\textstyle A}$ and ${\textstyle B}$ , respectively. The flux ${\textstyle q}$ includes both convective and diffusive terms; i.e. ${\textstyle q=u\phi -k{d\phi \over dx}}$ , where ${\textstyle \phi }$ is the transported variable, ${\textstyle v}$ is the velocity and ${\textstyle k}$ is the diffusitivity of the material.

Figure 1: Equilibrium of fluxes in a balance domain of finite size

Let us express now the fluxes ${\textstyle q_{A}}$ and ${\textstyle q_{B}}$ in terms of the flux at an arbitrary point ${\textstyle C}$ within the balance domain (Figure 1). Expanding ${\textstyle q_{A}}$ and ${\textstyle q_{B}}$ in Taylor series around point ${\textstyle C}$ up to second order terms gives

q_{A}=q_{C}-d_{1}{\frac {dq}{dx}}\vert _{C}+{\frac {d_{1}^{2}}{2}}{\frac {d^{2}q}{dx^{2}}}\vert _{C}+O(d_{1}^{3})\quad ,\quad q_{B}=q_{C}+d_{2}{\frac {dq}{dx}}\vert _{C}+{\frac {d_{2}^{2}}{2}}{\frac {d^{2}q}{dx^{2}}}\vert _{C}+O(d_{2}^{3})

(2)

Substituting eq.(2) into eq.(1) gives after simplification

(3)

where ${\textstyle h=d_{1}-d_{2}}$ and all derivatives are computed at point ${\textstyle C}$ .

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

Distance ${\textstyle h}$ in eq.(3) can be interpreted as a free parameter depending, of course, on the location of point ${\textstyle C}$ (note that ${\textstyle -d\leq h\leq d}$ ). However, the fact that eq.(3) is the exact balance equation (up to second order terms) for any 1D domain of finite size and that the position of point ${\textstyle C}$ is arbitrary, can be used to derive numerical schemes with enhanced properties simply by computing the characteristic length parameter from an adequate “optimality” rule.

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

(4)

with

(5)

where ${\textstyle Q}$ is the external source. For consistency a “finite” form of the Neumann boundary condition should be used. This can be readily obtained by invoking balance of fluxes in a domain of finite size next to the boundary ${\textstyle \Gamma _{q}}$ where the external (diffusive) flux is prescribed to a value ${\textstyle q_{p}}$ . The modified Neumann boundary condition can be written as [32]

(6)

The definition of the problem is completed with the standard Dirichlet condition prescribing the value of ${\textstyle \phi }$ at the boundary ${\textstyle \Gamma _{\phi }}$ .

The underlined terms in Eqs.(5) and (7) introduce the necessary stabilization in the discrete solution of the problem using whatever numerical scheme. For details see [32-41].

The starting point in the next section are the FIC equation for a viscous incompressible fluid written in a lagrangian frame of reference.

3 FIC EQUATIONS FOR VISCOUS INCOMPRESSIBLE FLOW USING A LAGRANGIAN DESCRIPTION

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

Momentum

(7)

Mass balance

(8)

where

(9)

(10)

Above ${\textstyle u_{i}}$ is the velocity along the ith global axis, ${\textstyle \rho }$ is the (constant) density of the fluid, ${\textstyle p}$ is the absolute pressure (defined positive in compression), ${\textstyle b_{i}}$ are the body forces and ${\textstyle s_{ij}}$ are the viscous deviatoric stresses related to the viscosity ${\textstyle \mu }$ by the standard expression

(11)

where ${\textstyle \delta _{ij}}$ is the Kronecker delta and the strain rates ${\textstyle {\dot {\varepsilon }}_{ij}}$ are

(12)

The FIC boundary conditions are

(13)

(14)

and the initial condition is ${\textstyle u_{j}=u_{j}^{0}}$ for ${\textstyle t=t_{0}}$ .

In Eqs.(13) and (14) ${\textstyle t_{i}}$ and ${\textstyle u_{j}^{p}}$ are surface tractions and prescribed displacements on the boundaries ${\textstyle \Gamma _{t}}$ and ${\textstyle \Gamma _{u}}$ , respectively, ${\textstyle n_{j}}$ are the components of the unit normal vector to the boundary and ${\textstyle \sigma _{ij}}$ are the total stresses given by ${\textstyle \sigma _{ij}=s_{ij}-\delta _{ij}p}$ . The sign in front the stabilization term in Eq.(13) is positive due to the definition of ${\textstyle r_{m_{i}}}$ in Eq.(9).

The ${\textstyle h_{i}'s}$ 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 [32].

Eqs.(7)–(14) are the starting point for deriving stabilized finite element methods for solving the incompressible Navier-Stokes equations in a lagrangian frame of reference using equal order interpolation for the velocity and pressure variables [37-39]. Application of the FIC formulations to meshless analysis of fluid flow problems using the finite point method can be found in [40,41].

3.1 Stabilized integral forms

From Eq.(7) it can be obtained (taking into account Eqs.(9) and (11))

{\partial r_{d} \over \partial x_{i}}=-{1 \over a_{i}}\left[{\bar {r}}_{m_{i}}-{h_{j} \over 2}{\partial r_{m_{i}} \over \partial x_{j}}\right]-{\rho u_{i}h_{k} \over 2a_{i}}{\partial r_{d} \over \partial x_{k}}\quad i,j=1,n_{d}\quad ,\quad k\not =i

(15)

where

(16.a)

and

(16.b)

Substituting Eq.(15) into Eq.(8) and retaining the terms involving the derivatives of ${\textstyle r_{m_{i}}}$ with respect to ${\textstyle x_{i}}$ only, leads to the following expression for the stabilized mass balance equation

(17)

with

(18)

The ${\textstyle \tau _{i}}$ 's in Eq.(17) are termed in the stabilization literature intrinsic time parameters. Similar values for ${\textstyle \tau _{i}}$ (usually ${\textstyle \tau _{i}=\tau }$ is taken) are used in other works from ad-hoc extensions of the 1D advective-diffusive problem [12-31]. It is remarkable that the intrinsic time parameters have been deduced here from the general FIC formulation and this shows the possibilities of the method.

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

Momentum

(19)

Mass balance

(20)

where ${\textstyle \delta u_{i}}$ and ${\textstyle q}$ are arbitrary weighting functions representing virtual velocity and virtual pressure fields. Integration by parts of the ${\textstyle r_{m_{i}}}$ terms leads to

(21.a)

(21.b)

The third integral in Eq.(21.a) is expressed as a sum of the element contributions to allow for discontinuities in the derivatives of ${\textstyle r_{m_{i}}}$ along the element interfaces.

Also in Eq.(21.a) we will neglect hereonwards the third integral by assuming that ${\textstyle r_{m_{i}}}$ is negligible on the boundaries. The deviatoric stresses and the pressure terms in the first integral of Eq.(21.a) are integrated by parts in the usual manner. The resulting equations are

Momentum

(22)

Mass balance

(23)

We note that in Eq.(22) we use now the standard form of the convective operator for incompressible flows (i.e. neglecting the contribution from the volumetric strain rate ${\textstyle \displaystyle {\partial u_{i} \over \partial x_{i}}}$ ). Also in Eq.(22)

3.2 Convective and pressure gradient projections

The computation of the residual terms can be simplified if we introduce now the pressure gradient projections ${\textstyle \pi _{i}}$ , defined as

(24)

We can express ${\textstyle r_{m_{i}}}$ in Eq.(23) in terms of the ${\textstyle \pi _{i}}$ 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 ${\textstyle r_{m_{i}}}$ vanishes (in average sense) for the form given by Eq.(24). This gives the final system of governing equation as:

(25)

(26)

(27)

with ${\textstyle i,j,k=1,n_{d}}$ . In Eqs.(27) ${\textstyle \delta \pi _{i}}$ are appropriate weighting functions and the ${\textstyle \tau _{i}}$ weights are introduced for convenience.

4 FINITE ELEMENT DISCRETIZATION

We choose ${\textstyle C^{\circ }}$ continuous linear interpolations of the velocities, the pressure and the pressure gradient projections ${\textstyle \pi _{i}}$ over three node triangles (2D) and four node tetrahedra (3D). The linear interpolations are written as

(28)

where ${\textstyle n=3}$ (4) for triangles (tetrahedra), ${\textstyle {\bar {(\cdot )}}^{j}}$ denotes nodal variables and ${\textstyle N_{j}}$ are the linear shape functions [1].

Substituting the approximations (28) into Eqs.(25–27) and choosing the Galerking form with ${\textstyle \delta u_{i}=q=\delta \pi _{i}=N_{i}}$ leads to following system of discretized equations

(29.a)

(29.b)

(29.c)

where the element contributions are given by (for 2D problems)

{\begin{array}{l}\displaystyle M_{ij}=\int _{\Omega ^{e}}\rho N_{i}N_{j}d\Omega \quad ,\quad \displaystyle {\boldsymbol {K}}_{ij}=\int _{\Omega ^{e}}{\boldsymbol {B}}_{i}^{T}{\boldsymbol {D}}{\boldsymbol {B}}_{j}d\Omega \\\\\displaystyle {\boldsymbol {G}}_{ij}=\int _{\Omega ^{e}}({\boldsymbol {\nabla }}N_{i})N_{j}d\Omega \quad ,\quad {\boldsymbol {\nabla }}=\left[{\partial  \over \partial x_{1}},{\partial  \over \partial x_{2}}\right]^{T}\\\\\displaystyle L_{ij}=\int _{\Omega ^{e}}{\boldsymbol {\nabla }}^{T}N_{i}[\tau ]{\boldsymbol {\nabla }}N_{j}d\Omega \quad ,\quad [\tau ]=\left[{\begin{matrix}\tau _{1}&0\\0&\tau _{2}\\\end{matrix}}\right]\\\\\displaystyle {\boldsymbol {Q}}=[{\boldsymbol {Q}}^{1},{\boldsymbol {Q}}^{2}]\quad ,\quad \displaystyle Q_{ij}^{k}=\int _{\Omega ^{e}}\tau _{k}{\partial N_{i} \over \partial x_{k}}N_{j}d\Omega \\\\\displaystyle {\hat {\boldsymbol {M}}}=\left[{\begin{matrix}{\hat {\boldsymbol {M}}}^{1}&{\boldsymbol {0}}\\{\boldsymbol {0}}&{\hat {\boldsymbol {M}}}^{2}\\\end{matrix}}\right]\quad ,\quad {\hat {\boldsymbol {M}}}_{ij}=\int _{\Omega ^{e}}\tau _{k}N_{i}N_{j}d\Omega \\\\\displaystyle {\boldsymbol {f}}_{i}=\int _{\Omega ^{e}}N_{i}{\boldsymbol {b}}d\Omega +\int _{\Gamma ^{e}}N_{i}{\boldsymbol {t}}d\Gamma \quad ,\quad {\boldsymbol {b}}=[b_{1},b_{2}]^{T}\quad ,\quad {\boldsymbol {t}}=[t_{1},t_{2}]^{T}\end{array}}

(30)

with ${\textstyle i,j=1,n}$ and ${\textstyle k,l=1,n_{d}}$ .

In above ${\textstyle {\boldsymbol {B}}_{i}}$ is the standard strain rate matrix and ${\textstyle {\boldsymbol {D}}}$ the deviatoric constitutive matrix (assuming ${\textstyle \displaystyle {\partial u_{i} \over \partial x_{i}}=0}$ ). For 2D problems

{\boldsymbol {B}}_{i}=\left[{\begin{matrix}\displaystyle {\partial N_{i} \over \partial x_{1}}&0\\0&\displaystyle {\partial N_{i} \over \partial x_{2}}\\\displaystyle {\partial N_{i} \over \partial x_{2}}&\displaystyle {\partial N_{i} \over \partial x_{1}}\\\end{matrix}}\right]\quad ,\quad {\boldsymbol {D}}=\mu \left[{\begin{matrix}2&0&0\\0&2&0\\0&0&1\\\end{matrix}}\right]

(31)

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

\left[{\begin{matrix}{\boldsymbol {K}}&-{\boldsymbol {G}}&{\boldsymbol {0}}\\-{\boldsymbol {G}}^{T}&-{\boldsymbol {L}}&-{\boldsymbol {Q}}\\{\boldsymbol {0}}&-{\boldsymbol {Q}}^{T}&-{\hat {\boldsymbol {M}}}\\\end{matrix}}\right]\left\{{\begin{matrix}{\bar {\boldsymbol {u}}}\\{\bar {\boldsymbol {p}}}\\{\bar {\boldsymbol {\pi }}}\\\end{matrix}}\right\}=\left\{{\begin{matrix}{\boldsymbol {f}}\\{\boldsymbol {0}}\\{\boldsymbol {0}}\\\end{matrix}}\right\}

(32)

The system is symmetric and always positive definite and therefore leads to a non singular solution. We note that this property holds for any interpolation function chosen for ${\textstyle {\bar {\boldsymbol {u}}},{\bar {\boldsymbol {p}}}}$ and ${\textstyle {\bar {\boldsymbol {\pi }}}}$ , therefore overcoming the Babuŝka-Brezzi (BB) restrictions [1].

A reduced velocity-pressure formulation can be obtained by eliminating the pressure gradient projection variables ${\textstyle {\bar {\boldsymbol {\pi }}}}$ from the last equation to give

\left[{\begin{matrix}{\boldsymbol {K}}&-{\boldsymbol {G}}\\-{\boldsymbol {G}}^{T}&-({\boldsymbol {L}}-{\boldsymbol {Q}}{\hat {\boldsymbol {M}}}^{-1}{\boldsymbol {Q}}^{T})\\\end{matrix}}\right]\left\{{\begin{matrix}{\bar {\boldsymbol {u}}}\\{\bar {\boldsymbol {p}}}\\\end{matrix}}\right\}=\left\{{\begin{matrix}{\boldsymbol {f}}\\{\boldsymbol {0}}\\\end{matrix}}\right\}

(33)

The reduction process is simplified by using a diagonal form of matrix ${\textstyle {\hat {\boldsymbol {M}}}}$ . Obviously above reduction is also applicable to the transient case.

4.1 Quasi-implicit scheme

The solution process can be advanced in time in a (quasi-nearly) implicit iterative manner using the following scheme.

Step 1

(34)

Step 2

(35)

Step 3

(36)

where ${\textstyle 0\leq \theta _{i}\leq 1}$ .

Steps 1–3 are repeated until converged values of the velocities, ${\textstyle u_{i}}$ the pressure ${\textstyle p}$ and the pressure gradient projection variables ${\textstyle \pi _{i}}$ is obtained. The process within a time step increment is completed with steps 4 and 5 described below.

Step 4

Update the mesh nodes in a lagrangian manner as

(37)

Step 5

Generate a new mesh. This can be effectively performed using the procedure described in Section 5. Indeed the mesh regeneration can take place after a prescribed number of time steps or when the nodal displacements induce significant distorsions in some element shapes.

Details of the treatment of the boundary conditions in the lagrangian flow formulation can be found in [51,52].

Some practical remarks

In Eqs.(34–36) ${\textstyle {\bar {(\cdot )}}^{n,i}}$ denote nodal values at the ${\textstyle n}$ th time step and the ith iteration. Note that ${\textstyle {\boldsymbol {A}}^{n+\theta _{1},i-1}\equiv {\boldsymbol {A}}({\bar {\boldsymbol {u}}}^{n+\theta _{1},i-1})}$ etc. Also ${\textstyle (\cdot )^{n+\theta _{i},0}\equiv (\cdot )^{n}}$ for the computations in step 1 at the onset of the iterations.

Steps 1 and 3 can be solved explicitely by choosing a lumped (diagonal) form of matrices ${\textstyle M}$ and ${\textstyle {\hat {\boldsymbol {M}}}}$ . In this manner the main computational cost is the solution of step 2 involving the inverse of a Laplacian matrix. This can be solved very effectively using an iterative method such as the conjugate gradient method or similar.

For ${\textstyle \theta _{i}\not =0}$ the iterative proces is unavoidable. The iterations follow until convergence is reached in an adequate error norm in terms of the velocity and pressure variables, or the residuals ${\textstyle r_{m_{i}}}$ and ${\textstyle r_{d}}$ . Indeed some ot the ${\textstyle \theta _{i}}$ 's in Eqs.(34)–(36) can be made equal to zero. Note that for ${\textstyle \theta _{2}=0}$ the algorithm is unconditionally unstable. A particularity interesting and simple semi-implicit form is obtained by making ${\textstyle \theta _{1}=\theta _{2}=\theta _{3}=0}$ . Now all steps can be solved explicitely with exception of Step 2 for the pressure, which still requires the solution of a simultaneous system of equations.

Convergence of this solution scheme is however difficult for some problems. An enhanced version of the algorithm can be obtained by simply adding the term ${\textstyle {\hat {\boldsymbol {L}}}({\bar {\boldsymbol {p}}}^{n+1,i}-{\bar {\boldsymbol {p}}}^{n+1,i-1})}$ where ${\textstyle {\hat {L}}_{ij}=\Delta t\int _{\Omega ^{e}}{\boldsymbol {\nabla }}^{T}N_{i}{\boldsymbol {\nabla }}N_{j}d\Omega }$ to the equation for the computation of the pressure in the second step. The new term acts as a preconditioner of the pressure equation given now by

(38)

Note that the added term vanishes for the converged solution (i.e. when ${\textstyle {\bar {\boldsymbol {p}}}^{n+1,i}={\bar {\boldsymbol {p}}}^{n+1,i-1}}$ ).

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

4.2 Fractional step method

An effective algorithm can be obtained by splitting the pressure from the momentum equations as follows

(39.a)

(39.b)

In Eq.(39.a) ${\textstyle \alpha }$ is a variable taking values equal to zero or one. For ${\textstyle \alpha =0}$ , ${\textstyle \delta p\equiv p^{n+1}}$ and for ${\textstyle \alpha =1}$ , ${\textstyle \delta p=\Delta p}$ . Note that in both cases the sum of Eqs.(39.a) and (39.b) gives the time discretization of the momentum equations with the pressures computed at ${\textstyle t^{n+1}}$ . The value of ${\textstyle {\bar {\boldsymbol {u}}}^{n+1}}$ from Eq.(39.b) is substituted now into Eq.(29.b) to give

(40.a)

The product ${\textstyle {\boldsymbol {G}}^{T}{\boldsymbol {M}}^{-1}{\boldsymbol {G}}}$ can be approximated by a laplacian matrix, i.e.

(40.b)

A semi-implicit algorithm can thus be derived as follows.

Step 1 Compute ${\textstyle {\boldsymbol {u}}^{*}}$ explicitely from Eq.(39.a) with ${\textstyle {\boldsymbol {M}}={\boldsymbol {M}}_{d}}$ where subscript ${\textstyle d}$ denotes hereonwards a diagonal matrix.

Step 2 Compute ${\textstyle \delta {\bar {\boldsymbol {p}}}}$ and ${\textstyle {\boldsymbol {p}}^{n+1}}$ from Eq.(40.a) as

(41.a)

(41.b)

Step 3 Compute ${\textstyle {\bar {\boldsymbol {u}}}^{n+1}}$ explicitely from Eq.(39.b) with ${\textstyle {\boldsymbol {M}}={\boldsymbol {M}}_{d}}$

Step 4 Compute ${\textstyle {\bar {\boldsymbol {\pi }}}^{n+1}}$ explicitely from Eq.(36) as

(42)

Steps 5 and 6 coincide with steps 4 and 5 of Section 4.1.

This algorithm has an additional step than the iterative scheme of Section 4.1. The advantage is that now Steps 1 and 2 can be fully linearized by choosing ${\textstyle \theta _{1}=\theta _{2}=\theta _{3}=0}$ . Also the equation for the pressure variables in Step 2 has improved stabilization properties due to the additional laplacian matrix ${\textstyle {\hat {\boldsymbol {L}}}}$ .

The boundary conditions are applied as follows. No condition is applied in the computation of the fractional velocities ${\textstyle {\boldsymbol {u}}^{*}}$ in Eq.(39.a). The prescribed velocities at the boundary are applied when solving for ${\textstyle {\bar {\boldsymbol {u}}}^{n+1}}$ in step 3. The prescribed pressures at the boundary are imposed by making zero the pressure increments at the relevant boundary nodes and making ${\textstyle {\bar {\boldsymbol {p}}}^{n}}$ equal to the prescribed pressure values.

5 GENERATION OF A NEW MESH

One of the key points for success of the lagrangian flow formulation here described is the fast regeneration of a mesh at every time step on the basis of the position of the nodes in the space domain. In our work the mesh is generated using the so called extended Delaunay tesselation (EDT) presented in [44,47]. The EDT allows to generate meshes of elements with arbitrary polyhedrical shapes (combining triangles, quadrilaterals and other polygons in 2D and tetrahedra, hexahedra and arbitrary polyhedra in 3D) in a computing time of order ${\textstyle n}$ , being ${\textstyle n}$ the total number of nodes in the mesh. The shape functions of the elements can be simply obtained using the so called non-sibsonian interpolations. Details of the mesh generation procedure can be found in [44-47].

Once the new mesh has been generated at each time step the numerical solution is found using the finite element algorithm described in the paper.

The combination of elements with different geometrical shapes in the same mesh is one of the innovative aspects of the lagrangian formulation presented here.

6 IDENTIFICATION OF BOUNDARY SURFACES

One of the main problems in mesh generation is the correct definition of the boundary domain. Sometimes, boundary nodes are explicitly defined as special nodes, which are different from internal nodes. In other cases, the total set of nodes is the only information available and the algorithm must recognize the boundary nodes. Such is the case in lagrangian flow methods in which, at each time step, new nodal positions are obtained and the boundary-surface must be recognized using the new positions of the nodes.

The use of the Extended Delaunay partition makes it easier to recognize boundary nodes.

Considering that the nodes follow a variable ${\textstyle h(x)}$ distribution, where ${\textstyle h(x)}$ is the minimum distance between two nodes, the following criterion has been used:

All nodes on an empty sphere with a radius bigger than, are considered as boundary nodes (See Figure 2).

Thus, ${\textstyle \alpha }$ is a parameter close to, but greater than one. Note that this criterion is coincident with the Alpha Shape concept [47].

Figure 2: Contour recognition: Empty circles with radius

define the boundary particles.

Once a decision has been made concerning which of the nodes are on the boundaries, the boundary surface must be defined. It is well known that in 3-D problems the surface fitting a number of nodes is not unique. For instance, four boundary nodes on the same sphere may define two different boundary surfaces, a concave one and convex one.

In this work, the boundary surface is defined with all the polyhedral surfaces having all their nodes on the boundary and belonging to just one polyhedron.

The correct boundary surface may be important to define the correct normal external to the surface. Furthermore; in weak forms (Galerkin) such as those used here a correct evaluation of the volume domain is also important. Nevertheless, it must be noted that in the criterion proposed above, the error in the boundary surface definition is proportional to ${\textstyle h}$ . This is the error order accepted in a numerical method for a given node distribution. The only way to obtain more accurate boundary surface definition is by decreasing the distance between the nodes.

7 COMPUTATION OF THE CHARACTERISTIC LENGTHS

The evaluation of the stabilization parameters is one of the crucial issues in stabilized methods. Most of existing methods use expressions which are direct extensions of the values obtained for the simplest 1D case. In this work we will accept the so called SUPG assumption, i.e. to admit that vector ${\textstyle {\boldsymbol {h}}}$ has the direction of the velocity field [32-37]. This gives

(43)

where ${\textstyle u=\vert {\boldsymbol {u}}\vert }$ and ${\textstyle h_{s}}$ as the “streamline” characteristic length

(44)

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

A more consistent evaluation of ${\textstyle h}$ based on a diminishing residual technique can be found in [37].

8 EXAMPLES

All examples presented here have been obtained withthe fractional step algorithm described in Section 4.2.

8.1 Collapse of a water column

The first problem solved with the lagrangian formulation is the study of the collapse of a water column. This problem was solved by Koshizu and Oka [43] both experimentally and numerically. It has became a classical example to test the validation of the lagrangian formulation in fluid flows. The water is initially located on the left supported by a removable board. The collapse starts at time ${\textstyle t=0}$ , when the removable board is removed. Viscosity and surface tension are neglected. Figure 3 shows the point positions at different time steps. The dark points represent the free-surface detected with the algorithms described in Section 5. The internal points are gray and the fixed points are black.

The water is running on the bottom wall until, near ${\textstyle 0.3}$ sec, it impinges on the right vertical wall. Breaking waves appear at ${\textstyle 0.6}$ sec. Around ${\textstyle 1}$ sec. the water reaches the left wall. Agreement with the experimental results of [43] both in the shape of the free surface as well as in the time evolution are excellent.

Figure 3: Water column collapse at different time steps.

Figure 3 shows the point positions at different time steps. The darker points represent the free-surface detected.

The water is running on the bottom wall until, near ${\textstyle 0.3}$ sec, it impinges on the right vertical wall. Breaking waves appear at ${\textstyle 0.6}$ sec. Around ${\textstyle t=1}$ sec. the main water wave reaches the left wall again Agreement with the experimental results of ref. [43] both in the shape of the free surface as well as in the time development are excellent.

The 3D solution of the same problem is shown in Figure 4.




Figure 4: Dam collapse in a 3D domain. Different time step.

8.2 Semi-submerged rotating water mill

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

Figure 5: Rotating water mill.

8.3 Sloshing problems

The simple problem of the free oscillation of an incompressible liquid in a container is considered next. Numerical solutions for this problem can be found in several references [44]. This problem is interesting because there is an analytical solution for small amplitudes. Figure 6 shows a schematic of the problem, and represents also the point distribution in the initial position. The dark points represent the fixed points in which the velocity is fixed to zero.

Figure 6: Sloshing. Initial point distribution.

Figure 7: Sloshing: Comparison of the numerical and analytical solution.

Figure 7 shows the variation in time of the amplitude compared with the analytical results for the near in viscid case. Little numerical viscosity is observed on the phase wave and amplitude in spite of the relative poor point distribution.

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 can be separated from the fluid domain due to their large velocity. Figure 8 shows the numerical results obtained with the method presented in this paper for larger sloshing amplitudes. Breaking waves as well as separation effects can be seen on the free-surface. This particular and very complicated effect is apparently well represented by this model.

Figure 8: Sloshing: Different time step for large amplitudes.

Figure 9: Sloshing: Different time step for 3D domains.

In order to test the potentiality of the method in a 3D domain, the same sloshing problem was solved as a 3D problem. Figure 9 show the different point position 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 sphere representation is only used in order to improve the visualization of the numerical results.

8.4 Wave breaking on a beach

A simulation of the propagation of a water wave and its breaking due to shoaling over a plane slope is presented next. This example was numerically studied in [44] with a lagrangian formulation using directly the standard Finite Element Method with remeshing. There is also an analytical solution for a simplified approximation that is used for comparisons [45]. Figure 10 shows the initial point distribution and Figure 11 a comparison with the analytical free-surface at different time step. The geometry of the problem as well as a discussion of the analytical solution may be found in [44].

Figure 10: Wave breaking on a beach. Initial geometry and point positions.

Initially (Figures 11.a. and 11.b.) the wave travels over a constant depth bottom towards the slope with no ostensible change of shape. Strongly non-linear effects appear when the wave hits the slope (Fig. 11.c.). The crest of the wave accelerates until it reaches the shore (Fig. 11.d.). At this time the comparisons with the analytical solution are in agreement only in the wave position. The shape of the wave obtained with the numerical solution is totally different. The reason is that the analytical solution gives symmetrical shape waves, which are not physical, before the breaking process. Subsequently, a water jet is formed at the crest plunge making the breaking wave (Figs. 11.e. and 11.f.) and coming in contact with the nearly still surface of the water ahead. In Ref. 44 the evaluation is stopped before this contact point. Using the methodology proposed in this paper, the analysis may be continued until the end. In Figs. 11.g. and 11.h. the wave finally hits a lateral wall (introduced in the model to stop the lateral effects) producing drop separations, and then coming back toward to the left as a new wave.

Figure 11: Wave breaking on a beach. Comparison with analytical results at different time steps. Top: Numerical solution. Bottom: Analytical solution.

The ability of the model to accurately simulate the various stages of the wave breaking is noteworthy.

Nevertheless, a 2D domain is an easy case, and may be solved acceptably with any mesh generator. The true problems appears in a 3D domain, where the mesh generation is complicated with presence of slivers an other geometric mesh generation problems. In order to show the power of the tool presented, the same problem was solved in a 3D domain.

To transform the wave breaking described before in a true 3D problem, the initial position of the wave was introduced having an oblique angle with the beach line. In this way, a 3D effect appears. When the wave hits the slope, the crest of the wave accelerates differently accordingly with the depth, inducing the wave to correct its oblique position and break parallel to the beach. The results may be seen in Figure 12 for different time steps.

Figure 12: Breaking wave on a beach: Oblique wave on a 3D domain.

8.5 Solid floating on a free-surface

The following example shown schematically in Figure 13 represents a very interesting problem of fluid structure interaction when there is a weak interaction between the fluid and a large rigid deformation of the structure. In this case, there is also a free-surface problem, representing a schematic case of see-keeping in ship hydrodynamics.

Figure 13: Solid floating on a free-surface. Initial geometry and point distribution.

The example shows an initially stationary recipient with a floating piece of wood in which a wave is produced on the left side. The wave intercepts the wood piece producing a breaking wave and moving the floating wood. In this example the solid was represented by very viscous flows with a viscosity parameter order ten times greater than the water viscosity. Figure 14 shows the pressure contours and the free-surface position for different time steps.

Figure 14: Solid floating on a free-surface. Pressure contours and free-surface positions for different time step.

8.6 Solid cube falling in a recipient with water

This last example is also a case of fluid structure interaction. The solid is initially totally free and is falling down into a recipient with a fluid. In this example, the solid was modeled as a boundary condition for the fluid. Once the pressure and the viscous forces have been evaluated in the fluid, the solid is accelerated using Newton law. The solid has a mass and a gravity force concentrate in its gravity center. The solid is considered to be ligh compared to the liquid weight.

At the beginning the solid falls free due to the gravity forces (Figure 15). Once in contact with the water free-surface, the alpha-shape method recognizes the different boundary contours. The pressure and viscous forces are evaluated in all the domain and in particular on the solid contour. The flow forces introduce a negative acceleration in the vertical direction until , once the solid is completely inside the water, the vertical velocity becomes zero. Then, Arquimides forces bring the solid up to the free-surface. It is interesting to observe that there is a rotation of the solid. The reason is that the center of the floating forces is higher in the rotated position than in the initial ones.

Figure 15: Solid cube falling into a recipient with water. Free- surface position at different time steps.

9 CONCLUSIONS

The finite calculus expression of the fluid mechanics equations written in a lagrangian frame of reference is a good starting point for deriving stabilized finite element method for solving a variety of fluid flow problems. Both monolithic and fractional step algorithms with intrinsic stabilization properties can be readily derived as shown here. The lagrangian formulation allows to solve in an effective manner fluid flow problems involving large motions of the free surface and complex fluid-structure interactions.

ACKNOWLEDGEMENTS

Thanks are also given to Mr. Romain Aubry for many useful discussions.

References

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

[2] C. Hirsch, Numerical computation of internal and external flow, J. Wiley, Vol. 1 1988, Vol. 2, 1990.

[3] A.J. Chorin, “A numerical solution for solving incompressible viscous flow problems” J. Comp. Phys., 2, 12–26, 1967.

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

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

[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–42, 1996.

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

[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–143, 1977.

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

[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–46, 1996.

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

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

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

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

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

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

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

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

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

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

[21] J. Donea, “A Taylor-Galerkin method for convective transport problems”, Int. J. Num. Meth. Engng., 20, 101–119, 1984.

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

[23] O. Pironneau, “On the transport-diffusion algorithm and its applications to the Navier-Stokes equations”, Numer. Math., 38, 309, 1982.

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

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

[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–112, 2002.

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

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

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

[30] R. Codina, “Stabilized finite element approximation of transient incompressible flows using orthogonal subscales”, Comput. Methods Appl. Mech. Engrg., 191, 4295–4321, 2002.

[31] J. Donea and A. Huerta, “Finite element method for flow problems”, J. Wiley, 2003.

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

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

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

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

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

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

[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–660, (2001).

[39] E. Oñate, “Possibilities of finite calculus in computational mechanics”, Submitted to Int. J. Num. Meth. Engng., 2002.

[40] E. Oñate and S. Idelsohn, A mesh free finite point method for advective-diffusive transport and fluid flow problems,, Computational Mechanics, 21, 283–292, 1988.

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

[42] J. García and E. Oñate, “An unstructured finite element solver for ship hydrodynamic problems”, in J. Appl. Mech., 70, 18–26, January 2003.

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

[44] S.R. Idelsohn, E. Oñate, N. Calvo and F. del Pin, “The meshless finite element method”, in Int. J. Num. Meth. Engng., 58, 4, 2003.

[45] S.R. Idelsohn, E. Oñate, F. Del Pin and N. Calvo, “Lagrangian formulation: the only way to solve some free-surface fluid mechanics problems”, Fith World Congress on Computational Mechanics, Mang HA, Rammerstorfer FG and Eberhardsteiner J. (eds), July 7–12, Viena, Austria, 2002, Web...

[46] S.R. Idelsohn, E. Oñate and F. Del Pin, “A lagrangian meshless finite element method applied to fluid-structure interaction problems”, in Computer and Structures, 81, 655–671, 2003.

[47] S.R. Idelsohn, N. Calvo and E. Oñate. Polyhedrization of an arbitrary point set. Comput. Method Appl. Mech. Engng.. Accepted for publication 2003.

[48] H. Edelsbrunner and E.P. Mucke. Three dimensional alpha shapes. ºit ACM Trans. Graphics, 13, 43–72, 1999.

[49] S. Koshizuka and Y. Oka. Moving particle semi-implicit method for fragmentation of incompressible fluid. Nuclear Engineering Science, 123, 421–434, 1996.

[50] R. Radovitzki and M. Ortiz. Lagrangian finite element analysis of Newtonian flows. Int. J. Num. Meth. Engng., 43, 607–619, 1998.

[51] E.V. Laitone. The second approximation to cnoidal waves. Journal of Fluid Mechanics, 9, 430, 1960.

Latest revision as of 17:07, 10 December 2019

Abstract

1 INTRODUCTION

2 THE FINITE CALCULUS METHOD

3 FIC EQUATIONS FOR VISCOUS INCOMPRESSIBLE FLOW USING A LAGRANGIAN DESCRIPTION

3.1 Stabilized integral forms

3.2 Convective and pressure gradient projections

4 FINITE ELEMENT DISCRETIZATION

4.1 Quasi-implicit scheme

4.2 Fractional step method

5 GENERATION OF A NEW MESH

6 IDENTIFICATION OF BOUNDARY SURFACES

7 COMPUTATION OF THE CHARACTERISTIC LENGTHS

8 EXAMPLES

8.1 Collapse of a water column

8.2 Semi-submerged rotating water mill

8.3 Sloshing problems

8.4 Wave breaking on a beach

8.5 Solid floating on a free-surface

8.6 Solid cube falling in a recipient with water

9 CONCLUSIONS

ACKNOWLEDGEMENTS

References

Document information

Document Score

Share this document

Keywords

claim authorship

@@ Line 36: / Line 36: @@
 We present a general formulation for incompressible fluid flow analysis using the finite element method (FEM) and a lagrangian description. The necessary stabilization for dealing with the incompressibility condition is introduced via the so called finite calculus (FIC) method. Both a quasi-implicit algorithm and a fractional step scheme are described. Examples of application of the lagrangian flow description are presented.
-'''keywords''' Lagrangian formulation, incompressible fluid flow, finite calculus, finite element method.
+'''Keywords:''' Lagrangian formulation, incompressible fluid flow, finite calculus, finite element method.
 ==1 INTRODUCTION==
-The development of efficient and robust numerical methods such as the finite element, finite volume and finite difference methods for analysis of incompressible flows has been a subject of intensive research in last decades <span id='citeF-1'></span>[[#cite-1|1]]. Much effort has been spent in developing the so called stabilized numerical methods overcoming the two main sources of instability in incompressible flow analysis, namely those originated by the high values of the convective terms and those induced by the difficulty in satisfying the incompressibility conditions.
+The development of efficient and robust numerical methods such as the finite element, finite volume and finite difference methods for analysis of incompressible flows has been a subject of intensive research in last decades<span id='citeF-1'></span><span id='citeF-1'></span> [[#cite-1|[1]],[[#cite-2|2]]]. Much effort has been spent in developing the so called stabilized numerical methods overcoming the two main sources of instability in incompressible flow analysis, namely those originated by the high values of the convective terms and those induced by the difficulty in satisfying the incompressibility conditions.
 The problems originated by the convective terms dissapear if a lagrangian description is used. In the lagrangian formulation the motion of each individual flow particle is followed as in solid mechanics problems. The difficulty is tranferred now to the problem if adequately moving the mesh nodes. Indeed for large mesh motions remeshing is necessary after flow time steps.
-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 pressure segregation schemes <span id='citeF-3'></span>[[#cite-3|3]], artificial compressibility methods <span id='citeF-4'></span>[[#cite-4|4]] and preconditioning techniques <span id='citeF-7'></span>[[#cite-7|7]]. Other FEM schemes with good stabilization properties for the 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]]. More recently a general class of Galerkin FEM has been developed where the standard Galerkin variational form is extended with adequate residual-based terms in order to achieve a stabilized numerical scheme. Among the many methods of this kind 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 technique <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]].
+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 pressure segregation schemes <span id='citeF-3'></span>[[#cite-3|[3]]], artificial compressibility methods <span id='citeF-4'></span><span id='citeF-5'></span><span id='citeF-6'></span>[[#cite-4|[4]]-[[#cite-6|6]]] and preconditioning techniques <span id='citeF-7'></span>[[#cite-7|[7]]]. Other FEM schemes with good stabilization properties for the 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><span id='citeF-8'></span>[[#cite-2|[2]],[[#cite-8|8]]]. More recently a general class of Galerkin FEM has been developed where the standard Galerkin variational form is extended with adequate residual-based terms in order to achieve a stabilized numerical scheme. Among the many methods of this kind we can name the Streamline Upwind Petrov Galerkin (SUPG) method <span id='citeF-1'></span><span id='citeF-9'></span><span id='citeF-18'></span>[[#cite-1|[1]],[[#cite-9|9]]-[[#cite-18|18]]] the Galerkin Least Square (GLS) method <span id='citeF-19'></span><span id='citeF-20'></span>[[#cite-19|[19]],[[#cite-20|20]]], the Taylor-Galerkin method <span id='citeF-21'></span>[[#cite-21|[21]]], the Characteristic Galerkin method <span id='citeF-22'></span><span id='citeF-24'></span>[[#cite-22|[22]]-[[#cite-23|24]]] and its variant the Characteristic Based Split (CBS) method <span id='citeF-25'></span><span id='citeF-26'></span>[[#cite-25|[25]],[[#cite-25|26]]] the pressure gradient operator technique <span id='citeF-27'></span>[[#cite-27|[27]]] and the Subgrid Scale (SS) method <span id='citeF-28'></span><span id='citeF-29'></span><span id='citeF-30'></span>[[#cite-28|[28]]-[[#cite-30|30]]]. A good review of these methods can be found in <span id='citeF-31'></span>[[#cite-31|[31]]].
-In this paper a stabilized finite element formulation for incompressible lagrangian flows is derived in a different manner. The starting point are the modified governing differential equations of the fluid flow problem formulated via a finite calculus (FIC) approach <span id='citeF-32'></span>[[#cite-32|32]]. The FIC method is based in invoking the balance of fluxes in a fluid domain of finite size. This introduces naturally additional terms in the classical differential equations of infinitesimal fluid mechanics which are a function of the balance domain dimensions. The new terms in the modified governing equations provide the necessary stabilization to the discrete equations obtained via the standard Galerkin finite element method <span id='citeF-33'></span>[[#cite-33|33]].
+In this paper a stabilized finite element formulation for incompressible lagrangian flows is derived in a different manner. The starting point are the modified governing differential equations of the fluid flow problem formulated via a finite calculus (FIC) approach <span id='citeF-32'></span>[[#cite-32|[32]]]. The FIC method is based in invoking the balance of fluxes in a fluid domain of finite size. This introduces naturally additional terms in the classical differential equations of infinitesimal fluid mechanics which are a function of the balance domain dimensions. The new terms in the modified governing equations provide the necessary stabilization to the discrete equations obtained via the standard Galerkin finite element method <span id='citeF-33'></span><span id='citeF-39'></span>[[#cite-33|[33]]-[[#cite-33|39]]].
 The layout of the chapter is the following. In the next section, the main concepts of the FIC approach are introduced via a simple 1D convection-diffusion model problem. Then the basic FIC equations for incompressible lagrangian flow problems are presented. The finite element discretization is introduced and the resulting matrix formulation is detailed. Both monolithic and fractional step schemes for the transient solution are presented.
@@ Line 56: / Line 56: @@
 ==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
+We will consider a convection-diffusion problem in a 1D domain <math display="inline">\Omega </math> of length <math display="inline">L</math>. The equation of balance of fluxes in a subdomain of size <math display="inline">d</math> belonging to <math display="inline">\Omega </math> (Figure [[#img-1|1]]) is written as
+<span id="eq-1"></span>
 {| class="formulaSCP" style="width: 100%; text-align: left;"
 |-
@@ Line 73: / Line 74: @@
 {| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
 |-
-|[[Image:Draft_Samper_800412471_9634_img-1.JPG|516px|Equilibrium of fluxes in a  balance domain of finite size]]
+|style="padding: 10px;"|[[File:Draft_Samper_187181919_2981_img-1.JPG|516px|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
+| colspan="1" style="padding: 10px;"| '''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
+Let us express now the fluxes <math display="inline">q_A</math> and <math display="inline">q_B</math> in terms of the flux at an arbitrary point <math display="inline">C</math> within the balance domain (Figure [[#img-1|1]]). Expanding <math display="inline">q_A</math> and <math display="inline">q_B</math> in Taylor series around point <math display="inline">C</math> up to second order terms gives
+<span id="eq-2"></span>
 {| class="formulaSCP" style="width: 100%; text-align: left;"
 |-
@@ Line 90: / Line 92: @@
 |}
-Substituting eq.(2) into eq.(1) gives after simplification
+Substituting eq.([[#eq-2|2]]) into eq.([[#eq-1|1]]) gives after simplification
+<span id="eq-3"></span>
 {| class="formulaSCP" style="width: 100%; text-align: left;"
 |-
@@ Line 97: / Line 100: @@
 {| 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="text-align: center;" | <math>\frac{dq}{dx}-\underline{\frac{h}{2} \frac{d^2q}{dx^2}}=0 </math>
 |}
 | style="width: 5px;text-align: right;white-space: nowrap;" | (3)
@@ Line 104: / Line 107: @@
 where <math display="inline">h=d_1-d_2</math> and all derivatives are computed at point <math display="inline">C</math>.
-Standard calculus theory assumes  that the domain <math display="inline">d</math> is of infinitesimal size and the resulting balance equation is simply <math display="inline">{dq\over dx}=0</math>. We will relax this assumption and allow the balance domain to have a ''finite size''. The new balance equation (3) incorporates now the underlined term which introduces the ''characteristic length'' <math display="inline">h</math>. Obviously, accounting for higher order terms in eq.(2) would lead to new terms in eq.(3) involving higher powers of <math display="inline">h</math>.
+Standard calculus theory assumes  that the domain <math display="inline">d</math> is of infinitesimal size and the resulting balance equation is simply <math display="inline">{dq\over dx}=0</math>. We will relax this assumption and allow the balance domain to have a ''finite size''. The new balance equation ([[#eq-3|3]]) incorporates now the underlined term which introduces the ''characteristic length'' <math display="inline">h</math>. Obviously, accounting for higher order terms in eq.([[#eq-2|2]]) would lead to new terms in eq.([[#eq-3|3]]) involving higher powers of <math display="inline">h</math>.
-Distance <math display="inline">h</math> in eq.(3) can be interpreted as a free parameter depending, of course, on the location of point <math display="inline">C</math> (note that <math display="inline">-d\le h \le d</math>). However, the fact that eq.(3) is the exact balance equation (up to second order terms) for any 1D domain of finite size and that the position of point <math display="inline">C</math> is arbitrary, can be used to derive numerical schemes with enhanced properties simply by computing the characteristic length parameter from an adequate “optimality” rule.
+Distance <math display="inline">h</math> in eq.([[#eq-3|3]]) can be interpreted as a free parameter depending, of course, on the location of point <math display="inline">C</math> (note that <math display="inline">-d\le h \le d</math>). However, the fact that eq.([[#eq-3|3]]) is the exact balance equation (up to second order terms) for any 1D domain of finite size and that the position of point <math display="inline">C</math> is arbitrary, can be used to derive numerical schemes with enhanced properties simply by computing the characteristic length parameter from an adequate “optimality” rule.
-Equation (3) can be extended to account for source effects. The full stabilized equation can be then written in compact form as
+Equation ([[#eq-3|3]]) can be extended to account for source effects. The full stabilized equation can be then written in compact form as
+<span id="eq-4"></span>
 {| class="formulaSCP" style="width: 100%; text-align: left;"
 |-
@@ Line 115: / Line 119: @@
 {| 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="text-align: center;" | <math>r- \underline{\frac{h}{2} \frac{dr}{dx}}=0 </math>
 |}
 | style="width: 5px;text-align: right;white-space: nowrap;" | (4)
@@ Line 122: / Line 126: @@
 with
+<span id="eq-5"></span>
 {| class="formulaSCP" style="width: 100%; text-align: left;"
 |-
@@ Line 132: / Line 137: @@
 |}
-where <math display="inline">Q</math> is the external source. For consistency a “finite” form of the Neumann boundary condition should be used. This can be readily obtained by invoking balance of fluxes in a domain of finite size next to the boundary <math display="inline">\Gamma _q</math> where the external (diffusive) flux is prescribed to a value <math display="inline">q_p</math>. The modified Neumann boundary condition can be written as <span id='citeF-32'></span>[[#cite-32|32]]
+where <math display="inline">Q</math> is the external source. For consistency a “finite” form of the Neumann boundary condition should be used. This can be readily obtained by invoking balance of fluxes in a domain of finite size next to the boundary <math display="inline">\Gamma _q</math> where the external (diffusive) flux is prescribed to a value <math display="inline">q_p</math>. The modified Neumann boundary condition can be written as <span id='citeF-32'></span>[[#cite-32|[32]]]
+<span id="eq-6"></span>
 {| class="formulaSCP" style="width: 100%; text-align: left;"
 |-
@@ Line 139: / Line 145: @@
 {| 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="text-align: center;" | <math>k \frac{d\phi }{dx}+ q_p- \underline{\frac{h}{2}r}=0 \quad \hbox{at } \, \Gamma _q </math>
 |}
 | style="width: 5px;text-align: right;white-space: nowrap;" | (6)
@@ Line 146: / Line 152: @@
 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.(5) and (7) introduce the necessary stabilization in the discrete solution of the problem using whatever numerical scheme. For details see <span id='citeF-32'></span>[[#cite-32|32]].
+The underlined terms in Eqs.([[#eq-5|5]]) and ([[#eq-7|7]]) introduce the necessary stabilization in the discrete solution of the problem using whatever numerical scheme. For details see <span id='citeF-32'></span>[[#cite-32|[32]]-[[#cite-41|41]]].
 The starting point in the next section are the FIC equation for a viscous incompressible fluid written in a lagrangian frame of reference.
@@ Line 154: / Line 160: @@
 The FIC governing equations for a viscous incompressible fluid can be written in a lagrangian frame as
-===Momentum===
+'''Momentum'''
+<span id="eq-7"></span>
 {| class="formulaSCP" style="width: 100%; text-align: left;"
 |-
@@ Line 161: / Line 168: @@
 {| 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  </math>
+| style="text-align: center;" | <math>r_{m_i} - \underline{{1\over 2} h_j{\partial r_{m_i} \over \partial x_j}}=0  </math>
 |}
 | style="width: 5px;text-align: right;white-space: nowrap;" | (7)
 |}
-===Mass balance===
+'''Mass balance'''
+<span id="eq-8"></span>
 {| class="formulaSCP" style="width: 100%; text-align: left;"
 |-
@@ Line 173: / Line 181: @@
 {| 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  </math>
+| style="text-align: center;" | <math>r_d - \underline{{1\over 2} h_j {\partial r_d \over \partial x_j}}=0  </math>
 |}
 | style="width: 5px;text-align: right;white-space: nowrap;" | (8)
@@ Line 180: / Line 188: @@
 where
+<span id="eq-9"></span>
 {| class="formulaSCP" style="width: 100%; text-align: left;"
 |-
@@ Line 186: / Line 195: @@
 |-
 | style="text-align: center;" | <math>r_{m_i} = \rho {\partial u_i \over \partial t} + {\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)
+|}
+<span id="eq-10"></span>
+{| class="formulaSCP" style="width: 100%; text-align: left;"
+|-
+|
+{| style="text-align: left; margin:auto;width: 100%;"
 |-
 | style="text-align: center;" | <math> r_d = {\partial u_i \over \partial x_i}\qquad i,j = 1, n_d </math>
+|}
 | style="width: 5px;text-align: right;white-space: nowrap;" | (10)
 |}
-|}
 Above <math display="inline">u_i</math> is the velocity along the ith global axis, <math display="inline">\rho </math> is the (constant) density of the fluid, <math display="inline">p</math> is the absolute pressure (defined positive in compression), <math display="inline">b_i</math> are the body forces and <math display="inline">s_{ij}</math> are the viscous deviatoric stresses related to the viscosity <math display="inline">\mu </math> by the standard expression
+<span id="eq-11"></span>
 {| class="formulaSCP" style="width: 100%; text-align: left;"
 |-
@@ Line 207: / Line 226: @@
 where <math display="inline">\delta _{ij}</math> is the Kronecker delta and the strain rates <math display="inline">\dot \varepsilon _{ij}</math> are
+<span id="eq-12"></span>
 {| class="formulaSCP" style="width: 100%; text-align: left;"
 |-
@@ Line 219: / Line 239: @@
 The FIC boundary conditions  are
+<span id="eq-13"></span>
 {| class="formulaSCP" style="width: 100%; text-align: left;"
 |-
@@ Line 224: / Line 245: @@
 {| 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="text-align: center;" | <math>n_j \sigma _{ij} -t_i + \underline{{1\over 2} h_j n_j r_{m_i}}=0 \quad \hbox{on }\Gamma _t </math>
 |}
 | style="width: 5px;text-align: right;white-space: nowrap;" | (13)
 |}
+<span id="eq-14"></span>
 {| class="formulaSCP" style="width: 100%; text-align: left;"
 |-
@@ Line 241: / Line 263: @@
 and the initial condition is <math display="inline">u_j =u_j^0</math> for <math display="inline">t=t_0</math>.
-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).
+In Eqs.([[#eq-13|13]]) and ([[#eq-14|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.([[#eq-13|13]]) is positive due to the definition of <math display="inline">r_{m_i}</math> in Eq.([[#eq-9|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]].
+The <math display="inline">h_i's</math> in above equations are characteristic lengths of the domain where balance of momentum and mass is enforced. In Eq.([[#eq-13|13]]) these lengths define the domain where equilibrium of boundary tractions is established <span id='citeF-32'></span>[[#cite-32|[32]]].
-Eqs.(7)&#8211;(14) are the starting  point for deriving stabilized finite element methods for solving the incompressible Navier-Stokes equations in a lagrangian frame of reference using equal order interpolation for the velocity and pressure variables <span id='citeF-37'></span>[[#cite-37|37]]. Application of the FIC formulations to meshless analysis of fluid flow problems using the finite point method can be found in <span id='citeF-40'></span>[[#cite-40|40]].
+Eqs.([[#eq-7|7]])&#8211;([[#eq-14|14]]) are the starting  point for deriving stabilized finite element methods for solving the incompressible Navier-Stokes equations in a lagrangian frame of reference using equal order interpolation for the velocity and pressure variables <span id='citeF-37'></span>[[#cite-37|[37]]-[[#cite-39|39]]]. Application of the FIC formulations to meshless analysis of fluid flow problems using the finite point method can be found in <span id='citeF-40'></span>[[#cite-40|[40]],[[#cite-41|41]]].
 ===3.1 Stabilized integral forms===
-From Eq.(7) it can be obtained (taking into account Eqs.(9) and (11))
+From Eq.([[#eq-7|7]]) it can be obtained (taking into account Eqs.([[#eq-9|9]]) and ([[#eq-11|11]]))
+<span id="eq-15"></span>
 {| class="formulaSCP" style="width: 100%; text-align: left;"
 |-
@@ Line 263: / Line 286: @@
 where
+<span id="eq-16.a"></span>
 {| class="formulaSCP" style="width: 100%; text-align: left;"
 |-
@@ Line 270: / Line 294: @@
 | style="text-align: center;" | <math>a_i = {2\mu \over 3} </math>
 |}
-| style="width: 5px;text-align: right;white-space: nowrap;" | (16a)
+| style="width: 5px;text-align: right;white-space: nowrap;" | (16.a)
 |}
 and
+<span id="eq-16.b"></span>
 {| class="formulaSCP" style="width: 100%; text-align: left;"
 |-
@@ Line 282: / Line 307: @@
 | style="text-align: center;" | <math>\bar r_{m_i}=\rho {\partial u_i \over \partial t} + {\partial p \over \partial x_i}-{\partial  \over \partial x_j} (2\mu \dot \varepsilon _{ij}) - b_i</math>
 |}
-| style="width: 5px;text-align: right;white-space: nowrap;" | (16b)
+| style="width: 5px;text-align: right;white-space: nowrap;" | (16.b)
 |}
-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 expression for the stabilized mass balance equation
+Substituting Eq.([[#eq-15|15]]) into Eq.([[#eq-8|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 expression for the stabilized mass balance equation
-{| class="formulaSCP" style="width: 100%; text-align: left;"
+<span id="eq-17"></span>
+{| class="formulaSCP" style=" width: 100%; text-align: left;"
 |-
 |
-{| style="text-align: left; margin:auto;width: 100%;"
+{| style="text-align: left; margin:auto;"
 |-
-| style="text-align: center;" | <math>\displaystyle{r_d - \sum \limits _{i=1}^{n_d} \tau _i {\partial r_{m_i} \over \partial x_i}=0} </math>
+| style="text-align: center;   border: 1px solid; padding: 10px;" | <math>\displaystyle{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)
@@ Line 299: / Line 325: @@
 with
+<span id="eq-18"></span>
 {| class="formulaSCP" style="width: 100%; text-align: left;"
 |-
@@ Line 309: / Line 336: @@
 |}
-The <math display="inline">\tau _i</math>'s in Eq.(17) are termed in the stabilization literature ''intrinsic time parameters''. Similar values for <math display="inline">\tau _i</math> (usually <math display="inline">\tau _i =\tau </math> is taken) are used in other works from ad-hoc extensions of the 1D advective-diffusive problem <span id='citeF-12'></span>[[#cite-12|12]]. It is remarkable that the intrinsic time parameters have been deduced here from the general FIC formulation and this shows the possibilities of the method.
+The <math display="inline">\tau _i</math>'s in Eq.([[#eq-17|17]]) are termed in the stabilization literature ''intrinsic time parameters''. Similar values for <math display="inline">\tau _i</math> (usually <math display="inline">\tau _i =\tau </math> is taken) are used in other works from ad-hoc extensions of the 1D advective-diffusive problem <span id='citeF-12'></span>[[#cite-12|[12]]-[[#cite-31|31]]]. It is remarkable that the intrinsic time parameters have been deduced here from the general FIC formulation and this shows the possibilities of the method.
-The weighted residual form of the momentum and mass balance equations (Eqs.(7) and (17)) is written as
+The weighted residual form of the momentum and mass balance equations (Eqs.([[#eq-7|7]]) and ([[#eq-17|17]])) is written as
-===Momentum===
+'''Momentum'''
+<span id="eq-19"></span>
 {| class="formulaSCP" style="width: 100%; text-align: left;"
 |-
@@ Line 325: / Line 353: @@
 |}
-===Mass balance===
+'''Mass balance'''
+<span id="eq-20"></span>
 {| class="formulaSCP" style="width: 100%; text-align: left;"
 |-
@@ Line 339: / Line 368: @@
 where <math display="inline">\delta u_i</math> and <math display="inline">q</math> are arbitrary weighting functions representing virtual velocity and virtual pressure fields. Integration by parts of the <math display="inline">r_{m_i}</math> terms  leads to
+<span id="eq-21.a"></span>
 {| class="formulaSCP" style="width: 100%; text-align: left;"
 |-
@@ Line 346: / Line 376: @@
 | style="text-align: center;" | <math>\int _\Omega \delta u_i r_{m_i} + \int _{\Gamma _t} \delta u_i (\sigma _{ij} n_j - t_i)d\Gamma + \sum \limits _e \int _{\Omega ^e} {h_j\over 2}{\partial \delta u_i \over \partial x_j} r_{m_i} d\Omega =0 </math>
 |}
-| style="width: 5px;text-align: right;white-space: nowrap;" | (21a)
+| style="width: 5px;text-align: right;white-space: nowrap;" | (21.a)
 |}
+<span id="eq-21.b"></span>
 {| class="formulaSCP" style="width: 100%; text-align: left;"
 |-
@@ Line 356: / Line 387: @@
 | 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)
+| style="width: 5px;text-align: right;white-space: nowrap;" | (21.b)
 |}
-The third integral in Eq.(21a) is expressed as a sum of the element contributions to allow for discontinuities in the derivatives of <math display="inline">r_{m_i}</math> along the element interfaces.
+The third integral in Eq.([[#eq-21.a|21.a]]) is expressed as a sum of the element contributions to allow for discontinuities in the derivatives of <math display="inline">r_{m_i}</math> along the element interfaces.
-Also in Eq.(21a) we will neglect hereonwards the third integral by assuming that <math display="inline">r_{m_i}</math> is negligible on the boundaries. The deviatoric stresses and the pressure terms in the first integral of Eq.(21a) are integrated by parts in the usual manner. The resulting equations are
+Also in Eq.([[#eq-21.a|21.a]]) we will neglect hereonwards the third integral by assuming that <math display="inline">r_{m_i}</math> is negligible on the boundaries. The deviatoric stresses and the pressure terms in the first integral of Eq.([[#eq-21.a|21.a]]) are integrated by parts in the usual manner. The resulting equations are
-===Momentum===
+'''Momentum'''
+<span id="eq-22"></span>
 {| class="formulaSCP" style="width: 100%; text-align: left;"
 |-
@@ Line 375: / Line 407: @@
 |}
-===Mass balance===
+'''Mass balance'''
+<span id="eq-23"></span>
 {| class="formulaSCP" style="width: 100%; text-align: left;"
 |-
@@ Line 387: / Line 420: @@
 |}
-We note that in Eq.(22) we use now the standard form of the convective operator for incompressible flows (i.e. neglecting the contribution from the volumetric strain rate <math display="inline">\displaystyle {\partial u_i \over \partial x_i}</math>). Also in Eq.(22)
+We note that in Eq.([[#eq-22|22]]) we use now the standard form of the convective operator for incompressible flows (i.e. neglecting the contribution from the volumetric strain rate <math display="inline">\displaystyle {\partial u_i \over \partial x_i}</math>). Also in Eq.([[#eq-22|22]])
 {| class="formulaSCP" style="width: 100%; text-align: left;"
@@ Line 402: / Line 435: @@
 The computation of the residual terms can be simplified if we introduce now the pressure gradient projections <math display="inline">\pi _i</math>, defined as
+<span id="eq-24"></span>
 {| class="formulaSCP" style="width: 100%; text-align: left;"
 |-
@@ Line 412: / Line 446: @@
 |}
-We can express <math display="inline">r_{m_i}</math> in  Eq.(23) in terms of the <math display="inline">\pi _i</math>  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 the form given by Eq.(24). This gives the final system of governing equation as:
+We can express <math display="inline">r_{m_i}</math> in  Eq.([[#eq-23|23]]) in terms of the <math display="inline">\pi _i</math>  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 the form given by Eq.([[#eq-24|24]]). This gives the final system of governing equation as:
+<span id="eq-25"></span>
 {| class="formulaSCP" style="width: 100%; text-align: left;"
 |-
@@ Line 424: / Line 459: @@
 |}
+<span id="eq-26"></span>
 {| class="formulaSCP" style="width: 100%; text-align: left;"
 |-
@@ Line 434: / Line 470: @@
 |}
+<span id="eq-27"></span>
 {| class="formulaSCP" style="width: 100%; text-align: left;"
 |-
@@ Line 444: / Line 481: @@
 |}
-with <math display="inline">i,j,k=1,n_d</math>.  In Eqs.(27) <math display="inline">\delta \pi _i</math> are appropriate weighting functions and the <math display="inline">\tau _i</math> weights are introduced for convenience.
+with <math display="inline">i,j,k=1,n_d</math>. In Eqs.([[#eq-27|27]]) <math display="inline">\delta \pi _i</math> are appropriate weighting functions and the <math display="inline">\tau _i</math> weights are introduced for convenience.
 ==4 FINITE ELEMENT DISCRETIZATION==
@@ Line 450: / Line 487: @@
 We choose <math display="inline">C^\circ </math> continuous linear interpolations of the velocities, the pressure and the pressure gradient projections <math display="inline">\pi _i</math> over three node triangles (2D) and four node tetrahedra (3D). The linear interpolations are written as
+<span id="eq-28"></span>
 {| class="formulaSCP" style="width: 100%; text-align: left;"
 |-
@@ Line 460: / Line 498: @@
 |}
-where <math display="inline">n=3</math> (4) for triangles (tetrahedra), <math display="inline">\bar {(\cdot )}^j</math> denotes nodal variables and <math display="inline">N_j</math> are the linear shape functions <span id='citeF-1'></span>[[#cite-1|1]].
+where <math display="inline">n=3</math> (4) for triangles (tetrahedra), <math display="inline">\bar {(\cdot )}^j</math> denotes nodal variables and <math display="inline">N_j</math> are the linear shape functions <span id='citeF-1'></span>[[#cite-1|[1]]].
-Substituting the approximations (28) into Eqs.(25&#8211;27) and choosing the Galerking form with <math display="inline">\delta u_i =q=\delta \pi _i =N_i</math> leads to following system of discretized equations
+Substituting the approximations ([[#eq-28|28]]) into Eqs.([[#eq-25|25]]&#8211;[[#eq-27|27]]) and choosing the Galerking form with <math display="inline">\delta u_i =q=\delta \pi _i =N_i</math> leads to following system of discretized equations
+<span id="eq-29.a"></span>
 {| class="formulaSCP" style="width: 100%; text-align: left;"
 |-
@@ Line 471: / Line 510: @@
 | style="text-align: center;" | <math>\displaystyle {\boldsymbol M}\dot{\bar {\boldsymbol u}} + {\boldsymbol K} \bar {\boldsymbol u} - {\boldsymbol G}\bar {\boldsymbol p} ={\boldsymbol f}</math>
 |}
-| style="width: 5px;text-align: right;white-space: nowrap;" | (29a)
+| style="width: 5px;text-align: right;white-space: nowrap;" | (29.a)
 |}
+<span id="eq-29.b"></span>
 {| class="formulaSCP" style="width: 100%; text-align: left;"
 |-
@@ Line 481: / Line 521: @@
 | style="text-align: center;" | <math>\displaystyle {\boldsymbol G}^T \bar {\boldsymbol u} + {\boldsymbol L}\bar {\boldsymbol p}+{\boldsymbol Q}\bar {\boldsymbol \pi }={\boldsymbol 0}</math>
 |}
-| style="width: 5px;text-align: right;white-space: nowrap;" | (29b)
+| style="width: 5px;text-align: right;white-space: nowrap;" | (29.b)
 |}
+<span id="eq-29.c"></span>
 {| class="formulaSCP" style="width: 100%; text-align: left;"
 |-
@@ Line 491: / Line 532: @@
 | 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;" | (29c)
+| style="width: 5px;text-align: right;white-space: nowrap;" | (29.c)
 |}
 where the element contributions are given by (for 2D problems)
+<span id="eq-30"></span>
 {| class="formulaSCP" style="width: 100%; text-align: left;"
 |-
@@ Line 503: / Line 545: @@
 | style="text-align: center;" | <math>\begin{array}{l}\displaystyle M_{ij}= \int _{\Omega ^e} \rho N_iN_j d\Omega \quad ,\quad \displaystyle {\boldsymbol K}_{ij}= \int _{\Omega ^e}{\boldsymbol B}_i^T {\boldsymbol D}{\boldsymbol B}_j d\Omega \\ \\ \displaystyle {\boldsymbol G}_{ij}= \int _{\Omega ^e}({\boldsymbol \nabla } N_i)N_j d\Omega \quad ,\quad {\boldsymbol \nabla }= \left[{\partial  \over \partial x_1},{\partial  \over \partial x_2}\right]^T \\ \\ \displaystyle L_{ij}=\int _{\Omega ^e} {\boldsymbol \nabla }^T N_i [\tau ] {\boldsymbol \nabla } N_j d\Omega \quad ,\quad [\tau ]= \left[\begin{matrix}\tau _1 &0\\ 0 & \tau _2\\\end{matrix}\right]\\ \\ \displaystyle {\boldsymbol Q}= [{\boldsymbol Q}^1,{\boldsymbol Q}^2] \quad ,\quad  \displaystyle Q_{ij}^k = \int _{\Omega ^e}\tau _k {\partial N_i \over \partial x_k} N_jd\Omega \\ \\ \displaystyle \hat {\boldsymbol M}= \left[\begin{matrix}\hat {\boldsymbol M}^1 &{\boldsymbol 0}\\  {\boldsymbol 0} & \hat {\boldsymbol M}^2\\\end{matrix}\right]\quad ,\quad  \hat {\boldsymbol M}_{ij} = \int _{\Omega ^e} \tau _k N_i N_j d\Omega \\ \\ \displaystyle {\boldsymbol f}_i = \int _{\Omega ^e} N_i {\boldsymbol b}d\Omega + \int _{\Gamma ^e}N_i {\boldsymbol t} d\Gamma \quad ,\quad  {\boldsymbol b}=[b_1,b_2]^T \quad ,\quad {\boldsymbol t}=[t_1,t_2]^T  \end{array} </math>
 |}
+| style="width: 5px;text-align: right;white-space: nowrap;" | (30)
 |}
@@ Line 509: / Line 552: @@
 In above <math display="inline">{\boldsymbol B}_i</math> is the standard strain rate matrix and <math display="inline">{\boldsymbol D}</math> the deviatoric constitutive matrix (assuming <math display="inline">\displaystyle {\partial u_i \over \partial x_i}=0</math>). For 2D problems
+<span id="eq-31"></span>
 {| class="formulaSCP" style="width: 100%; text-align: left;"
 |-
@@ Line 516: / Line 560: @@
 | style="text-align: center;" | <math>{\boldsymbol B}_i =\left[\begin{matrix}\displaystyle {\partial N_i \over \partial x_1}&0\\ 0 & \displaystyle {\partial N_i \over \partial x_2}\\ \displaystyle {\partial N_i \over \partial x_2} & \displaystyle {\partial N_i \over \partial x_1}\\\end{matrix}\right]\quad ,\quad {\boldsymbol D} =\mu \left[\begin{matrix}2&0&0\\ 0&2&0\\ 0&0&1\\\end{matrix}\right] </math>
 |}
-| style="width: 5px;text-align: right;white-space: nowrap;" | (30)
+| style="width: 5px;text-align: right;white-space: nowrap;" | (31)
 |}
-The steady-state form of Eqs.(29) can be expressed in matrix form as
+The steady-state form of Eqs.([[#eq-29|29]]) can be expressed in matrix form as
+<span id="eq-32"></span>
 {| class="formulaSCP" style="width: 100%; text-align: left;"
 |-
@@ Line 528: / Line 573: @@
 | style="text-align: center;" | <math>\left[\begin{matrix}{\boldsymbol K}&-{\boldsymbol G}&{\boldsymbol 0}\\  -{\boldsymbol G}^T & - {\boldsymbol L} &-{\boldsymbol Q}\\ {\boldsymbol 0}& -{\boldsymbol Q}^T&-\hat {\boldsymbol M}\\\end{matrix}\right]\left\{\begin{matrix}\bar {\boldsymbol u}\\ \bar {\boldsymbol p}\\ \bar {\boldsymbol \pi }\\\end{matrix}\right\}= \left\{\begin{matrix}{\boldsymbol f}\\ {\boldsymbol 0}\\{\boldsymbol 0}\\ \end{matrix}\right\} </math>
 |}
-| style="width: 5px;text-align: right;white-space: nowrap;" | (31)
+| style="width: 5px;text-align: right;white-space: nowrap;" | (32)
 |}
-The system is symmetric and always positive definite and therefore leads to a non singular solution. We note that this property holds for ''any interpolation function'' chosen for <math display="inline">\bar {\boldsymbol u},\bar {\boldsymbol p}</math> and <math display="inline">\bar {\boldsymbol \pi }</math>, therefore overcoming the Babuŝka-Brezzi (BB) restrictions <span id='citeF-1'></span>[[#cite-1|1]].
+The system is symmetric and always positive definite and therefore leads to a non singular solution. We note that this property holds for ''any interpolation function'' chosen for <math display="inline">\bar {\boldsymbol u},\bar {\boldsymbol p}</math> and <math display="inline">\bar {\boldsymbol \pi }</math>, therefore overcoming the Babuŝka-Brezzi (BB) restrictions <span id='citeF-1'></span>[[#cite-1|[1]]].
 A reduced velocity-pressure formulation can be obtained by eliminating the pressure gradient projection variables <math display="inline">\bar {\boldsymbol \pi }</math> from the last equation to give
+<span id="eq-33"></span>
 {| class="formulaSCP" style="width: 100%; text-align: left;"
 |-
@@ Line 542: / Line 588: @@
 | style="text-align: center;" | <math>\left[\begin{matrix}{\boldsymbol K}&-{\boldsymbol G}\\ -{\boldsymbol G}^T & - ({\boldsymbol L} -{\boldsymbol Q}\hat {\boldsymbol M}^{-1}{\boldsymbol Q}^T)\\\end{matrix}\right]\left\{\begin{matrix}\bar {\boldsymbol u}\\ \bar {\boldsymbol p}\\\end{matrix}\right\}= \left\{\begin{matrix}{\boldsymbol f}\\ {\boldsymbol 0}\\\end{matrix}\right\} </math>
 |}
-| style="width: 5px;text-align: right;white-space: nowrap;" | (32)
+| style="width: 5px;text-align: right;white-space: nowrap;" | (33)
 |}
@@ Line 551: / Line 597: @@
 The solution process can be advanced in time in a (quasi-nearly) implicit iterative manner using the following scheme.
-===Step 1===
+<span id="step1"></span>
+'''Step 1'''
+<span id="eq-34"></span>
 {| class="formulaSCP" style="width: 100%; text-align: left;"
 |-
@@ Line 560: / Line 607: @@
 | style="text-align: center;" | <math>\bar{\boldsymbol u}^{n+1,i} = \bar{\boldsymbol u}^n - \Delta t {\boldsymbol M}^{-1} [{\boldsymbol K} \bar{\boldsymbol u}^{n+\theta _1,i-1}- {\boldsymbol G}\bar{\boldsymbol p}^{n+\theta _2,i-1} -{\boldsymbol f}^{n+1}] </math>
 |}
-| style="width: 5px;text-align: right;white-space: nowrap;" | (33)
+| style="width: 5px;text-align: right;white-space: nowrap;" | (34)
 |}
-===Step 2===
+<span id="step2"></span>
+'''Step 2'''
+<span id="eq-35"></span>
 {| class="formulaSCP" style="width: 100%; text-align: left;"
 |-
@@ Line 572: / Line 620: @@
 | style="text-align: center;" | <math>\bar{\boldsymbol p}^{n+1,i}=-{\boldsymbol L}^{-1} [{\boldsymbol G}^T \bar {\boldsymbol u}^{n+1,i}+{\boldsymbol Q}\bar {\boldsymbol \pi }^{n+\theta _3,i-1}] </math>
 |}
-| style="width: 5px;text-align: right;white-space: nowrap;" | (34)
+| style="width: 5px;text-align: right;white-space: nowrap;" | (35)
 |}
-===Step 3===
+<span id="step3"></span>
+'''Step 3'''
+<span id="eq-36"></span>
 {| class="formulaSCP" style="width: 100%; text-align: left;"
 |-
@@ Line 584: / Line 633: @@
 | style="text-align: center;" | <math>\bar {\boldsymbol \pi }^{n+1,i}= - \hat {\boldsymbol M}^{-1}{\boldsymbol Q}^T \bar {\boldsymbol p}^{n+1,i} </math>
 |}
-| style="width: 5px;text-align: right;white-space: nowrap;" | (35)
+| style="width: 5px;text-align: right;white-space: nowrap;" | (36)
 |}
 where <math display="inline">0\le \theta _i\le 1</math>.
-Steps 1&#8211;3 are repeated until converged values of the velocities, <math display="inline">u_i</math> the pressure <math display="inline">p</math> and the pressure gradient projection variables <math display="inline">\pi _i</math> is obtained. The process within a time step increment is completed with steps 4 and 5 described below.
+Steps [[#step1|1]]&#8211;[[#step3|3]] are repeated until converged values of the velocities, <math display="inline">u_i</math> the pressure <math display="inline">p</math> and the pressure gradient projection variables <math display="inline">\pi _i</math> is obtained. The process within a time step increment is completed with steps [[#step4|4]] and [[#step5|5]] described below.
-===Step 4===
+<span id="step4"></span>
+'''Step 4'''
 Update the mesh nodes in a lagrangian manner as
+<span id="eq-37"></span>
 {| class="formulaSCP" style="width: 100%; text-align: left;"
 |-
@@ Line 602: / Line 652: @@
 | style="text-align: center;" | <math>{\boldsymbol x}_i^{n+1} = {\boldsymbol x}_i^{n}+\bar {\boldsymbol u}_i^{n+1} \Delta t </math>
 |}
-| style="width: 5px;text-align: right;white-space: nowrap;" | (36)
+| style="width: 5px;text-align: right;white-space: nowrap;" | (37)
 |}
-===Step 5===
+<span id="step5"></span>
+'''Step 5'''
 Generate a new mesh. This can be effectively performed using the procedure described in Section 5. Indeed the mesh regeneration can take place after a prescribed number of time steps or when  the nodal displacements induce significant distorsions in some element shapes.
-Details of the treatment of the boundary conditions in the lagrangian flow formulation can be found in <span id='citeF-45'></span>[[#cite-45|45]].
+Details of the treatment of the boundary conditions in the lagrangian flow formulation can be found in <span id='citeF-51'></span>[[#cite-51|[51]],[[#cite-52|52]]].
-<br/>
+'''''Some practical remarks'''''
-<math display="inline">\boldsymbol {Some~practical~remarks}</math>
+In Eqs.([[#eq-34|34]]&#8211;[[#eq-36|36]]) <math display="inline">\bar{(\cdot )}^{n,i}</math> denote nodal values at the <math display="inline">n</math>th time step and the ith iteration. Note that <math display="inline">{\boldsymbol A}^{n+\theta _1,i-1}\equiv {\boldsymbol A}(\bar {\boldsymbol u}^{n+\theta _1,i-1})</math> etc. Also <math display="inline">(\cdot )^{n+\theta _i,0}\equiv (\cdot )^n</math> for the computations in step [[#step1|1]] at the onset of the iterations.
-In Eqs.(34&#8211;36) <math display="inline">\bar{(\cdot )}^{n,i}</math> denote nodal values at the <math display="inline">n</math>th time step and the ith iteration. Note that <math display="inline">{\boldsymbol A}^{n+\theta _1,i-1}\equiv {\boldsymbol A}(\bar {\boldsymbol u}^{n+\theta _1,i-1})</math> etc. Also <math display="inline">(\cdot )^{n+\theta _i,0}\equiv (\cdot )^n</math> for the computations in step 1 at the onset of the iterations.
+Steps [[#step1|1]] and [[#step3|3]] can be solved explicitely by choosing a ''lumped (diagonal) form'' of matrices <math display="inline">M</math> and <math display="inline">\hat {\boldsymbol M}</math>. In this manner the main computational cost is the solution of step [[#step2|2]] involving the inverse of a Laplacian matrix. This can be solved very effectively using an iterative method such as the conjugate gradient method or similar.
-Steps 1 and 3 can be solved explicitely by choosing a ''lumped (diagonal) form'' of matrices '''M''' and <math display="inline">\hat {\boldsymbol M}</math>. In this manner the main computational cost is the solution of step 2 involving the inverse of a Laplacian matrix. This can be solved very effectively using an iterative method such as the conjugate gradient method or similar.
+For <math display="inline">\theta _i\not =0</math> the iterative proces is unavoidable. The iterations follow until convergence is reached in an adequate error norm in terms of the velocity and pressure variables, or the residuals <math display="inline">r_{m_i}</math> and <math display="inline">r_d</math>. Indeed some ot the <math display="inline">\theta _i</math>'s in Eqs.([[#eq-34|34]])&#8211;([[#eq-36|36]]) can be made equal to zero. Note that for <math display="inline">\theta _2 =0</math> the algorithm is unconditionally unstable. A particularity interesting and simple semi-implicit form is obtained by making <math display="inline">\theta _1 = \theta _2=\theta _3=0</math>. Now all steps can be solved explicitely with exception of Step [[#step2|2]] for the pressure, which still requires the solution of a simultaneous system of equations.
-For <math display="inline">\theta _i\not =0</math> the iterative proces is unavoidable. The iterations follow until convergence is reached in an adequate error norm in terms of the velocity and pressure variables, or the residuals <math display="inline">r_{m_i}</math> and <math display="inline">r_d</math>. Indeed some ot the <math display="inline">\theta _i</math>'s in Eqs.(34)&#8211;(36) can be made equal to zero. Note that for <math display="inline">\theta _2 =0</math> the algorithm is unconditionally unstable. A particularity interesting and simple semi-implicit form is obtained by making <math display="inline">\theta _1 = \theta _2=\theta _3=0</math>. Now all steps can be solved explicitely with exception of Step 2 for the pressure, which still requires the solution of a simultaneous system of equations.
 Convergence of this solution scheme is however difficult for some problems. An enhanced version of the algorithm can be obtained by simply adding the term <math display="inline">\hat {\boldsymbol  L} (\bar {\boldsymbol p}^{n+1,i} - \bar {\boldsymbol p}^{n+1,i-1})</math> where <math display="inline">\hat L_{ij} =\Delta t \int _{\Omega ^e} {\boldsymbol \nabla }^T N_i {\boldsymbol \nabla } N_j d\Omega </math> to the equation for the computation of the pressure in the second step. The new term acts as a preconditioner of the pressure equation given now by
+<span id="eq-38"></span>
 {| class="formulaSCP" style="width: 100%; text-align: left;"
 |-
@@ Line 630: / Line 680: @@
 | style="text-align: center;" | <math>\bar {\boldsymbol p}^{n+1,i}=-[{\boldsymbol L} + \hat {\boldsymbol L}]^{-1}[{\boldsymbol G}^T \bar {\boldsymbol u}^{n+1,i}+\hat {\boldsymbol L}\bar {\boldsymbol p}^{n+1,i-1} +{\boldsymbol Q}\bar {\boldsymbol \pi }^{n+\theta _3,i-1}] </math>
 |}
-| style="width: 5px;text-align: right;white-space: nowrap;" | (37)
+| style="width: 5px;text-align: right;white-space: nowrap;" | (38)
 |}
@@ Line 641: / Line 691: @@
 An effective algorithm can be obtained by splitting the pressure from the momentum equations as follows
+<span id="eq-39.a"></span>
 {| class="formulaSCP" style="width: 100%; text-align: left;"
 |-
@@ Line 648: / Line 699: @@
 | style="text-align: center;" | <math>\bar {\boldsymbol u}^* = \bar {\boldsymbol u}^n -\Delta t {\boldsymbol M}^{-1}[{\boldsymbol K} + \bar{\boldsymbol u}^{n+\theta _1}-\alpha {\boldsymbol G}\bar{\boldsymbol p}^{n} -{\boldsymbol f}^{n+1}]</math>
 |}
-| style="width: 5px;text-align: right;white-space: nowrap;" | (39a)
+| style="width: 5px;text-align: right;white-space: nowrap;" | (39.a)
 |}
+<span id="eq-39.b"></span>
 {| class="formulaSCP" style="width: 100%; text-align: left;"
 |-
@@ Line 658: / Line 710: @@
 | style="text-align: center;" | <math>\bar{\boldsymbol u}^{n+1}= \bar{\boldsymbol u}^* + \Delta t {\boldsymbol M}^{-1}{\boldsymbol G}\delta \bar{\boldsymbol p}</math>
 |}
-| style="width: 5px;text-align: right;white-space: nowrap;" | (39b)
+| style="width: 5px;text-align: right;white-space: nowrap;" | (39.b)
 |}
-In Eq.(39a) <math display="inline">\alpha </math> is a variable taking values equal to zero or one. For <math display="inline">\alpha =0</math>, <math display="inline">\delta p \equiv p^{n+1}</math> and for <math display="inline">\alpha =1</math>, <math display="inline">\delta p =\Delta p</math>. Note that in both cases  the sum of Eqs.(39a) and (39b) gives the time discretization of the momentum equations with the pressures computed at <math display="inline">t^{n+1}</math>. The value of <math display="inline">\bar {\boldsymbol u}^{n+1}</math> from Eq.(39b) is substituted now into Eq.(29b) to give
+In Eq.([[#eq-39.a|39.a]]) <math display="inline">\alpha </math> is a variable taking values equal to zero or one. For <math display="inline">\alpha =0</math>, <math display="inline">\delta p \equiv p^{n+1}</math> and for <math display="inline">\alpha =1</math>, <math display="inline">\delta p =\Delta p</math>. Note that in both cases  the sum of Eqs.([[#eq-39.a|39.a]]) and ([[#eq-39.b|39.b]]) gives the time discretization of the momentum equations with the pressures computed at <math display="inline">t^{n+1}</math>. The value of <math display="inline">\bar {\boldsymbol u}^{n+1}</math> from Eq.([[#eq-39.b|39.b]]) is substituted now into Eq.([[#eq-29.b|29.b]]) to give
+<span id="eq-40.a"></span>
 {| class="formulaSCP" style="width: 100%; text-align: left;"
 |-
@@ Line 670: / Line 723: @@
 | style="text-align: center;" | <math>{\boldsymbol G}^T\bar {\boldsymbol u}^* + \Delta t {\boldsymbol G}^T {\boldsymbol M}^{-1}{\boldsymbol G}\delta \bar{\boldsymbol p} + {\boldsymbol L} \bar{\boldsymbol p}^{n+1}+ {\boldsymbol Q}\bar {\boldsymbol \pi }^{n+\theta _3}=0</math>
 |}
-| style="width: 5px;text-align: right;white-space: nowrap;" | (40a)
+| style="width: 5px;text-align: right;white-space: nowrap;" | (40.a)
 |}
 The product <math display="inline">{\boldsymbol G}^T {\boldsymbol M}^{-1}{\boldsymbol G}</math> can be approximated by a laplacian matrix, i.e.
+<span id="eq-40.b"></span>
 {| class="formulaSCP" style="width: 100%; text-align: left;"
 |-
@@ Line 682: / Line 736: @@
 | style="text-align: center;" | <math>{\boldsymbol G}^T {\boldsymbol M}^{-1}{\boldsymbol G}=\hat {\boldsymbol L} \quad \hbox{with } \hat{\boldsymbol L} \simeq \int _{\Omega ^e} {1\over \rho } {\boldsymbol \nabla }^T N_i {\boldsymbol \nabla } N_j \,d\Omega  </math>
 |}
-| style="width: 5px;text-align: right;white-space: nowrap;" | (40b)
+| style="width: 5px;text-align: right;white-space: nowrap;" | (40.b)
 |}
-A semi-implicit algorithm can thus be derived as follows.<br/>
+A semi-implicit algorithm can thus be derived as follows.
-'''Step 1''' Compute <math display="inline">{\boldsymbol u}^*</math> explicitely from Eq.(39a) with <math display="inline">{\boldsymbol M}={\boldsymbol M}_d</math> where subscript <math display="inline">d</math> denotes hereonwards a diagonal matrix.
+<span id="step-1"></span>
+'''Step 1''' Compute <math display="inline">{\boldsymbol u}^*</math> explicitely from Eq.([[#eq-39.a|39.a]]) with <math display="inline">{\boldsymbol M}={\boldsymbol M}_d</math> where subscript <math display="inline">d</math> denotes hereonwards a diagonal matrix.
-<br/>
+<span id="step-2"></span>
+'''Step 2''' Compute <math display="inline">\delta \bar {\boldsymbol p}</math> and <math display="inline">{\boldsymbol p}^{n+1}</math> from Eq.([[#eq-40.a|40.a]]) as
-'''Step 2''' Compute <math display="inline">\delta \bar {\boldsymbol p}</math> and <math display="inline">{\boldsymbol p}^{n+1}</math> from Eq.(40a) as
+<span id="eq-41.a"></span>
 {| class="formulaSCP" style="width: 100%; text-align: left;"
 |-
@@ Line 700: / Line 755: @@
 | style="text-align: center;" | <math>\delta \bar {\boldsymbol p} =-({\boldsymbol L}+\Delta t \hat  {\boldsymbol L})^{-1} [{\boldsymbol G}^T\bar{\boldsymbol u}^* +{\boldsymbol Q}\bar {\boldsymbol \pi }^{n+\theta _4}+ \alpha {\boldsymbol L} \bar {\boldsymbol p}^n]</math>
 |}
-| style="width: 5px;text-align: right;white-space: nowrap;" | (41a)
+| style="width: 5px;text-align: right;white-space: nowrap;" | (41.a)
 |}
+<span id="eq-41.b"></span>
 {| class="formulaSCP" style="width: 100%; text-align: left;"
 |-
@@ Line 710: / Line 766: @@
 | style="text-align: center;" | <math>{\boldsymbol p}^{n+1} = {\boldsymbol p}^n + \delta {\boldsymbol p}</math>
 |}
-| style="width: 5px;text-align: right;white-space: nowrap;" | (41b)
+| style="width: 5px;text-align: right;white-space: nowrap;" | (41.b)
 |}
-'''Step 3''' Compute <math display="inline"> \bar{\boldsymbol u}^{n+1}</math> explicitely from Eq.(39b) with <math display="inline">{\boldsymbol M}={\boldsymbol M}_d</math>
+<span id="step-3"></span>
+'''Step 3''' Compute <math display="inline"> \bar{\boldsymbol u}^{n+1}</math> explicitely from Eq.([[#eq-39.b|39.b]]) with <math display="inline">{\boldsymbol M}={\boldsymbol M}_d</math>
-'''Step 4''' Compute <math display="inline">  \bar{\boldsymbol \pi }^{n+1}</math> explicitely from Eq.(36) as
+<span id="step-4"></span>
+'''Step 4''' Compute <math display="inline">  \bar{\boldsymbol \pi }^{n+1}</math> explicitely from Eq.([[#eq-36|36]]) as
+<span id="eq-42"></span>
 {| class="formulaSCP" style="width: 100%; text-align: left;"
 |-
@@ Line 727: / Line 786: @@
 |}
-Steps 5 and 6 coincide with steps 4 and 5 of Section 4.1.
+Steps 5 and 6 coincide with steps [[#step4|4]] and [[#step5|5]] of Section [[#4.1 Quasi-implicit scheme|4.1.]]
-This algorithm has an additional step than  the iterative scheme of Section 4.1. The advantage is that now Steps 1 and 2 can be fully linearized by choosing <math display="inline">\theta _1 = \theta _2 =\theta _3 =0</math>. Also the equation for the pressure variables in Step 2 has improved stabilization properties due to the additional laplacian matrix <math display="inline">\hat {\boldsymbol L}</math>.
+This algorithm has an additional step than  the iterative scheme of Section [[#4.1 Quasi-implicit scheme|4.1.]] The advantage is that now Steps [[#step-1|1]] and [[#step-2|2]] can be fully linearized by choosing <math display="inline">\theta _1 = \theta _2 =\theta _3 =0</math>. Also the equation for the pressure variables in Step [[#step-2|2]] has improved stabilization properties due to the additional laplacian matrix <math display="inline">\hat {\boldsymbol L}</math>.
-The boundary conditions are applied as follows. No condition is applied in the computation of the fractional velocities <math display="inline">{\boldsymbol u}^*</math> in Eq.(39a). The prescribed velocities at the boundary are applied when solving for <math display="inline">\bar{\boldsymbol u}^{n+1}</math> in step 3. The prescribed pressures at the boundary are imposed by making zero the pressure increments at the relevant boundary nodes and making <math display="inline">\bar{\boldsymbol p}^n</math> equal to the prescribed pressure values.
+The boundary conditions are applied as follows. No condition is applied in the computation of the fractional velocities <math display="inline">{\boldsymbol u}^*</math> in Eq.([[#eq-39.a|39.a]]). The prescribed velocities at the boundary are applied when solving for <math display="inline">\bar{\boldsymbol u}^{n+1}</math> in step [[#step-3|3]]. The prescribed pressures at the boundary are imposed by making zero the pressure increments at the relevant boundary nodes and making <math display="inline">\bar{\boldsymbol p}^n</math> equal to the prescribed pressure values.
 ==5 GENERATION OF A NEW MESH==
-One of the key points for success of the lagrangian flow formulation here described is the fast regeneration of a mesh at every time step on the basis of the position of the nodes in the space domain. In our work the mesh is generated using the so called extended Delaunay tesselation (EDT) presented in [44,47]. The EDT allows to generate meshes of elements with arbitrary polyhedrical shapes (combining 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>, being <math display="inline">n</math> the total number of nodes in the mesh. The shape functions of the elements can be simply obtained using the so called non-sibsonian interpolations. Details of the mesh generation procedure can be found in [44&#8211;47].
+One of the key points for success of the lagrangian flow formulation here described is the fast regeneration of a mesh at every time step on the basis of the position of the nodes in the space domain. In our work the mesh is generated using the so called extended Delaunay tesselation (EDT) presented in <span id='citeF-44'></span><span id='citeF-47'></span>[[#cite-44|[44]],[[#cite-47|47]]]. The EDT allows to generate meshes of elements with arbitrary polyhedrical shapes (combining 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>, being <math display="inline">n</math> the total number of nodes in the mesh. The shape functions of the elements can be simply obtained using the so called non-sibsonian interpolations. Details of the mesh generation procedure can be found in <span id='citeF-44'></span><span id='citeF-47'></span>[[#cite-44|[44]]-[[#cite-47|47]]].
 Once the new mesh has been generated at each time step the numerical solution is found using the finite element algorithm described in the paper.
@@ Line 749: / Line 808: @@
 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  bigger than, are considered as boundary nodes'' (See Figure 2).
+''All nodes on an empty sphere with a radius  bigger than, are considered as boundary nodes'' (See Figure [[#img-2|2]]).
-Thus, <math display="inline">\alpha </math>  is a parameter close to, but greater than one. Note that this criterion is coincident with the Alpha Shape concept <span id='citeF-47'></span>[[#cite-47|47]].
+Thus, <math display="inline">\alpha </math>  is a parameter close to, but greater than one. Note that this criterion is coincident with the Alpha Shape concept <span id='citeF-47'></span>[[#cite-47|[47]]].
 <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_Test_756558304-Fig2.png|351px|Contour recognition:  Empty circles with radius αh(x)  define the boundary particles.]]
+|style="padding: 10px;"|[[Image:Draft_Samper_187181919-Fig2.png|200px|Contour recognition:  Empty circles with radius αh(x)  define the boundary particles.]]
 |- style="text-align: center; font-size: 75%;"
-| colspan="1" | '''Figure 2:''' Contour recognition:  Empty circles with radius <math>\alpha h(x)</math>  define the boundary particles.
+| colspan="1" style="padding: 10px;"|  '''Figure 2:''' Contour recognition:  Empty circles with radius <math>\alpha h(x)</math>  define the boundary particles.
 |}
@@ Line 769: / Line 828: @@
 ==7 COMPUTATION OF THE CHARACTERISTIC LENGTHS==
-The evaluation of the stabilization parameters is one of the crucial issues in stabilized methods. Most of existing methods use expressions which are direct extensions of the values obtained for the simplest 1D case. In this work we will accept the so called SUPG assumption, i.e. to admit that vector <math display="inline">{\boldsymbol h}</math> has the direction of the velocity field <span id='citeF-32'></span>[[#cite-32|32]]. This gives
+The evaluation of the stabilization parameters is one of the crucial issues in stabilized methods. Most of existing methods use expressions which are direct extensions of the values obtained for the simplest 1D case. In this work we will accept the so called SUPG assumption, i.e. to admit that vector <math display="inline">{\boldsymbol h}</math> has the direction of the velocity field <span id='citeF-32'></span><span id='citeF-37'></span>[[#cite-32|[32]]-[[#cite-37|37]]]. This gives
+<span id="eq-43"></span>
 {| class="formulaSCP" style="width: 100%; text-align: left;"
 |-
@@ Line 783: / Line 843: @@
 where <math display="inline">u=\vert {\boldsymbol u}\vert </math> and <math display="inline">h_s</math> as the “streamline” characteristic length
+<span id="eq-44"></span>
 {| class="formulaSCP" style="width: 100%; text-align: left;"
 |-
@@ Line 788: / Line 849: @@
 {| style="text-align: left; margin:auto;width: 100%;"
 |-
-| style="text-align: center;" | <math>h_s=\max ({\boldsymbol l}^T_j {\boldsymbol u})/{u} \quad , \quad  j=1,n_s </math>
+| style="text-align: center;" | <math>h_s=\max \frac{({\boldsymbol l}^T_j {\boldsymbol u})}{u} \quad , \quad  j=1,n_s </math>
 |}
 | style="width: 5px;text-align: right;white-space: nowrap;" | (44)
@@ Line 795: / Line 856: @@
 where <math display="inline">{\boldsymbol l}_j</math> are the vectors defining the element sides (<math display="inline">n_s=6</math> for tetrahedra).
-A more consistent evaluation of '''h''' based on a diminishing residual technique can be found in <span id='citeF-37'></span>[[#cite-37|37]].
+A more consistent evaluation of <math display="inline">h</math> based on a diminishing residual technique can be found in <span id='citeF-37'></span>[[#cite-37|[37]]].
 ==8 EXAMPLES==
-All examples presented here have been obtained withthe fractional step algorithm described in Section 4.2.
+All examples presented here have been obtained withthe fractional step algorithm described in Section [[#4.2 Fractional step method|4.2.]]
 ===8.1 Collapse of a water column===
-The first problem solved with the lagrangian formulation is the study of the collapse of a water column. This problem was solved by Koshizu and Oka <span id='citeF-43'></span>[[#cite-43|43]] both experimentally and numerically. It has became a classical example to test the validation of the lagrangian formulation in fluid flows. The water is initially located on the left supported by a removable board. The collapse starts at time t = 0, when the removable board is removed. Viscosity and surface tension are neglected. Figure&nbsp;3 shows the point positions at different time steps. The dark points represent the free-surface detected with the algorithms described in Section 5. The internal points are gray and the fixed points are black.
+The first problem solved with the lagrangian formulation is the study of the collapse of a water column. This problem was solved by Koshizu and Oka <span id='citeF-43'></span>[[#cite-43|[43]]] both experimentally and numerically. It has became a classical example to test the validation of the lagrangian formulation in fluid flows. The water is initially located on the left supported by a removable board. The collapse starts at time <math display="inline">t = 0</math>, when the removable board is removed. Viscosity and surface tension are neglected. Figure&nbsp;[[#img-3|3]] shows the point positions at different time steps. The dark points represent the free-surface detected with the algorithms described in Section [[#5 GENERATION OF A NEW MESH|5]]. The internal points are gray and the fixed points are black.
-The water is running on the bottom wall until, near 0.3 sec, it impinges on the right vertical wall. Breaking waves appear at 0.6 sec. Around 1 sec. the water reaches the left wall. Agreement with the experimental results of <span id='citeF-43'></span>[[#cite-43|43]] both in the shape of the free surface as well as in the time evolution are excellent.
+The water is running on the bottom wall until, near <math display="inline">0.3</math> sec, it impinges on the right vertical wall. Breaking waves appear at <math display="inline">0.6</math> sec. Around <math display="inline">1</math> sec. the water reaches the left wall. Agreement with the experimental results of <span id='citeF-43'></span>[[#cite-43|[43]]] both in the shape of the free surface as well as in the time evolution are excellent.
 <div id='img-3'></div>
 {| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
 |-
-|[[Image:Draft_Samper_800412471_3321_img-3.JPG|600px|Water column collapse at different time steps.]]
+|style="padding: 10px;"| [[File:Draft_Samper_187181919_2382_img-3.JPG|400px|Water column collapse at different time steps.]]
-|- style="text-align: center; font-size: 75%;"
-| colspan="1" | '''Figure 3:''' Water column collapse at different time steps.
-|}
-<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_800412471_9131_img-3b.JPG|600px|cont.]]
+|[[File:Draft_Samper_187181919_6277_img-3b.JPG|400px|cont.]]
 |- style="text-align: center; font-size: 75%;"
-| colspan="1" | '''Figure 4:''' cont.
+| colspan="1" style="padding: 10px;"|  '''Figure 3:''' Water column collapse at different time steps.
 |}
-Figure 3 shows the point positions at different time steps. The darker points represent the free-surface detected.
+Figure [[#img-3|3]] shows the point positions at different time steps. The darker points represent the free-surface detected.
-The water is running on the bottom wall until, near 0.3 sec, it impinges on the right vertical wall. Breaking waves appear at 0.6 sec. Around <math display="inline">t=1</math> sec. the main water wave reaches the left wall again Agreement with the experimental results of ref. <span id='citeF-43'></span>[[#cite-43|43]] both in the shape of the free surface as well as in the time development are excellent.
+The water is running on the bottom wall until, near <math display="inline">0.3</math> sec, it impinges on the right vertical wall. Breaking waves appear at <math display="inline">0.6</math> sec. Around <math display="inline">t=1</math> sec. the main water wave reaches the left wall again Agreement with the experimental results of ref. <span id='citeF-43'></span>[[#cite-43|[43]]] both in the shape of the free surface as well as in the time development are excellent.
-The 3D solution of the same problem is shown in Figure 4.
+The 3D solution of the same problem is shown in Figure [[#img-4|4]].
-<div id='img-5'></div>
+<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_Test_756558304-Fig7_1.png|600px|]]
+|style="padding: 10px;"| [[Image:Draft_Samper_187181919-Fig7_1.png|400px|]]
-|[[Image:Draft_Test_756558304-Fig7_2.png|600px|]]
 |-
-| colspan="2"|[[Image:Draft_Test_756558304-Fig7_3.png|600px|Water column collapse in a 3D domain. Different time step.]]
+|[[Image:Draft_Samper_187181919-Fig7_2.png|400px|]]
+|-
+|[[Image:Draft_Samper_187181919-Fig7_3.png|400px|Water column collapse in a 3D domain. Different time step.]]
 |- style="text-align: center; font-size: 75%;"
-| colspan="2" | '''Figure 5:''' Water column collapse in a 3D domain. Different time step.
+| colspan="2" style="padding: 10px;"| '''Figure 4:''' Dam collapse in a 3D domain. Different time step.
 |}
 ===8.2 Semi-submerged rotating water mill===
-The second example is the analysis of a rotating water mill semi submerged in water. A schematic representation of a water mill is presented in Figure 5. The blades of the mill have an imposed rotating velocity, while the water is initially in a stationary and flat position. Fluid structure interactions with free-surfaces and water fragmentation are well reproduced in this example.
+The second example is the analysis of a rotating water mill semi submerged in water. A schematic representation of a water mill is presented in Figure [[#img-5|5]]. The blades of the mill have an imposed rotating velocity, while the water is initially in a stationary and flat position. Fluid structure interactions with free-surfaces and water fragmentation are well reproduced in this example.
-<div id='img-6'></div>
+<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_800412471_6628_img-5.JPG|600px|Rotating water mill.]]
+|style="padding: 10px;"| [[File:Draft_Samper_187181919_4263_img-5.JPG|600px|Rotating water mill.]]
 |- style="text-align: center; font-size: 75%;"
-| colspan="1" | '''Figure 6:''' Rotating water mill.
+| colspan="1" style="padding: 10px;"| '''Figure 5:''' Rotating water mill.
 |}
 ===8.3 Sloshing problems===
-The simple problem of the free oscillation of an incompressible liquid in a container is considered next. Numerical solutions for this problem can be found in several references <span id='citeF-44'></span>[[#cite-44|44]]. This problem is interesting because there is an analytical solution for small amplitudes. Figure 6 shows a schematic of the problem, and represents also the point distribution in the initial position. The dark points represent the fixed points in which the velocity is fixed to zero.
+The simple problem of the free oscillation of an incompressible liquid in a container is considered next. Numerical solutions for this problem can be found in several references <span id='citeF-44'></span>[[#cite-44|[44]]]. This problem is interesting because there is an analytical solution for small amplitudes. Figure [[#img-6|6]] shows a schematic of the problem, and represents also the point distribution in the initial position. The dark points represent the fixed points in which the velocity is fixed to zero.
-<div id='img-7'></div>
+<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_Test_756558304-Fig3.png|351px|Sloshing. Initial point distribution.]]
+|style="padding: 10px;"| [[Image:Draft_Samper_187181919-Fig3.png|252px|Sloshing. Initial point distribution.]]
 |- style="text-align: center; font-size: 75%;"
-| colspan="1" | '''Figure 7:''' Sloshing. Initial point distribution.
+| colspan="1" style="padding: 10px;"| '''Figure 6:''' Sloshing. Initial point distribution.
 |}
-<div id='img-8'></div>
+<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_Test_756558304-Fig4.png|351px|Sloshing: Comparison of the numerical and analytical solution.]]
+|style="padding: 10px;"| [[Image:Draft_Samper_187181919-Fig4.png|351px|Sloshing: Comparison of the numerical and analytical solution.]]
 |- style="text-align: center; font-size: 75%;"
-| colspan="1" | '''Figure 8:''' Sloshing: Comparison of the numerical and analytical solution.
+| colspan="1" style="padding: 10px;"| '''Figure 7:''' Sloshing: Comparison of the numerical and analytical solution.
 |}
-Figure 7 shows the variation in time of the amplitude compared with the analytical results for the near in viscid case.  Little numerical viscosity is observed on the phase wave and amplitude in spite of the relative poor point distribution.
+Figure [[#img-7|7]] shows the variation in time of the amplitude compared with the analytical results for the near in viscid case.  Little numerical viscosity is observed on the phase wave and amplitude in spite of the relative poor point distribution.
-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 can be separated from the fluid domain due to their large velocity.  Figure 8 shows the numerical results obtained with the method presented in this paper for larger sloshing amplitudes. Breaking waves as well as separation effects can be seen on the free-surface. This particular and very complicated effect is apparently well represented by this model.
+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 can be separated from the fluid domain due to their large velocity.  Figure [[#img-8|8]] shows the numerical results obtained with the method presented in this paper for larger sloshing amplitudes. Breaking waves as well as separation effects can be seen on the free-surface. This particular and very complicated effect is apparently well represented by this model.
-<div id='img-9'></div>
+<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_Test_756558304-Fig5.png|351px|Sloshing: Different time step for large amplitudes.]]
+|style="padding: 10px;"| [[Image:Draft_Samper_187181919-Fig5.png|351px|Sloshing: Different time step for large amplitudes.]]
 |- style="text-align: center; font-size: 75%;"
-| colspan="1" | '''Figure 9:''' Sloshing: Different time step for large amplitudes.
+| colspan="1" style="padding: 10px;"| '''Figure 8:''' Sloshing: Different time step for large amplitudes.
 |}
-<div id='img-10'></div>
+<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_Test_756558304-Fig6.png|351px|Sloshing: Different time step for 3D domains.]]
+|style="padding: 10px;"| [[Image:Draft_Samper_187181919-Fig6.png|351px|Sloshing: Different time step for 3D domains.]]
 |- style="text-align: center; font-size: 75%;"
-| colspan="1" | '''Figure 10:''' Sloshing: Different time step for 3D domains.
+| colspan="1" style="padding: 10px;"| '''Figure 9:''' Sloshing: Different time step for 3D domains.
 |}
-In order to test the potentiality of the method in a 3D domain, the same sloshing problem was solved as a 3D problem.  Figure 9 show the different point position 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 sphere representation is only used in order to improve the visualization of the numerical results.
+In order to test the potentiality of the method in a 3D domain, the same sloshing problem was solved as a 3D problem.  Figure [[#img-9|9]] show the different point position 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 sphere representation is only used in order to improve the visualization of the numerical results.
 ===8.4 Wave breaking on a beach===
-A simulation of the propagation of a water wave and its breaking due to shoaling over a plane slope is presented next. This example was numerically studied in <span id='citeF-44'></span>[[#cite-44|44]] with a lagrangian formulation using directly the standard Finite Element Method with remeshing. There is also an analytical solution for a simplified approximation that is used for  comparisons <span id='citeF-45'></span>[[#cite-45|45]]. Figure 10 shows the initial point distribution and Figure 11 a comparison with the analytical free-surface at different time step. The geometry of the problem as well as a discussion of the analytical solution may be found in <span id='citeF-44'></span>[[#cite-44|44]].
+A simulation of the propagation of a water wave and its breaking due to shoaling over a plane slope is presented next. This example was numerically studied in <span id='citeF-44'></span>[[#cite-44|[44]]] with a lagrangian formulation using directly the standard Finite Element Method with remeshing. There is also an analytical solution for a simplified approximation that is used for  comparisons <span id='citeF-45'></span>[[#cite-45|[45]]]. Figure [[#img-10|10]] shows the initial point distribution and Figure [[#img-11|11]] a comparison with the analytical free-surface at different time step. The geometry of the problem as well as a discussion of the analytical solution may be found in <span id='citeF-44'></span>[[#cite-44|[44]]].
-<div id='img-11'></div>
+<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_Test_756558304-Fig8.png|351px|Wave breaking on a beach. Initial geometry and point positions.]]
+|style="padding: 10px;"| [[Image:Draft_Samper_187181919-Fig8.png|351px|Wave breaking on a beach. Initial geometry and point positions.]]
 |- style="text-align: center; font-size: 75%;"
-| colspan="1" | '''Figure 11:''' Wave breaking on a beach. Initial geometry and point positions.
+| colspan="1" style="padding: 10px;"| '''Figure 10:''' Wave breaking on a beach. Initial geometry and point positions.
 |}
-Initially (Figures 11a. and 11b.) the wave travels over a constant depth bottom towards the slope with no ostensible change of shape. Strongly non-linear effects appear when the wave hits the slope (Fig. 11c.). The crest of the wave accelerates until it reaches the shore (Fig. 11d.). At this time the comparisons with the analytical solution are in agreement only in the wave position. The shape of the wave obtained with the numerical solution is totally different. The reason is that the analytical solution gives symmetrical shape waves, which are not physical, before the breaking process. Subsequently, a water jet is formed at the crest plunge making the breaking wave (Figs. 11e. and 11f.)  and coming in contact with the nearly still surface of the water ahead. In Ref. <span id='citeF-44'></span>[[#cite-44|44]] the evaluation is stopped before this contact point. Using the methodology proposed in this paper, the analysis may be continued until the  end.  In Figs. 11g. and 11h. the wave finally hits a lateral wall (introduced in the model to stop the lateral effects) producing drop separations, and then coming back toward to the left as a new wave.
+Initially (Figures [[#img-11|11.a.]] and [[#img-11|11.b.]]) the wave travels over a constant depth bottom towards the slope with no ostensible change of shape. Strongly non-linear effects appear when the wave hits the slope (Fig. [[#img-11|11.c.]]). The crest of the wave accelerates until it reaches the shore (Fig. [[#img-11|11.d.]]). At this time the comparisons with the analytical solution are in agreement only in the wave position. The shape of the wave obtained with the numerical solution is totally different. The reason is that the analytical solution gives symmetrical shape waves, which are not physical, before the breaking process. Subsequently, a water jet is formed at the crest plunge making the breaking wave (Figs. [[#img-11|11.e.]] and [[#img-11|11.f.]])  and coming in contact with the nearly still surface of the water ahead. In Ref. <span id='citeF-44'></span>[[#cite-44|44]] the evaluation is stopped before this contact point. Using the methodology proposed in this paper, the analysis may be continued until the  end.  In Figs. [[#img-11|11.g.]] and [[#img-11|11.h.]] the wave finally hits a lateral wall (introduced in the model to stop the lateral effects) producing drop separations, and then coming back toward to the left as a new wave.
-<div id='img-12'></div>
+<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_Test_756558304-Fig9.png|351px|Wave breaking on a beach. Comparison with analytical results at different time steps. Top: Numerical solution. Bottom: Analytical solution.]]
+|style="padding: 10px;"| [[Image:Draft_Samper_187181919-Fig9.png|454px|Wave breaking on a beach. Comparison with analytical results at different time steps. Top: Numerical solution. Bottom: Analytical solution.]]
 |- style="text-align: center; font-size: 75%;"
-| colspan="1" | '''Figure 12:''' Wave breaking on a beach. Comparison with analytical results at different time steps. Top: Numerical solution. Bottom: Analytical solution.
+| colspan="1" style="padding: 10px;"| '''Figure 11:''' Wave breaking on a beach. Comparison with analytical results at different time steps. Top: Numerical solution. Bottom: Analytical solution.
 |}
@@ Line 920: / Line 976: @@
 Nevertheless, a 2D domain is an easy case, and may be solved acceptably with any mesh generator. The true problems appears in a 3D domain, where the mesh generation is complicated with presence of slivers an other geometric mesh generation problems. In order to show the power of the tool presented, the same problem was solved in a 3D domain.
-To transform the wave breaking described before in a true 3D problem, the initial position of the wave was introduced having an oblique angle with the beach line. In this way, a 3D effect appears. When the wave  hits the slope, the crest of the wave accelerates differently accordingly with the depth, inducing the wave to correct its oblique position and break parallel to the beach. The results may be seen in Figure 12 for different time steps.
+To transform the wave breaking described before in a true 3D problem, the initial position of the wave was introduced having an oblique angle with the beach line. In this way, a 3D effect appears. When the wave  hits the slope, the crest of the wave accelerates differently accordingly with the depth, inducing the wave to correct its oblique position and break parallel to the beach. The results may be seen in Figure [[#img-12|12]] for different time steps.
-<div id='img-13'></div>
+<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_Test_756558304-Fig10.png|351px|Breaking wave on a beach: Oblique wave on a 3D domain.]]
+| [[Image:Draft_Samper_187181919-Fig10.png|351px|Breaking wave on a beach: Oblique wave on a 3D domain.]]
 |- style="text-align: center; font-size: 75%;"
-| colspan="1" | '''Figure 13:''' Breaking wave on a beach: Oblique wave on a 3D domain.
+| colspan="1" style="padding: 10px;"| '''Figure 12:''' Breaking wave on a beach: Oblique wave on a 3D domain.
 |}
 ===8.5 Solid floating on a free-surface===
-The following example shown schematically in Figure 13 represents a very interesting problem of fluid structure interaction when there is a weak interaction between the fluid and a large rigid deformation of the structure. In this case, there is also a free-surface problem, representing a schematic case of see-keeping in ship hydrodynamics.
+The following example shown schematically in Figure [[#img-13|13]] represents a very interesting problem of fluid structure interaction when there is a weak interaction between the fluid and a large rigid deformation of the structure. In this case, there is also a free-surface problem, representing a schematic case of see-keeping in ship hydrodynamics.
-<div id='img-14'></div>
+<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_Test_756558304-Fig11.png|351px|Solid floating on a free-surface. Initial geometry and point distribution.]]
+|style="padding: 10px;"| [[Image:Draft_Samper_187181919-Fig11.png|351px|Solid floating on a free-surface. Initial geometry and point distribution.]]
 |- style="text-align: center; font-size: 75%;"
-| colspan="1" | '''Figure 14:''' Solid floating on a free-surface. Initial geometry and point distribution.
+| colspan="1" style="padding: 10px;"| '''Figure 13:''' Solid floating on a free-surface. Initial geometry and point distribution.
 |}
-The example shows an initially stationary recipient with a floating piece of wood in which a wave is produced on the left side. The wave intercepts the wood piece producing a breaking wave and moving the floating wood. In this example the solid was represented by very viscous flows with a viscosity parameter order ten times greater than the water viscosity. Figure 14 shows the pressure contours and the free-surface position for different time steps.
+The example shows an initially stationary recipient with a floating piece of wood in which a wave is produced on the left side. The wave intercepts the wood piece producing a breaking wave and moving the floating wood. In this example the solid was represented by very viscous flows with a viscosity parameter order ten times greater than the water viscosity. Figure [[#img-14|14]] shows the pressure contours and the free-surface position for different time steps.
-<div id='img-15'></div>
+<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_Test_756558304-Fig12.png|351px|Solid floating on a free-surface. Pressure contours and free-surface positions for different time step.]]
+|style="padding: 10px;"| [[Image:Draft_Samper_187181919-Fig12.png|404px|Solid floating on a free-surface. Pressure contours and free-surface positions for different time step.]]
 |- style="text-align: center; font-size: 75%;"
-| colspan="1" | '''Figure 15:''' Solid floating on a free-surface. Pressure contours and free-surface positions for different time step.
+| colspan="1" style="padding: 10px;"| '''Figure 14:''' Solid floating on a free-surface. Pressure contours and free-surface positions for different time step.
 |}
@@ Line 956: / Line 1,012: @@
 This last example is also a case of fluid structure interaction. The solid is initially totally free and is falling down into a recipient with a fluid. In this example, the solid was modeled as a boundary condition for the fluid. Once the pressure and the viscous forces have been evaluated in the fluid, the solid is accelerated using Newton law. The solid has a mass and a gravity force concentrate in its gravity center.  The solid is considered to be ligh compared to the liquid weight.
-At the beginning the solid falls free due to the gravity forces (Figure 15). Once in contact with the water free-surface, the alpha-shape method recognizes the different boundary contours.  The pressure  and viscous forces are evaluated in all the domain and in particular on the solid contour. The flow forces introduce a negative acceleration in the vertical direction until , once the solid is completely inside the water, the vertical  velocity  becomes zero. Then,  Arquimides  forces bring the solid up to the free-surface. It is interesting to observe that there is a rotation of the solid. The reason is that the center of the floating forces is higher in the rotated position than in the initial ones.
+At the beginning the solid falls free due to the gravity forces (Figure [[#img-15|15]]). Once in contact with the water free-surface, the alpha-shape method recognizes the different boundary contours.  The pressure  and viscous forces are evaluated in all the domain and in particular on the solid contour. The flow forces introduce a negative acceleration in the vertical direction until , once the solid is completely inside the water, the vertical  velocity  becomes zero. Then,  Arquimides  forces bring the solid up to the free-surface. It is interesting to observe that there is a rotation of the solid. The reason is that the center of the floating forces is higher in the rotated position than in the initial ones.
-<div id='img-16'></div>
+<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_Test_756558304-Fig13.png|351px|Solid cube falling into a recipient with water. Free- surface position at different time steps.]]
+|style="padding: 10px;"| [[Image:Draft_Samper_187181919-Fig13.png|401px|Solid cube falling into a recipient with water. Free- surface position at different time steps.]]
 |- style="text-align: center; font-size: 75%;"
-| colspan="1" | '''Figure 16:''' Solid cube falling into a recipient with water. Free- surface position at different time steps.
+| colspan="1" style="padding: 10px;"| '''Figure 15:''' Solid cube falling into a recipient with water. Free- surface position at different time steps.
 |}
@@ Line 974: / Line 1,030: @@
 Thanks are also given to Mr. Romain Aubry for many useful discussions.
-===BIBLIOGRAPHY===
+==References==
 <div id="cite-1"></div>
-'''[[#citeF-1|[1]]]''' O.C. Zienkiewicz  and  R.C. Taylor,''The finite element method'', 5th Edition, 3 Volumes, Butterworth&#8211;Heinemann, 2000.
+[[#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.
+[[#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>
-'''[[#citeF-3|[3]]]''' A.J. Chorin,  “A numerical solution for solving incompressible viscous flow problems” ''J. Comp. Phys.'', '''2''', 12&#8211;26, 1967.
+[[#citeF-3|[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.
+[[#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.
+[[#citeF-5|[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.
+[[#citeF-6|[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.
+[[#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.
+[[#citeF-8|[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.
+[[#citeF-9|[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.
+[[#citeF-10|[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.
+[[#citeF-11|[11]]] A. Brooks and T.J.R. Hughes, “Streamline upwind/Petrov-Galerkin  formulation for convection dominated flows with particular emphasis  on the incompressible Navier-Stokes equations”, ''Comput. Methods Appl.  Mech. Engrg'', 32,  199&#8211;259, 1982.
 <div id="cite-12"></div>
-'''[[#citeF-12|[12]]]''' T.J.R. Hughes and M. Mallet, “A new finite element formulations  for computational fluid dynamics: III. The generalized streamline operator  for multidimensional advective-diffusive systems”, ''Comput Methods Appl.  Mech. Engrg.'', '''58''', pp. 305&#8211;328,  1986.
+[[#citeF-12|[12]]] T.J.R. Hughes and M. Mallet, “A new finite element formulations  for computational fluid dynamics: III. The generalized streamline operator  for multidimensional advective-diffusive systems”, ''Comput Methods Appl.  Mech. Engrg.'', 58, pp. 305&#8211;328,  1986.
 <div id="cite-13"></div>
-'''[13]''' P. Hansbo and a. Szepessy, “A velocity-pressure streamline diffusion finite element method for the incompressible Navier-Stokes equations”, ''Comput. Methods Appl. Mech. Engrg.'', '''84''', 175&#8211;192, 1990.
+[[#citeF-13|[13]]] P. Hansbo and a. Szepessy, “A velocity-pressure streamline diffusion finite element method for the incompressible Navier-Stokes equations”, ''Comput. Methods Appl. Mech. Engrg.'', 84, 175&#8211;192, 1990.
 <div id="cite-14"></div>
-'''[14]''' T.J.R. Hughes, L.P. Franca and M. Balestra, “A new finite element formulation for computational fluid dynamics. V Circumventing the Babuska-Brezzi condition: A stable Petrov-Galerkin formulation of the Stokes problem accomodating equal order interpolations”, ''Comput. Methods Appl. Mech. Engrg.'', '''59''', 85&#8211;89, 1986.
+[[#citeF-14|[14]]] T.J.R. Hughes, L.P. Franca and M. Balestra, “A new finite element formulation for computational fluid dynamics. V Circumventing the Babuska-Brezzi condition: A stable Petrov-Galerkin formulation of the Stokes problem accomodating equal order interpolations”, ''Comput. Methods Appl. Mech. Engrg.'', 59, 85&#8211;89, 1986.
 <div id="cite-15"></div>
-'''[15]''' L.P. Franca and S.L. Frey, “Stabilized finite element methods:  II. The incompressible Navier-Stokes equations”, ''Comput. Method Appl.  Mech. Engrg.'', Vol. '''99''', pp. 209&#8211;233, 1992.
+[[#citeF-15|[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.
+[[#citeF-16|[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.
+[[#citeF-17|[17]]]  M.A. Cruchaga and E. Oñate, “A finite element formulation for incompressible flow problems using a generalized streamline operator”, ''Comput. Methods in Appl. Mech. Engrg.'', 143, 49&#8211;67, 1997.
 <div id="cite-18"></div>
-'''[18]'''   M.A. Cruchaga and E. Oñate, “A generalized streamline finite element approach for the analysis of incompressible flow problems including moving surfaces”, ''Comput. Methods in Appl. Mech. Engrg.'', '''173''', 241&#8211;255, 1999.
+[[#citeF-18|[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.
+[[#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.
+[[#citeF-20|[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.
+[[#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.
+[[#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.
+[[#citeF-23|[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.
+[[#citeF-24|[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.
+[[#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.
+[[#citeF-26|[26]]]  R. Codina and O.C. Zienkiewicz, “CBS versus GLS stabilization of the incompressible Navier-Stokes equations and the role of the time step as stabilization parameter”, ''Communications in Numerical Methods in Engineering'', 18 (2), 99&#8211;112, 2002.
 <div id="cite-27"></div>
-'''[[#citeF-27|[27]]]'''  R. Codina and J. Blasco, “Stabilized finite element method for the transient Navier-Stokes equations based on a pressure gradient operator”,  ''Comput. Methods in Appl. Mech. Engrg.'', '''182''', 277&#8211;301, 2000.
+[[#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.
+[[#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.
+[[#citeF-29|[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.
+[[#citeF-30|[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.
+[[#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).
+[[#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).
+[[#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).
+[[#citeF-34|[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.
+[[#citeF-35|[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.
+[[#citeF-36|[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.
+[[#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>
-'''[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).
+[[#citeF-38|[38]]]  E. Oñate and J. García, “A finite element method for  fluid-structure interaction with surface waves using a finite calculus formulation”, ''Comput. MethodsAppl. Mech. Engrg.'', 191, 635&#8211;660, (2001).
 <div id="cite-39"></div>
-'''[39]'''  E. Oñate, “Possibilities of finite calculus in computational mechanics”, Submitted to ''Int. J. Num. Meth. Engng.'', 2002.
+[[#citeF-39|[39]]]  E. Oñate, “Possibilities of finite calculus in computational mechanics”, Submitted to ''Int. J. Num. Meth. Engng.'', 2002.
 <div id="cite-40"></div>
-'''[[#citeF-40|[40]]]''' E. Oñate and S. Idelsohn, A mesh free finite point method for advective-diffusive transport and fluid flow problems,, ''Computational Mechanics'', 21, 283&#8211;292, 1988.
+[[#citeF-40|[40]]] E. Oñate and S. Idelsohn, A mesh free finite point method for advective-diffusive transport and fluid flow problems,, ''Computational Mechanics'', 21, 283&#8211;292, 1988.
 <div id="cite-41"></div>
-'''[41]'''  E. Oñate, C. Sacco and  S. Idelsohn, “A finite point method for incompressible flow problems”, ''Computing and Visualization in Science'', '''2''', 67&#8211;75, 2000.
+[[#citeF-41|[41]]]  E. Oñate, C. Sacco and  S. Idelsohn, “A finite point method for incompressible flow problems”, ''Computing and Visualization in Science'', 2, 67&#8211;75, 2000.
 <div id="cite-42"></div>
-'''[42]''' J. García and E. Oñate, “An unstructured finite element solver for ship hydrodynamic problems”, in ''J. Appl. Mech.'', '''70''', 18&#8211;26, January 2003.
+[[#citeF-42|[42]]] J. García and E. Oñate, “An unstructured finite element solver for ship hydrodynamic problems”, in ''J. Appl. Mech.'', 70, 18&#8211;26, January 2003.
 <div id="cite-43"></div>
-'''[[#citeF-43|[43]]]'''  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.
+[[#citeF-43|[43]]]  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-44"></div>
-'''[[#citeF-44|[44]]]''' S.R. Idelsohn, E. Oñate, N. Calvo and F. del Pin, “The meshless finite element method”,  in ''Int. J. Num. Meth. Engng.'', '''58''', '''4''', 2003.
+[[#citeF-44|[44]]] S.R. Idelsohn, E. Oñate, N. Calvo and F. del Pin, “The meshless finite element method”,  in ''Int. J. Num. Meth. Engng.'', 58, 4, 2003.
 <div id="cite-45"></div>
-'''[[#citeF-45|[45]]]''' S.R. Idelsohn, E. Oñate, F. Del Pin and N. Calvo, “Lagrangian formulation: the only way to solve some free-surface fluid mechanics problems”, ''Fith World Congress on Computational Mechanics'', Mang HA, Rammerstorfer FG and Eberhardsteiner J. (eds), July 7&#8211;12, Viena, Austria, 2002, Web...
+[[#citeF-45|[45]]] S.R. Idelsohn, E. Oñate, F. Del Pin and N. Calvo, “Lagrangian formulation: the only way to solve some free-surface fluid mechanics problems”, ''Fith World Congress on Computational Mechanics'', Mang HA, Rammerstorfer FG and Eberhardsteiner J. (eds), July 7&#8211;12, Viena, Austria, 2002, Web...
 <div id="cite-46"></div>
-'''[46]''' S.R. Idelsohn, E. Oñate and F. Del Pin, “A lagrangian meshless finite element method applied to fluid-structure interaction problems”, in ''Computer and Structures'', '''81''', 655&#8211;671, 2003.
+[[#citeF-46|[46]]] S.R. Idelsohn, E. Oñate and F. Del Pin, “A lagrangian meshless finite element method applied to fluid-structure interaction problems”, in ''Computer and Structures'', 81, 655&#8211;671, 2003.
 <div id="cite-47"></div>
-'''[[#citeF-47|[47]]]'''  S.R. Idelsohn, N. Calvo and E. Oñate. Polyhedrization of an arbitrary point set. ''Comput. Method Appl. Mech. Engng.''. Accepted for publication 2003.
+[[#citeF-47|[47]]]  S.R. Idelsohn, N. Calvo and E. Oñate. Polyhedrization of an arbitrary point set. ''Comput. Method Appl. Mech. Engng.''. Accepted for publication 2003.
 <div id="cite-48"></div>
-'''[48]'''  H. Edelsbrunner and E.P. Mucke. Three dimensional alpha shapes. ºit ACM Trans. Graphics, '''13''', 43&#8211;72, 1999.
+[[#citeF-48|[48]]]  H. Edelsbrunner and E.P. Mucke. Three dimensional alpha shapes. ºit ACM Trans. Graphics, 13, 43&#8211;72, 1999.
 <div id="cite-49"></div>
-'''[49]'''  S. Koshizuka and Y. Oka. Moving particle semi-implicit method for fragmentation of incompressible fluid. ''Nuclear Engineering Science'', '''123''', 421&#8211;434, 1996.
+[[#citeF-49|[49]]]  S. Koshizuka and Y. Oka. Moving particle semi-implicit method for fragmentation of incompressible fluid. ''Nuclear Engineering Science'', 123, 421&#8211;434, 1996.
 <div id="cite-50"></div>
-'''[50]'''  R. Radovitzki and M. Ortiz. Lagrangian finite element analysis of Newtonian flows. ''Int. J. Num. Meth. Engng.'', '''43''', 607&#8211;619, 1998.
+[[#citeF-50|[50]]]  R. Radovitzki and M. Ortiz. Lagrangian finite element analysis of Newtonian flows. ''Int. J. Num. Meth. Engng.'', 43, 607&#8211;619, 1998.
 <div id="cite-51"></div>
-'''[51]'''  E.V. Laitone. The second approximation to cnoidal waves. ''Journal of Fluid Mechanics'', '''9''', 430, 1960.
+[[#citeF-51|[51]]]  E.V. Laitone. The second approximation to cnoidal waves. ''Journal of Fluid Mechanics'', 9, 430, 1960.