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

|}

|}

​

Equation [[#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.

​

Equations [[#2.2 Continuum equations|(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)]], 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.[[#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.[[#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;" | (19)

|}

​

In Eq. [[#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|380px|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 [[#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|338px|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|338px|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|410px|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|320px|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|320px|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:

​

'''Step 1)''' Convect the particles using the streamlines}<br />

'''Step 2)''' Project information on the mesh, with  <math>q</math> fluid particles and <math>r</math> solid particles

​

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>     Common~~variables \qquad \qquad  \qquad Solid~~Variables   </math>

|-

| style="text-align: center;padding-top:10px;" | <math>    \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>Fluid~~variables \qquad \qquad \qquad \nu_{solid} = \frac{\sum _{i}^{r} \nu_i }{ \sum _{i}^{r} 1}</math>

|-

| style="text-align: center;" | <math>\mu_{fluid} = \frac{\sum_{i}^{q} \mu_i }{ \sum_{i}^{q} 1} \qquad \qquad\qquad \sigma^{'n}_{solid} = \frac{\sum_{i}^{r} {\sigma '}^n_i }{ \sum _{i}^{r} 1}         </math>

|-

| style="text-align: center;" | <math>\rho _{fluid} = \frac{ \sum _{i}^{q} \rho _i }{ \sum _{i}^{q} 1} \qquad\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>

|}

|}

​

'''Step 3)''' Detect interface elements<br />

'''Step 4)''' Estimate pressure enrichments<br />

'''Step 5)''' Assemble the system with enrichments<br />

'''Step 6)''' Solve the system of equations<br />

'''Step 7)''' Recover condensed DoFs<br />

'''Step 8)''' Update the velocity, pressure and stress at the particles

​

==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|440px|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|50px|rayleigh_re256_1,0]]

| style="padding-left:5px;"|[[Image:Draft_Samper_820077714-rayleigh_re256_1,5.png|50px|rayleigh_re256_1,5]]

| style="padding-left:5px;"|[[Image:Draft_Samper_820077714-rayleigh_re256_1,75.png|50px|rayleigh_re256_1,75]]

| style="padding-left:5px;"|[[Image:Draft_Samper_820077714-rayleigh_re256_2,0.png|50px|rayleigh_re256_2,0]]

| style="padding-left:5px;"|[[Image:Draft_Samper_820077714-rayleigh_re256_2,25.png|50px|rayleigh_re256_2,25]]

| style="padding-left:5px;"|[[Image:Draft_Samper_820077714-rayleigh_re256_2,5.png|50px|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: