Lagrangian formulation for fluid-structure interaction problems using the particle finite element method

==Abstract==

We present a general formulation for analysis of fluid-structure interaction problems using the particle finite element method (PFEM). The key feature of the PFEM is the use of a Lagrangian description to model the motion of nodes (particles) in both the fluid and the structure domains. Nodes are thus viewed as particles which can freely move and even separate from the main analysis domain representing, for instance,  the effect of water drops. A mesh connects the nodes defining the discretized domain where the governing equations, expressed in an integral from, are solved as in the standard FEM. The necessary stabilization for dealing with the incompressibility of the fluid is introduced via the finite calculus (FIC) method. A fractional step scheme for the transient coupled fluid-structure solution is described. Examples of application of the PFEM  to solve a number of fluid-structure interaction problems involving large motions of the free surface and splashing of waves are presented.

'''Keyword''': Lagrangian formulation, fluid-structure, particle finite element.

==1 Introduction==

There is an increasing interest in the development of robust and efficient numerical methods for analysis of engineering problems involving the interaction of fluids and structures accounting for large motions of the fluid free surface and the existence of  fully  or partially submerged bodies. Examples of this kind are common  in  ship hydrodynamics, off-shore structures, spill-ways in dams, free surface channel flows, liquid containers, stirring reactors, mould filling processes, etc.

The movement of solids in fluids is usually analyzed with the finite element method (FEM) <span id='citeF-20'></span>[[#cite-20|[20]]] using the so called arbitrary Lagrangian-Eulerian (ALE) formulation <span id='citeF-3'></span>[[#cite-3|[3]]]. In the ALE approach the movement of the fluid particles is decoupled from that of the mesh nodes. Hence the relative velocity between mesh nodes and particles is used as the convective velocity in the momentum equations.

Typical difficulties of FSI analysis using the FEM with both the Eulerian and ALE formulation include the treatment of the convective terms and the incompressibility constraint in the fluid equations, the modelling and tracking of the free surface in the fluid, the transfer of information between the fluid and solid domains via the contact interfaces, the modelling of wave splashing, the possibility to deal with large rigid body motions of the structure within the fluid domain, the efficient updating of the finite element meshes for both the structure and the fluid, etc.

Most of these problems disappear if a ''Lagrangian description'' is used to formulate the governing equations of both the solid and the fluid domain. In the Lagrangian formulation the motion of the individual particles are followed and, consequently, nodes in a finite element mesh can be viewed as moving “particles”. Hence, the motion of the  mesh discretizing the total domain (including both the fluid and solid parts) is followed during the transient solution.

In this paper we present an overview of a particular class of Lagrangian formulation developed by the authors to solve problems involving the interaction between fluids and solids in a unified manner. The method, called the ''particle finite element method'' (PFEM), treats the mesh nodes in the fluid and solid domains as  particles which can freely move and even separate from the main fluid domain representing, for instance, the effect of water drops. A finite element mesh connects the nodes defining the discretized domain where the governing equations are solved in the standard FEM fashion. The PFEM is the natural evolution of recent work of the authors for  the solution of FSI problems using Lagrangian finite element and meshless methods <span id='citeF-1'></span>[[#cite-1|[1]],<span id='citeF-7'></span>[[#cite-7|7]]-<span id='citeF-10'></span>[[#cite-10|10]],<span id='citeF-17'></span>[[#cite-17|17]]-<span id='citeF-19'></span>[[#cite-19|19]]].

An obvious advantage of the Lagrangian formulation is that the convective terms disappear from the fluid equations. The difficulty is however transferred to the problem of adequately (and efficiently) moving the mesh nodes. Indeed for large mesh motions remeshing may be a frequent necessity along the time solution. We use an innovative mesh regeneration procedure blending elements of different shapes using an extended Delaunay tesselation <span id='citeF-7'></span>[[#cite-7|[7]]-<span id='citeF-9'></span>[[#cite-9|9]]]. Furthermore, this special polyhedral finite element needs special shape functions. In this paper, meshless finite element (MFEM) shape functions have been used <span id='citeF-7'></span>[[#cite-7|[7]]].

The need to properly treat the incompressibility condition in the fluid still remains in the Lagrangian formulation. The use of standard finite element interpolations may lead to a volumetric locking defect unless some precautions are taken <span id='citeF-3'></span>[[#cite-3|[3]],<span id='citeF-20'></span>[[#cite-20|20]]].  In our work the stabilization via a finite calculus (FIC) procedure has been chosen <span id='citeF-12'></span>[[#cite-12|[12]]]. Recent applications of the FIC method for incompressible flow analysis using linear triangles and tetrahedra are reported in <span id='citeF-5'></span>[[#cite-5|[5]],<span id='citeF-12'></span>[[#cite-12|12]],<span id='citeF-13'></span>[[#cite-13|13]],<span id='citeF-14'></span>[[#cite-14|14]],<span id='citeF-15'></span>[[#cite-15|15]],<span id='citeF-18'></span>[[#cite-18|18]],<span id='citeF-19'></span>[[#cite-19|19]]].

The layout of the paper is the following. In the next section the basic ideas of the PFEM are outlined. Next the basic equation for an incompressible flow using a Lagrangian description and the FIC formulation are presented. Then a fractional step scheme for the transient solution via standard finite element procedures is described. Details of the treatment of the coupled FSI problem are given. The procedures for mesh generation and for identification of the free surface nodes are briefly outlined. Finally, the efficiency of the ''particle finite element method'' (PFEM) is shown in its application to a number of FSI problems involving large flow motions, surface waves, moving bodies. etc.

==2 The basis of the Particle Finite Element Method==

Let us consider a domain containing both fluid and solid subdomains. The moving fluid particles interact with the solid boundaries thereby inducing the deformation of the solid which in turn affects the flow motion and, therefore, the problem is  fully coupled.

In the PFEM approach presented here, both the fluid and the solid domains are modelled using an ''updated Lagrangian formulation''. That is, all variables in the fluid and solid domains are assumed to be known in the current configuration at time t. The new set of variables in both domains are sought for in the next or updated configuration at time <math>t+\Delta t</math> (Figure [[#img-1|1]]). The finite element method (FEM) is used to solve the continuum equations in both domains. Hence a mesh discretizing these domains must be generated in order to solve the governing equations for both the fluid and solid problems in the standard FEM fashion. We note once more that the nodes discretizing the fluid and solid domains can be viewed as material particles which motion is tracked during the transient solution. This is useful to model the separation of water particles from the main fluid domain and to follow their subsequent motion as individual particles with an initial velocity and subject to gravity forces.

It is important to note once more that each particle is a material point characterized
by the density of the solid or fluid domain to which it belongs. The mass of a
given domain is obtained by integrating the density at the different material points
over the domain.

The quality of the numerical solution  depends on the discretization chosen as in the standard FEM. Adaptive mesh refinement techniques can be used to improve the solution in zones where large motions of the fluid or the structure occur.

<span id='img-1'></span>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[File:Draft_Samper_662632188_1293_im1.JPG|400px]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figura 1''': Updated lagrangian description for a continuum containing a fluid and a solid domain
|}

===2.1 Basic steps of the PFEM===

For clarity purposes we will define the ''collection or cloud of nodes'' (<math>C</math>) pertaining to the fluid and solid domains, the volume (<math>V</math>) defining the analysis domain for the fluid and the solid and the mesh (<math>M</math>) discretizing both domains.
A typical solution with the PFEM involves the following steps.

<ol>
<div id="step-1"></div>
<li>The starting point at each time step is the cloud of points in the fluid and solid domains. For instance <math>^nC</math> denotes the cloud at time <math>t = tn</math> (Figure [[#img-2|2]]).</li>

<li>Identify the boundaries for both the fluid and solid domains defining the analysis domain <math>^nV</math> in the fluid and the solid. This is an essential step as some boundaries (such as the free surface in fluids) may be severely distorted during the solution process including separation and re-entering of nodes. The Alpha Shape method <div id="citeF-4"></div>[[#cite-4|4]] is used for the boundary definition (see Section [[#7 Identification of Boundary Surfaces|7]]). </li>

<li>Discretize the fluid and solid domains with a finite element mesh <math>^nM</math>. In our work we use an innovative mesh generation scheme based on the extended Delaunay tesselation (Section [[#6 Generation of a new mesh|6]])<div id="citeF-7"></div><div id="citeF-8"></div><div id="citeF-10"></div>[[#cite-7|[7]],[[#cite-8|8]],[[#cite-10|10]]]</li>
<div id="step-4"></div>
<li>Solve the coupled Lagrangian equations of motion for the fluid and the solid domains. Compute the relevant state variables in both domains at the next (updated) configuration for <math>t+\Delta t</math>: velocities, pressure and viscous stresses in the fluid and displacements, stresses and strains in the solid. An overview of the coupled FSI algorithm is given in the next section. </li>

<li>Move the mesh nodes to a new position <math>^{n+1}C</math> where <math>n + 1</math> denotes the time <math>t_n+\Delta t</math>, in terms of the time increment size. This step is typically a consequence of the solution process of step [[#step-4|4]]. </li>

<li>Go back to step [[#step-1|1]] and repeat the solution process for the next time step. </li>

</ol>

<span id='img-2'></span>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[File:Draft_Samper_662632188_8173_img2.JPG|400px]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figura 2''': Sequence of steps to update a “cloud” of nodes from time <math>n (t = t_n)</math> to time <math>n + 1 (t = t_n + \Delta t)</math>
|}

===2.2 Overview of the coupled FSI algoritm===

Figure [[#img-3|3]] shows a typical domain V with external boundaries <math>\Gamma_v</math> and <math>\Gamma_t</math> where the velocity and the surface tractions are prescribed, respectively. The domain <math>V</math> is formed by fluid (<math>V_F</math>) and solid (<math>V_S</math>) subdomains. Both subdomains interact at a common boundary <math>\Gamma_{FS}</math> where the surface tractions and the kinematic variables (displacements, velocities and acelerations) are the same for both subdomains. Note that both set of variables (the surface tractions and the kinematic variables) are equivalent in the equilibrium configuration.

<span id='img-3'></span>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[File:Draft_Samper_662632188_8253_img3.JPG|400px]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figura 3''': Split of the analysis domain <math>V</math> into fluid and solid subdomains. Equality of surface tractions and kinematic variables at the common interface
|}

Let us define <math>^tS</math> and <math>^tF</math> the set of variables defining the kinematics and the stress-strain fields at the solid and fluid domains at time <math>t</math>, respectively, i.e.

<span id='eq-1'></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>^tS := [^tx_S,^tu_S,^tv_S,^ta_S,^t\epsilon_S,^t\sigma_S,\cdot \cdot \cdot]^T</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (1)
|}
<span id='eq-2'></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>^tF := [^tx_F,^tu_F,^tv_F,^ta_F,^t\dot{\epsilon}_F,^t\sigma_F,\cdot \cdot \cdot]^T</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (2)
|}

where <math>x</math> is the nodal coordinate vector, <math>u</math>, <math>v</math> and <math>a</math> are the vector of displacements, velocities and accelerations, respectively, <math>\epsilon</math>, <math>\dot{\epsilon}</math> and <math>\sigma</math> are the strain vector, the strainrate (or rate of deformation) vectors and the Cauchy stress vector, respectively and <math>F</math> and <math>S</math> denote the variables in the fluid and solid domains, respectively. In the discretized problem, a bar over these variables will denote nodal values. 

The coupled fluid-structure interaction (FSI) problem of Figure [[#img-3|3]] will be solved using the following conceptual scheme:

# We assume that the variables in the solid and fluid domains at time <math>t (^tS</math> and <math>^tF)</math> are known.
# Solve for the variables at the solid domain at time <math>t+\Delta t (^{t+\Delta t}S)</math> under prescribed surface tractions at the fluid-solid boundary <math>\Gamma_{FS}</math>.
# Solve for the variables at the fluid domain at time <math>t+\Delta t (^{t+\Delta t}F)</math> under prescribed surface tractions at the external boundary <math>\Gamma_t</math> and prescribed velocities at the external and internal boundaries <math>\Gamma_V</math> and <math>\Gamma_{FS}</math>, respectively.

Iterate between 2 and 3 until convergence.

The variables at the solid domain <math>^{t+\Delta t}S</math> are found via the integration of the dynamic equations of motion in the solid written as

<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>M_Sa_S + g_S - f_S = 0</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (3)
|}

where <math>M_S</math>, <math>g_S</math> and <math>f_S</math> denote the mass matrix, the internal node force vector and the external nodal force vector at the solid domain. The time integration of Eq.([[#eq-3|3]]) is performed using a standard Newmark method. An incremental iterative scheme is implemented within each time step to account for non linear geometrical and material effects.

The FEM solution of the variables in the (incompressible) fluid domain implies solving the momentum and incompressibility equations. As mentioned above this is not such as simple problem as the incompressibility condition limits the choice of the FE approximations for the velocity and pressure to overcome the well known divstability conditions <span id='citeF-20'></span>[[#cite-20|[20]]]. In our work we use a stabilized FEM based on the Finite Calculus approach which allows to use a linear approximation for the velocity and pressure variables. Details of the FEM/FIC formulation used are given in the next section.

Figure [[#img-4|4]] shows a typical example of a PFEM solution in 2D. The pictures correspond to the analysis of the problem of breakage of a water column <span id='citeF-10'></span>[[#cite-10|[10]],<span id='citeF-13'></span>[[#cite-13|13]]]. Figure [[#img-4|4a]] shows the initial grid of four node rectangles discretizing the fluid domain and the solid walls. Boundary nodes identified with the Alpha-Shape method have been marked with a circle. Figures [[#img-4|4b]] and [[#img-4|4c]] show the mesh for the solution at two later times.

<span id='img-4'></span>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[File:Draft_Samper_662632188_7522_img4.JPG]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figura 4''': Breakage of a water column. (a) Discretization of the fluid domain and the solid walls. Boundary nodes are marked with circles. (b) and (c) Mesh in the fluid and solid domains at two different times.
|}

==3 FIC/FEM formulation for a Lagrangian incompressible fluid==

The standard infinitesimal  equations for a viscous incompressible fluid can be written in a  Lagrangian frame as <span id='citeF-11'></span>[[#cite-11|[11]],<span id='citeF-20'></span>[[#cite-20|20]]].

'''Momentum'''
<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>r_{m_i} =0 \quad \hbox{in } \Omega  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (4)
|}

'''Mass balance'''
<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>r_d =0 \quad \hbox{in } \Omega  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (5)
|}

where
<span id='eq-6'></span><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>r_{m_i} = \rho {\partial v_i \over \partial t} + {\partial \sigma _{ij} \over \partial x_j}-b_i\quad ,\quad \sigma _{ji}=\sigma _{ij}</math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (6)
|-
| style="text-align: center;" | <math> r_d = {\partial v_i \over \partial x_i}\qquad i,j = 1, n_d </math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (7)
|}
|}

Above <math display="inline">n_d</math> is the number of space dimensions, <math display="inline">v_i</math> is the velocity along the ith global axis (<math display="inline">v_i = {\partial u_i \over \partial t}</math>, where <math display="inline">u_i</math> is the ''i''th displacement), <math display="inline">\rho </math> is the (constant) density of the fluid, <math display="inline">b_i</math> are the body forces, <math display="inline">\sigma _{ij}</math> are the total stresses given by <math display="inline">\sigma _{ij}=s_{ij}-\delta _{ij}p</math>, <math display="inline">p</math> is the absolute pressure (defined positive in compression) and <math display="inline">s_{ij}</math> are the viscous deviatoric stresses related to the viscosity <math display="inline">\mu </math> by the standard expression
<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>s_{ij}=2\mu \left(\dot \varepsilon _{ij} - \delta _{ij} {1\over 3} {\partial v_k \over \partial x_k}\right) </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (8)
|}

where <math display="inline">\delta _{ij}</math> is the Kronecker delta and the strain rates <math display="inline">\dot \varepsilon _{ij}</math> are
<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>\dot \varepsilon _{ij}={1\over 2} \left({\partial v_i \over \partial x_j}+{\partial v_j \over \partial x_i}\right) </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (9)
|}

In the above all variables are defined at the current time <math display="inline">t</math> (current configuration).

In our work we will solve a ''modified set of governing'' equations derived using a finite calculus (FIC) formulation. The FIC governing equations are <span id='citeF-11'></span>[[#cite-11|[11]]-<span id='citeF-13'></span>[[#cite-13|13]],<span id='citeF-15'></span>[[#cite-15|15]]].

'''Momentum'''
<span id='eq-10'></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>r_{m_i} - \underline{{1\over 2} h_j{\partial r_{m_i} \over \partial x_j}}=0 \qquad \hbox{in } \Omega  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (10)
|}

'''Mass balance'''
<span id='eq-11'></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>r_d - \underline{{1\over 2} h_j {\partial r_d \over \partial x_j}}=0 \qquad \hbox{in }\Omega  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (11)
|}

The problem definition is completed with the following boundary conditions
<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>n_j \sigma _{ij} -t_i + \underline{{1\over 2} h_j n_j r_{m_i}}=0 \quad \hbox{on }\Gamma _t </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (12)
|}
<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>v_j - v_j^p =0 \quad \hbox{on }\Gamma _v </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (13)
|}

and the initial condition is <math display="inline">v_j =v_j^0</math> for <math display="inline">t=t_0</math>. The standard summation convention for repeated indexes is assumed unless otherwise specified.

In Eqs.([[#eq-12|12]]) and ([[#eq-13|13]]) <math display="inline">t_i</math> and <math display="inline">v_j^p</math> are surface tractions and prescribed velocities on the boundaries <math display="inline">\Gamma _t</math> and <math display="inline">\Gamma _v</math>, respectively, <math display="inline">n_j</math> are the components of the unit normal vector to the boundary.

The <math display="inline">h_i's</math> in above equations are ''characteristic lengths'' of the domain where balance of momentum and mass is enforced. In Eq.([[#eq-9|9]]) these lengths define the domain where equilibrium of boundary tractions is established. Details of the derivation of Eqs.([[#eq-10|10]])-([[#eq-13|13]]) can be found in <span id='citeF-11'></span>[[#cite-11|[11]]-<span id='citeF-13'></span>[[#cite-13|13]]].

Eqs.([[#eq-10|10]])-([[#eq-13|13]]) are the starting  point for deriving stabilized finite element methods to solve the incompressible Navier-Stokes equations in a Lagrangian frame of reference using equal order interpolation for the velocity and pressure variables <span id='citeF-1'></span>[[#cite-1|[1]]<span id='citeF-6'>,</span>[[#cite-6|6]]-<span id='citeF-8'></span>[[#cite-8|8]],<span id='citeF-10'></span>[[#cite-10|10]],<span id='citeF-17'></span>[[#cite-17|17]]]. Application of the FIC formulation to finite element and meshless analysis of fluid flow problems can be found in <span id='citeF-5'></span>[[#cite-5|[5]]<span id='citeF-12'>,</span>[[#cite-12|12]],<span id='citeF-13'></span>[[#cite-13|13]],<span id='citeF-14'></span>[[#cite-14|14]],<span id='citeF-15'></span>[[#cite-15|15]],<span id='citeF-17'></span>[[#cite-17|17]],<span id='citeF-18'></span>[[#cite-18|18]]].

===3.1 Transformation of the Mass Balance Equation. Integral Governing Equations===

The underlined term in Eq.([[#eq-8|8]]) can be expressed in terms of the momentum equations. The new expression for the mass balance equation is <span id='citeF-12'></span>[[#cite-12|[12]],<span id='citeF-19'></span>[[#cite-19|19]]]
<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>r_d - \sum \limits _{i=1}^{n_d} \tau _i {\partial r_{m_i} \over \partial x_i}=0 </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (14)
|}

with
<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>\tau _i = {3h_i^2\over 8\mu } </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (15)
|}

At this stage it is no longer necessary to retain the stabilization terms in the momentum equations. These terms are critical in Eulerian formulations to stabilize the numerical solution for high values of the convective terms. In the Lagrangian formulation the convective terms dissappear from the momentum equations and the FIC terms in these  equations are just useful to derive the form of the mass balance equation given by Eq.([[#eq-11|11]]) and can be disregarded there onwards. Consistently, the stabilization terms are also neglected in the Neuman boundary conditions (eqs.[[#eq-9|9]]).

The weighted residual expression of the final form of the momentum and mass balance equations can  be written as
<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>\int _{VF} \delta v_i r_{m_i} dV + \int _{\Gamma _t} \delta v_i (n_j\sigma _{ij} - t_i) d\Gamma =0 </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (16)
|}
<span id='eq-17'></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\int _{VF} q \left[r_d - \sum \limits _{i=1}^{n_d} \tau _i {\partial r_{m_i} \over \partial x_i}\right]dV =0 </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (17)
|}

where <math display="inline">\delta v_i</math> and <math display="inline">q</math> are arbitrary weighting functions equivalent to virtual velocity and virtual pressure fields.

The <math display="inline">r_{m_i}</math> term  in Eq.([[#eq-14|14]]) and the deviatoric stresses and the pressure terms within <math display="inline">r_{m_i}</math> in Eq.([[#eq-13|13]]) are integrated by parts to give
<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>\int _{VF} \left[\delta v_i\rho  {\partial v_i \over \partial t}  + \delta \dot \varepsilon _{ij}(s_{ij}- \delta _{ij}p )\right]d\Omega -  \int _{VF} \delta v_i b_i d\Omega - \int _{\Gamma _t} \delta v_i t_id\Gamma  =0 </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (18)
|}
<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>\int _{VF} q {\partial v_i \over \partial x_i} dV + \int _{VF} \left[\sum \limits _{i=1}^{n_d} \tau _i {\partial q \over \partial x_i} r_{m_i}\right]dV =0 </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (19)
|}

In Eq.([[#eq-18|18]]) <math display="inline">\delta \dot \varepsilon _{ij}</math> are virtual strain rates. Note that the boundary term resulting from the integration by parts of <math display="inline">r_{m_i}</math> in Eq.([[#eq-19|19]]) has been neglected  as the influence of this term in the numerical solution has been found to be negligible.

===3.2 Pressure Gradient Projection===

The computation of the residual terms in Eq.([[#eq-19|19]]) can be simplified if we introduce  the pressure gradient projections <math display="inline">\pi _i</math>, defined as
<span id='eq-20'></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\pi _i = r_{m_i} - {\partial p \over \partial x_i} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (20)
|}

We  express now <math display="inline">r_{m_i}</math> in  Eq.([[#eq-20|20]]) in terms of the <math display="inline">\pi _i</math>  which then become additional variables. The system of integral equations is therefore augmented in the necessary number of  equations by imposing that the residual <math display="inline">r_{m_i}</math> vanishes within the analysis domain (in an average sense). This gives the final system of governing equation as:
<span id='eq-21'></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\int _{VF} \left[\delta v_i\rho {\partial v_i \over \partial t} + \delta \dot \varepsilon _{ij}(s_{ij}- \delta _{ij}p )\right]dV -  \int _{VF} \delta v_i b_i dV - \int _{\Gamma _t} \delta v_i t_id\Gamma =0 </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (21)
|}
<span id='eq-22'></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\int _{VF} q {\partial v_i \over \partial x_i} dV + \int _{VF} \sum \limits _{i=1}^{n_d} \tau _i {\partial q \over \partial x_i} \left({\partial p \over \partial x_i}+\pi _i\right)dV =0 </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (22)
|}
<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>\int _{VF} \delta \pi _i \tau _i \left({\partial p \over \partial x_i}+\pi _i\right)dV =0\qquad \hbox{no sum in }i </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (23)
|}

with <math display="inline">i,j,k=1,n_d</math>.  In Eqs.([[#eq-23|23]]) <math display="inline">\delta \pi _i</math> are appropriate weighting functions and the <math display="inline">\tau _i</math> weights are introduced for symmetry reasons.

===3.3 Finite Element Discretization===

We choose equal order <math display="inline">C^\circ </math> continuous  interpolations of the velocities, the pressure and the pressure gradient projections <math display="inline">\pi _i</math> over each element with <math display="inline">n</math> nodes. The interpolations are written as
<span id='eq-24'></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>v_i = \sum \limits _{j=1}^n N_j \bar v_i^j \quad , \quad p = \sum \limits _{j=1}^n N_j \bar p^j\quad , \quad \pi _i = \sum \limits _{j=1}^n N_j \bar \pi _i^j </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (24)
|}

where <math display="inline">\bar {(\cdot )}^j</math> denotes nodal variables and <math display="inline">N_j</math> are the  shape functions <span id='citeF-20'></span>[[#cite-20|[20]]]. More details of the mesh discretization process and the choice of shape functions are given in Section [[#7 Identification of Boundary Surfaces|7]].

Substituting the approximations ([[#eq-24|24]]) into Eqs.([[#eq-22|22]]-[[#eq-23|23]]) and choosing a Galerkin form with <math display="inline">\delta v_i =q=\delta \pi _i =N_i</math> leads to the following system of discretized equations
<span id='eq-25'></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>{\boldsymbol M}\dot{\bar{\boldsymbol v}} + {\boldsymbol K} \bar {\boldsymbol v} - {\boldsymbol G} \bar {\boldsymbol p}- {\boldsymbol f}={\boldsymbol 0}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (25a)
|}
<span id='eq-25b'></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\displaystyle {\boldsymbol G}^T \bar {\boldsymbol v} + {\boldsymbol L}\bar {\boldsymbol p}+{\boldsymbol Q}\bar {\boldsymbol \pi }={\boldsymbol 0}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (25b)
|}
<span id='eq-25c'></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\displaystyle {\boldsymbol Q}^T \bar {\boldsymbol p} + \hat {\boldsymbol M}\bar {\boldsymbol \pi }={\boldsymbol 0}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (25c)
|}

The matrices and vectors in Eqs.([[#eq-25|25]]) are assembled from the element contributions given by (for 2D problems)

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\displaystyle {\boldsymbol M}_{ij}= \int _{V_F^e} \rho{\boldsymbol  N}_i {\boldsymbol N}_j dV \quad ,\quad  {\boldsymbol K}_{ij} =\int _{V_F^e} {\boldsymbol B}_i^T {\boldsymbol D} {\boldsymbol B}_j dV </math>
|-
| style="text-align: center;" | <math> {\boldsymbol D} =\mu \left[\begin{matrix}2 &0&0\\ 0&2&0\\ 0&0&1\\\end{matrix}\right]\quad ,\quad  {\boldsymbol B}_i =\left[\begin{matrix}\displaystyle {\partial N_i \over \partial x_1}&0\\ 0 & \displaystyle {\partial N_i \over \partial x_2}\\ \displaystyle {\partial N_i \over \partial x_2} & \displaystyle {\partial N_i \over \partial x_1}\\\end{matrix}\right]</math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (26)
|-
| style="text-align: center;" | <math> \displaystyle L_{ij}= \int _{V_F^e} \tau _k {\partial N_i \over \partial x_k} {\partial N_j  \over \partial x_k}dV \quad ,\quad \displaystyle {\boldsymbol Q}= [{\boldsymbol Q}^1,{\boldsymbol Q}^2] \quad ,\quad  \displaystyle Q_{ij}^k = \int _{V_F^e}\tau _k {\partial N_i \over \partial x_k} N_jdV</math>
|-
| style="text-align: center;" | 
|-
| style="text-align: center;" | <math> \displaystyle \hat {\boldsymbol M}= \left[\begin{matrix}\hat {\boldsymbol M}^1 &{\boldsymbol 0}\\  {\boldsymbol 0} & \hat {\boldsymbol M}^2\\\end{matrix}\right]\quad ,\quad  \hat {\boldsymbol M}^k_{ij} = \int _{V_F^e} \tau _k N_i N_j dV  \quad ,\quad \displaystyle {\boldsymbol G}_{ij}= \int _{V_F^e} {\boldsymbol B}_i^T {\boldsymbol m} N_j dV </math>
|-
| style="text-align: center;" | <math> \displaystyle {\boldsymbol f}_i = \int _{V_F^e} N_i {\boldsymbol b}dV + \int _{\Gamma ^e_V}N_i {\boldsymbol t} d\Gamma \quad ,\quad  {\boldsymbol b}=[b_1,b_2]^T \quad ,\quad {\boldsymbol t}=[t_1,t_2]^T  </math>
|}
|}

with <math display="inline">i,j=1,n</math> and <math display="inline">k,l=1,2</math>.

In above <math>B</math> is the strain rate matrix and <math display="inline">{\boldsymbol m} = [1,1,0]^T</math> for 2D problems.

===3.4 Fractional Step Method for Fluid-Structure Interaction Analysis===

The starting point of the iterative algorithm are the variables at time n in the fluid domain (<math>^nF</math>). The sought variables are the variables at time <math>n+1 (^{n+1}F)</math>. For the sake of clarity we will skip the upper left index n + 1 for all variables, i.e.
<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>^{n+1}\bar{x}\equiv  \bar{x}; \quad ^{n+1}\bar{p}\equiv  \bar{p}; \quad ^{n+1}\bar{\pi}\equiv  \bar{\pi}; \quad  ^{n+1}\bar{x}\equiv  \bar{x}; \quad etc</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (27)
|}


A simple iterative algorithm is obtained by splitting the pressure from the momentum equations as follows
<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>\bar {\boldsymbol v}^* = \bar {\boldsymbol v}^n -\Delta t {\boldsymbol M}^{-1}[{\boldsymbol K} {\boldsymbol v}^{n+1,j-1} -{\boldsymbol G} \mathbf{p}^n - {\boldsymbol f}^{n+1}]</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (28)
|}
<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>\bar{\boldsymbol v}^{n+1,j}= \bar{\boldsymbol v}^* + \Delta t {\boldsymbol M}^{-1}{\boldsymbol G}\delta \bar{\boldsymbol p}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (29)
|}

where <math display="inline">\delta \bar{\boldsymbol p}</math> denotes a pressure increment. In above equations and in the following the left upper index <math display="inline">n</math> refers to values in the current application <math>^nV</math> whereas the right index <math display="inline">j</math> denotes an iteration number within each time step.

The value of <math display="inline">\bar {\boldsymbol v}^{n+1,j}</math> from Eq.([[#eq-28|28]]) is substituted now into Eq.([[#eq-25b|25b]]) to give
<span id='eq-30'></span>
<span id='eq-30a'></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>{\boldsymbol G}^T\bar {\boldsymbol v}^* + \Delta t {\boldsymbol S}\delta \bar{\boldsymbol p} + {\boldsymbol L} \bar{\boldsymbol p}^{j+1}+ {\boldsymbol Q}\bar {\boldsymbol \pi }^{j}={\boldsymbol 0}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (30a)
|}

where
<span id='eq-30b'></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>{\boldsymbol S} = {\boldsymbol G}^T {\boldsymbol M}^{-1}{\boldsymbol G} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (30b)
|}

Typically matrix <math>S</math> is computed using a diagonal matrix <math display="inline">{\boldsymbol M} = {\boldsymbol    M}_d</math>, where the subscript <math display="inline">d</math> denotes hereonward a diagonal matrix.

An alternative is to approximate matrix <math>S</math> by a Laplacian matrix. This reduces considerably the bandwith of <math>S</math>. The disadvantage is that the pressure increment must be then prescribed  on the free surface and this reduces the accuracy in the satisfaction of the incompressibility condition in these regions (see Remark [[#remark-1|1]]).

A semi-implicit algorithm can  be derived as follows. For each iteration:<br/>
<span id='step-1'></span>
'''Step 1''' Compute <math display="inline">{\boldsymbol v}^*</math>  from Eq.([[#eq-28|28]]) with <math display="inline">{\boldsymbol M}={\boldsymbol    M}_d</math>. For the first iteration <math display="inline">(\bar{v}^1, \bar{p}^1, \bar{\pi}^1, \bar{x}^1)\equiv(^n\bar{v}, ^n\bar{p}, ^n\bar{\pi}, ^n\bar{x})</math>.<br/>
<span id='step-2'></span>
'''Step 2''' Compute <math display="inline">\delta \bar {\boldsymbol p}</math> and <math display="inline">{\boldsymbol p}^{j+1}</math> from Eq.([[#eq-30|30a]]) as
<span id='eq-31'></span>
<span id='eq-31a'></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\delta \bar {\boldsymbol p} =-({\boldsymbol L}+\Delta t  {\boldsymbol S})^{-1} [{\boldsymbol G}^T\bar{\boldsymbol v}^* +{\boldsymbol Q}\bar {\boldsymbol \pi }^{j}+ {\boldsymbol L} \bar {\boldsymbol p}^n]</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (31a)
|}

The pressure <math display="inline">\bar {\boldsymbol p}^{n+1,j}</math> is computed as follows
<span id='eq-31b'></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\bar {\boldsymbol p}^{n+1,j} = \bar{\boldsymbol p}^n + \delta \bar{\boldsymbol p} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (31b)
|}
<span id='step-3'></span>
'''Step 3''' Compute <math display="inline"> \bar{\boldsymbol v}^{n+1,j}</math>  from Eq.([[#eq-29|29]]) with <math display="inline">{\boldsymbol M}={\boldsymbol M}_d</math>
<span id='step-4'></span>
'''Step 4''' Compute <math display="inline">  \bar{\boldsymbol \pi }^{n+1,j}</math>  from Eq.([[#eq-25c|25c]]) as
<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>\bar{\boldsymbol \pi }^{n+1,j}=- \hat {\boldsymbol M}_d^{-1} {\boldsymbol Q}^T \bar {\boldsymbol p}^{n+1,j} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (32)
|}
<span id='step-5'></span>
'''Step 5''' Update the coordinates of the mesh nodes. From the definition of the velocity <math display="inline">vi={\partial u_i\over \partial t}</math> it is deduced.
<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>x^{j+1}_i= ^n{x}_i+\bar{v}^{j+1}_i\Delta t</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (33)
|}
<span id='step-6'></span>
'''Step 6''' Check the convergence of the velocity and pressure fields. If convergence is achieved move to the next time step, otherwise return to step [[#step-1|1]] for the next iteration with <math display="inline">j \gets 1 + j</math>.

Note that solution of steps [[#step-1|1]], [[#step-3|3]] and [[#step-4|4]] does not require the solution of  a system of equations as a diagonal form is chosen for <math display="inline">\boldsymbol M</math> and <math display="inline">\hat {\boldsymbol M}</math>.

In the examples presented in the paper the time increment size has been chosen as
<span id='eq-34'></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\Delta t =\min (\Delta t_i ) \quad \hbox{with}\quad \Delta t_i ={h_i^{\min } \over \vert {\boldsymbol v}\vert } </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (34)
|}

where <math display="inline">h_i^{\min }</math> is the distance between node <math display="inline">i</math> and the closest node in the mesh.

Although not explicitely mentioned all matrices and vectors in Eqs.([[#eq-27|27]])-([[#eq-31|31]]) are computed at the  configuration <math display="inline">^{n+1}V_F</math>. This  means that the integration domain changes for each iteration and, hence,  all the terms involving space derivatives  must be updated at each iteration. An alternative is to refer the integrations domain at each time step to <math display="inline">^nV_F</math>. The jacobian matrix is needed in this case to transform the space derivatives and the differencial of volume from <math display="inline">^{n+1}V_F</math> to <math display="inline">^nV_F</math> at each iteration.

The boundary conditions are applied as follows. No condition is applied for the computation of the fractional velocities <math display="inline">{\boldsymbol v}^*</math> in Eq.([[#25a|25]]). The prescribed velocities at the boundary are applied when solving for <math display="inline">\bar{\boldsymbol v}^{j+1}</math> in step [[#step-3|3]].

<span id='remark-1'></span>
'''Remark 1'''. The form of <math>S</math> in Eq.([[#eq-30b|30b]]) avoids the need to prescribing the pressure at the boundary nodes. It is however recommended to fix the pressure at a point in the analysis domain so as to ensure the correct definition of the pressure field for all problems. An alternative procedure is to approximate  <math>S</math> by a Laplacian matrix <span id='citeF-10'></span>[[#cite-10|[10]],<span id='citeF-13'></span>[[#cite-13|13]]]. The pressure increment on the free surface must be prescribed in this case to the value <math display="inline">2\mu {\partial \Delta \bar v_n\over    \partial n}</math> where <math display="inline">\Delta \bar v_n={\boldsymbol n}^T \Delta \bar{\boldsymbol v}</math> is the velocity increment along the normal direction to the boundary and <math display="inline">\Delta \bar {\boldsymbol v}= \bar{\boldsymbol v}^{n+1,j}- \bar{\boldsymbol v}^n</math>. We have found however that the form of  <math>S</math> given by Eq.([[#eq-30b|30b]]) allows to satisfy better the incompressibility condition in the vecinity of the free surfaces, thereby leading to smaller   volume losses in transient problems involving large motions of the free surface.

==4 Staggered scheme for the FSI problem==
The solution for the variables in the solid and fluid domains at the updated configuration <math>^{n+1}F</math>, <math>^{n+1}s</math> is found using the staggered scheme shown in Box [[#box-1|1]].
<div id='box-1'></div>
{| class="floating_imageSCP" style="text-align: center; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[File:Draft_Samper_662632188_5116_box1.JPG|500px]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Box 1:''' Staggered scheme for the FSI problem
|}

==5 Treatment of contact between the fluid and a fixed boundary==

The motion of the solid is governed by the action of the fluid flow forces induced by the pressure and the viscous stresses acting at the solid boundary, as mentioned above.

The  condition of prescribed velocities  at the solid boundaries in the PFEM are  applied in strong form to the boundary nodes. These nodes might belong to fixed external boundaries or to moving boundaries linked to the interacting solids. Contact between water particles and the solid boundaries is accounted for by the incompressibility condition which ''naturally prevents the penetration of the water nodes into the solid boundaries'' (Figure [[#img-5|5]]). This simple way to treat the water-wall contact is another attractive feature of the PFEM formulation.

<span id='img-5'></span>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[File:Draft_Samper_662632188_7050_img5.JPG|500px]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figura 5''': Automatic treatment of contact condition at the fluid-wall interface
|}

==6 Generation of a new mesh==

One of the key points for the success of the Lagrangian flow formulation  described here is the fast regeneration of a mesh at every time step on the basis of the position of the nodes in the space domain. In our work the mesh is generated using the so called extended Delaunay tesselation (EDT) presented in <span id='citeF-7'></span>[[#cite-7|[7]],<span id='citeF-9'></span>[[#cite-9|9]],<span id='citeF-10'></span>[[#cite-10|10]]]. The EDT allows one to generate non standard meshes combining elements of  arbitrary polyhedrical shapes  (triangles, quadrilaterals and other polygons in 2D and tetrahedra, hexahedra and arbitrary polyhedra in 3D) in a computing time of order <math display="inline">n</math>, where <math display="inline">n</math> is the total number of nodes in the mesh (Figure [[#img-6|6]]). The <math display="inline">C^\circ </math> continuous shape functions of the elements can be simply obtained using the so called meshless finite element interpolation (MFEM). Details of the mesh generation procedure and the derivation of the MFEM shape functions can be found in <span id='citeF-7'></span>[[#cite-7|[7]],<span id='citeF-9'></span>[[#cite-9|9]],<span id='citeF-10'></span>[[#cite-10|10]]].

Once the new mesh has been generated  the numerical solution is found at each time step using the fractional step algorithm described in the previous section.

<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_662632188-Figure3.png|350px|Generation of non standard meshes combining different polygons (in 2D) and polyhedra (in 3D) using the extended Delaunay technique.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 6:''' Generation of non standard meshes combining different polygons (in 2D) and polyhedra (in 3D) using the extended Delaunay technique.
|}

==7 Identification of Boundary Surfaces==

One of the main tasks  in the PFEM is the correct definition of the boundary domain. Sometimes, boundary nodes are explicitly identified  differently from internal nodes. In other cases, the total set of nodes is the only information available and the algorithm must recognize the boundary nodes.

The  extended Delaunay partition makes it easier to recognize boundary nodes. Considering that the nodes follow a variable <math display="inline">h(x)</math>  distribution, where <math display="inline">h(x)</math> is typically the minimum distance between two nodes, the following criterion has been used. ''All nodes on an empty sphere with a radius  greater than <math>\alpha h</math>, are considered as boundary nodes''. In practice <math display="inline">\alpha </math>  is a parameter close to, but greater than one. This criterion is coincident with the Alpha Shape concept <span id='citeF-4'></span>[[#cite-4|[4]]]. Figure [[#img-7|7]] shows example of the boundary recognition using the Alpha Shape technique.

Once a decision has been made concerning which  nodes are on the boundaries, the boundary surface is defined by  all the polyhedral surfaces (or polygons in 2D) having all their nodes on the boundary and belonging to just one polyhedron.

The correct boundary surface is important to define the  normal external to the surface. Furthermore, in weak forms (Galerkin) such as those used here a correct evaluation of the volume domain is also important. In the criterion proposed above, the error in the boundary surface definition is proportional to <math display="inline">h</math> which is an acceptable error.

The method described also allows one to identify isolated fluid particles outside the main fluid domain. These particles are treated as part of the external boundary where the pressure is fixed to the atmospheric value (Figure [[#img-7|7]]). We recall that each particle is a material point characterized by the density of the solid or fluid domain to which it belongs. The mass which is lost when a boundary element is eliminated due to departure of a node (<math display="inline">a</math> particle) from the main analysis domain is again regained when the “flying” node falls donw and a new boundary element is created by the Alpha Shape algorithm when the lost node is at a distance less than <math display="inline">\alpha h</math> from the boundary.

A practical application of the method for identifying free surface particles is shown in Figure [[#img-8|8]]. The example corresponds to the analysis of the motion of a fluid within an oscillating ellipsoidal container.

<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_662632188-Figure5.png|400px|Identification of individual particles (or a group of particles) starting from a given collection of nodes.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 7:''' Identification of individual particles (or a group of particles) starting from a given collection of nodes.
|}

<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_662632188-Figure6b.png|400px|Motion of a liquid within an oscillating container. Position of the liquid particles at two different times. ]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 8:''' Motion of a liquid within an oscillating container. Position of the liquid particles at two different times. 
|}

===7.1 Contact between the solid-solid interfaces===

The contact between two solid interfaces can be simply treated by introducing a layer of ''contact elements'' between the two interacting solid interfaces. This layer is automatically created during the mesh generation by prescribing a minimum distance between two solid boundaries. If the distance exceeds the minimum value then the generated elements are treated as fluid elements. Otherwise the elements are treated as contact elements where a relationship between the tangential and normal forces and the corresponding displacement is introduced so as to model elastic and frictional contact in the normal and tangential directions, respectively (Figure [[#img-9|9]]).

This algorithm has proven to be very effective and it allows to identify and model complex frictional contact conditions between two or more interacting solids in an extremely simple manner. Of course the accuracy on this contact model depends of the critical distance above mentioned.

<div id='img-9'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[File:Draft_Samper_662632188_2696_img-9.JPG|400px|Contact conditions at a solid-solid interface]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 9:''' Contact conditions at a solid-solid interface
|}

==8 Examples==

The examples chosen show the applicability of the PFEM to solve problems involving large fluid motions and FSI situations.

===8.1 Rigid objects filling into water===

The analysis of the motion of submerged or floating objects in water is of great interest in many areas of harbour and coastal engineering and naval architecture among others.

Figure [[#img-10|10]] shows the penetration and evolution of a cube and a cylinder in a container with water. The colours denote the different sizes of the elements at several times. In order to increase the accuracy of the FSI problem smaller size elements have been generated in the vicinity of the moving bodies during their motion (Figure [[#img-11|11]]).

<div id='img-10'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[File:Draft_Samper_662632188_3160_img-10.JPG|400px|2D simulation of the penetration and evolution of a cube and a cylinder in a water container. The colours denote the different sizes of the elements at several times.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 10:''' 2D simulation of the penetration and evolution of a cube and a cylinder in a water container. The colours denote the different sizes of the elements at several times.
|}

<div id='img-11'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|
[[File:Draft_Samper_662632188_9623_img-11.JPG|400px|Detail of element sizes during the motion of a rigid cylinder within   a water container.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 11:''' Detail of element sizes during the motion of a rigid cylinder within   a water container.
|}

===8.2 Impact of water streams on rigid structures===

Figure [[#img-12|12]] shows the simulation of the propagation of a wave generated in a test canal impacting on a vertical wall. Figure [[#img-13|13]] shows a comparison of the experimental and numerical values of the water pressure on the vertical wall. The overall agreement is noticeable.

Figure [[#img-14|14]] shows an example of a wave breaking within a prismatic container including a vertical cylinder. Finally Figure [[#img-15|15]] shows the impact of a wave on a vertical column sustained by four pillars. The objective of this example was to model the impact of a water stream on a bridge pier accounting for the foundation effects.

<div id='img-12'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[File:Draft_Samper_662632188_6304_img-12.JPG|400px|Propagation of a wave generated in a test canal impacting with a   vertical wall.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 12:''' Propagation of a wave generated in a test canal impacting with a   vertical wall.
|}

<div id='img-13'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[File:Draft_Samper_662632188_9901_img-13.JPG|600px|Comparison of experimental and numerical water pressures at the   vertical wall for the example of Figure 7]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 13:''' Comparison of experimental and numerical water pressures at the   vertical wall for the example of Figure [[#img-7|7]].
|}

<div id='img-14'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|
[[File:Draft_Samper_662632188_1676_img-14.JPG|600px|Evolution of a water column within a prismatic container including a   vertical cylinder.]]
|- style="text-align: center; font-size: 75%;"
| colspan="2" | '''Figure 14:''' Evolution of a water column within a prismatic container including a   vertical cylinder.
|}

<div id='img-15'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|
[[File:Draft_Samper_662632188_5434_img-15.JPG|600px|Impact of a wave on a prismatic column on a slab sustained by four   pillars.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 15:''' Impact of a wave on a prismatic column on a slab sustained by four   pillars.
|}

===8.3 Dragging of objects by water streams===

Figure [[#img-16|16]] shows the effect of a wave impacting on a cube representing a vehicle. This situation is typical in  flooding and Tsunami situations.

<div id='img-16'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[File:Draft_Samper_662632188_9287_img-16.JPG|450px|Dragging of a cubic object by a water stream.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 16:''' Dragging of a cubic object by a water stream.
|}

===8.4 Impact of sea waves on breakwater and piers===

Figure [[#img-17|17]] shows the simulation of the impact of a wave generated in an experimental flume on a collection of water motion in the vicinity of the rocks represent a breakwater. Details of the  water-rock interaction are shown in Figure [[#img-18|18]].

Figure [[#img-19|19]] shows a 3D analysis of a similar problem. Figure [[#img-20|20]] shows the 3D simulation of the interaction of a wave with a vertical pier formed by a collection of reinforced concrete cylinders.

The last examples shown in Figures [[#img-21|21]] and [[#img-22|22]] evidence the potential of the PFEM to solve 3D problems involving complex interactions between water  and solid objects. Figure [[#img-21|21]] shows the simulation of the falling of two tetrapods in a water container. Finally, Figure [[#img-22|22]] shows the motion of a collection of ten tetrapods placed in a slope under an incident wave.

<div id='img-17'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[File:Draft_Samper_662632188_6864_img-17.JPG|500px|Generation and impact of a wave on a collection of   rocks in a breakwater.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 17:''' Generation and impact of a wave on a collection of   rocks in a breakwater.
|}

<div id='img-18'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[File:Draft_Samper_662632188_5375_img-18.JPG|600px|Detail of the impact of a wave on a breakwater. The arrows indicate   the water force on the rocks at different instants.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 18:''' Detail of the impact of a wave on a breakwater. The arrows indicate   the water force on the rocks at different instants.
|}

<div id='img-19'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[File:Draft_Samper_662632188_4507_img-19.JPG|500px|3D simulation of the impact of a wave on a collection of rocks in a breakwater.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 19:''' 3D simulation of the impact of a wave on a collection of rocks in a breakwater.
|}

<div id='img-20'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[File:Draft_Samper_662632188_5106_img-20.JPG|500px|Interaction of a wave with a vertical pier formed by   reinforced concrete cylinders.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 20:''' Interaction of a wave with a vertical pier formed by   reinforced concrete cylinders.
|}

<div id='img-21'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[File:Draft_Samper_662632188_1027_img-21.JPG|500px|Motion of two tetrapods falling in a water container.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 21:''' Motion of two tetrapods falling in a water container.
|}

<div id='img-22'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[File:Draft_Samper_662632188_4754_img-22a.JPG|600px|]]
|-
|[[File:Draft_Samper_662632188_3107_img-22b.JPG|600px|Motion of ten tetrapods on a slope under an incident wave.]]
|- style="text-align: center; font-size: 75%;"
| colspan="2" | '''Figure 22:''' Motion of ten tetrapods on a slope under an incident wave.
|}

==10 Conclusions==

The particle finite element method (PFEM) is ideal to treat problems involving fluids with free surface and submerged or floating structures within a unified Lagrangian finite element framework. Problems such as the analysis of fluid-structure interaction, large motion of  fluid or solid particles, surface waves, water splashing, separation of water drops, etc. can be easily solved with the PFEM. The success of the method lies in the accurate and efficient solution of the equations of an incompressible fluid and of solid dynamics using a stabilized finite element method  allowing the use of low order elements with equal order interpolation for all the variables. Other essential solution ingredients are the efficient regeneration of the finite element mesh using an extended Delaunay tesselation, the meshless finite element interpolation (MFEM) and the identification of the boundary nodes using an Alpha Shape type technique. The examples presented have shown the potential of the PFEM for solving a wide class of practical FSI problems.

==Acknowledgements==

Thanks are given to Dr. F. Del Pin, Dr. N. Calvo and Ms. M. de Mier for many useful suggestions.

==References==

<div id="cite-1"></div>
[[#citeF-1|[1]]] Aubry R, Idelsohn SR, Oñate E (2005) Particle finite element   method in fluid mechanics including thermal convection-diffusion. Computer   & Structures 83(17-18):1459&#8211;1475

<div id="cite-2"></div>
[[#citeF-2|[2]]]  Codina R, Zienkiewicz OC (2002) CBS versus GLS stabilization of   the incompressible Navier-Stokes equations and the role of the time step as   stabilization parameter. Communications in Numerical Methods in Engineering   18(2):99&#8211;112

<div id="cite-3"></div>
[[#citeF-3|[3]]]  Donea J, Huerta A (2003) Finite element method for flow problems. J. Wiley.

<div id="cite-4"></div>
[[#citeF-4|[4]]]  Edelsbrunner H, Mucke EP (1999) Three dimensional alpha shapes. ACM   Trans. Graphics 13:43&#8211;72

<div id="cite-5"></div>
[[#citeF-5|[5]]]  García J,  Oñate E (2003) An unstructured finite element   solver for ship hydrodynamic problems. J. Appl. Mech. 70:18&#8211;26 January

<div id="cite-6"></div>
[[#citeF-6|[6]]] Idelsohn SR, Oñate E, Del Pin F,  Calvo N (2002) Lagrangian formulation: the only way to solve some free-surface fluid mechanics problems. Fith World Congress on Computational Mechanics, Mang HA, Rammerstorfer FG, Eberhardsteiner J (eds), July 7&#8211;12, Viena, Austria

<div id="cite-7"></div>
[[#citeF-7|[7]]]  Idelsohn SR, Oñate E, Calvo N, Del Pin F (2003a) The meshless finite element method.  Int. J. Num. Meth. Engng.  58(6):893&#8211;912

<div id="cite-8"></div>
[[#citeF-8|[8]]] Idelsohn SR, Oñate E, Del Pin F (2003b) A lagrangian meshless finite element method applied to fluid-structure interaction problems. Computer and Structures 81:655&#8211;671

<div id="cite-9"></div>
[[#citeF-9|[9]]] Idelsohn SR, Calvo N, Oñate E (2003c) Polyhedrization of an arbitrary point set. Comput. Method Appl. Mech. Engng.  192(22-24):2649&#8211;2668

<div id="cite-10"></div>
[[#citeF-10|[10]]]   Idelsohn SR, Oñate E, Del Pin F (2004) The particle finite   element method: a powerful tool to solve incompressible flows with   free-surfaces and breaking waves. Int. J. Num. Meth. Engng. 61:964-989

<div id="cite-11"></div>
[[#citeF-11|[11]]] Oñate E  (1998) Derivation of stabilized equations for advective-diffusive transport and fluid flow problems. Comput. Meth. Appl. Mech. Engng. 151:233&#8211;267

<div id="cite-12"></div>
[[#citeF-12|[12]]] Oñate E, (2000) A stabilized finite element method for incompressible viscous flows using a finite increment calculus formulation. Comp. Meth. Appl. Mech. Engng. 182(1&#8211;2):355&#8211;370

<div id="cite-13"></div>
[[#citeF-13|[13]]]  Oñate E (2004) Possibilities of finite calculus in computational mechanics. Int. J. Num. Meth. Engng. 60(1):255&#8211;281

<div id="cite-14"></div>
[[#citeF-14|[14]]]  Oñate E,  Idelsohn SR (1998) A mesh free finite point method for advective-diffusive transport and fluid flow problems. Computational Mechanics 21:283&#8211;292

<div id="cite-15"></div>
[[#citeF-15|[15]]]  Oñate E, García J (2001) A finite element method for  fluid-structure interaction with surface waves using a finite calculus formulation. Comput. Meth. Appl. Mech. Engrg. 191:635&#8211;660

<div id="cite-16"></div>
[[#citeF-16|[16]]] Oñate E, Sacco C,  Idelsohn SR (2000) A finite point method for   incompressible flow problems.  Comput. Visual. in Science 2:67&#8211;75

<div id="cite-17"></div>
[[#citeF-17|[17]]]  Oñate E, Idelsohn SR, Del Pin F (2003) Lagrangian formulation for   incompressible fluids using finite calculus and the finite element   method. Numerical Methods for Scientific Computing Variational Problems and   Applications, Y Kuznetsov, P Neittanmaki, O Pironneau (Eds.), CIMNE, Barcelona

<div id="cite-18"></div>
[[#citeF-18|[18]]] Oñate E, García J, Idelsohn SR (2004a) Ship   hydrodynamics. Encyclopedia of Computational Mechanics, E Stein, R de Borst,   T.J.R. Hughes (Eds), J. Wiley.

<div id="cite-19"></div>
[[#citeF-19|[19]]] Oñate E, Idelsohn SR, Del Pin F, Aubry R (2004b) The particle   finite element method. An overview. Int. J. Comput. Methods 1(2):267-307

<div id="cite-20"></div>
[[#citeF-20|[20]]]  Zienkiewicz OC, Taylor RL, Nithiarasu P (2006) The finite element   method for fluid dynamics,   Elsevier

<div id="cite-21"></div>
[[#citeF-21|[21]]]  Zienkiewicz OC, Taylor RL (2005) The finite element method for   solid and structural mechanics. Elsevier