Advances in the particle finite element method for the analysis of fluid-multibody interaction and bed erosion in free surface flows

<!-- metadata commented in wiki content
==ADVANCES IN THE PARTICLE FINITE ELEMENT  METHOD FOR THE ANALYSIS  OF FLUID-MULTIBODY INTERACTION AND BED EROSION IN FREE SURFACE FLOWS  ==

'''Eugenio Oñate, Sergio R. Idelsohn<math>^{\scriptstyle *}</math>, Miguel A. Celigueta

'''and'''   Riccardo Rossi

'''

{|
|-
|International Center for Numerical Methods in Engineering (CIMNE) 
|-
| Technical University of Catalonia 
|-
| Campus Norte UPC, 08034 Barcelona, Spain 
|-
| Email: [mailto:onate@cimne.upc.edu onate@cimne.upc.edu] 
|-
| Web page: [http://www.cimne.com http://www.cimne.com]
|-
| <math>^{\scriptstyle *}</math> ICREA Research Professor at CIMNE
|-
| 
|}
-->

==Abstract==

We present some advances in the formulation of the Particle Finite   Element Method (PFEM) for solving complex fluid-structure   interaction problems with free surface waves. In particular, we present   extensions of   the PFEM for the analysis of the interaction between a collection of bodies in water   allowing for frictional contact conditions at the fluid-solid and   solid-solid interfaces via mesh generation. An  algorithm to treat bed   erosion in free surface flows is   also presented. Examples of application of the PFEM  to solve a number of   fluid-multibody interaction problems involving splashing of waves, large   motions of floating and submerged bodies and bed erosion situations   are given.

'''keywords''' Lagrangian formulation, fluid-structure interaction, particle finite element method, bed erosion, free surface flows

==1 INTRODUCTION==

The analysis of 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 which interact among themselves is of big relevance in many areas of engineering. Examples are common in ship hydrodynamics, off-shore and harbour structures, spill-ways in dams, free surface channel flows, environmental flows, liquid containers, stirring reactors, mould filling processes, etc.

Typical difficulties of fluid-multibody interaction  analysis in free surface flows 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 the moving solid domains via the contact interfaces, the modeling of wave splashing, the possibility to deal with large  motions of the bodies within the fluid domain, the efficient updating of the finite element meshes for both the structure and the fluid, etc. For a comprehensive list of references in FEM for fluid flow problems see <span id='citeF-5'></span><span id='citeF-34'></span>[[#cite-5|[5,34]]] and the references there included. A survey of recent works in fluid-structure interaction analysis can be found in <span id='citeF-16'></span>[[#cite-16|[16]]], <span id='citeF-25'></span>[[#cite-25|[25]]] and <span id='citeF-32'></span>[[#cite-32|[32]]].

Most of the above problems disappear if a ''Lagrangian description'' is used to formulate the governing equations of both the solid and the fluid domains. 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 material points (hereforth called “particles”). Hence, the motion of the  mesh discretizing the total domain (including both the fluid and solid parts) is followed during the transient solution.

The authors have successfully developed in previous works a particular class of Lagrangian formulation for solving problems involving complex interaction between fluids and solids. 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 using a stabilized FEM based in the Finite Calculus (FIC) approach. An 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.  We use a  mesh regeneration procedure blending elements of different shapes using an extended Delaunay tesselation with special shape functions <span id='citeF-9'></span><span id='citeF-11'></span>[[#cite-9|[9,11]]]. The theory and applications of the PFEM are reported in <span id='citeF-2'></span><span id='citeF-4'></span><span id='citeF-9'></span><span id='citeF-10'></span><span id='citeF-12'></span><span id='citeF-13'></span><span id='citeF-24'></span><span id='citeF-25'></span><span id='citeF-26'></span><span id='citeF-28'></span><span id='citeF-29'></span>[[#cite-2|[2,4,9,10,12,13,24,25,26,28,29]]].

The aim of this paper is to describe two recent advances of the PFEM: a) the analysis of the interaction between a collection of bodies which are floating and/or submerged in the fluid, and b) the modeling of bed erosion in open channel flows. Both problems are of great relevance in many areas of civil, marine and naval engineering, among others. It is shown in the paper that the PFEM provides a general analysis methodology for treat such a complex problems in a simple and efficient manner.

The layout of the paper is the following. In the next section the key ideas of the PFEM are outlined. Next the basic equations for an incompressible flow using a Lagrangian description and the FIC formulation are presented. Then a fractional step scheme for the transient solution is briefly described. Details of the treatment of the coupled FSI problem are given. The methods for mesh generation and for identification of the free surface nodes are outlined. The procedure for treating at mesh generation level the contact conditions at fluid-wall interfaces and the frictional contact interaction between moving solids is explained. A methodology  for modeling bed erosion due to fluid forces is described. Finally, the efficiency of the PFEM is shown in its application to a number of problems involving large flow motions, surface waves, moving bodies in water and bed erosion.

==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  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 <math display="inline">t</math>. The new set of variables in both domains are sought for in the ''next or updated configuration'' at time <math display="inline">t+\Delta t</math> (Figure 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. Recall that the nodes discretizing the fluid and solid domains are treated as ''material particles'' which motion is tracked during the transient solution. This is useful to model the separation of fluid particles from the main fluid domain in a splashing wave, or soil particles in a bed erosion problem, and to follow their subsequent motion as individual particles with a known density, an initial acceleration and velocity and subject to gravity forces. 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.

<div id='img-1'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[File:Draft_Samper_677613731_7270_Fig1.png|500px|Updated lagrangian description for a continuum containing a fluid and   a solid domain]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 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 (C)'' pertaining to the fluid and solid domains, the ''volume (V)'' defining the analysis domain for the fluid and the solid and the ''mesh (M)'' discretizing both domains.

A typical solution with the PFEM involves the following steps.

<ol>

<li>The starting point at each time step is the cloud of points in the fluid and solid   domains. For instance <math display="inline">{}^nC</math> denotes the cloud at time <math display="inline">t=t_n</math> (Figure 2). </li>

<li>Identify the  boundaries for both the fluid and solid domains defining   the analysis domain <math display="inline">{}^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, including separation   and re-entering of nodes. The Alpha Shape method   <span id='citeF-6'></span>[[#cite-6|[6]]] is used for the boundary definition   (Section 5). </li>

<li>Discretize the fluid and solid domains with a finite element mesh <math display="inline">{}^nM</math>. In our work we use an innovative mesh generation scheme based on the extended Delaunay tesselation (Section 4) <span id='citeF-9'></span><span id='citeF-10'></span><span id='citeF-12'></span>[[#cite-9|[9,10,12]]]. </li>

<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 display="inline">t+\Delta t</math>: velocities, pressure   and viscous stresses in the fluid and   displacements, stresses and strains in the solid. </li>

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

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

</ol>

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

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

Figure 3 shows a typical domain <math display="inline">V</math> with external boundaries <math display="inline">\Gamma _V</math> and <math display="inline">\Gamma _t</math> where the velocity and the surface tractions are prescribed, respectively. The domain <math display="inline">V</math> is formed by fluid (<math display="inline">V_F</math>) and solid (<math display="inline">V_S</math>) subdomains (i.e. <math display="inline">V=V_F \cup V_S</math>). Both subdomains interact at a common boundary <math display="inline">\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.

<div id='img-3'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[File:Draft_Samper_677613731_6179_fig3a.png|400px|]]
|[[File:Draft_Samper_677613731_4213_fig3b.png|400px|Split of the analysis domain V into fluid and solid   subdomains. Equality of surface tractions and kinematic variables at the   common interface]]
|- style="text-align: center; font-size: 75%;"
| colspan="2" | '''Figure 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 display="inline">{}^t S</math> and <math display="inline">{}^t F</math> the set of variables defining the kinematics and the stress-strain fields at the solid and fluid domains at time <math display="inline">t</math>, respectively, i.e.

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>{}^t S  :=  [{}^t {x}_s,{}^t {u}_s, {}^t {v}_s, {}^t {a}_s, {}^t {\boldsymbol \varepsilon }_s, {}^t {\boldsymbol \sigma }_s,\cdots ]^T</math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (1)
|-
| style="text-align: center;" | <math> {}^t F  :=  [{}^t {x}_F,{}^t {u}_F, {}^t {v}_F, {}^t {a}_F, {}^t \dot{\boldsymbol \varepsilon }_F, {}^t {\boldsymbol \sigma }_F,\cdots ]^T </math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (2)
|}
|}

where <math display="inline">{x}</math> is the nodal coordinate vector, '''u''', '''v''' and '''a''' are the vector of displacements, velocities and accelerations, respectively, <math display="inline">{\boldsymbol \varepsilon }, \dot{\boldsymbol \varepsilon }</math> and <math display="inline">{\boldsymbol \sigma }</math> are the strain vector, the strain-rate (or rate of deformation) vectors and the Cauchy stress vector, respectively and <math display="inline">F</math> and <math display="inline">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 3 is solved, in this work, using the following ''strongly coupled staggered scheme'':

<ol>

<li>We assume that the variables in the solid and fluid domains at time <math display="inline">t</math>   (<math display="inline">{}^t S</math> and <math display="inline">{}^t F</math>) are known. </li>
<li>Solve for the variables at the solid domain at time <math display="inline">t+\Delta t</math>   (<math display="inline">{}^{t+\Delta t}{S}</math>) under prescribed surface tractions at the   fluid-solid boundary <math display="inline">\Gamma _{FS}</math>. The boundary conditions at   the part of the external boundary intersecting the domain are   the standard ones in solid mechanics. </li>
<li>Solve for the variables at the fluid domain at time <math display="inline">t+\Delta t</math>   (<math display="inline">{}^{t+\Delta t}{F}</math>) under prescribed surface tractions at the   external boundary <math display="inline">\Gamma _{t}</math> and prescribed velocities at the external and   internal boundaries <math display="inline">\Gamma _{V}</math> and <math display="inline">\Gamma _{FS}</math>, respectively. </li>
<li>Iterate between 1 and 2 until convergence. </li>

</ol>

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

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>{M}_s {a}_s +{g}_s -  {f}_s = {0} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (3)
|}

where <math display="inline">{M}_s, {g}_s</math> and <math display="inline">{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.(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 <math display="inline">div</math>-stability condition <span id='citeF-5'></span><span id='citeF-34'></span>[[#cite-5|[5,34]]]. In our work we use a stabilized mixed FEM based on the Finite Calculus (FIC) approach which allows for a linear approximation for the velocity and pressure variables. Details of the FEM/FIC formulation  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-12'></span><span id='citeF-26'></span>[[#cite-12|[12,26]]]. Figure [[#img-4|4]]a shows the initial grid of four node rectangles discretizing the fluid domain and the solid walls.  Figures [[#img-4|4]]b and [[#img-4|4]]c show the mesh for the solution at two later times.

<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_677613731-Figure2a_b.png|600px|]]
|-
|[[File:Draft_Samper_677613731_7019_fig4.png|400px|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 domain at two different times.]]
|- style="text-align: center; font-size: 75%;"
| colspan="2" | '''Figure 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 domain 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-17'></span><span id='citeF-34'></span>[[#cite-17|[17,34]]].

'''Momentum'''

{| 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 } V_F </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (4)
|}

'''Mass balance'''

{| 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 } V_F </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (5)
|}

where

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

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

{| 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-17'></span><span id='citeF-18'></span><span id='citeF-19'></span><span id='citeF-21'></span>[[#cite-17|[17,18,19,21]]].

'''Momentum'''

{| 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}}{1\over 2} h_j{\partial r_{m_i} \over \partial x_j}=0 \qquad \hbox{in } V_F </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (10)
|}

'''Mass balance'''

{| 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}}{1\over 2} h_j {\partial r_d \over \partial x_j}=0 \qquad \hbox{in }V_F </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (11)
|}

The problem definition is completed with the following boundary conditions

{| 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}}{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)
|}

{| 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.(12) and (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. Recall that <math display="inline">\Gamma _v</math> includes both the external boundary and the internal boundary <math display="inline">\Gamma _{F_S}</math> (Figure 3).

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.(12) these lengths define the domain where equilibrium of boundary tractions is established. We note that at the discretized level, the <math display="inline">h_i's</math> become of the order of a typical element or grid dimension  <span id='citeF-17'></span><span id='citeF-18'></span><span id='citeF-19'></span>[[#cite-17|[17,18,19]]].

Eqs.(10)&#8211;(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-2'></span><span id='citeF-8'></span><span id='citeF-9'></span><span id='citeF-10'></span><span id='citeF-12'></span><span id='citeF-24'></span>[[#cite-2|[2,8,9,10,12,24]]].  Application of the FIC formulation to finite element and meshless analysis of fluid flow problems can be found in <span id='citeF-7'></span><span id='citeF-18'></span><span id='citeF-19'></span><span id='citeF-20'></span><span id='citeF-21'></span><span id='citeF-23'></span><span id='citeF-25'></span><span id='citeF-27'></span>[[#cite-7|[7,18,19,20,21,23,25,27]]].

===3.1 Transformation of the mass balance equation. Integral governing equations===

The underlined term in Eq.(11) can be expressed in terms of the momentum equations. The new expression for the mass balance equation is <span id='citeF-18'></span><span id='citeF-26'></span>[[#cite-18|[18,26]]]

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

{| 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)
|}

In our work we have taken the characteristic distances <math display="inline">h_i</math> to be constant within each element and equal to a typical element dimension computed as <math display="inline">h^e = [V^e]^m</math> where <math display="inline">V^e</math> is the element domain and <math display="inline">m=1/2</math> for 2D problems and <math display="inline">m=1/3</math> for 3D problems.

At this stage it is no longer necessary to retain the stabilization terms in the momentum equations and the traction boundary conditions (Eqs.(10) and (12)). These terms are  critical in Eulerian formulations to stabilize the numerical solution for high values of the convective terms <span id='citeF-18'></span><span id='citeF-21'></span><span id='citeF-27'></span><span id='citeF-28'></span>[[#cite-18|[18,21,27,28]]].

The weighted residual expression of the final form of the momentum and mass balance equations is written as

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\int _{V_F} \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)
|}

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\int _{V_F} 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.(17) and the deviatoric stresses and the pressure terms within <math display="inline">r_{m_i}</math> in Eq.(16) are integrated by parts to give

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\int _{V_F} \left[\delta v_i\rho  {\partial v_i \over \partial t}  + \delta \dot \varepsilon _{ij}(s_{ij}- \delta _{ij}p )\right]dV -  \int _{V_F} \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)
|}

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\int _{V_F} q {\partial v_i \over \partial x_i} dV + \int _{V_F} \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.(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.(19) has been neglected in this work. Retaining this term has been recently found to be advantageous for enhancing the satisfaction of the incompressibility condition in FEM predictor-corrector schemes for incompressible fluid flow analysis <span id='citeF-30'></span>[[#cite-30|[30]]].

===3.2 Pressure gradient projection===

The computation of the residual terms in Eq.(19) is simplified if we introduce the pressure gradient projections <math display="inline">\pi _i</math>, defined as

{| 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 <math display="inline">r_{m_i}</math> in  Eq.(19) in terms of the <math display="inline">\pi _i</math>  which then become additional variables. The system of integral equations is now 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:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\int _{V_F} \left[\delta v_i\rho {\partial v_i \over \partial t} + \delta \dot \varepsilon _{ij}(s_{ij}- \delta _{ij}p )\right]dV -  \int _{V_F} \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)
|}

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\int _{V_F} q {\partial v_i \over \partial x_i} dV + \int _{V_F} \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)
|}

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\int _{V_F}\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.(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

{| 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-34'></span>[[#cite-34|[34]]]. More details of the mesh discretization process and the choice of shape functions are given in Section 4.

Substituting the approximations (24) into Eqs.(21&#8211;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

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>{M}\dot{\bar{v}} + {K} \bar {v} - {G} \bar {p}- {   f}={0} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (25.a)
|}

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\displaystyle {G}^T \bar {v} + {L}\bar {p}+{Q}\bar {\boldsymbol    \pi }={0} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (25.b)
|}

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\displaystyle {Q}^T \bar {p} + \hat {M}\bar {\boldsymbol \pi }={0} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (25.c)
|}

The form of the element matrices and vectors in Eqs.(25) can be found in <span id='citeF-28'></span>[[#cite-28|[28]]].

===3.4 Fractional step algorithm for the fluid variables===

The starting point of the iterative algorithm are the variables at time <math display="inline">n</math> in the fluid domain (<math display="inline">{}^n F</math>). The sought variables are the variables at time <math display="inline">n+1</math> (<math display="inline">{}^{n+1} F</math>). For the sake of clarity we will skip the upper left index <math display="inline">n+1</math> for all variables, i.e.

{| 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}\,\, ;\,\, {}^{n+1} \bar {p}\equiv \bar  {p}\,\, ;\,\, {}^{n+1} \bar {\boldsymbol \pi }\equiv \bar {\boldsymbol \pi }\,\, ;\,\, {}^{n+1} \bar {x}\equiv \bar {x}\,\, ; \hbox{etc.} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (26)
|}

A simple  iterative algorithm is obtained by splitting the pressure from the momentum equations as follows

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\bar {v}^* = {}^n\bar {v} -\Delta t {M}^{-1}[{K} \bar {   v}^{j} -{G}{}^n \mathbf{p} - {f}] </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (27)
|}

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\bar{v}^{j+1}= \bar{v}^* + \Delta t {M}^{-1}{G}\delta  \bar{p} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (28)
|}

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

The value of <math display="inline">\bar {v}^{j+1}</math> from Eq.(29) is substituted now into Eq.(25b) to give

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>{G}^T\bar {v}^* + \Delta t {S}\delta \bar{p} + {L} \bar{p}^{j+1}+ {Q}\bar {\boldsymbol \pi }^{j}={0}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (29.a)
|}

where

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>{S} = {G}^T {M}^{-1}{G} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (29.b)
|}

Typically matrix '''S''' is computed using a diagonal matrix <math display="inline">{M} = {   M}_d</math>, where the subscript <math display="inline">d</math> denotes hereonward a diagonal matrix. We note   that the form of '''S''' in Eq.(29b) avoids the need for prescribing the   pressure at the boundary nodes.

An alternative is to approximate  '''S''' by a Laplacian matrix. This reduces considerably the bandwidth of '''S'''. 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. These problems are overcome however by retaining the residual term <math display="inline">r_{m_i}</math> in the boundary integral resulting from the integration by parts of Eq.(17) <span id='citeF-30'></span>[[#cite-30|[30]]]. In this work however the form of matrix '''S''' given by Eq.(29a) has been used.

A semi-implicit algorithm can  be derived as follows. For each iteration:<br/>

'''Step 1''' Compute <math display="inline">{v}^*</math>  from Eq.(27) with <math display="inline">{ M}={   M}_d</math>. For the first iteration

<math display="inline">(\bar {v}^1,\bar {p}^1,  \bar {\boldsymbol \pi }^1, \bar {x}^1 )\equiv ({}^n\bar {v}, {}^n\bar {p},  {}^n\bar  {\boldsymbol \pi }, {}^n\bar {x})</math><br/>

'''Step 2''' Compute <math display="inline">\delta  \bar {p}</math> and <math display="inline">{p}^{j+1}</math> from Eq.(29a) as

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\delta \bar {p} =-({L}+\Delta t  {S})^{-1} [{G}^T\bar{v}^* +{Q}\bar {\boldsymbol \pi }^{j}+ {L} \bar {p}^j] </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (30.a)
|}

The pressure <math display="inline">\bar {p}^{n+1,j}</math> is computed as

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\bar {p}^{j+1} = \bar{p}^j + \delta \bar{p}^j </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (30.b)
|}

'''Step 3''' Compute <math display="inline"> \bar{v}^{j+1}</math>  from Eq.(28) with <math display="inline">{M}={M}_d</math>

'''Step 4''' Compute <math display="inline">  \bar{\boldsymbol \pi }^{j+1}</math>  from Eq.(25c) as

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\bar{\boldsymbol \pi }^{j+1}=- \hat {M}_d^{-1} {Q}^T \bar {p}^{j+1} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (31)
|}

'''Step 5''' Update the coordinates of the mesh nodes as

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>{x}_i^{j+1} = {}^n {x}_i +\bar {v}_i^{j+1} \Delta t </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (32)
|}

'''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 1 for the next iteration with <math display="inline">j\leftarrow j+1 </math>.

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

In the examples presented in the paper the time increment size has been chosen as

{| 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 {v}\vert } </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (33)
|}

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.(25) are computed at the updated 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.

The boundary conditions are applied as follows. No condition is applied for the computation of the fractional velocities <math display="inline">{v}^*</math> in Eq.(27). The prescribed velocities at the boundary are applied when solving for <math display="inline">\bar{   v}^{j+1}</math> in step 3. As mentioned earlier there is no need for prescribing the   pressure at the boundary if the form of Eq.(29b) is chosen for '''S'''.

Box 1 shows a summary of the staggered scheme used for the solution for the variables in the solid and fluid domain at the updated configuration (<math display="inline">{}^{n+1}F,{}^{n+1}S</math>).

{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Samper_677613731-BoxI.png|451px|]]
|}
<div class="center" style="font-size: 75%;">
'''Box 1'''. Staggered scheme for the FSI problem (see also Figure 3)</div>


<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_677613731-polygon.png|400px|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 5:''' Generation of non standard meshes combining different polygons (in 2D) and polyhedra (in 3D) using the extended Delaunay technique.
|}

==4 GENERATION OF A NEW MESH==

One of the key points for the success of the PFEM  is the fast regeneration of a mesh at every time step on the basis of the position of the nodes in the space domain. Indeed, any fast meshing algorithm can be used for this purpose. In our work the mesh is generated at each time step using the so called extended Delaunay tesselation (EDT) presented in <span id='citeF-9'></span><span id='citeF-11'></span><span id='citeF-12'></span>[[#cite-9|[9,11,12]]]. 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-5|5]]). 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). In our work the simpler linear C<math display="inline">^\circ </math> interpolation has been chosen. Details of the mesh generation procedure and the derivation of the linear MFEM shape functions can be found in <span id='citeF-9'></span><span id='citeF-11'></span><span id='citeF-12'></span>[[#cite-9|[9,11,12]]].

Figure [[#img-6|6]] shows the evolution of the CPU time required for generating the mesh, for solving the system of equations and for assembling such a system in terms of the number of nodes. the numbers correspond to the solution of a 3D flow in an open channel with the PFEM. The figure shows the CPU time in seconds for each time step of the algorithm of Section 3.4. It is clearly seen that the CPU time required for meshing grows linearly with the number of nodes, as expected. Note also that the CPU time for solving the equations exceeds that required for meshing as the number of nodes increases. This situation has been found in all the problems solved with the PFEM. As a general rule for large 3D problems meshing consumes around 30% of the total CPU time for each time step, while the solution of the equations and the assembling of the system consume approximately 40% and 20% of the  CPU time for each time step, respectively. These figures prove that the generation of the mesh has an acceptable cost in the PFEM solution. An improvement of the mesh generation process will in any case help to reducing the computational cost.

<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_677613731-etiquetas.png|600px|3D flow problem solved with the PFEM. CPU time for meshing, assembling   and solving the system of equations at each time step in terms of the number of nodes]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 6:''' 3D flow problem solved with the PFEM. CPU time for meshing, assembling   and solving the system of equations at each time step in terms of the number of nodes
|}

==5 IDENTIFICATION OF BOUNDARY SURFACES==

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

In our work we use an extended Delaunay partition for recognizing 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. Values of <math display="inline">\alpha </math> ranging between 1.3 and 1.5 have been found to be optimal in all examples analyzed. This criterion is coincident with the Alpha Shape concept <span id='citeF-6'></span>[[#cite-6|[6]]]. Figure [[#img-7|7]] shows an 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 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. 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 (a particle) from the main analysis domain is again regained when the “flying” node falls down and a new boundary element is created by the Alpha Shape algorithm (Figures 2 and [[#img-7|7]]).

<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_677613731-Identification.png|600px|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.
|}

The boundary recognition method above described is also useful for detecting contact conditions between the fluid domain and a fixed boundary, as well as between different solids interacting with each other. The contact detection procedure is detailed in Section 6.

In order to show the quality of the boundary recognition approach, the following simple example has been performed. A square fluid domain of 0.25m<math display="inline">^2</math> is at a stationary position within a recipient (Figure [[#img-8|8]]). Then, as time evolves, the fluid falls down into the lower part of the recipient due to gravity effects. At the end of the process the total volume of the fluid within the recipient must be the same as that of the initial square domain. It must be noted that during the different time steps, the fluid has completely different free-surfaces including waves, breaking waves and fluid fragmentation zones.

The meshes used have average element sizes of 0.05m, 0.025m and 0.01m each which correspond to a total initial number of particles of 161, 552 and 3105 each. A value of <math display="inline">\alpha =1.4</math> for the Alpha-Shape method was used for the three analyses. Figure [[#img-8|8]] shows the initial position of the fluid domain and one of the three  meshes used for the analysis. Figure [[#img-9|9]] shows the fluid domain at different time steps.

<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_677613731-mesh.png|450px|Fluid domain following into a recipient. Initial position. Fine mesh of 3105 nodes (element size of 0.01m)]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 8:''' Fluid domain following into a recipient. Initial position. Fine mesh of 3105 nodes (element size of 0.01m)
|}

<div id='img-9a'></div>
<div id='img-9b'></div>
<div id='img-9c'></div>
<div id='img-9d'></div>
<div id='img-9e'></div>
<div id='img-9f'></div>
<div id='img-9g'></div>
<div id='img-9h'></div>
<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_677613731-time0_2.png|600px|t=0.2s]]
|[[Image:Draft_Samper_677613731-time0_4.png|600px|t=0.4s]]
|- style="text-align: center; font-size: 75%;"
| (a) <math>t=0.2</math>s
| (b) <math>t=0.4</math>s
|-
|[[Image:Draft_Samper_677613731-time0_6.png|481px|t=0.6s]]
|[[Image:Draft_Samper_677613731-time0_8.png|491px|t=0.8s]]
|- style="text-align: center; font-size: 75%;"
| (c) <math>t=0.6</math>s
| (d) <math>t=0.8</math>s
|-
|[[Image:Draft_Samper_677613731-time1.png|591px|t=1.0s]]
|[[Image:Draft_Samper_677613731-time2.png|600px|t=2.0s]]
|- style="text-align: center; font-size: 75%;"
| (e) <math>t=1.0</math>s
| (f) <math>t=2.0</math>s
|-
|[[Image:Draft_Samper_677613731-time3.png|600px|t=3.0s]]
|[[Image:Draft_Samper_677613731-time4.png|551px|t=4.0s]]
|- style="text-align: center; font-size: 75%;"
| (g) <math>t=3.0</math>s
| (h) <math>t=4.0</math>s
|- style="text-align: center; font-size: 75%;"
| colspan="2" | '''Figure 9:''' Positions of the fluid domain at different time steps.
|}

This simple example is interesting to show the quality of the boundary identification procedure. Another aim is to evaluate the volume variation from the incompressibility point of view, as well as the preservation of the total volume of the fluid due to possible errors in the boundary recognition using the Alpha-Shape method. Figure [[#img-10|10]] shows the total fluid volume during the different time steps for the three different meshes. The charge of volume is insignificant for the fine mesh and becomes larger but acceptable for the coarse meshes.

<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_677613731_3992_volume_history.jpg|500px|Total volume change as a function of time for different meshes.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 10:''' Total volume change as a function of time for different meshes.
|}

It must be noted that the main difference between the PFEM and the classical FEM is just the remeshing technique and the evaluation of the boundary position at each time step. The rest of the steps in the computation are coincident with those of the classical FEM.  This simple example shows that in spite of the permanent remeshing and the  evaluation of the boundary position via the Alpha-Shape method, the total fluid mass is preserved. We note however, that a good selection of the <math display="inline">\alpha </math> parameter is essential for the good behaviour of the boundary recognition process. Examples showing the accuracy of the PFEM for fixed boundary problems can be found in [2].

==6 TREATMENT OF CONTACT  CONDITIONS IN THE PFEM==

===6.1 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 common boundary <math display="inline">\Gamma _{FS}</math>, as mentioned above.

The  condition of prescribed velocities  at the fixed 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 the fluid particles and the fixed  boundaries is accounted for by the incompressibility condition which ''naturally prevents    the fluid nodes to penetrate into the solid boundaries'' (Figure [[#img-11|11]]). This simple way to treat the fluid-wall contact at mesh generation level is a distinct and attractive feature of the PFEM formulation.

<div id='img-11'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Samper_677613731-Figure6_1.png|600px|]]
|[[Image:Draft_Samper_677613731-Figure8.png|600px|Automatic treatment of contact conditions at the fluid-wall interface]]
|- style="text-align: center; font-size: 75%;"
| colspan="2" | '''Figure 11:''' Automatic treatment of contact conditions at the fluid-wall interface
|}

===6.2 Contact between solid-solid interfaces===

The contact between two solid interfaces is simply treated by introducing a layer of ''contact elements'' between the two interacting solid interfaces. This layer is ''automatically created during the mesh   generation step'' by prescribing a minimum distance (<math display="inline">h_c</math>) between two solid boundaries. If the distance exceeds the minimum value (<math display="inline">h_c</math>) 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 effects in the normal and tangential directions, respectively (Figure [[#img-12|12]]).

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

This contact algorithm can also be used effectively to model frictional contact conditions between rigid or elastic solids in standard structural mechanics applications. Figures [[#img-13|13]]&#8211;[[#img-16|16]] show examples of application of the contact algorithm to the  bumping of a ball falling in a container, the failure of a domino set, the failure of an arch formed by a collection of stone blocks under a seismic loading and the motion of five tetrapods as they  fall and slip over an inclined plane, respectively. The images in Figures [[#img-13|13]] and [[#img-17|17]] show explicitely  the  layer of contact elements which controls the accuracy of the contact algorithm.

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

<div id='img-13'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Samper_677613731-bola.png|451px|Bumping of a ball within a container. The layer of contact   elements is shown at each contact instant]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 13:''' Bumping of a ball within a container. The layer of contact   elements is shown at each contact instant
|}

<div id='img-14'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Samper_677613731-domino.png|451px|Failure of a domino set. The distance between the domino chips   shows the contact element layer]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 14:''' Failure of a domino set. The distance between the domino chips   shows the contact element layer
|}

<div id='img-15'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Samper_677613731-arch.png|451px|Failure of an arch   formed by  stone blocks under seismic loading]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 15:''' Failure of an arch   formed by  stone blocks under seismic loading
|}

<div id='img-16'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Samper_677613731-tetrapods-six.png|451px|Motion of five tetrapods on an inclined plane]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 16:''' Motion of five tetrapods on an inclined plane
|}

<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_677613731_8331_tetrapods5fall.jpg|400px|Detail of five tetrapods on an inclined plane. The layer of elements modeling the   frictional contact conditions is shown]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 17:''' Detail of five tetrapods on an inclined plane. The layer of elements modeling the   frictional contact conditions is shown
|}

==7 MODELING OF BED EROSION==

Prediction of bed erosion and sediment transport in open channel flows are  important tasks in many areas of river and environmental engineering. Bed erosion can lead to instabilities of the river basin slopes. It can also undermine the foundation of bridge piles thereby favouring structural failure. Modeling of bed erosion is also relevant for predicting the evolution of surface material dragged in earth dams in overspill situations. Bed erosion is one of the main causes of environmental damage in floods.

Bed erosion models are traditionally based on a relationship between the rate of erosion and the shear stress level <span id='citeF-14'></span><span id='citeF-33'></span>[[#cite-14|[14,33]]]. The effect of water velocity on soil erosion was studied in <span id='citeF-31'></span>[[#cite-31|[31]]]. In a recent work we have proposed an extension of the PFEM to model bed erosion <span id='citeF-29'></span>[[#cite-29|[29]]]. The erosion model is based on the frictional work  at the bed surface originated by the shear stresses in the fluid. The resulting erosion model resembles Archard law typically used for modeling abrasive wear in surfaces under frictional contact conditions <span id='citeF-1'></span><span id='citeF-22'></span>[[#cite-1|[1,22]]].

The algorithm for modeling the  erosion of soil/rock particles at the fluid bed is the following:

<ol>

<li>Compute at every point of the bed surface the resultant tangential   stress <math display="inline">\tau </math> induced by the fluid motion. In 3D problems <math display="inline">\tau =   (\tau _{s}^2 + \tau _{t})^2</math> where <math display="inline">\tau _s</math> and <math display="inline">\tau _t</math> are the tangential stresses   in the plane defined by the normal direction <math display="inline">{n}</math> at the bed   node. The value of <math display="inline">\tau </math> for 2D problems can be estimated as follows:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>

\tau _t =\mu \gamma _t</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (34.a)
|}</li>

with

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>

\gamma _t ={1\over 2}{\partial v_t\over \partial n} ={v_t^k\over 2h_k}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (34.b)
|}

where <math display="inline">v_t^k</math> is the modulus of the tangential velocity at the node <math display="inline">k</math>  and <math display="inline">h_k</math> is a prescribed distance along the normal of the bed node <math display="inline">k</math>. Typically <math display="inline">h_k</math> is of the order of magnitude of the smallest fluid element adjacent to node <math display="inline">k</math> (Figure [[#img-18|18]]).

<li>Compute the frictional work originated by the tangential stresses at the   bed surface as

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>

W_f =\int _\circ ^t \tau _t \gamma _t\, dt = \int _\circ ^t {\mu \over 4} \left({v_t^k\over h_k}\right)^2 dt </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (35)
|}</li>

Eq.(35) is integrated in time using a simple scheme as

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>

{}^n W_f ={}^{n-1} W_f + \tau _t \gamma _t\, \Delta t </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (36)
|}

<li>The onset of erosion at a bed point occurs when <math display="inline">{}^nW_f</math> exceeds a critical   threshold value <math display="inline">W_c</math> defined empirically   according to the specific properties of the bed material. </li>

<li>If <math display="inline">{}^nW_f > W_c</math> at a bed node, then the node is detached from the bed   region and it is allowed to move with the fluid flow, i.e. it becomes a fluid node. As a consequence, the   mass of the patch of bed elements surrounding the bed node vanishes in the   bed domain and it is transferred to the new   fluid node. This mass is subsequently transported with the fluid. Conservation of   mass of the bed particles within the fluid is guaranteed by changing the   density of the new fluid node so that the mass of the suspended sediment   traveling with the fluid  equals the  mass originally assigned to the bed   node. Recall that the mass assigned to a node is computed by multiplying the   node density by the tributary domain of the node. </li>

<li>Sediment deposition can be modeled by an inverse process to that described   in the previous step. Hence, a suspended node adjacent to the bed surface with a   velocity below a threshold value is assigned to the bed surface. This   automatically leads to the generation of   new bed elements adjacent to the boundary of the bed region. The original   mass of the bed region is recovered by adjusting the density of   the newly generated bed elements. </li>

</ol>

Figure [[#img-18|18]] shows an schematic view of the bed erosion algorithm proposed.

<div id='img-18'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Samper_677613731-Erosion.png|600px|Modeling of bed erosion by dragging of bed material]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 18:''' Modeling of bed erosion by dragging of bed material
|}

==8 FSI EXAMPLES==

The examples chosen show the applicability of the PFEM to solve problems involving large  motions of the free surface, fluid-multibody interactions and bed erosion.

===8.1 Rigid objects falling 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-19|19]] shows the penetration and evolution of a cube and a cylinder of rigid shape 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-20|20]]).

<div id='img-19'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Samper_677613731-Figure5.png|600px|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 19:''' 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-20'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Samper_677613731-Figure6.png|600px|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 20:''' 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-21|21]] shows an example of a wave breaking within a prismatic container including a vertical cylinder. Figure [[#img-22|22]] 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-21'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Samper_677613731-diap9.png|600px|]]
|-
|[[Image:Draft_Samper_677613731-diap10.png|600px|Evolution of a water column within a prismatic container including a   vertical cylinder.]]
|- style="text-align: center; font-size: 75%;"
| '''Figure 21:''' Evolution of a water column within a prismatic container including a   vertical cylinder.
|}

<div id='img-22'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Samper_677613731-diap17.png|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 22:''' 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-23|23]] shows the effect of a wave impacting on a rigid cube representing a vehicle. This situation is typical in  flooding and Tsunami situations. Note the layer of contact elements modeling the frictional contact conditions between the cube and the bottom surface.

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

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

Figure [[#img-24|24]] shows the 3D simulation of the impact of a wave generated in an experimental flume on a collection of  rigid rocks representing a breakwater. Details of the  water-rock interaction are shown in Figure [[#img-25|25]].

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

<div id='img-25'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Samper_677613731-Fig13.png|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 25:''' Detail of the impact of a wave on a breakwater. The arrows indicate   the water force on the rocks at different instants.
|}


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

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

Figure [[#img-30|30]] shows a detail of the complex three-dimensional interactions between the water particles and the tetrapods and between the tetrapods themselves, which can be easily modeled with the PFEM.

<div id='img-26'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Samper_677613731-Fig14.png|451px|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 26:''' 3D simulation of the impact of a wave on a collection of rocks in a breakwater.
|}

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

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

<div id='img-29'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[File:Draft_Samper_677613731_2572_tetrapods_slope_1.png]]
|style="padding-left:10px;"|[[File:Draft_Samper_677613731_4274_tetrapods_slope_2.png]]
|-
|[[File:Draft_Samper_677613731_6831_tetrapods_slope_3.png]]
|style="padding-left:10px;"|[[File:Draft_Samper_677613731_9144_tetrapods_slope_4.png]]
|-
|[[File:Draft_Samper_677613731_3839_tetrapods_slope_5.png]]
|style="padding-left:10px;"|[[File:Draft_Samper_677613731_1159_tetrapods_slope_6.png]]
|- style="text-align: center; font-size: 75%;"
| colspan="2" | '''Figure 29:''' Motion of ten tetrapods on a slope under an incident wave.
|}

<div id='img-30'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Samper_677613731-Figure24a.png|600px|]]
|[[Image:Draft_Samper_677613731-Figure24b.png|600px|Detail of the motion of ten tetrapods on a slope under an incident   wave. The figure shows the complex interactions between the water particles   and the tetrapods.]]
|- style="text-align: center; font-size: 75%;"
| colspan="2" | '''Figure 30:''' Detail of the motion of ten tetrapods on a slope under an incident   wave. The figure shows the complex interactions between the water particles   and the tetrapods.
|}

===8.5 Erosion of a 3D earth dam due to an overspill stream===

We present finally a simple, schematic, but very illustrative example showing the potential of the PFEM to model bed erosion in free surface flows.

The  example represents  the erosion of an earth dam under a water stream running over the dam top. A schematic geometry of the dam has been chosen to simplify the computations. Sediment deposition is not considered in the solution. The images of Figure [[#img-31|31]]  show the progressive erosion of the dam  until the whole dam is dragged out by the fluid flow.

Other applications of the PFEM to bed erosion problems can be found in <span id='citeF-29'></span>[[#cite-29|[29]]].

<div id='img-31'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Samper_677613731-erosion_3D.png|451px|Erosion of a 3D earth dam due to an overspill stream.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 31:''' Erosion of a 3D earth dam due to an overspill stream.
|}

==9 CONCLUSIONS==

The particle finite element method (PFEM) is ideal to treat problems involving fluids with free surfaces and submerged or floating structures and bodies within a unified Lagrangian finite element framework. Problems such as fluid-structure interaction, large motion of  fluid or solid particles, surface waves, water splashing, separation of water drops, frictional  contact situations between fluid-solid and solid-solid interfaces, bed erosion, 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 an updated Lagrangian formulation and 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 identification of the boundary nodes using an Alpha-Shape type technique and the simple algorithm to treat frictional contact conditions at fluid-solid and solid-solid interfaces via mesh generation. The examples presented have shown the great potential of the PFEM for solving a wide class of practical FSI problems in engineering. Examples of validation of the PFEM results with data from experimental tests are reported in <span id='citeF-15'></span>[[#cite-15|[15]]].

==ACKNOWLEDGEMENTS==

Thanks are given to Mrs. M. de Mier  for many useful suggestions. This research was partially supported by project SEDUREC of the Consolider Programme of the Ministerio de Educación y Ciencia of Spain.

==REFERENCES==

<div id="cite-1"></div>
'''[[#citeF-1|[1]]]'''   J.F. Archard, Contact and rubbing of flat surfaces,      J. Appl. Phys.  24(8) (1953) 981&#8211;988.

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

<div id="cite-3"></div>
'''[3]'''   R. Codina, O.C. Zienkiewicz, 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  (2002)   18(2) (2002) 99&#8211;112.

<div id="cite-4"></div>
'''[[#citeF-4|[4]]]'''  F. Del Pin, S.R. Idelsohn, E. Oñate, R. Aubry, The ALE/Lagrangian particle finite element method: A new approach to computation of free-surface flows and fluid-object interactions. Computers & Fluids  36  (2007) 27&#8211;38.

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

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

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

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

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

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

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

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

<div id="cite-13"></div>
'''[[#citeF-13|[13]]]'''  S.R. Idelsohn, E. Oñate, F. Del Pin, N. Calvo, Fluid-structure   interaction using the particle finite element   method, Comput. Meth. Appl. Mech. Engng. 195 (2006) 2100-2113.

<div id="cite-14"></div>
'''[[#citeF-14|[14]]]'''  A. Kovacs, G. Parker,  A new vectorial bedload formulation and   its application to the time evolution of straight river channels, J. Fluid   Mech. 267 (1994) 153&#8211;183.

<div id="cite-15"></div>
'''[[#citeF-15|[15]]]'''  L. Larese,  R. Rossi, E. Oñate, S.R. Idelsohn,  Validation of the Particle Finite Element Method (PFEM) for simulation of free surface flows, submitted to Engineering Computations, March (2007).

<div id="cite-16"></div>
'''[[#citeF-16|[16]]]'''  R. Ohayon, Fluid-structure interaction problem, in: E. Stein, R. de Borst, T.J.R. Hugues (Eds.), Enciclopedia of Computatinal Mechanics, Vol. 2, (J. Wiley, 2004) 683&#8211;694.

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

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

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

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

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

<div id="cite-22"></div>
'''[[#citeF-22|[22]]]'''  E. Oñate, J. Rojek, Combination of discrete element and   finite element method for dynamic analysis of geomechanic   problems,  Comput. Meth. Appl. Mech. Engrg. 193 (2004) 3087&#8211;3128.

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

<div id="cite-24"></div>
'''[[#citeF-24|[24]]]'''  E. Oñate, S.R. Idelsohn, F. Del Pin,  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 (2003).

<div id="cite-25"></div>
'''[[#citeF-25|[25]]]'''  E. Oñate, J. García, S.R. Idelsohn, Ship   hydrodynamics, in: E. Stein,   R. de Borst, T.J.R. Hughes (Eds), Encyclopedia of Computational Mechanics, Vol 3,   (J. Wiley 2004a) 579&#8211;610.

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

<div id="cite-27"></div>
'''[[#citeF-27|[27]]]'''  E. Oñate, A. Valls, J. García,  FIC/FEM formulation   with matrix stabilizing terms for incompressible flows at low and high   Reynold's numbers, Computational Mechanics 38 (4-5) (2006a) 440-455.

<div id="cite-28"></div>
'''[[#citeF-28|[28]]]'''  E. Oñate, J. García,  S.R. Idelsohn, F. Del Pin,  FIC   formulations for finite element analysis of incompressible flows. Eulerian,   ALE and Lagrangian approaches,  Comput. Meth. Appl. Mech. Engng.   195 (23-24) (2006b) 3001-3037.

<div id="cite-29"></div>
'''[[#citeF-29|[29]]]'''  E. Oñate, M.A. Celigueta, S.R. Idelsohn, Modeling bed erosion in   free surface flows by the Particle Finite Element Method, Acta   Geotechnia 1 (4) (2006c) 237-252.

<div id="cite-30"></div>
'''[[#citeF-30|[30]]]'''  E. Oñate, S.R. Idelsohn, R. Rossi, Enhanced FIC-FEM formulation for incompressible flows, Research Report, CIMNE Barcelona, March 2007.

<div id="cite-31"></div>
'''[[#citeF-31|[31]]]'''  D.B. Parker, T.G. Michel, J.L. Smith,  Compaction and water   velocity effects on soil erosion in shallow flow,  J. Irrigation and Drainage   Engineering 121 (1995) 170&#8211;178.

<div id="cite-32"></div>
'''[[#citeF-32|[32]]]'''  T.E. Tezduyar, Finite element method for fluid dynamics with moving boundaries and interface, in: E. Stein, R. de Borst, T.J.R. Hugues (Eds.), Enciclopedia of Computatinal Mechanics, Vol. 3, (J. Wiley, 2004) 545&#8211;578.

<div id="cite-33"></div>
'''[[#citeF-33|[33]]]'''  C.F. Wan, R. Fell, Investigation of erosion of soils in   embankment dams, J. Geotechnical and Geoenvironmental Engineering   130 (2004) 373&#8211;380.

<div id="cite-34"></div>
'''[[#citeF-34|[34]]]'''  O.C. Zienkiewicz, R.L. Taylor, P. Nithiarasu, The finite element   method for fluid dynamics,   Elsevier, (2006).

<div id="cite-35"></div>
'''[35]'''  O.C. Zienkiewicz, R.L. Taylor,   The finite element method for   solid and structural mechanics,  Elsevier, (2005).