==A GENERAL STABILIZED FORMULATION FOR INCOMPRESSIBLE FLUID FLOW USING FINITE CALCULUS AND THE FINITE ELEMENT METHOD==

​

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

​

{|

|-

|International Centre for Numerical Methods in Engineering

|-

| Universidad Politécnica de Cataluña 

|-

| Gran Capitán s/n, 08034 Barcelona, Spain 

|-

| E–mail: [mailto:onate@cimne.upc.es onate@cimne.upc.es]

|-

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

|}

-->

​

==Abstract==

​

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

​

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

​

==1 INTRODUCTION==

​

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

​

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

​

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

​

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

​

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

​

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

​

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

​

==2 THE FINITE CALCULUS METHOD==

​

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

​

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

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

|-

|[[Image:Draft_Samper_187067199-majesus2.png|350px|Equilibrium of fluxes in a  balance domain of finite size]]

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

| colspan="1" | '''Figure 1:''' Equilibrium of fluxes in a  balance domain of finite size

|}

​

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

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

|-

| 

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

|-

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

|}

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

|}

​

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

​

Let us express now the fluxes <math display="inline">q_A</math> and <math display="inline">q_B</math> in terms of the flux at an arbitrary point <math display="inline">C</math> within the balance domain (Figure 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.

​

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

​

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

|-

| 

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

|-

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

|}

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

|}

​

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

​

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

|}

​

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;" | (6)

|}

​

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

​

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

|}

​

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

​

The underlined terms in Eqs.(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 time dimension can be  introduced in the FIC method by considering the balance equation in a space-time slab domain <span id='citeF-32'></span>[[#cite-32|32]]. Quite generally the FIC equation can be written for any problem in mechanics 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>r_i - \underline{\frac{h_j}{2}{\partial r_i \over \partial x_j}}\frac{h_j}{2}{\partial r_i \over \partial x_j}-\underline{\frac{\delta }{2}{\partial r_i \over \partial t}}\frac{\delta }{2}{\partial r_i \over \partial t}=0 \quad , \begin{array}{l}i=1,n_b\\ j=1,n_d \end{array} </math>

|}

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

|}

​

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

​

For example, in the case of the convection-diffusion problem <math display="inline">n_b=1</math>, Eq.(8) is particularized 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_j}{2}{\partial r \over \partial x_j}}\frac{h_j}{2}{\partial r \over \partial x_j}-\underline{\frac{\delta }{2}{\partial r \over \partial t}}\frac{\delta }{2}{\partial r \over \partial t}=0 \qquad ,\quad j=1,n_d </math>

|}

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

|}

​

with <math display="inline">r:=-\left(\displaystyle {\partial \phi  \over \partial t}+u_i {\partial \phi  \over \partial x_i}\right)+ \displaystyle \frac{d}{dx_k}\left(k \frac{d\phi }{dx_k}\right)+ Q</math>.  The modified Neumann boundary conditions in the FIC formulation can be expressed in the general case as

​

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

|-

| 

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

|-

| style="text-align: center;" | <math>q_{ij}n_j - t_i - \underline{h_jn_j r_i}h_jn_j r_i =0\quad \hbox{on } \Gamma _q \quad i=1,n_d </math>

|}

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

|}

​

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

​

In Eqs.(8)–(10) we have underlined once more the terms introduced by the FIC approach which are essential for deriving stabilized numerical formulations.

​

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

​

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

​

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

​

===Momentum===

​

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

|}

​

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

|}

​

where

​

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

|-

| 

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

|-

| style="text-align: center;" | <math>r_{m_i} = \rho \left({\partial u_i \over \partial t}+{\partial (u_iu_j) \over \partial x_j}\right)+ {\partial p \over \partial x_i}- {\partial s_{ij} \over \partial x_j}-b_i</math>

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

|-

| style="text-align: center;" | <math> r_d = {\partial u_i \over \partial x_i}\qquad i,j = 1, n_d </math>

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

|}

|}

​

Above <math display="inline">u_i</math> is the velocity along the ith global axis, <math display="inline">\rho </math> is the (constant) density of the fluid, <math display="inline">p</math> is the absolute pressure (defined positive in compression), <math display="inline">b_i</math> are the body forces and <math display="inline">s_{ij}</math> are the viscous deviatoric stresses related to the viscosity <math display="inline">\mu </math> by the standard expression

​

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

|-

| 

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

|-

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

|}

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

|}

​

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

​

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

|-

| 

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

|-

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

|}

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

|}

​

The FIC boundary conditions  are

​

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

|}

​

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

|-

| 

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

|-

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

|}

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

|}

​

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

​

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

​

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

​

Eqs.(11)&#8211;(18) are the starting  point for deriving stabilized finite element methods for solving the incompressible Navier-Stokes equations using equal order interpolation for the velocity and pressure variables <span id='citeF-37'></span>[[#cite-37|37]]. 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]].

​

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

​

===3.1 Stabilized integral forms===

​

From Eq.(11) it can be obtained (taking into account Eq.(15))

​

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

|-

| 

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

|-

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

|}

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

|}

​

where

​

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

|-

| 

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

|-

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

|}

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

|}

​

and

​

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

|-

| 

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

|-

| style="text-align: center;" | <math>\bar r_{m_i}=\rho \left({\partial u_i \over \partial t}+u_j {\partial u_i \over \partial x_j}\right)+ {\partial p \over \partial x_i}-{\partial  \over \partial x_j} (2\mu \dot \varepsilon _{ij}) - b_i</math>

|}

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

|}

​

Substituting Eq.(19) into Eq.(12) and retaining the terms involving the derivatives of <math display="inline">r_{m_i}</math> with respect to <math display="inline">x_i</math> only, leads to the following expression for the stabilized mass balance equation

​

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

|-

| 

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

|-

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

|}

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

|}

​

with

​

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

|-

| 

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

|-

| style="text-align: center;" | <math>\tau _i = \left({8\mu \over 3h_i^2}+{2\rho u_i\over h_i}\right)^{-1} </math>

|}

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

|}

​

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

|}

​

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

|}

​

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

​

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

|-

| 

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

|-

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

|}

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

|}

​

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

|-

| 

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

|-

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

|}

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

|}

​

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

​

Also in Eq.(25b) we will neglect hereonwards the third integral by assuming that <math display="inline">r_{m_i}</math> is negligible on the boundaries. The deviatoric stresses and the pressure terms in the first integral of Eq.(25a) 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 \left({\partial u_i \over \partial t}+u_j {\partial {u_i} \over \partial x_j}\right)+ \delta \dot \varepsilon _{ij}(s_{ij}- \delta _{ij}p )\right]d\Omega -  \int _{\Omega } \delta u_i b_i d\Omega - \int _{\Gamma _t} \delta u_i t_id\Gamma +</math>

|-

| style="text-align: center;" | <math> + \sum \limits _e \int _{\Omega ^e} {h_j\over 2}{\partial \delta u_i \over \partial x_j} r_{m_i} d\Omega =0 </math>

|}

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

|}

​

===Mass balance===

​

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

|-

| 

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

|-

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

|}

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

|}

​

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

​

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

|-

| 

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

|-

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

|}

|}

​

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

​

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

​

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

|-

| 

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

|-

| style="text-align: center;" | <math>\begin{array}{l}\displaystyle c_i = r_{m_i} -\rho u_j {\partial u_i \over \partial x_j}\\ \displaystyle \pi _i = r_{m_i} - {\partial p \over \partial x_i} \end{array} </math>

|}

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

|}

​

We can express <math display="inline">r_{m_i}</math> in  Eqs.(26) and (27) in terms of <math display="inline">c_i</math> and <math display="inline">\pi _i</math>, respectively which then become additional variables. The system of integral equations is now augmented in the necessary number of additional equations by imposing that the residual <math display="inline">r_{m_i}</math> vanishes (in average sense) for both forms given by Eqs.(28). This gives the final system of governing equation as:

​

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

|-

| 

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

|-

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

|-

| style="text-align: center;" | <math> + \sum \limits _e \int _{\Omega ^e} {h_k\over 2}{\partial (\delta u_i) \over \partial x_k} \left(\rho u_j {\partial {u_i} \over \partial x_j} + c_i\right)d\Omega =0 </math>

|}

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

|}

​

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

|-

| 

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

|-

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

|}

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

|}

​

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

|-

| 

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

|-

| style="text-align: center;" | <math>\int _\Omega \delta c_i \rho \left(\rho u_j {\partial {u_i} \over \partial x_j} + c_i\right)d\Omega =0 \qquad \hbox{no sum in }i </math>

|}

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

|}

​

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

|-

| 

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

|-

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

|}

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

|}

​

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

​

==4 FINITE ELEMENT DISCRETIZATION==

​

We choose <math display="inline">C^\circ </math> continuous linear interpolations of the velocities, the pressure, the convection projections <math display="inline">c_i</math> and the pressure gradient projections <math display="inline">\pi _i</math> over three node triangles (2D) and four node tetrahedra (3D). The linear interpolations are written as

​

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

|-

| 

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

|-

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

|}

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

|}

​

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

​

Substituting the approximations (33) into Eqs.(29&#8211;32) and choosing the Galerking form with <math display="inline">\delta u_i =q=\delta c_i=\delta \pi _i =N_i</math> leads to following system of discretized equations

​

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

|-

| 

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

|-

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

|}

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

|}

​

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

|-

| 

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

|-

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

|}

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

|}

​

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

|-

| 

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

|-

| style="text-align: center;" | <math>\displaystyle \hat {\boldsymbol C}\bar{\boldsymbol u}+ {\boldsymbol M}\bar {\boldsymbol c}={\boldsymbol 0}</math>

|}

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

|}

​

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

|-

| 

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

|-

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

|}

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

|}

​

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

​

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

|-

| 

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

|-

| style="text-align: center;" | <math>\begin{array}{l}\displaystyle M_{ij}= \int _{\Omega ^e} \rho N_iN_j d\Omega \quad ,\quad A_{ij}= \int _{\Omega ^e} N_i \rho {\boldsymbol u}^T {\boldsymbol \nabla } N_j d\Omega \\ \\ \displaystyle {\boldsymbol K}_{ij}= \int _{\Omega ^e}{\boldsymbol B}_i^T {\boldsymbol D}{\boldsymbol B}_j d\Omega \quad ,\quad \hat {\boldsymbol K}_{ij}= \int _{\Omega ^e} \rho {h_ku_l\over 2}\displaystyle {\partial N_i \over \partial x_k} \displaystyle {\partial N_j \over \partial x_l}d\Omega \\ \\ \displaystyle {\boldsymbol G}_{ij}= \int _{\Omega ^e}({\boldsymbol \nabla } N_i)N_j d\Omega \quad ,\quad {\boldsymbol \nabla }= \left[{\partial  \over \partial x_1},{\partial  \over \partial x_2}\right]^T \quad ,\quad  {C}_{ij} = \int _{\Omega ^e}{h_k \over 2} \displaystyle {\partial N_i \over \partial x_k}N_j d\Omega \, d\Omega \\ \\ \displaystyle L_{ij}=\int _{\Omega ^e} {\boldsymbol \nabla }^T N_i [\tau ] {\boldsymbol \nabla } N_j d\Omega \quad ,\quad [\tau ]= \left[\begin{matrix}\tau _1 &0\\ 0 & \tau _2\\\end{matrix}\right]\quad ,\quad  \hat {C}_{ij}=\int _{\Omega ^e} N_i \rho ^2 u_k \displaystyle {\partial N_j \over \partial x_k}d\Omega \\ \\ \displaystyle {\boldsymbol Q}= [{\boldsymbol Q}^1,{\boldsymbol Q}^2] \quad ,\quad  \displaystyle Q_{ij}^k = \int _{\Omega ^e}\tau _k {\partial N_i \over \partial x_k} N_jd\Omega \\ \\ \displaystyle \hat {\boldsymbol M}= \left[\begin{matrix}\hat {\boldsymbol M}^1 &{\boldsymbol 0}\\  {\boldsymbol 0} & \hat {\boldsymbol M}^2\\\end{matrix}\right]\quad ,\quad  \hat {\boldsymbol M}_{ij} = \int _{\Omega ^e} \tau _k N_i N_j d\Omega \\ \\ \displaystyle {\boldsymbol f}_i = \int _{\Omega ^e} N_i {\boldsymbol b}d\Omega + \int _{\Gamma ^e}N_i {\boldsymbol t} d\Gamma \quad ,\quad  {\boldsymbol b}=[b_1,b_2]^T \quad ,\quad {\boldsymbol t}=[t_1,t_2]^T \end{array} </math>

|}

|}

​

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;" | (35)

|}

​

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

​

===4.1 Quasi-implicit scheme===

​

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

​

===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 A}^{n+\theta _1,i-1}+  {\boldsymbol K} + \hat {\boldsymbol K}^{n+\theta _1,i-1})\bar{\boldsymbol u}^{n+\theta _1,i-1}-</math>

|-

| style="text-align: center;" | <math> - {\boldsymbol G}\bar{\boldsymbol p}^{n+\theta _2,i-1} + {\boldsymbol C} \bar {\boldsymbol c}^{n+\theta _3,i-1}-{\boldsymbol f}^{n+1}] </math>

|}

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

|}

​

===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 _4,i-1}] </math>

|}

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

|}

​

===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 c}^{n+1,i}= - {\boldsymbol M}^{-1} \hat {\boldsymbol C}^{n+1,i}\bar {\boldsymbol u}^{n+1,i} </math>

|}

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

|}

​

===Step 4===

​

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

|}

​

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

​

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

​

Steps 1, 3 and 4 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.(37)&#8211;(40) can be made equal to zero. Note that for <math display="inline">\theta _2 =0</math> the algorithm is unconditionally unstable. A particularity interesting and simple semi-implicit form is obtained by making <math display="inline">\theta _1 = \theta _3=\theta _4=0</math>. Now all steps can be solved explicitely with exception of Step 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 _4,i-1}] </math>

|}

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

|}

​

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

​

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

​

===4.2 Fractional step method===

​

An alternative algorithm can be obtained by splitting the pressure from the momentum equations as follow

​

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

|-

| 

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

|-

| style="text-align: center;" | <math>\bar {\boldsymbol u}^* = \bar {\boldsymbol u}^n -\Delta t {\boldsymbol M}^{-1}[({\boldsymbol A}^{n+\theta _1}+  {\boldsymbol K} + \hat {\boldsymbol K}^{n+\theta _1})\bar{\boldsymbol u}^{n+\theta _1}-\alpha {\boldsymbol G}\bar{\boldsymbol p}^{n} + {\boldsymbol C} \bar {\boldsymbol c}^{n+\theta _3}-{\boldsymbol f}^{n+1}]</math>

|}

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

|}

​

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

|-

| 

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

|-

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

|}

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

|}

​

In Eq.(42a) <math display="inline">\alpha </math> is a variable taking values equal to zero or one. For <math display="inline">\alpha =0</math>, <math display="inline">\delta p \equiv p^{n+1}</math> and for <math display="inline">\alpha =1</math>, <math display="inline">\delta p =\Delta p</math>. Note that in both cases the sum of Eqs.(42a) and (42b) gives the time discretization of the momentum equations with the pressures computed at <math display="inline">t^{n+1}</math>. The value of <math display="inline">\bar {\boldsymbol u}^{n+1}</math> from Eq.(42b) is substituted now into Eq.(34b) 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 _4}=0</math>

|}

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

|}

​

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

​

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

|-

| 

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

|-

| style="text-align: center;" | <math>{\boldsymbol G}^T {\boldsymbol M}^{-1}{\boldsymbol G}=\hat {\boldsymbol L} \quad \hbox{with } \hat{\boldsymbol L} \simeq \int _{\Omega ^e} {1\over \rho } {\boldsymbol \nabla }^T N_i {\boldsymbol \nabla } N_j \,d\Omega  </math>

|}

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

|}

​

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

​

'''Step 1''' Compute <math display="inline">{\boldsymbol u}^*</math> explicitely from Eq.(41a) 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> from Eq.(43a) 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;" | (44)

|}

​

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

​

<br/>

​

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

​

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

|-

| 

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

|-

| style="text-align: center;" | <math>\bar{\boldsymbol c}^{n+1}=-{\boldsymbol M}_d^{-1}\hat {\boldsymbol C}^{n+1}\bar{\boldsymbol u}^{n+1} </math>

|}

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

|}

​

'''Step 5''' Compute <math display="inline">  \bar{\boldsymbol \pi }^{n+1}</math> explicitely from Eq.(40) 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;" | (46)

|}

​

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

​

===4.3 Stokes flow===

​

The formulation for a Stokes flow can be readily obtained simply by neglecting the convective terms in the general Navier-Stokes formulation. This also implies neglecting the convective stabilization terms in the momentum equations and, consequently, the convective projection variables are not larger necessary. Also the intrinsic time parameters <math display="inline">\tau _i</math> take now the simpler form (see Eq.(22)):

​

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

|-

| 

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

|-

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

|}

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

|}

​

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

​

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

|-

| 

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

|-

| style="text-align: center;" | <math>\begin{array}{l}\displaystyle {\boldsymbol M}\dot{\bar {\boldsymbol u}} + {\boldsymbol K}\bar{\boldsymbol u} - {\boldsymbol G}\bar {\boldsymbol p}={\boldsymbol f}\\  \displaystyle {\boldsymbol G}^T \bar {\boldsymbol u} + {\boldsymbol L}\bar {\boldsymbol p}+{\boldsymbol Q}\bar {\boldsymbol \pi }={\boldsymbol 0}\\ \displaystyle {\boldsymbol Q}^T \bar {\boldsymbol p} + \hat {\boldsymbol M}\bar {\boldsymbol \pi }={\boldsymbol 0} \end{array} </math>

|}

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

|}

​

The iterative algorithm of Section 4.1 can now be implemented. Convergence is now faster due to the absence of the non linear convective terms in the momentum equation.

​

The steady-state form of Eqs.(48) 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;" | (49)

|}

​

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

​

A reduced velocity-pressure formulation can be obtained by eliminating the pressure gradient projection variables <math display="inline">\bar {\boldsymbol \pi }</math> from the last equation to give

​

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

|-

| 

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

|-

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

|}

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

|}

​

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

​

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

​

===5.1 General coupled solution scheme===

​

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

​

===Step  A1.  Solve  for  the  movement  of  the  structure  due  to  the  fluid  flow  forces===

​

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

​

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

|-

| 

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

|-

| style="text-align: center;" | <math>{\boldsymbol M}_s \ddot {\boldsymbol d}+ {\boldsymbol K}_s {\boldsymbol d}={\boldsymbol f}_{ext} </math>

|}

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

|}

​

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

​

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

​

===Step A2.Computethenew position of the mesh nodes===

​

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

​

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

​

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

​

Chiandussi, Bugeda and Oñate <span id='citeF-45'></span>[[#cite-45|45]] have proposed a simple method for  the movement of mesh nodes ensuring minimum element distortion. The method is based on the iterative solution of a fictious linear elastic problem on the mesh domain. In order to minimize the mesh deformation the “elastic” properties of each mesh element are appropiately selected so that elements suffering greater movements are stiffer. The basis of the method is given below.

​

Let us consider an elastic domain with arbitrary homogeneous isotropic elastic properties characterized by the Young modulus <math display="inline">\bar E</math> and the Poisson coefficient <math display="inline">\nu </math>. Once a discretized finite element problem has been solved using, for instance, standard <math display="inline">C^o</math> linear triangles (in 2D) or linear tetraedra (in 3D), the principal stresses <math display="inline">{}^1\sigma _i</math> at the center of each element can be obtained. For 3D problems

​

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

|-

| 

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

|-

| style="text-align: center;" | <math>{}^1\sigma _i ={\bar E (1-\nu )\over ( 1+\nu )(1-2\nu )}\left[\varepsilon _i + {\nu \over  1-\nu } (\varepsilon _j + \varepsilon _k)\right]\qquad i,j,k=1,2,3 ~~i\not =j\not = k </math>

|}

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

|}

​

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

​

Let us assume now that a uniform strain field <math display="inline">\varepsilon _i =\bar \varepsilon </math> throughout the mesh is sought. The principal stresses are then given by

​

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

|-

| 

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

|-

| style="text-align: center;" | <math>{}^2\sigma _i =  E \bar \varepsilon \qquad i=1,2,3 </math>

|}

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

|}

​

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

​

A number of criteria can be now used to find the value of <math display="inline">E</math>. An effective approach found in <span id='citeF-45'></span>[[#cite-45|45]] is to make equal the element strain energy densities in both analysis. Thus (for <math display="inline">\nu =0</math>)

​

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

|-

| 

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

|-

| style="text-align: center;" | <math>U_1 =  {{}^1\sigma _i \varepsilon _i\over 2} = {\bar E [(\varepsilon _1^2+\varepsilon _2^2 + \varepsilon _3^2)]\over 2} </math>

|}

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

|}

​

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

|-

| 

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

|-

| style="text-align: center;" | <math>U_2 =  {{}^2\sigma _i \varepsilon _i \over 2}= {3E\bar \varepsilon ^2\over 2}</math>

|}

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

|}

​

Equaling eqs.(54a) and (54b) gives the sought Young modulus <math display="inline">E</math> as

​

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

|-

| 

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