Unified Lagrangian formulation for elastic solids and incompressible fluids. application to fluid-structure interaction problems via the PFEM



==Abstract==

We present a general Lagrangian formulation for treating elastic   solids and quasi/fully incompressible fluids in a unified form. The   formulation allows to treat solid and fluid subdomains in a unified manner in   fluid-structure interaction (FSI) situations. In our work the FSI problem is   solved via the Particle Finite Element Method (PFEM). The PFEM is an   effective technique for modeling complex interactions between floating and   submerged bodies and free surface flows, accounting for splashing of waves,   large motions of the bodies and frictional contact conditions. Applications of   the unified Lagrangian formulation to a number of FSI problems are given.

'''Keywords''': Lagrangian formulation, fluid-structure, particle finite element method

==1 INTRODUCTION==

Fluid-structure interaction (FSI) problems are typically solved with a staggered scheme [23] by which the relevant variables at the fluid and solid subdomains are separately and sequentially (and iteratively) solved using as boundary conditions at the common fluid-solid boundaries the velocities (for the fluid problem) and the surface tractions (for the solid problems). The staggered scheme is ideal for using existing finite element codes, initially developed for fluid dynamics and solid mechanics problems, and the computing effect is mainly focused on the interfacing of the relevant data between the common fluid and solid boundaries. Indeed the FSI problem typically involves the motion of mesh nodes in both the fluid and solid domains. As a consequence, an arbitrary Lagrangian-Eulerian (ALE) formulation [24] is used to model the governing equations for the fluid, while a standard Lagrangian formulation is used for the equations of the solid part.

In our work we propose a different route for solving FSI problems. Our goal is to solve the equations for both the fluid and solid domains using a  unified lagrangian formulation. This basically means that the analysis domain, containing both fluid and solid subdomains which interact with each other, is seen as a single continuum domain with different material properties assigned to each of the interacting subdomains (i.e. the fluid and solid regions). This allows to make no distinction between fluids and solids for the numerical solution and a single computer code can be used for solving the FSI problem.

The governing equations for the fluid and solid domains (in a lagrangian frame of reference) are discretized and solved with the Particle Finite Element Method (PFEM). The PFEM  treats the mesh nodes in the fluid and solid domains as moving material points (henceforth called 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 motion of the  mesh discretizing the total domain (including both the fluid and solid parts) is followed during the transient solution. 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  [1,7,8,10,18,19,20]

In summary, the key ingredients of the unified formulation presented in this paper are: a) the definition of a unified constitutive equation for the fluid and solid materials, b) the use of a Lagrangian description to model the kinematic of both fluid and solid domains, and c) the use of the Particle Finite Element Method (PFEM) for redefinition of the domain boundaries and the treatment of frictional contact conditions.

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 [7,9]. Furthermore, this special polyhedral finite element needs special shape functions. In this paper, meshless finite element (MFEM) shape functions have been used [7]. Another possibility is the use of a standard mesh generator with sliver elimination [25]. Nevertheless, standard mesh generator are som times too expensive in computational cost. An efficient mesh generator as such presented in ref. [9] is a key issue in a Lagrangian formulation.

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 elastic solid using a hypo-elastic approximation are presented.  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 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 [[#img-1|1]]). </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 process including separation and   re-entering of nodes. This allows to model splashing of waves.   The Alpha Shape method [4]  is used for the boundary definition  (see Section 6). </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 [9]. </li>

<li>Solve the  Lagrangian equations of motion for  the fluid and the   solid domains using the unified formulation proposed in this work. 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. Indeed this step requires   an iterative scheme as large motion of both the fluid and solid domain may   occur during the time step. Also material non linearities may appear in the   solid although this situation is not considered in this work. </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. </li>

</ol>

<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_577315608_7676_Fig2.jpg|500px|Sequence of steps to update a “cloud” of nodes from time n   (t=tₙ)  to   time n+1 (t=tₙ+∆t)  ]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 1:''' Sequence of steps to update a “cloud” of nodes from time <math>n</math>   (<math>t=t_n</math>)  to   time <math>n+1</math> (<math>t=t_n +\Delta t</math>)  
|}

==3 CONSTITUTIVE EQUATIONS TO BE USED FOR BOTH FLUID AND SOLID MATERIALS==

===3.1 Constitutive equations for hypo-elastic solids===

Let a material with a hypo- elastic constitutive equation like:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>L(\tau _{ij}) =\lambda _s \delta _{ij} d_{ll}+2\mu _s d_{ij} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (1)
|}

where <math display="inline">\tau _{ij} = J\sigma _{ij}</math> is the Kirchhoff stress tensor, <math display="inline">d_{ij}=\displaystyle \frac{1}{2} \left(\displaystyle {\partial V_i \over \partial x_j}+{\partial V_j \over \partial x_i}\right)</math>  is the rate of deformation tensor, <math display="inline">V_i</math> the velocity along the <math display="inline">i</math>th axis, <math display="inline">\sigma _{ij}</math>  the Cauchy stress tensor, <math display="inline">\lambda _s</math> and <math display="inline">\mu _s</math> the Lamé  parameters, <math display="inline">J = {det} (F_{ij})</math>  the Jacobian, being <math display="inline">F_{ij}=\displaystyle {\partial x_i \over \partial X_j}={\partial u_i \over \partial X_j}+\delta _{ij}</math> the deformation gradient tensor, <math display="inline">u_i</math> the <math display="inline">i</math>th displacement component and <math display="inline">L(\tau _{ij})=F\dot S F^T</math>, the time Lie derivative , with <math display="inline">S=F^{-1} \tau F^{-T} =F^{-1}\sigma F^{-T} J</math>, the second Piola-Kirchhoff stress tensor.

Dividing the strain and the stress in their deviatoric (<math display="inline">^\prime </math>) and the volumetric parts

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>d_{ij}= d_{ij}' + {d_{ll}\over 3} \delta _{ij} \quad ;\quad  \tau _{ij}=\tau _{ij}' + {\tau _{ll}\over 3} \delta _{ij} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (2)
|}

Then

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>L(\tau _{ij}) =\left(\lambda _s + {2\mu _s \over 3}\right)\delta _{ij} d_{ll}+2\mu _s d_{ij}' </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (3)
|}

This may be split as

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>L(\tau _{ij}')  = 2\mu _s d_{ij}'</math>
|-
| style="text-align: center;" | <math> L\left({\tau _{ll}\over 3} \right)= \left(\lambda _s + {2\mu _s \over 3}\right)d_{ll} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (4)
|}

The volumetric strain rate and the pressure will be written as

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\varepsilon _v =tr (d_{ij}) = d_{ll}\quad \hbox{and}\quad p =tr \left(-{\sigma _{ij}\over 3}  \right)= - {\sigma _{ll}\over 3} = - {\tau _{ll}\over 3J} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (5)
|}

Approaching the derivative in (4) by a finite time step

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\left(F {\Delta S\over \Delta t} F^T\right)^\prime  =2\mu _s d'   \quad \hbox{and}\quad  tr \left(F   {\Delta S\over 3\Delta t} F^T\right)= \left( \lambda _s + {2\mu _s \over      3}\right)\varepsilon _v </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (6)
|}

which may be written as a function of the Cauchy stress:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\sigma _{ij}^{'(n+1)} =\hat \sigma _{ij}^{'(n)}+ {2\Delta t\mu _s\over J}{ d'}_ {ij}\quad \hbox{and}\quad p^{n+1}-\hat p^n = - {\Delta t\over J} \left( \lambda _s + {2\mu _s \over    3}\right)\varepsilon _v </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (7)
|}

where:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" |<math>\begin{array}{l}
\hat\sigma^{(n)}= \displaystyle \frac{F S^n F^T}{J} \quad ; \quad  \hat\sigma^{'(n)}=\displaystyle \frac{( F   S^n F^T)^\prime}{J}\\[.3cm] 
\displaystyle \hat p^{(n)} = \displaystyle \frac{- tr \left[F \left(S^n\over 3  \right)F^T   \right]}{J} \quad ;\quad   p^{n+1}=\displaystyle \frac{-tr \left[F \left(\frac{S^{n+1}}{3}  \right)F^T \right]}{J}\end{array}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (8)
|}

In the previous equations <math display="inline">\hat \sigma ^{'(n)}</math> and <math display="inline">\hat p^{(n)}</math>  represent the deviatoric stress and the pressure at the beginning of the time step <math display="inline">(n)</math>, but referred to the final time step configuration <math display="inline">(n+1)</math>.

Finally for the hypo-elastic material the constitutive relations may be written in terms of the following three expressions:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\sigma _{ij}^{n+1} = \hat \sigma _{ij}^{'(n)}+ {2\Delta t\mu _s\over J}{ d'}_ {ij} - p^{n+1}\delta _{ij}</math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (9)
|-
| style="text-align: center;" | <math> \sigma _{ij}^{n+1} = \hat \sigma _{ij}^{(n)}+ {2\Delta t\mu _s\over J}{ d'}_{ij}+{\Delta t\over J} \left( \lambda _s + {2\mu _s \over    3}\right) \varepsilon _v\delta _{ij}</math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (10)
|-
| style="text-align: center;" | <math> \sigma _{ij}^{n+1} = \hat \sigma _{ij}^{(n)}+ {2\Delta t\mu _s\over J} d_{ij}+{\Delta t\over J} \lambda _s \varepsilon _v\delta _{ij} </math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (11)
|}
|}

===3.2 Constitutive relations for incompressible and quasi-incompressible fluids===

For Newtonian fluid flows the standard constitutive relations are:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\sigma _{ij}^{n+1} = 2\mu {_L d'}_{ij}- p^{n+1}\delta _{ij} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (12)
|}

where <math display="inline">\mu _L</math> is the viscosity.

For quasi-incompressible flows, the volumetric train rate may be written as a strain function of the sound speed by:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>{Dp\over Dt}\simeq {\Delta p\over \Delta t} = -C^2 \rho \varepsilon _v \Rightarrow p^{n+1} =p^n - \Delta t C^2 \rho \varepsilon _v =p^n - \Delta t \kappa _l \varepsilon _v </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (13)
|}

where <math display="inline">C</math> is the sound speed and <math display="inline">\kappa _l  = C^2 \rho </math>.

Then, the Cauchy stress tensor may be written in function of the volumetric strain rate by:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\sigma _{ij}^{n+1} = 2\mu {_L d'}_{ij}- p^n \delta _{ij} + \Delta t \kappa _l \varepsilon _v \delta _{ij} = 2\mu _L d_{ij} - p^n \delta _{ij} + \left(\Delta t \kappa _l - {2\mu _l\over 3} \right)\varepsilon _v \delta _{ij} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (14)
|}

From (12&#8211;14) Newtonian constitutive relation for incompressible or near incompressible flow may be written in one of the  following three manners:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\sigma _{ij}^{n+1} = 2\mu {_L d'}_{ij}- p^{n+1}\delta _{ij}\quad \hbox{ (for incompressible)}</math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (15)
|-
| style="text-align: center;" | <math> \sigma _{ij}^{n+1} = 2\mu {_L d'}_{ij}- p^n \delta _{ij} + \Delta t \kappa _l \varepsilon _v \delta _{ij}\,\, \hbox{ (for quasi-incompressible)}</math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (16)
|-
| style="text-align: center;" | <math> \sigma _{ij}^{n+1} = 2\mu _L d_{ij}- p^n \delta _{ij}+\left(\Delta t \kappa _l - {2\mu _l\over 3} \right)\varepsilon _v \delta _{ij}\,\, \hbox{ (for quasi-incompressible)}\qquad  </math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (17)
|}
|}

===3.3 Unique constitutive equation for fluids and solids===

In the following, only a unified constitutive equation will be used for both the elastic solid and the incompressible or nearly incompressible fluid:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\sigma _{ij}^{n+1} = \hat \sigma _{ij}^{(n)}+ 2\mu {d'}_{ij}- \Delta  p^{n+1}\delta _{ij}\quad \hbox{ (for incompressible)}</math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (18)
|-
| style="text-align: center;" | <math> \sigma _{ij}^{n+1} = \hat \sigma _{ij}^{(n)}+ 2\mu {d'}_{ij}+\Delta t \kappa \varepsilon _v \delta _{ij}\quad \hbox{ (for quasi-incompressible)}</math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (19)
|-
| style="text-align: center;" | <math> \sigma _{ij}^{n+1} = \hat \sigma _{ij}^{(n)}+ 2\mu d_{ij} + \lambda \varepsilon _v \delta _{ij}\quad \hbox{ (for quasi-incompressible)} </math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (20)
|}
|}

with the  definitions for <math display="inline">\hat \sigma _{ij}^{(n)}</math>, <math display="inline">\mu </math>, <math display="inline">\lambda </math>  and <math display="inline">\kappa </math> given in Table 1.


{|  class="floating_tableSCP wikitable" style="text-align: center; margin: 1em auto;min-width:50%;"
|+ style="font-size: 75%;" |Table. 1 Definition of physical parameters
|- style="border-top: 2px solid;"
| style="border-left: 2px solid;border-right: 2px solid;" |  '''Solid''' 
| style="border-left: 2px solid;border-right: 2px solid;" | '''Fluid''' 
| style="border-left: 2px solid;border-right: 2px solid;" | '''Unique'''

|- style="border-top: 2px solid;"
| style="border-left: 2px solid;border-right: 2px solid;" | <math>\hat \sigma _{ij}^{(n)}</math>
| style="border-left: 2px solid;border-right: 2px solid;" | <math>-p^n \delta _{ij}</math>
| style="border-left: 2px solid;border-right: 2px solid;" | <math>\hat \sigma _{ij}^{(n)}</math>


|- style="border-top: 2px solid;"
| style="border-left: 2px solid;border-right: 2px solid;" | <math>\displaystyle {\Delta t\mu _s\over J}</math>
| style="border-left: 2px solid;border-right: 2px solid;" | <math>\mu _L</math>
| style="border-left: 2px solid;border-right: 2px solid;" | <math>\mu </math>


|- style="border-top: 2px solid;"
| style="border-left: 2px solid;border-right: 2px solid;" | <math>\displaystyle {\Delta t\lambda _s\over J}</math>
| style="border-left: 2px solid;border-right: 2px solid;" | <math> \Delta t \kappa _l - \displaystyle {2\mu _l\over 3}</math>
| style="border-left: 2px solid;border-right: 2px solid;" | <math>\lambda </math>


|- style="border-top: 2px solid;"
| style="border-left: 2px solid;border-right: 2px solid;" | <math>\displaystyle {1\over J} \left( \lambda _s + {2\mu _s \over    3}\right)</math>
| style="border-left: 2px solid;border-right: 2px solid;" | <math>\kappa _l</math>
| style="border-left: 2px solid;border-right: 2px solid;" | <math>\kappa </math>


|}

Eq.(18) will be used only in such cases in which all the domain or a part of the domain is totally incompressible, while Eq.(19)  or (20) will be used in such cases in which all the domain may be considered as compressible or quasi-incompressible.

Then, depending of the material the following definitions will be used:

<br/>

''a) For the fluid''

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\mu = \mu _l</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>\hat \sigma _{ij}^{'(n)}\equiv 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>\hat \sigma _{ij}^{(n)} = -p^n \delta _{ij} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (23)
|}

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\kappa = \kappa _l =C^2 \rho _l</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (24)
|}

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\lambda = \lambda _l = \Delta t\kappa \varepsilon _v </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (25)
|}

Totally incompressible flows means <math display="inline">C=\infty </math> and then <math display="inline">\kappa = \kappa _l =\infty </math>.<br/>

''b) For the solid part:''

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\mu = {\Delta t\over J}\mu _s = {\Delta t\over J}{E\over 2(1+\nu )} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (26)
|}

where <math display="inline">E</math> is the Young modulus, <math display="inline">\nu </math> the  Poison coefficient and <math display="inline">\lambda _s</math> and <math display="inline">\mu _s</math> the Lamé parameters.

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\lambda = \lambda _s = {\nu E\over (1+\nu )(1-2\nu )} </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>\kappa = \kappa _s = {1\over J} \left( \lambda _s + {2\mu _s \over    3}\right) </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (28)
|}

===3.4 Momentum conservation equations===

The standard infinitesimal momentum conservation equation may be written in a Lagrangian frame as:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\rho a_i = \rho {DV_i\over Dt} = {\partial \sigma _{ij} \over \partial x_j}+\rho f_i </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (29)
|}

The equations are completed with the pressure-volumetric strain equations

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\varepsilon _v  = -{\Delta p\over \kappa \Delta t} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (30)
|}

and the boundary conditions:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\sigma _{ni} =\bar \sigma _{ni}\quad \hbox{on } \Gamma _\sigma  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (31)
|}

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>V_i = \bar V_i \quad \hbox{on } \Gamma _V </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (32)
|}

It must be noted that the term <math display="inline">{DV_i\over Dt}</math>, represents the time material derivative (lagrangian) of the velocity <math display="inline"> {DV_i\over Dt}={V_i^{n+1}-V_i^n\over \Delta t}</math> where <math display="inline">V_i^{n+1}</math> represents the velocity at time <math display="inline">(n+1)</math> in the position <math display="inline">x^{n+1}_i</math>. The convective term, normally included in the fluid mechanics equations, are not explicitly present in the lagrangian formulation.

A weighted residual method will be used to approximate above equations:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\int _V W_i \left(\rho {DV_i\over Dt}-{\partial \sigma _{ij} \over \partial x_j}-\rho f_i\right) dV + \int _{\Gamma _\sigma } W_i (\sigma _{ni}-\bar \sigma _{ni} ) d\Gamma =0</math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (33)
|-
| style="text-align: center;" | <math>  \int _V W_p \left(-\varepsilon _v - {\Delta p\over \kappa \Delta t}\right) dV=0 </math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (34)
|}
|}

In weak form:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\int _V \left(W_i \rho {DV_i\over Dt}+ {\partial W_i \over \partial x_j}\sigma _{ij} -W_i \rho    f_i\right) dV - \int _{\Gamma _\sigma } W_i \bar \sigma _{ni} d\Gamma =0</math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (35)
|-
| style="text-align: center;" | <math>  \int _V W_p \left(-\varepsilon _v - {\Delta p\over \kappa \Delta t}\right) dV=0 </math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (36)
|}
|}

==4 TEMPORAL INTEGRATION OF GOVERNING EQUATIONS==

Let

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>a_i^{n+1/2} = {V_i^{n+1}-V_i^n\over \Delta t}= {a_i^{n+1}+a_i^n\over    2}\Rightarrow a_i^{n+1} = 2\left({V_i^{n+1}-V_i^n\over \Delta t}  \right)- a_i^n </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (37)
|}

Then, at time <math display="inline">(n+1)</math> the weak form of the weighted residual equation becomes:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\int _V \left(W_i 2\rho \left({V_i^{n+1}-V_i^n\over \Delta t}\right)- W_i \rho   a_i^n + {\partial W_i \over \partial x_j}\sigma _{ij}^{n+1} -W_i \rho f_i^{n+1}\right)dV- </math>
|-
| style="text-align: center;" | <math> \qquad - \int _{\Gamma _\sigma } W_i \bar \sigma _{ni} d\Gamma =0</math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (38)
|-
| style="text-align: center;" | <math> \int _V W_p \left(-\varepsilon _v^{n+1} - {\Delta p\over \kappa \Delta t}\right)dV=0 </math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (39)
|}
|}

In the following, the upper index <math display="inline">n+1</math> will be dropped resulting in:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\int _V \left(W_i 2\rho \left({V_i-V_i^n\over \Delta t}\right)- W_i \rho   a_i^n + {\partial W_i \over \partial x_j}\sigma _{ij} -W_i \rho f_i\right)dV -  \!\int _{\Gamma _\sigma } W_i \bar \sigma _{ni} d\Gamma =0 </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (40)
|}

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\int _V W_p \left(-\varepsilon _v - {\Delta p\over \kappa \Delta t}\right)dV=0 </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (41)
|}

===4.1 Case in which all the domain may be considered as  quasi-incompressible===

In this case, the constitutive equations  are the described in Eqs.(19) or (20). The pressure is not a state variable:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\int _V \left(W_i 2 \rho \left({V_i-V_i^n\over \Delta t}\right)- W_i \rho   a_i^n + {\partial W_i \over \partial x_j} (\hat \sigma _{ij}^{(n)} +2\mu d_{ij} + \lambda  \varepsilon _v \delta _{ij})-W_i \rho f_i \right)dV</math>
|-
| style="text-align: center;" | <math>  - \int _{\Gamma _\sigma } W_i\bar \sigma _{ni} d\Gamma =0 </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (42)
|}

or

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\int _V \left[W_i (2\rho )V_i + \Delta t {\partial W_i \over \partial x_j}2\mu d_{ij} + \Delta t   {\partial W_i \over \partial x_i} \lambda \varepsilon _v \right]dV=</math>
|-
| style="text-align: center;" | <math> \int _V \left[W_i (2\rho )V_i^n + \Delta t W_i \rho a_i^n + \Delta t W_i \rho    f_i -\Delta t {\partial W_i \over \partial x_j}\hat \sigma _{ij}^{(n)}\right]dV + \Delta t \int _{\Gamma _\sigma } \bar \sigma _{ni} d\Gamma  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (43)
|}

===''Spatial discretization''===

Each of the velocity components <math display="inline">V_i</math> is interpolated using the MFEM shape functions<math display="inline">^7</math>  as:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>V_i = N^T Q_i </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (44)
|}

where <math display="inline">N^T</math>  are the MFEM shape functions and <math display="inline">Q_i</math> a vector containing the nodal values of the velocity components.

For Galerkin residual approximations, the arbitrary weighting functions <math display="inline">W_i</math>  are:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>[W_1,W_2,W_3]= \left[\begin{matrix}N &0&0\\ 0 &N&0\\ 0&0&N \end{matrix}\right] </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (45)
|}

The equation to be solved in matrix form becomes:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>(M_{ij} +\Delta t K_{ij}) Q_j = G_i^n </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (46)
|}

with

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\begin{array}{l}K_{ij} = K_{ij}^1 + K_{ij}^2 + K_{ij}^3\\[.2cm] M_{ii} = \int _V  N (2\rho )N^T dV \quad ;\quad K_{ii}^1 = \int _V \displaystyle {\partial N \over \partial x_j} (\mu ) \displaystyle {\partial N^T \over \partial x_j}dV\\[.2cm] M_{ij} = K_{ij}^1 =0 \qquad \forall i\not = j\\[.2cm] K_{ij}^2 = \int _V \displaystyle {\partial N \over \partial x_j} (\mu ) \displaystyle {\partial N^T \over \partial x_i}dV\\[.2cm] K_{ij}^3 = \int _V \displaystyle {\partial N \over \partial x_i} (\lambda ) \displaystyle {\partial N^T \over \partial x_j}dV\\[.2cm]  G_i^n = \int _V\left[N (2\rho ) V_i^n + N (\Delta t\rho )a^n - \Delta t \displaystyle     {\partial N \over \partial x_j}\hat \sigma ^n_{ij} + N \Delta t\rho f_i \right]dV +  \int _{\Gamma _\sigma } N \Delta t \bar \sigma _{ni} d\Gamma   \end{array} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (47)
|}

===4.2 Case in which all the domain or a part of it must be considered as incompressible===

In this case, the only possibility is to use a mixed formulation Eq.(18) using the velocity and the pressure as unknown variables. The weighted residual equation remains:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\int _V \left(W_i 2\rho \left({V_i-V_i^n\over \Delta t}\right)- W_i \rho   a_i^n + {\partial W_i \over \partial x_j} (\hat \sigma _{ij}^{(n)} +2\mu {d'}_{ij}-\Delta p  \delta _{ij} ) - W_i \rho f_i\right)dV -</math>
|-
| style="text-align: center;" | <math> \qquad - \int _{\Gamma _\sigma } \bar \sigma _{ni}  d\Gamma =0</math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (48)
|-
| style="text-align: center;" | <math> \int _V W_p \left(- \varepsilon _v - {\Delta p\over \kappa \Delta t}\right)dV =0 </math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (49)
|}
|}

or also

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\int _V \left[W_i (2\rho )V_i + \Delta t {\partial W_i \over \partial x_j}2\mu {d'}_{ij} - \Delta t   {\partial W_i \over \partial x_i} \Delta p \right]dV</math>
|-
| style="text-align: center;" | <math> =\!\!\int _V \left[W_i (2\rho )V_i^n + \Delta t W_i \rho a_i^n - \Delta t   {\partial W_i \over \partial x_j} \hat \sigma _{ij}^{(n)} + \Delta t    W_i \rho    f_i\right]dV + \Delta t \int _{\Gamma _\sigma } \bar \sigma _{ni} d\Gamma \quad \qquad    </math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (50)
|-
| style="text-align: center;" | <math>  \int _V \left[\Delta t W_p {\partial V_i \over \partial x_i}+W_p {\Delta p\over    \kappa } \right]dV=0 </math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (51)
|}
|}

Taking into account the definition of the deviatoric strain rate:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\int _V \left[W_i (2\rho )V_i + \Delta t {\partial W_i \over \partial x_j}\mu    {\partial V_i \over \partial x_i}+\Delta {\partial W_i \over \partial x_j} \mu {\partial V_j \over \partial x_i}-\Delta t   {\partial W_i \over \partial x_j}{2\mu \over 3} \varepsilon _v \delta _{ij}- \Delta t   {\partial W_i \over \partial x_i}\Delta p \right]dV</math>
|-
| style="text-align: center;" | <math> = \!\!\int _V \left[W_i (2\rho )V_i^n + \Delta t W_i \rho a_i^n - \Delta t   {\partial W_i \over \partial x_j} \hat \sigma _{ij}^{(n)} + \Delta t    W_i \rho    f_i\right]dV + \Delta t \int _{\Gamma _\sigma } \bar \sigma _{ni} d\Gamma </math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (52)
|-
| style="text-align: center;" | <math> \int _V \left[\Delta t W_p {\partial V_i \over \partial x_i}+W_p {\Delta p\over    \kappa } \right]dV=0 </math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (53)
|}
|}

==5 SPATIAL DISCRETIZATION==

Both the velocity components and the pressure increment are discretized by

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>V_i = N^T Q_i\quad \hbox{and}\quad \Delta p = N_p \Delta P </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (54)
|}

The equation to be solved in matrix form becomes:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\left\{\left[\begin{matrix}M_{ij} & 0\\ 0 & - {1\over \kappa }       M_{pp}\end{matrix}\right]+\Delta t \left[\begin{matrix}K_{ij} & -       B_{ip}\\- B_{pj} & 0\end{matrix}\right]\right\}\left[\begin{matrix}Q_i\\ \Delta P\end{matrix}\right]= \left[\begin{matrix}G_i^n\\ 0\end{matrix}\right] </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (55)
|}

with

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\begin{array}{l}K_{ij} = K_{ij}^1 + K_{ij}^2 + K_{ij}^3\\[.2cm] M_{ii} = \int _V  N (2\rho )N^T dV \quad ;\quad K_{ii}^1 = \int _V \displaystyle {\partial N \over \partial x_j} (\mu ) \displaystyle {\partial N^T \over \partial x_j}dV \quad ;\quad M_{pp} = \int _V N_p N_p^TdV\\[.2cm] M_{ij} = K_{ij}^1 =0 \qquad \forall i\not = j\\[.2cm] K_{ij}^2 = \int _V \displaystyle {\partial N \over \partial x_j} (\mu ) \displaystyle {\partial N^T \over \partial x_i}dV\\[.2cm] K_{ij}^3 = \int _V \displaystyle {\partial N \over \partial x_i} (-{2\mu \over 3} ) \displaystyle  {\partial N^T \over \partial x_j}dV\\[.2cm] B_{ip} = \int _V \displaystyle {\partial N \over \partial x_i}N_p^TdV \quad ;\quad B_{pi}= \int _V N_p {\partial N^T \over \partial x_i}dV\\[.2cm]  G_i^n = \int _V\left[N (2\rho ) v_i^n + N (\Delta t\rho )a^n - \Delta t \displaystyle     {\partial N \over \partial x_j}\hat \sigma ^n_{ij} + N \Delta t\rho f_i \right]dV +  \int _{\Gamma _\sigma } N \Delta t \bar \sigma _{ni} d\Gamma   \end{array} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (56)
|}

It must be noted that this equation must be stabilized in order to avoid wiggles in the pressure solution due to the lack of the Babuska-Brezzi conditions. In this paper a Finite Calculus (FIC) formulation will be used to stabilize the solution [12,13,14,16].

At the end of each time step, the Cauchy stress <math display="inline">\sigma _{ij}^{n+1}</math>  are evaluated and saved for the next time step. At the beginning of each time step, the previous Cauchy stress are  considered as the second Piola-Kirchhoff stress for the present step and they are  evaluated by

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>S_{ij}^n \Leftarrow \sigma _{ij}^{n+1}</math>
|-
| style="text-align: center;" | <math> \hat \sigma ^n = (FS^n F) /J </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (57)
|}

where <math display="inline">F</math> is the deformation gradient tensor and <math display="inline">J</math> the <math display="inline">\det (F)</math> defined in Section 3.1.

==6 GENERATION OF A NEW MESH AND IDENTIFICATION OF THE BOUNDARY SURFACES==

One of the key points for the success of the Lagrangian  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 [7&#8211;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. 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 [7&#8211;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.

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 [4].

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  to the surface. Moreover, in weak forms (Galerkin) such as those used here a correct evaluation of the volume domain is  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. 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 when the lost node is at a distance less than <math display="inline">\alpha h</math> from the boundary.

==7 EXAMPLES==

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

Only 2D quasi-incompressible materials have been used in all the examples shown in the following sections.

The use of quasi-incompressible formulation to approximate an incompressible flow has been largely used in the literature. The advantage of this approximation is that no stabilization is necessary to obtain smooth solutions. In our unique lagrangian formulation for both, the elastic solid and the incompressible flow, the advantage is more evident because both domains: the solid and the fluid may be solved identical in this case. The only difference between the solid and the fluid is the constitutive equation but both do not need the pressure as a state variable.

There is some type of fluid-structure interaction problems, named gravitational problems in which the introduction of a speed of sound much smaller than the real one do not disturb too much the results. The large majority of free-surface problem where the free-surface is in contact with the atmospheric pressure are inside this kind of gravitational problems and for this reason, the introduction of the real speed of sound is unnecessary and we can take advantage of this factor. The reader is referred to previous publications of the authors [10,1] to see 2D and 3D examples of FSI between totally incompressible fluid flows and rigid solids.

In all the examples shown bellow, the mesh in the solid domain is generated only once and then the nodes are moved without re-meshing as in a classical finite element method. The PFEM approximation described before is only used in the fluid domain. In this way, the stresses <math display="inline">\hat \sigma _{ij}^{(n)}</math> from the previous time step remain at the element level for the next time step. Nevertheless the pressure values for the fluid domain are evaluated at the node (particles) to be transmitted to next time step with a new mesh.

===7.1 Elastic solid oscillation of a cantilever beam===

The unified formulation presented has been widely used in fluid dynamics problems using PFEM [6,10,18,20]. Lecturers are asked to read previous references for PFEM validation on fluid mechanics problems. Nevertheless, for structural problems the solution with PFEM is new. For this reason a simple a pure structural dynamic example has been tested. The free vibration of a cantilever beam subjected to a suddenly applied shear stress across the other end boundary face is shown in Figure [[#img-2|2]].  A plane strain beam with length <math display="inline">l=20</math>m, height <math display="inline">h= 5</math>, <math display="inline">E=4.00</math>GPa, Poisson's ratio = 0.32, <math display="inline">G=1.51</math>GPaN/cm<math display="inline">^2</math>  and density =1450kg/m<math display="inline">^3</math> was discretized using 1331 equal space particles. The total shear stress was equal to 1MPa with the upper and lower faces being traction free.

<div id='img-2'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Samper_577315608-barra_final.png|600px|Cantilever beam under a shear stress at the end length: Geometry]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 2:''' Cantilever beam under a shear stress at the end length: Geometry
|}

This example was presented in [27] and compare with an analytical solution for the 1D beam theory with a correction for longitudinal shear deformation.

The maximum displacement of the beam as time function is represented in Figure [[#img-3|3]] and compare with the analytical solution. Different time step were used in order to test de numerical diffusion of the method. Courant number equal 0.5 and 50 were tested. We can conclude that the approximation used in the solid part of the domain has a small numerical diffusion and that time step order of Courant number between 0.5 and 50 are enough to have excellent results.

<div id='img-3'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Samper_577315608-desplaz-barra.png|600px|Cantilever beam under a shear stress at the end length:  Maximum displacement for Courant numbers 0.5 and 50 ]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 3:''' Cantilever beam under a shear stress at the end length:  Maximum displacement for Courant numbers 0.5 and 50 
|}

===7.2 Impact of sea waves on solid object===

The collapse of a water column illustrated in Figure [[#img-4|4]] is calculated using the present formulation and  compared with experiment results obtained by S. Koshizuka [11]. Figure [[#img-5|5]] shows the experimental and the numerical results at different characteristic time step.

At time 0.1 sec, the right surfaces of the water start the disturbance due to the obstacle. At time 0.2 sec, the water surface is completely disturbed by the obstacle. The results compare very well with the experimental results. At 0.3 sec. the collapsed water crashes to the right wall. At 0.4 sec, the water goes up along the right wall with separations and several drops. Finally, at 1,0 sec, the water along the right wall falls down and a new breaking wave will soon occur on the left wall.

<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_577315608-dam.png|600px|Initial geometry of the water column]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 4:''' Initial geometry of the water column
|}

<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_577315608-Fig11a.png|360px|]]
|[[Image:Draft_Samper_577315608-dam011.png|300px|]]
|-
|[[Image:Draft_Samper_577315608-Fig11c.png|360px|]]
|[[Image:Draft_Samper_577315608-dam021.png|300px|]]
|-
|[[Image:Draft_Samper_577315608-Fig11e.png|360px|]]
|[[Image:Draft_Samper_577315608-dam031.png|300px|]]
|-
|[[Image:Draft_Samper_577315608-Fig11g.png|360px|]]
|[[Image:Draft_Samper_577315608-dam041.png|300px|]]
|-
|[[Image:Draft_Samper_577315608-Fig11i.png|360px|]]
|[[Image:Draft_Samper_577315608-dam101.png|300px|Comparison between experimental and numerical results of the collapse of a water column with an obstacle]]
|- style="text-align: center; font-size: 75%;"
| colspan="2" | '''Figure 5:''' Comparison between experimental and numerical results of the collapse of a water column with an obstacle
|}

After 0.4 sec a large deviation between the experimental results and the numerical one is observed. The reason is that the experimental results were performed in recipient with water in contact with the external air. After 0.4 sec, an air bubble is formed with the jet of water and the recipient. This incompressible air bubble changes significantly the experimental results.

In order to improve the numerical results, the same problem was solved including  the recipient water and air. Figure [[#img-6|6]] shows the new results in which the air is represented by a fluid with their corresponding density and viscosity. The results are now in a better agreement with the experimental ones and show the powerful of the PFEM to handle large differences in fluid physical characteristic.

<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_577315608-Figure4.png|600px|Comparison between experimental and numerical results of the collapse of a water column with an obstacle. Case with water and air.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 6:''' Comparison between experimental and numerical results of the collapse of a water column with an obstacle. Case with water and air.
|}

Figure [[#img-7|7]] represents the same problem but with a elastic obstacle with density <math display="inline">\rho = 2500 [kg/m^3]</math>, Young's modulus <math display="inline">E=10^6 [kg/s^2m]</math> and Poisson's ratio <math display="inline">\nu =0</math>. The geometry of the more slender obstacle is of width <math display="inline">h= 1.2 [cm]</math> and height <math display="inline">20/3\cdot h</math>. This example was also solved in Ref. 26 with a staggered and level-set method for the free surface flow. A unified formulation for incompressible flow and for the flexible obstacle has been used in this paper. No experimental results have been found for this example, but the shape of the deformation of the elastic beam as well as the free surface perturbation seems to be in agreement with the physics of the problem.

The left upper corner of the solid first gets a defection to the left when the water acts on its lower part and moves to the right while the water rises. It obtains its maximum deflection (mark (a) in Figure [[#img-8|8]]) when the water jet passes the top and is fully attached to the left side of the solid. Finally, the impact of the fluid causes the thin solid to start oscillating (b). This oscillation  is damped (c) by the water located left and right of the wall. The results obtained in Ref. 26 are also presented in Figure [[#img-8|8]].

<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_577315608-Figure_collapse2.png|600px|Collapse of a water column with an elastic obstacle]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 7:''' Collapse of a water column with an elastic obstacle
|}

<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_577315608-Grafico5.png|600px|Collapse of a water column with an elastic obstacle: History of   the displacement and comparison with Ref. 26]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 8:''' Collapse of a water column with an elastic obstacle: History of   the displacement and comparison with Ref. 26
|}

In order to increase the interaction between the fluid and the elastic solid, a similar problem to the previous one has been solved but using a taller elastic beam built-in on the bottom of the recipient. The obstacle with density <math display="inline">\rho = 7800 [kg/m^3]</math>, Young's modulus <math display="inline">E=2.1\cdot 10^6 [kg/s^2m]</math> and Poisson's ratio <math display="inline">\nu =0.3</math>. Geometry is illustrated in Figure [[#img-9|9]].

<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_577315608-barra.png|540px|Initial geometry of a high elastic bar under a lateral wave]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 9:''' Initial geometry of a high elastic bar under a lateral wave
|}

Now, the system is simulated for <math display="inline">2[s]</math> with a constant time step <math display="inline">= 0.001 [s]</math>. The hypoelastic solid is deflected by the impulse of the fluid wave (Figure [[#img-10|10]]).

<div id='img-10'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Samper_577315608-elastic_beam.png|600px|A high elastic beam built-in on the bottom of the recipient]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 10:''' A high elastic beam built-in on the bottom of the recipient
|}

===7.3 Elastic objects falling or floating in water===

The analysis of the motion of submerged or floating object in water is of great interest in many areas of harbour and costal engineering and naval architecture among others. In the following figures a series of examples will be described in order to show the potential of the PFEM to handle  fully coupled elastic solid and viscous fluid flows in problems in which the position of the body depends totally on the internal forces produces by the fluid. The presence of free surfaces introduces another complexity to the problem which is solved easily by PFEM.

Figure [[#img-11|11]] shows the penetration and motion of an elastic cube in a container with water. The difference in the density of the elastic body with the water results that the body will float on the free-surface. The density for the body used was 5% bigger than fluid. The colours denote the pressure distribution at different time instants. In this case, the speed of sound used for the water was <math display="inline">C= 10 [m/s]</math>.The small value for the sound speed may be observed in some time fluctuation of the pressure values, which for this kind of gravitational problems, do not affect considerably the geometric results.

Figure [[#img-12|12]] shows the evolution of an elastic floating plate on a free-surface with breaking waves. Here again a unified fluid-solid formulation has been used. The floating plate has a density 5% bigger than the fluid. For this reason the plate becomes totally submerged in the water.

<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_577315608-elastic_cube.png|600px|Elastic cube falling into a container with water]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 11:''' Elastic cube falling into a container with water
|}

<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_577315608-elastic-plate.png|600px|Dragging an elastic plate in a water container ]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 12:''' Dragging an elastic plate in a water container 
|}

==8 CONCLUSIONS==

A unified Lagrangian finite element framework in conjunction with the Particle Finite Element Method is ideal to treat problems involving fluids with free surfaces and submerged or floating structures. Problems such as 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 an updated Lagrangian formulation. Other essential solution ingredients are the efficient regeneration of the finite element mesh using an extended Delaunay tesselation, the meshless finite element interpolation (MFEM),  the identification of the boundary nodes using an Alpha Shape type technique and the simple algorithm to treat 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.

==ACKNOWLEDGEMENTS==

The authors thank Prof. J. Oliver and Dr. R. Rossi for many useful discussions. This work 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>
'''[1]'''   R. Aubry, S.R. Idelsohn, E. Oñate.  Particle finite element   method in fluid mechanics including thermal convection-diffusion. ''Computer   & Structures '', '''83''' (17-18), 1459&#8211;1475, (2005).

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

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

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

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

<div id="cite-6"></div>
'''[6]'''  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'', Mang HA,   Rammerstorfer FG, Eberhardsteiner J (eds), July 7&#8211;12, Viena, Austria, (2002).

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

<div id="cite-8"></div>
'''[8]'''  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''', 655&#8211;671, (2003b).

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

<div id="cite-10"></div>
'''[10]'''  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''', 964-989, (2004).

<div id="cite-11"></div>
'''[11]'''  S. Koshizuka, H. Tamko, Y. Oka, A particle method for   incompressible viscous flow with fluid fragmentation. ''Comp. Fluid Dynamic   Journal'', '''4''', (1), 29&#8211;46, (1995).

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

<div id="cite-13"></div>
'''[13]'''  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), 355&#8211;370, (2000).

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

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

<div id="cite-16"></div>
'''[16]'''   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''', 635&#8211;660, (2001).

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

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

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

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

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

<div id="cite-23"></div>
'''[23]'''  Ch. Farhat, G. Kristoffer, V. Zee, Ph. Geozine. Probably second-order time accurate loosely-coupled solution algorithm for transient nonlinear computational aeroelasticity. ''Comput. Meth. Appl. Mech. Engrg.'', '''195'''(17&#8211;18), 1973&#8211;2001, (2006).

<div id="cite-24"></div>
'''[24]'''  Ch. Foster, W.A. Wall, E. Ramm. Artificial added mass instabiliting in sequential staggered coupling of nonlinear structures and incompressible flows. ''Comput. Meth. Appl. Mech. Engrg.'', '''196''', 1278&#8211;1293, (2007).

<div id="cite-25"></div>
'''[25]'''  R. Lohner. A parallel advancing front grid generation scheme. ''AIAA'', 00-1005, (2000).

<div id="cite-26"></div>
'''[26]'''  E. Walhorn, A. Kolke, B. Hubner, D. Dinkler. Fluid-structure coupling with monolithic model involving free surface flows. ''Computers & Structures'', '''83''', 2100&#8211;2111, (2005).

<div id="cite-27"></div>
'''[27]'''  C.J. Greenshields, H.G. Weller. A unified formulation for continuum  mechanics applied to fluid-structure interaction in flexible  tubes. ''Int. J. Num. Meth. Engng.'', '''64''', 1575&#8211;1593. (2005).