Unified updated Lagrangian formulation for the analysis of quasi and fully incompressible fluids and solids and their interaction via a partitioned scheme and the PFEM

==Abstract==

In this work, a unified updated Lagrangian formulation for solving fluid-structure interaction (FSI) problems is derived. The mixed velocity-pressure formulation for hypoelastic solids and quasi and fully incompressible Newtonian fluids is obtained as an extension of the velocity formulation derived for a general continuum. The space discretization for the fluid domain is performed via the Particle Finite Element Method (PFEM), where for the solid domain a standard FEM is used. Linear interpolation is used for both the velocity and the pressure fields. The global FSI problem is solved using a Gauss-Seidel iterative scheme.  The required stabilization for dealing with incompressible situations is given by an enhanced formulation of the Finite Calculus (FIC) technique <span id='citeF-22'></span>[[#cite-22|[22]]]. The efficiency of the proposed strategy is tested by solving benchmark FSI problems.

'''keywords''' Unified Updated Lagrangian Formulation, Quasi and Fully Incompressible, Partitioned Scheme, Particle Finite Element Method

==1 Introduction==

The great interest of the computational mechanics community in solving fluid-structure interaction (FSI) problems is due to two main reasons. On the one hand the analytical solution of these problems is generally impossible to obtain and experimental tests can not be performed, or they have an excessive cost. On the other hand, FSI problems are very attractive for their multidisciplinarity: they occur in many fields, such as civil, aerospace and nuclear engineering, or biology.

In the last decades, many numerical methods have been proposed for solving FSI problems. The most common classification distinguishes the ''monolithic'' and the ''staggered'' (or ''partitioned'') approaches. The first group includes those strategies in which the fluid and the structure are solved within the same global linear system (<span id='citeF-13'></span>[[#cite-13|[13]]], <span id='citeF-20'></span>[[#cite-20|[20]]] and <span id='citeF-26'></span>[[#cite-26|[26]]], among others). Instead, staggered methods solve the FSI problem via an iterative loop where the fluid solver and solid solvers are activated alternatively. The fluid and the solid domains interact tranfering the Dirichlet and Neumann conditions through the interface (<span id='citeF-11'></span>[[#cite-11|[11]]], <span id='citeF-28'></span>[[#cite-28|[28]]] and <span id='citeF-9'></span>[[#cite-9|[9]]]).

In this paper, a unified updated Lagrangian finite element formulation for solving FSI problems is derived and discussed. The present method belongs to the monolithic category of FSI strategies. The aim of the work is to give a general procedure for dealing with solid and fluid mechanics problems indifferently and for modeling their interaction in a natural way. With this purpose, the derivation of the formulation is developed so that the changes required by a specific material are minimized and introduced only at the end of the derivation. The materials considered in this work are hypoelastic solids and Newtonian fluids. However, it will be shown that the formulation can be easily extended to other types of materials, such as incompressible or elasto-plastic solids among others. The mixed velocity pressure formulation emanates from the velocity formulation and both cases are derived first for a general continuum. Only later the formulations are adapted to specific materials. The finite element approximation is performed by interpolating both the velocity and the pressure fields using linear shape functions.

Concerning the mixed formulation, the global system of algebraic equations is solved via an iterative partitioned Gauss-Seidel scheme. This means that first the momentum equations are solved in terms of the unknown velocity increments using the pressures computed at the previous iteration and then the unknown pressures are computed using the updated velocities. This procedure is the key point of the unified formulation. In fact, it allows to treat hypoelastic solids and quasi-incompressible fluids in a similar way. In particular, in the momentum equations, the structure of the tangent matrix and the terms of the right hand side are the same for both materials.

The main differences in the treatment of fluids and solids are in the incompressibility of the fluid and in the distortion of its mesh. In order to overcome the latter difficulty, the Particle Finite Element Method has been adopted for the analysis of the fluid domain. The PFEM is a Lagrangian formulation based on an Alpha-Shape <span id='citeF-10'></span>[[#cite-10|[10]]]  remeshing procedure that allows one  to deal with a fine and regular mesh all over the duration of the analysis. Many scientific publications have shown the efficiency of the PFEM in the simulation of free-surface fluids (<span id='citeF-16'></span>[[#cite-16|[16]]] and  <span id='citeF-17'></span>[[#cite-17|[17]]] among others) and FSI problems (<span id='citeF-18'></span>[[#cite-18|[18]]] and  <span id='citeF-23'></span>[[#cite-23|[23]]] among others).

In regard to the fluid incompressibility (or quasi-incompressibility) a specific stabilization is required because the linear interpolations chosen for the velocity and the pressure fields do not fulfil the so-called <math display="inline">{LBB  inf-sup  condition}</math>  <span id='citeF-4'></span>[[#cite-4|[4]]]. In this work, the mass balance equation in the fluid is stabilized using the updated formulation of Finite Increment Calculus (FIC) technique <span id='citeF-22'></span>[[#cite-22|[22]]]. This choice is motivated by the small intrusively of this method (its stabilization terms affect the continuity equation only) and by its mass preservation features <span id='citeF-22'></span>[[#cite-22|[22]]].

The FSI solver is based on the mixed velocity-pressure formulations for hypoelastic solids and Newtonian fluids. Fluids and solids are computed by using the same unknown variables (velocities and pressure), the same framework (Lagrangian) and the same solving scheme (Gauss-Seidel). All this simplifies the coupling for solving FSI problems allowing the use of a monolithic scheme; <math display="inline">{}</math> fluids and solids can be easily solved using a similar system of equations ensuring a strong coupling automatically. The FSI coupling is performed in a natural way and it essentially consists in realizing two tasks: to detect the interface creating the coupling elements (elements that share at least one interface node) and to assembly properly the global system.

The outline of paper is the following. First the velocity formulation in the updated Lagrangian framework is derived for a general continuum. Then the unknown pressures are introduced and the mixed velocity-pressure formulation is presented. In the third section, the formulations are adapted to specific materials. First, both the velocity and the mixed formulations are specified for the case of hypoelastic solids. Then the mixed formulation is adapted for the analysis of quasi-incompressible free surface fluids. Some reminders remarks about the stabilization procedure and the Particle Finite Element Method are also given. In the fifth section, the coupling algorithm for solving FSI problems is explained. Particular attention is devoted to the detection of the interface and the assembly of the global system. In the next section, some representative problems are solved and discussed in order to test the efficiency of the proposed FSI solution strategy. In particular, the collapse of a water column on a deformable membrane and the water entry of a nylon sphere are analyzed. Finally, the main conclusions are given.

==2 Velocity formulation==

In this section, the velocity formulation for solving transient problems for a general continuum is derived. The governing equations are the linear momentum equations and they are derived in the updated Lagrangian (UL) framework. The essential feature of a Lagrangian description is that the independent variables are the Lagrangian coordinates <span id='citeF-3'></span>[[#cite-3|[3]]]. For this reason, a typical Eulerian measure, as the Cauchy stress tensor, can be used in a Lagrangian framework if it is expressed in function of the Lagrangian coordinates. In the UL formulation used in this work the governing equations are integrated over the unknown configuration <math display="inline">{\Omega }</math> (the so-called updated configuration). As a consequence, the space derivatives for the UL description are computed with respect to the spatial coordinates.

===2.1 From local form to the spatial semi-discretization===

In this section the spatial semi-discretization of the linear momentum equations is derived.

For a general continuum, the local form of the linear momentum equations using the updated Lagrangian description reads:

<span id="eq-1"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>\rho (\hbox{X},t) \frac{ \partial \hbox{v}(\hbox{X},t)}{\partial t}-{\partial \sigma (\hbox{X},t) \over \partial \hbox{x}}-\hbox{b}(\hbox{X},t)=0   \quad \quad in   \Omega \times (0,T)   </math>
|}
| style="width: 5px;text-align: right;" | (1)
|}

where <math display="inline">\rho </math> is the density of the material, <math display="inline">v</math> are the velocities, <math display="inline">\sigma </math> is the Cauchy stress tensor and <math display="inline">b</math> is the body force. The variables within the brackets are the independent variables. In particular,  '''X''' are the Lagrangian  or material coordinates, '''x''' the Eulerian or spatial coordinates and <math display="inline">t</math> is the time. For simplicity, in the following the independent variables will be not specified.

The set of governing equations is completed by the following conditions at the Dirichlet (<math display="inline">\Gamma _v</math>) and Neumann (<math display="inline">\Gamma _t</math>) boundaries:

<span id="eq-2"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>v_i -v_i^p =0 \qquad \hbox{on }\Gamma _v </math>
|}
| style="width: 5px;text-align: right;" | (2)
|}

<span id="eq-3"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>\sigma _{ij}n_j -t_i^p =0 \qquad \hbox{on }\Gamma _t </math>
|}
| style="width: 5px;text-align: right;" | (3)
|}

where <math display="inline">v_i^p</math> and <math display="inline">t_i^p</math>, <math display="inline">i=1,...,n_s</math> are the prescribed velocities and the prescribed tractions, respectively.

The spaces for the trial and test functions are defined, respectively, as:

<span id="eq-4"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>v_ i \in {U},  \quad \quad {U}=\{ v_ i |v_ i  \in C^0,   v_ i=v_i^p    on  \Gamma _v \}   </math>
|}
| style="width: 5px;text-align: right;" | (4)
|}

<span id="eq-5"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>w_ i \in {U_0}, \quad \quad {U_0}=\{ w_ i | w_ i  \in C^0,   w_ i=0    on  \Gamma _v \}   </math>
|}
| style="width: 5px;text-align: right;" | (5)
|}

Multiplying  Eqs.([[#eq-1|1]]) by the test functions and integrating over the updated configuration domain, the following global integral form is obtained:

<span id="eq-6"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>\int _\Omega w_i \left(\rho \frac{\partial v_i}{\partial t}- {\partial \sigma _{ij} \over \partial x_j}-b_i\right)d\Omega =0 </math>
|}
| style="width: 5px;text-align: right;" | (6)
|}

Integrating by parts  the term involving <math display="inline">\sigma _{ij}</math> in Eq.([[#eq-6|6]]) and using the Neumann boundary conditions ([[#eq-3|3]]) yields the weak variational form of the momentum equations as:

<span id="eq-7"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>\int _\Omega w_i \rho \frac{\partial v_i}{\partial t} d\Omega + \int _\Omega \frac{\partial w_i}{\partial x_j} \sigma _{ij} d\Omega - \int _\Omega w_i b_i d\Omega - \int _{\Gamma _t} w_i t_i^p d\Gamma =0  </math>
|}
| style="width: 5px;text-align: right;" | (7)
|}

Eq.([[#eq-7|7]]) is the standard form of the principle of virtual power <span id='citeF-3'></span>[[#cite-3|[3]]].

The spatial discretization is introduced using the classical FEM-Galerkin prodedure. Hence both the trial and the test functions are interpolated in the space by means of the same shape functions <math display="inline">{N}</math>.

<span id="eq-8"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>v_i = \sum \limits _{I=1}^n N_I \bar v_{iI} \quad ,\quad w_i = \sum \limits _{I=1}^n N_I \bar w_{iI}     </math>
|}
| style="width: 5px;text-align: right;" | (8)
|}

where <math display="inline">n=3/4</math> for 2D/3D problems is the number of the nodes of the finite element, <math display="inline">\bar{(\cdot )}</math> denotes a nodal value, the capital subscript specifies the node and the lower case subscript represents the cartesian direction. In this work, the velocities have been interpolated using linear shape functions.

Using the arbitrariness of the test functions and introducing the spatial discretization ([[#eq-8|8]]) into  Eq.([[#eq-7|7]]), the spatial semi-discretized form of the momentum equations in the UL framework reads:

<span id="eq-9"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>\underbrace{  \int _\Omega N_I \rho d\Omega  \dot{v_i}  }_{{{W}^{k}_{Ii}}}  +  \underbrace{  \int _\Omega \frac{\partial N_I}{\partial x_j} \sigma _{ij} d\Omega  }_{{{W}^{int}_{Ii}}} = \underbrace{  \int _\Omega N_I b_i d\Omega + \int _{\Gamma _t} N_I t_i^p d\Gamma  }_{{{W}^{ext}_{Ii}}}  </math>
|}
| style="width: 5px;text-align: right;" | (9)
|}

For convenience, the semi-discretized form of the momentum equations in the total Lagrangian (TL) framework is also presented here. This is written as <span id='citeF-3'></span>[[#cite-3|[3]]]:

<span id="eq-10"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>\underbrace{  \int _{\Omega _0} N_I \rho _0 d\Omega  \dot{v_i} }_{{^{TL}{W}^{k}_{Ii}}}+ \underbrace{  \int _{\Omega _0} \frac{\partial N_I}{\partial X_j} P_{ij} d\Omega }_{{^{TL}{W}^{int}_{Ii}}} = \underbrace{  \int _{\Omega _0} N_I b_i d\Omega + \int _{\Gamma ^0_t} N_I t_i^{0p} d\Gamma  }_{{^{TL}{W}^{ext}_{Ii}}} </math>
|}
| style="width: 5px;text-align: right;" | (10)
|}

where <math display="inline">{P}</math> is the first Piola-Kirchhoff stress tensor, or the nominal stress tensor, and all the variables with the subscript <math display="inline">{(\cdot )}_0</math> refer to the last known configuration. Note that Eq.([[#eq-10|10]]) can be obtained from Eq.([[#eq-9|9]]) by properly performing pull back transformations on all its terms <span id='citeF-3'></span>[[#cite-3|[3]]].

For the sake of clarity in the notation, the terms referred to the TL description are denoted with the exponent <math display="inline">{^{TL}(\cdot )}</math>. If not specified, the variables belong to the UL description.

===2.2 Time integration===

In this work, the kinematic variables  have been integrated in time using a second order scheme. In particular, the implicit Newmark's integration rule has been adopted. For the general case, the Newmark's integration rule states that accelerations and displacements are computed as:

<span id="eq-11"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>{\dot{v}}_{n+1}</math>
| style="text-align: right;" | <math>= \frac{1}{\gamma \Delta t} \left(v_{n+1} -v_{n}\right)- \frac{1-\gamma }{\gamma }{\dot{v}}_{n}  </math>
|-
| style="text-align: center;" | <math> u_{n+1}</math>
| style="text-align: right;" | <math>=u_{n} + {\Delta t}\frac{\gamma - \beta }{\gamma }v_{n} + \Delta t\frac{\beta }{\gamma } v_{n+1}+ {\Delta t^2}\frac{\gamma{-} 2\beta }{2\gamma }{\dot{v}}_{n} </math>
|}
| style="width: 5px;text-align: right;" | (11)
|}

Where <math display="inline">{\beta }</math> and <math display="inline">{\gamma }</math> are the so-called Newmark's parameters <span id='citeF-3'></span>[[#cite-3|[3]]]. A time integration scheme is unconditionally stable if the following relation holds:

<span id="eq-12"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>\gamma \ge 2 \beta \ge \frac{1}{2} </math>
|}
| style="width: 5px;text-align: right;" | (12)
|}

In the present work, the time integration scheme is implicit and the Newmark's parameters chosen are <math display="inline">{\beta = \frac{1}{4}}</math> and <math display="inline">{\gamma = \frac{1}{2}}</math> in order to fulfill relation ([[#eq-12|12]]).

Replacing the numerical values of the constants in Eq.([[#eq-11|11]]) yields:

<span id="eq-13"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>{\dot{v}}_{n+1}= \frac{2}{ \Delta t} \left(v_{n+1} -v_{n}\right)- {\dot{v}}_{n}  </math>
|}
| style="width: 5px;text-align: right;" | (13)
|}

<span id="eq-14"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>u_{n+1}=u_{n} + \frac{\Delta t}{2}\left(v_{n+1}+v_{n}\right)  </math>
|}
| style="width: 5px;text-align: right;" | (14)
|}

===2.3 Linearization===

Although the problem is setted out in the UL framework, the linearization of the momentum equations is performed first on the TL  semi-discretized form ([[#eq-10|10]]). The updated Lagrangian linearized form is subsequently obtained by performing a push-forward transformation on the total Lagrangian form. The reason for this is that the derivation of the tangent matrix in the total Lagrangian framework is easier. In fact, in Eq.([[#eq-10|10]]) the only variable that depends on time is the nominal stress <math display="inline">{P}</math>, while in the updated Lagrangian form ([[#eq-9|9]]) the time-dependent variables are the updated domain <math display="inline">{\Omega }</math>, the Cauchy stress tensor <math display="inline">{\sigma }</math> and the spatial derivatives <math display="inline">{   \partial N/ \partial x}</math>. For the sake of clarity, the linearization of the internal and kinematic work terms will be performed separately. Thus, in the following section  First the internal components of the tangent matrix are derived.

====2.3.1 Internal components of the tangent matrix====

From  Eq.([[#eq-10|10]]) the internal work in the TL description is defined as:

<span id="eq-15"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>^{TL}W^{int}_{Ii}=\int _{\Omega _0} \frac{\partial N_I}{\partial X_j} P_{ij} d\Omega  </math>
|}
| style="width: 5px;text-align: right;" | (15)
|}

In order to obtain the tangent matrix, the material time derivative of ([[#eq-15|15]]) is computed as:

<span id="eq-16"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>^{TL}\dot{W}^{int}_{Ii}=\int _{\Omega _0} \frac{\partial N_I}{\partial X_j} \dot{P}_{ij} d\Omega  </math>
|}
| style="width: 5px;text-align: right;" | (16)
|}

The material time derivative of the material part of the internal work is now discretized to obtain  its increment in a time interval <math display="inline">{\Delta t}</math>:

<span id="eq-17"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>\frac{\Delta ^{TL}W^{int}_{Ii}}{\Delta t}=\int _{\Omega _0} \frac{\partial N_I}{\partial X_j} \dot{P}_{ij} d\Omega  </math>
|}
| style="width: 5px;text-align: right;" | (17)
|}

From Eq.([[#eq-17|17]]) we deduce:

<span id="eq-18"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>\Delta ^{TL}W^{int}_{Ii}=\int _{\Omega _0} \frac{\partial N_I}{\partial X_j} \Delta t \dot{P}_{ij} d\Omega  </math>
|}
| style="width: 5px;text-align: right;" | (18)
|}

The first Piola-Kirchhoff stress tensor <math display="inline">{P}</math> is not typically used because it is not symmetric and its rate is a non-objective measure. For these reasons, the second Piola-Kirchhoff stress tensor <math display="inline">{S}</math> and its rate are  preferred quantities in the TL framework. The mentioned stress rate measures are related each other through the following relation:

<span id="eq-19"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>\dot{P}_{ij}=\dot{S}_{ir} F^T_{rj} + S_{ir} \dot{F}^T_{rj}  </math>
|}
| style="width: 5px;text-align: right;" | (19)
|}

where <math display="inline">{F}</math> is the so-called deformation gradient defined as:

<span id="eq-20"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>{F}_{ij}= \frac{\partial x_i}{\partial X_j}  </math>
|}
| style="width: 5px;text-align: right;" | (20)
|}

Substituting Eq.([[#eq-19|19]]) into ([[#eq-18|18]]) yields:

<span id="eq-21"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>\Delta ^{TL}{W}^{int}_{Ii}= \underbrace{\int _{\Omega _0} \frac{\partial N_I}{\partial X_j} F_{ir} \Delta t \dot{S}_{jr}d\Omega } _{{^{TL}\dot{W}^{m}_{Ii}}} +  \underbrace{\int _{\Omega _0} \frac{\partial N_I}{\partial X_j}  S_{ir} \Delta t \dot{F}^T_{rj} d\Omega } _{{^{TL}\dot{W}^{g}_{Ii}}}     </math>
|}
| style="width: 5px;text-align: right;" | (21)
|}

Eq.([[#eq-17|17]]) shows that the internal power can be split into the material and the geometric parts, <math display="inline">{\dot{W}^{m}}</math> and <math display="inline">{\dot{W}^{g}}</math>, respectively. The former accounts for the material response through the rate of the second Piola-Kirchhoff stress tensor, the second term is the initial stress term that contains the information of the updated stress field. It is to be noted that, up to now, no constitutive relations have been introduced and the present derivation holds for a general continuum. 

<math>{{Material  tangent   matrix}}</math>

For the derivation of the material tangent matrix, the constitutive relation between the stress and the strain measures needs to be defined. In order to maintain the formulation as general as possible, the stress rate measure is related to the deformation rate through a tangent modulus tensor, such that

<span id="eq-22"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>\dot{S}_{ij}=C_{ijkl} \dot{E}_{kl}=C_{ijkl}  \frac{\partial N_J}{\partial X_s}  F_{kl}  \bar  v_{sJ}  </math>
|}
| style="width: 5px;text-align: right;" | (22)
|}

where <math display="inline">{C}</math> is a fourth-order tangent moduli tensor and <math display="inline">{E}</math> is the Green-Lagrange strain tensor and the following substitution has been made:

<span id="eq-23"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>\dot{E}_{kl}= \frac{\partial N_J}{\partial X_s}  F_{kl}  \bar  v_{sJ}  </math>
|}
| style="width: 5px;text-align: right;" | (23)
|}

One may note that in literature (<span id='citeF-3'></span>[[#cite-3|[3]]],<span id='citeF-2'></span>[[#cite-2|[2]]]), often the term <math display="inline">{\frac{\partial N_I}{\partial X}  F}</math> is grouped in the matrix <math display="inline">{B^I_{0}}</math> defined in two dimension as: 


{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| style="text-align: center;" |
<math>\begin{array}{l} \\       {B^I_{0}}= \left[                      \begin{array}{cccccc}                        \frac{\partial N_I}{\partial {X}}\frac{\partial x}{\partial {X}} & \frac{\partial N_I}{\partial {X}}\frac{\partial y}{\partial {X}} \\  \frac{\partial N_I}{\partial {Y}}\frac{\partial x}{\partial {Y}} & \frac{\partial N_I}{\partial {Y}}\frac{\partial y}{\partial {Y}} \\  \frac{\partial N_I}{\partial {X}}\frac{\partial x}{\partial {Y}}+\frac{\partial N_I}{\partial {Y}}\frac{\partial x}{\partial {X}} &  \frac{\partial N_I}{\partial {X}}\frac{\partial y}{\partial {Y}}+\frac{\partial N_I}{\partial {Y}}\frac{\partial y}{\partial {X}} \\                      \end{array} \right] \\ \\ \end{array}</math>
|}


In this work, for convenience, the matrix <math display="inline">{B^I_{0}}</math> is not used. As it will be shown in the following sections, Eq.([[#eq-22|22]]) can represent  both a Kirchhoff solid material and a quasi-incompressible fluid. If different constitutive laws are used, Eq.([[#eq-22|22]]) should be modified accordingly in order to derive the material part of the tangent matrix.

Substituting Eq.([[#eq-22|22]]) into <math display="inline">{\dot{W}^{m}}</math> of  Eq.([[#eq-17|17]]) yields

<span id="eq-24"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>\Delta ^{TL}{W}^{m}_{Ii}=\int _{\Omega _0} \frac{\partial N_I}{\partial X_j}  F_{rj} \Delta t C_{ijkl}  \frac{\partial N_J}{\partial X_s}  F_{kl}  d\Omega   \bar   v_{sJ}  </math>
|}
| style="width: 5px;text-align: right;" | (24)
|}

where <math display="inline">{v_{sJ}}</math> is the <math display="inline">{s}</math>-component of the velocity of node <math display="inline">{J}</math>.

From Eq.([[#eq-24|24]]), the material tangent matrix in the TL description can be computed as

<span id="eq-25"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>^{TL}{K}^{m}_{IJis}=\int _{\Omega _0} \frac{\partial N_I}{\partial X_j}  F_{rj} \Delta t C_{ijkl}  \frac{\partial N_J}{\partial X_s}  F_{kl}  d\Omega  </math>
|}
| style="width: 5px;text-align: right;" | (25)
|}

The material tangent matrix for the UL framework can be derived by applying a push-forward transformation on each term of Eq.([[#eq-25|25]]). The following relations hold

<span id="eq-26"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>{\Omega _0}=\frac{\Omega }{J}   and   d{\Omega _0}=\frac{ d \Omega }{J}   </math>
|}
| style="width: 5px;text-align: right;" | (26)
|}

<span id="eq-27"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>\frac{\partial N_I}{\partial X_j}=\frac{\partial N_I}{\partial x_k} F_{kj}  </math>
|}
| style="width: 5px;text-align: right;" | (27)
|}

<span id="eq-28"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>C_{ijkl}=F^{-1}_{mi} F^{-1}_{nj} F^{-1}_{ok} F^{-1}_{pl} c^{\tau }_{mnop} = F^{-1}_{mi} F^{-1}_{nj} F^{-1}_{ok} F^{-1}_{pl} c^{\sigma }_{mnop} J  </math>
|}
| style="width: 5px;text-align: right;" | (28)
|}

where <math display="inline">{c^{\tau }}</math> is the tangent moduli corresponding to the material time derivative of the Kirchhoff stress tensor and <math display="inline">{c^{ \sigma }}</math> is the tangent moduli for the rate of the Cauchy stress. The rate of the Cauchy stress tensor is related to the rate of deformation <math display="inline">{D}</math> through the fourth-order tensor <math display="inline">{C}</math> in the following way

<span id="eq-29"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>\sigma ^{\bigtriangledown }=c^{\sigma } : D  </math>
|}
| style="width: 5px;text-align: right;" | (29)
|}

Substituting Eqs.([[#eq-26|26]]-[[#eq-28|28]]) into ([[#eq-25|25]]), the material tangent matrix for the UL description is obtained as

<span id="eq-30"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>K^{m}_{IJis}=\int _{\Omega } \frac{\partial N_I}{\partial x_j} \Delta t  c^{\sigma }_{ijkl}  \frac{\partial N_J}{\partial x_s}    d\Omega  </math>
|}
| style="width: 5px;text-align: right;" | (30)
|}

<math>{Geometric  tangent   matrix}</math>

The geometric tangent matrix for the UL framework is here derived by using the same procedure adopted for the material components. Hence, first the linearization is performed on the TL form and then the UL tangent matrix is obtained by performing the required transformation over the TL terms.

From Eq.([[#eq-17|17]])

<span id="eq-31"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>\Delta ^{TL}{W}^{g}_{Ii}=\int _{\Omega _0} \frac{\partial N_I}{\partial X_j} \Delta t  S_{ir}  \dot{F}^T_{rj} d\Omega   </math>
|}
| style="width: 5px;text-align: right;" | (31)
|}

where the rate of the deformation gradient is defined as

<span id="eq-32"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>\dot{F}_{ij}= \frac{\partial N_J}{\partial X_i}  \bar v_{Jj}  </math>
|}
| style="width: 5px;text-align: right;" | (32)
|}

Substituting Eq.([[#eq-32|32]]) into ([[#eq-31|31]]), the geometric components of the internal power in the TL description can be written as

<span id="eq-33"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>\Delta ^{TL}{W}^{g}_{Ii}=\int _{\Omega _0} \frac{\partial N_I}{\partial X_j}  \Delta t S_{ir} \frac{\partial N_J}{\partial X_r}   d\Omega   \bar  v_{Jj}  </math>
|}
| style="width: 5px;text-align: right;" | (33)
|}

The geometric tangent matrix reads:

<span id="eq-34"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>^{TL}K^{g}_{IJij}=\int _{\Omega _0} \frac{\partial N_I}{\partial X_j} \Delta t   S_{ir} \frac{\partial N_J}{\partial X_r} d\Omega   </math>
|}
| style="width: 5px;text-align: right;" | (34)
|}

In order to recover the UL form, the Piola identity has to be recalled, <math display="inline">{}</math>

<span id="eq-35"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>{S}= {F}^{-1} {\sigma } {F}^{-T} J  </math>
|}
| style="width: 5px;text-align: right;" | (35)
|}

The geometric tangent matrix in the UL framework is obtained by substituting Eqs.([[#eq-26|26]]), ([[#eq-27|27]]) and ([[#eq-35|35]]) into ([[#eq-34|34]]). This leads to

<span id="eq-36"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>K^{g}_{IJij}=\int _{\Omega } \frac{\partial N_I}{\partial x_j} \Delta t  {\sigma }_{ir} \frac{\partial N_J}{\partial x_r} d\Omega   </math>
|}
| style="width: 5px;text-align: right;" | (36)
|}

====2.3.2 Kinematic components of the tangent matrix====

The kinematic components of the tangent matrix in the UL description can be derived directly from the UL kinematic term <math display="inline">{{W}^{k}_{Ii}}</math> of Eq.([[#eq-9|9]]).

<span id="eq-37"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>W^{k}_{Ii}= \int _\Omega N_I \rho d\Omega  \dot{v}_i   </math>
|}
| style="width: 5px;text-align: right;" | (37)
|}

Eq.([[#eq-37|37]]) has to be discretized on time with the purpose of replacing the accelerations with the velocities using the time integration described in Eq.([[#eq-13|13]]).

Introducing Eq.([[#eq-13|13]]) into ([[#eq-37|37]]) and differentiating with respect of the unknown velocities <math display="inline">{v_{n+1}}</math>, the kinematic components of the tangent matrix are obtained as:

<span id="eq-38"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>K^{k}_{IJ}= \int _\Omega N_I \frac{2\rho }{\Delta t}  N_J d\Omega   </math>
|}
| style="width: 5px;text-align: right;" | (38)
|}

The tangent matrix is computed as the sum of the internal and the kinematic components ([[#eq-30|30]]), ([[#eq-36|36]]) and  ([[#eq-38|38]]) as

<span id="eq-39"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>K= K^{m}+ K^{g} +K^{k}   </math>
|}
| style="width: 5px;text-align: right;" | (39)
|}

===2.4 Residual form and incremental solution scheme===

With the aim of deriving the incremental solution scheme, Eq.([[#eq-9|9]]) has to be rewritten in a residual form. Using Eq.([[#eq-13|13]]) and shifting all the terms to the left hand side, yields:

<span id="eq-40"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>R^{n+1}_{Ii}= \int _\Omega N_I \rho  N_J d\Omega   \bar{\dot{v}}^{n+1}_{Ji} +  \int _\Omega \frac{\partial N_I}{\partial x_j} \sigma ^{n+1}_{ij} d\Omega  -  \int _\Omega N_I b^{n+1}_i d\Omega  </math>
|-
| style="text-align: center;" | <math> - \int _{\Gamma _t} N_I t_i^{p,n+1} d\Gamma =0  </math>
|}
| style="width: 5px;text-align: right;" | (40)
|}

where <math display="inline">{R^{n+1}_{Ii}}</math> is the residual of the momentum equations referred to node <math display="inline">{I}</math> and the cartesian direction <math display="inline">{i}</math>. One can note that in Eq.([[#eq-40|40]]), as for Eq.([[#eq-36|36]]), the Cauchy stress tensor still appears. This tensor will be expressed as a function of the nodal unknowns of the problem in the following sections dedicated to the constitutive laws.

In Box 1 the iterative solution incremental scheme of the velocity formulation for a generic time increment <math display="inline">{}</math> of duration <math display="inline">{\Delta t}</math> is described.


<div class="center" style="font-size: 75%;">[[File:Draft_Samper_438761226_3662_Box1.png]]

'''Box 1'''. Iterative incremental solution scheme for the velocity formulation.
 </div>

==3 Mixed velocity-pressure formulation==

In solid and fluid mechanics, there are problems where a mixed formulation is recommendable, or even necessary. In fluid dynamics, the mixed formulation is required to apply properly the incompressibility constraint and to guarantee the stability of the problem. In solid mechanics, the mixed formulation represents the most reasonable choice to circumvent the locking problem in rubber-type materials, or in cases where plasticity or fracture is induced (<span id='citeF-4'></span>[[#cite-4|[4]]],<span id='citeF-6'></span>[[#cite-6|[6]]],<span id='citeF-7'></span>[[#cite-7|[7]]],<span id='citeF-24'></span>[[#cite-24|[24]]]). Also, for modelling FSI problems with standard compressible solids involved, the mixed formulation may represent an useful choice <span id='citeF-14'></span>[[#cite-14|[14]]]. In fact using the same variables for the analysis of the fluid and the solid represents both an important facility and a computational advantage: fluid and solid can be solved through an unified code and the risk of ill-conditioning the global problem is significally reduced. For all these reasons, the velocity formulation is here extended to the mixed velocity-pressure problems, considering the pressure as an additional unknown. Furthermore, in order to obtain a well-posed problem, the continuity equation is introduced in the solution scheme.

The whole problem is solved using a Gauss-Seidel partitioned iterative scheme. Hence, first the momentum equations are solved in terms of increments of the velocities and including the (known) pressures of the previous iteration in the residual, then the continuity equation is solved for the pressure using the updated velocities computed from the momentum equations. It will be shown that using this not intrusive scheme, it is possible to take advantage of most of the velocity-formulation derived so far.  In particular, the incremental velocity scheme for the momentum equations (Box 1) and the structure of the tangent matrix  ([[#eq-39|39]]) are still applicable. In this work, the same linear interpolation has been used for the velocity and the pressure fields. It is well known that, for incompressible (or quasi-incompressible) problems, this combination does not fulfill the <math display="inline">{inf-sup}</math> condition <span id='citeF-4'></span>[[#cite-4|[4]]] and a stabilization method is required.

===3.1 Splitting of the stress measures===

The Cauchy stress tensor is rewritten as the sum of its deviatoric and hydrostatic parts as

<span id="eq-41"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>\sigma =\sigma ' - pI </math>
|}
| style="width: 5px;text-align: right;" | (41)
|}

where <math display="inline">{\sigma '}</math> is the deviatoric part of the Cauchy stress, <math display="inline">{p}</math> is the pressure  and <math display="inline">{I}</math> is the identity tensor.

In terms of stress rate, that would be:

<span id="eq-42"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>\sigma ^{\bigtriangledown }=  \sigma ' ^{\bigtriangledown } - p^{\bigtriangledown }I    </math>
|}
| style="width: 5px;text-align: right;" | (42)
|}

where

<span id="eq-43"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>\sigma ' ^{\bigtriangledown } =  c^{\sigma '} : D^{dev}  </math>
|}
| style="width: 5px;text-align: right;" | (43)
|}

where <math display="inline">{{{\sigma } ' ^{\bigtriangledown }}}</math>  is the deviatoric part of the Cauchy stress rate tensor, <math display="inline">{p^{\bigtriangledown }}</math> is the rate of  the pressure, <math display="inline">{c^{\sigma '} }</math> is the tangent moduli for <math display="inline">{{{\sigma } ' ^{\bigtriangledown }}}</math> and <math display="inline">{D^{dev} }</math> is the deviatoric part of the deformation rate tensor. The pressure is related to the volumetric part of the rate of deformation by the following relation:

<span id="eq-44"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>\frac{1}{\kappa }\frac{\partial p}{\partial t}=  {D}^{v}  </math>
|}
| style="width: 5px;text-align: right;" | (44)
|}

where <math display="inline">{\kappa }</math> is the bulk modulus of the material and <math display="inline">{{D}^{v}= trace(D)}</math> is the volumetric strain rate.   Eq.([[#eq-44|44]]) is the closure equation for the mixed velocity-pressure formulation.

===3.2 Time integration===

As regards the variation on time of the pressure, a linear scheme has been adopted. That is  because in fluid dynamics the lagrangian mesh suffers large displacements and the nodal pressure information of the previous time steps has to be recovered by means of an adequate interpolation. This operation introduces approximation errors in the scheme that increase beyond the current time step. This also has a huge computational cost because the information of the previous steps has to be stored. With the procedure here presented only the current time step is needed and only the first variation in time of pressure has to be saved. Hence, the first and the second variations in time of the pressure are respectively:

<span id="eq-45"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>{^{n+1}\dot{p}}= \frac{ ^{n+1}{p} - {^{n}p} }{\Delta t}  </math>
|}
| style="width: 5px;text-align: right;" | (45)
|}

<span id="eq-46"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>{^{n+1}\ddot{p}}= \frac{^{n+1} {p} - {^{n}p} }{{\Delta t}^2}-\frac{^{n}\dot{p} }{\Delta t}   </math>
|}
| style="width: 5px;text-align: right;" | (46)
|}

Using  Eq.([[#eq-45|45]]), the time discretization of Eq.([[#eq-44|44]]) within the time interval <math display="inline">{[n,n+1]}</math> of duration <math display="inline">{\Delta t}</math> leads to:

<span id="eq-47"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>^{n+1}p=^{n}p + \Delta t \kappa ^{n+1}{D}^{v}  </math>
|}
| style="width: 5px;text-align: right;" | (47)
|}

===3.3 Fully discretized form and incremental solution scheme===

As already mentioned, the pressure field is interpolated with the same linear shape functions used for the velocity field (Eqs.([[#eq-8|8]])). So:

<span id="eq-48"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>p = \sum \limits _{I=1}^n N_I \bar p_{I}     </math>
|}
| style="width: 5px;text-align: right;" | (48)
|}

The Galerking-FEM approximation of the global form of Eq.([[#eq-47|47]]) leads to the following matrix form:

<span id="eq-49"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: right;" | <math>-Q^T {^{n+1}{\bar{v} }} + \frac{1}{\Delta t}M_1 {^{n+1}{\bar{p} }} </math>
| <math>= {^{n}g}  </math>
|}
| style="width: 5px;text-align: right;" | (49)
|}

where:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>{M}_{1_{IJ}} =\int _{\Omega ^e} N_I \frac{1}{\kappa } N_J d\Omega    </math>
|}
| style="width: 5px;text-align: right;" | (50)
|}

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>\hbox{Q}_{IJ}  =\int _{\Omega ^e} \hbox{B}_I^T \hbox{m} N_J d\Omega </math>
|-
| style="text-align: center;" | 
|}
| style="width: 5px;text-align: right;" | (51)
|}

<span id="eq-52"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>^{n}g_I=\int _{\Omega ^n}  N_I \frac{1}{\kappa \Delta t}  N_J d{\Omega }  {^{n}\bar p}_J   </math>
|}
| style="width: 5px;text-align: right;" | (52)
|}

with
<div class="center">
<math>\displaystyle{B}^e_I = \left[\begin{matrix} \displaystyle {\partial N_I \over \partial x} &0&0\\ \displaystyle{0}& \displaystyle {\partial N_I \over \partial y}&0\\ \displaystyle{0}&0&\displaystyle {\partial N_I \over \partial z}\\ \displaystyle {\partial N_I \over \partial y}&\displaystyle {\partial N_I \over \partial x}&0\\[.25cm] \displaystyle {\partial N_I \over \partial z}&0&\displaystyle {\partial N_I \over \partial x}\\[.25cm] \displaystyle{0}&\displaystyle {\partial N_I \over \partial z}&\displaystyle {\partial N_I \over \partial y}          \end{matrix}  \right]  \quad \quad \quad and \quad  \quad  \quad {m}=[1,1,1,0,0,0]^T </math>
</div>

Finally, for ensuring the coupling of the linear momentum equations with Eq.([[#eq-49|49]]), the pressure must be expressed explicitly  in the residual of the momentum equations (Eqs.([[#eq-40|40]])). For this reason, in  Eq.([[#eq-40|40]]) the Cauchy stress tensor is split into its deviatoric and pressure parts via Eq.([[#eq-41|41]]). The resulting residual form is

<span id="eq-53"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>R^{n+1}_{Ii}= \int _\Omega N_I \rho  N_J d\Omega   \bar{\dot{v}}^{n+1}_{Ji}+  \int _\Omega \frac{\partial N_I}{\partial x_j} \sigma ^{n+1}_{ij} d\Omega - \int _\Omega \frac{\partial N_I}{\partial x_j}\delta _{ij} N_J d\Omega \bar  p^{n+1}_J+ </math>
|-
| style="text-align: center;" | <math> -  \int _\Omega N_I b^{n+1}_i d\Omega - \int _{\Gamma _t} N_I t_i^{p,n+1} d\Gamma =0  </math>
|}
| style="width: 5px;text-align: right;" | (53)
|}

In Box 2, the iterative solution incremental solution scheme for a generic continuum using the mixed velocity-pressure formulation  is shown for a time increment <math display="inline">{}</math>.

<div class="center" style="font-size: 75%;">
[[File:Draft_Samper_438761226_7536_Box2.png]]

'''Box 2.'''  Iterative solution scheme for a generic continuum using the mixed velocity-pressure for-
mulation.</div>

==4 Constitutive laws==

In the formulation derived so far the constitutive model has not been specified. The only requirement stated so far is that the relation between the rate of the Cauchy stress and the rate of deformation is linear (Eq.([[#eq-29|29]])). This means that the constitutive relation must be rate-independent, incrementally linear and reversible <span id='citeF-25'></span>[[#cite-25|[25]]].For small increments the relation between the stress and strain increments is linear and they are recovered upon unloading. However, for large deformations, the energy is not necessarily conserved and the work done in a closed deformation path is not necessarily zero <span id='citeF-3'></span>[[#cite-3|[3]]].

As it has already been pointed out,  the Cauchy stress tensor, or its deviatoric part, still appear explicitly ([[#eq-36|36]], [[#eq-40|40]]) in the formulation. In the present section, the formulation will be adapted to specific constitutive laws and the stress tensor will be defined in terms of the nodal velocities.

First, it will be shown that both the velocity and the mixed velocity-pressure formulations hold for hypoelastic materials. Then, the case of quasi-incompressible Newtonian fluids will be studied. Due to the (quasi) incompressibility constraint, the fluid problem can be solved only using the mixed velocity-pressure formulation. Furthermore, because linear shape functions have been used for the interpolation of both the velocity and pressure fields, the <math display="inline">{inf-sup}</math> condition <span id='citeF-4'></span>[[#cite-4|[4]]] is not satisfied and the solution needs to be stabilized.

For both constitutive relations, the complete form of the tangent matrix will be derived and the incremental solution scheme will be explained in detail.

===4.1 Hypoelastic material===

The direct relation between the rate of stress and the rate of deformation is the main feature of hypoelastic material laws. In a large class of these, this relation is linear. The fundamental requirement for the stress rate measures is that  they must be objective, or, equally, frame-invariant. If not, many inconveniences can appear; for instance, rigid rotations can originate undesiderable states of stress. Many objective measures for the stress rates are available; the most common ones are the Truesdell and Jaumann measures of the Cauchy stress rate.

In this section, both the velocity and the mixed formulations are adapted for hypoelastic solids. 

<math>{Velocity   formulation}</math>

The Truesdell and Jaumann measures of the Cauchy stress rate are defined, respectively, as

<span id="eq-54"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>\sigma ^{\bigtriangledown \tau }=C^{\sigma \tau } : D  </math>
|}
| style="width: 5px;text-align: right;" | (54)
|}

<span id="eq-55"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>\sigma ^{\bigtriangledown J}=C^{\sigma J} : D  </math>
|}
| style="width: 5px;text-align: right;" | (55)
|}

where <math display="inline">{C^{\sigma \tau }}</math> and <math display="inline">{C^{\sigma J}}</math> are the Truesdell and Jaumann tangent moduli and the relation between the two tensors is

<span id="eq-56"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>C^{\sigma \tau }=C^{\sigma J} - C^*  </math>
|}
| style="width: 5px;text-align: right;" | (56)
|}

where:

<span id="eq-57"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>C^{\sigma J}_{ijkl}= \lambda \delta _{ij} \delta _{kl} + \mu \left(\delta _{ik} \delta _{jl} +  \delta _{il} \delta _{kj} \right) </math>
|}
| style="width: 5px;text-align: right;" | (57)
|}

<span id="eq-58"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>C^{*}_{ijkl}= \frac{1}{2} \left(\delta _{ik} \sigma _{jl} +  \delta _{il} \sigma _{jk} +  \delta _{jk} \sigma _{il} + \delta _{jl} \sigma _{ik} \right) - \delta _{kl} \sigma _{ij}  </math>
|}
| style="width: 5px;text-align: right;" | (58)
|}

where <math display="inline">{\lambda }</math> and <math display="inline">{\mu }</math> are the Lamé constants.

Eqs.([[#eq-54|54]]-[[#eq-58|58]]) show that, contrary to the Truesdell measure, the Jaumann tangent moduli tensor <math display="inline">{C^{\sigma J}}</math> is symmetric. This feature represents the main advantage of the Jaumann measure and it explains why it is the largest used measure of the Cauchy stress rate. However, the Jaumann stress rate represents an approximation of the Truesdell's one. The solution given by the latter is more accurate especially for shear-dominant problems <span id='citeF-3'></span>[[#cite-3|[3]]].

The material rate for the Cauchy stress tensor can be computed as:

<span id="eq-59"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>\dot{\sigma }=\sigma ^{\bigtriangledown \tau }+ \Omega ^{\bigtriangledown \tau }  </math>
|}
| style="width: 5px;text-align: right;" | (59)
|}

<span id="eq-60"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>\dot{\sigma }=\sigma ^{\bigtriangledown J}+ \Omega ^{\bigtriangledown J}  </math>
|}
| style="width: 5px;text-align: right;" | (60)
|}

where <math display="inline">{\Omega ^{\bigtriangledown \tau }}</math> and <math display="inline">{\Omega ^{\bigtriangledown J}}</math> are the part of the material rate of the Cauchy stress due to the rotations for the Truesdell and the Jaumann measures, respectively. They are defined as follows:

<span id="eq-61"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>\Omega ^{\bigtriangledown \tau }=-W\cdot \sigma  -\sigma \cdot W^T    </math>
|}
| style="width: 5px;text-align: right;" | (61)
|}

<span id="eq-62"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>\Omega ^{\bigtriangledown J}=W\cdot \sigma   +\sigma \cdot W^T    </math>
|}
| style="width: 5px;text-align: right;" | (62)
|}

where <math display="inline">{ W}</math> is the spin tensor defined as

<span id="eq-63"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>W_{ij}=\frac{1}{2} \left(L_{ij}-L_{ji}\right)= \frac{1}{2} \left(\frac{\partial v_i }{\partial x_j }-\frac{\partial v_j }{\partial x_i }\right) </math>
|}
| style="width: 5px;text-align: right;" | (63)
|}

Discretizing in time  and introducing the constitutive laws into Eqs.([[#eq-59|59]] and [[#eq-60|60]]), for the time step increment <math display="inline">{}</math>  the following relations hold

<span id="eq-64"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>\frac{{^{n+1}\sigma }- {^{n}\sigma }}{\Delta t}= {^{n+1}C}^{\sigma \tau } : {^{n+1}D} -{^{n+1}W}\cdot{^{n+1}\sigma }-{^{n+1}\sigma }\cdot{^{n+1} W^T}    </math>
|}
| style="width: 5px;text-align: right;" | (64)
|}

<span id="eq-65"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>\frac{ {^{n+1}\sigma }- {^{n}\sigma }}{\Delta t}= C^{\sigma J} : {^{n+1}D} + {^{n+1}W}\cdot {^{n+1}\sigma }+{^{n+1}\sigma }\cdot{^{n+1} W}^T    </math>
|}
| style="width: 5px;text-align: right;" | (65)
|}

Relations ([[#eq-64|64]]) and ([[#eq-65|65]]) are not linear and the Cauchy stress tensor can not be computed explicitly. Hence, the stress tensor has to computed iteratively within the loop of the incremental solution scheme  described in Box 1. Once the Cauchy stress tensor has been computed, it has to be introduced into the geometric part of the tangent matrix ([[#eq-36|36]])  and into the residual term ([[#eq-40|40]]). Concerning the material part ([[#eq-30|30]]), the tangent moduli tensor has to be replaced at each iteration only if the Truesdell stress measure is used. In fact, the Jaumann tangent moduli tensor does not change with time and it can be computed only once at the beginning of the analysis.

In Box 3, the iterative solution incremental solution scheme for hypoelastic solids using the velocity formulation  is described for a generic time increment <math display="inline">{}</math>. 

<math>{Mixed  velocity - pressure   formulation}</math>

In order to solve an hypoelastic solid with the mixed velocity-pressure formulation, the modifications over the velocity scheme described in Section [[#3 Mixed velocity-pressure formulation|3]] have to be introduced.

The deviatoric part of the Truesdell and Jaumann stress rate measure are:

<span id="eq-66"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>\sigma '^{\bigtriangledown \tau }=C^{\sigma '\tau } : D^{dev}  </math>
|}
| style="width: 5px;text-align: right;" | (66)
|}

<span id="eq-67"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>\sigma '^{\bigtriangledown J}=C^{\sigma ' J} : D^{dev}  </math>
|}
| style="width: 5px;text-align: right;" | (67)
|}

where <math display="inline">{C^{\sigma ' \tau }}</math> and <math display="inline">{C^{\sigma ' J}}</math> are defined as

<span id="eq-68"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>C^{\sigma ' \tau }=C^{\sigma 'J} - C'^*  </math>
|}
| style="width: 5px;text-align: right;" | (68)
|}

with:

<span id="eq-69"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>C^{\sigma ' J}_{ijkl}= \mu \left(\delta _{ik} \delta _{jl} +  \delta _{il} \delta _{kj} \right) </math>
|}
| style="width: 5px;text-align: right;" | (69)
|}

<span id="eq-70"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>{C'}^{*}_{ijkl}= \frac{1}{2} \left(\delta _{ik} \sigma '_{jl} +  \delta _{il} \sigma '_{jk} +  \delta _{jk} \sigma '_{il} + \delta _{jl} \sigma '_{ik} \right) - \delta _{kl} \sigma '_{ij}  </math>
|}
| style="width: 5px;text-align: right;" | (70)
|}

Using the time discretization described in Eqs.([[#eq-59|59]]-[[#eq-65|65]]) the material rate for the Cauchy stress tensor can be computed as:

<span id="eq-71"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>\frac{{^{n+1}\sigma '}- {^{n}\sigma '}}{\Delta t}= {^{n+1}C}^{\sigma ' \tau } : {^{n+1}D^{dev}} +\Omega '^{\bigtriangledown \tau }  </math>
|}
| style="width: 5px;text-align: right;" | (71)
|}

<span id="eq-72"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>\frac{ {^{n+1}\sigma '}- {^{n}\sigma '}}{\Delta t}= C^{\sigma ' J} : {^{n+1}D^{dev}} +\Omega '^{\bigtriangledown J}    </math>
|}
| style="width: 5px;text-align: right;" | (72)
|}

with:

<span id="eq-73"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>\Omega '^{\bigtriangledown \tau } = -{^{n+1}W}\cdot{^{n+1}\sigma '}-{^{n+1}\sigma '}\cdot{^{n+1} W^T}    </math>
|}
| style="width: 5px;text-align: right;" | (73)
|}

<span id="eq-74"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>\Omega '^{\bigtriangledown J} = {^{n+1}W}\cdot {^{n+1}\sigma '}+{^{n+1}\sigma '}\cdot{^{n+1} W}^T    </math>
|}
| style="width: 5px;text-align: right;" | (74)
|}

As for Eqs. ([[#eq-64|64]]) and ([[#eq-65|65]]), the computation of the deviatoric part of the Cauchy stress tensor  (Eqs. ([[#eq-71|71]]) and ([[#eq-72|72]])) is performed iteratively within the global solution scheme.

In Box 4, the iterative solution incremental scheme for hypoelastic solids using the mixed velocity-pressure formulation  is described for a generic time increment <math display="inline">{}</math>.

Note that the linear momentum and the continuity equations can be easily decoupled. For this purpose, the residual of the linear momentum equations is written in terms of the Cauchy stress and this tensor is computed using the velocity only and not as the sum of its deviatoric part and the pressure. In this case, a velocity formulation is recovered because the continuity equation is used to compute the pressures only.

===4.2 Quasi-incompressible Newtonian fluid===

The formulation is here adapted to the case of quasi-incompressible Newtonian fluids. Considering a quasi-incompressible fluid is equivalent to solving a Navier-Stokes problem where the divergence of the velocity is not required to vanish in the continuity equation. At the same time, this compressibility is such small that the variation with time of the density is considered null as for a fully incompressible material. As already mentioned, in order to avoid the numerical instabilities, for this type of material only the mixed velocity-pressure formulation can be used with confidence. In this work, the fluid is discretized using the Particle Finite Element Method (PFEM) (www.cimne.com/pfem) <span id='citeF-5'></span><span id='citeF-17'></span><span id='citeF-19'></span>[[#cite-5|[5,17,19]]]. The essential features of the PFEM will be given in Section [[#4.2.3 About the Particle Finite Element Method (PFEM)|4.2.3]]. The same linear interpolation has been used for the velocity and pressure fields. It is well known that this combination does not fulfill the <math display="inline">{inf-sup}</math> condition <span id='citeF-4'></span>[[#cite-4|[4]]] and the solution needs to be stabilized. In this work, the required stabilization is provided by the Finite Calculus (FIC) technique presented in <span id='citeF-22'></span>[[#cite-22|[22]]]. The stabilization terms only affect the continuity equation. In conclusion, the momentum equation is modified only by introducing the Newtonian constitutive law into the internal component of the tangent matrix.

First of all the constitutive law for a Newtonian fluid is recalled as

<span id="eq-75"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>\sigma =\sigma ' - pI =2 \mu D^{dev} - pI </math>
|}
| style="width: 5px;text-align: right;" | (75)
|}

where <math display="inline">{\mu }</math> is the fluid viscosity.

For Newtonian fluids the governing equations are Eqs.([[#eq-1|1]]) and ([[#eq-44|44]]). Note that, if a fully-incompressible material is considered, <math display="inline">{\kappa \rightarrow \infty }</math> and the standard form of the incompressibility equation is recovered from  Eq.([[#eq-44|44]])  (<math display="inline">{ {D}^{v}}=0</math>).

Considering a time interval <math display="inline">{[n,n+1]}</math> and substituting Eq.([[#eq-47|47]]) into Eq.([[#eq-41|41]]) yields:

<span id="eq-76"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>^{n+1}\sigma = \left(2 \mu I^{dev} - \Delta t \kappa I\otimes I  \right): D - ^{n}p I  </math>
|}
| style="width: 5px;text-align: right;" | (76)
|}

where:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>I^{dev}=I-\frac{1}{3}I\otimes I </math>
|}
| style="width: 5px;text-align: right;" | (77)
|}

For convenience, Eq.([[#eq-76|76]]) is rewritten in the following form:

<span id="eq-78"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>^{n+1}\Delta \sigma =^{n+1}\sigma  - ^{n}\sigma =\left( C^{d} +  C^{\kappa } \right): D </math>
|}
| style="width: 5px;text-align: right;" | (78)
|}

where the following substitutions have been done:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>^{n}\sigma _{ij}= - ^{n}p I </math>
|}
| style="width: 5px;text-align: right;" | (79)
|}

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>C^{d}=2 \mu I^{dev} </math>
|}
| style="width: 5px;text-align: right;" | (80)
|}

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>C^{\kappa }= - \Delta t \kappa I\otimes I </math>
|}
| style="width: 5px;text-align: right;" | (81)
|}

The goal is to obtain a relation between the measure of rate of stress and the rate of deformation similar to the one introduced in the material part of the tangent matrix (Eq.([[#eq-29|29]])). As shown in Eq.([[#eq-78|78]]), in fluids the rate of deformation is related through the constitutive parameter to Cauchy stress and not to rate of the Cauchy stress, as for Eq.([[#eq-29|29]]). For this reason, in fluids the concept of objectivity of the stress rate measures is not as important as for solids. Rigid rotations do not cause any state of stress, because the Cauchy stress is expressed in term of the rate of deformation. For these reasons, the rate of Cauchy stress can be simply defined using Eq.([[#eq-78|78]]) as:

<span id="eq-82"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>\dot \sigma =\frac{\Delta  \sigma }{\Delta t}= \left(\frac{ C^{d}}{\Delta t} + \frac{C^{\kappa }}{\Delta t}   \right): D  </math>
|}
| style="width: 5px;text-align: right;" | (82)
|}

Note that  Eq.([[#eq-82|82]]) has the same structure as Eq.([[#eq-29|29]]). For convenience, the constitutive matrix has been split in the deviatoric and volumetric part. As a consequence, for a quasi incompressible Newtonian fluid, the material part of the tangent matrix is written as the sum of two submatrices defined as follows:

<span id="eq-83"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>K^{m}_{NF}=K^{\mu } +K^{\kappa }  </math>
|}
| style="width: 5px;text-align: right;" | (83)
|}

where <math display="inline">{K^{\mu }}</math> and <math display="inline">{K^{\kappa }}</math> are defined for a generic finite element <math display="inline">{e}</math> and the pair of nodes <math display="inline">{I,J}</math> as:

<span id="eq-84"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>K^{\mu }_{IJ}=\int _{\Omega ^e} B^{eT}_I  C^{\mu }  B^e_J    d\Omega  </math>
|}
| style="width: 5px;text-align: right;" | (84)
|}

<span id="eq-85"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>K^{\kappa }_{IJ}=\int _{\Omega ^e} B^{eT}_I  m \kappa  m^T B^e_J    d\Omega  </math>
|}
| style="width: 5px;text-align: right;" | (85)
|}

<math display="inline">\begin{array}{l} \\   \hbox{with}       {C^{\mu }}=\mu \left[                      \begin{array}{cccccc}                        2/3 & -1/3 & -1/3 & 0 & 0 & 0 \\                         &   2/3 & -1/3 &0  &0  & 0 \\                         &  &  2/3 & 0 & 0 & 0 \\                         &  &  & 1 & 0 & 0 \\                        {Sym.} &  &  &  & 1 & 0 \\                         &  &  &  &  & 1 \\                      \end{array} \right] \\ \\ \end{array}</math>

The volumetric part of the tangent matrix can compromise the conditioning of the linear system because its terms are orders of magnitude larger than the viscous part. In order to prevent the numerical instabilities originated by the ill-conditioning of the tangent matrix, a reduced pseudo-bulk modulus can be used in the expression of <math display="inline">{K^{\kappa }}</math> without altering the numerical results <span id='citeF-12'></span>[[#cite-12|[12]]].

Concerning the material part of the tangent matrix, the form of Eq.([[#eq-36|36]]) holds also for quasi-incompressible Newtonian fluids. Now the Cauchy stress tensor should be simply replaced with the expression obtained using Eq.([[#eq-76|76]]).

====4.2.1 Stabilization procedure using FIC====

The interpolation orders of the velocity and pressure fields do not fullfil the so-called <math display="inline">{LBB  inf-sup  condition}</math>  <span id='citeF-4'></span>[[#cite-4|[4]]]. Consequently the numerical scheme needs to be stabilized. The required stabilization is introduced via  the Finite Calculus (FIC) technique presented in <span id='citeF-22'></span>[[#cite-22|[22]]]. In the mentioned work, a FIC-based stabilized finite element formulation for quasi-incompressible Lagrangian fluids is presented and successfully applied to the analysis of several free-surface flow problems. The formulation has excellent mass preservation features. The details of this stabilized formulation lie outside the objectives of this work and can be found in  <span id='citeF-22'></span>[[#cite-22|[22]]]. As already mentioned, this stabilization technique only affects the continuity equation (Eq.([[#eq-44|44]])) and not the linear momentum equations  (Eq.([[#eq-9|9]])).

The fully-discretized form of the stabilized mass conservation equation reads

<span id="eq-86"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>\left(\frac{1}{ \Delta t} \hbox{M}_{1} +   \frac{1}{ {\Delta t}^2} \hbox{M}_{2} + {L}+{M}_b \right){^{n+1}\bar {p}} = {Q}^T {^{n+1}\bar{v}}+{^ng}^S </math>
|}
| style="width: 5px;text-align: right;" | (86)
|}

where:

<span id="eq-87"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>^n{g}^S = {f}_p + \frac{1}{ \Delta t} \hbox{M}_{1} {^{n} \bar {p}} + \frac{1}{ \Delta t} \hbox{M}_{2}  \left(\frac{ ^{n}\bar {p}}{\Delta t} + {^n \dot{\bar{p}}} \right)  </math>
|}
| style="width: 5px;text-align: right;" | (87)
|}

with:


{| class="wikitable" style="text-align: left; margin: 1em auto;"
|-
| <math>\begin{array}{l} \displaystyle \displaystyle   {M}_{2_{IJ}} =\int _{\Omega ^e} \frac{\tau }{\kappa }N_I N_J d\Omega ~ {where}    \displaystyle ~\tau = \left(\frac{8\mu }{h^2}+ \frac{2\rho }{\delta } \right)^{-1}\\[.4cm] \displaystyle{M}_{b_{IJ}}=\int _{\Gamma _t} \frac{2\tau }{h_N} N_I N_J d\Gamma ~    ;   ~ L_{IJ}= \int _{\Omega ^e} \tau \left(\bigtriangledown ^T N_I\right)\bigtriangledown N_J d\Omega  ;\\[.4cm] \displaystyle f_{p_{I}}=\int _{\Gamma _t}q\tau N_I \left[\frac{Dv_N}{Dt}-\frac{2}{h_N} (2\mu D_N -t_N)\right]d\Gamma - \int _{\Omega ^e} \tau  {\bigtriangledown ^TN_I{b}}d\Omega  \end{array}</math>


|}

where <math display="inline">{\tau }</math> is the stabilization parameters and <math display="inline">{h}</math> and <math display="inline">{\delta }</math> are characteristic distances in space and time, respectively. In practice, <math display="inline">{h}</math> and <math display="inline">{\delta }</math> have the order of magnitude of the element size and the time step increment, respectively. The subscripts <math display="inline">{(\cdot )_N}</math> denote the normal component of the variable.

====4.2.2 Incremental solution scheme====

For treating Newtonian quasi-incompresible fluids, the iterative scheme for the mixed formulation described in Box 2 has to be modified slightly for taking into account the features of the material described at the beginning of this section. In particular, the material part tangent matrix in the momentum equations has two contributions (Eqs.([[#eq-84|84]]) and ([[#eq-85|85]])), the continuity equation must be considered in its stabilized form (Eq.([[#eq-86|86]])) and, obviously, the stress measures are computed using the Newtonian constitutive law (Eq.([[#eq-41|41]])).

The incremental iterative solution scheme for a quasi incompressible Newtonian fluid is described for a generic time increment <math display="inline">{}</math> in Box 5.

====4.2.3 About the Particle Finite Element Method (PFEM)====

The Particle Finite Element Method is a particular class of Lagrangian finite element formulation. Its efficacy has been proved in many scientific publications, as  <span id='citeF-17'></span>[[#cite-17|[17]]], <span id='citeF-21'></span>[[#cite-21|[21]]]  and many others. PFEM treats the mesh nodes in the fluid and solid domains as particles which can freely move and even separate from the main fluid domain representing, for instance, the effect of water drops. A mesh connects the nodes discretizing the domain where the governing equations are solved using a stabilized FEM.  These characteristic features makes this method the ideal instrument to model and simulate free surface flows.

Consider a global domain containing a fluid and a solid subdomains. In the PFEM both subdomains are modelled using an ''updated'' ''Lagrangian formulation''. The finite element method (FEM) is used to solve the equations of continuum mechanics for each of the subdomains. Hence a mesh discretizing these domains is generated in order to solve the governing equations for each subdomain in the standard FEM fashion.

For clarity purposes,  the '' collection or cloud of nodes'' pertaining to the fluid and solid domains is defined'' C'', the ''volume'' defining the analysis domain for the fluid and the solid is termed''V''and the ''mesh'' discretizing both domains  is called ''M''.

The solution steps within a time step are as follows:

<div id='img-1'></div>
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:draft_Samper_438761226-CloudPoints.png|600px|Sequence of steps to update a “cloud” of nodes representing a domain containing a fluid and a solid part from time n   (t=ⁿt)  to   time n+2 (t=ⁿt +2∆t)  ]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 1:''' Sequence of steps to update a “cloud” of nodes representing a domain containing a fluid and a solid part from time <math>n</math>   (<math>t=^n t</math>)  to   time <math>n+2</math> (<math>t=^n t +2\Delta t</math>)  
|}

<ol>

<li>The starting point at each time step is the cloud of points in the fluid and solid domains. For instance <math display="inline">^{n} C</math> denotes the cloud at time <math display="inline">t=^n t </math> Fig.([[#img-1|1]]). </li>

<li>Identify the boundaries of both the fluid and solid domains defining the analysis domain <math display="inline">^{n} V</math> in the fluid and the solid. This is an essential step as some boundaries (such as the free surface in fluids) may be severely distorted during the solution, including separation and re-entering of nodes. The Alpha Shape method <span id='citeF-10'></span>[[#cite-10|[10]]] is used for the boundary definition. </li>

<li> Discretize the fluid and solid domains with a finite element mesh<math display="inline">^{n} M.</math> In this work, an efficient mesh generation scheme based on the extended Delaunay tesselation <span id='citeF-15'></span>[[#cite-15|[15]]] is used. </li>

<li> Solve the Lagrangian equations of motion for the overall continuum. Compute the state variables in at the next (updated) configuration for <math display="inline">t+\Delta t</math>: velocities, pressure and viscous stresses in the fluid and displacements, stresses and strains in the solid. </li>

<li>  Move the mesh nodes to a new position <math display="inline">^{n+1} C</math> where ''n''+1 denotes the time <math display="inline">t^n +\Delta t</math>, in terms of the time increment size. This step is typically a consequence of the solution process of step 4. </li>

<li> Go back to step 1 and repeat the solution for the next time step to obtain<math display="inline">^{n+2} C</math>  </li>

</ol>

It is to note that the key differences between the PFEM and the classical FEM are the remeshing technique and the identification of the domain boundary at each time step.

The quality of the numerical solution depends on the discretization chosen as in the standard FEM. Adaptive mesh refinement techniques can be used to improve the solution.

==5 Coupling strategy for fluid-structure interaction (FSI) analysis==

In this section the coupling strategy for solving FSI problems using the unified formulation is described and tested. The algorithm uses the same mixed velocity-pressure formulations for hypoelastic solids and Newtonian fluids. As explained in the previous sections, in the unified formulation solids and fluids are treated practically in the same way. Apart from the material parameters, the analysis of  newtonian fluids and hypoelastic solids with the PFEM formulation here presented differs only for the stabilization and the remeshing step that are implemented for the fluid only. Solids and fluids have many common points: they use the same set of the same unknown variables (velocities and pressure), the same framework (Lagrangian) and the same solution scheme (Gauss-Seidel). All these common features simplify the solution of the FSI problems with a monolithic scheme. In summary fluids and solids can be easily solved within a unique system of equations.

The main advantage of the monolithic scheme versus the staggered approach, is that it ensures a strong coupling automatically. In fact in the staggered schemes, a few iterations may be necessary in order to obtain a similar degree of coupling as for a monolithic scheme.

On the other hand, a monolithic scheme can lead to ill-conditioned linear systems because the order of magnitude of the terms corresponding to the fluid and the solid parts of the domain may be very different. In this work, this drawback is overcome via the segregated scheme and the unified formulation. In fact, the partitioning of the governing equations separates the velocity and pressure unknowns and the matrices referred to them. In this way the size of the matrices is reduced and its conditioning is improved as the risk of mixing in the same matrix terms with different order of  magnitude is prevented. Moreover, in the unified formulation the solid and the fluid use the same nodal unknowns. As consequence, the fluid and solid matrices terms differ only in the material constants and in for the variable transformation.

By using  the unified formulation and the PFEM, FSI problems can be solved in a natural way. It essentially consists in performing two tasks: to detect the interface for creating the coupling elements (elements that share at least one interface node) and to assembly properly the global system.

As explained in the previous sections, the classical FEM is used for the solid while the PFEM is adopted for the fluid. Hence, the solid domain keeps the same grid during all over the analysis, whereas a remeshing of the fluid domain is performed whenever its discretization becomes excessively distorted as it is typically done in the PFEM. The nodes over the boundaries of the solid are the interface nodes, and they may be part of the fluid domain too. The solid boundaries are detected by exploiting the capability of the PFEM to match contours. In practise, the fluid domain detect the solid interface in the same way it recognizes its fixed contours. This represents one of the key features of the PFEM. The technique uses an extended Delaunay partition for recognizing boundary nodes <span id='citeF-15'></span>[[#cite-15|[15]]]. For all the nodes a characteristic length ''h'', typically the minimum distance between two nodes, is assigned. All the nodes on an empty sphere with a radius greater than a critical value are considered as boundary nodes. The critical value is  defined as <math display="inline">{h_{crit}=\alpha h}</math>, where <math display="inline">\alpha </math> is a parameter close to one (typically, values of <math display="inline">\alpha </math> ranging between 1.3 and 1.5 are chosen). This criterion is coincident with the Alpha Shape concept <span id='citeF-10'></span>[[#cite-10|[10]]]. In this way, the fluid boundaries are recognized exactly. In the case of FSI problems, this test is realized also on the solid interface particles. If the fluid domain is sufficiently close to the interface, at least, one interface node will have the variable ''h'' smaller than <math display="inline">{h_{crit}}</math>. If this occurs, an element joining the fluid domain and the solid interface will be built, whereas the two domains keep separated. Fig. [[#img-2|2]] shows a graphic representation of this technique. In order to avoid that particles of fluid pass through the solid boundary, the solid discretization in the region close to the interface should have a size similar to that used for the fluid mesh.

<div id='img-2'></div>
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 80%;max-width: 100%;"
|-
|[[Image:draft_Samper_438761226-pfemContact.png|480px|Detection of the interface using PFEM]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 2:''' Detection of the interface using PFEM
|}

When at least one coupling element is created, there will be some nodes on the interface that belong to both the solid and the fluid domains. In regard to the DOFs of these nodes, the velocities are the same as for the fluid and the solid (in this way the Dirichlet  condition is prescribed strongly) while the pressure DOFs must be duplicated. In fact, along the interface the Neumann transmission condition has to be imposed (in weak form) on the normal component of the Cauchy stress and not on the pressure. Hence only for the momentum equations, the interface nodes contribute to the global linear system the sum of the contributions of the solid and the fluid elements they are sharing, as represented in Fig.[[#img-3|3]].

<div id='img-3'></div>
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:draft_Samper_438761226-FSImatrix.png|600px|Graphic representation of domain contributions to the momentum equations global system]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 3:''' Graphic representation of domain contributions to the momentum equations global system
|}

==6 Numerical examples==

===6.1 Collapse of a water column on a deformable elastic membrane===

The problem is illustrated in Fig.([[#img-4|4]]) and it was introduced by Walhorn <math display="inline">{et  al.}</math> in <span id='citeF-27'></span>[[#cite-27|[27]]]. The water column collapses by instantaneously removing the vertical wall retaining the water. This originates the flow of water within the tank, the formation of a jet after the water stream hits the rigid step and the subsequent sloshing of the fluid as it impacts a highly deformable elastic membrane. The membrane bends and starts oscillating under the effect of its inertial forces and the impact with the water stream.

<div id='img-4'></div>
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:draft_Samper_438761226-FSIinput.png|600px|Collapse of a water column on a deformable membrane: starting geometry and problem data.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 4:''' Collapse of a water column on a deformable membrane: starting geometry and problem data.
|}

In Figs.([[#img-5|5]],[[#img-6|6]]) some representative snapshots of 2D and 3D numerical simulation results are illustrated. The results obtained with the present formulation have been compared with the ones computed by other formulations presented in scientific publications <span id='citeF-13'></span>[[#cite-13|[13]]], <span id='citeF-14'></span>[[#cite-14|[14]]] and <span id='citeF-8'></span>[[#cite-8|[8]]]. In the graph of Fig.([[#img-7|7]]) the time evolution of the horizontal deflection of the right top corner is illustrated. The diagram shows that the present formulation agrees well with the results found in the literature. In the graph of Fig.([[#img-8|8]]) the numerical results obtained with the 2D and 3D simulations are compared.

<div id='img-5a'></div>
<div id='img-5b'></div>
<div id='img-5c'></div>
<div id='img-5d'></div>
<div id='img-5'></div>
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:draft_Samper_438761226-fsi2D0235.png|400px|t = 0.235 st = 0.365 s]]
|[[Image:draft_Samper_438761226-fsi2D0365.png|400px|t = 0.515 s]]
|- style="text-align: center; font-size: 75%;"
| (a) t = 0.235 s
| (b) t = 0.515 s
|-
|[[Image:draft_Samper_438761226-fsi2D0515.png|400px|t = 0.595s]]
|[[Image:draft_Samper_438761226-fsi2D0595.png|400px|Collapse of a water column on a deformable membrane: snapshots at different instants of the 2D simulation.]]
|- style="text-align: center; font-size: 75%;"
| (c) t = 0.595s
|- style="text-align: center; font-size: 75%;"
| colspan="2" | '''Figure 5:''' Collapse of a water column on a deformable membrane: snapshots at different instants of the 2D simulation.
|}

<div id='img-6'></div>
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:draft_Samper_438761226-FSIDB3B023.png|400px|t = 0.23 st = 0.37 s]]
|[[Image:draft_Samper_438761226-FSIDB3B037.png|400px|t = 0.49 s]]
|- style="text-align: center; font-size: 75%;"
| t = 0.23 s
| t = 0.49 s
|-
|[[Image:draft_Samper_438761226-FSIDB3B049.png|400px|t = 0.89 s]]
|[[Image:draft_Samper_438761226-FSIDB3B089.png|400px|Collapse of a water column on a deformable membrane: snapshots of different instants of the 3D simulation.]]
|- style="text-align: center; font-size: 75%;"
| t = 0.89 s
|- style="text-align: center; font-size: 75%;"
| colspan="2" | '''Figure 6:''' Collapse of a water column on a deformable membrane: snapshots of different instants of the 3D simulation.
|}

<div id='img-7'></div>
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:draft_Samper_438761226-FSI2Dcomparison.png|600px|Collapse of a water column on a deformable membrane: horizontal deflection of the right top corner on time. Comparison with numerical results obtained by other formulations. Curves  '1', '2', '3' correspond to <span id='citeF-27'></span>[[#cite-27|[27]]], <span id='citeF-14'></span>[[#cite-14|[14]]] and <span id='citeF-8'></span>[[#cite-8|[8]]] respectively]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 7:''' Collapse of a water column on a deformable membrane: horizontal deflection of the right top corner on time. Comparison with numerical results obtained by other formulations. Curves  '1', '2', '3' correspond to <span id='citeF-27'></span>[[#cite-27|[27]]], <span id='citeF-14'></span>[[#cite-14|[14]]] and <span id='citeF-8'></span>[[#cite-8|[8]]] respectively
|}

<div id='img-8'></div>
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 85%;max-width: 100%;"
|-
|[[Image:draft_Samper_438761226-FSIdbComparison23D.png|510px|Collapse of a water column on a deformable membrane: horizontal deflection of the right top corner on time. Comparison between 2D and 3D analyses.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 8:''' Collapse of a water column on a deformable membrane: horizontal deflection of the right top corner on time. Comparison between 2D and 3D analyses.
|}

===6.2 Water entry of a nylon sphere===

The problem was presented by Aristoff <math display="inline">{et  al.}</math> in <span id='citeF-1'></span>[[#cite-1|[1]]]. In the mentioned work the experimental results of the water entry of spheres of different materials are studied. In this section, the case of a nylon sphere is analyzed. The numerical results given by the unified formulation are compared with the results of the laboratory test. The sphere has a diameter of 2.54 cm and impacts the water in the tank with a vertical velocity of 2.17 <math display="inline">{m/s}</math>. The density of nylon is 1140 <math display="inline">{kg/m^3}</math> and its Young modulus and Poisson coefficient are 3 <math display="inline">{GPa}</math> and 0.2, respectively. Water has been simulated considering a density of 1000 <math display="inline">{kg/m^3}</math>, a dynamic viscosity 0.00089 <math display="inline">{Pa \cdot s}</math> and bulk modulus 2.15 <math display="inline">{GPa}</math>. In order to simulate  this problem correctly, a very fine mesh was necessary. For this reason the whole domain has been discretized using 1059924 tetrahedra. The time step used for the analysis is <math display="inline">{\Delta t =10^{-4} s}</math>.

In Fig.([[#img-9|9]]) the numerical results are compared with the experimental ones published in  <span id='citeF-1'></span>[[#cite-1|[1]]]. The comparison shows the good agreement of the results obtained with the present formulation with the laboratory experience.

<div id='img-9'></div>
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:draft_Samper_438761226-aristoffReal1.png|400px|]]
|[[Image:draft_Samper_438761226-aristoffNumerical1.png|400px|]]
|-
|[[Image:draft_Samper_438761226-aristoffReal2.png|400px|]]
|[[Image:draft_Samper_438761226-aristoffNumerical2.png|400px|Water entry of a nylon sphere: comparison between experimental and numerical results.]]
|- style="text-align: center; font-size: 75%;"
| colspan="2" | '''Figure 9:''' Water entry of a nylon sphere: comparison between experimental and numerical results.
|}

==7 Conclusions==

A unified updated Lagrangian finite element formulation for solving FSI problems has been derived, discussed and tested. The method belongs to the monolithic strategies for FSI problems because both the solid and the fluid domains are solved with a general system of equations. First a velocity formulation for a general continuum has been derived. After that, the mixed velocity-pressure formulation has been obtained. After, both formulations have been adapted for specific materials. In this paper only hypoelastic solids and quasi incompressible fluids have been studied. However the formulation holds for other type of material, such as incompressible and elasto-plastic solids, among others.

It has been shown that the essential differences in the treatise of the fluid and the solid domain are the different meshing techniques and the inclusion of the stabilization terms in the mass balance equation of the fluid. In this work, the stabilized formulation emanates from the Finite Calculus (FIC) technique presented in <span id='citeF-22'></span>[[#cite-22|[22]]]. Concerning the meshing, the Particle Finite Element Method is adopted for the fluid analysis, while for the discretization of the solid domain the classical FEM is used.

With the unified mixed formulation, both fluid and solid domains are solved in a Lagrangian framework via a Gauss-Seidel solution scheme using  the velocities and the pressure as the unknown variables. All these common points facilitate the coupling for solving FSI problems. In fact, it has been shown that the coupling essentially consists in detecting the interface between the fluid and the solid domains and assembling the global system. It has been shown the importance of the PFEM for the former task.

Finally, the unified formulation has been tested solving different benchmark problems. In particular, the collapse of a water column on a deformable membrane and the water entry of a nylon sphere have been numerically studied in 2D and 3D. The comparison of the numerical results with the ones published in the literature and with scientific tests is fully satisfactory and evidences the efficiency and the accuracy of the proposed method.

===BIBLIOGRAPHY===

<div id="cite-1"></div>
'''[[#citeF-1|[1]]]'''  J.&nbsp;Aristoff, T.T. Truscott, A.H. Techet, and J.W.M Bush. The water entry of decelerating spheres. ''Physics of Fluids'', 22:032102, 2010.

<div id="cite-2"></div>
'''[[#citeF-2|[2]]]'''  K.J. Bathe. ''Finite Element Procedures''. Prentice-Hall, New Jersey, 1996.

<div id="cite-3"></div>
'''[[#citeF-3|[3]]]'''  T.&nbsp;Belytschko, W.K. Liu, and B.&nbsp;Moran. ''Nonlinear Finite Elements For Continua And Structures''. John Wiley & Sons, New York, 2000.

<div id="cite-4"></div>
'''[[#citeF-4|[4]]]'''  F.&nbsp;Brezzi. On the existence, uniqueness and approximation of saddle-point   problems arising from lagrange multipliers. ''Revue francaise d'automatique, informatique, recherche   opérationnelle. Série rouge. Analyse numérique'', 8(R-2):129&#8211;151,   November 1974.

<div id="cite-5"></div>
'''[[#citeF-5|[5]]]'''  J.M. Carbonell, E.&nbsp;Oñate, and B.&nbsp;Suarez. Modeling of ground excavation with the particle finite-element   method. ''Journal of Engineering Mechanics'', 136:455&#8211;463, 2010.

<div id="cite-6"></div>
'''[[#citeF-6|[6]]]'''  M.&nbsp;Cervera, M.&nbsp;Chiumenti, and R.&nbsp;Codina. Mixed stabilized finite element methods in nonlinear solid mechanics:   Part i: Formulation. ''Computer Methods In Applied Mechanics And Engineering'',   199:2559&#8211;2570, 2010.

<div id="cite-7"></div>
'''[[#citeF-7|[7]]]'''  M.&nbsp;Chiumenti, Q.&nbsp;Valverde, C.&nbsp;Agelet&nbsp;De Saracibar, and M.Cervera. A stabilized formulation for incompressible elasticity using linear   displacement and pressure interpolations. ''Computer Methods In Applied Mechanics And Engineering'',   191:5253&#8211;5264, 2002.

<div id="cite-8"></div>
'''[[#citeF-8|[8]]]'''  M.&nbsp;Cremonesi. ''Doctoral thesis: A Lagrangian Finite Element Method for the   Interaction Between Flexible Structures and Free Surfaces Fluid Flows''. 2010.

<div id="cite-9"></div>
'''[[#citeF-9|[9]]]'''  M.&nbsp;Cremonesi, A.&nbsp;Frangi, and U.&nbsp;Perego. A lagrangian finite element approach for the analysis of   fluid–structure interaction problems. ''International Journal of Numerical Methods in Engineering'',   84:610&#8211;630, 2010.

<div id="cite-10"></div>
'''[[#citeF-10|[10]]]'''  H.&nbsp;Edelsbrunner and E.P. Mucke. Three dimensional alpha shapes. ''ACM Trans Graphics'', 13:43–72, 1999.

<div id="cite-11"></div>
'''[[#citeF-11|[11]]]'''  C.A. Felippa and K.C. Park. Staggered transient analysis procedures for coupled mechanical   systems: Formulation. ''Computer Methods in Applied Mechanics and Engineering'',   24:61&#8211;111, 1980.

<div id="cite-12"></div>
'''[[#citeF-12|[12]]]'''  A.&nbsp;Franci, E.&nbsp;Oñate, and J.M. Carbonell. On the effect of the tangent bulk stiffness matrix in the analysis of   free surface lagrangian flows using pfem. ''Researh Report CIMNE PI402Computational Mechanics'', Submitted to   International Journal for Numerical Methods in Engineering.

<div id="cite-13"></div>
'''[[#citeF-13|[13]]]'''  B.&nbsp;Hubner, E.&nbsp;Walhorn, and D.Dinkler. A monolithic approach to fluid-structure interaction using space-time   finite elements. ''Computer Methods in Applied Mechanics and Engineering'',   193:2087&#8211;2104, 2004.

<div id="cite-14"></div>
'''[[#citeF-14|[14]]]'''  S.R. Idelshon, J.&nbsp;Marti, A.&nbsp;Limache, and E.&nbsp;Oñate. Unified lagrangian formulation for elastic solids and incompressible   fluids: Applications to fluid-structure interaction problems via the pfem. ''Computer Methods In Applied Mechanics And Engineering'',   197:1762&#8211;1776, February 2008.

<div id="cite-15"></div>
'''[[#citeF-15|[15]]]'''  S.&nbsp;Idelsohn, N.&nbsp;Calvo, and E.&nbsp;Oñate. Polyhedrization of an arbitrary point set. ''Computer Methods in Applied Mechanics and Engineering'',   92(22–24):2649–2668, 2003.

<div id="cite-16"></div>
'''[[#citeF-16|[16]]]'''  S.R. Idelsohn, M.Mier-Torrecilla, and E.&nbsp;Oñate. Multi-fluid flows with the particle finite element method. ''Computer methods in applied mechanics and engineering'',   198:2750&#8211;2767, 2007.

<div id="cite-17"></div>
'''[[#citeF-17|[17]]]'''  S.R. Idelsohn, E.&nbsp;Oñate, and F.&nbsp;Del Pin. The particle finite element method: a powerful tool to solve   incompressible flows with free-surfaces and breaking waves. ''International Journal for Numerical Methods in Engineering'',   61:964&#8211;989, 2004.

<div id="cite-18"></div>
'''[[#citeF-18|[18]]]'''  S.R. Idelsohn, E.&nbsp;Oñate, F.&nbsp;Del Pin, and N.&nbsp;Calvo. Fluid-structure interaction using the particle finite element method. ''Computer methods in applied mechanics and engineering'',   195:2100&#8211;2113, 2006.

<div id="cite-19"></div>
'''[[#citeF-19|[19]]]'''  A.&nbsp;Larese, R.&nbsp;Rossi, E.&nbsp;Oñate, and S.R. Idelsohn. Validation of the particle finite element method (pfem) for   simulation of free surface flows. ''International Journal for Computer-Aided Engineering and   Software'', 25:385&#8211;425, 2008.

<div id="cite-20"></div>
'''[[#citeF-20|[20]]]'''  C.&nbsp;Michler, S.J. Hulshoff, and E.H. Van Brummelenand R.&nbsp;De Borst. A monolithic approach to fluid-structure interaction. ''Computers and Fluids'', 33:839&#8211;848, 2004.

<div id="cite-21"></div>
'''[[#citeF-21|[21]]]'''  E.&nbsp;Oñate, M.A. Celigueta, S.R. Idelsohn, F.&nbsp;Salazar, and B.&nbsp;Suarez. Possibilities of the particle finite element method for   fluid–soil–structure interaction problems. ''Computation mechanics'', 48:307&#8211;318, 2011.

<div id="cite-22"></div>
'''[[#citeF-22|[22]]]'''  E.&nbsp;Oñate, A.&nbsp;Franci, and J.M. Carbonell. Lagrangian formulation for finite element analysis of   quasi-incompressible fluids with reduced mass losses. ''International Journal for Numerical Methods in Fluids'', 74.

<div id="cite-23"></div>
'''[[#citeF-23|[23]]]'''  E.&nbsp;Oñate, A.&nbsp;Franci, and J.M. Carbonell. A particle finite element method for analysis of industrial forming   processes. ''Computational Mechanics'', DOI: 10.1007/s00466-014-1016-2.

<div id="cite-24"></div>
'''[[#citeF-24|[24]]]'''  E.&nbsp;Oñate, J.&nbsp;Rojek, R.L. Taylor, and O.C. Zienkiewicz. Finite calculus formulation for incompressble solids using linear   triangles and tetrahedra. ''International Journal For Numerical Methods In Engineering'',   59:1473&#8211;1500, November 2004.

<div id="cite-25"></div>
'''[[#citeF-25|[25]]]'''  W.&nbsp;Prager. ''Introduction to Mechanics of Continua''. Ginn and Company, Boston, 1961.

<div id="cite-26"></div>
'''[[#citeF-26|[26]]]'''  P.B. Ryzhakov, R.&nbsp;Rossi, S.R. Idelsohn, and E.&nbsp;Oñate. A monolithic lagrangian approach for fluid-structure interaction   problems. ''Computational Mechanics'', 46:883&#8211;899, 2010.

<div id="cite-27"></div>
'''[[#citeF-27|[27]]]'''  E.&nbsp;Walhorn, A.&nbsp;Kolke&nbsp;B. Hubner, and D.Dinkler. Fluid-structure coupling within a monolithic model involving free   surface flows. ''Computer <math>{&}</math> Structures Methods in Applied Mechanics and   Engineering'', 83 (25-26):2100&#8211;2111, 2005.

<div id="cite-28"></div>
'''[[#citeF-28|[28]]]'''  C.&nbsp;Wood, A.J. Gil, O.Hassan, and J.Bonet. Partitioned block gauss-seidel coupling for dynamic fluid-strucure   interaction. ''Computers and Structures'', 88:1367&#8211;1382, 2010.