Revision as of 09:17, 7 June 2019 (view source) Cinmemj (talk | contribs) ← Older edit	Latest revision as of 09:17, 7 June 2019 (view source) Cinmemj (talk | contribs) m (Cinmemj moved page Draft Samper 970825501 to Onate Carbonell 2014a)
(No difference)

---

Latest revision as of 09:17, 7 June 2019

Published in Computational Mechanics, Vol. 54 (6), pp. 1583-1596, 2014
DOI: 10.1007/s00466-014-1078-1

Abstract

We present a mixed velocity-pressure finite element formulation for solving the updated Lagrangian equations for quasi and fully incompressible fluids. Details of the governing equations for the conservation of momentum and mass are given in both differential and variational form. The finite element interpolation uses an equal order approximation for the velocity and pressure unknowns. The procedure for obtaining stabilized FEM solutions is outlined. The solution in time of the discretized governing conservation equations using an incremental iterative segregated scheme is described. The linearization of these equations and the derivation of the corresponding tangent stiffness matrices is detailed. Other iterative schemes for the direct computation of the nodal velocities and pressures at the updated configuration are presented. The advantages and disadvantages of choosing the current or the updated configuration as the reference configuration in the Lagrangian formulation are discussed.

Keywords: Updated Lagrangian formulation, incompressible fluids, finite element method, incremental solution, tangent matrix, mixed formulation, stabilized method

1 INTRODUCTION

In Lagrangian analysis procedures, the motion of the particles of a deforming body is followed in time. In Eulerian formulations, on the other hand, attention is focused in the motion of the material through a stationary control volume. Lagrangian methods have been traditionally used for the numerical analysis of solids and structures, while Eulerian methods are typical in computational fluid dynamics [2,4,5,36,40,41].

Despite this evidence, the use of a Lagrangian description for solving fluid dynamics problems using the finite element method (FEM) [10,41] and different meshless and mesh-based particle-based numerical techniques [6,7,12],[16]–[20],[25,26] [29]–[35],[37,38] has received much attention in recent years. Of particular interest are numerical procedures, such as the Particle Finite Element Method (PFEM) [16,25,26,29,31], that combine the advantage of particle-based procedures with the formalism and accuracy of the FEM.

Despite the increasing interest in the Lagrangian approach for solving the equations of fluid mechanics using the FEM, there are few references to the derivation of incremental iterative solution schemes using linearized forms of the discretized Lagrangian equations for fluid flow problems. An early attempt in this direction was reported by Radovitzky and Ortiz [34] who derived the tangent matrix for the FEM Lagrangian analysis of compressible flows using a Newton-Rapsohn type iterative scheme.

The goal of this paper is to present a mixed velocity-pressure formulation for the finite element analysis of quasi and fully incompressible fluids using an updated Lagrangian formulation. We advocate an equal order interpolation for the velocity and pressure variables which invariably requires using an adequate stabilization scheme for the mass balance equation in order to obtain accurate numerical results. In the paper we derive in some detail both the continuum and discretized (FEM) forms of the equations governing the motion of the fluid in the updated Lagrangian description. An incremental Newton-Raphson type iterative staggered scheme for the solution in time of the non linear discretized equations is detailed. The derivations are carried out first using the current configuration as the reference configuration in the Lagrangian description of the motion. The expression of the different matrices and vectors involved in the incremental iterative scheme is given. The particular form of the FEM equations when the updated configuration is taken as the reference configuration is presented. The direct transient solution of the nodal velocities and pressures using monolithic and staggered schemes is also presented for completeness.

The chapter concludes with a discussion of the computational advantages and disadvantages of choosing the current or the the updated configuration as the reference configuration in the Lagrangian description.

In the next section the basic concepts of the motion of a fluid are briefly revisited. These concepts are standard and can be found in many books on computational solid and fluid mechanics and fluid-structure interaction, among others [2,3,4,5,13,36]. This introductory section aims to set up the updated Lagrangian framework where the governing equations for the fluid are written and subsequently solved with the FEM using a mixed velocity-pressure formulation.

2 MOTION. DISPLACEMENT, VELOCITY AND ACCELERATION

We will consider a domain containing a fluid which evolves in time due to external and internal forces and prescribed velocities from an initial configuration at time ${\textstyle t={}^{0}t}$ (typically ${\textstyle {}^{0}t=0}$) to a current configuration at time ${\textstyle t={}^{n}t}$. The fluid volume ${\textstyle V}$ and its boundary ${\textstyle \Gamma }$ at the initial and current configurations are denoted as (${\textstyle {}^{0}V,{}^{0}\Gamma }$) and (${\textstyle {}^{n}V,{}^{n}\Gamma }$), respectively. The goal is to find the domain that the fluid occupies, as well as the velocities, strain rates and stresses (the deviatoric stresses and the pressure) in the so-called updated configuration at time ${\textstyle {}^{n+1}t={}^{n}t+\Delta t}$ (Figure 1). In the following a left superindex denotes the configuration of the fluid domain where the variable is computed.

Figure 1: Initial (${\displaystyle t={}^{0}t}$), current (${\displaystyle t={}^{n}t}$) and updated (${\displaystyle t={}^{n+1}t}$) configurations of a fluid domain ${\displaystyle V}$ with Neumann (${\displaystyle \Gamma _{t}}$) and Dirichlet (${\displaystyle \Gamma _{v}}$) boundaries. Trajectory of a material point ${\displaystyle i}$ and velocity (${\displaystyle {v}_{i}}$) and displacement (${\displaystyle {u}_{i}}$) vectors of the point at each configuration. ${\displaystyle {}^{n}{u}}$, ${\displaystyle {}^{n+1},{u}}$ and ${\displaystyle \Delta {u}}$ denote schematically the trajectories of the overall domain between two configurations.

The motion of the fluid domain is described by

${\displaystyle {}^{t}{x}={\boldsymbol {\phi }}({}^{0}{x},t)\quad {\hbox{or}}\quad {}^{t}x_{i}=\phi _{i}({}^{0}{x},t)}$

(1)

where ${\textstyle {}^{t}{x}}$ is the position of the material point ${\textstyle {}^{0}{x}}$ laying on the initial configuration at time ${\textstyle t}$. The coordinates ${\textstyle {}^{t}{x}}$ give the spatial position of a point and are called spatial or Eulerian coordinates. The function ${\textstyle {\boldsymbol {\phi }}}$ maps the initial configuration into the current configuration at time ${\textstyle t}$. The position of the points in the current and initial configurations at time ${\textstyle t={}^{n}t}$ and ${\textstyle t={}^{0}t}$, respectively are expressed by

${\displaystyle {\begin{array}{l}{}^{n}{x}={\boldsymbol {\phi }}({}^{0}{x},{}^{n}t)\quad {\hbox{or}}\quad {}^{n}x_{i}=\phi _{i}({}^{n}{x},{}^{n}t)\\{}^{0}{x}={\boldsymbol {\phi }}({}^{0}{x},{}^{0}t)\quad {\hbox{or}}\quad {}^{0}x_{i}=\phi _{i}({}^{0}{x},{}^{0}t)\end{array}}}$

(2)

Thus the mapping ${\textstyle {\boldsymbol {\phi }}({}^{0}{x},{}^{0}t)}$ is the identity mapping.

In the Lagrangian description (also called the material description) the independent variables are the material coordinates ${\textstyle {}^{0}{x}}$ of the point on the initial configuration and the time ${\textstyle t}$. In the Eulerian description, on the other hand, the independent coordinates are the spatial coordinates ${\textstyle {}^{t}{x}}$ and the time ${\textstyle t}$ [4,5].

The displacement of a material point is given by the difference between its position at time ${\textstyle t}$ and its original position, so

${\displaystyle {u}({}^{0}{x},t)={\boldsymbol {\phi }}({}^{0}{x},t)-{\boldsymbol {\phi }}({}^{0}{x},{}^{0}t)={}^{t}{x}-{}^{0}{x}}$

(3)

where ${\textstyle {u}=[u_{1},u_{2},u_{3}]^{T}}$ is the displacement vector of a point.

For ${\textstyle t={}^{n}t}$ we have

${\displaystyle {}^{n}{u}\equiv {u}({}^{0}{x},{}^{n}t)={}^{n}{x}-{}^{0}{x}}$

(4)

The velocity vector is the rate of change of the position vector for a material point, i.e. the time derivative of ${\textstyle \phi ({}^{0}{x},t)}$ with ${\textstyle {}^{0}{x}}$ held constant (also called the material time derivative or total derivative). The velocity vector is usually written as (using Eq.(3) and noting that the coordinates ${\textstyle {}^{0}{x}}$ are fixed)

${\displaystyle {}^{t}{v}=[{}^{t}v_{1},{}^{t}v_{2},{}^{t}v_{3}]^{T}\equiv {v}({}^{0}{x},t)={\frac {\partial {\boldsymbol {\phi }}({}^{0}{x},t)}{\partial t}}={\frac {\partial {u}({}^{0}{x},t)}{\partial t}}\equiv {\dot {u}}}$

(5)

The velocity vector of a material point in the current configuration is written as

${\displaystyle {}^{n}{v}=[{}^{n}v_{1},{}^{n}v_{2},{}^{n}v_{3}]^{T}\equiv {v}({}^{0}{x},{}^{n}t)}$

(6)

The acceleration vector is the rate of change of the velocity vector of a material point (i.e. the material time derivative of the velocity vector) and it can be written as

${\displaystyle {}^{t}{a}={a}({}^{0}{x},t)\equiv {\frac {\partial {v}({}^{0}{x},t)}{\partial t}}={\frac {\partial ^{2}{u}({}^{0}{x},t)}{\partial t^{2}}}={\ddot {u}}}$

(7)

and

${\displaystyle {}^{n}{a}=[{}^{n}a_{1},{}^{n}a_{2},{}^{n}a_{3}]^{T}={a}({}^{0}{x},{}^{n}t)}$

(8)

Eq.(7) is the material form of the acceleration. Note the difference of Eq.(7) with the expression of the acceleration written in the Eulerian description which involves the convective terms [2,4,5,13,36,41].

In the total Lagrangian description the various equations and variables are referred to the initial configuration which is taken as the reference configuration. In the updated Lagrangian description, either the current or the updated configurations can be taken as the reference configuration [2,4,5,13,36].

For simplicity, in the first part of this work the current configuration will be taken as the reference configuration for the derivation of the incremental FEM equations. The reason is that the current configuration remain fixed during the iterative solution process. The particularization for the case when the updated configuration is taken as the reference configuration is presented in the last part of the paper.

The displacement increment vector ${\textstyle \Delta {u}({}^{0}{x},t)=[\Delta u_{1},\Delta u_{2},\Delta u_{3}]^{T}}$ of a material point between the updated and the current configurations is defined as

${\displaystyle \Delta {u}\equiv \Delta {u}({}^{0}{x},t)={\boldsymbol {\phi }}({}^{0}{x},{}^{n+1}t)-{\boldsymbol {\phi }}({}^{0}{x},{}^{n}t)={}^{n+1}{x}-{}^{n}{x}}$

(9)

The displacement increment of a material point can be computed from the velocity as

${\displaystyle \Delta {u}=\int _{{}^{n}t}^{{}^{n+1}t}{}^{\tau }{v}d\tau }$

(10)

where ${\textstyle {}^{\tau }{v}}$ is the velocity vector of the points laying on the trajectory followed by the material point between times ${\textstyle {}^{n}t}$ and ${\textstyle {}^{n+1}t}$ (Figure 1). The integral of Eq.(10) can be approximated in a number of ways (see Remark 4).

3 MOMENTUM EQUATIONS AND BOUNDARY CONDITIONS

Let us assume that at time ${\textstyle {}^{n+1}t}$ the fluid occupies a volume ${\textstyle {}^{n+1}V}$ with a boundary ${\textstyle {}^{n+1}\Gamma }$. The fluid is subjected to body forces ${\textstyle {}^{n+1}b_{i}}$ acting over the volume ${\textstyle {}^{n+1}V}$ and surface tractions ${\textstyle {}^{n+1}t_{i}^{p}}$ acting on the part of the boundary termed as ${\textstyle {}^{n+1}\Gamma _{t}}$, where index ${\textstyle i}$ denotes the component of the force along the ${\textstyle i}$th cartesian axis (Figure 1).

The equations of internal equilibrium between the applied body forces and the Cauchy stresses ${\textstyle \sigma _{ij}}$ in the fluid are expressed by the momentum equations written in the updated configuration as [2,4,5,13,36]

${\displaystyle {}^{n+1}r_{m_{i}}:={}^{n+1}\rho ~{}^{n+1}a_{i}-{\partial \sigma _{ij} \over \partial x_{j}}-{}^{n+1}b_{i}=0\quad {\hbox{in }}{}^{n+1}V\quad ,\quad i,j=1,\cdots ,n_{s}}$

(11)

where ${\textstyle {}^{n+1}r_{m_{i}}}$ is the ${\textstyle i}$th momentum equation, ${\textstyle {}^{n+1}\rho }$, ${\textstyle {}^{n+1}v_{i}}$ and ${\textstyle {}^{n+1}a_{i}}$ are the fluid density and the ${\textstyle i}$th component of the velocity and the acceleration at time ${\textstyle {}^{n+1}t}$ and ${\textstyle n_{s}}$ is the number of space dimensions (${\textstyle n_{s}=3}$ for three dimensional (3D) problems). Note that ${\textstyle \sigma _{ij}}$ is always referred to the updated configuration, i.e. ${\textstyle \sigma _{ij}\equiv {}^{n+1}\sigma _{ij}}$.

In Eq.(11) and in the following, the standard summation convention for repeated indices is adopted, unless otherwise specified.

The boundary conditions are

${\displaystyle {}^{n+1}v_{i}-{}^{n+1}v_{i}^{p}=0\quad {\hbox{on }}{}^{n+1}\Gamma _{v}}$

(12)

${\displaystyle \sigma _{ij}n_{j}-{}^{n+1}t_{i}^{p}=0\quad {\hbox{on }}{}^{n+1}\Gamma _{t}}$

(13)

where ${\textstyle {}^{n+1}v_{i}^{p}}$ is the prescribed value of the ${\textstyle i}$th velocity component at the external boundary ${\textstyle {}^{n+1}\Gamma _{v}}$ with ${\textstyle {}^{n+1}\Gamma _{v}\cup {}^{n+1}\Gamma _{t}={}^{n+1}\Gamma }$. In Eq.(13) ${\textstyle n_{j}}$ is the ${\textstyle j}$th component of the unit normal to the boundary ${\textstyle {}^{n+1}\Gamma _{t}}$ where the prescribed surface tractions ${\textstyle {}^{n+1}t_{i}^{p}}$ are applied.

The goal is to obtain the geometry of the updated configuration ${\textstyle {}^{n+1}V}$, as well as the velocities and stresses in ${\textstyle {}^{n+1}V}$ that satisfy Eqs.(11)–(13).

4 PRINCIPLE OF VIRTUAL POWER IN THE UPDATED CONFIGURATION

The principle of virtual power (PVP) can be written in the updated configuration as [2,4,5]

${\displaystyle {}^{n+1}\delta A+{}^{n+1}\delta U-{}^{n+1}\delta W=0}$

(14)

where ${\textstyle {}^{n+1}\delta A}$, ${\textstyle {}^{n+1}\delta U}$ and ${\textstyle {}^{n+1}\delta W}$ are the virtual powers due to the acceleration terms, the stresses and the external loads acting on ${\textstyle {}^{n+1}V}$, respectively given by [2,4,5]

${\displaystyle {}^{n+1}\delta A=\int _{{}^{n+1}V}\delta {v}^{T}{~}^{n+1}\rho ~{}^{n+1}{a}~d~{}^{n+1}V}$ (15) ${\displaystyle {}^{n+1}\delta U=\int _{{}^{n+1}V}\left\{\delta {D}\right\}^{T}\left\{{\boldsymbol {\sigma }}\right\}d~{}^{n+1}V}$ (16) ${\displaystyle {}^{n+1}\delta W=\int _{{}^{n+1}V}\delta {v}^{T}{~}^{n+1}{b}~d{~}^{n+1}V+\int _{{}^{n+1}\Gamma _{t}}\delta {v}^{T}{~}^{n+1}{t}~d~{}^{n+1}\Gamma }$ (17)

In Eq.(15)–(16) ${\textstyle \delta {v}}$, ${\textstyle {b}}$ and ${\textstyle {t}}$ are the virtual velocity vector, the body forces vector and the surface tractions vector, respectively defined for 3D problems as

${\displaystyle \delta {v}=[\delta v_{1},\delta v_{2},\delta v_{3}]^{T}~~,~~{b}=[b_{1},b_{2},b_{3}]^{T}~~,~~{t}=[t_{1},t_{2},t_{3}]^{T}}$

(18)

In Eq.(16) ${\textstyle {\boldsymbol {\sigma }}}$ is the Cauchy stress tensor and ${\textstyle \delta {D}}$ is the virtual rate of deformation tensor defined as

${\displaystyle \delta {D}_{ij}={\frac {1}{2}}\left({\partial \delta v_{i} \over \partial {~}^{n+1}x_{j}}+{\partial \delta v_{j} \over \partial {~}^{n+1}x_{i}}\right)}$

(19)

where ${\textstyle \delta v_{i}}$ is the ${\textstyle i}$th component of the virtual velocity.

In the following ${\textstyle \left\{{A}\right\}}$, where A is a symmetric tensor, will denote the vector representation of A. Hence, in Eq.(16) ${\textstyle \left\{{\boldsymbol {\sigma }}\right\}}$ and ${\textstyle \left\{\delta {D}\right\}}$ are the Cauchy stress vector and the rate of deformation vector, respectively defined in Voigt notation [4,5] as

${\displaystyle \left\{{\boldsymbol {\sigma }}\right\}=\left[\sigma _{11},\sigma _{22},\sigma _{33},\sigma _{12},\sigma _{13},\sigma _{23}\right]^{T}}$ (20.a) ${\displaystyle \left\{\delta {D}\right\}=\left[\delta D_{11},\delta D_{22},\delta D_{33},\delta D_{12},\delta D_{13},\delta D_{23}\right]^{T}}$ (20.b)

Similarly, ${\textstyle \left[{A}\right]}$ will denote hereonwards the matrix form of tensor ${\textstyle A}$.

Remark 1. The PVP can be obtained by applying the standard weighted residual method to the governing equations (11) and (13) using the virtual velocities ${\textstyle \delta v_{i}}$ as weighting functions [4].

Remark 2. Tensor ${\textstyle \delta {D}}$ can be interpreted as the variation of the rate of deformation tensor ${\textstyle D}$ defined in terms of the velocities at the updated configuration as

${\displaystyle {D}_{ij}={\frac {1}{2}}\left({\partial {~}^{n+1}v_{i} \over \partial {~}^{n+1}x_{j}}+{\partial {~}^{n+1}v_{j} \over \partial {~}^{n+1}x_{i}}\right)}$

(21)

5 TRANSFORMATION TO THE CURRENT CONFIGURATION. LAGRANGIAN STRESS AND STRAIN MEASURES

In the following sections we will use the current configuration as the reference configuration for computing the internal virtual due to the acceleration terms and the stresses, as well as for subsequently performing the linearization of its discretized form via the FEM. The alternative of using the updated configuration as the reference configuration is presented in Section 12.

After the standard transformations of continuum mechanics [2,4,5,13,36] the internal virtual power at the updated configuration due to the acceleration terms and the stresses can be expressed in terms of material parameters, integrals, strain rates and stress measures evaluated at the current configuration ${\textstyle {}^{n}V}$ as

${\displaystyle {}^{n+1}\delta A=\int _{{}^{n}V}\delta {v}^{T}{~}^{n}\rho {~}^{n+1}{a}~d{~}^{n}V}$

(22.a)

${\displaystyle {}^{n+1}\delta U=\int _{{}^{n}V}\left\{\delta {\dot {E}}\right\}^{T}\left\{{S}\right\}d{}^{n}V}$

(22.b)

where ${\textstyle \delta {\dot {E}}}$ and ${\textstyle {S}}$ are the material virtual strain rate tensor and the second Piola-Kirchhoff stress tensor, respectively. The relationship between the material strain rate tensor ${\textstyle {\dot {E}}}$ and the rate of deformation tensor ${\textstyle {D}}$ and between the stress tensors ${\textstyle {\boldsymbol {\sigma }}}$ and ${\textstyle {S}}$ is [2,4,5,13,36]

${\displaystyle {\dot {E}}={F}^{T}{D}{F}\quad ,\quad {S}=J{F}^{-1}{\boldsymbol {\sigma }}{F}^{-T}}$

(23)

where F is the deformation gradient and ${\textstyle J}$ is the Jacobian, respectively defined as

${\displaystyle F_{ij}={\frac {\partial {~}^{n+1}{x}_{i}}{\partial {~}^{n}{x}_{j}}}\quad ,\quad J=\vert {F}\vert }$

(24)

From the expression of ${\textstyle {\dot {E}}}$ of Eq.(23) we can deduce

${\displaystyle \delta {\dot {E}}={F}^{T}\delta {D}{F}}$

(25)

Remark 3. The material strain rate tensor can also be obtained from the time derivative of the Green-Lagrange strain tensor ${\textstyle {E}}$ as [2,4,5]

${\displaystyle {\dot {E}}={\frac {d}{dt}}({E})={\frac {d}{dt}}\left({\frac {1}{2}}\left({C}-{I}\right)\right)={\frac {1}{2}}{\frac {d}{dt}}{C}={\frac {1}{2}}{\frac {d}{dt}}({F}^{T}{F})={\frac {1}{2}}({F}^{T}{\dot {F}}+{\dot {F}}^{T}{F})}$

(26)

where ${\textstyle {C}={F}^{T}{F}}$ is the right Cauchy-Green tensor and ${\textstyle {\dot {F}}={\frac {d}{dt}}({F})}$. The expression of ${\textstyle \delta {\dot {E}}}$ is obtained by writing the variation of ${\textstyle {\dot {E}}}$ with respect to the velocities as

${\displaystyle \delta {\dot {E}}={\frac {1}{2}}\left({F}^{T}\delta {\dot {F}}+\delta {\dot {F}}^{T}{F}\right)}$

(27)

A useful explicit expression of the material strain rate tensor components in terms of the velocities ${\textstyle {}^{n+1}v_{i}}$ and the displacement increments ${\textstyle \Delta u_{i}}$ is

${\displaystyle {\dot {E}}_{ij}={\frac {1}{2}}\left({}_{n}^{n+1}v_{i,j}+{}_{n}^{n+1}v_{j,i}+{}_{n}^{n+1}v_{k,i}~{}_{n}\Delta u_{k,j}+{}_{n}^{n+1}v_{k,j}~{}_{n}\Delta u_{k,i}\right)}$

(28.a)

with

${\displaystyle {}_{n}^{n+1}v_{i,j}={\partial {~}^{n+1}v_{i} \over \partial {~}^{n}x_{j}}\quad ,\quad ~{}_{n}\Delta u_{i,j}={\partial \Delta u_{i} \over \partial {~}^{n}x_{j}}}$

(28.b)

From Eqs.(5) we deduce

${\displaystyle \delta {\dot {E}}_{ij}={\frac {1}{2}}\left({}_{n}\delta v_{i,j}~+{}_{n}\delta v_{j,i}+{}_{n}\delta v_{k,i}~{}_{n}\Delta u_{k,j}+{}_{n}\delta v_{k,j}~{}_{n}\Delta u_{k,i}\right)}$

(29.a)

with

${\displaystyle {}_{n}\delta v_{i,j}={\frac {\partial \delta v_{i}}{\partial {}^{n}x_{j}}}}$

(29.b)

6 SPLIT OF THE INTERNAL VIRTUAL POWER INTO DEVIATORIC AND VOLUMETRIC COMPONENTS

The Cauchy stress tensor can be split in its deviatoric component ${\textstyle {\boldsymbol {\sigma }}'}$ and the hydrostatic pressure component ${\textstyle p}$ as

${\displaystyle {\boldsymbol {\sigma }}={\boldsymbol {\sigma }}'+p{I}_{3}}$

(30)

where ${\textstyle {I}_{3}}$ is the ${\textstyle 3\times 3}$ identity matrix.

Note that in incompressible fluid mechanics the pressure ${\textstyle p}$ is an independent variable. Also, unless otherwise specified we will assume that ${\textstyle p}$ is the pressure at the updated configuration (i.e. ${\textstyle p={~}^{n+1}p}$).

Substituting Eq.(30) into the internal virtual power expression in Eq.(16) gives

${\displaystyle {}^{n+1}\delta U=\int _{{}^{n+1}V}\left\{\delta {D}\right\}^{T}\left(\left\{{\boldsymbol {\sigma }}'\right\}+{m}p\right)d{~}^{n+1}V=\left(\int _{{}^{n+1}V}\left\{\delta {D}\right\}^{T}\left\{{\boldsymbol {\sigma }}'\right\}+\delta D_{v}p\right)d{~}^{n+1}V}$

(31)

where ${\textstyle D_{v}}$ is the volumetric strain rate given by

${\displaystyle D_{v}=D_{ii}={m}^{T}\left\{{D}\right\}\quad {\hbox{with}}\quad {m}=[1,1,1,0,0,0]^{T}}$

(32)

It can be proven that [4,5,13,36]

${\displaystyle D_{v}={\dot {E}}_{v}\quad {\hbox{with}}\quad {\dot {E}}_{v}=\left\{{\dot {E}}\right\}^{T}\left\{{C}^{-1}\right\}}$

(33)

where ${\textstyle {\dot {E}}_{v}}$ is the volumetric material strain rate,

${\displaystyle \left\{{E}'\right\}=\left[{\dot {E}}_{11},{\dot {E}}_{22},{\dot {E}}_{33},2{\dot {E}}_{12},2{\dot {E}}_{13},2{\dot {E}}_{23}\right]^{T}}$

(34.a)

${\displaystyle \left\{{C}^{-1}\right\}=\left[C_{11}^{-1},C_{22}^{-1},C_{33}^{-1},C_{12}^{-1},C_{13}^{-1},C_{23}^{-1}\right]^{T}}$

(34.b)

and ${\textstyle C_{ij}^{-1}}$ is the element ${\textstyle ij}$ of tensor ${\textstyle {C}^{-1}}$.

The internal virtual power expression of Eq.(31) can be written in the current configuration as follows.

Substituting the Cauchy stress split of Eq.(30) into the expression of the second Piola-Kirchhoff stress tensor of Eq.(23) gives

${\displaystyle {S}=J{F}^{-1}{\boldsymbol {\sigma }}'{F}^{-T}+pJ{C}^{-1}={S}'+pJ{C}^{-1}}$

(35)

where

${\displaystyle {S}'=J{F}^{-1}{\boldsymbol {\sigma }}'{F}^{-T}}$

(36)

is the deviatoric second Piola-Kirchhoff stress tensor. ${\textstyle {S}'}$ is usually called the (true) deviatoric component of ${\textstyle S}$ [4,5,13,36].

Substituting Eq.(35) into (22.b) gives

${\displaystyle {}^{n+1}\delta U=\int _{{}^{n}V}\left(\left\{\delta {\dot {E}}\right\}^{T}\left\{{S}'\right\}+J\delta {\dot {E}}_{v}p\right)d{}^{n}V}$

(37)

with

${\displaystyle \left\{{S}'\right\}{=[S'}_{11}{,S'}_{22}{,S'}_{33}{,S'}_{12}{,S'}_{13}{,S'}_{23}]^{T}\quad ,\quad \delta {\dot {E}}_{v}=\left\{\delta {\dot {E}}\right\}^{T}\left\{{C}^{-1}\right\}}$

(38)

Eq.(37) can also be obtained from Eq.(31) using the relationship between ${\textstyle {D},~{\boldsymbol {\sigma }}'}$ and ${\textstyle D_{v}}$ with ${\textstyle {\dot {E}}}$, ${\textstyle {S}'}$ and ${\textstyle {\dot {E}}_{v}}$, respectively and the expression ${\textstyle d^{n+1}V=Jd^{n}V}$ [2,4,5].

Substituting Eqs.(15), (22.a) and (37) into (14), the PVP in the updated configuration can be written in terms of the pressure, the deviatoric second Piola-Kirchhoff stresses and the virtual Green-Lagrange strains computed in the current configuration as

${\displaystyle {\begin{array}{r}\displaystyle \int _{{}^{n}V}\delta {v}^{T}{~}^{n}\rho {~}^{n+1}{a}d{~}^{n}V+\int _{{}^{n}V}\left(\left\{\delta {\dot {E}}\right\}^{T}\left\{{S}'\right\}+J\delta {\dot {E}}_{v}p\right)d{~}^{n}V-\\-\displaystyle \int _{{}^{n+1}V}\delta {v}^{T}{~}^{n+1}{b}~d{~}^{n+1}V-\int _{{}^{n+1}\Gamma _{t}}\delta {v}^{T}{~}^{n+1}{t}~d{~}^{n+1}\Gamma =0\end{array}}}$

(39)

The vector form of tensors ${\textstyle \delta {\dot {E}}}$ and ${\textstyle {S}}$ in Eq.(39) is

${\displaystyle \left\{\delta {\dot {E}}\right\}=[\delta {\dot {E}}_{11},\delta {\dot {E}}_{22},\delta {\dot {E}}_{33},2\delta {\dot {E}}_{12},2\delta {\dot {E}}_{13},2\delta {\dot {E}}_{23}]^{T}}$

(40.a)

${\displaystyle \left\{{S}\right\}=[S_{11},S_{22},S_{33},S_{12},S_{13},S_{23}]^{T}}$

(40.b)

Note that the contribution of the external forces in Eq.(39) is computed at the updated configuration and this requires using the correct expression for the density and the surface tractions on ${\textstyle {}^{n+1}V}$ (see also Remark 5).

7 CONSTITUTIVE RELATIONSHIP FOR THE DEVIATORIC STRESSES

The relationship between the deviatoric Cauchy stresses ${\textstyle {\sigma '}_{ij}}$ and the rates of deformation ${\textstyle D_{ij}}$ for a Newtonian fluid is written as

${\displaystyle {\sigma '}_{ij}={c}_{ijkl}{D'}_{kl}={c}_{ijkl}\left(D_{kl}-{\frac {1}{3}}D_{v}\delta _{kl}\right)}$

(41)

where ${\textstyle {D'}_{kl}}$ are the deviatoric rates of deformation. The rates of deformation ${\textstyle D_{ij}}$ are obtained in terms of the velocities by Eq.(21).

The components of the fourth order constitutive tensor ${\textstyle {c}}$ in the updated configuration are

${\displaystyle {c}_{ijkl}=\mu \left(\delta _{ik}\delta _{jl}+\delta _{il}\delta _{jk}\right)}$

(42)

where ${\textstyle \mu }$ is the viscosity of the fluid.

Eq.(41) can be written in vector form as

${\displaystyle \left\{{\boldsymbol {\sigma }}'\right\}=\left[{c}\right]\left\{{D}\right\}}$

(43)

with

${\displaystyle {\begin{array}{l}\left\{{\boldsymbol {\sigma }}'\right\}=\left[\sigma '_{11},\sigma '_{22},\sigma '_{33},\sigma '_{12},\sigma '_{13},\sigma '_{23}\right]^{T}\\\left\{{D}\right\}=\left[D_{11},D_{22},D_{33},2D_{12},2D_{13},2D_{23}\right]^{T}\end{array}}}$

(44)

and the constitutive matrix ${\textstyle [{c}]}$ is given by (for 3D problems)

${\displaystyle \left[{c}\right]=\mu \left[\displaystyle {\begin{matrix}\displaystyle {\frac {4}{3}}&\displaystyle -{\frac {2}{3}}&\displaystyle -{\frac {2}{3}}&0&0&0\\[.25cm]&\displaystyle {\frac {4}{3}}&\displaystyle -{\frac {2}{3}}&0&0&0\\[.25cm]&&\displaystyle {\frac {4}{3}}&0&0&0\\{\hbox{ Sym.}}&&&1&0&0\\&&&&1&0\\&&&&&1\end{matrix}}\right]}$

(45)

The constitutive equation (41) can be transformed to the current configuration to yield the relationship between the deviatoric second Piola-Kirchhoff stresses and the material strain rates as

${\displaystyle {S'}_{ij}={\cal {C}}_{ijkl}{\dot {E}}_{kl}}$

(46)

The constitutive tensor in the current configuration ${\textstyle \left[{\boldsymbol {\cal {C}}}\right]}$ is obtained from its counterpart in the updated configuration as [36]

${\displaystyle {\cal {C}}_{ijkl}=F_{Ai}^{-1}F_{Bj}^{-1}F_{Ck}^{-1}F_{Dl}^{-1}c_{ABCD}}$

(47)

The vector form of Eq.(46) is written as

${\displaystyle \left\{\mathbf {S} '\right\}=\left[{\boldsymbol {\cal {C}}}\right]\left\{{\dot {E}}\right\}}$

(48)

Matrix ${\textstyle \left[{\boldsymbol {\cal {C}}}\right]}$ is obtained by applying Voigt rule to the terms of tensor ${\textstyle {\boldsymbol {\cal {C}}}}$ [4,5].

8 THE MASS BALANCE EQUATION

The mass balance equation in the updated configuration is written for a quasi-incompressible fluid as

${\displaystyle {}^{n+1}r_{v}:=-{\frac {1}{{}^{n+1}\rho c^{2}}}{\frac {\partial p}{\partial t}}+D_{v}=0\quad {\hbox{in }}{}^{n+1}V}$

(49)

where ${\textstyle c}$ is the speed of sound in the fluid. For a fully incompressible fluid ${\textstyle c=\infty }$ and Eq.(49) reduces to the standard form ${\textstyle D_{v}=0}$. The quasi-incompressible form will be retained here and this will allow us to account for the effect of the (small) compressibility in most fluids.

Eq.(49) can be written in terms of the bulk modulus of the fluid ${\textstyle \kappa }$ defined as ${\textstyle \kappa =\rho c^{2}}$. For convenience we will retain the form of Eq.(49).

The variational form of the mass balance equation is obtained via the standard weighted residual method [41] as

${\displaystyle \int _{{}^{n+1}V}q\left(-{\frac {1}{{}^{n+1}\rho c^{2}}}{\frac {\partial p}{\partial t}}+D_{v}\right)d{~}^{n+1}V=0}$

(50)

where ${\textstyle q}$ are appropriate test functions with dimensions of pressure [4,5,41].

The integral expression (50) can be written in the current configuration using Eq.(33) as

${\displaystyle \displaystyle \int _{{}^{n}V}q\left(-{\frac {J^{2}}{{}^{n}\rho c^{2}}}{\frac {\partial p}{\partial t}}+J{\dot {E}}_{v}\right)d{~}^{n}V=0}$

(51)

In the derivation of Eq.(51) we have used the standard expressions [4,5]

${\displaystyle {}^{n}\rho ={}^{n+1}\rho J\quad {\hbox{and}}\quad d{~}^{n+1}V=Jd{~}^{n}V}$

(52)

Eqs.(39) and (51) together with the constitutive relationship (47) and the boundary conditions (12) complete the set of governing variational equations for a fluid in the updated Lagrangian description. The solution of these equations with the FEM is described in the next section.

9 FINITE ELEMENT DISCRETIZATION

We interpolate the velocities and the pressure in terms of their nodal values in the standard finite element fashion [28,39,41]. For a mesh of ${\textstyle n}$-noded ${\textstyle C^{\circ }}$ continuous elements we can write

${\displaystyle {}^{n+1}{v}={N}_{v}{}^{n+1}{\bar {v}}\quad ,\quad p={N}_{p}{}^{n+1}{\bar {p}}}$

(53)

where

${\displaystyle {\bar {v}}=\left\{{\begin{matrix}{\bar {v}}^{1}\\{\bar {v}}^{2}\\\vdots \\{\bar {v}}^{N}\end{matrix}}\right\}\quad {\hbox{with }}{\bar {v}}^{i}=[{\bar {v}}_{1}^{i},{\bar {v}}_{2}^{i},{\bar {v}}_{3}^{i}]^{T}\quad ,\quad {\bar {p}}=[{\bar {p}}_{1},{\bar {p}}_{2},\cdots ,{\bar {p}}_{N}]^{T}}$

(54.a)

where ${\textstyle N}$ is the total number of nodes in the mesh, ${\textstyle {\bar {v}}_{j}^{i}}$ and ${\textstyle {\bar {p}}^{i}}$ are ${\textstyle j}$th velocity component and the pressure unknowns for node ${\textstyle i}$,

${\displaystyle {\begin{array}{l}{N}_{v}=[{N}_{1}^{v},{N}_{2}^{v},\cdots ,{N}_{N}^{v}]\quad ,\quad {N}_{i}^{v}=N_{i}^{v}{I}_{3}\\{N}_{p}=[{N}_{1}^{p},{N}_{2}^{p},\cdots ,{N}_{N}^{p}]\end{array}}}$

(54.b)

In Eq.(54.b) ${\textstyle {N}_{i}^{v}}$ and ${\textstyle N_{i}^{p}}$ are the global shape functions [28,39] for the velocity and pressure interpolations for node ${\textstyle i}$, respectively.

The velocity and pressure increments and the virtual velocities are interpolated in terms of their nodal values in the same fashion as in Eq.(53).

The strain rate vector ${\textstyle \left\{{\dot {E}}\right\}}$ and its virtual expression ${\textstyle \left\{\delta {\dot {E}}\right\}}$ are respectively expressed in terms of the nodal velocities ${\textstyle {\bar {v}}}$ and their virtual values ${\textstyle \delta {\bar {v}}}$ via Eqs.(28.a), (29.a) and (53) as

${\displaystyle \left\{{\dot {E}}\right\}={B}{}^{n+1}{\bar {v}}\quad ,\quad \left\{\delta {\dot {E}}\right\}={B}\delta {\bar {v}}}$

(55)

The actual and virtual volumetric material strain rates are related to the virtual velocities as

${\displaystyle {\dot {E}}_{v}=\left\{{C}^{-T}\right\}{B}{~}^{n+1}{\bar {v}}\quad ,\quad \delta {\dot {E}}_{v}=\delta {\bar {v}}^{T}{B}^{T}\left\{{C}^{-1}\right\}}$

(56)

In the above equations ${\textstyle {B}}$ is the material strain rate matrix given by

${\displaystyle {B}=[{B}_{1},{B}_{2},\cdots ,{B}_{N}]^{T}\quad ,\quad {B}_{i}={}_{n}{B}_{i}^{0}+{B}_{i}^{L}}$

(57)

where

${\displaystyle {\begin{array}{l}{}_{n}{B}_{i}^{0}=\left[{\begin{matrix}{}_{n}N_{i,1}^{v}&0&0\\0&{}_{n}N_{i,2}^{v}&0\\0&0&{}_{n}N_{i,3}^{v}\\{}_{n}N_{i,2}^{v}&{}_{n}N_{i,1}^{v}&0\\{}_{n}N_{i,3}^{v}&0&{}_{n}N_{i,1}^{v}\\0&{}_{n}N_{i,3}^{v}&{}_{n}N_{i,2}^{v}\end{matrix}}\right]~~{\hbox{and}}~~{B}_{i}^{L}={}_{n}{B}_{i}^{0}{L}^{T}~~{\hbox{with}}~~{L}=\left[{\begin{matrix}l_{11}&l_{12}&l_{13}\\l_{21}&l_{22}&l_{23}\\l_{31}&l_{32}&l_{33}\end{matrix}}\right]\end{array}}}$

(58.a)

are the linear and non linear counterparts of the material strain rate matrix, respectively. The components of ${\textstyle {}_{n}{B}_{i}^{0}}$ and ${\textstyle {L}}$ are

${\displaystyle {}_{n}N_{i,j}^{v}={\partial N_{i}^{v} \over \partial {}^{n}x_{j}}\quad ,\quad l_{ij}=\sum \limits _{k=1}^{n}{}_{n}N_{k,j}^{v}\Delta u_{k}}$

(58.b)

The deviatoric second Piola-Kirchhoff stresses are related to the nodal velocities via Eqs.(47) and (55) as

${\displaystyle \left\{\mathbf {S} '\right\}=\left[{\boldsymbol {\cal {C}}}\right]\left\{{\dot {E}}\right\}=\left[{\boldsymbol {\cal {C}}}\right]{B}{~}^{n+1}{\bar {v}}}$

(59)

Remark 4. The displacement increment ${\textstyle \Delta u_{i}}$ can be computed in terms of the velocity in a number of ways. A popular choice is

${\displaystyle \Delta u_{i}=\Delta t{~}^{n+\alpha }v_{i}}$

(60.a)

where ${\textstyle {~}^{n+\alpha }v_{i}}$ is the velocity at ${\textstyle t=t_{n}+\alpha \Delta t}$ where ${\textstyle \alpha }$ is a positive number (${\textstyle 0\leq \alpha \leq 1}$). The value of ${\textstyle {~}^{n+\alpha }v_{i}}$ is typically computed as

${\displaystyle {~}^{n+\alpha }v_{i}=(1-\alpha )^{n}v_{i}+\alpha {}^{n+1}v_{i}}$

(60.b)

9.1 Discretized form of the PVP

Substituting Eqs.(53), (55) and (56) into the PVP (Eq.(39)) we obtain, after simplifying the virtual velocities

${\displaystyle {}^{n+1}{\bar {r}}_{m}:=\int _{{}^{n}V}{}^{n}\rho {N}_{v}^{T}{N}_{v}{\dot {\bar {v}}}d{~}^{n}V+\int _{{}^{n}V}{B}^{T}\left[\left\{\mathbf {S} '\right\}+J\left\{{C}^{-1}\right\}{p}\right]d{~}^{n}V-}$ ${\displaystyle \int _{{}^{n+1}V}{N}_{v}^{T}{~}^{n+1}{b}~d{~}^{n+1}V-\int _{{}^{n+1}\Gamma }{N}_{v}^{T}{~}^{n+1}{t}~d{~}^{n+1}\Gamma ={0}}$

(61)

where ${\textstyle {}^{n+1}{\bar {r}}_{m}}$ is the residual vector of the discretized momentum equations in the updated configuration and ${\textstyle {\dot {\bar {v}}}\equiv {\partial {}^{n+1}{\bar {v}} \over \partial t}}$ is the nodal acceleration vector.

Eq.(61) can be written in a more compact form as

${\displaystyle \displaystyle {}^{n+1}{\bar {r}}_{m}:={M}_{v}{\dot {\bar {v}}}+{}^{n+1}{g}_{v}+{}^{n+1}{g}_{p}-{}^{n+1}{f}_{m}}$

(62)

where

${\displaystyle {M}_{v}=\!\!\int _{{}^{n}V}{}^{n}\rho {N}_{v}^{T}{N}_{v}d~{}^{n}V~,~{}^{n+1}{g}_{v}=\!\!\int _{{}^{n}V}{B}^{T}{S}'d{~}^{n}V~,~{}^{n+1}{g}_{p}=\!\!\int _{{}^{n}V}{B}^{T}J\left\{{C}^{-1}\right\}pd{~}^{n}V}$ ${\displaystyle {}^{n+1}{f}_{m}=\!\!\int _{{}^{n+1}V}{N}_{v}^{T}{~}^{n+1}~{b}~d{~}^{n+1}V+\!\!\int _{{}^{n+1}\Gamma }{N}_{v}^{T}{~}^{n+1}{t}~d{~}^{n+1}V}$

(63)

In Eq.(62) ${\textstyle {}^{n+1}{g}_{v}}$ and ${\textstyle {}^{n+1}{g}_{p}}$ are internal force vectors contributed by the deviatoric second Piola-Kirchhoff stresses and the pressure, respectively and ${\textstyle {}^{n+1}{f}_{m}}$ is the external force vector.

Recall that the internal force vectors at the current configuration are obtained in terms of expressions at the updated configuration.

Remark 5. The computation of the body force vector ${\textstyle {}^{n+1}{b}}$ in Eq.(63) for the case of selfweight requires computing the density at the updated configuration as ${\textstyle {}^{n+1}\rho ={\frac {1}{J}}{}^{n}\rho }$. We also note that the surface tractions ${\textstyle {}^{n+1}{t}}$ are applied on the boundary of the updated configuration ${\textstyle {}^{n+1}\Gamma }$, which requires the accurate identification of this boundary.

The acceleration term in Eq.(62) can be approximated in a number of manners. A backward Euler scheme gives

${\displaystyle {}^{n+1}{\bar {r}}_{m}:={M}_{v}{\frac {{}^{n+1}{\bar {v}}-{}^{n}{\bar {v}}}{\Delta t}}+{}^{n+1}{g}_{v}+{}^{n+1}{g}_{p}-{}^{n+1}{f}_{m}=0}$

(64)

9.2 Discretization of the mass conservation equation

The arbitrary pressure test functions ${\textstyle q}$ are interpolated in the same fashion as the pressure as

${\displaystyle q={N}_{p}{\bar {q}}={\bar {q}}^{T}{N}_{p}^{T}}$

(65)

where vector ${\textstyle {\bar {q}}}$ contains the nodal values of ${\textstyle q}$.

Substituting the expression of ${\textstyle {\dot {E}}_{v}}$ in term of ${\textstyle {}^{n+1}{\bar {v}}}$ from Eq.(56) and Eq.(65) in the variational form of the mass balance equation (Eq.(51)) we obtain, after eliminating the pressure test functions ${\textstyle {\bar {q}}}$

${\displaystyle \displaystyle {}^{n+1}{\bar {r}}_{v}:=-{M}_{p}{\dot {\bar {p}}}+{Q}^{T}{~}^{n+1}{\bar {v}}={0}}$

(66)

where ${\textstyle {}^{n+1}{\bar {r}}_{v}}$ is the residual vector of the discretized mass conservation equation, ${\textstyle {\dot {\bar {p}}}={\partial {~}^{n+1}{\bar {p}} \over \partial t}}$ and the terms of ${\textstyle {M}_{p}}$ and ${\textstyle {}^{n}{Q}}$ are

${\displaystyle M_{p_{ij}}=\int _{{}^{n}V}{\frac {J^{2}}{^{n}\rho c^{2}}}N_{i}^{p}N_{j}^{p}d{~}^{n}V~~,~~{Q}=\int _{{}^{n}V}{B}^{T}\left\{{C}^{-1}\right\}{N}_{p}Jd{~}^{n}V}$

(67)

9.3 Stabilization of the mass conservation equation

It is well known that the FEM solutions for the fully incompressible limit are unstable for some particular forms of the approximation for the velocities and the pressure [10,39,41]. This is the case, for instance, when an equal order interpolation is chosen for the velocity components and the pressure (i.e. ${\textstyle N_{i}^{v}=N_{i}^{p}}$). This problem can be overcome by using finite element approximations for ${\textstyle {v}}$ and ${\textstyle p}$ satisfying the so-called LBB conditions [10,39,41], or else by introducing stabilization terms into the discretized mass balance equation [10,41]. In this work the later approach is chosen for obtaining stabilized numerical solutions.

In essence all stabilized formulations modify the discretized form of the mass balance equation as

${\displaystyle \displaystyle {}^{n+1}{\hat {r}}_{v}:={}^{n+1}{\bar {r}}_{v}+^{n+1}{\bar {r}}_{s}=0}$

(68)

where ${\textstyle {}^{n+1}{\bar {r}}_{v}}$ was defined in Eq.(66) and ${\textstyle {}^{n+1}{\bar {r}}_{s}}$ is a stabilization mass balance vector which expression can be typically written as

${\displaystyle {}^{n+1}{\bar {r}}_{s}=-{\boldsymbol {\cal {S}}}{~}^{n+1}{\bar {p}}+{}^{n+1}{f}_{s}}$

(69)

where ${\textstyle {\boldsymbol {\cal {S}}}}$ is a stabilization matrix and ${\textstyle {}^{n+1}{f}_{s}}$ is a force vector that depends on the nodal velocities and pressures. The form of ${\textstyle {\boldsymbol {\cal {S}}}}$ and ${\textstyle {f}_{s}}$ is different for each stabilization method. Typically, ${\textstyle {\boldsymbol {\cal {S}}}}$ and ${\textstyle {}^{n+1}{f}_{s}}$ contain integrals over the volume and the Neumann boundary of the analysis domain and are both a function of the so-called stabilization parameters which expression also depends on the particular stabilization method chosen [3,8,9,10,15,31,37,41].

The derivation of a stabilized mixed velocity-pressure formulation for incompressible fluids exceeds the objectives of this paper. Details can be obtained in the references cited in the previous paragraphs.

A particular stabilized Lagrangian formulation with excellent mass conservation features based on the Finite Calculus (FIC) theory [21]–[24],[27,30] can be found in [31].

10 INCREMENTAL SOLUTION FOR THE NODAL VELOCITIES AND PRESSURES

We will derive next an incremental iterative procedure for solving the discretized form of the equations for conservation of linear momentum and mass conservation in a segregated form. This requires the linearization of Eqs.(61) and (68) with respect to the nodal velocity and pressure unknowns, respectively. The linearization procedure takes advantage from the fact that the reference configuration (i.e. the current configuration ${\textstyle {}^{n}V}$) remains fixed during the linearization process.

10.1 Linearization of the discretized momentum equations with respect to the nodal velocities

We linearize the residual vector ${\textstyle {}^{n+1}{\bar {r}}_{m}}$ of Eq.(61) with respect to the nodal velocities ${\textstyle {\bar {v}}}$ as

${\displaystyle \delta _{\bar {v}}{}^{n+1}{\bar {r}}_{m}({x},{\bar {v}},{\bar {p}})={\frac {d}{d\epsilon }}{\Bigg |}_{\epsilon =0}{}^{n+1}{\bar {r}}_{m}({x},{\bar {v}}+\epsilon d{\bar {v}},{\bar {p}})}$

(70)

Expression (70) is the directional derivative of the residual vector ${\textstyle {}^{n+1}{\bar {r}}_{m}}$ at a point ${\textstyle {x}}$ in the direction of the nodal velocity vector ${\textstyle d{\bar {v}}}$ (hereafter called nodal velocity pseudo-increment vector). The same definition applies to the directional derivative of a matrix or a scalar depending on the space coordinate ${\textstyle {x}}$, the nodal velocities ${\textstyle {\bar {v}}}$ and the nodal pressures ${\textstyle {\bar {p}}}$ [4,5,36].

Using the expression of vector ${\textstyle {}^{n+1}{\bar {r}}_{m}}$ of Eq.(62) in Eq.(70) and neglecting the changes of vector ${\textstyle {}^{n+1}{f}_{m}}$ with the velocity, gives

${\displaystyle \delta _{\bar {v}}{}^{n+1}{\bar {r}}_{m}={\frac {1}{\Delta t}}{M}_{v}d{\bar {v}}+\int _{{}^{n}V}{\Big \{}{B}^{T}{\Big [}\delta _{\bar {v}}\left\{\mathbf {S} '\right\}+\left(\delta _{\bar {v}}J\left\{{C}^{-1}\right\}+{J}\delta _{\bar {v}}\left\{{C}^{-1}\right\}\right)p+}$ ${\displaystyle +J\left\{{C}^{-1}\right\}\delta _{\bar {v}}p{\Big ]}+\delta _{\bar {v}}{B}^{T}\left\{\mathbf {S} '\right\}{\Big \}}d{~}^{n}V}$

(71)

After some algebra detailed in the Appendix, the following linearized form of ${\textstyle {}^{n+1}{\bar {r}}_{m}}$ with respect to the nodal velocities is found

${\displaystyle \delta _{v}{}^{n+1}{\bar {r}}_{m}=\left({\frac {1}{\Delta t}}{M}_{v}+{K}_{c}+{K}_{\sigma }\right)d{\bar {v}}+{Q}{\delta }_{\bar {v}}^{n+1}{\bar {p}}}$

(72)

where matrices ${\textstyle {M}_{v}}$ and ${\textstyle {}^{n}{Q}}$ were given in Eq.(67) and ${\textstyle {}^{n}{K}_{c}}$ and ${\textstyle {K}_{\sigma }}$ are the constitutive and initial stress matrices, respectively. The expression of these matrices is

${\displaystyle {K}_{c}=\!\!\int _{{}^{n}V}{B}^{T}[{\boldsymbol {\cal {C}}}_{T}]{B}d{~}^{n}V~~,~~{K}_{\sigma }=\Delta t\!\int _{{}^{n}V}{G}^{T}{\hat {S}}'{G}d{~}^{n}V}$

(73)

The form of matrices ${\textstyle {G}}$ and ${\textstyle {\hat {S}}'}$ is given in the Appendix.

The expression of the tangent constitutive tensor ${\textstyle {\boldsymbol {\cal {C}}}_{T}}$ is (see Appendix)

${\displaystyle {\boldsymbol {\cal {C}}}_{T}={\boldsymbol {\cal {C}}}+J\Delta t~p({C}^{-1}\otimes {C}^{-1}-2{\boldsymbol {\cal {I}}})}$

(74)

Tensor ${\textstyle {\boldsymbol {\cal {C}}}_{T}}$ contains contributions from the constitutive tensor ${\textstyle {\boldsymbol {\cal {C}}}}$ of Eq.(46), the time step and the pressure. This form of ${\textstyle {\boldsymbol {\cal {C}}}_{T}}$ is very similar to that used for incompressible continua [4,5,13,36]. It can be shown that tensor ${\textstyle {\boldsymbol {\cal {I}}}}$ is symmetric (see Appendix) and, hence, the tangent constitutive tensor ${\textstyle {\boldsymbol {\cal {C}}}_{T}}$ is symmetric [4,5,13,36].

10.2 Linearization of the nodal pressures with respect to the nodal velocities. Derivation of the tangent matrix

From Eq.(66) we can obtain the directional derivative of the nodal pressure variables in the direction of the nodal velocity increments. Using a trapezoidal rule for approximating the time derivative term gives

${\displaystyle {M}_{p}{\frac {1}{\theta \Delta t}}({~}^{n+1}{\bar {p}}-{}^{n+\theta }{\bar {p}})+{Q}^{T}{\bar {v}}={0}}$

(75)

where ${\textstyle \theta }$ is a time parameter such that ${\textstyle 0<\theta \leq 1}$. A value ${\textstyle \theta =1}$ defines a backward Euler scheme [39,41].

From Eq.(75) it is straightforward to obtain

${\displaystyle {\delta }_{\bar {v}}{~}^{n+1}{\bar {p}}=\theta \Delta t{M}_{p}^{-1}{Q}^{T}d{\bar {v}}}$

(76)

Substituting Eq.(76) into (72) gives the final linearized form of the momentum equations in terms of the nodal velocities increments as

${\displaystyle {\delta }_{\bar {v}}{}^{n+1}{\bar {r}}_{m}={K}_{T}^{i}d{\bar {v}}}$

(77)

where ${\textstyle (\cdot )^{i}}$ denotes values at the ${\textstyle i}$th iteration and the tangent matrix ${\textstyle ^{n}{K}_{T}}$ is

${\displaystyle \displaystyle {K}_{T}={\frac {1}{\Delta t}}{M}_{v}+{K}_{c}+{K}_{\sigma }+{K}_{p}}$

(78)

with ${\textstyle {M}_{v}}$, ${\textstyle {K}_{c}}$ and ${\textstyle {K}_{\sigma }}$ defined in Eqs.(63) and (72) and the tangent “bulk” matrix ${\textstyle {K}_{p}}$ is

${\displaystyle {K}_{p}=\theta \Delta t{Q}{M}_{p}^{-1}{Q}^{T}}$

(79)

In practice, ${\textstyle {K}_{p}}$ is computed using the diagonal form of ${\textstyle {M}_{p}}$, i.e.

${\displaystyle {K}_{p}=\theta \Delta t{Q}{M}_{pD}^{-1}{\bar {Q}}^{T}}$

(80)

where ${\textstyle {M}_{pD}={diag}({M}_{p})}$.

The incremental form of the momentum equations can written as

${\displaystyle {}^{n+1}{\bar {r}}_{m}\simeq {}^{n+1}{\bar {r}}_{m}^{i}+{\delta }_{\bar {v}}{~}^{n+1}{\bar {r}}_{m}={}^{n+1}{\bar {r}}_{m}^{i}+{K}_{T}^{i}d{\bar {v}}=0}$

(81)

Solution of Eq.(81) yields the nodal velocity increments ${\textstyle d{\bar {v}}}$ for the ${\textstyle i}$th iteration.

Remark 6. For fully incompressible fluids (${\textstyle c=\infty }$) a large but finite value of ${\textstyle c}$ is used in practice. This allows to eliminate the pressure DOFs in the momentum equations via Eq.((76)).

Remark 7. The tangent “bulk” stiffness matrix ${\textstyle \mathbf {K} _{p}}$ in the tangent matrix ${\textstyle \mathbf {K} _{T}}$ accounts for the changes of the pressure due to the velocity. Including matrix ${\textstyle \mathbf {K} _{p}}$ in ${\textstyle \mathbf {K} _{T}}$ has proven to be essential for the fast convergence, mass preservation and overall accuracy of the iterative solution in all cases [11]. The ${\textstyle \theta }$ parameter in ${\textstyle {K}_{p}}$ (Eq.(79)) has the role of preventing the ill-conditioning of ${\textstyle \mathbf {K} _{T}}$ for very large values of the speed of sound in the fluid that lead to a dominant role of the terms of the tangent bulk stiffness matrix ${\textstyle {K}_{p}}$. Clearly, the value of the terms of ${\textstyle {K}_{p}}$ can also be limited by reducing the time step size. This, however, leads to an increase in the overall cost of the computations [11]. An adequate selection of ${\textstyle \theta }$ improves the convergence of the iterative process and leads to a more accurate numerical solution with reduced mass loss, even for large time steps [11]. These considerations, however, do not affect the value of ${\textstyle c}$ within matrix ${\textstyle {M}_{p}}$ in Eq.(67) that vanishes for the fully incompressible case (${\textstyle c=\infty }$).

10.3 Linearization of the stabilized mass conservation equation with respect to the nodal pressures

The stabilized mass conservation equation (68) is linearized with respect to the nodal pressure pseudo-increment vector ${\textstyle d{\bar {p}}}$ using Eqs.(66), (68) and (69) as

${\displaystyle \delta _{\bar {p}}{}^{n+1}{\hat {r}}_{v}({x},{\bar {v}},{\bar {p}})={\frac {d}{d\epsilon }}{\Bigg |}_{\epsilon =0}{}^{n+1}{\hat {r}}_{v}({x},{\bar {v}},{\bar {p}}+\epsilon d{\bar {p}})}$

(82)

Using Eqs.(67)–(68) and a backward Euler scheme for approximating the time derivative of the nodal pressure in Eq.(67) gives

${\displaystyle \delta _{\bar {p}}{~}^{n+1}{\hat {r}}_{v}=-\left({\frac {1}{\Delta t}}{M}_{p}+{\boldsymbol {\cal {S}}}\right)d{\bar {p}}}$

(83)

In the derivation of Eq.(83) we have neglected the pressure dependence of the terms of matrix ${\textstyle {\boldsymbol {\cal {S}}}}$ and ${\textstyle {}^{n+1}{f}_{s}}$ in Eq.(69).

The incremental form of the mass balance equation is therefore

${\displaystyle {}^{n+1}{\hat {r}}_{v}\simeq {}^{n+1}{\hat {r}}_{v}^{i}+\delta _{\bar {p}}{~}^{n+1}{\hat {r}}_{v}={}^{n+1}{\hat {r}}_{v}^{i}-\left({\frac {1}{\Delta t}}{M}_{p}+{\boldsymbol {\cal {S}}}^{i}\right)d{\bar {p}}=0}$

(84)

Solution of Eq.(84) yields the nodal pressure pseudo-increment vector ${\textstyle d{\bar {p}}}$ at the ${\textstyle i}$th iteration.

10.4 Incremental solution of the discretized equations

An incremental Newton-Raphson type iterative solution scheme for the stabilized velocity-pressure formulation is as follows.

Within a time increment ${\textstyle [n,n+1]}$:

• Initialize variables: ${\textstyle ({}^{n+1}{x},{}^{n+1}{\bar {v}}^{1},{}^{n+1}{\bar {p}}^{1},{}^{n+1}{\bar {r}}_{m},{}^{n+1}{\hat {r}}_{v})\equiv ({}^{n}{x},{}^{n}{\bar {v}},{}^{n}{\bar {p}},{}^{n}{\bar {r}}_{m},{}^{n}{\hat {r}}_{v})}$.
• Iteration loop = ${\textstyle i=1,\cdots ,NITER}$

For each iteration

1. Compute the nodal velocity pseudo-increments , ${\textstyle d{\bar {v}}}$ (from Eq.(81)) from

   ${\displaystyle {K}_{T}^{i}d{\bar {v}}=-{}^{n+1}{\bar {r}}_{m}^{i}({}^{n+1}{\bar {v}}^{i},{}^{n+1}{\bar {p}}^{i})</li>}$

   (85)

2. Update the nodal velocities :

   ${\displaystyle ^{n+1}{\bar {v}}^{i+1}=^{n+1}{\bar {v}}^{i}+d{\bar {v}}</li>}$

   (86)

3. Compute the nodal pressure pseudo-increments , ${\textstyle d{\bar {p}}}$. From Eq.(84),

   ${\displaystyle \left({\frac {1}{\Delta t}}{M}_{p}+{\boldsymbol {\cal {S}}}^{i}\right)d{\bar {p}}=^{n+1}{\hat {r}}_{v}^{i}({}^{n+1}{\bar {v}}^{i+1},{}^{n+1}{\bar {p}}^{i})</li>}$

   (87)

4. Update the nodal pressures :

   ${\displaystyle ^{n+1}{\bar {p}}^{i+1}=^{n+1}{\bar {p}}^{i}+d{\bar {p}}</li>}$

   (88)

5. Update the nodal displacement increments ${\textstyle \Delta {u}}$ using Eq.(60.a) and the approximate value for the nodal velocities ${\textstyle {}^{n+1}{\bar {v}}^{i+1}}$.
6. Update the nodal coordinates in the updated configuration as

   ${\displaystyle {}^{n+1}{x}^{i+1}={}^{n+1}{x}^{i}+\Delta {u}</li>}$

   (89)

7. Compute the material derivative of the Green strains ${\displaystyle {\dot {E}}_{ij}}$ and the deviatoric second Piola-Kirchhoff stresses ${\displaystyle {S'}_{ij}}$ (from Eqs.(55) and (45)).
8. Compute the momentum and mass balance residuals : ${\textstyle ^{n+1}{\bar {r}}_{m}^{i+1}}$ and ${\textstyle ^{n+1}{\hat {r}}_{v}^{i+1}}$
9. Check convergence

   ${\displaystyle {\begin{array}{l}\Vert ^{n+1}{\bar {r}}_{m}^{i+1}\Vert \leq e_{m}\Vert {}^{n}{f}_{v}\Vert \\[.5cm]\Vert ^{n+1}{\hat {r}}_{v}^{i+1}\Vert \leq e_{v}{~}^{n}V</li>\end{array}}}$

   (90)

   where ${\textstyle \Vert \cdot \Vert }$ denotes the quadratic norm and ${\textstyle e_{m}}$ and ${\textstyle e_{v}}$ are prescribed error tolerances for the discretized residuals of the momentum and mass balance equations, respectively.

   If both conditions (90) are satisfied then make ${\textstyle n\leftarrow n+1}$ and proceed to the next time increment. Otherwise, make the iteration counter ${\textstyle i\leftarrow i+1}$ and repeat Steps 1–9.

Remark 8. An alternative convergence criteria based on the nodal velocities and pressures can be defined as

${\displaystyle \Vert ^{n+1}{\bar {v}}^{i+1}-{}^{n+1}{\bar {v}}^{i}\Vert \leq e_{\bar {v}}\Vert ^{n+1}{\bar {v}}^{i}\Vert }$ (91.a) ${\displaystyle \Vert ^{n+1}{\bar {p}}^{i+1}-^{n+1}{\bar {p}}^{i}\Vert \leq e_{\bar {p}}\Vert ^{n+1}{\bar {p}}^{i}\Vert }$ (91.b)

where ${\textstyle e_{\bar {v}}}$ and ${\textstyle e_{\bar {p}}}$ are error tolerances.

Remark 9. The nodal velocity and pressure increment vectors ${\textstyle \Delta {\bar {v}}}$ and ${\textstyle \Delta {\bar {p}}}$ can be computed at the end of each time step as

${\displaystyle \Delta {\bar {v}}={}^{n+1}{\bar {v}}-{}^{n}{\bar {v}}\quad {\hbox{and}}\quad \Delta {\bar {p}}={}^{n+1}{\bar {p}}-{}^{n}{\bar {p}}}$

(92)

where ${\textstyle {}^{n+1}{\bar {v}}}$ and ${\textstyle {}^{n+1}{\bar {p}}}$ denote the converged values at the end of the iteration loop. Clearly, ${\textstyle \Delta {\bar {v}}}$ and ${\textstyle \Delta {\bar {p}}}$ can also be obtained by adding up the pseudo-increment vectors ${\textstyle d{\bar {v}}}$ and ${\textstyle d{\bar {p}}}$ within the iterative solution.

Remark 10. The nodal pressures ${\textstyle ^{n+1}{\bar {p}}^{i+1}}$ can be directly obtained from Eq.(68), after substitution of Eqs.(66) and (69), as

${\displaystyle \left({\frac {1}{\Delta t}}{M}_{p}+{\boldsymbol {\cal {S}}}^{i}\right){}^{n+1}{\bar {p}}^{i+1}={Q}^{T}{~}^{n+1}{\bar {v}}^{i+1}+{\frac {1}{\Delta t}}{M}_{p}{~}^{n}{\bar {p}}+{}^{n+1}{f}_{s}}$

(93)

The rest of the iterative solution scheme is similar to that described above. Eq.(93) substitutes Eqs.(87) and (88) and the convergence of the nodal pressures is verified by Eq.(91.b).

Remark 11. For a fully incompressible fluid ${\textstyle c=\infty }$ and matrix ${\textstyle {M}_{p}=0}$ in Eqs.(67), (87) and (93).

The problem can still be accurately solved in this case using an adequate expression for the stabilization matrix ${\textstyle {\boldsymbol {\cal {S}}}}$. The form of ${\textstyle {\boldsymbol {\cal {S}}}}$ presented in [31] includes a Laplacian term over the whole domain ${\textstyle \Omega }$ and an integral along the Neumann boundary ${\textstyle \Gamma _{t}}$. The boundary term in ${\textstyle {\boldsymbol {\cal {S}}}}$ avoids the need for prescribing the pressure on the domain boundary. If a standard Laplacian form is chosen for ${\textstyle {\boldsymbol {\cal {S}}}}$, then the value of the pressure has to be prescribed in strong form at some boundary points in order to obtain good results [41].

11 DIRECT ITERATIVE SOLUTION OF THE NODAL VELOCITIES AND PRESSURES

Substituting the FEM approximation for the velocities and the pressure (53) into Eq.(64) and assuming a Newtonian fluid with a constitutive equation given by Eq.(47), the following matrix expression of the discretized momentum equations can be obtained

${\displaystyle \displaystyle {}^{n+1}{\bar {r}}_{m}:={M}_{v}{\dot {\bar {v}}}+{K}_{v}{~}^{n+1}{\bar {v}}+{Q}{~}^{n+1}{\bar {p}}-{}^{n+1}{f}_{m}=0}$

(94)

where ${\textstyle {M}_{v}}$, ${\textstyle {}^{n}{Q}}$ and ${\textstyle {}^{n+1}{f}_{m}}$ have been defined previously and

${\displaystyle {K}_{v}=\int _{{}^{n}V}{B}^{T}\left[{\boldsymbol {\cal {C}}}\right]{B}d{~}^{n}V}$

(95)

Combining Eq.(94) with the stabilized mass balance equation (68) and using Eq.(66) gives the following matrix system for the nodal velocities and pressures

${\displaystyle \left[{\begin{matrix}{M}_{v}&{0}\\{0}&-{M}_{p}\end{matrix}}\right]\left\{{\begin{matrix}{\dot {\bar {v}}}\\{\dot {\bar {p}}}\end{matrix}}\right\}+\left[{\begin{matrix}{K}_{v}&{Q}\\{Q}^{T}&-{\boldsymbol {\cal {S}}}\end{matrix}}\right]\left\{{\begin{matrix}{}^{n+1}{\bar {v}}\\{}^{n+1}{\bar {p}}\end{matrix}}\right\}-\left\{{\begin{matrix}{}^{n+1}{f}_{m}\\{}^{n+1}{f}_{s}\end{matrix}}\right\}={0}}$

(96)

Eq.(96) can be written in a compact (monolithic) form as

${\displaystyle {M}{\dot {\bar {a}}}+{K}{~}^{n+1}{\bar {a}}-{}^{n+1}{f}=0}$

(97)

where

${\displaystyle {M}=\left[{\begin{matrix}{M}_{v}&{0}\\{0}&-{M}_{p}\end{matrix}}\right]~~,~~{K}=\left[{\begin{matrix}{K}_{v}&{Q}\\{Q}^{T}&-{\boldsymbol {\cal {S}}}\end{matrix}}\right]~~,~~{~}^{n+1}{\bar {a}}=\left\{{\begin{matrix}{}^{n+1}{\bar {v}}\\{}^{n+1}{\bar {p}}\end{matrix}}\right\}~~,~~{}^{n+1}{f}=\left\{{\begin{matrix}{}^{n+1}{f}_{m}\\{}^{n+1}{f}_{s}\end{matrix}}\right\}}$

(98)

11.1 Monolithic solution scheme

Eq.(98) is the basis for deriving monolithic iterative time integration schemes for directly computing the nodal velocities and pressure at the updated configuration. For instance, the standard backward Euler scheme gives

${\displaystyle {H}^{i}{~}^{n+1}{\bar {a}}^{i+1}={}^{n+1}{f}^{i}+{\frac {1}{\Delta t}}{M}{}^{n}{\bar {a}}}$

(99)

with ${\textstyle {H}^{i}={\frac {1}{\Delta t}}{M}+{K}_{v}^{i}}$.

The values of ${\textstyle {}^{n+1}{\bar {a}}}$ can be directly found by solving Eq.(99) iterative.

The non linearity of matrix ${\textstyle {}^{n}{K}_{v}}$ emanates from the non linear terms in the material strain rate matrix ${\textstyle {B}^{L}}$ involving the displacement increments (Eqs.(9)).

11.2 Segregated solution scheme

The nodal velocities and pressures at the updated configuration can be also computed starting from Eq.(97) using a segregated iterative scheme as follows

${\displaystyle \left[{\frac {1}{\Delta t}}{M}_{v}+{K}_{v}^{i}\right]{~}^{n+1}{\bar {v}}^{i+1}={}^{n+1}{f}_{m}-{Q}{~}^{n+1}{p}^{i}}$

(100)

${\displaystyle \left[{\frac {1}{\Delta t}}{M}_{p}+{\boldsymbol {\cal {S}}}^{i}\right]{~}^{n+1}{\bar {p}}^{i+1}={~}^{n+1}{f}_{s}+{Q}^{T}{~}^{n+1}{\bar {v}}^{i+1}+{\frac {1}{\Delta t}}{M}_{p}{~}^{n}{\bar {p}}}$

(101)

Eqs.(100) and (101) are solved sequentially and iteratively until convergence for the nodal velocities and pressures is found.

The above iterative scheme can be considered as a simplification of the more accurate incremental iterative segregated scheme described in Section 10.4.

Remark 12. A variant of the iteration matrix in the left hand side of Eq.(100) can be used by adding the bulk stiffness matrix ${\textstyle {K}_{p}}$ to matrix ${\textstyle {H}^{i}}$ [11,31].

12 PARTICULARIZATION FOR THE UPDATED CONFIGURATION

The finite element formulation presented in the previous sections can be particularized for the case when the updated configuration ${\textstyle {~}^{n+1}{V}}$ is chosen as the reference configuration.

The particularization is simple if we note that now ${\textstyle \Delta u_{i}=0}$. Hence, ${\textstyle {L}=0}$ and ${\textstyle {B}^{L}={0}}$ (see Eq.(58.a)). Also all the space derivatives are taken with respect to the updated coordinates ${\textstyle {~}^{n+1}x_{i}}$ and the integrals are carried out in the updated configuration ${\textstyle {~}^{n+1}V}$. In addition, the tangent constitutive tensor ${\textstyle {\boldsymbol {\cal {C}}}_{T}={c}}$. The expressions of the relevant matrices and vectors for this case are given in Box 1.

13 PROS AND CONS OF USING THE CURRENT OR THE UPDATED CONFIGURATION AS THE REFERENCE CONFIGURATION

We have shown in the previous sections that either the current or the updated configuration can be indistinctly used as the reference configuration for the finite element analysis of quasi and fully incompressible fluids using a Lagrangian formulation. Both choices imply solving a non linear system of equations. However, the nature of the non-linearity is different for each case.

If the current configuration is chosen as the reference configuration, all integrals in the tangent matrix are performed on the known configuration ${\textstyle {}^{n}V}$ which remains fixed during the iterative solution process. The non-linearity affects, however, the non linear terms of the material strain rate matrix (i.e. ${\textstyle {B}^{L}}$) and also the constitutive tensor ${\textstyle {\boldsymbol {\cal {C}}}}$ that involves the deformation gradient. A simplification of the tangent stiffness matrix for this case can be introduced by neglecting ${\textstyle {B}^{L}}$ in the material strain rate matrix and assuming that ${\textstyle {F}={C}={I}_{3}}$ and, hence, ${\textstyle {\boldsymbol {\cal {C}}}_{T}={c}}$.

If, on the other hand, the updated configuration is chosen as the reference configuration, then all the integrals in the tangent matrix are computed in the unknown configuration ${\textstyle {}^{n+1}V}$, which should be updated at each iteration. On the other hand, the expression for the material strain rate matrix is linear and also ${\textstyle {\boldsymbol {\cal {C}}}_{T}={c}}$. A simplification of the iterative process can be introduced by keeping the tangent matrix constant for a fixed number of iterations.

Indeed the above mentioned simplifications can affect the convergence rate of the iterative solution and should be implemented with care.

Which reference configuration should be chosen can be problem dependent and, certainly, the choice will affect the organization of the computer program and its efficiency. What should be kept in mind is that the final solution, i.e. the geometry of the updated configuration and the velocities and stresses on ${\textstyle {}^{n+1}V}$, should be identical in both cases.

In previous works of the authors with the Particle Finite Element Method (PFEM) the current configuration ${\textstyle {}^{n}V}$ was typically chosen as the reference configuration with the simplified form for the tangent matrix as explained above and also neglecting the initial stress matrix terms in ${\textstyle {K}_{T}}$ [16,17,25,26,29,31]. Recent experiences indicate that using the updated configuration ${\textstyle {}^{n+1}V}$ as the reference configuration can be advantageous in many Lagrangian flow problems. The topic is still open for research and hopefully this paper will be useful for choosing the adequate FEM expressions for each case.

14 CONCLUDING REMARKS

We have presented a mixed velocity-pressure finite element formulation for solving the updated Lagrangian equations for quasi and fully incompressible fluids. The finite element interpolation uses an equal order approximation for the velocity and pressure unknowns. We have detailed a number of iterative algorithms for solving the non-linear stabilized FEM equations for the velocities and the pressure at the nodes using incremental and direct solution schemes. The algorithms presented are useful for the study of Lagrangian flows, as well as for solving fluid-structure-interaction problems using a unified Lagrangian finite element formulation for modelling both the fluid and the structure [11,17].

The choice of the current or the updated configurations as the reference Lagrangian configuration is still an open topic. Researchers interested in Lagrangian CFD procedures will find in this paper the equations to be used for each case, from which simplifications or further computational refinements can be made.

ACKNOWLEDGEMENTS

This work was partially supported by the Advanced Grant project SAFECON of the European Research Council (ERC).

References

[1] S. Badia, R. Codina, Unified stabilized finite element formulation for the Stokes and the Darcy problems. SIAM J. on Numerical Analysis, 47, 1971–2000, 2009.

[2] K.J. Bathe, Finite element procedures. Prentice-Hall, 1996.

[3] Y. Bazilevs, K Takizawa, T.E Tezduyar, Computational Fluid-Structure Interaction: Methods and Applications, Wiley, 2013.

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

[5] J. Bonet and R.D. Wood, Non linear continuum mechanics for finite element analysis. Wiley, 2nd Edition, 2008.

[6] J.U. Brackbill, D.B. Kothe, H.M. Ruppel, FLIP: A low-dissipation, particle-in-cell method for fluid flow. Computer Physics Communications, 48, 25–38, 1988.

[7] D. Burgess, D. Sulsky, J.U. Brackbill, Mass matrix formulation of the FLIP particle-in-cell method. Journal of Computational Physics, 103 (1), 1–15, 1992.

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

[9] R. Codina, H. Coppola-Owen, P. Nithiarasu and C. Liu, Numerical comparison of CBS and SGS as stabilization techniques for the incompressible Navier-Stokes equations. Int. J. Numer. Meth. Engng., 66, 1672–1689, 2006.

[10] J. Donea and A. Huerta, Finite Element Methods for Flow Problems. Wiley, 2003.

[11] A. Franci, E. Oñate, J.M. Carbonell, Unified Lagrangian formulation for analysis of fluid-structure interaction problems. Research Report PI-400, CIMNE, Barcelona 2013.

[12] F.H. Harlow, The particle-in-cell computing method for fluid dynamics. Methods Comput. Phys., 3, 219, 1963.

[13] G.A. Holzaphel, Non linear solid mechanics. Wiley, 2000.

[14] A. Huerta, Y. Vidal, J. Bonet, Updated Lagrangian formulation for corrected Smooth Particle Hydrodynamics. Int. J. of Computational Methods, 3 (4), 383–399, 2006.

[15] T.J.R. Hughes, G. Scovazzi, L.P. Franca, Multiscale and stabilized methods. Encyclopedia of Comput. Mechanics. E. Stein, R. de Borst and T.J.R. Hughes (Eds.), Vol. 3 Fluids, 5–60, Wiley, 2004.

[16] 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. Num. Meth. Eng., 61(7), 964–989, 2004.

[17] S.R. Idelsohn, J. Marti, A. Limache, E. Oñate Unified Lagrangian formulation for elastic solids and incompressible fluids: Application to fluid-structure interaction problems via the PFEM. Comput. Meth. Appl. Mech. Engrg., 197 (19-20), 1762–1776, 2008.

[18] S. Li, W.K. Liu, Meshfree and particle methods and their applications. Appl. Mech. Rev., 55, 1, 2002.

[19] W.K. Liu, Y. Chen, S. Jun, J.S. Chen, T. Belytschko, C. Pan, R.A. Uras, C.T. Chang, Overview and applications of the Reproducing Kernel Particle Methods. Arch Comput Methods Eng., 3 (1), 3–80, 1996.

[20] M.B. Liu and G.R. Liu, Smoothed particle hydrodynamics (SPH): an overview and recent developments. Arch Comput Methods Eng., Vol. 17, 25–-76, 2010.

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

[22] 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), 355–370, 2000.

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

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

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

[26] E. Oñate, J. García, S.R. Idelsohn, F. Del Pin, FIC formulations for finite element analysis of incompressible flows. Eulerian, ALE and Lagrangian approaches. Comput. Meth. Appl. Mech. Engrg. 195(23-24), 3001–3037, 2006.

[27] 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, 609–637, 2007.

[28] E. Oñate, Structural analysis with the finite element method. Linear statics. Volume 1. Basis and Solids. CIMNE-Springer, 2009.

[29] 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. Computational Mechanics, 48(3), 307–318, 2011.

[30] E. Oñate, S.R. Idelsohn, C. Felippa, Consistent pressure Laplacian stabilization for incompressible continua via higher-order finite calculus. Int. J. Numer. Meth. Engng, 87 (1-5), 2011, 171–195.

[31] E. Oñate, A. Franci, J.M. Carbonell, Lagrangian formulation for finite element analysis of incompressible fluids with reduced mass losses. Int. J. Numer. Meth. Fluids, 2013. DOI:0.1002/fld.3870.

[32] E. Oñate, P. Nadukandi, S.R. Idelsohn, P1/P0+ elements for incompressible flows with discontinuous material properties. Comput. Meth. Appl. Mech. Engrg., 271, 185–209, 2014, .

[33] N.A. Patankar, D.D. Joseph, Lagrangian numerical simulation of particulate flows. Int. J. of Multiphase Flow, 27 (10), 1685–-1706, 2001.

[34] R. Radovitzky and M. Ortiz, Lagrangian finite element analysis of newtonian fluid flows. Int. J. Numer. Meth. Engng., 43, 607–619, 1998.

[35] B. Ramaswamy and M. Kawahara, Lagrangian finite element analysis applied to viscous free surface fluid flow. Int. J. Num. Meth. Fluids, 7, 953–984, 1986.

[36] P. Wriggers, Non linear finite element methods. Springer, 2008.

[37] T.E. Tezduyar, Finite elements for flow problems with moving boundaries and interfaces. Arch. for Comput. Methods Eng., 8, 83–130, 2001.

[38] D.Z. Zhang, Q. Zou, W.B. VanderHeyden, X. Ma, Material point method applied to multiphase flows. Journal of Computational Physics, 227 (6), 3159-3173, 2008.

[39] O.C. Zienkiewicz, R.L. Taylor, J.Z. Zhu, The Finite Element Method. Vol. 1 The Basis. Elsevier, 6th Edition, 2005.

[40] O.C. Zienkiewicz, R.L. Taylor, The Finite Element Method. Vol. 2 Solid and Structural Mechanics. Elsevier, 6th Edition, 2005.

[41] O.C. Zienkiewicz, R.L. Taylor and P. Nithiarasu, The Finite Element Method. Vol. 3 Fluid Dynamics. Elsevier, 6th Edition, 2005.

APPENDIX A. LINEARIZATION OF THE MOMENTUM EQUATIONS WITH RESPECT TO THE NODAL VELOCITIES

Using the expression of ${\textstyle {}^{n+1}{\bar {r}}_{m}}$ of Eq.(62) and neglecting the changes of the external vector ${\textstyle {}^{n+1}{\bar {f}}_{m}}$ with the velocity (accounting for these changes is possible and will lead to additional terms in the tangent matrix), we can write

${\displaystyle \delta _{\bar {v}}{}^{n+1}{\bar {r}}_{m}={\frac {1}{\Delta t}}{M}_{v}d{\bar {v}}+\int _{{}^{n}V}{\Big \{}{B}^{T}{\Big [}\delta _{\bar {v}}\left\{\mathbf {S} '\right\}+\left(\delta _{\bar {v}}J\left\{{C}^{-1}\right\}+{J}\delta _{\bar {v}}\left\{{C}^{-1}\right\}\right)p+}$ ${\displaystyle +J\left\{{C}^{-1}\right\}\delta _{\bar {v}}p{\Big ]}+\delta _{\bar {v}}{B}^{T}\left\{\mathbf {S} '\right\}{\Big \}}d{~}^{n}V}$

(A.1)

Introducing Eqs.(55) and (59) into (A.1) gives after some algebra [4,5,36]

${\displaystyle \delta _{\bar {v}}\left\{\mathbf {S} '\right\}=[{\boldsymbol {\cal {C}}}]\delta _{\bar {v}}\left\{{\dot {E}}\right\}=[{\boldsymbol {\cal {C}}}]{B}d{\bar {v}}}$

(A.2)

${\displaystyle \delta _{\bar {v}}J\left\{{C}^{-1}\right\}=J{C}^{-1}\otimes {C}^{-1}\delta _{\bar {v}}\left\{{E}\right\}=J\Delta t{C}^{-1}\otimes {C}^{-1}{B}d{\bar {v}}}$

(A.3)

${\displaystyle J\delta _{\bar {v}}\left\{{C}^{-1}\right\}=-2J{\boldsymbol {\cal {I}}}\delta _{\bar {v}}\left\{{E}\right\}=-2J\Delta t{\boldsymbol {\cal {I}}}{B}d{\bar {v}}}$

(A.4)

where

${\displaystyle {\boldsymbol {\cal {I}}}_{ijkl}={\frac {1}{2}}\left[({C}^{-1})_{ik}({C}^{-1})_{jl}-({C}^{-1})_{il}({C}^{-1})_{jk}\right]}$

(A.5)

It can be shown that tensor ${\textstyle {\boldsymbol {\cal {I}}}}$ is symmetric [4,5].

On the other hand,

${\displaystyle \delta _{\bar {v}}{B}^{T}\left\{{S}'\right\}=\Delta t{G}^{T}{\hat {S}}'{G}d{\bar {v}}}$

(A.6a)

with

${\displaystyle {\begin{array}{l}{G}=\left[{\begin{matrix}{\bar {G}}&{\bar {0}}&{\bar {0}}\\{\bar {0}}&{\bar {G}}&{\bar {0}}\\{\bar {0}}&{\bar {0}}&{\bar {G}}\end{matrix}}\right]~~,~~{\bar {G}}=\left[{\begin{matrix}{}_{n}N_{1,1}^{v}&0&0&{}_{n}N_{2,1}^{v}&0&0&\cdots &{}_{n}N_{n,1}^{v}\\{}_{n}N_{1,2}^{v}&0&0&{}_{n}N_{2,2}^{v}&0&0&\cdots &{}_{n}N_{n,2}^{v}\\{}_{n}N_{1,3}^{v}&0&0&{}_{n}N_{2,3}^{v}&0&0&\cdots &{}_{n}N_{n,3}^{v}\end{matrix}}\right]\\[1cm]{\hat {S}}'=\left[{\begin{matrix}{S}'&{0}&{0}\\{0}&{S}'&{0}\\{0}&{0}&{S}'\end{matrix}}\right]~~,~~{0}=\left[{\begin{matrix}0&0&0\\0&0&0\\0&0&0\end{matrix}}\right]~~,~~{\bar {0}}=\left\{{\begin{matrix}0\\0\\0\end{matrix}}\right\}\end{array}}}$

(A.6b)

with ${\textstyle {}_{n}N_{i,j}^{v}}$ defined in Eq.(58.b).

In the derivation of Eqs.(A.3), (A.4) and (A.6a) we have assumed that ${\textstyle d(\Delta u_{i})=du_{i}=\Delta t~d({}^{n+\alpha }v_{i})=\Delta t~dv_{i}}$. This relationship follows from Eq.(60.a).

Substituting Eqs.(A.2)–(A.4) and (A.6a) and the interpolation for the pressure (Eq.(53)) into (A.1) yields the linearized form of the residual vector of the discretized momentum equations as

${\displaystyle \delta _{v}{}^{n+1}{\bar {r}}_{m}={\frac {1}{\Delta t}}{M}_{v}d{\bar {v}}+\left\{\int _{{}^{n}V}\left[{B}^{T}[{\boldsymbol {\cal {C}}}_{T}]{B}+{G}^{T}{\hat {S}}{G}\right]d{~}^{n}V\right\}d{\bar {v}}+}$ ${\displaystyle +\left\{\int _{{}^{n}V}{B}^{T}\left\{{C}^{-1}\right\}{N}_{p}Jd{~}^{n}V\right\}{\delta }_{\bar {v}}^{n+1}{\bar {p}}=\left({\frac {1}{\Delta t}}{M}_{v}+{K}_{c}+{K}_{\sigma }\right)d{\bar {v}}+{Q}{\delta }_{\bar {v}}^{n+1}{\bar {p}}}$

(A.7)

where matrices ${\textstyle {M}_{v}}$ and ${\textstyle {Q}}$ were given in Eq.(67) and ${\textstyle {K}_{c}}$ and ${\textstyle {K}_{\sigma }}$ are the constitutive and initial stress matrices, respectively. The expression of these matrices is

${\displaystyle {K}_{c}=\!\!\int _{{}^{n}V}{B}^{T}[{\boldsymbol {\cal {C}}}_{T}]{B}d{~}^{n}V~~,~~{K}_{\sigma }=\!\!\int _{{}^{n}V}{G}^{T}{\hat {S}}'{G}d{~}^{n}V}$

(A.8)

The tangent constitutive tensor ${\textstyle {\boldsymbol {\cal {C}}}_{T}}$ is deduced from Eqs.(A.1)-(A.4) as

${\displaystyle {\boldsymbol {\cal {C}}}_{T}={\boldsymbol {\cal {C}}}+J\Delta tp({C}^{-1}\otimes {C}^{-1}-2{\boldsymbol {\cal {I}}})}$

(A.9)

It is straightforward to show that tensor ${\textstyle {\boldsymbol {\cal {C}}}_{T}}$ is symmetric.