A unified monolithic approach for multi-fluid flows and Fluid-Structure Interaction using the Particle Finite Element Method with fixed mesh

<!-- metadata commented in wiki content
==A unified monolithic approach for multi-fluid flows and Fluid-Structure Interaction using the Particle Finite Element Method with fixed mesh==

'''P. Becker<math>^1</math>,         S.R. Idelsohn<math>^{1,2}</math>, E. Oñate<math>^{1,3}</math> '''
-->

==Abstract==

This paper describes a strategy to solve multi-fluid and Fluid-Structure Interaction (FSI) problems using Lagrangian particles combined with a fixed Finite Element (FE) mesh . Our approach is an extension of the fluid-only PFEM-2 <span id='citeF-1'></span>[[#cite-1|[1]]]<span id='citeF-2'></span>[[#cite-2|[2]]] which uses explicit integration over the streamlines to improve accuracy. As a result, the convective term does not appear in the set of equations solved on the fixed mesh. Enrichments in the pressure field are used to improve the description of the interface between phases.

'''Keywords''': Multi-fluids FSI Fixed mesh Lagrangian particles Unified approach

==1 Introduction==

Simulation of Fluid-Structure Interaction (FSI) problems is an area of growing interest due to the several application fields in which it is required. Nowadays all structures, from buildings to airplanes and boats are optimized to minimize weight, leading to a more flexible behaviour in air or water. Furthermore, new areas are always being explored, biomedical research being a fast-growing and stimulating field.

A common approach to treat FSI problems is the staggered or weakly coupled algorithm. This method consists on computing separately the fluid and the solid domain, without ensuring force balance at each time step. The popularity of this method lies in its speed and possibility to use specialized solvers for the structure and the fluid, a desired property since both subdomains are generally treated in different frameworks. While for the structural solver most algorithms employ a Lagrangian reference framework, for the fluid either Eulerian or Arbitrary Lagrangian-Eulerian (ALE) frameworks are used.

On the other hand, in strongly coupled algorithms the solution is considered to converge only when force balance is achieved at the interfaces. To couple the fluid and structural subdomains, monolithic and partitioned strategies have been developed. While the first method consists on assembling and solving all the equations in a single system, the partitioned algorithm solves separately the fluid and solid subdomains, relying on an iterative process to achieve convergence.

A starting point to reduce the complexity of the problem is selecting the right reference framework. An appropriate choice leads to the simplification of complex terms and higher accuracy in the results due to better approximations. Following this line, the Lagrangian framework offers important advantages due to the simplicity of the equations. However, since material points move and the configuration changes continuously in time, it is necessary to couple this set of equations with a strategy that solves for the movement of material points.   In <span id='citeF-3'></span>[[#cite-3|[3]]] a strongly coupled monolithic strategy was proposed for the PFEM to treat simultaneously the fluid and solid phases within the same framework. Using a Lagrangian approach, it was possible to write the momentum equations for both solids and fluids together and seamlessly, with only minor differences due to the stress memory of solids. This allowed us to solve all the equations within a unified computational frame, significantly improving accuracy. On the other hand, to deal with the large deformations experimented by the fluid phase, a remeshing algorithm was implemented. This way distorted meshes were avoided, but with the extra computational cost associated to the Delaunay triangulation meshing stage.

In this work we propose a new method which uses the same discretization and solving technique for all the phases present in the domain. This strategy is an extension of the Particle Finite Element Method-second generation (PFEM-2) developed for multi-fluids in <span id='citeF-2'></span>[[#cite-2|[2]]]. The PFEM-2 is based on the use of Lagrangian particles with no connectivities to convect all the material properties combined with a fixed mesh, leading to a Lagrangian formulation that does not rely on deforming the spatial mesh nor remeshing.

Since all the material points of the domain are represented by Lagrangian particles, it is possible to convect a large number of variables at very little cost. Furthermore, the convection stage is greatly improved by following the streamlines of the velocity in an explicit way, therefore improving accuracy over first order explicit methods but keeping the computational cost low. Once the particles have been convected, the variables are projected into the mesh, boundaries are identified where material properties change and the Lagrangian system is solved in the fixed mesh. Furthermore, by improving the definition of the interface using enrichment shape functions, the physics of the problem is very closely simulated despite that the mesh does not match the interfaces.

This paper is structured as follows. First, the monolithic strategy of the PFEM-2 is described. In the subsequent sections the strategy for the treatment of the FSI problem is described and finally validation examples are presented.

==2 Monolithic strategy for multi-fluids and FSI==

===2.1 Mixed Framework: The PFEM-2===

The current model is based on a mixed Eulerian-Lagran-gian framework, combining advantages of both methods. This is achieved by using a set of particles combined with a fixed FEM mesh in a solving strategy known as PFEM-2 <span id='citeF-1'></span>[[#cite-1|[1]]]. The method is obtained by writing, for a given particle <math display="inline">p</math>, the following definition for the position and velocity in a given timestep <math display="inline">n+1</math>, where <math display="inline">\mathbf{V}</math> is the velocity and <math display="inline">\mathbf{a}</math> is the acceleration

<span id="eq-1.a"></span>
<span id="eq-1.b"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\mathbf{x}^{n+1}_p = \mathbf{x}^{n}_p + \int _{n}^{n+1} \mathbf{V^t}   dt  </math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (1.a)
|-
| style="text-align: center;" | <math> \mathbf{V}^{n+1}_p = \mathbf{V}^{n}_p + \int _{n}^{n+1} \mathbf{a^t}   dt  </math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (1.b)
|}
|}

By using an explicit strategy in [[#eq-1.a|1.a]], the convective term is completely uncoupled from the momentum equations and Lagrangian equations are obtained. The main advantage of using Lagrangian particles is that the convection is obtained by simply moving the particles across the space and, therefore, the system to be solved does not include the convective term. This set of equations is calculated on the mesh, although it is possible to include partial contributions from the particles to improve accuracy <span id='citeF-1'></span>[[#cite-1|[1]]]. Then the complete scheme required to solve a step becomes:


[[File:Draft_Samper_820077714_7976_fig0.png]]


The particles used in this scheme do not represent a fixed amount of mass, but rather material points with certain properties and velocity. This allows for a variable number of particles per element depending on the zone, to ensure a better accuracy on selected areas. It must be noted that in the algorithm presented in this paper, the particles are only used to transport information (solve the convective term). However, in certain cases where the viscosity is low and there is a single fluid, it is possible to solve partially the momentum equations [[#eq-1.b|1.b]] in the particles, as explained in <span id='citeF-1'></span>[[#cite-1|[1]]]. Even if this strategy leads to higher accuracy in the cases analysed in that article, it lacks the generality that is required for the simulation of FSI problems.

Figure [[#img-1|1]]  shows the streamline integration for a single particle. This step is purely explicit , which is achieved by convecting the particles using the velocity from the previous time step in a Picard <span id='citeF-4'></span>[[#cite-4|[4]]] iteration fashion. As for the force integration (if added), information is also gathered from the last step, so the Lagrangian particle solution step remains purely explicit, allowing for a fast computation on each particle that is trivially parallelizable and, hence, fast.

<div id='img-1'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Samper_820077714-streamline_integration2.png|410px|PFEM-2 streamline integration]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 1:''' PFEM-2 streamline integration
|}

 ''Remark:'' It is possible to imagine this explicit convection as the first step of a non linear iteration process. If an iterative process would be used until convergence, once the new velocity has been calculated, each particle should have its coordinates resetted to the initial position and the streamline integration should be performed again to repeat all the tasks. However, practical experience suggests that the first iteration provides results that are accurate enough, since the numerical data is very close to the experimental results even for very large time steps.


Once all the particles have been convected, information is projected into the mesh. The interface is drawn by determining the exact place where the material properties change. This is defined in the mesh by an auxiliary pseudo level-set scalar function <math display="inline"> \varphi </math> that is zero at these lines(2D) or surfaces (3D). A correct location of the interface is critical, since the finite element mesh is enriched at the interface as it will be explained later, which allows for a more accurate solution. Figure [[#img-2|2]] shows the position of the interface for a given distribution of particles. The projection kernel to transfer the information from the particles to the mesh is a critical part of the strategy. Since no information is stored in the mesh, this kernel must be accurate enough as to guarantee that the projected field of the variables(in the mesh) resembles the field of the original domain (the particles). Failing to do so would mean that a different problem is being solved at each domain and therefore degrading accuracy. <div id='img-2'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Samper_820077714-particles_interface.png|250px|Interface detection using the particles. φ= 0  at the red line. ]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 2:''' Interface detection using the particles. <math>\varphi = 0</math>  at the red line. 
|}

Having detected the interface, the system of equation is assembled and solved implicitly. Since the sound speed is infinite in incompressible fluids, the complete set of equations cannot be calculated explicitly.

Finally, once the full system of equations  has been solved, corrections are passed to the particles and the cycle is restarted. Details of the convection stage and the procedure for solving the system of equations is explained in the following sections.

===2.2 Continuum equations===

As explained in the last section, the use of particles allows us to omit the convective term in the momentum equations. The system of equations to be solved in the physical domain <math display="inline">\Omega </math> is obtained by coupling the momentum balance for a general material with a compressibility/incompressibility equation. Together, these two equations provide exact solutions for problems in which thermal terms are considered to be uncoupled from the mechanical equations

<span id="eq-2.a"></span>
<span id="eq-2.b"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\displaystyle \rho \frac{D  ( {\mathbf{V}} )}{D t}  = \nabla \sigma   + \rho \mathbf{g} \qquad \hbox{in}  \Omega   </math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (2.a)
|-
| style="text-align: center;" | <math> \displaystyle \frac{D  p}{D t} + \kappa \nabla .   {\mathbf{V}} = 0  \qquad  \quad   \hbox{in}  \Omega    </math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (2.b)
|}
|}

Equations [[#eq-2.a|2.a]] are complemented by the standard boundary conditions of prescribed velocities and prescribed tractions at the Dirichlet (<math display="inline">\Gamma _u</math>) boundary conditions and Neumann (<math display="inline">\Gamma _t</math>) boundary conditions, respectively.

Eqquations [[#2.2 Continuum equations|2.2]] hold for any material, both fluids and solids, where <math display="inline">\rho </math> is the density, <math display="inline">\kappa </math> the compressibility modulus, <math display="inline">\mathbf{V}</math> the velocity vector, <math display="inline">\sigma </math> the stress tensor and <math display="inline">\mathbf{g}</math> the mass force vector. In this work fluids are treated as incompressible, meaning that <math display="inline">\kappa = \infty </math> and, therefore, the second equation simplifies for fluids to <math display="inline">\nabla .   {\mathbf{V}} = 0</math>.

====2.2.1 Continuum equations for fluids====

For the particular case of fluids, the stress tensor depends only on the current strain rate and the pressure as <math display="inline">\sigma _f = 2 \mu   \left(\nabla {\mathbf{V}}^S  \right) - \nabla p </math>. Since the aim of this work is to deal with problems of relatively low velocities <math display="inline">  \| \mathbf{V} \| <0.3 c </math> (where <math display="inline">c</math> is speed of sound), the fluid will be considered to be incompressible. Replacing into [[#2.2 Continuum equations|2.2]], the equations for fluids yields

<span id="eq-3"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\rho  \frac{D  ( {\mathbf{V}} )}{D t}    = \nabla \left[2 \mu   \left(\frac{\nabla {\mathbf{V}}+ (\nabla {\mathbf{V}})^T}{2}  \right) \right] - \nabla p + \rho \mathbf{g} </math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (3)
|-
| style="text-align: center;" | <math>  \nabla .  {\mathbf{V}}  = 0  </math>
|}
|}

It must be noted that the numerical formulation of this set of equations must be stabilized due to the infinite sound speed caused by the incompressibility constraint.

====2.2.2 Continuum equations for solids====

The hypoelastic model provides an excellent constitutive model for solids due to its simplicity and direct application for velocity formulations. Hypoelastic materials are those whose stress rate can be defined simply using the rate of deformation tensor <math display="inline">\mathbf{d}</math> as:

<span id="eq-4"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\dot \sigma  = f (\sigma  , \mathbf{d}) \qquad \hbox{where} \quad d_{ij} = \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;white-space: nowrap;" | (4)
|}

where the upper dot <math display="inline">\dot{( \cdot )})</math> represents the time derivative.

The constitutive relationship for hypoelastic materials defines the rate of stress rather than the stress itself. This is particularly useful in the velocity formulation since the rate stress function <math display="inline">f</math> depends on the rate of deformation <math display="inline">\mathbf{d}</math> (depending on the velocity) instead of the actual deformation. However, the Cauchy stress tensor <math display="inline">\sigma </math> is not objective (it is affected by rigid body rotations). Therefore, another stress measure is required. Using the Truesdell stress rate <math display="inline">\sigma  ^ {\nabla }</math> and eliminating the stress dependency of the model for simplicity, the stress rate relationship becomes:

<span id="eq-5"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\sigma  ^ {\nabla }  = f (\mathbf{d})  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (5)
|}

The Truesdell rate is related to the Second Piola-Kirchhof (PK2) Stress tensor <math display="inline">\mathbf{S}</math> and provides the exact stress rate for a given deformation as

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

Using the tensor of elastic moduli for the Truesdell stress rate <math display="inline">\mathbf{C}^{\tau }</math> yields

<span id="eq-7"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\sigma  ^ {\nabla } = \mathbf{C} ^ {\tau } : \mathbf{d}   </math>
|-
| style="text-align: center;" | <math>   \sigma  ^ {\nabla } = ( \lambda \mathbf{I} \otimes \mathbf{I} + 2 \mu   \mathbf{I}) :  \mathbf{d} -  \mathbf{d}  \cdot \sigma  -  \sigma   \cdot \mathbf{d} + \sigma  \otimes \mathbf{I}   </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (7)
|}

where <math display="inline">\lambda </math> and <math display="inline">\mu </math> are the Lame parameters. As it can be seen in equation [[#eq-7|7]] the constitutive expression is not linear since <math display="inline">C^\tau </math> depends on <math display="inline">\sigma </math>. To overcome this inconvenience, the Jaumann stress rate <math display="inline">\sigma  ^ {\nabla J}</math> <span id='citeF-5'></span>[[#cite-5|[5]]] is used, which simplifies the stress rate into a deformation  plus a rotation, avoiding all the non-linear terms. After some algebra, this yields

<span id="eq-8"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\sigma  ^ {\nabla J} = ( \lambda \mathbf{I} \otimes \mathbf{I} + 2 \mu   \mathbf{I}) :  \mathbf{d}   </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (8)
|}

Approximating the Truessdel stress rate by the Jaumann stress rate gives, combining equations [[#eq-8|8]] and  [[#eq-6|6]]:

<span id="eq-9"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\mathbf{F}^{(n+1)} \dot{\mathbf{S}}^{(n+1)} \mathbf{F}^{T(n+1)} = J^{(n+1)}  ( \lambda \mathbf{I} \otimes \mathbf{I} + 2 \mu   \mathbf{I}) :  \mathbf{d}^{n+1}  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (9)
|}

We use finite differences in time to obtain <math display="inline">\dot{\mathbf{S}}</math>

<span id="eq-10.a"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\mathbf{F}^{(n+1)} \frac{\mathbf{S}^{(n+1)} - \mathbf{S}^{(n)} }{\Delta t} \mathbf{F}^{T(n+1)} = </math>
|-
| style="text-align: center;" | <math> J^{(n+1)}  ( \lambda \mathbf{I} \otimes \mathbf{I} + 2 \mu   \mathbf{I}) :  \mathbf{d}^{n+1}  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (10.a)
|}

Hence,

<span id="eq-10.b"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\mathbf{F}^{(n+1)} \frac{\mathbf{S}^{(n+1)}}{\Delta t} \mathbf{F}^{T(n+1)} =  J^{(n+1)}  ( \lambda \mathbf{I} \otimes \mathbf{I} + 2 \mu   \mathbf{I}) :  \mathbf{d}^{n+1} + </math>
|-
| style="text-align: center;" | <math>   \mathbf{F}^{(n+1)} \frac{\mathbf{S}^{(n)} }{\Delta t} \mathbf{F}^{T(n+1)}  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (10.b)
|}

Recalling the relation between the PK2 stress and the Cauchy stress (<math display="inline"> J \sigma  = \mathbf{F S F}^T</math>) Equation [[#eq-10.b|10.b]] can be expressed as

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>J^{(n+1)} \frac{\sigma ^{(n+1)}}{\Delta t} = J^{(n+1)}  ( \lambda \mathbf{I} \otimes \mathbf{I} + 2 \mu   \mathbf{I}) :  \mathbf{d}^{n+1} + </math>
|-
| style="text-align: center;" | <math>   \mathbf{F}^{(n+1)} \frac{\mathbf{S}^{(n)} }{\Delta t} \mathbf{F}^{T(n+1)} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (11)
|}

Resetting the reference framework at each time step as the last known configuration, i.e. <math display="inline">\mathbf{x}^{(n+1)}:=\mathbf{X}^{(n)}</math>, the Cauchy stress for the previous time step directly becomes the PK2 stress , i.e. <math display="inline">\mathbf{S}^{(n)}:=\sigma ^{(n)}</math>. Finally the Cauchy stresses in the new configuration are obtained by

<span id="eq-12"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\sigma ^{(n+1)} =  J^{-1 (n+1)}  \mathbf{F}^{(n+1)}  \sigma ^{(n)} \mathbf{F}^{T(n+1)}  + </math>
|-
| style="text-align: center;" | <math>\Delta t  ( \lambda \mathbf{I} \otimes \mathbf{I} + 2 \mu   \mathbf{I}) :  \mathbf{d}^{n+1} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (12)
|}

or

<span id="eq-13"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\sigma ^{(n+1)} =   \sigma ^{(n)}_{n+1}  + \Delta t  ( \lambda  tr(\mathbf{d}^{n+1}) \mathbf{I} + 2 \mu   \mathbf{d}^{n+1} ) </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (13)
|}

where the subindex <math display="inline">{n+1}</math> implies that the stresses at time <math display="inline">t^n</math> have been updated to the new configuration.

To mimic the formulation of velocities and pressure used in fluids, it is useful to decompose the stresses into the pressure <math display="inline">p= - \frac{1}{3} tr(\sigma )</math> caused by the volumetric deformation <math display="inline">\epsilon _v=\frac{1}{3} tr(\mathbf{d})</math> and the deviatoric part <math display="inline">\sigma '</math> caused by the deviatoric strain <math display="inline">\mathbf{d}'=\mathbf{d} - \epsilon _v \mathbf{I}</math>. Introducing this splitting in [[#eq-13|13]] yields

<span id="eq-14"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\sigma ^{(n+1)} =   \sigma ^{(n)}_{n+1}  + \Delta t \left[(\lambda +\frac{2}{3} \mu )  tr(\mathbf{d}^{n+1}) \mathbf{I} + 2 \mu  ({\mathbf{d}'}^{n+1}) \right] </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (14)
|}

Finally,

<span id="eq-15"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\sigma ^{(n+1)} =   \sigma ^{(n)}_{n+1}  + \Delta p \mathbf{I} + \Delta \sigma'  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (15)
|}

where

<span id="eq-16"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\Delta p = \Delta t   (\lambda +\frac{2}{3} \mu )  tr(\mathbf{d}^{n+1}) </math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (16)
|-
| style="text-align: center;" | <math>  \Delta \sigma ' = 2 \Delta t  \mu   ({\mathbf{d}'}^{n+1})  </math>
|}
|}

Equations [[#eq-15|15]] and [[#eq-16|16]] provide the expressions for updating the stresses at each time step. A more detailed explanation can be found in <span id='citeF-3'></span>[[#cite-3|[3]]] and <span id='citeF-6'></span>[[#cite-6|[6]]].

==3 Fixed mesh domain==

===3.1 Spatial discretization===

The chosen spatial discretization procedure is the Finite Element Method (FEM) <span id='citeF-7'></span>[[#cite-7|[7]]]. It consists on dividing the domain into elements, whose geometry is defined by nodes. The unknown variables are the values at the nodes and the solution is interpolated from the nodal values inside each element. In this work linear shape functions will be used for all the variables. Three-noded triangles (2D) and 4-noded tetrahedra will be used for all the examples shown in this work. Using the FEM discretization and the test functions <math display="inline">\mathbf{w}</math> for the velocity and <math display="inline">q</math> for the pressure, the weak form <span id='citeF-7'></span>[[#cite-7|[7]]] of the system  [[#2.2 Continuum equations|2.2]] becomes

<span id="eq-17.a"></span>
<span id="eq-17.b"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\left({\mathbf{w}},  \rho \frac{D  ( {\mathbf{V}} )}{D t}  \right)_ \Omega   = \left({\mathbf{w}}, \rho \mathbf{g}  \right)_ \Omega +  \left({\mathbf{w}}, \mathbf{t}^p  \right)_ {\Gamma _t} - \left(\nabla \mathbf{w} , \sigma  \right)_ \Omega   </math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (17.a)
|-
| style="text-align: center;" | <math>  \left(q , \nabla . ( {\mathbf{V}}^{n+1} ) \right)_ \Omega = 0  </math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (17.b)
|}
|}

Note that the stress term is integrated by parts to avoid second space derivatives on the shape functions.

In [[#eq-17.a|17.a]], <math display="inline">\mathbf{t}^p</math> are the prescribed tractions on the Neumann boundaries <math display="inline">\Gamma _t</math>. On the other hand, we have omitted the boundary terms <math display="inline"> \left(\nabla \mathbf{w} , \mathbf{n}  \sigma  \right)_ \Gamma   </math> that appear due to the integration by parts. These terms are the boundary traction terms at the edges <math display="inline">\Gamma </math> of all the elements in the mesh, with the normal vectors <math display="inline">\mathbf{n}</math> pointing outwards. Omitting these terms ensures that traction stress continuity is fulfilled on the inter-element boundaries, where <math display="inline">\Gamma _1</math> and <math display="inline">\Gamma _2</math> are the boundaries of two adjacent elements, i.e.

<span id="eq-18"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\left( \mathbf{w} , \mathbf{n}_1  \sigma _1 \right)_ {\Gamma _1} - \left( \mathbf{w} , \mathbf{n}_1  \sigma _2 \right)_ {\Gamma _2} = 0 </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (18)
|}

Ensuring that [[#eq-18|18]] is satisfied is of particular interest in this work since the pressure field has to be discontinous between interfaces due to the different  material properties. Fullfilling this requirement by omitting this boundary term provides a straightforward yet accurate approach, even if discontinuities occur in the pressure field.

Using a Galerkin method (with <math display="inline">\mathbf{w}=\mathbf{N}</math>, where <math display="inline">\mathbf{N}</math> is the shape functions matrix) the set of algebraic equations from [[#3.1 Spatial discretization|3.1]] can be written as

<span id="eq-19"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\begin{bmatrix}{\mathbf{K}}(\mu )+ \frac{{\mathbf{M}}(\rho )}{\Delta t} & {\mathbf{G}}\\ {\mathbf{G}}^T & 0 \end{bmatrix} \begin{bmatrix}{\mathbf{\bar V}}^{n+1} \\ \bar p^{n+1} \end{bmatrix} = \begin{bmatrix}\frac{1}{\Delta t}{\mathbf{M  \bar V}}^n + \mathbf{Mg} \\ 0 \end{bmatrix} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (22)
|}

In Equations [[#eq-19|19]], first order finite differences are used in the acceleration term to advance in time and <math display="inline"> {\bar{(\cdot )}} </math> denotes nodal variables. The sub-matrices of the system are the standard ones obtained after using the Galerkin FEM for fluid problems <span id='citeF-8'></span>[[#cite-8|[8]]].

===3.2 Stabilized matrix form for the fluid phase===

A FEM discretization directly implemented in the form [[#eq-19|19]] using the same shape functions for the velocity and the pressure would be unstable in the pressure since it does not satisfy the inf-sup condition <span id='citeF-9'></span>[[#cite-9|[9]]]. Moreover, even if different functions were to be used to avoid this restriction, the system of equations would not be suited for several  linear solvers since it contains zeros in the diagonal terms of the pressure equation as seen in the matrix form [[#eq-19|19]].

In order to overcome this limitation, stabilization terms can be added to the set of equations [[#eq-19|19]]. As stated before, it is necessary to stabilize the pressure equation for the fluid problem only. While there are several methods to do so, only two will be mentioned in this work: the Pressure-Stabilizing/Petrov-Galerkin (PSPG) method <span id='citeF-10'></span>[[#cite-10|[10]]] and the Orthogonal Subgrid Scale (OSS) stabilization  <span id='citeF-11'></span>[[#cite-11|[11]]]. The reason behind this is that both methods add a Laplacian to the LHS of the system of equations [[#eq-19|19]], thus preserving the symmetry of the system.

In Equations [[#eq-23|23]] the final matrix system is presented for the PSPG method. Using the OSS method with an explicit projection stage would lead to a similar expression but with the explicit projection of the pressure in the RHS.

<span id="eq-23"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\begin{bmatrix}{\mathbf{K}}(\mu )+ \frac{1}{\Delta t} {\mathbf{M}}(\rho ) & {\mathbf{G}}\\ {\mathbf{G}}^T  & \tau L(\frac{1}{\rho }) \end{bmatrix} \begin{bmatrix}{\mathbf{\bar V}}^{n+1} \\ \bar p^{n+1} \end{bmatrix} = \begin{bmatrix}\frac{1}{\Delta t}{\mathbf{M \bar V}}^n + \mathbf{Mg} \\ \tau {\mathbf{G}}^T \mathbf{g} \end{bmatrix} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (26)
|}

In Equations [[#eq-23|23]], <math display="inline">\mathbf{L}(\frac{1}{\rho })</math> is the Laplacian of the pressure. Another alternative to stabilize the equations is the Finite Increment Calculus (FIC) procedure <span id='citeF-6'></span>[[#cite-6|[6]]].The FIC method avoids the need to prescribe the pressure in free surface solvers using staggered schemes. A mixed FEM Lagrangian formulation for treating both quasi and fully incompressible fluids as well as FSI problems using the FIC method is presented in <span id='citeF-6'></span>[[#cite-6|[6]]].

The system [[#eq-23|23]] can be solved by choosing an appropriate stabilization parameter <math display="inline">\tau </math>. In this work <math display="inline">\tau ={\left(1/ \Delta t + 4 \mu / h^2 + 2 V^n / h \right)} ^{-1} </math> following the work in <span id='citeF-12'></span>[[#cite-12|[12]]]<span id='citeF-10'></span>[[#cite-10|[10]]]. The dependency on <math display="inline">\Delta t</math> was added to avoid <math display="inline">\tau  \rightarrow \infty </math> when the velocity and the viscosity are zero <span id='citeF-10'></span>[[#cite-10|[10]]].

===3.3 Matrix form for the solid phase===

Unlike fluids, solids will not be considered incompressible. The spatial discretization and variables will be the same as for the fluids formulation: linear elements for both the pressure and the velocity. As for the deviatoric stress, it will not be an independent variable, but a state variable depending on the velocity unknowns.

Recalling the general system of equations [[#2.2 Continuum equations|2.2]] and splitting the stresses into its  deviatoric and volumetric components, the system of governing equations reads:

<span id="eq-27"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\rho  \frac{D  ( {\mathbf{V}} )}{D t}  = \nabla \left(\sigma ' - p  \mathbf{I}  \right)+ \rho \mathbf{g} \qquad \hbox{in}  \Omega  </math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (27)
|-
| style="text-align: center;" | <math> \frac{D  p}{D t} + \kappa \nabla .   {\mathbf{V}} = 0 \qquad \hbox{in}  \Omega  </math>
|}
|}

The weak form of the problem is

<span id="eq-28"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\left(\mathbf{w} , \rho  \frac{D  ( {\mathbf{V}} )}{D t} \right)_ \Omega = - \left(\nabla \mathbf{w} ,  \sigma ' - p  \mathbf{I}   \right)_ \Omega  </math>
|-
| style="text-align: right;" | <math>  + \left(\mathbf{w} , \rho \mathbf{g} \right)_ \Omega + \left(\mathbf{w} , \mathbf{t}^p \right)_ {\Gamma _t} </math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (28)
|-
| style="text-align: center;" | <math> \left(q , \frac{D  p}{D t} + \kappa \nabla .   {\mathbf{V}} \right)_ \Omega = 0  </math>
|}
|}

A first approximation in time for the velocity would lead to a stable but too diffusive algorithm. To avoid this shortcoming, the Newmark's Beta method was selected. Parameters were set for constant acceleration. The resulting expression for the velocity integration becomes:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\qquad \qquad 2 \rho  \frac{ \mathbf{V}^{n+1}}{\Delta t} = 2 \rho  \frac{ \mathbf{V}^{n}}{\Delta t} -  \rho  \mathbf{ \dot V}^{n}  </math>
|}
|}

Using finite differences in time for the pressure, dividing the second equation by the compressibility modulus <math display="inline">\kappa </math> and recalling the definition of the stresses [[#eq-15|15]] yields:

<span id="eq-29"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\left(\mathbf{w} , \frac{\rho \mathbf{V}^{n+1} }{\Delta t} \right)_{\Omega } - \left(\nabla \mathbf{w},    \Delta \sigma ' - \Delta p  \mathbf{I}  \right)_{\Omega }   = \left(\mathbf{w} ,\frac{\rho \mathbf{V}^{n} }{\Delta t} + \rho \mathbf{g} \right)_{\Omega }  </math>
|-
| style="text-align: center;" | <math>  + \left(\mathbf{w} , \mathbf{t}^p \right)_ {\Gamma _t} - \left(\nabla \mathbf{w},  \Delta {\sigma '}^n_{n+1} - (p^n)  \mathbf{I}   \right)_{\Omega }  </math>
|-
| style="text-align: center;" | <math> \qquad \qquad \qquad \quad \left(q , \frac{p^{n+1}}{\kappa \Delta t} + \nabla .   {\mathbf{V}} \right)_{\Omega } = \left(q , \frac{p^{n}}{\kappa \Delta t} \right)_{\Omega }   </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (29)
|}

The rate of defomations is calculated implicitly to avoid instabilities as

<span id="eq-30"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\qquad \qquad \mathbf{d} = \mathbf{B}   \mathbf{V}^{n+1} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (30)
|}

where <math display="inline">\mathbf{B}</math> is the standard strain-rate velocity matrix [[#eq-33|33]]. Finally, writing the equations in matrix form, the system can be expressed as:

<span id="eq-31"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\begin{bmatrix}\left(\frac{2}{\Delta t} \mathbf{M}(\rho ) + \Delta t \mathbf{K} \right)   & \mathbf{G} \\ \mathbf{G}^T & \frac{1}{\Delta t} \mathbf{M}(\frac{1}{\kappa }) \end{bmatrix} \begin{bmatrix}\mathbf{\bar V}^{n+1} \\ \bar p^{n+1} \end{bmatrix} = </math>
|-
| style="text-align: right;" | <math> \begin{bmatrix}\mathbf{M}(\rho )  ( \frac{2}{\Delta t} \mathbf{\bar V}^{n} - \mathbf{ \bar {\dot V}}^{n} + \mathbf{g} )+  \mathbf{G} {\sigma '}^n_{n+1} \\ \frac{1}{\Delta t} \mathbf{M}(\frac{1}{\kappa }) \bar p^n \end{bmatrix} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (31)
|}

In general all the matrices have the same form used for the fluid elements, except for the stiffness matrix <math display="inline">\mathbf{K}</math>. Its expression is:

<span id="eq-32"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\qquad \qquad \mathbf{K} =  \int _{\Omega } \mathbf{B^T C'  B} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (32)
|}

where <math display="inline">\mathbf{C'}</math> is the constitutive matrix for the deviatoric stresses. The expressions for <math display="inline">\mathbf{B}</math> and <math display="inline">\mathbf{C'}</math> for 3D problems are

<span id="eq-33"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\mathbf{B} = \begin{bmatrix}\frac{\partial N}{\partial x}    & 0 & 0 \\ 0    & \frac{\partial N}{\partial y} & 0 \\ 0    & 0 & \frac{\partial N}{\partial z} \\ \frac{\partial N}{\partial y} & \frac{\partial N}{\partial x} & 0 \\ 0 & \frac{\partial N}{\partial z} & \frac{\partial N}{\partial y} \\ \frac{\partial N}{\partial z} & 0 &  \frac{\partial N}{\partial x} \end{bmatrix}   \mathbf{C'} = \mu  \begin{bmatrix}4/3    & -2/3 & -2/3 & 0 & 0 & 0 \\ -2/3 & 4/3    & -2/3 & 0 & 0 & 0 \\ -2/3 & -2/3 &    4/3 & 0 & 0 & 0 \\   0  &   0  &   0  & 1 & 0 & 0 \\   0  &   0  &   0  & 0 & 1 & 0 \\   0  &   0  &   0  & 0 & 0 & 1 \end{bmatrix} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (43)
|}

It must be noted that even for 2D simulations, the terms in <math display="inline">\mathbf{B}</math> and <math display="inline">\mathbf{C'}</math> do not change. The only modification is the elimination of the strains in the third dimension.

The system of equations [[#eq-31|31]] can be solved without the use of stabilized formulations, as long as the solid material is not incompressible. Otherwise, the pressure mass matrix becomes zero and therefore stabilization is required. In this case the same stabilization technique used for the fluid elements can be used.

An additional hypothesis has to be made in order to linearise the system. The transformation of the stresses in the previous configuration to the new configuration requires the updated velocity <math display="inline">\mathbf{V}^{n+1}</math>, making the system non-linear. The simplification used here consists on using the previous step velocity, therefore assuming small accelerations for the solids. This hypothesis is similar to the one used for the convective term, which showed good accuracy in all the cases tested.

===3.4 Special considerations on the interface===

Despite the advantages that the FEM discretization provides, it bears a severe limitation when the variables undergo abrupt changes that the chosen FEM space  is unable to reproduce. As an example,when there is a sharp change of the density in a  two-fluid problem, the hydrostatic condition under the gravitational force leads to two different values for the pressure gradient. This is specially critical when the fluids have very different densities (i.e. air-water), where the term <math display="inline"> \nabla p =  \rho g </math>  changes abruptly at the interface.

<div id='img-3'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Samper_820077714-HYDROSTATIC.png|180px|Pressure distribution for the hydrostatic case]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 3:''' Pressure distribution for the hydrostatic case
|}

Figure [[#img-3|3]] shows the exact pressure distribution for the hydrostatic case of water and air and the one obtained using three linear elements. It is clear that this solution is very poor and in practice it would lead to an incorrect behavior of the interface and mass losses <span id='citeF-13'></span>[[#cite-13|[13]]].

When dealing with FSI problems, this issue is even more critical. The reason for this is again the strong discontinuity in the material properties, in this case the stiffness (or the viscosity). A simple example to illustrate this is a box of a  solid material surrounded by air under gravitational force. In this case the pressure distribution would take the shape illustrated in Figure&nbsp;[[#img-4|4]].

<div id='img-4'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Samper_820077714-box_in_air.png|138px|Pressure distribution of a solid body surrounded by air under gravitational force]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 4:''' Pressure distribution of a solid body surrounded by air under gravitational force
|}

The fulfilment of the transmission conditions requires that the tractions along the interface must be continuous. Neglecting the air pressure, this condition is satisfied by <math display="inline">\sigma _{xx} = 0</math> at the left and right boundaries of the solid. However, the implemented solid model makes use of the mean pressure in the solid, which means that (taking <math display="inline">\sigma _{zz}=0</math>) the pressure in the solid becomes

<span id="eq-44"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\qquad \qquad \qquad p_{solid} = \frac{1}{3} \sigma _{yy} \ne p_{air} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (44)
|}

The pressure given by [[#eq-44|44]] is clearly discontinuous at the interface. If a linear (or any continuous) element was located over the sharp boundary, the discrete pressure field would detect a non-existent gradient that would create normal velocities in the fluid, as seen in Figure [[#img-5|5]]. Since this would violate fluid incompressibility, the solver usually converges to a solution where the last solid node has a very low pressure value and, therefore, the structure is more flexible than it should be.

<div id='img-5'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Samper_820077714-discretized_box_pressure.png|138px|Close-up of a linear element located at the interface showing the interpolated pressure field]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 5:''' Close-up of a linear element located at the interface showing the interpolated pressure field
|}

===3.5 Raising the limitations: pressure enrichments===

To overcome the problems of standard finite elements, the space must be modified to allow for a more accurate reproduction of the solution. A first alternative would be remeshing the zone crossed by the interface. This would lead to the exact solution but would require adding new degrees of freedom (DoFs) to the new discretized geometry. This strategy, despite being possible, is computationally expensive since the system must be resized and new memory space must be allocated in the memory, task that is likely to be slower than the solution of the system itself.

An alternative to a new mesh is enrichening the FEM space. Enriching consist on creating new DoFs on each side of the interface elements and then statically condensing the new unknowns <span id='citeF-14'></span>[[#cite-14|[14]]].

In <span id='citeF-15'></span>[[#cite-15|[15]]] it is shown that the solution improves when more flexibility is given to the pressure field at the interface. Following this line, three enrichments shape functions plus special shape functions to replace the standard pressure field are used. Figure [[#img-6|6]] shows that a total of 6 functions are used, with the goal of fully uncoupling the pressure from both sides. To do so, first the standard shape functions in the interface elements are replaced by their discontinuous counterparts. These functions are basically the same as the original ones, but the integrals are calculated only with the contribution of the partitions whose sign matches the sign of the node. On the other hand, when partitions and nodes have different signs, the contribution is added to the enrichments functions <math display="inline">i^*,j^*</math> or <math display="inline">k^*</math>. In other words, the original shape functions are split in two independent functions across the interface. This allows more freedom in the pressure field, while retaining the partition of unity <span id='citeF-7'></span>[[#cite-7|[7]]].

<div id='img-6'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Samper_820077714-enrichments_ricc3.png|210px|Pressure shape functions: Top: Replacement functions. Bottom: Enrichments]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 6:''' Pressure shape functions: Top: Replacement functions. Bottom: Enrichments
|}

Stabilizing the pressure equation in the fluid phase is equivalent to adding artificial diffusion to this variable. Therefore, adding stabilization to the interface elements would inevitably lead to an inaccurate behavior of the solid phase due to dissipation of the stresses in the interface region. For this reason no stabilization is used in these elements. However, as stated previously, doing so would cause instabilities due to the non-fulfilment of the inf-sup condition <span id='citeF-11'></span>[[#cite-11|[11]]]. To override this problem, both sides of the interface are considered to be compressible (using the compressibility modulus of the solid phase). While this hypothesis is not correct, it leads to small errors only since the fluid volume is changed only when there is a sharp change in the pressure from one timestep to the next. This strategy yields almost exact results, avoiding pressure diffusion across the interface with minimal errors.

'''Remark:''' One of the advantages of this set of discontinuous functions is that the partition of unity is conserved. But most importantly, this also implies that the reconstruction of the pressure <math display="inline">p_n</math> after convection is easier to perform accurately. The strategy currently implemented consists on the following steps. First, a standard nodal projection algorithm is used and, therefore, the values of the replacement shape functions are set. Having done this, the enrichment shape functions are defined as the mean of the real nodal values for the previous timestep <math display="inline">n</math>. This ensures that the prediction for the gradient in the normal direction is minimized, thus reducing spurious velocities.

====3.5.1 Adding the new DoFs to the system of equations====

The enrichment functions can be added to the monolithic strategy following a standard condensation procedure <span id='citeF-14'></span>[[#cite-14|[14]]]. Using the subscript <math display="inline">N</math> for the standard shape functions and <math display="inline">*</math> for the enrichment functions, the extended system becomes

<span id="eq-45"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\left[ \begin{array}{cc|cc}{\mathbf{K}}_{NN}+ {\mathbf{M}_N}  & {\mathbf{G}}_{NN} & \mathbf{G}_{N*} \\ {\mathbf{G}_{NN}}^T  & \mathbf{M}_N(\frac{1}{\kappa }) & 0 \\  {\mathbf{G}_{N*}}^T  & 0 & \mathbf{M}_*(\frac{1}{\kappa }) \end{array} \right] \left[ \begin{array}{c}{\mathbf{\bar V_N}}^{n+1}  \\ \bar p_N^{n+1} \\  \bar p_*^{n+1} \end{array} \right] </math>
|-
| style="text-align: center;" | <math> = \left[ \begin{array}{c}\frac{1}{\Delta t} {\mathbf{M_N \bar V_N }}^n + \mathbf{M_N g} \\ {\mathbf{M}_N(\frac{1}{\kappa }) {\bar p}^n_N } \\  {\mathbf{M}_*(\frac{1}{\kappa }) {\bar p}^n_* } \end{array} \right] </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (51)
|}

Note that the density <math display="inline">\rho </math> is omitted in the velocity mass matrices due to lack of space. Also, since all mass matrices (both velocity and pressure) are multiplied by <math display="inline">\frac{1}{\Delta t}</math>, this is also omitted for the same reason. To simplify the notation, the system is written as:

<span id="eq-52"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\qquad  \left[ \begin{array}{c|c}\mathbf{A_{NN}} & \mathbf{A_{N*}} \\  \mathbf{A_{*N}} & \mathbf{A_{**}} \end{array} \right] \left[ \begin{array}{c}\mathbf{X_{N}} \\  \mathbf{X_{*}} \end{array} \right] = \left[ \begin{array}{c}\mathbf{F_{N}} \\  \mathbf{F_{*}} \end{array} \right] </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (55)
|}

Condensing the enrichment DoFs the initial size of the system is recovered, as:

<span id="eq-56"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\qquad \qquad \qquad \mathbf{\tilde{A}} \mathbf{{X}_N}= \mathbf{\tilde{F}} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (56)
|}

where

<span id="eq-57"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\qquad \qquad \qquad \mathbf{\tilde{A}} = \mathbf{{A}_{NN}} - \mathbf{{A}_{N*}} \mathbf{{A}^{-1}_{**}} \mathbf{{A}_{*N}} </math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (57)
|-
| style="text-align: center;" | <math>  \qquad \qquad \qquad \mathbf{\tilde{F}} = \mathbf{{F}_{N}} - \mathbf{{A}_{N*}} \mathbf{{A}^{-1}_{**}} \mathbf{{F}_{*}} \mathbf{\tilde{F}} </math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (58)
|}
|}

==4 Lagrangian particles domain==

===4.1 Streamline integration===

The solution strategy described in the previous section requires an accurate and robust tool to account for the convective term. Since the aim of the PFEM-2 is to achieve large time steps with small errors, it is not enough to simply move the particles with the last velocity multiplied by the time step (<math display="inline">\delta \mathbf{x} =  {V}^n \cdot \Delta _t</math>).

The main goal is to obtain as much information as possible from the previous time step to calculate the new position of the particles. The streamline integration method used in this work follows the same strategy of <span id='citeF-1'></span>[[#cite-1|[1]]] for single fluids. For a particle <math display="inline">p</math>, its exact position at time <math display="inline">n+1</math> is computed as [[#eq-59|59]]:

<span id="eq-59"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\mathbf{x}^{n+1}_p = \mathbf{x}^{n}_p + \int _{n}^{n+1} \mathbf{V^t}   dt </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (59)
|}

Unfortunately, the velocity is not known continuously in time but only at discrete time steps( <math display="inline">n-1</math>,<math display="inline">n</math>, <math display="inline">n+1</math> ... ). Approximating the particle trajectory using the velocity of the previous time step to obtain an explicit scheme, the new position of the particle becomes:

<span id="eq-60"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\mathbf{x}^{n+1}_p \approx \mathbf{x}^{n}_p + \int _{n}^{n+1} \mathbf{V^n}   dt </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (60)
|}

While this approximation may be valid for quasi-stationary problems, it can lead to severe errors in transient situations. In <span id='citeF-1'></span>[[#cite-1|[1]]] it was shown that this strategy provides good results for single fluids. However in multi-fluids and FSI problems the difference in the properties, in particular the density, makes this approximation unsuitable for some cases. When the density jump is important, one fluid dominates over the other and their behaviour is completely different (ie: air and water). As an example, a breaking wave is represented in Figure [[#img-7|7]]. Following the air streamlines to predict the water particle trajectory would lead to large errors.

<div id='img-7'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Samper_820077714-streamline_wave2.png|120px|Breaking wave with two fluids. Blue line: Streamline at tⁿ;   Green line: Actual particle trajectory]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 7:''' Breaking wave with two fluids. Blue line: Streamline at <math>t^n</math>;   Green line: Actual particle trajectory
|}

To overcome this problem, the streamline integration must be deactivated when a particle of the heavier fluid transits across a region of the lighter fluid. In this case, the gravity will be taken as the only acting force, leading to a parabolic (projectile-type) motion, where the last known velocity of the particle is taken as the starting data. Setting the distance function <math display="inline">\varphi \le 0</math> for the denser fluid regions and <math display="inline">\varphi > 0 </math> for the lighter fluid regions, the updated position of the denser fluid particles becomes

<span id="eq-61"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\mathbf{x}^{n+1}_p \approx \Bigg\{{ \begin{matrix}   \mathbf{x}^{n}_p + \int _{n}^{n+1} \mathbf{V^n}   dt \qquad \qquad \hbox{if}   \varphi _{x_p^t} \le 0  \\   \mathbf{x}^{n}_p + \int _{n}^{n+1}  \mathbf{V^n_p} + g\cdot t   dt \quad \hbox{if}   \varphi _{x_p^t} > 0 \end{matrix}  } </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (62)
|}

It is important to note that this method is purely explicit and, therefore, no further corrections are done for computing the convective term. Even though an iterative process could be performed to improve the accuracy, in all the examples solved with this strategy the errors were small, despite the large time steps chosen. This is specially interesting since the PFEM-2 was compared against other algorithms that make use of implicit strategies and yet they have to rely on smaller time steps to obtain similar results. Based on these tests it was decided to keep only the explicit computation of the convective term, since it provides an excellent balance between computation time and accuracy.

===4.2 Projection kernel===

Once the particles have been convected to their new position, <math display="inline">\mathbf{x}_{n+1}^p</math>, it is necessary to project the information into the mesh in order to solve the Lagrangian system of equations. The simplest projection kernel consists on directly using the shape functions of the elements. As an example, the definition of the interface for two materials is explained next.

Each particle has an associated material and the interface should be located where the material properties change. The distribution of particles inside each element can define complicated curves that become impossible to manage with simple shape functions due to the large number of particles. As an example, in Figure [[#img-8|8]] some particles are in the "wrong" side of the interface. For this reason the instantaneous local interface inside each element is simply defined by a line (or a plane in 3D) taking into account a weighted average using the shape functions. An instantaneous , pseudo level-set function <math display="inline">\varphi </math> is created that defines the interface where it is valued zero. To calculate the value of <math display="inline">\varphi </math> at each <math display="inline">j</math>th node, all the <math display="inline">i</math> particles in the neighbour elements to <math display="inline">j</math> are used. The weighting function  is defined as the standard shape function of the element for the <math display="inline">j</math>th node in the position of each <math display="inline">i</math>th  particle. Once the contribution of all the particles has been added, the location of the interface is computed as

<span id="eq-63"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\varphi _ j = \frac{ \sum _{i}^{n} {sign}_i N_i^j}{ \sum _{i}^{n} N_i^j} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (63)
|}

where <math display="inline">sign_i</math> denotes the sign of the particle, positive (<math display="inline">+1</math>) for the lighter fluid and negative (<math display="inline">-1</math>) for the denser fluid.

<div id='img-8'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Samper_820077714-projection_kernel.png|120px|Contribution of i particle to the j node]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 8:''' Contribution of <math>i</math> particle to the <math>j</math> node
|}

The same procedure is  for all the variables, replacing the <math display="inline">sign_i</math> by the desired variable of the particle, either scalar, vectorial or tensorial. Using this kernel, the projected velocity can be written as

<span id="eq-64"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\hat { \hat { \mathbf{V}  } } _ j ^{n+1} = \frac{ \sum _{i}^{n} \mathbf{V}_i N_i^j}{ \sum _{i}^{n} N_i^j} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (64)
|}

The <math display="inline"> \hat{\hat{(\cdot )}} </math> in the notation implies that it is the first approximation of the velocity, taking into account only the convecting term and neglecting all the other contributions. This velocity will be used as the starting point for the finite element fixed mesh solver

It is worth noting that the element shape functions <math display="inline">N_i</math> are not the only possible kernel. As an example, the function proposed in <span id='citeF-16'></span>[[#cite-16|[16]]] for the  Smoothed Particle Hydrodynamics method assigns a higher weight to those particles that are closer to the node. Setting <math display="inline">h</math> as the element size and <math display="inline">u=R/h</math>, where <math display="inline">R</math> is the distance from the particle to the node, the Wendland kernel is written as

<span id="eq-65"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>W_i^j = \frac{7}{4 \pi h^2}(1+2u)(1-u/2)^4; </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (65)
|}

After the information of the Lagrangian domain has been projected and the mesh computations have been performed, the velocity of the particles has to be updated. At this stage it is important not to replace the particle velocity with the updated velocity but only a correction. The reason for this is simple; no matter how good the kernel algorithm may be, it is never exact. Replacing the unknowns in the particles would destroy valuable data, leading to excessive artificial diffusion in the strategy. Based on this fact, the update stage on the particle will only transfer the variation in the velocity (or any other variable) calculated in the mesh, that is:

<span id="eq-66"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\delta \mathbf{V}  _ j ^{n+1} =   \mathbf{V}  _ j ^{n+1}  - \hat { \hat { \mathbf{V}  } } _ j ^{n+1} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (66)
|}

===4.3 Computational issues and code speed===

The code used in this work was developed in the C++  programming environment KRATOS <span id='citeF-17'></span>[[#cite-17|[17]]]. KRATOS is an open-source multi-physics software framework mainly designed for the FEM, which eases the tasks that are common to all FEM codes, such as the assembly of the system of equations or implementing an efficient solver. Ensuring low execution times at the particles convection stage is critical since the number of particles is 10 to 20 times larger than the number of elements in the mesh for a given domain, meaning that millions of particles have to be displaced at each time step for large problems. The optimization and parallelization of this stage lead to a code that spends approximately half of the execution time on particle tasks and the other half on the solution of the implicit system of equations in the mesh.

The current implementation of this work is capable of solving two phase flows accurately and with execution times similar, or even faster, than industry standard codes such as OpenFoam, as shown in <span id='citeF-1'></span>[[#cite-1|[1]]] for one phase flows and in <span id='citeF-2'></span>[[#cite-2|[2]]] for multi-fluids.

==5 Full PFEM-2 algorithm==

Coupling the strategy described in the previous sections, the PFEM-2 solver for multi-fluids and FSI is obtained. The solution steps within a time step can be summarized as follows:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\hbox{'''Step 1)''' Convect the particles using the streamlines}    </math>
|-
| style="text-align: center;" | <math>   \hbox{ '''Step 2)''' Project information on the mesh,}    </math>
|-
| style="text-align: center;" | <math>     \hbox{ with  }q\hbox{ fluid particles and }r\hbox{ solid particles}    </math>
|-
| style="text-align: center;" | <math>  \qquad   Common  variables  \qquad \qquad Solid  Variables   </math>
|-
| style="text-align: center;" | <math>   \qquad \hat { \hat { \mathbf{V}  } } _ j ^{n+1} = \frac{ \sum _{i}^{r+q} \mathbf{V}_i N_i^j}{ \sum _{i}^{r+q} N_i^j} \qquad \qquad p^n _ j = \frac{ \sum _{i}^{r} p_i N_i^j}{ \sum _{i}^{r} N_i^j}     </math>
|-
| style="text-align: center;" | <math>          \qquad   Fluid  variables \qquad \qquad \qquad \nu _{solid} = \frac{ \sum _{i}^{r} nu _i }{ \sum _{i}^{r} 1}             </math>
|-
| style="text-align: center;" | <math>         \qquad \qquad \mu _{fluid} = \frac{ \sum _{i}^{q} \mu _i }{ \sum _{i}^{q} 1} \qquad \qquad {\sigma '}n __{solid} = \frac{ \sum _{i}^{r} {\sigma '}^n__i }{ \sum _{i}^{r} 1}         </math>
|-
| style="text-align: center;" | <math>          \qquad \qquad \rho _{fluid} = \frac{ \sum _{i}^{q} \rho _i }{ \sum _{i}^{q} 1} \qquad \qquad E _{solid} = \frac{ \sum _{i}^{r} E _i }{ \sum _{i}^{r} 1} </math>
|-
| style="text-align: center;" | <math>       \qquad \qquad \qquad \qquad \qquad \qquad \qquad     \rho _{solid} = \frac{ \sum _{i}^{r} \rho _i }{ \sum _{i}^{r} 1}      </math>
|}
|}

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\hbox{'''Step 3)''' Detect interface elements}    </math>
|-
| style="text-align: center;" | <math>   \hbox{ '''Step 4)''' Estimate pressure enrichments}     </math>
|-
| style="text-align: center;" | <math>   \hbox{ '''Step 5)''' Assemble the system with enrichments}    </math>
|-
| style="text-align: center;" | <math>   \hbox{ '''Step 6)''' Solve the system of equations}    </math>
|-
| style="text-align: center;" | <math>   \hbox{ '''Step 7)''' Recover condensed DoFs}    </math>
|-
| style="text-align: center;" | <math>   \hbox{ '''Step 8)''' Update the velocity, pressure and stress}    </math>
|-
| style="text-align: center;" | <math>   \qquad  \qquad \hbox{ at the particles. }    </math>
|}
|}

==6 Numerical Examples==

===6.1 Rayleigh-Taylor instability===

The classical Rayleigh-Taylor instability benchmark for multi-fluids was tested to verify the PFEM-2 solver. Despite there are no analytical results or experimental data in 2D to check the error, it is possible to compare the solution against other numerical results. He et al. <span id='citeF-18'></span>[[#cite-18|[18]]] presented results for a Reynolds number <math display="inline">Re=256</math> using a Lattice Boltzmann scheme. With the same initial conditions, results for <math display="inline">Re=1000</math> can be found in <span id='citeF-19'></span>[[#cite-19|[19]]].  The geometry of the two-fluid problem can be seen Figure [[#img-9|9]]. We have used a mesh of 163000 3-noded triangles. The parameters of the simulation are:

<math display="inline">\rho _ q = 3  kg/m^3 \qquad  \nu _ q = \frac{\sqrt{g}  \cdot  d^{3/2} }{Re}   m^2/s</math>

<math display="inline">\rho _ p = 1  kg/m^3 \qquad  \nu _ p = \frac{\sqrt{g}  \cdot  d^{3/2} }{Re}  m^2/s</math>

<math display="inline">g _ y =  - 10  m/s^2</math>

<math display="inline">\delta t = 0.01  s</math>

Initial conditions:

<math display="inline">V=0  in  \Omega </math>

<math display="inline"> \varphi =  - y - \delta _0 (cos(2 \pi (x)- \pi )+1 )+2.0</math>

<math display="inline">\delta _0 = 0.1</math>

<div id='img-9'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Samper_820077714-rayleigh_geometry2.png|100px|Rayleigh-Taylor instability geometry]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 9:''' Rayleigh-Taylor instability geometry
|}

Figure [[#img-11|11]] shows snapshots of the two-fluid problem at different time steps. For the Reynolds number analysed, the results obtained match the ones found in <span id='citeF-19'></span>[[#cite-19|[19]]]<span id='citeF-18'></span>[[#cite-18|[18]]]. The non-dimensional time <math display="inline">t_{adim}=t \sqrt{g/2}</math> is used. Finally the position of the advancing front of the denser fluid (centre) and the rising of the bubbles (walls) are compared to the results of <span id='citeF-19'></span>[[#cite-19|[19]]]. Again good matching is observed.

<div id='img-10'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Samper_820077714-rayleigh_table.png|240px|Vertical position of the advancing front and bubbles in the Rayleigh-Taylor instability for Re=1000 <span id='citeF-19'></span>[[#cite-19|[19]]]]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 10:''' Vertical position of the advancing front and bubbles in the Rayleigh-Taylor instability for <math>Re=1000</math> <span id='citeF-19'></span>[[#cite-19|[19]]]
|}

<div id='img-11'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|
{|  style="text-align: center; margin: 1em auto;min-width:50%;width:100%;"
|-
| [[File:Draft_Samper_820077714_4972_rayleigh_re256_1,0.png|100px|rayleigh_re256_1,0]]
| style="padding-left:5px;"|[[Image:Draft_Samper_820077714-rayleigh_re256_1,5.png|100px|rayleigh_re256_1,5]]
| style="padding-left:5px;"|[[Image:Draft_Samper_820077714-rayleigh_re256_1,75.png|100px|rayleigh_re256_1,75]]
| style="padding-left:5px;"|[[Image:Draft_Samper_820077714-rayleigh_re256_2,0.png|100px|rayleigh_re256_2,0]]
| style="padding-left:5px;"|[[Image:Draft_Samper_820077714-rayleigh_re256_2,25.png|100px|rayleigh_re256_2,25]]
| style="padding-left:5px;"|[[Image:Draft_Samper_820077714-rayleigh_re256_2,5.png|100px|rayleigh_re256_2,5]]
|}
|-
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 11:''' Snapshots for the Rayleigh-Taylor instability two-fluid problem at <math>t_{adim}=1.0  ,  1.5  ,  1.75  ,  2.0  ,  2.25  ,   2.5</math>. <math>Re=256</math>
|}

===6.2 Static cantilever beam immersed in fluid===

To test the accuracy of the solver and the potential of the enrichment functions to uncouple the pressure, the following numerical example was  designed. A beam under small deformations is surrounded by a very light fluid. Inertial terms are neglected to test the static solution, leading to a Stokes problem in the fluid. The viscosity is set as small as possible to allow for a faster convergence to the statical solution while maintaining stability. The geometry is depicted in Figure [[#img-12|12]]. The material properties are:

<math display="inline">\rho _q = 1  kg/m^3 \qquad \rho _p = 0.001  kg/m^3</math>

<math display="inline">E = 10^7 Pa \qquad    \mu _p = 0.01 Pa.s </math>

<math display="inline">\nu = 0 </math>

<div id='img-12'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Samper_820077714-fluid_cantilever_geometry.png|240px|Cantilever beam with fluid geometry]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 12:''' Cantilever beam with fluid geometry
|}

The beam dimensions are <math display="inline">L=10</math> and <math display="inline">H=1</math>. A vertical mass force of <math display="inline">q = 1[m/s^2]</math> is applied to the whole domain. According to Timoshenko beam theory and correcting the force to account for the buoyancy effect, the vertical displacement <math display="inline">w</math> at the tip should be

<span id="eq-67"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>w_{max} =  (\rho _q - \rho _p) \left(\frac{q  L^4}{8  E  I} + \frac{1}{2} \frac{q L^2}{10/12   h  \mu } \right) = 0.001511 </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (67)
|}

Structured meshes with decreasing element size were created to check convergence. Mesh sizes from <math display="inline">h=1</math> to <math display="inline">h=0.01</math> were used. To test the influence of the projection algorithm, the simulation was run with and without projection. In the latter case, the nodal pressure is kept as a nodal variable from one time-step to the next instead of being obtained from the particles using the projection kernel.

In order to test the behaviour of the splitted elements, the interface between the solid phase and the fluid phase is located  in all the meshes at approximately the middle of the elements, as shown in Figure [[#img-13|13]]. In this case the enrichment shape functions in the pressure are imperative due to the strong discontinuity of the material properties, leading to a jump in the pressure field. If this variable was not enriched to allow for discontinuities, a strong pressure gradient would appear at the interface and, therefore, the resulting velocity field would be inaccurate. It must be noted that Figure [[#img-13|13]] does not correctly represent the pressure at the interface elements, since the field is printed using linear interpolation, while in fact it is discontinuous due to the enrichments.

<div id='img-13'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Samper_820077714-solid_interface.png|228px|Beam interface and pressure distribution]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 13:''' Beam interface and pressure distribution
|}

Figure  [[#img-14|14]] clearly shows the detrimental effect of the projection stage. While the algorithm without projection is able to maintain the quadratic convergence of the solid-only case, the results with projection show a convergence slope of only <math display="inline">3/2</math>.

<div id='img-14'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Samper_820077714-errors_fluid_cantilever2.png|240px|Convergence for the cantilever beam in fluid]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 14:''' Convergence for the cantilever beam in fluid
|}

Despite the solver with the projection stage provides a less accurate solution, it leads to acceptable errors for all the mesh sizes chosen. Even for coarse meshes of only three elements along the thickness, the errors are still low (at around 10 % ).

===6.3 Turek FSI benchmark===

To test the PFEM-2 algorithm in a fully coupled scenario, the Turek FSI benchmark <span id='citeF-20'></span>[[#cite-20|[20]]], was simulated. This problem is typically used as a reference to validate new FSI strategies and it is particularly challenging for staggered schemes due to the similar properties of the fluid and the solid. On the other hand, the monolithic solver used in this work is well suited for this benchmark. Unlike the previous examples presented in this work (in which the PSPG stabilization was used) in this example the OSS stabilization is employed, with the explicit projection of the pressure gradient in the RHS. This is done to avoid having zeros in the RHS of the pressure equation of Equation [[#eq-23|23]] due to lack of gravity, which would provide excessive numerical diffusion.

The problem models a 2D rigid cylinder attached to an elastic beam immersed in a fluid stream. The dimensions and problem parameters can be found in <span id='citeF-20'></span>[[#cite-20|[20]]]. In this paper a fixed mesh of 48500 3-noded triangles was used. The frequency and amplitude of the oscillations are analyzed to test the accuracy. For the case studied we use <math display="inline">Re=200</math>, with a parabolic inlet velocity profile and the following parameters:

<math display="inline">\rho _{fluid} = 1000  kg/m^3 \qquad  \mu _{fluid}  = 1 Pa.s</math>

<math display="inline">\rho _{solid} = 1000  kg/m^3 \qquad E_{solid} = 1.4 MPa </math>

<math display="inline">\nu _{solid} = 0.4 \qquad \qquad \qquad \delta t = 0.005  s</math>

In Figure [[#img-15|15]] the velocity contour is plotted after the stationary regime has been achieved. Finally, in Figure [[#img-16|16]] the vertical displacement of the tip can be observed. In Table 1 the results are analysed. Compared to Turek's solution, both the amplitude and the frequency are 10 % below the reference value. This result is considered acceptable since only ten 3-noded triangular elements have been used to discretize the beam thickness.


{|  class="floating_tableSCP wikitable" style="text-align: left; margin: 1em auto;min-width:50%;"
|+ style="font-size: 75%;" |<span id='table-1'></span>Table. 1 PFEM2 and Turek results <span id='citeF-20'></span>[[#cite-20|[20]]]
|- style="border-top: 2px solid;"
| style="border-left: 2px solid;border-right: 2px solid;" |     
| style="text-align: center;border-left: 2px solid;" | Frequency 
| style="text-align: right;border-right: 2px solid;" | Amplitude 
|- style="border-top: 2px solid;"
| style="border-left: 2px solid;border-right: 2px solid;" |  Turek 
| style="text-align: center;border-left: 2px solid;" | 5.3  
| style="text-align: right;border-right: 2px solid;" | 0.0707 
|- style="border-bottom: 2px solid;"
| style="border-left: 2px solid;border-right: 2px solid;" | PFEM2 
| style="text-align: center;border-left: 2px solid;" | 4.8  
| style="text-align: right;border-right: 2px solid;" | 0.0610 

|}

<div id='img-15'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Samper_820077714-turek_snapshot.png|180px|Turek FSI velocity contour]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 15:''' Turek FSI velocity contour
|}

<div id='img-16'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Samper_820077714-graph_turek.png|240px|Tip displacement of the Turek FSI]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 16:''' Tip displacement of the Turek FSI
|}

===6.4 3D benchmark. Dynamic cantilever beam ===

The last example tests the PFEM-2 algorithm in 3D dynamic problems. The same geometry of the previous 2D cantilever beam was selected but restoring the inertial terms and taking the width <math display="inline">b=1</math> to obtain a square section. The geometry can be seen in Figure [[#img-17|17]], with a vertical mass force to initiate the movement. The mesh size is <math display="inline">h=0.125</math> in the surroundings of the beam and <math display="inline">h=0.6</math> in the distant fluid areas, for a total of 283500 4-noded tetrahedra and more than four million particles. The time step <math display="inline">\Delta t</math> is set for <math display="inline">0.001</math> and the properties of the solid are the same as the 2D static cantilever beam example. The subindexes <math display="inline">q</math> and <math display="inline">p</math> are used for the solid and fluid properties respectively

<math display="inline"> \qquad \rho _q = 1  kg/m^3 \qquad E = 10^7 Pa \qquad \nu = 0 </math>   <div id='img-17'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Samper_820077714-3d_beam_geometry2.png|210px|Tridimensional Cantilever geometry]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 17:''' Tridimensional Cantilever geometry
|}

For cases in which the density of the fluid is negligible, it is possible to calculate the analytical vibration frequency in vacuum with the  following expression:

<span id="eq-68"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\omega _{vac} = \sqrt{\frac{C_n^4  E  I}{L^4  \rho _q}}   </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (68)
|}

where the subscript <math display="inline">n</math> denotes the mode and <math display="inline">C_n</math> are the positive roots of the following equation

<span id="eq-69"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>1 + cos (C_n)  cosh (C_n) = 0.   </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (69)
|}

Taking the first vibration mode, <math display="inline">C_1 = 1.8751</math> and the frequency of the beam becomes:

<span id="eq-70"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\omega _{vac} = \sqrt{\frac{C_1^4  E  I}{L^4  \rho _q}} = 32.097  rad/s  \rightarrow  5.108  Hz   </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (70)
|}

The frequency observed in the simulation for low fluid densities (<math display="inline"><0.001kg/m^3</math>) is <math display="inline">5.175 Hz</math>, in good agreement with the theoretical value. The slightly higher value is caused by the low number of elements in the thickness, which leads to a stiffer behavior.

To test the influence of the fluid density, according to Van Eysden <span id='citeF-21'></span>[[#cite-21|[21]]] it is possible to estimate the lower vibration modes for rectangular section beams immersed in a fluid using the following formula:

<span id="eq-71"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\omega _{fluid}=\omega _{vac} \left(1 + \frac{ \pi \rho _p b }{4 \rho _q H } \Gamma _f \right)^{-1/2}    </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (71)
|}

where <math display="inline">H</math> and <math display="inline">b</math> are the height and width of the beam and <math display="inline">\Gamma _f</math> is the hydrodynamic function which depends on the viscosity and the cross section of the beam. For inviscid flows and square sections <math display="inline">\Gamma _f=1.513</math> <span id='citeF-22'></span>[[#cite-22|[22]]].

Fluid densities <math display="inline">\rho _p</math> from <math display="inline">0.0001 kg/m^3 </math> to <math display="inline">0.99 kg/m^3 </math> were selected to validate the algorithm. The results can be seen in Figure [[#img-18|18]], toghether with the comparison against Equation  [[#eq-71|71]]. Analyzing this graph, it is clear that the numerical solution diverges from the expected curve as the fluid density increases. This is caused by the limitations of the standard linear elements used for the velocity field, which lead to interface elements that are assembled with averaged stiffness properties between the fluid and the solid and, therefore, the added mass effect is larger than it should be.   <div id='img-18'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Samper_820077714-beam_oscilations_3d.png|270px|Oscillation frequency of the 3D cantillever beam for different fluid densities]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 18:''' Oscillation frequency of the 3D cantillever beam for different fluid densities
|}

To analyze the convergence of the method with the mesh size, the case with a fluid density <math display="inline">\rho _f=0.5 kg/m^3</math> was selected. Mesh sizes from <math display="inline">h=0.33</math> to <math display="inline">h=0.08</math> were chosen (3 to 12 elements in the thickness of the beam). In all  cases the interface was located in the middle of the elements to test the worst  scenario.

<div id='img-19'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Samper_820077714-beam_oscilations_convergence.png|240px|Mesh convergence for the oscillation frequency of the 3D cantilever beam]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 19:''' Mesh convergence for the oscillation frequency of the 3D cantilever beam
|}

A quadratic convergence with the mesh size can be observed in Figure [[#img-19|19]], indicating that the interaction between the fluid and the beam is accurately simulated.

==7 Conclusions and remarks==

A numerical method to solve monolithically multi-fluids and FSI problems was developed. The strategy is based on the PFEM-2 technique originally formulated in <span id='citeF-1'></span>[[#cite-1|[1]]] for homogeneous fluids, which allows for fast computations while maintaining accuracy.

The enhanced PFEM-2 algorithm makes use of Lagrangian particles and a fixed mesh. Several enrichment degrees of freedom (DoFs) are added in the mesh to improve the pressure field at the interfaces. Adding extra DoFs allows for a more accurate capturing of the solution where the material properties change abruptly, since this leads to discontinuities of the pressure or its gradients. On the other hand, to avoid resizing the system of equations, the new variables are statically condensed and, hence, the computation times are not significantly increased.

The 3D example presented shows that the standard velocity field might not be accurate for coarse meshes due to the artificial added mass effect caused by the FEM linear shape functions.  As a possible solution, enrichening the velocity variables to improve the definition of the interface would reduce the errors without the need for refining the mesh.

More advanced constitutive models are also needed to account for a more realistic simulation of complex phenomena such as debris flows, both at the solid stage and the rheological behaviour of the flow. This would allow to simulate all the phenomena in this type of problem, from initiation to deposition of the eroded material using a single and accurate tool.

==8 Acknowledgements==

This research has been partly funded by  the People Programme (Marie Curie Actions) of the European Union's Seventh Framework Programme FP7/2007-2013/ under REA grant agreement n° 289911. This work was  also supported by the ERC Advanced Grant ‘REALTIME’ project (AdG-2009325), the ERC Advanced Grant ‘SAFECON’ project (AdG-26752) and  the HFLUIDS project of the National RTD Plan of the Spanish Ministry of Science and Innovation (BIA2010-15880).

The authors would also like to thank Dr. Riccardo Rossi for suggesting the special shape functions used in this work.

==BIBLIOGRAPHY==

<div id="cite-1"></div>
'''[[#citeF-1|[1]]]''' Idelsohn, S. and Nigro, N. and Gimenez, J. and Rossi, R. and Marti, J. (2013) "A fast and accurate method to solve the incompressible Navier-Stokes equations", Volume 30. Emerald Group Publishing Limited. Engineering Computations 2 2&#8211;2

<div id="cite-2"></div>
'''[[#citeF-2|[2]]]''' Idelsohn, Sergio R and Marti, Julio and Becker, Pablo and Oñate, Eugenio. (2014) "Analysis of multifluid flows with large time steps using the particle finite element method". Wiley Online Library. Int. J. Numer. Meth. Fluids

<div id="cite-3"></div>
'''[[#citeF-3|[3]]]''' Idelsohn, Sergio R and Marti, Julio and Limache, A and Oñate, Eugenio. (2008) "Unified Lagrangian formulation for elastic solids and incompressible fluids: application to fluid&#8211;structure interaction problems via the PFEM", Volume 197. Elsevier. Comput. Method. Appl. M. 19 1762&#8211;1776

<div id="cite-4"></div>
'''[[#citeF-4|[4]]]''' DeBlois, Bruce M. (1997) "Linearizing convection terms in the Navier-Stokes equations", Volume 143. Elsevier. Comput. Method. Appl. M. 3 289&#8211;297

<div id="cite-5"></div>
'''[[#citeF-5|[5]]]''' Belytschko, Ted and Liu, Wing Kam and Moran, Brian and Elkhodary, Khalil. (2013) "Nonlinear Finite Elements for Continua and Structures". John Wiley & Sons

<div id="cite-6"></div>
'''[[#citeF-6|[6]]]''' Oñate, Eugenio and Franci, Alessandro and Carbonell, Josep M. (2014) "A particle finite element method for analysis of industrial forming processes", Volume 54. Springer. Comput. Mech. 1&#8211;23

<div id="cite-7"></div>
'''[[#citeF-7|[7]]]''' Zienkiewicz, O.C. and Taylor, R.L. and Zhu, J. Z. (2006) "The Finite Element Method: Its Basic and Fundamentals", Volume 1. Elsevier, 6th Edition

<div id="cite-8"></div>
'''[[#citeF-8|[8]]]''' Donea, J. and Huerta, H. (2003) "Finite Element Methods for Flow Problems". John Wiley and Sons

<div id="cite-9"></div>
'''[[#citeF-9|[9]]]''' Codina, Ramon and Badia, Santiago. (2006) "On some pressure segregation methods of fractional-step type for the finite element approximation of incompressible flow problems", Volume 195. Elsevier. Comput. Method. Appl. M. 23 2900&#8211;2918

<div id="cite-10"></div>
'''[[#citeF-10|[10]]]''' Tezduyar, Tayfun E. (1992) "Stabilized finite element formulations for incompressible flow computations", Volume 28. Elsevier. Adv. Appl. Mech. 1&#8211;44

<div id="cite-11"></div>
'''[[#citeF-11|[11]]]''' Codina, Ramon. (2000) "Stabilization of incompressibility and convection through orthogonal sub-scales in finite element methods", Volume 190. Elsevier. Comput. Method. Appl. M. 13 1579&#8211;1599

<div id="cite-12"></div>
'''[[#citeF-12|[12]]]''' Oñate, Eugenio and Franci, Alessandro and Carbonell, Josep M. (2014) "Lagrangian formulation for finite element analysis of quasi-incompressible fluids with reduced mass losses", Volume 74. Wiley Online Library. Int. J. Numer. Meth. Fluids 10 699&#8211;731

<div id="cite-13"></div>
'''[[#citeF-13|[13]]]''' Coppola, H. (2009) "A finite element model for free surface and two fluid flows on fixed meshes". Universitat Politecnica de Catalunya

<div id="cite-14"></div>
'''[[#citeF-14|[14]]]''' Felippa, C. A. (2004) "Introduction to finite element methods". University of Colorado, Boulder, http://www. colorado. edu/engineering/CAS/courses. d/IFEM. d

<div id="cite-15"></div>
'''[[#citeF-15|[15]]]''' Ausas, Roberto F and Buscaglia, Gustavo C and Idelsohn, Sergio R. (2012) "A new enrichment space for the treatment of discontinuous pressures in multi-fluid flows", Volume 70. Wiley Online Library. Int. J. Numer. Meth. Fluids 7 829&#8211;850

<div id="cite-16"></div>
'''[[#citeF-16|[16]]]''' Wendland, Holger. (1995) "Piecewise polynomial, positive definite and compactly supported radial functions of minimal degree", Volume 4. Springer. Adv. Comput. Math. 1 389&#8211;396

<div id="cite-17"></div>
'''[[#citeF-17|[17]]]''' Dadvand, Pooyan and Rossi, Riccardo and Oñate, Eugenio. (2010) "An object-oriented environment for developing finite element codes for multi-disciplinary applications", Volume 17. Springer. Arch. Comput. Method E. 3 253&#8211;297

<div id="cite-18"></div>
'''[[#citeF-18|[18]]]''' He, Xiaoyi and Chen, Shiyi and Zhang, Raoyang. (1999) "A lattice Boltzmann scheme for incompressible multiphase flow and its application in simulation of Rayleigh&#8211;Taylor instability", Volume 152. Elsevier. J. Comput. Phys. 2 642&#8211;663

<div id="cite-19"></div>
'''[[#citeF-19|[19]]]''' Guermond, J-L and Quartapelle, L. (2000) "A projection FEM for variable density incompressible flows", Volume 165. Elsevier. J. Comput. Phys. 1 167&#8211;188

<div id="cite-20"></div>
'''[[#citeF-20|[20]]]''' Turek, Stefan and Hron, Jaroslav. (2006) "Proposal for numerical benchmarking of fluid-structure interaction between an elastic object and laminar incompressible flow". Springer

<div id="cite-21"></div>
'''[[#citeF-21|[21]]]''' Van Eysden, Cornelis A and Sader, John E. (2006) "Resonant frequencies of a rectangular cantilever beam immersed in a fluid", Volume 100. AIP Publishing. J. Appl. Phys. 11 114916

<div id="cite-22"></div>
'''[[#citeF-22|[22]]]''' Brumley, Douglas R and Willcox, Michelle and Sader, John E. (2010) "Oscillation of cylinders of rectangular cross section immersed in fluid", Volume 22. AIP Publishing. Phys. Fluids (1994-present) 5 052001