==LAGRANGIAN FORMULATION FORINCOMPRESSIBLE FLUIDS  USING FINITE CALCULUS AND THE FINITE ELEMENT METHOD==

​

'''E. Oñate<math>^1</math> 

S.R. Idelsohn<math>^{1,2}</math> and F. Del Pin<math>^2</math>  '''

​

{|

|-

| <math>^1</math> International Center for Numerical

|-

| Methods in Engineering (CIMNE)

|-

| Universidad Politecnica de Cataluna

|-

| Campus Norte UPC, 08034 Barcelona, Spain 

|-

| Email: [mailto:onate@cimne.upc.es onate@cimne.upc.es] 

|-

| Web page: [http://www.cimne.upc.es http://www.cimne.upc.es]   

<math>^2</math> CIMEC

|-

| Universidad Nacional del Litoral

|-

| Güemes 3450

|-

| 3000 Santa Fe

|-

| Argentina

|-

| Email: [mailto:sergio@ceride.gov.ar sergio@ceride.gov.ar]  

|}

-->

​

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

​

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

​

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 <math display="inline">\Omega </math> of length <math display="inline">L</math>. The equation of balance of fluxes in a subdomain of size <math display="inline">d</math> belonging to <math display="inline">\Omega </math> (Figure 1) is written as

​

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>q_A - q_B=0 </math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | (1)

|}

​

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

​

<div id='img-1'></div>

{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"

|-

|[[Image:Draft_Samper_800412471_9634_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

|}

​

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

​

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>q_A= q_C - d_1 \frac{dq}{d x}\vert _C+ \frac{d^2_1}{2}\frac{d^2q}{dx^2}\vert _C + O(d^3_1)\quad ,\quad  q_B= q_C + d_2 \frac{dq}{d x}\vert _C+\frac{d^2_2}{2}\frac{d^2q}{dx^2}\vert _C + O(d^3_2) </math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | (2)

|}

​

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

​

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>\frac{dq}{dx}-\underline{\frac{h}{2} \frac{d^2q}{dx^2}}\frac{h}{2} \frac{d^2q}{dx^2}=0 </math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | (3)

|}

​

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

​

Standard calculus theory assumes  that the domain <math display="inline">d</math> is of infinitesimal size and the resulting balance equation is simply <math display="inline">{dq\over dx}=0</math>. We will relax this assumption and allow the balance domain to have a ''finite size''. The new balance equation (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>.

​

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.

​

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

​

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>r- \underline{\frac{h}{2} \frac{dr}{dx}}\frac{h}{2} \frac{dr}{dx}=0 </math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | (4)

|}

​

with

​

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>r = -u \frac{d\phi }{dx}+ \frac{d}{dx}\left(k \frac{d\phi }{dx}\right)+ Q </math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | (5)

|}

​

where <math display="inline">Q</math> is the external source. 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]]

​

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>k \frac{d\phi }{dx}+ q_p- \underline{\frac{h}{2}r}\frac{h}{2}r=0 \quad \hbox{at } \, \Gamma _q </math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | (6)

|}

​

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

​

The underlined terms in Eqs.(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 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===

​

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>r_{m_i} - \underline{{1\over 2} h_j{\partial r_{m_i} \over \partial x_j}}{1\over 2} h_j{\partial r_{m_i} \over \partial x_j}=0  </math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | (7)

|}

​

===Mass balance===

​

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>r_d - \underline{{1\over 2} h_j {\partial r_d \over \partial x_j}}{1\over 2} h_j {\partial r_d \over \partial x_j}=0  </math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | (8)

|}

​

where

​

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>r_{m_i} = \rho {\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)

|-

| 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

​

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>s_{ij}=2\mu \left(\dot \varepsilon _{ij} - \delta _{ij} {1\over 3} {\partial u_k \over \partial x_k}\right) </math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | (11)

|}

​

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

​

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

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

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | (12)

|}

​

The FIC boundary conditions  are

​

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>n_j \sigma _{ij} -t_i + \underline{{1\over 2} h_j n_j r_{m_i}}{1\over 2} h_j n_j r_{m_i}=0 \quad \hbox{on }\Gamma _t </math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | (13)

|}

​

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>u_j - u_j^p =0 \quad \hbox{on }\Gamma _u </math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | (14)

|}

​

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

​

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

​

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

​

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

​

===3.1 Stabilized integral forms===

​

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

​

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>{\partial r_d \over \partial x_i}=-{1\over a_i} \left[\bar r_{m_i}-{h_j\over 2} {\partial r_{m_i} \over \partial x_j}\right]-  {\rho u_i h_k\over 2a_i} {\partial r_d \over \partial x_k}\quad i,j=1,n_d\quad ,\quad  k\not =i </math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | (15)

|}

​

where

​

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>a_i = {2\mu \over 3} </math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | (16a)

|}

​

and

​

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>\bar r_{m_i}=\rho {\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)

|}

​

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

​

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| 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="width: 5px;text-align: right;white-space: nowrap;" | (17)

|}

​

with

​

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>\tau _i = {3h_i^2\over 8\mu } </math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | (18)

|}

​

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 weighted residual form of the momentum and mass balance equations (Eqs.(7) and (17)) is written as

​

===Momentum===

​

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>\int _\Omega \delta u_i \left[r_{m_i} - {h_j\over 2} {\partial r_{m_i} \over \partial x_j}\right]+ \int _{\Gamma _t} \delta u_i (\sigma _{ij} n_j - t_i + {h_j \over 2} n_j r_{m_i}) d\Gamma =0 </math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | (19)

|}

​

===Mass balance===

​

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>\int _\Omega q \left[r_d - \sum \limits _{i=1}^{n_d} \tau _i {\partial r_{m_i} \over \partial x_i}\right]d\Omega =0 </math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | (20)

|}

​

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

​

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>\int _\Omega \delta u_i r_{m_i} + \int _{\Gamma _t} \delta u_i (\sigma _{ij} n_j - t_i)d\Gamma + \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)

|}

​

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>\int _\Omega q r_d d\Omega + \int _\Omega \left[\sum \limits _{i=1}^{n_d}\tau _i {\partial q \over \partial x_i}r_{m_i} \right]d\Omega - \int _\Gamma \left[\sum \limits _{i=1}^{n_d} q \tau _i n_i r_{m_i}\right]d\Gamma =0</math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | (21b)

|}

​

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.

​

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

​

===Momentum===

​

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>\int _\Omega \left[\delta u_i\rho  {\partial u_i \over \partial t}  + \delta \dot \varepsilon _{ij}(s_{ij}- \delta _{ij}p )\right]d\Omega -  \int _{\Omega } \delta u_i b_i d\Omega - \int _{\Gamma _t} \delta u_i t_id\Gamma  =0 </math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | (22)

|}

​

===Mass balance===

​

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>\int _\Omega q {\partial u_i \over \partial x_i} d\Omega + \int _\Omega \left[\sum \limits _{i=1}^{n_d} \tau _i {\partial q \over \partial x_i} r_{m_i}\right]d\Omega =0 </math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | (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 <math display="inline">\displaystyle {\partial u_i \over \partial x_i}</math>). Also in Eq.(22)

​

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

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

|}

|}

​

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

​

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

​

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>\pi _i = r_{m_i} - {\partial p \over \partial x_i} </math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | (24)

|}

​

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:

​

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>\int _\Omega \left[\delta u_i\rho {\partial u_i \over \partial t} + \delta \dot \varepsilon _{ij}(s_{ij}- \delta _{ij}p )\right]d\Omega -  \int _{\Omega } \delta u_i b_i d\Omega - \int _{\Gamma _t} \delta u_i t_id\Gamma =0 </math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | (25)

|}

​

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>\int _\Omega q {\partial u_i \over \partial x_i} d\Omega + \int _\Omega \sum \limits _{i=1}^{n_d} \tau _i {\partial q \over \partial x_i} \left({\partial p \over \partial x_i}+\pi _i\right)d\Omega =0 </math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | (26)

|}

​

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>\int _\Omega \delta \pi _i \tau _i \left({\partial p \over \partial x_i}+\pi _i\right)d\Omega =0\qquad \hbox{no sum in }i </math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | (27)

|}

​

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.

​

==4 FINITE ELEMENT DISCRETIZATION==

​

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

​

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>u_i = \sum \limits _{j=1}^n N_j \bar u_i^j \quad , \quad p = \sum \limits _{j=1}^n N_j \bar p^j\quad , \quad \pi _i = \sum \limits _{j=1}^n N_j \bar \pi _i^j </math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | (28)

|}

​

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

​

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>\displaystyle {\boldsymbol M}\dot{\bar {\boldsymbol u}} + {\boldsymbol K} \bar {\boldsymbol u} - {\boldsymbol G}\bar {\boldsymbol p} ={\boldsymbol f}</math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | (29a)

|}

​

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>\displaystyle {\boldsymbol G}^T \bar {\boldsymbol u} + {\boldsymbol L}\bar {\boldsymbol p}+{\boldsymbol Q}\bar {\boldsymbol \pi }={\boldsymbol 0}</math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | (29b)

|}

​

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>\displaystyle {\boldsymbol Q}^T \bar {\boldsymbol p} + \hat {\boldsymbol M}\bar {\boldsymbol \pi }={\boldsymbol 0}</math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | (29c)

|}

​

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

​

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>\begin{array}{l}\displaystyle M_{ij}= \int _{\Omega ^e} \rho N_iN_j d\Omega \quad ,\quad \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>

|}

|}

​

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

​

In above <math display="inline">{\boldsymbol B}_i</math> is the standard strain rate matrix and <math display="inline">{\boldsymbol D}</math> the deviatoric constitutive matrix (assuming <math display="inline">\displaystyle {\partial u_i \over \partial x_i}=0</math>). For 2D problems

​

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>{\boldsymbol B}_i =\left[\begin{matrix}\displaystyle {\partial N_i \over \partial x_1}&0\\ 0 & \displaystyle {\partial N_i \over \partial x_2}\\ \displaystyle {\partial N_i \over \partial x_2} & \displaystyle {\partial N_i \over \partial x_1}\\\end{matrix}\right]\quad ,\quad {\boldsymbol D} =\mu \left[\begin{matrix}2&0&0\\ 0&2&0\\ 0&0&1\\\end{matrix}\right] </math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | (30)

|}

​

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

​

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>\left[\begin{matrix}{\boldsymbol K}&-{\boldsymbol G}&{\boldsymbol 0}\\  -{\boldsymbol G}^T & - {\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)

|}

​

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

​

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>\left[\begin{matrix}{\boldsymbol K}&-{\boldsymbol G}\\ -{\boldsymbol G}^T & - ({\boldsymbol L} -{\boldsymbol Q}\hat {\boldsymbol M}^{-1}{\boldsymbol Q}^T)\\\end{matrix}\right]\left\{\begin{matrix}\bar {\boldsymbol u}\\ \bar {\boldsymbol p}\\\end{matrix}\right\}= \left\{\begin{matrix}{\boldsymbol f}\\ {\boldsymbol 0}\\\end{matrix}\right\} </math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | (32)

|}

​

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

​

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

​

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>\bar{\boldsymbol u}^{n+1,i} = \bar{\boldsymbol u}^n - \Delta t {\boldsymbol M}^{-1} [{\boldsymbol 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)

|}

​

===Step 2===

​

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>\bar{\boldsymbol p}^{n+1,i}=-{\boldsymbol L}^{-1} [{\boldsymbol G}^T \bar {\boldsymbol u}^{n+1,i}+{\boldsymbol Q}\bar {\boldsymbol \pi }^{n+\theta _3,i-1}] </math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | (34)

|}

​

===Step 3===

​

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>\bar {\boldsymbol \pi }^{n+1,i}= - \hat {\boldsymbol M}^{-1}{\boldsymbol Q}^T \bar {\boldsymbol p}^{n+1,i} </math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | (35)

|}

​

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.

​

===Step 4===

​

Update the mesh nodes in a lagrangian manner as

​

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>{\boldsymbol x}_i^{n+1} = {\boldsymbol x}_i^{n}+\bar {\boldsymbol u}_i^{n+1} \Delta t </math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | (36)

|}

​

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

​

<br/>

​

<math display="inline">\boldsymbol {Some~practical~remarks}</math>

​

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

​

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>\bar {\boldsymbol p}^{n+1,i}=-[{\boldsymbol L} + \hat {\boldsymbol L}]^{-1}[{\boldsymbol G}^T \bar {\boldsymbol u}^{n+1,i}+\hat {\boldsymbol L}\bar {\boldsymbol p}^{n+1,i-1} +{\boldsymbol Q}\bar {\boldsymbol \pi }^{n+\theta _3,i-1}] </math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | (37)

|}

​

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

​

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

​

===4.2 Fractional step method===

​

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

​

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>\bar {\boldsymbol u}^* = \bar {\boldsymbol u}^n -\Delta t {\boldsymbol M}^{-1}[{\boldsymbol 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)

|}

​

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>\bar{\boldsymbol u}^{n+1}= \bar{\boldsymbol u}^* + \Delta t {\boldsymbol M}^{-1}{\boldsymbol G}\delta \bar{\boldsymbol p}</math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | (39b)

|}

​

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

​

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>{\boldsymbol G}^T\bar {\boldsymbol u}^* + \Delta t {\boldsymbol G}^T {\boldsymbol M}^{-1}{\boldsymbol G}\delta \bar{\boldsymbol p} + {\boldsymbol L} \bar{\boldsymbol p}^{n+1}+ {\boldsymbol Q}\bar {\boldsymbol \pi }^{n+\theta _3}=0</math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | (40a)

|}

​

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

​

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>{\boldsymbol G}^T {\boldsymbol M}^{-1}{\boldsymbol G}=\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)

|}

​

A semi-implicit algorithm can thus be derived as follows.<br/>

​

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

​

<br/>

​

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

​

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>\delta \bar {\boldsymbol p} =-({\boldsymbol L}+\Delta t \hat  {\boldsymbol L})^{-1} [{\boldsymbol G}^T\bar{\boldsymbol u}^* +{\boldsymbol Q}\bar {\boldsymbol \pi }^{n+\theta _4}+ \alpha {\boldsymbol L} \bar {\boldsymbol p}^n]</math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | (41a)

|}

​

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

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

|}

​

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

​

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

​

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>\bar{\boldsymbol \pi }^{n+1}=- \hat {\boldsymbol M}_d^{-1} {\boldsymbol Q}^T \bar {\boldsymbol p}^{n+1} </math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | (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 <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>.

​

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.

​

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

​

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

​

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

|}

​

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 <math display="inline">h</math>. 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 <math display="inline">{\boldsymbol h}</math> has the direction of the velocity field <span id='citeF-32'></span>[[#cite-32|32]]. This gives

​

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>{\boldsymbol h}=h_s {{\boldsymbol u}\over {u}} </math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | (43)

|}

​

where <math display="inline">u=\vert {\boldsymbol u}\vert </math> and <math display="inline">h_s</math> as the “streamline” characteristic length

​

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>h_s=\max ({\boldsymbol l}^T_j {\boldsymbol u})/{u} \quad , \quad  j=1,n_s </math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | (44)

|}

​

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

​

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

​

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

|- style="text-align: center; font-size: 75%;"

| colspan="1" | '''Figure 4:''' cont.

|}

​

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 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 3D solution of the same problem is shown in Figure 4.

​

<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_Test_756558304-Fig7_1.png|600px|]]

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

|- style="text-align: center; font-size: 75%;"

| colspan="2" | '''Figure 5:''' Water column 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.

​

<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-Fig4ArtBathe.png|600px|Rotating water mill.]]

|- style="text-align: center; font-size: 75%;"

| colspan="1" | '''Figure 6:''' 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.

​

<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-Fig3.png|351px|Sloshing. Initial point distribution.]]

|- style="text-align: center; font-size: 75%;"

| colspan="1" | '''Figure 7:''' Sloshing. Initial point distribution.

|}

​

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

|}

​

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.

​

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

|}

​

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

|}

​

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