P1/P0+ elements for incompressible flows with discontinuous material properties

Revision as of 15:33, 5 October 2018

Published in Comput. Meth. Appl. Mech. Engng., Vol. 293, pp. 191-206, 2015

Abstract

We present a 3-noded triangle and a 4-noded tetrahedra with a continuous linear velocity and a discontinuous linear pressure field formed by the sum of an unknown constant pressure field and a prescribed linear field that satisfies the steady state momentum equations for a constant body force. The elements are termed P1/P0+ as the “effective” pressure field is linear, although the unknown pressure field is piecewise constant within each element. The elements have an excellent behaviour for incompressible viscous flow problems with discontinuous material properties formulated in either Eulerian or Lagrangian descriptions. The necessary numerical stabilization for dealing with the inf-sup condition imposed by the incompressibility constraint and high convective effects (in Eulerian flows) is introduced via the Finite Calculus (FIC) approach. For the sake of clarity, the element derivation is presented first for the simpler Stokes equations written in the standard Eulerian frame. The extension of the formulation to the Navier-Stokes equations written in the Eulerian and Lagrangian frameworks is straightforward and is presented in the second part of the paper.

The efficiency and accuracy of the new P1/P0+ triangle is verified by solving a set of incompressible multifluid flow problems using a Lagrangian approach and a classical Eulerian description. The excellent performance of the new triangular element in terms of mass conservation and general accuracy for analysis of fluids with discontinuous material properties is highlighted.

Keywords: P1/P0+ elements, Incompressible flow, Discontinuous material properties, Multifluids

1 INTRODUCTION

Preservation of mass is a great challenge in the numerical study of incompressible flow problems. Mass losses can be introduced by the so-called stabilization terms which are typically added to the discretized forms of the momentum and mass balance equations in order to account for large values of the convective acceleration terms in the momentum equations in high Reynolds number flows and to satisfy the inf-sup condition imposed by the incompressibility constraint when equal order interpolation of the velocities and the pressure are used in mixed finite element methods [4,11,39].

Mass loss also typically occurs in the numerical study of free-surface incompressible flow problems using Eulerian and Lagrangian descriptions. In both cases, the inaccuracy in capturing the free surface can lead to considerable mass losses unless special numerical schemes are used [22].

An additional source of mass loss occurs in the numerical analysis of the so-called multifluid problems with discontinuous changes in the viscosity and/or the density in parts of the domain. Most numerical methods have difficulties for accurately capturing the jumps in the pressure and/or the pressure gradient at the interfaces between the different fluids [20,21,25].

The numerical study of multifluids via the FEM and similar computational techniques has been the subject of much research in last decades. Several authors have proposed alternative stabilized FEM procedures for accounting for the discontinuity in the pressure (and/or the pressure gradient) at the interface of fluids with different viscosity (and/or pressure). Among these formulations we note those based in injecting a discontinuous pressure field within the appropriate elements [2] and those based on using a stabilized formulation based on the introduction of stabilization terms including the jumps in the pressure and the viscous terms at the element boundaries [3,5,6,7,9,12,14,17,23].

In this work we present a new 3-noded triangle (and the 4-noded tetrahedron counterpart) with a continuous linear velocity and a discontinuous linear pressure field formed by the sum of an unknown constant pressure field and a prescribed linear field that satisfies the steady state momentum equations for a constant body force. The so-called P1/P0+ elements have an excellent behaviour for incompressible viscous flow problems formulated in Eulerian and Lagrangian descriptions. For the sake of clarity, the elements derivation is presented first for the simpler Stokes equations written in the standard Eulerian frame. The extension of the formulation to the Navier-Stokes equations written in the Eulerian and Lagrangian frameworks is straightforward and is presented in the second part of the paper. A motivation of this work is the study of multifluids problems using a Lagrangian formulation via the Particle Finite Element Method (PFEM) [10],[18]–[21],[25,32],[36]–[38] and similar procedures.

The success of the P1/P0+ formulation lays in the consistent derivation of a residual-based expression of the mass balance equation using the Finite Calculus (FIC) technique. The FIC approach in mechanics is based on expressing the equations of balance of mass and momentum in a domain of finite size and retaining higher order terms in the Taylor series expansion typically used for expressing the change in the transported variables within the balance domain. In addition to the standard terms of infinitesimal theory, the FIC forms of the balance equations contain derivatives of the classical differential equations in mechanics multiplied by characteristic distances in space and time. Examples of stabilized FIC-FEM formulations in fluid and solid mechanics can be found in [26]–[31],[33]–[35]. In our work we use the FIC forms of the mass balance equation in space and time for obtaining a variational residual form useful for finite element analysis.

The discretized stabilized variational form for the mass balance equation using a P1/P0+ approximation for the velocity/pressure variables involves the jumps in the viscous stresses, the pressure, the surface tractions and the acceleration term in the normal direction to each side (or face) of the elements. These stabilization terms resemble those proposed by other authors for similar fluid flow problems [3,9,12,14,17,23]. The method presented in this paper yields a consistent and extended form of the stabilization terms that has proven to give a superior behaviour in terms of mass conservation and overall stability in multifluids problems.

The lay-out of the paper is the following. In the next section we present the basic equations for an incompressible Stokes fluid. Next we derive the stabilized variational FIC form of the mass balance equation. Then the P1/P0+ finite element discretization for the 3-noded triangle and the 4-noded tetrahedron is presented and the key matrices and vectors of the discretized system of equations are given. Details of the solution of the FEM equations are given.

The extension of the general stabilized FIC-FEM formulation to Navier-Stokes flows using a standard Eulerian approach and a Lagrangian description is outlined.

The efficiency and accuracy of the new P1/P0+ triangle is verified by solving a set of transient incompressible multifluid flow problems using a Lagrangian approach and a steady state problem via a classical Eulerian description. The excellent performance of the P1/P0+ triangle in terms of mass conservation and general accuracy for fluid flow problems with discontinuous material properties is highlighted.

2 BASIC EQUATIONS

We write the governing equations for an incompressible Stokes flow problem as follows.

Momentum equations

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): r_{m_i}:= \rho \frac{\partial v_i}{\partial t}-{\partial \sigma _{ij} \over \partial x_j}-b_i=0\quad , \quad i,j=1,n_s \quad \hbox{in } \Omega

(1)

In Eq.(1), Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle v_i}

and Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle b_i}
are the velocity and body force components  along the Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle i}

th Cartesian axis, Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle n_s}

is the number of space dimension (i.e. Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle n_s=3}
for a 3D problem), Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \Omega }
is the analysis domain and Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \sigma _{ij}}
are the Cauchy stresses that are split in the deviatoric (Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle s_{ij}}

) and pressure (Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle p} ) components as

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): \sigma _{ij} = s_{ij}+ p \delta _{ij}

(2)

where Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \delta _{ij}}

is the Kronecker delta. Note that the pressure is assumed to be positive for a tension state.

Summation of terms with repeated indices is assumed in Eq.(1) and in the following, unless otherwise specified.

Constitutive equation and volumetric strain rates

The relationship between the deviatoric stresses and the strain rates has the standard form for a Newtonian fluid,

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): s_{ij} = 2\mu \left(\varepsilon _{ij} - {1 \over 3}\varepsilon _v \delta _{ij}\right)

(3)

where the strain rates Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \varepsilon _{ij}}

are related to the velocities by

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): \varepsilon _{ij} = {1 \over 2}\left({\partial v_i \over \partial x_j} + {\partial v_j \over \partial x_i} \right)

(4)

In Eq.(4) Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \varepsilon _v}

is the volumetric strain rate defined as

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): \varepsilon _v= \varepsilon _{ii}

(5)

Substituting Eqs.(2) and (4) into (1) gives a useful form of the momentum equations as

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): \rho \frac{\partial v_i}{\partial t} - {\partial \over \partial x_j}(2\mu \varepsilon _{ij})+{\partial \over \partial x_i}\left(\frac{2}{3}\mu \varepsilon _v\right)- {\partial p \over \partial x_i}-b_i =0\quad , \quad i,j=1,n_s

(6)

Boundary conditions

The boundary conditions at the Dirichlet (Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \Gamma _v} ) and Neumann/traction (Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \Gamma _t} ) boundaries are

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): v_i -v_i^p =0 \qquad \hbox{on }\Gamma _v

(7.a)

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): \sigma _{ij}n_j -t_i^p =0 \qquad \hbox{on }\Gamma _t

(7.b)

where Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle v_i^p}

and Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle t_i^p}

, Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle i=1,n_s}

are the prescribed velocities and prescribed tractions on the Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \Gamma _v}
and Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \Gamma _t}
boundaries, respectively.

Mass balance equation

The mass balance equation for an incompressible fluid is written as

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): \varepsilon _v =0 \qquad \hbox{in }\Omega

(8)

Stabilized FEM techniques are needed for solving the general momentum equations in fluid mechanics when they are written in the Eulerian description. This is due to the effect of the convective acceleration terms that lead to loss of stability of standard Galerkin FEM. For Stokes flows (or for fluid flows formulated in the Lagrangian description) the convective terms vanish from the momentum equations and, consequently, the equations can be solved with the Galerkin FEM.

The problem, however, remains for obtaining stable solution for incompressible flows when an equal order interpolation is used for the velocities and the pressure, as it is the case for the element derived in the paper. This situation violates the so called inf-sup (or LBB) condition and, hence, there is a need to use stabilization techniques for solving the mass balance equation with the FEM [11,39].

3 STABILIZED FIC FORMS OF THE MASS BALANCE EQUATION

In our work the stabilized form of the mass balance equation is obtained using the Finite Calculus (FIC) procedure [26]–[31],[33]–[35]. We will use both the second order form of the mass balance equation for an incompressible fluid obtained using the FIC method in space, as well as the FIC form of the mass balance equation in time. These forms have the following expressions.

Second order FIC mass balance equation in space

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): \varepsilon _v + \frac{h_i^2}{12} \frac{\partial ^2 \varepsilon _v}{\partial x^2_i}=0\qquad \hbox{in }\Omega \qquad i=1,n_s

(9)

FIC mass balance equation in time

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): \varepsilon _v + \frac{\delta }{2} \frac{\partial \varepsilon _v}{\partial t}=0 \qquad \hbox{in }\Omega

(10)

Eq.(9) is obtained by expressing the balance of mass in a rectangular domain of finite size (with space dimensions Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle h_1\times h_2}

for 2D problems) and retaining higher order terms in the Taylor series expansions than those typically used in the infinitesimal theory for expressing the change of mass along the sides of the balance domain. The characteristic lengths Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle h_i}
are related to the finite element sizes in the discretized problem, as it will be explained later.

Eq.(10), on the other hand, is obtained by expressing the balance of mass in a domain of infinitesimal size in space and of finite dimension in time, where Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \delta }

is a characteristic time value.

In Eqs.(9) and (10) the terms involving Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle h_i}

and Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \delta }
play the role of stabilization terms respectively. The form of the Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle h_i}
and Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \delta }
parameters will be defined later. Note that for Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle h_i\to 0}
and Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \delta \to 0}
the standard expression of the mass balance equation (8) is recovered in all cases.

The derivation of Eqs.(9) and (10) for a 1D problem are shown in [37].

4 FIC FORM OF THE MASS BALANCE EQUATION IN TERMS OF THE MOMENTUM EQUATIONS

We will derive next a more useful FIC form of the stabilized mass balance equation expressed in terms of the momentum equations.

From the momentum equations (6) we obtain (neglecting the space changes of the viscosity Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \mu }

in the term involving Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \varepsilon _v}

)

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): \frac{2}{3}\mu \frac{\partial \varepsilon _v}{\partial x_i} = - \rho \frac{\partial v_i}{\partial t} + {\partial \over \partial x_j} (2\mu \varepsilon _{ij}) + {\partial p \over \partial x_i}+b_i = -\rho {\partial v_i \over \partial t} + \hat r_{m_i} \qquad i,j=1,n_s

(11)

From the last equation we deduce

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): \frac{\partial \varepsilon _v}{\partial x_i} = \frac{3}{2\mu } \left[-\rho \frac{\partial v_i}{\partial t} + \hat r_{m_i}\right]\qquad i=1,n_s

(12)

In the last two equations Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \hat r_{m_i}}

is a static momentum term defined as

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): \hat r_{m_i}= {\partial \hat \sigma _{ij} \over \partial x_j}+b_i \quad \hbox{with}\quad \hat \sigma _{ij}= 2\mu \varepsilon _{ij}+p \delta _{ij}

(13)

Let us introduce Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \frac{\partial \varepsilon _v}{\partial x_i}}

from Eq.(12) intro Eq.(9). This gives, after small algebra

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): \varepsilon _v + \frac{h_i^2}{12}\frac{\partial }{\partial x_i} \left(\frac{\partial \varepsilon _v}{\partial x_i}\right)= \rho \varepsilon _v + \frac{h_i^2}{8 \mu } \frac{\partial }{\partial x_i} \left(-\rho \frac{\partial v_i}{\partial t}+\hat r_{m_i}\right)\qquad i=1,n_s

(14)

Observing the term involving the time derivative of Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle v_i}

in Eq.(14) gives

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): \frac{\partial }{\partial x_i} \left(-\rho \frac{\partial v_i}{\partial t}\right)=-\rho \frac{\partial }{\partial t} \left({\partial v_i \over \partial x_i}\right)=-\rho \frac{\partial \varepsilon _v}{\partial t}

(15)

Substituting Eq.(15) into (14) gives

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): \varepsilon _v + \frac{h_i^2}{8 \mu }\left(-\rho \frac{\partial \varepsilon _v}{\partial t}+ {\partial \hat r_{m_i} \over \partial x_i} \right)=0

(16)

On the other hand, from Eq.(10) we deduce

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): -\frac{\partial \varepsilon _v}{\partial t}= \frac{2}{\delta } \varepsilon _v

(17)

Substituting Eq.(17) into (16) gives

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): \varepsilon _v + \frac{h_i^2}{8\mu }\left(\frac{2\rho }{\delta } \varepsilon _v + {\partial \hat r_{m_i} \over \partial x_i} \right)=0

(18)

In the following we will assume Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle h_i = h}

where Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle h}
is a characteristic length that will be related to a typical average dimension of each finite element in the mesh. Multiplying Eq.(18) by Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle  \frac{8\mu }{h^2}}
gives, after grouping some terms,

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): \displaystyle \varepsilon _v +\tau {\partial \hat r_{m_i} \over \partial x_i} =0

(19)

where Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \tau }

is a stabilization parameter given by

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): \tau = \left(\frac{8\mu }{h^2}+ \frac{2\rho }{\delta } \right)^{-1}

(20)

Eq.(19) is the FIC form of the stabilized mass balance equation. This equation will be taken as the starting point for deriving the stabilized FIC-FEM formulation. Note again that for Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \tau =0} , the standard form of the incompressibility condition (8) is obtained.

5 VARIATIONAL EQUATIONS

5.1 Variational expression of the momentum equation

Multiplying Eq.(1) by arbitrary test functions Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle w_i ({x})}

(with dimension of velocity) and integrating over the analysis domain Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \Omega }
gives the standard weighted residual form of the momentum equations as

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): \int _\Omega w_i \left(\rho {\partial v_i \over \partial t}- {\partial \sigma _{ij} \over \partial x_j}-b_i\right)d\Omega =0

(21)

Integrating by parts the term involving Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \sigma _{ij}}

and making use of the Neumann  boundary conditions (7.b), yields the standard principle of virtual power as [4,39]

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): \int _\Omega w_i \rho {\partial v_i \over \partial t} d\Omega + \int _\Omega \delta \varepsilon _{ij} \sigma _{ij} d\Omega - \int _\Omega w_i b_i d\Omega - \int _{\Gamma _t} w_i t_i^p d\Gamma =0

(22)

where Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \delta \varepsilon _{ij}= \frac{1}{2}\left({\partial w_i \over \partial x_j}+{\partial w_j \over \partial x_i}\right)}

is an arbitrary (virtual) strain rate field.

Substituting the expression of the stresses from Eq.(2) into (22) gives

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): \int _\Omega \! w_i \rho {\partial v_i \over \partial t} d\Omega + \!\int _\Omega \left[\delta \varepsilon _{ij} 2\mu \left(\!\varepsilon _{ij} - \frac{1}{3}\varepsilon _{kk} \delta _{ij}\!\right)+ \delta \varepsilon _{v}p\right]d\Omega - \!\int _\Omega w_i b_i d\Omega - \!\int _{\Gamma _t} w_i t_i^p d\Gamma =0

(23)

Eq.(23) can be written in matrix form as

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): \displaystyle \int _\Omega {\bf w}^T \rho {\partial {\bf v} \over \partial t} d\Omega + \!\int _\Omega \delta {\boldsymbol \varepsilon }^T {\bf D} {\boldsymbol \varepsilon } d\Omega + \!\int _\Omega \delta {\boldsymbol \varepsilon }^T {\bf m} p - \int _\Omega {\bf w}^T {b}d\Omega - \!\int _{\Gamma _t} \mathbf{\bf w}^T \mathbf{\bf t}^p d\Gamma =0

(24)

In Eq.(24) Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle {\bf w},{\bf v},{\boldsymbol \varepsilon }}

are vectors containing the weighting functions, the velocities and the strain rates, respectively; Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle {\bf b}}
and Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \mathbf{\bf t}^p}
are body force and surface tractions vectors, respectively; Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle {\bf D}}
is the viscous constitutive matrix and Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle {\bf m}}
is an auxiliary vector. These vectors are defined as (for 3D problems)

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): \begin{array}{l}\mathbf{w} =[w_1,w_2,w_3]^T\quad ,\quad \mathbf{v} =[v_1,v_2,v_3]^T\\[.5cm] {\boldsymbol \varepsilon }= [\varepsilon _{11},\varepsilon _{22}, \varepsilon _{33},2\varepsilon _{12},2\varepsilon _{13}, 2\varepsilon _{23}]^T \quad ,\quad \mathbf{\bf b} =[b_1,b_2,b_3]^T \quad ,\quad \mathbf{\bf t}^p =[t_1^p,t_2^p,t_3^p]^T\\[.5cm] {\bf D}=\mu \left[ \begin{array}{cccccc}\frac{4}{3} & -\frac{2}{3} &-\frac{2}{3} & 0 & 0 & 0 \\ & \frac{4}{3} & -\frac{2}{3} &0 &0 & 0 \\ & & \frac{4}{3} & 0 & 0 & 0 \\ & & & 1 & 0 & 0 \\ {Sym.} & & & & 1 & 0 \\ & & & & & 1 \\ \end{array} \right]\qquad ,\qquad {\bf m}=[1,1,1,0,0,0]^T \end{array}

(25)

The 2D form of above expressions is straightforward. For instance, the 2D expression of matrix Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle {\bf D}}

is obtained by deleting the rows and columns 3, 5 and 6 in Eq.(25).

Remark 1. From the definition of Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle {\bf m}}

and Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle {\bf D}}
and Eqs.(2), (3) and (5)  we deduce

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\boldsymbol \sigma } = {\bf s} + {\bf m}p \quad ,\quad {\bf s}= {\bf D}{\boldsymbol \varepsilon } \quad \hbox{and}\quad \varepsilon _v= {\bf m}^T {\boldsymbol \varepsilon }

(26)

where Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle {\boldsymbol \sigma } =[\sigma _{11},\sigma _{22},\sigma _{33}, \sigma _{12},\sigma _{13}, \sigma _{31}]^T,~{s}=[s_{11},s_{22},s_{33}, s_{12},s_{13}, s_{23}]^T}

are the stress and deviatoric stress vectors, respectively.

5.2 Variational expression of the stabilized mass balance equation

Multiplying equation (19) by arbitrary test functions Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle q(\mathbf{x})}

(with dimension of pressure) defined over the analysis domain Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \Omega }
and integrating over Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \Omega }
gives

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): \int _\Omega q \varepsilon _v d\Omega + \int _\Omega q \tau {\partial \hat r_{m_i} \over \partial x_i} d\Omega =0

(27)

Let us integrate by parts the second integral in Eq.(27) over a mesh of Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle N_e}

elements each one  with area Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \Omega ^e }
and boundary Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \Gamma ^e}

. This yields

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): \sum \limits _{e=1}^{N_e} \left\{\int _{\Omega ^e} q \varepsilon _v d\Omega - \int _{\Omega ^e} \tau {\partial q \over \partial x_i} \hat{r}_{m_{i}} d \Omega +\int _{\Gamma ^e} q \tau n_i \hat{r}_{m_{i}} d \Gamma \right\}=0

(28)

where Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle n_i}

are the components of the unit normal vector to the element boundary Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \Gamma ^e}

. Note that Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \Gamma ^e}

is a side for a 3-noded triangle and a face for a 4-noded tetrahedron.

Remark 2. In Eq.(28), the second integral over the domain Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \Omega ^e}

involving Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle {\partial q \over \partial x_i}}
is zero for a piecewise constant approximation of the weighting function Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle q}
over the  elements. This situation occurs for the P1/P0+ elements presented in this work.

From the definition of Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \hat{r}_{m_{i}}}

of Eqs.(13) and (1) we deduce

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): \hat{r}_{m_{i}} = \rho {\partial v_i \over \partial t}+ \frac{2\mu }{3}{\partial \varepsilon _v \over \partial x_i}

(29)

Multiplying by the normal components gives

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): \hat{r}_{m_{i}}n_i = \rho {\partial v_n \over \partial t}+\frac{2\mu }{3}{\partial \varepsilon _v \over \partial n}

(30)

where Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle v_n =v_in_i}

is the projection of the velocity in the direction of the normal to the element boundary and Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle {\partial \varepsilon _v \over \partial n}={\partial \varepsilon _v \over \partial x_i}n_i}
is the derivative of Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \varepsilon _v}
in the normal direction.

The term Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \mu {\partial \varepsilon _v \over \partial n}}

at an element boundary Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \Gamma _{ea}}
connecting elements Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle e}
and Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle a}
can be computed as

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): \frac{2}{3}\mu {\partial \varepsilon _v \over \partial n} = \frac{2}{l_{ea}} \left(\frac{2}{3}\mu ^+ \varepsilon _v^+ - \frac{2}{3}\mu ^- \varepsilon _v^- \right)

(31)

where Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle (\cdot )^+}

and Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle (\cdot )^-}
denote the values of the relevant magnitude at the external and internal points adjacent to the element boundary, respectively and Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle l_{ea}}
is the characteristic length of the boundary. For triangles Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle l_{ea}}
is taken as the side lenght  (Figure 1). For 4-noded tetrahedra Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle l_{ea} = 2 (\Omega _{ea})^{1/2}}
where Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \Omega _{ea}}
is the area of the triangular face connecting elements Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle e}
and Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle a}

.

(a)	(b)

(c)

Figure 1: (a) Patch of four triangles associated to a central 3-noded triangle Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): e . (b) Definition of Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\partial \varepsilon _v \over \partial n} at the side Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): ij of element Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): e . (c) Equilibrium of tractions at the side Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): ij adjacent to elements Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): e and Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): a . (d) Jump of Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): \frac{2}{3}\mu \varepsilon _v across a side of element Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): e

Eq.(31) is just a simple (yet effective) model for computing Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \mu {\partial \varepsilon _v \over \partial n}}

at the element boundary. Indeed, more sophisticated procedures can be used [9].

From the Neumann boundary conditions (7b) and Eq.(31) we deduce

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): \frac{2}{3}\mu ^+ \varepsilon _v^+ = 2\mu ^+ {\partial v_n^+ \over \partial n} + p^+ - t^+_n \quad \hbox{ at }\Gamma _{ea}^+	(32.a)
Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): \frac{2}{3}\mu ^- \varepsilon _v^- = 2 \mu ^- {\partial v_n^- \over \partial n} + p^- - t^-_n \quad \hbox{ at }\Gamma _{ea}^-	(32.b)

where indices Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle +}

and – denote values at the external and internal sides Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \Gamma _{ea}^+}
and Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \Gamma _{ea}^-}
of the element boundary, respectively. In Eqs.(5.2) Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle t_n^+}
and Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle t_n^-}
are the normal tractions respectively acting at Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \Gamma _{ea}^+}
and Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \Gamma _{ea}^-}
(Figure 1c).

Substituting Eqs.(5.2) into (31) gives

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): \frac{2}{3}\mu {\partial \varepsilon _v \over \partial n} = \frac{2}{l_{ea}} \left[\!\!\left[2\mu {\partial v_n \over \partial n}+p -t_n \right]\!\!\right]_{\Gamma ^e}

(33)

where Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \left[\!\left[a \right]\!\right]}

denotes the jump of the magnitude Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle a}
across the element boundary Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \Gamma ^e}

, i.e. Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \left[\!\left[a \right]\!\right]= a^+-a^-} .

From the equilibrium of tractions at the element boundaries we obtain

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): \left[\!\left[t_n \right]\!\right]= t_n^+ - t_n^- = t_n^p = \left[{t}^p \right]^T {n}^+

(34)

where Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle t_n^p}

is the external normal traction acting on the element side in the direction of the  normal to Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \Gamma ^+_{ea}}
(Figure 1). Clearly for unloaded boundaries Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle t_n^p=0}
and, consequently, Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \left[\!\left[t_n \right]\!\right]=0}

.

We note that above derivations are independent of the choice of Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \Gamma _{ea}^+}

and Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \Gamma _{ea}^-}

.

Substituting Eq.(33) into (30) and this one into (28) gives the final stabilized variational form of the mass balance equation as

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): \sum \limits _{e=1}^{N_e} \left\{\int _{\Omega ^e} q \varepsilon _v d\Omega - \int _{\Omega ^e} \tau {\partial q \over \partial x_i} \hat r_{m_i} d\Omega + \int _{\Gamma ^e} \frac{2\tau }{l^e}q \left(\frac{\rho l^e}{2} {\partial v_n \over \partial t} +\left[\!\!\left[2\mu {\partial v_n \over \partial n}+p -t_n \right]\!\!\right]\right)d\Gamma \right\}=0

(35)

where Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle l^e}

are the lengths of the element sides.

For element boundaries laying on the Dirichlet boundary Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \Gamma _v} , Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \hat r_{m_i} n_i =0}

and the boundary term vanishes from Eq.(35). This is justified by the fact that both Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle {\partial v_n \over \partial t}}
and Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle {\partial \varepsilon _v \over \partial n}}
are zero at Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \Gamma _v}

.

Remark 3. At element boundaries laying on an external Neumann boundary Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \Gamma _t} , Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle t^+_n=0} , Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle {n}^+=-{n}^- = {n}^e}

and, hence, Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \left[\!\left[t_n \right]\!\right]=-t_n^p}
at Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \Gamma _t}

.

Remark 4. The boundary integral in Eq.(35) resembles the jump stabilization terms across element sides proposed by different authors for fluid flow problems [3,9,12,14,17,23]. It is remarkable that Eq.(35) emanates naturally from the FIC formulation. Also, the form of Eq.(35) introduces the effect of the acceleration in the normal direction to the element side, as well as the effect of the external surface tractions acting on the element side. These two terms have proven to be important for the improved mass conservation and overall stability of the numerical solution of the multifluid problems solved in this paper and also for free surface homogeneous fluid flows [37].

Remark 5. For steady state problems the acceleration term Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \rho {\partial v_n \over \partial t}}

vanishes from the boundary integral in Eq.(35).

6 FEM DISCRETIZATION

We discretize the analysis domain into a mesh of Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle N_e}

finite elements in the standard manner. In our work we will choose simple 3-noded triangles (for 2D problems) and 4-noded tetrahedra (for 3D problems). The velocity is linearly interpolated in terms of the nodal values, while a discontinuous linear interpolation for the pressure over each element is chosen, i.e.

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): v_i = \sum \limits _{j=1}^n N_j \bar v_i^j \quad ,\quad p = \bar p^e + \tilde p({\bf x})

(36)

where Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle N_i}

are the standard linear shape functions for simplicial elements, Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle n}
is the number of element nodes (Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle n=3/4}
for 2D/3D problems), Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \bar{v}_i^j}
denotes the value of the Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle i}

th velocity component for the Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle j} th node of an element, Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \bar p^e}

is a constant pressure field over each element and Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle  \tilde {p}({\bf x})}
is a discontinuous linear pressure field chosen as (Figure 2)

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): \tilde {p}({\bf x}) = ({\bf x}_{c} - {\bf x})^{T} {\bf b}

(37)

where Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle {x}_{c}}

are the coordinates of the element midpoint.

Figure 2: One dimensional representation of the discontinuous pressure field. (a) Constant element pressure (Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): \bar p^e

) and linear pressure field (Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): \tilde p ). (b) Pressure field Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): p=\bar p^e + \tilde p

over the element. (c) Discontinuous linear pressure field over a three 1D element patch

Note that Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \tilde {p}({x})}

satisfies the steady state form of the momentum equation for a linear velocity field and  constant body forces Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle b_i}

, i.e.

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\partial \tilde p \over \partial x_i}+b_i=0

(38)

The actual unknowns of the problem are the nodal velocities Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \bar v_i^j}

and the constant pressure field Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \bar p^e}
for each element (Figure 3).

We will choose now Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle N_e}

piecewise constant unit weighting functions Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle q^e ({\bf x})}
 so that Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle q^e ({\bf x})=1}
and Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle  {\partial q^e \over \partial x_i}=0 }
if Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle  {x} \in \Omega ^e}
and Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle q^e ({\bf x}) =0}
if Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle  {x} \not \in \Omega ^e}

. With these assumptions the variational form (35) for the stabilized mass balance equation simplifies for 3-noded triangles to

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): \int _{\Omega ^e} \varepsilon _v d\Omega + \sum \limits _{i=1}^{3} 2 \tau \left( \frac{\rho l_i^e}{2} {\partial v_n \over \partial t}+ \left[\!\!\left[2\mu {\partial v_n \over \partial n}+p -t_n\right]\!\!\right] \right) =0 \quad ,\quad e=1,N_e

(39)

where the sum in the second term extends over the three sides of each triangular element and Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle l_i^e}

is the length of the Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle i}

th element side.

Figure 3: Nodal velocities and element pressure variables at a patch of four 3-noded triangles

Remark 6. For 4-noded tetrahedra the sum in Eq.(39) extends over the four faces of the tetrahedra and Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle l_i^e}

is a characteristic distance for the Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle i}

th face with area Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \Omega _i^e}

computed as Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle l_i^e = 2(\Omega _i^e)^{1/2}}
where Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \Omega _i^e}
is the area of the face.

Substituting the approximation (36) into Eqs.(39) and (24) and choosing a Galerkin form with Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle w_j =N_j}

gives the discretized expression of the momentum and (stabilized) mass balance equations as

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\bf M} {\dot{\bar{\bf v}}} + {\bf K}\bar{\bf v}+{\bf Q}\bar {\bf p}- {\bf f}_v={\bf 0}

(40.a)

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\bf Q}^T \bar{\bf v} + {\bf S}\bar{\bf p}-{\bf f}_p={\bf 0}

(40.b)

with

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): \bar{\bf v} = \left\{\begin{matrix}\bar{\bf v}^1\\\bar{\bf v}^2\\\vdots \\ \bar{\bf v}^N\end{matrix}\right\}~~,~~ \bar{\bf v}^j = [\bar{\bf v}^j_1,\bar{\bf v}^j_2,\bar{\bf v}^j_3]^T \quad \hbox{and } \bar{\bf p} = [\bar{p}^1,\bar{p}^2,\cdots , \bar{p}^{N_e}]^{T}

(40.c)

where Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle N}

is the number of nodes in the mesh. The different matrices and vectors in Eqs.(6a,b) are assembled from the element contributions given in Box 1 for 3-noded triangles.

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): \begin{array}{l} \displaystyle \textbf{M}_{ij}^e =\int_{\Omega^e}\rho \textbf{N}_i^T \textbf{N}_jd\Omega ~,~\textbf{K}^e_{ij} =\int_{\Omega^e} \textbf{B}_i^T \textbf{D}\textbf{B}_j d\Omega ~,~\textbf{Q}^e_{ij} =\int_{\Omega^e} \textbf{B}_i^T \textbf{m}d\Omega\\[.4cm] \displaystyle {S}_{ee} = - 2 \left(\tau^{ea} + \tau^{eb}+\tau^{ec}\right)\quad ,\quad {S}_{eb} = 2\tau^{eb}~,~{S}_{ec} = 2\tau^{ec}~,~{S}_{ea} = 2\tau^{ea}\\[.4cm] \displaystyle{\bf B}_i = \left[ \begin{matrix} \displaystyle\frac{N_i}{x_1} &0\\ 0& \displaystyle\frac{N_i}{x_2}\\ \displaystyle\frac{N_i}{x_2}&\displaystyle\frac{N_i}{x_1}\\[.25cm] \end{matrix} \right]\quad , \quad {\bf N}_i = \left[ \begin{array}{ccc} N_i & 0 \\ 0 & N_i \\ \end{array} \right] \quad \quad , \quad {\bf m} = [1,1,0]^T\\\\ \displaystyle \textbf{f}^e_{v_i}= \int_{\Omega^e}\textbf{N}_i^T \textbf{b} d\Omega + \int_{\Gamma_t} \textbf{N}_i^T {\bf t}d\Gamma\\[.4cm] \displaystyle f_{p_e}=-\sum\limits_{r=a,b,c} 2\tau^{er} \left( \frac{\rho l_{er}}{2} \frac{v_n}{t}+ \left[\!\!\left[ 2\mu \frac{v_n}{n}+ ({\bf x}_c^T - {\bf x}_{er}^T) {\bf b}-{t}_n \right]\!\!\right]\right) \end{array}

where Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): l_{er}

is the length of the side connecting elements e and r and Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\bf x}_{er}
 is the mid-point of the side. For continuous body forces between elements Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): \left[\!\left[{\bf x}_{er}^T {\bf b} \right]\!\right]=0

The expression of Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): \tau^{er}

is given in Eqs.(45) and (46).

Box 1. Element matrices and vectors in Eqs.(\ref{eq36}) for 3-noded triangles (Figure 1)

Remark 7. The FEM approximation (36) yields the strain rates and the stresses within an element in terms of the nodal velocities and pressure as

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): \begin{array}{c}{\boldsymbol \varepsilon } = \sum \limits _{j=1}^n \mathbf{B}_j \bar {\bf v}^j \quad ,\quad \varepsilon _v = \mathbf{m}^T \sum \limits _{j=1}^n \mathbf{B}_j \bar {\bf v}^j\\ \displaystyle \mathbf{s} = \mathbf{D} \sum \limits _{j=1}^n \mathbf{B}_j \bar {\bf v}^j \quad ,\quad {\boldsymbol \sigma } = \mathbf{D} \sum \limits _{j=1}^n \mathbf{B}_j \bar {\bf v}^j + \mathbf{m} \bar p^e\end{array}

(41)

where Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \mathbf{B}_j}

is given in Box 1 for 2D problems.

7 SOLUTION OF THE DISCRETIZED EQUATIONS

The discretized form in time of the system of Eqs.(6) is expressed as

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\bf M} \frac{{}^{n+1}\bar{\bf v}-{}^{n}\bar{\bf v}}{\Delta t} + {\bf K} {}^{n+1}\bar{\bf v}+{\bf Q}{}^{n+1}\bar {\bf p}-{\bf f}_v={\bf 0}

(42.a)

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\bf Q}^T {}^{n+1}\bar{\bf v}+ {\bf S} {}^{n+1}\bar{\bf p}-{\bf f}_p={\bf 0}

(42.b)

A symmetric monolithic form of Eq.(7) can be written as

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): \begin{bmatrix}\frac{1}{\Delta t}{\bf M}+ {\bf K}& {\bf Q} \\[.2cm] {\bf Q}^T & {\bf S} \end{bmatrix} \left\{\begin{matrix}{}^{n+1}\bar{\bf v} \\[.2cm]{}^{n+1}\bar {\bf p}\\[.2cm]\end{matrix}\right\}=\left\{\begin{matrix}{\bf f}_v+ \frac{1}{\Delta t} {\bf M} {}^{n}\bar{\bf v}\\[.2cm] {\bf f}_p \end{matrix}\right\}

(43)

The steady state form of Eq.(43) is simply

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): \begin{bmatrix}{\bf K}&{\bf Q}\\[.2cm]{\bf Q}^T & {\bf S} \end{bmatrix} \left\{\begin{matrix}\bar{v} \\[.2cm]\bar {\bf p} \end{matrix}\right\}= \left\{\begin{matrix}{\bf f}_v\\[.2cm]{\bf f}_p \end{matrix}\right\}

(44)

Note that the terms involving the normal velocity to the element, the normal tractions to the side and the term emanating form the side discontinuous linear pressure field Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \tilde p}

have been incorporated into the force vector Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle {\bf f}_p}
in Eq.(40.b). An obvious alternative will be to add the contribution of these terms to matrices Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \bf Q}
and Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \bf S}
in the l.h.s. of this equation.

Remark 8. The stabilization parameter Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \tau }

of Eq.(20) is computed at each element boundary Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle ij}
connecting elements Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle i}
and Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle j}
as

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): \tau = \tau ^{ij}= \left(\frac{8\mu _{ij}}{l^2_{ij}}+\frac{2\rho _{ij}}{\Delta t} \right)^{-1}

(45.a)

where Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle l_{ij}}

is a characteristic length of the element boundary. For triangles Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle l_{ij}}
is taken as the side length. For 4-noded tetrahedra Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle l_{ij}=2 (\Omega _{ij})^{1/2}}

, where Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \Omega _{ij}}

is the area of the face connecting elements Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle i}
and Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle j}

. The material parameters Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \mu _{ij}}

and Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \rho _{ij}}
are computed as

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): \mu _{ij} = \frac{1}{2} (\mu _i + \mu _j) \quad ,\quad \rho _{ij}= \frac{1}{2} (\rho _i + \rho _j)

(45.b)

where indices Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle i}

and Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle j}
denote the element where the material parameter is computed.

Remark 9. For the steady state case the stabilization parameter of Eq.(45.a) is computed as

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): \tau ^{ij}= \left(\frac{8\mu _{ij}}{l^2_{ij}}+\frac{2\rho _{ij} \vert{v}_{ij}\vert }{l_{ij}} \right)^{-1}

(46)

with Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle {v}_{ij}}

being the velocity vector of the mid-point of the boundary connecting elements Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle i}
and Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle j}

.

The above expressions of the stabilization parameter are similar to those used by many stabilization methods [11,39].

8 EXTENSION TO THE NAVIER-STOKES EQUATIONS

The momentum equations for Navier-Stokes flows are written in the Eulerian framework as

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): r_{m_i}:= \rho \left(\frac{\partial v_i}{\partial t}+v_j{\partial v_{i} \over \partial x_j} \right)- {\partial \sigma _{ij} \over \partial x_j}-b_i=0

(47)

Note that the difference with Eq.(1) is that the convective acceleration terms are now accounted for in Eq.(47).

The constitutive equation and the boundary conditions are identical to Eqs.(2)–(4) and (7), respectively.

The mass balance equation is given by Eq.(8). The stabilized form is given by Eq.(19).

The FEM solution of the Navier-Stokes equation for incompressible fluids requires a stabilized numerical method that can capture the internal layers introduced by the convective acceleration terms, as well as for satisfying the inf-sup condition introduced by the incompressibility constraint. Different stabilization procedures have been proposed in the past two decades. An overview can be found in [11,39]. In our work we use the FIC approach for obtaining stabilized FEM solutions for the Navier-Stokes equations [26]–[31],[33]–[35].

In the FIC method the momentum equations are derived in a balance domain of finite size. This yields a non-local form of the equations. For the Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle i} th momentum equation we have

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): r_{m_i} - \frac{h_{ij}}{2} \frac{\partial r_{m_i}}{\partial x_j}=0 \quad i,j=1,n_s \hbox{ (sum in }j\hbox{ only)}

(48)

where Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle h_{ij}}

are characteristic distances that define the balance domain [28]–[30].

For consistency reasons, in the FIC method the Neumann boundary conditions expressing balance of tractions at the boundary Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \Gamma _t}

are also derived in a finite domain adjacent to the boundary. The modified Neumann boundary conditions are written as

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): \sigma _{ij}n_h - t_i^p + \frac{h_{in}}{2} r_{m_i} =0 \quad i,j=1,n_s \quad \hbox{at } \Gamma _t

(49)

where Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle h_{in} = h_{ij} n_j} .

Applying the standard weighted residual procedure to Eq.(48) we obtain

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): \int _\Omega w_i r_{m_i} d\Omega - \int _\Omega w_i \frac{h_{ij}}{2} \frac{\partial r_{m_i}}{\partial x_j}d\Omega=0

(50)

where Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle w_i}

are standard weighting functions.

Integrating by parts the second integral in (50) assuming that Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle w_i=0}

on Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \Gamma _v}
and using the modified Neumann boundary conditions (49) we obtain

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): \int _\Omega w_i r_{m_i} d\Omega + \sum \limits _{e=1}^{N_e} \int _{\Omega ^e} \frac{h_{ij}}{2} \frac{\partial w_i}{\partial x_j}r_{m_i}d\Omega=0

(51)

The characteristic distances can be defined in a number of ways. For instance choosing

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): h_{ij}=h^e \frac{v_i}{\vert {v}\vert }

(52)

where Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle h^e}

is a characteristic distance for element Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle e}

, reproduces the standard SUPG method [11,39].

The effect of sharp internal gradients in the solution can be introduced by defining Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle h_{ij}}

as

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): h_{ij}=h^e_v\frac{v_i}{\vert {v}\vert } + h_g^e \frac{1}{\vert {v}\vert } {\partial v_{i} \over \partial x_j}

(53)

The first term in the r.h.s. of Eq.(53) introduces the SUPG stabilization while the second term introduces a shock-capturing type of stabilization. In Eq.(53) Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle h^e_v}

and Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle h_g^e}
are appropriate characteristic distances for 1D elements Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle h_v^e = h_g^e = l^e}

. The simplest choice for simplicial 2D and 3D elements is Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle h^e_v=h_g^e= 2 [\Omega ^e]^{1/n_s}} .

As for the treatment of the mass balance equation, the same procedure explained in Section 5.2 has been applied. Hence, the stabilized mass balance expression is identical to Eq.(19). Also, the corresponding variational form of the mass balance equation coincides with Eq.(35).

The final system of discretized equations using P1/P0+ elements is identical to Eq.(6). For 3-noded triangles, all matrices and vectors coincide with the expressions of Box 1 except Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle {K}}

and Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle {f}_v}
which are now given by (for each element)

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\bf K}^e ={\bf K}^e_0 + {\bf K}^e_v + {\bf K}^e_s \quad ,\quad {\bf f}^e_v ={\bf f}^e_{v_0} + {\bf f}^e_s

(54)

In Eq.(54) Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle {\bf K}^e_0}

is the viscous contribution  coinciding with the expression of Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle {\bf K}^e}
of Box 1 and Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle  {\bf K}^e_v}
and Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle  {\bf K}^e_s}
are the convective and stabilization contributions for the element given by

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\mathbf{K}}_{v_{ij}}^e = \int _{\Omega ^e} \rho {N}_i ({v}^T {\boldsymbol \nabla } N_j){I}_{ns} d\Omega \quad ,\quad {K}^e_{s_{ij}}= \int _{\Omega ^e} \frac{1}{2} \mathbf{N}_i \mathbf{G}^T \mathbf{H} (\mathbf{\nabla } {v}^T)^T \mathbf{N}_jd\Omega

(55.a)

where Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle {\bf I}_{ns}}

is the Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle n_s \times n_s}
unit matrix and

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): \displaystyle{\bf G} = \left[\begin{matrix}{\boldsymbol \nabla } & {0} & {0}\\ {0}&{\boldsymbol \nabla } & {0}\\ {0} & {0}& {\boldsymbol \nabla } \end{matrix} \right]~,~ \mathbf{H} = \left[\begin{matrix}{\bf h}_1 & {\bf 0}& {\bf 0}\\ {\bf 0}& {\bf h}_2 & {\bf 0}\\ {\bf 0} & {\bf 0}& {\bf h}_3 \end{matrix} \right]~,~ \displaystyle {\bf h}_i = \left\{\begin{matrix}h_{i_1}\\h_{i_2}\\h_{i_3} \end{matrix} \right\}~,~ {\boldsymbol \nabla } = \left\{\begin{matrix}\displaystyle {\partial \over \partial x_1}\\ \displaystyle {\partial \over \partial x_2}\\ \displaystyle {\partial \over \partial x_3} \end{matrix} \right\}

(55.b)

In Eq.(54) Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle {\bf f}^e_{v_0}}

coincides with the expression of Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle {\bf f}^e_v}
in Box 1 and Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle  {\bf f}^e_s}
is the contribution from the stabilization terms given by

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\bf f}^e_{s_i} = \int _{\Omega ^e} \frac{1}{2} \mathbf{N}_i^T \mathbf{G}^T \mathbf{H} \mathbf{b} d\Omega

(56)

The general form of Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle h_{ij}}

of Eq.(53) can be used in (56) to introduce SUPG and/or shock capturing effects in the stabilized formulation, as appropriate.

For the problems we are solving in this work the velocity field is free from internal layers. Therefore, good results are obtained with the standard Galerkin form of the momentum equations neglecting the stabilization term. The full stabilized formulation expressed by Eq.(51), with Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle h_{ij}}

given by Eq.(53), should however be used for the general solution of Navier-Stokes problems involving high gradients of the velocity.

9 PARTICULARIZATION FOR LAGRANGIAN FLOWS

The formulation presented in Sections 2–6 for Stokes flows in an Eulerian framework can be easily particularized for analysis of a viscous fluid flow using a Lagrangian description of the motion [4].

The relevant change is the definition of the acceleration term in the momentum equations. In the updated Lagrangian formulation these are written as [4]

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): r_{m_i} =\rho \frac{Dv_i}{Dt}- {\partial \sigma _{ij} \over \partial {}^{n+1} x_j}- {}^{n+1}b_i =0\quad i,j = 1,n_s

(57)

where Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \frac{Dv_i}{Dt}}

is the material derivative of the Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle i}

th velocity component Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle v_i ({}^n {\bf x})}

of a material point with coordinates Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle {}^n {\bf x}}
at time Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle t=t_n}

. Also, the Cauchy stresses Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \sigma _{ij}} , the body forces Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle b_i}

and the coordinates Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle {}^{n+1}x_i}
are referred to the updated configuration at time Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle t_{n+1}}
[4].

The material derivative of the velocity in the Lagrangian formulation is typically approximated as

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): \frac{Dv_i}{Dt} = \frac{{}^{n+1}v_i({}^{n+1}{\bf x})-{}^nv_i({}^n{\bf x})}{\Delta t}

(58)

where Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle {}^nv_i({}^n{\bf x})}

denotes the value of the Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle i}

th velocity component of the material point Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle {}^n {x}}

at time Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle t_n}

, etc.

Remark 10. The material derivative in an Eulerian description is expressed in a fixed control domain as

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): \frac{Dv_i}{Dt} = {\partial v_i \over \partial t} +v_j {\partial v_i \over \partial x_j}

(59)

Eq.(59) is the form used in the Eulerian formulation of the Navier-Stokes equations (see Eq.(47)). Neglecting the convective acceleration term in Eq.(59) yields the expression used in Eq.(1) for Stokes flows.

A difficulty in the Lagrangian formulation of flow problems is the need for tracking the motion of the material points in space and time. The authors have developed in recent years a particular class of Lagrangian method for analysis of fluid flow problems termed the Particle Finite Element Method (PFEM, www.cimne.com/pfem). The PFEM treats the nodes in a continuum as virtual particles that can freely move and even separate from the main domain representing, for instance, the effect of water drops in a splashing fluid or soil particles in an excavation problem. A finite element mesh connects the nodes defining the discretized domain where the governing equations are solved using the FEM. An advantage of the Lagrangian formulation in the PFEM is that the convective terms disappear from the fluid equations. The difficulty is however transferred to the problem of adequately (and efficiently) moving the mesh node using a fast mesh regeneration procedure at each time step. The theory and applications of the PFEM are reported in [1,8],[18]–[21],[25,32],[36]–[38].

The particularization of the formulation of the P1/P0+ triangle presented in Sections 2–6 to the Lagrangian analysis of viscous flows is straightforward and simply implies computing the acceleration term with the expression given by Eq.(58). The rest of the terms are identical to those given in those sections.

In this work we have applied the Lagrangian formulation of the P1/P0+ triangle to the study of multifluid problems with relatively small changes of the geometry of the fluid domain. The application of the P1/P0+ triangle and tetrahedra to a wider class of Lagrangian flow problems using the PFEM will be the subject of further work of the authors.

10 EXAMPLES

In this section we present the solutions obtained by the FIC method for four benchmark problems. The first two examples are steady state problems formulated in the Eulerian framework. The last two examples are transient problems formulated in the Lagrangian framework. In all these examples, the pressure solution has a strong discontinuity across the interface where the fluid viscosities are relevant and distinct. These pressure jumps occur either due to a relevant normal derivative of the velocity or in its absence due to a prescribed over-pressure force at the interface.


(a) Eulerian extrusion problem.	(b) Zhong's problem.
Figure 4: Description of the steady-state examples. (a) fixed-mesh (Eulerian) extrusion problem and (b) the Zhong's problem.

10.1 Fixed-mesh (Eulerian) extrusion problem

This example deals with a flow entering a rectangular 2D domain from the left boundary with a unit normal velocity and moving to the right where it finds an impermeable slip wall that deviates the flow upwards and downwards [20]. The 2D domain occupies the region: Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle (x,y) \in [0, 1] \times [-0.5, 0.5]}

(Figure 4a). The upper-half and lower-half of the domain are occupied by two immiscible fluids which have the same density (Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \rho _{1} = \rho _{2} = \rho }

), but with distinct viscosities (Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \mu _{1} \neq \mu _{2}} ). This problem has an analytical solution which can be expressed as follows.

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\bf v}(x,y) = \begin{bmatrix}1 - x \\ y \end{bmatrix} , \quad p(x,y) = \left\{ \begin{array}{ll}\rho \left[x - \dfrac{1}{2}(x^{2} + y^{2}) -\dfrac{7}{24}\right]+ (\mu _{1} - \mu _{2}) & \forall y > 0 \\ \rho \left[x - \dfrac{1}{2}(x^{2} + y^{2}) -\dfrac{7}{24}\right]- (\mu _{1} - \mu _{2}) & \forall y < 0 \end{array} \right.

(60)

As the pressure solution is known up to a constant, in Eq. (60) the later is chosen such that the pressure has a zero mean value over the domain. At the interface (Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle y = 0} ), due to a jump in the viscosities and a relevant directional derivative of the velocity along its normal, the pressure field has a jump across it.

Remark 11. This is a steady state example involving the solution of the Navier-Stokes equations solved in the Eulerian framework (Section 8). Generally, the presence of the convective term in the momentum equation triggers both global and local numerical instabilities which needs to be stabilized. However, these instabilities are manifested only when the solution of the continuous problem involves layers. In examples like the current one velocity field does not contain layers. Hence, in this problem, we can obtain a good numerical approximation to the velocity even without convective stabilization.


(a) Obtained solution Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): p^h structured mesh.	(b) Best approximation Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\cal P}^h_{L^2}p structured mesh.

(c) Obtained solution Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): p^h unstructured mesh.	(d) Best approximation Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\cal P}^h_{L^2}p unstructured mesh.
Figure 5: Fixed-mesh (Eulerian) extrusion problem. The obtained pressure solutions compared with the corresponding best approximations. The material properties are chosen as Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): (\rho _{1}, \mu _{1}) = (5,5) and Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): (\rho _{2}, \mu _{2}) = (5,1) .

Figure 5a illustrates the pressure solution Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle p^{h}}

obtained by the P1-P0+ triangle. The domain is discretized using a symmetric structured mesh of Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle 2 \times 30 \times 30}
elements. The material properties are chosen to be: Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \rho = 5}

, Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \mu _{1} = 5}

and Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \mu _{2} = 1}

. The velocity boundary conditions are chosen to be of Dirichlet (essential) type, taking values given by the expression in Eq. (60). The indeterminacy of the pressure solution is removed by imposing the zero-mean pressure condition [11]. The nonlinear convective terms in the Navier-Stokes equations are linearized using the Newton–Raphson method (Appendix A) and they converge in just two iterations for a tolerance of Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle 10^{-6}} .

Figure 5b illustrates the best approximation to the exact pressure solution Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle p}

given in Eq. (60) from the discrete pressure space Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle Q^{h}}
of piecewise constant pressure values over the elements of the considered mesh. The best approximation Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \mathcal{P}^{h}_{\hbox{L}^{2}}p}
is chosen from Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle Q^{h}}
using the Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \hbox{L}^{2}}
norm (Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \| \star \| := \sqrt{\int _{\Omega } (\star )^{2} \mbox{d}\Omega }}

) as the metric. Thus, Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \mathcal{P}^{h}_{\hbox{L}^{2}}p}

is the Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \hbox{L}^{2}}
projection of Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle p}
onto Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle Q^{h}}

. It can be obtained as follows: Find Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \mathcal{P}^{h}_{\hbox{L}^{2}}p \in Q^{h}} , such that

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): \int _{\Omega }{q^{h}(p-\mathcal{P}^{h}_{\hbox{L}^{2}}p) \mbox{d}\Omega } = 0 \forall q^{h} \in Q^{h} \Rightarrow {\|{p-\mathcal{P}^{h}_{\hbox{L}^{2}}p}\| } \leq {\|{p-p^{h}}\| } \forall p^{h} \in Q^{h}

(61)

Figure 5c illustrates the constant pressure solution Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle p^{h}}

(termed Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \bar p^e}
in Section 6) obtained by the new  P1-P0+ triangle and using an unstructured mesh. The unstructured mesh is generated by small random perturbations of the interior nodes (excluding the nodes on the interface) of the structured mesh followed by a Delaunay tessellation. Figure 5d illustrates the best approximation Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \mathcal{P}^{h}_{\hbox{L}^{2}}p}
obtained using the same mesh.

(a) Fixed-mesh (Eulerian) extrusion problem.

(b) Zhong's problem.

Figure 6: Error convergence in the Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): \hbox{L}^{2}

norm on mesh refinement for Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): p^{h}
and Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): (\boldsymbol \nabla  \cdot {v}^{h})

.

Figure 6a shows the convergence of the error Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle (p - p^{h})}

in the Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \hbox{L}^{2}}
norm with respect to mesh refinement. This error is compared with the best approximation error Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle (p - \mathcal{P}^{h}_{\hbox{L}^{2}}p)}

. To measure the convergence rate the domain is discretized by a sequence of symmetrical structured meshes of Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle 2 \times n \times n}

elements, using Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle n \in \{ 20,24,28,32,36,40,44,48,52,56\} }

. The relative error in the pressure solutions are measured as follows.

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): e^{h}_{p\hbox{L}^{2}} = \frac{\| p - p^{h}\| }{\| p\| } = \frac{\sqrt{\int _{\Omega } (p-p^{h})^{2} \mbox{d}\Omega }}{\sqrt{\int _{\Omega } p^{2} \mbox{d}\Omega }}

(62)

As expected for a piecewise constant pressure approximation, the convergence rate is found to be one for both solutions. Further, the error line of the numerical solution Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle p^{h}}

is found to be close to the error line of the best approximation Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \mathcal{P}^{h}_{\hbox{L}^{2}}p}

, indicating the good accuracy of the obtained solution. Additionally, in Figure 6a the convergence in the Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \hbox{L}^{2}}

norm of the divergence of the velocity is also shown. As the exact solution is solenoidal, this error is measured as follows.

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): \| \boldsymbol \nabla \cdot {\bf v}^{h} \| = \sqrt{\int _{\Omega } (\boldsymbol \nabla \cdot {\bf v}^{h})^{2} \mbox{d}\Omega }

(63)

where Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle {\bf v}^{h}}

is the velocity solution obtained with the P1/P0+ triangle. The convergence rate for Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \| \boldsymbol \nabla  \cdot {\bf v}^{h} \| }
is found to be close to Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle 1.5}

, which is more than the expected first-order rate. This super-convergence for Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \| \boldsymbol \nabla \cdot {\bf v}^{h} \| }

could be due to the smooth (linear, cf. Eq. (60)) solution profile for the velocity.

10.2 Zhong's problem: buoyant Stokes flow with a columnar viscosity structure

This example deals with the buoyancy-driven Stokes flow of a fluid with a columnar viscosity structure and was first presented in [40]. It has been previously used as a benchmark in [13,24]. It is an idealization of the mantle convection and the plate dynamics that occur in the subduction zones of the Earth. The temperature of the subducted plates are much colder than the ambient mantle which in turn leads to sharp lateral variations in the viscosity. The 2D domain (cf. Figure 4b) is a unit square: Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle (x,y) \in [0, 1] \times [0, 1]} . The left-half and right-half of the domain have viscosities given by Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \mu _{1} = 1}

and Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \mu _{2} = 10^6}

, respectively. The density and the acceleration due to gravity are specified as Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \rho (x,y) = - \cos (\pi x) \sin (\pi y)}

and Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle {\bf g} = (0,-1)}

, respectively. Free-slip boundary conditions are imposed everywhere on the domain boundary.

The evaluation ¹ of the analytical solution is described in [40] using the method of separation of variables (used earlier in [16]) and the propagator matrix techniques [15]. The constant to fix the pressure is chosen such that the latter has a zero mean value over the domain. Due to a jump in the viscosities at the viscosity boundary (Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle x = 0.5} ) and a non-zero directional derivative of the velocity along its normal, the pressure field acquires a jump across it. The expression for the pressure jump at the viscosity boundary can be written as follows.

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): [[p (x=0.5,y) ]] = \frac{2 |{g}| (\mu _{2}-\mu _{1})^2 \cosh ^2(\pi /2) \cos (\pi y)}{(\mu _{1}+\mu _{2})^2 \sinh ^2(\pi ) - \pi ^2 (\mu _{2}-\mu _{1})^2}

(64)


(a) Obtained solution Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): p^h structured mesh.	(b) Best approximation Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\cal P}^h_{L^2}p structured mesh.

(c) Obtained solution Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): p^h unstructured mesh.	(d) Best approximation Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\cal P}^h_{L^2}p unstructured mesh.
Figure 7: Zhong's problem: buoyancy-driven Stokes flow with a columnar viscosity structure in a square domain with slip walls. The density, the acceleration due to gravity and the viscosities are specified as Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): \rho = -\cos (\pi x) \sin (\pi y) , Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {g} = (0,-1) , Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): \mu _{1} = 1 and Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): \mu _{2} = 10^6 , respectively.

Figures 7a and 7c show the pressure solution Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle p^h}

obtained by the P1-P0+ triangle using a structured mesh (Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle 2 \times 30 \times 30}
elements) and an unstructured mesh (Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle 2 \times 30 \times 30}
elements), respectively. The corresponding best approximations in the Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \hbox{L}^{2}}
norm Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \mathcal{P}^{h}_{\hbox{L}^{2}}p}
are shown in Figures  7b and 7d, respectively.

Figure 6b shows the convergence of the error Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle (p - p^{h})}

in the Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \hbox{L}^{2}}
norm with respect to mesh refinement. This error is compared with the best approximation error Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle (p - \mathcal{P}^{h}_{\hbox{L}^{2}}p)}

. As expected for a piecewise constant pressure approximation, the convergence rate is found to be one for both solutions. Further, the error line for the numerical solution Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle p^{h}}

is found to be close to the error line of the best approximation Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \mathcal{P}^{h}_{\hbox{L}^{2}}p}

, indicating the good accuracy of the obtained solution. Additionally, in Figure 6b the convergence in the Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \hbox{L}^{2}}

norm of the divergence of the velocity is also shown.  As expected for a piecewise linear velocity approximation, a first-order convergence rate is found for Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \| \boldsymbol \nabla  \cdot {\bf v}^{h} \| }

.

(¹) The algebraic work involved to arrive at a closed form analytical expression of the solution is very tedious which includes the symbolic inversion of a Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): 4 \times 4

matrix. Thus, the analytical solution is often computed numerically for a specified problem data.

10.3 Moving-mesh (Lagrangian) extrusion problem

This is a transient example where the initial configuration is a rectangular 2D domain occupied by two immiscible fluids contained (due to the action of gravity) in a region bounded with sufficiently high walls (Figure 8a). Both the walls and the floor are considered to be slippery. The floor and the left wall are fixed while the right wall is pushed towards left at a constant prescribed velocity Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle V_{p}} . This is a free surface problem with an interface between the two immiscible fluids and the shape of the domain varies due to the action of the right wall. The numerical solution is sought using a Lagrangian formulation of the problem (Section 9).


(a) Lagrangian) extrusion problem.	(b) Lagrangian example with an unstable interface.
Figure 8: Description of the transient examples. (a) moving-mesh (Lagrangian) extrusion problem and (b) Lagrangian example with fixed-walls and an unstable interface: the interface initially has a serrated profile.


(a) initial state (t = 0); structured mesh.	(b) final state (t = 0); structured mesh.

(c) initial state (t = 0); unstructured mesh.	(d) final state (t = 0); unstructured mesh.
Figure 9: Moving-mesh (Lagrangian) extrusion problem: the initial (Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): t = 0 ) and final (Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): t = 2 ) configurations are independent to the choice of the mesh-type.

The fluid viscosities are chosen as Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \mu _{1} = 1}

and Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \mu _{2} = 10}

. The fluid densities are chosen as Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \rho _{1} = 1}

and Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \rho _{2} = 5}

. The acceleration due to gravity is taken as Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle [0,-10] \hbox{m/s}^{2}}

and the prescribed velocity Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle V_{p} = 0.1 \hbox{m/s}}

. The 2D domain has an initial length Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle l_{1} = 0.8 \hbox{m}}

and height Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle l_{2} = 0.4 \hbox{m}}

. The upper-half of the domain is occupied by the fluid with properties Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \rho _{1}, \mu _{1}} . For these data, the dynamics is dominated by the viscous term and hence in the exact solution the inertial effects can be neglected.

As a shorthand notation, we denote by Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle {}^{0}{\bf x}}

and Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle {}^{t}{\bf x}}
as the positions of an arbitrary fluid particle at times Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle t = 0}
and Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle t}

, respectively. Similar notation is used for the rest of the dependent variables in the Lagrangian description, e.g. Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle {}^{0}{\bf v}}

and Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle ^{t}{\bf v}}

. Using this notation and the expressions Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle L_{1} := (l_{1} - V_{p}t)}

and Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle L_{2} := (l_{1}l_{2}/L_{1})}

, the exact solution to this problem can be expressed as follows.

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {}^{t}{\bf x} = \begin{bmatrix}(L_{1}/l_{1}) & 0 \\ 0 & (L_{2}/l_{2}) \end{bmatrix} {}^{0}{\bf v},\quad {}^{t}{\bf v} := \frac{\mbox{d}}{\mbox{d}t}({\bf x}^{t}) = \frac{v_{p}}{L_{1}} \begin{bmatrix}-1 & 0 \\ 0 & 1 \end{bmatrix} {}^{t}{\bf v}

(65)

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {}^{t}{\bf a} := \frac{\mbox{d}}{\mbox{d}t}({}^{t}{\bf v}) = \left(\frac{V_{p}}{L_{1}}\right)^{2} \begin{bmatrix}0 & 0 \\ 0 & 2 \end{bmatrix} {}^{t}{\bf x},\quad {}^{t}\boldsymbol \nabla ({}^{t}{\bf v}) := \frac{\mbox{d}}{\mbox{d}{}^{t}{\bf x}}({}^{t}{\bf v}) = \frac{V_{p}}{L_{1}} \begin{bmatrix}-1 & 0 \\ 0 & 1 \end{bmatrix}

(66)

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): p ({\bf x},t) = \left\{ \begin{array}{ll}\rho _{1} g (L_{2}-{}^{t}Y) & \hbox{ for} L_{2} \geq ^{t}Y > \dfrac{L_{2}}{2} \\ \rho _{1} g\dfrac{L_{2}}{2} + 2(\mu _{2}-\mu _{1})\dfrac{V_{p}}{L_{1}} + \rho _{2}g\left(\dfrac{L_{2}}{2} - ^{t}Y\right)& \hbox{ for} \dfrac{L_{2}}{2} > ^{t}Y \geq 0 \end{array} \right.

(67)

where Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle ^{t}{\bf a}}

denotes the acceleration of the fluid particle and Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle ^{t}Y}
is the y-coordinate of the trajectory Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle ^{t}{x}}

. The constant in the pressure solution given by Eq. (67) is chosen such that it takes a value zero at the free surface. In the numerical solution, no boundary conditions are imposed for the velocity at the free surface. On the rest of the boundary, Dirichlet (essential) conditions are imposed for the velocity. As the velocity boundary conditions are not exclusively of the Dirichlet type, no further conditions need to be imposed for the pressure. For the sake of simplicity the discretized (Lagrangian) equations are linearized using the Picard (Fixed-point) method. The configuration Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle ^{t}{x}}

is updated at every iteration. The time integration is done using the Newmark method with the choice Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \theta = 1/2}
and Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \beta = 1/4}
(zero dissipation; second-order accuracy). The time increment is chosen as Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \Delta t = 0.1}

s (Appendix B).

Figure 9b illustrates the final configuration (at time Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle t = 2 \hbox{s}} ) obtained with the P1-P0+ triangle using a symmetric structured mesh of Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle 2 \times 24 \times 12}

elements. The initial configuration is shown in Figure 9a. Figures 9c and 9d show the initial and final configurations when using an unstructured mesh. For the initial domain area to be conserved the interface and the free surface should raise to a height of  ${\textstyle 0.2667{\hbox{m}}}$  and  ${\textstyle 0.5333{\hbox{m}}}$ , respectively. This is indeed observed in Figures 9b and 9d showing the good area conservation properties for the P1/P0+ triangle.


(a) Obtained solution $p^{2,h}$ ; structured mesh.	(b) Best approximation ${\cal {P}}_{L^{2}}^{h}p^{2}$ ; structured mesh.

(c) Obtained solution $p^{2,h}$ ; unstructured mesh.	(d) Best approximation ${\cal {P}}_{L^{2}}^{h}p^{2}$ ; unstructured mesh.
Figure 10: Moving-mesh (Lagrangian) extrusion problem. The obtained pressure solutions at time $t=2$ are compared with the corresponding best approximations. The material properties are chosen as $(\rho _{1},\mu _{1})=(1,1)$ and $(\rho _{2},\mu _{2})=(5,10)$ .

Figure 10a shows the pressure solution ${\textstyle p^{2,h}}$ obtained using the structured mesh. As the two fluids have distinct densities, it leads to a jump in both the pressure and its gradient at the interface. Figure 10b shows the ${\textstyle {\hbox{L}}^{2}}$ projection of the exact pressure ${\textstyle p^{2}}$ given in Eq. (67) onto the discrete pressure space ${\textstyle Q^{h}}$ and is denoted as ${\textstyle {\mathcal {P}}_{{\hbox{L}}^{2}}^{h}p^{2}}$ . For the sake of comparison, the constants in the discrete pressure solutions ${\textstyle p^{2,h}}$ and ${\textstyle {\mathcal {P}}_{{\hbox{L}}^{2}}^{h}p^{2}}$ are chosen a posteriori such that the minimum value on the free surface is zero. Likewise, Figures 10c and 10d show the pressure solutions ${\textstyle p^{2,h}}$ and ${\textstyle {\mathcal {P}}_{{\hbox{L}}^{2}}^{h}p^{2}}$ obtained on the unstructured mesh.


(a) Obtained solution $p^{2,h}$ ; structured mesh.	(b) Best approximation ${\cal {P}}_{L^{2}}^{h}p^{2}$ ; structured mesh.

(c) Obtained solution $p^{2,h}$ ; unstructured mesh.	(d) Best approximation ${\cal {P}}_{L^{2}}^{h}p^{2}$ ; unstructured mesh.
Figure 11: Moving-mesh (Lagrangian) extrusion problem: $p^{2,h}$ and ${\mathcal {P}}_{{\hbox{L}}^{2}}^{h}p^{2}$ projected onto the $x-z$ plane. The material properties are chosen as $(\rho _{1},\mu _{1})=(1,1)$ and $(\rho _{2},\mu _{2})=(5,10)$ .

Figure 11 presents the same results as shown in Figure 10 but with the solutions projected onto the ${\textstyle x-z}$ plane. In this view, the inter-element jumps in the obtained solution ${\textstyle p^{2,h}}$ can be better compared with those in the best approximation ${\textstyle {\mathcal {P}}_{{\hbox{L}}^{2}}^{h}p^{2}}$ . We see that on both meshes the obtained solution ${\textstyle p^{2,h}}$ is in good agreement with the best approximation ${\textstyle {\mathcal {P}}_{{\hbox{L}}^{2}}^{h}p^{2}}$ .


(a) Effective pressure $p^{2,h}$ ; structured mesh.	(b) Effective pressure $p^{2,h}$ ; unstructured mesh.
Figure 12: Moving-mesh (Lagrangian) extrusion problem: The effective pressure solutions at time $t=2$ which includes a prescribed linear field defined over each element. The material properties are chosen as $(\rho _{1},\mu _{1})=(1,1)$ and $(\rho _{2},\mu _{2})=(5,10)$ .

Figures 12a and 12b illustrate the effective pressure solution obtained with the P1/P0+ triangle using the structured and the unstructured meshes, respectively. Recall that the effective pressure solution includes a prescribed linear field defined over each element given by Eq.(37). A nearly identical figure is obtained by plotting the exact pressure solution given by Eq.(67), thus indicating the good accuracy of the method.

Finally, in Figure 13 we present the time evolution of the error in the area for all the cases discussed in this section. The relative error in the area is measured in three different ways

e_{\Omega }^{h}:={\dfrac {\displaystyle \sum _{e\in \Omega }({}^{0}\Omega _{e}-{}^{t}\Omega _{e})}{\displaystyle \sum _{e\in \Omega }{}^{0}\Omega _{e}}},\quad e_{\Omega _{1}}^{h}:={\dfrac {\displaystyle \sum _{e\in \Omega _{1}}({}^{0}\Omega _{e}-{}^{t}\Omega _{e})}{\displaystyle \sum _{e\in \Omega _{1}}{}^{0}\Omega _{e}}},\quad e_{\Omega _{2}}^{h}:={\dfrac {\displaystyle \sum _{e\in \Omega _{2}}({}^{0}\Omega _{e}-{}^{t}\Omega _{e})}{\displaystyle \sum _{e\in \Omega _{2}}{}^{0}\Omega _{e}}}

(68)

where ${\textstyle {}^{0}\Omega _{e}}$ and ${\textstyle {}^{t}\Omega _{e}}$ represent the initial and current (at time ${\textstyle t}$ ) areas of an element with index ${\textstyle e}$ . ${\textstyle \Omega _{1}}$ and ${\textstyle \Omega _{2}}$ denote the regions occupied by the two immiscible fluids, respectively and ${\textstyle \Omega :=\Omega _{1}\cup \Omega _{2}}$ is the domain. The indices ${\textstyle 1}$ and ${\textstyle 2}$ refer to the upper-half and lower-half regions of the domain. ${\textstyle e_{\Omega }^{h}}$ is the relative error in the total area of both fluids. ${\textstyle e_{\Omega _{1}}^{h}}$ and ${\textstyle e_{\Omega _{2}}^{h}}$ are the relative errors in the total area of the fluids occupying the regions ${\textstyle \Omega _{1}}$ and ${\textstyle \Omega _{2}}$ , respectively. For the simulated time period ( ${\textstyle 2}$ s), the maximum relative error is below ${\textstyle 3.25\times 10^{-4}}$ for ${\textstyle e_{\Omega }^{h}}$ , ${\textstyle e_{\Omega _{1}}^{h}}$ and ${\textstyle e_{\Omega _{2}}^{h}}$ .


(a) Area errors using a structured mesh.	(b) Area errors using an unstructured mesh.
Figure 13: Moving-mesh (Lagrangian) extrusion problem: time evolution of the area errors. The material properties are chosen as $(\rho _{1},\mu _{1})=(1,1)$ and $(\rho _{2},\mu _{2})=(5,10)$ .

10.4 Lagrangian example with fixed-walls and an unstable interface

Finally, we present another transient example where the two-fluid interface initially has a serrated profile (thus unstable) and becomes straight in due course. This configuration (Figure 8b) is obtained by meshing the domain with an unstructured mesh and identifying the elements that occupy the upper-half of the domain to represent a fluid with lower density and the remaining elements to represent the other fluid with a higher density.

To observe the interface evolution in a reasonable time period, fluid-pairs with low viscosities are chosen. Further, to avoid sloshing-like effects that arise due to the low viscosities, all the walls of the domain are fixed and are considered slippery. Thus, there are only buoyancy driven low-speed motions in the fluid due to its initial unstable configuration.

Due to the small viscosities and low-speed motions, the pressure has a negligible jump across the interface. In order to have a visible pressure jump, we prescribe a fictitious over-pressure force at the interface given by the expression ${\textstyle {t}^{\hbox{I}}=-\gamma {n}}$ . Here ${\textstyle {n}}$ represents the outward normal to the interface with respect to the heavier fluid and ${\textstyle \gamma }$ represents the prescribed jump in the pressure value across the interface.

This problem is solved considering two fluid-pairs. For the first fluid-pair, the material properties of the fluids on top and bottom are taken as ${\textstyle (\rho _{1},\mu _{1})=(1,10^{-6})}$ and ${\textstyle (\rho _{2},\mu _{2})=(5,10^{-5})}$ , respectively. At the interface, the over-pressure force is taken as ${\textstyle \gamma =5}$ . The second fluid-pair consists of sunflower oil on the top and water at the bottom. The material properties for sunflower oil and water at ${\textstyle 25^{\circ }{\hbox{C}}}$ are (in International units) ${\textstyle (\rho _{1},\mu _{1})=(919,4.9\times 10^{-2})}$ and ${\textstyle (\rho _{2},\mu _{2})=(997,8.9\times 10^{-4})}$ , respectively. At the oil–water interface, the over-pressure force is taken as ${\textstyle \gamma =1000}$ . The initial configuration is the same for both fluid-pairs, cf. figure 8b. The time integration is done by the Newmark method with the choice ${\textstyle \theta =1/2}$ and ${\textstyle \beta =1/4}$ (zero-dissipation; second-order accuracy). The time increment is chosen as ${\textstyle \Delta t=0.01}$ s. As no re-meshing is done, the simulation is stopped after ${\textstyle 100}$ time steps, i.e. a simulation time of ${\textstyle 1}$ s, to avoid numerical difficulties associated with severe mesh deformations.


(a) First fluid-pair interface length evolution.	(b) Oil-Water interface length evolution.

(c) First fluid-pair domain at t = 0.11s.	(d) Oil-Water domain at t = 0.55s.
Figure 14: Lagrangian example with fixed-walls and an unstable interface. The time evolution of the interface length and the domain when the former first attains a minimum are shown. The left column corresponds to a fluid-pair: $(\rho _{1},\mu _{1})=(1,10^{-6})$ , $(\rho _{2},\mu _{2})=(5,10^{-5})$ and $\gamma =5$ . The right column corresponds to oil–water interface: $(\rho _{1},\mu _{1})=(919,4.9\times 10^{-2})$ and $(\rho _{2},\mu _{2})=(997,8.9\times 10^{-4})$ and $\gamma =1000$ .


(a) Obtained solution $P^{0.11,h}$ ; 3D view.	(b) $P^{0.11,h}$ projected onto the $x-z$ plane.

(c) Obtained solution $P^{0.55,h}$ ; 3D view.	(d) $P^{0.55,h}$ projected onto the $x-z$ plane.
Figure 15: Lagrangian example with fixed-walls and an unstable interface. The obtained pressure solutions when the interface lengths first attain a minimum are shown. The top row corresponds to a fluid-pair: $(\rho _{1},\mu _{1})=(1,10^{-6})$ , $(\rho _{2},\mu _{2})=(5,10^{-5})$ and $\gamma =5$ . The bottom row corresponds to oil–water interface: $(\rho _{1},\mu _{1})=(919,4.9\times 10^{-2})$ and $(\rho _{2},\mu _{2})=(997,8.9\times 10^{-4})$ and $\gamma =1000$ .

Due to the small viscosities, the profile of the interface will have an oscillatory behavior (possibly damped due to momentum distribution in all directions) and a steady state can be expected in due course. We use the total length of the interface as a metric to measure the “distance” of the perturbed state of the interface from the stable horizontal profile. For an oscillating interface profile, the interface length will have an oscillatory behavior and will take a minimum value equal to ${\textstyle 0.8}$ (horizontal length of the domain) whenever the interface profile is horizontal.

Within the simulated time period, the P1-P0+ triangle is able to reproduce this oscillatory behavior of the interface for the first fluid-pair, cf. Figure 14a. The interface length first attains a minimum value of ${\textstyle 0.80014}$ in eleven time steps (at ${\textstyle t=0.11}$ s). The interface at this instant attains a horizontal profile as shown in Figure 14c. However, near the end of the simulation we observe an unphysical increase in the interface length. This could be due to the numerical difficulties associated to mesh-distortion (recall that no re-meshing is done here).

As the viscosities in the second fluid-pair (oil–water) are three orders of magnitude greater than those for the first fluid-pair, it is reasonable to expect a greater time scale for the former. Figure 14b shows that the oil–water interface length gradually reduces (i.e. the interface flattens) and takes a minimum value of ${\textstyle 0.80842}$ at ${\textstyle t=0.55}$ s. The interface at this instant attains a nearly horizontal profile as shown in Figure 14d. Past ${\textstyle t=0.55}$ s, the interface length gradually returns back to its initial value. However, the trend near the end of the simulation shows an unphysical increase in the interface length, which again could be due to the numerical difficulties associated to mesh-distortion.

Figure 15 shows the pressure solutions obtained for the considered fluid-pairs when their interfaces first attain a minimum length. We observe jumps in the pressure across the interface equivalent to the prescribed over-pressure forces, i.e. ${\textstyle \gamma =5}$ for the first fluid-pair (Figure 15a) and ${\textstyle \gamma =1000}$ for the oil–water fluid-pair (Figure 15c). Figures 15b and 15d present the same result as shown in Figures 15a and 15c, respectively, but with the solutions projected onto the ${\textstyle x-z}$ plane. In this view, all the inter-element jumps in the obtained pressure solutions are better seen.


(a) Signed relative error in areas.	(b) Absolute relative error (log-scale) in areas.

(c) Signed relative error in areas.	(d) Absolute relative error (log-scale) in areas.
Figure 16: Lagrangian example with fixed-walls and an unstable interface: time evolution of the area errors are shown. The top row corresponds to a fluid-pair: $(\rho _{1},\mu _{1})=(1,10^{-6})$ , $(\rho _{2},\mu _{2})=(5,10^{-5})$ and $\gamma =5$ . The bottom row corresponds to oil–water interface: $(\rho _{1},\mu _{1})=(919,4.9\times 10^{-2})$ and $(\rho _{2},\mu _{2})=(997,8.9\times 10^{-4})$ and $\gamma =1000$ .

The time evolution of the area-errors is shown in Figure 16 (both signed and absolute relative errors). The area-errors are measured in the same way as it is done earlier in section 10.3, i.e. using the expressions given in Eq. (68). For the simulated time period and the first fluid-pair, we find ${\textstyle |e_{\Omega _{1}}^{h}|<8\times 10^{-4}}$ , ${\textstyle |e_{\Omega _{2}}^{h}|<8\times 10^{-4}}$ and ${\textstyle |e_{\Omega }^{h}|<10^{-8}}$ . Likewise for the oil–water fluid-pair we find ${\textstyle |e_{\Omega _{1}}^{h}|<4\times 10^{-4}}$ , ${\textstyle |e_{\Omega _{2}}^{h}|<4\times 10^{-4}}$ and ${\textstyle |e_{\Omega }^{h}|<10^{-8}}$ . These results indicate the good area preservation property of the P1-P0+ triangle.

11 CONCLUDING REMARKS

We have presented a 3-noded triangle and a 4-noded tetrahedra with a continuous linear velocity and a discontinuous linear pressure field formed by the sum of an unknown constant pressure field and a prescribed linear field that satisfies the steady state momentum equations for a constant body force. In the so called P1/P0+ elements the “effective” pressure field is linear, although the unknown pressure field is piecewise constant within each element. The P1/P0+ elements have shown an excellent behaviour in terms of accuracy and mass conservation for incompressible viscous flow problems with discontinuous material properties formulated in either Eulerian or Lagrangian descriptions.

ACKNOWLEDGEMENTS

This research was supported by project HFLUIDS of the Ministerio de Educación y Ciencia of Spain and projects SAFECON and REALTIME of the European Research Council.

The authors are grateful to Prof. Shijie Zhong and Dr. Thibault Duretz for sharing their codes to evaluate the analytical solution of the buoyant Stokes flow example with a columnar viscosity structure.

APPENDIX A. LINEARIZATION OF THE NAVIER-STOKES EQUATIONS

Consider the steady Navier–Stokes equations described in an Eulerian framework. The stabilized momentum and mass balance equations discretized in space by the FEM lead to the following system of equations,

{\begin{bmatrix}(\mathbf {K} _{0}+\mathbf {K} _{v}+\mathbf {K} _{s})&\mathbf {Q} \\\mathbf {Q} ^{\hbox{T}}&\mathbf {S} \end{bmatrix}}{\begin{Bmatrix}{\bar {\bf {v}}}\\{\bar {\bf {p}}}\end{Bmatrix}}={\begin{Bmatrix}{\bf {f}}_{v}\\{\bf {f}}_{p}\end{Bmatrix}}

(A.1)

It is assumed that the stabilization terms in matrices ${\textstyle \mathbf {K} _{s}}$ and ${\textstyle \mathbf {S} }$ and vector ${\textstyle {\bf {f}}_{p}}$ are linear, i.e. they do not depend on the solution. The only nonlinear term that remains is, therefore, matrix ${\textstyle \mathbf {K} _{v}}$ obtained from the Galerkin FEM discretization of the convective term (Eq.(55a)).

The Newton–Raphson linearization of Eq.(1) gives the following system of equations.

{\begin{bmatrix}\mathbf {H} &\mathbf {Q} \\\mathbf {Q} ^{\hbox{T}}&\mathbf {S} \end{bmatrix}}{\begin{Bmatrix}\delta {\bar {\bf {v}}}\\\delta {\bar {\bf {p}}}\end{Bmatrix}}={\begin{Bmatrix}{\bf {f}}_{v}-(\mathbf {K} _{0}+\mathbf {K} _{v}+\mathbf {K} _{s}){\bar {\bf {v}}}^{m}-\mathbf {Q} {\bar {\bf {p}}}^{m}\\{\bf {f}}_{p}-\mathbf {Q} ^{\hbox{T}}{\bar {\bf {v}}}^{m}-\mathbf {S} _{p}{\bar {\bf {p}}}^{m}\end{Bmatrix}}

(A.2)

with

(A.3)

where matrices ${\textstyle \mathbf {K} _{0}^{e},~\mathbf {K} _{v}^{e}}$ and ${\textstyle \mathbf {K} _{s}^{e}}$ are given in Eqs.(54) and (8) and

(A.4)

In Eq.(2) index ${\textstyle m}$ denotes the iterations. Hence,

(A.5)

The Picard linearization of Eq.(1) yields the following system of equations

{\begin{bmatrix}(\mathbf {K} +\mathbf {K} _{v}^{m}+\mathbf {K} _{s})&\mathbf {Q} \\\mathbf {Q} ^{\hbox{T}}&\mathbf {S} \end{bmatrix}}{\begin{Bmatrix}{\bar {\bf {v}}}^{m+1}\\{\bar {\bf {p}}}^{m+1}\end{Bmatrix}}={\begin{Bmatrix}{\bf {f}}_{v}\\{\bf {f}}_{p}\end{Bmatrix}}

(A.6)

APPENDIX B. TIME INTEGRATION IN THE LAGRANGIAN FRAMEWORK

The momentum and the stabilized mass balance equations discretized in space by the FEM leads to the following system of equations in the Lagrangian framework

$\mathbf {M} {\bar {\bf {v}}}+\mathbf {K} {\bar {\bf {v}}}+\mathbf {Q} {\bar {\bf {p}}}={\bf {f}}_{v}$	(B.1)
$\mathbf {Q} ^{\hbox{T}}{\bar {\bf {v}}}+\mathbf {S} {\bar {\bf {p}}}={\bf {f}}_{p}$	(B.2)

where, the matrices and vectors with sub-scripts are the ones that appear due to the stabilization of the mass balance equation. Let ${\textstyle {\bar {a}}={\partial {\bar {v}} \over \partial t}}$ denote the vector of unknown nodal accelerations. The momentum balance equation, i.e. Eq.(B.1), is a dynamic equation. On the contrary, the stabilized mass balance equation, i.e. Eq.(B.2), is a quasi-static equation.

Using the Newmark algorithm the balance equations given in Eqs.(B.1) and (B.2) are integrated in time ( ${\textstyle \{{}^{n}\mathbf {x} ,{}^{n}\mathbf {v} ,{}^{n}\mathbf {a} ,{}^{n}{\bar {p}}\}\rightarrow \{{}^{n+1}{x},{}^{n+1}{\bar {v}},{}^{n+1}{\bar {a}},{}^{n+1}{\bar {p}}\}}$ ) as follows

${\begin{bmatrix}{\dfrac {1}{\Delta t}}\mathbf {M} +\theta ~{}^{n+1}\mathbf {K} &\theta ~{}^{n+1}\mathbf {Q} \\\theta \mathbf {Q} ^{\hbox{T}}&\theta ~{}^{n+1}\mathbf {S} \end{bmatrix}}{\begin{Bmatrix}{}^{n+1}{\bar {\bf {v}}}\\{}^{n+1}{\bar {\bf {p}}}\end{Bmatrix}}={\begin{Bmatrix}\theta ~{}^{n+1}{\bf {f}}_{v}+{}^{n}{\bf {r}}_{m}\\\theta ~{}^{n+1}{\bf {f}}_{p}\end{Bmatrix}}$
${}^{n}{\bf {r}}_{m}:={\frac {1}{\Delta t}}\mathbf {M} {}^{n}{\bar {\bf {v}}}-(1-\theta )\left[{}^{n}\mathbf {K} {}^{n}{\bar {\bf {v}}}+{}^{n}\mathbf {Q} {}^{n}{\bar {\bf {p}}}-{}^{n}{\bf {f}}_{v}\right]$	(6)
${}^{n+1}{\bar {\bf {a}}}={\frac {1}{\theta \Delta t}}\left({\bar {\bf {v}}}^{n+1}-{}^{n}{\bar {\bf {v}}}\right)-{\frac {(1-\theta )}{\theta }}{}^{n}{\bf {a}}$
${}^{n+1}{\bf {x}}={}^{n}{\bf {x}}+{\Delta t}{}^{n}{\bar {\bf {v}}}+{\frac {\Delta t^{2}}{2}}\left[(1-2\beta ){}^{n}{\bar {\bf {a}}}+2\beta {}^{n+1}{\bar {\bf {a}}}\right]$

REFERENCES

[1] R. Aubry, S.R. Idelsohn, E. Oñate, Particle finite element method in fluid mechanics including thermal convection–diffusion, Comput. Struct. 83 (2004) 1459–-75.

[2] R.F. Ausas , G.C. Buscaglia, S.R. Idelsohn, A new enrichment space for the treatment of discontinuous pressures in multi-fluid flows, Int. J. Numer. Meth. Fluids 70 (7) (2012) 829–850.

[3] R. Araya, G.R. Barrechea, F. Valentin, Stabilized finite element methods based on multiscale enrichment for the Stokes problem, Siam J. Numer. Anal. 44 (1) (2006) 322–-348.

[4] T. Belytschko, W.K. Liu, B. Moran, Non linear finite element for continua and structures, 2d Edition, Wiley, 2013.

[5] E. Burman, M. Fernández, P. Hansbo, Edge stabilizatio n for the incompressible Navier–Stokes equations: a continuous interior penalty finite element method, Tech. Report RR-5349, INRIA, Le Chesnay, France, 2004.

[6] E. Burman, P. Hansbo, A unified stabilized method for Stokes’ and Darcy equations, Tech. Report 2002-15, Chalmers Finite Element Center, Göteborg, Sweden, 2002.

[7] A. Caiazzo, M.A. Fernández, V. Martin, Analysis of a stabilized finite element method for fluid flows through a porous interface, Applied Mathematics Letters 24 (2011) 2124–2127.

[8] J.M. Carbonell, E. Oñate, B. Suárez, Modelling of tunnelling processes and cutting tool wear with the Particle Finite Element Method (PFEM), Accepted in Comput. Mech. (2013) DOI:10.1007/s00466-013-0835-x.

[9] R. Codina, S. Badia, On the design of discontinuous Galerkin methods for elliptic problems based on hybrid formulation, Comput. Methods Appl. Mech. Engrg. (2013) http://dx.doi.org/10.1016/j.cma.2013.05.004.

[10] F. Del Pin, S.R. Idelsohn, E. Oñate, R. Aubry, The ALE/Lagrangian Particle Finite Element Method: A new approach to computation of free-surface flows and fluid–object interactions, Computers & Fluids 36 (1) (2007) 27–38.

[11] J. Donea, A. Huerta, Finite element method for flow problems, J. Wiley & Sons, 2003, pp. 362.

[12] J. Douglas, J. Wang, An absolutely stabilized finite element method for the Stokes problem, Math. Comp. 52 (186) (1989) 495–508.

[13] T. Duretz, D. A. May, T. V. Gerya, P. J. Tackley, Discretization errors and free surface stabilization in the finite difference and marker-in-cell method for applied geodynamics: A numerical study, Geochemistry, Geophysics, Geosystems 12 (7) (2011).

[14] M.A. Fernández, J.-F. Gerbeau, V. Martin, Numerical simulation of blood flows through a porous interface, ESAIM: Mathematical Modelling and Numerical Analysis 42 (2008) 961–990.

[15] F. Gilbert, G. E. Backus, Propagator matrices in elastic wave and vibration problems, Geophysics 31 (2) (1966) 326–332.

[16] B. H. Hager, R. J. O'Connell, A simple global model of plate dynamics and mantle convection, Journal of Geophysical Research 86 (B6) (1981) 4843–4867.

[17] T.J.R. Hughes, L.P. Franca, A new finite element formulation for computational fluid dynamics: VII. The Stokes problem with various well-posed boundary conditions: Symmetric formulations that converge for all velocity/pressure spaces, Comput. Methods Appl. Mech. Engrg. 65 (1987) 85–-96.

[18] S.R. Idelsohn, E. Oñate, F. Del Pin, The particle finite element method: a powerful tool to solve incompressible flows with free-surfaces and breaking waves, Int. J. Numer. Methods Engrg. 61 (2004) 964–989.

[19] S.R. Idelsohn, E. Oñate, F. Del Pin, Fluid-structure interaction with the Particle Finite Element Method, Comput. Methods Appl. Mech. Engrg. 195 (2006) 2100–2113.

[20] S.R. Idelsohn, M. Mier-Torrecilla, E. Oñate, Multi-fluid flows with the Particle Finite Element Method, Comput Methods Appl Mech Engrg. 198 (2009) 2750–2767.

[21] S.R. Idelsohn, M. Mier-Torrecilla, R. Nigro, E. Oñate, On the analysis of heterogeneous fluids with jumps in the viscosity using a discontinuous pressure field, Comput. Mech. 46 (1) (2010) 115–124

[22] S.R. Idelsohn, E. Oñate, The challenge of mass conservation in the solution of free-surface flows with the fractional-step method: Problems and solutions, Int. J. Numer. Meth. Biomed. Engng. 26 (2010) 1313–-1330.

[23] N. Kechar, D. Silvester, Analysis of a locally stabilized mixed finite element method for the Stokes problem, Math. Comp. 58 (1992) 1–-10.

[24] M. Kronbichler, T. Heister, W. Bangerth, High accuracy mantle convection simulation through modern numerical methods, Geophysical Journal International 191 (1) (2012) 12–29.

[25] M. Mier-Torrecilla, S.R. Idelsohn, E. Oñate, Advances in the simulation of multi-fluid flows with the particle finite element method. Application to bubble dynamics, Int. J. Numer. Meth. Fluids 67 (2011) 1516–-1539.

[26] E. Oñate, Derivation of stabilized equations for advective-diffusive transport and fluid flow problems, Comput. Meth. Appl. Mech. Engng. 151 (1998) 233–267.

[27] E. Oñate, A stabilized finite element method for incompressible viscous flows using a finite increment calculus formulation, Comput Methods Appl Mech Engrg. 182 (1–2) (2000) 355–370.

[28] E. Oñate, J. García, A finite element method for fluid-structure interaction with surface waves using a finite calculus formulation, Comput. Meth. Appl. Mech. Engrg. 191 (2001) 635–660.

[29] E. Oñate, Multiscale computational analysis in mechanics using finite calculus: an introduction, Comput. Meth. Appl. Mech. Engrg. 192(28-30) (2003) 3043–3059.

[30] E. Oñate, Possibilities of finite calculus in computational mechanics, Int. J. Num. Meth. Engng. 60 (1) (2004) 255–281.

[31] E. Oñate, J. Rojek, R.L. Taylor, O.C. Zienkiewicz, Finite calculus formulation for incompressible solids using linear triangles and tetrahedra, Int. J. Numer. Meth. Engng. 59 (11) (2004a) 1473–1500.

[32] E. Oñate, S.R. Idelsohn, F. Del Pin, R. Aubry, The particle finite element method. An overview, Int. J. Comput. Methods 1 (2) (2004b) 267–307.

[33] E. Oñate, A. Valls, J. García, FIC/FEM formulation with matrix stabilizing terms for incompressible flows at low and high Reynold's numbers, Comput. Mech. 38 (4-5) (2006) 440–455.

[34] E. Oñate, J. García, S.R. Idelsohn, F. Del Pin, Finite calculus formulation for finite element analysis of incompressible flows. Eulerian, ALE and Lagrangian approaches, Comput. Meth. Appl. Mech. Engng. 195 (2006b) 3001–3037.

[35] E. Oñate, A. Valls, J. García, Computation of turbulent flows using a finite calculus-finite element formulation, Int. J. Numer. Meth. Engng. 54 (2007) 609–637.

[36] E. Oñate, M.A. Celigueta, S.R. Idelsohn, F. Salazar, B. Suárez, Possibilities of the particle finite element method for fluid-soil-structure interaction problems, Comput. Mech. 48(3) (2011) 307–318.

[37] E. Oñate, A. Franci, J.M. Carbonell, A FIC-based stabilized Lagrangian formulation with negligible mass losses for incompressible fluids. Application to free-surface flows using PFEM, Publication CIMNE, No. PI394. Submitted to Int. J. Num. Meth. Fluids, 2013.

[38] P. Ryzhakov, R. Rossi, S.R. Idelsohn, E. Oñate, A monolithic Lagrangian approach for fluid-structure interaction problems, Comput. Mech. 46 (6) 2010 883–899.

[39] O.C. Zienkiewicz, R.L. Taylor, P. Nithiarasu, The finite element method for fluid dynamics, Elsevier, 2006.

[40] S. Zhong. Analytic solutions for Stokes’ flow with lateral variations in viscosity, Geophysical Journal International 124 (1) 1996 18–28.

@@ Line 1,649: / Line 1,649: @@
 |}
-==References==
+==REFERENCES==
 <div id="cite-1"></div>