​

'''A. Franci<math>^1</math>, E. Oñate<math>^{1,2}</math>, J.M. Carbonell<math>^{1,2}</math>'''

​

{|

|-

<math>^1</math> Centre Internacional de Metodes Numerics en Enginyeria (CIMNE)

|-

| Campus Norte UPC, 08034 Barcelona, Spain

|-

| [http://www.cimne.com www.cimne.com], Email: <code>onate,falessandro[mailto:@cimne.upc.edu @cimne.upc.edu]</code> 

|-

| <math>^2</math> Universitat Politecnica de Catalunya (UPC), <code>[mailto:cpuigbo@cimne.upc.edu cpuigbo@cimne.upc.edu]</code>

|}

​

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

​

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

​

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

|}

​

<math>\begin{array}{l} \\ \hbox{with}\\   \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 \\ \\ \end{array}</math>

​

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

​

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