Advances in the particle finite element method (PFEM) for solving coupled problems in engineering

==Abstract==

We present some developments in the formulation of the Particle Finite   Element Method (PFEM) for analysis of complex coupled problems on fluid and solid   mechanics in engineering accounting for fluid-structure interaction and coupled thermal effects, material degradation and surface wear. The   PFEM uses  an updated Lagrangian description to   model the motion of nodes (particles) in both the fluid and the structure   domains. Nodes are viewed as material points  which can freely move and even   separate from the main analysis domain representing, for instance,  the   effect of water drops. A mesh connects the nodes defining the discretized   domain where the governing equations are   solved, as in the standard FEM. The necessary stabilization for dealing with   the incompressibility of the fluid is introduced via the finite calculus   (FIC) method. An incremental iterative scheme for the solution of the non   linear transient coupled fluid-structure problem is described. The procedure for modelling frictional contact conditions at fluid-solid and   solid-solid interfaces via mesh generation  are   described. A simple algorithm to treat soil erosion in  fluid beds is   presented. An straight forward extension of the PFEM  to model excavation processes and wear of rock cutting tools is described. Examples of application of the PFEM  to solve a wide number of   coupled problems in engineering such as the effect of large waves on breakwaters and bridges, the large   motions of floating and submerged bodies, bed erosion in open channel flows, the wear of rock cutting tools during excavation and tunneling  and the melting,     dripping and burning of polymers in fire situations are presented.

==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-9'></span><span id='citeF-49'></span>[[#cite-9|[9]],[[#cite-49|49]]] and the references there included. A survey of recent works in fluid-structure interaction (FSI) analysis can be found in <span id='citeF-26'></span><span id='citeF-35'></span><span id='citeF-47'></span><span id='citeF-49'></span>[[#cite-26|[26]],[[#cite-35|35]],[[#cite-47|47]],[[#cite-49|49]]].

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 the past years a particular class of Lagrangian formulation for problems involving complex interactions between fluids and solids. The so called ''particle finite element method'' ([[http://www.cimne.com/pfem   PFEM]]), treats the  nodes in the fluid and solid domains as particles which can freely move and even separate from the main fluid (or solid) domain representing, for instance, the effect of water drops. A mesh connects the nodes discretizing the domain where the governing equations are solved using a stabilized FEM.

The FEM solution in the (incompressible) fluid domain implies solving the momentum and incompressibility equations. This is not a 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-9'></span><span id='citeF-49'></span>[[#cite-9|[9]],[[#cite-49|49]]]. 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.

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-13'></span><span id='citeF-15'></span>[[#cite-13|[13]],[[#cite-15|15]]]. The theory and applications of the PFEM are reported in <span id='citeF-2'></span><span id='citeF-8'></span><span id='citeF-13'></span><span id='citeF-14'></span><span id='citeF-16'></span><span id='citeF-17'></span><span id='citeF-34'></span><span id='citeF-35'></span><span id='citeF-36'></span><span id='citeF-38'></span><span id='citeF-40'></span><span id='citeF-44'></span><span id='citeF-45'></span><span id='citeF-46'></span>[[#cite-2|[2]],[[#cite-8|8]],[[#cite-13|13]],[[#cite-14|14]],[[#cite-16|16]],[[#cite-17|17]],[[#cite-34|34]],[[#cite-35|35]],[[#cite-36|36]],[[#cite-38|38]],[[#cite-40|40]],[[#cite-44|44]],[[#cite-45|45]],[[#cite-46|46]]].

The PFEM has been recently extended to model the frictional interaction between water and solids, as well as between deformable solids accounting for surface wear situations. Successful applications of the PFEM in this field include the modeling of bed erosion in free surface flows <span id='citeF-40'></span>[[#cite-40|[40]]], the simulation of excavation and tunneling problems and the study of wear in rock cutting tools <span id='citeF-5'></span><span id='citeF-6'></span>[[#cite-5|[5]],[[#cite-6|6]]].

Yet another successful application of the PFEM is the study of how objects melt, drip and burn in presence of fire. The solution of this complex FSI problem requires solving the equations of a coupled thermal-flow in a multifluid environment including an appropriate combustion model and taking into account the large deformations and eventual  loss of mass in the burning object <span id='citeF-24'></span><span id='citeF-41'></span><span id='citeF-46'></span>[[#cite-24|[24]],[[#cite-41|41]],[[#cite-46|46]]].

The aim of this paper is to describe recent advances of the PFEM for a) the the interaction between a collection of bodies which are fixed, floating and/or submerged in a fluid, b) the soil erosion in open channel flows, c) the wear of rock cutting tools and their performance during excavation and tunneling processes and d) the melting, dripping and burning of polymer objects in fire situations. All these problems are of great relevance in many areas of engineering. It is shown that the PFEM provides a general analysis methodology for treat such  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 thermal flow using a Lagrangian description and the FIC formulation are presented. Then an algorithm for the transient solution is briefly described. The treatment of the coupled FSI problem and 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. The extension of this erosion technique to model excavation in soil/rock and wear of rock cutting tools with the PFEM is presented. The potencial of the PFEM is shown in its application to FSI problems involving large flow motions, surface waves, moving bodies in water and bed erosion. Other examples shown include the application of PFEM to excavation and tunneling problems and to the melting, dripping and burning of polymers in fire situations.

==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 [[#img-1|1]]). The finite element method (FEM) is used to solve the continuum equations in both domains. Hence a mesh discretizing these domains must be generated in order to solve the governing equations for both the fluid and solid problems in the standard FEM fashion. 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%;"
|-
|[[Image:Draft_Samper_357825070-Figure1a.png|350px|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.

<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_357825070-Figure2a.png|400px|Sequence of steps to update a “cloud” of nodes representing a domain containing a fluid and a solid part 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 representing a domain containing a fluid and a solid part from time <math>n</math>   (<math>t=t_n</math>)  to   time <math>n+2</math> (<math>t=t_n +2\Delta t</math>)  
|}

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-2|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-10'></span>[[#cite-10|[10]]] is used for the boundary definition   (Section [[#5 Generation of a new mesh|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 Overview of the coupled FSI algoritm|4]]) <span id='citeF-13'></span><span id='citeF-14'></span><span id='citeF-16'></span>[[#cite-13|[13]],[[#cite-14|14]],[[#cite-16|16]]]. </li>

<li>Solve the coupled Lagrangian equations of motion for  the fluid and the   solid domains. Compute the state variables in both domains at the   next (updated) configuration for <math display="inline">t+\Delta t</math>: velocities, pressure,   viscous stresses and temperature in the fluid and   displacements, stresses, strains and temperature 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 [[#img-2|2]]. </li>

</ol>

Figure [[#img-3|3]] shows another conceptual example of application of the PFEM to modelling the melting and dripping of a polymer object under a heat source <math display="inline">q</math> acting at a boundary.

<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_357825070-Cloud_nodes.png|350px|Sequence of steps to update in time a “cloud” of nodes representing a polymer object under thermal loads represented by a prescribed boundary heat flux q. Crossed circles denote prescribed nodes at the boundary. The same figure explains the application of the PFEM to modelling a rock cutting problem]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 3:''' Sequence of steps to update in time a “cloud” of nodes representing a polymer object under thermal loads represented by a prescribed boundary heat flux <math>q</math>. Crossed circles denote prescribed nodes at the boundary. The same figure explains the application of the PFEM to modelling a rock cutting problem
|}

<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_357825070-Figure2a_b.png|400px|]]
|-
|[[Image:Draft_Samper_357825070-Fig2c.png|500px|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="1" | '''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
|}

Figure [[#img-3|3]] can be also used to explain the application of the PFEM to rock cutting problems. In those cases <math display="inline">q</math> represents the forces of the rock cutting tool acting on a rock mass represented by the cloud of points. The figure shows the detachment of the rock mass during the cutting process.

Figure [[#img-4|4]] shows a typical example of a PFEM solution of a free surface flow problem in 2D. The images correspond to the analysis of the problem of breakage of a water column <span id='citeF-16'></span><span id='citeF-36'></span>[[#cite-16|[16]],[[#cite-36|36]]]. Figure [[#img-4|4a]] shows the initial grid of four-noded rectangles discretizing the fluid domain and the solid walls.  Figures [[#img-4|4b]] and [[#img-4|4c]] show the deformed mesh at two later times.

==3 FIC/FEM formulation for a lagrangian incompressible thermal fluid==

===3.1 Governing equations===

The key equations to be solved in the incompressible thermal flow problem, written in the Lagrangian frame of reference, are the following:

'''Momentum'''

<span id="eq-1"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\rho {\partial v_i \over \partial t}={\partial \sigma _{ij} \over \partial x_j}+b_i\qquad \hbox{in } \Omega   </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (1)
|}

'''Mass balance'''

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

'''Heat transport'''

<span id="eq-3"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\rho c {\partial T \over \partial t}= {\partial  \over \partial x_i} \left(k {\partial T \over \partial x_i}\right)+Q \qquad \hbox{in } \Omega   </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (3)
|}

In above equations <math display="inline">v_{i}</math> is the velocity along the ''i''th global (cartesian) axis, <math display="inline">T</math> is the temperature, <math display="inline">\rho </math>, <math display="inline">c</math> and <math display="inline">k_i</math> are the density (assumed constant), the specific heat and the conductivity of the material, respectively, <math display="inline">b_i</math> and <math display="inline">Q</math> are the body forces and the heat source per unit mass, respectively and <math display="inline">\sigma _{ij}</math> are the (Cauchy) stresses related to the velocities by the standard constitutive equation (for incompressible Newtonian material)
<span id="eq-4"></span>
<span id="eq-4.a"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\sigma _{ij} =s _{ij} - p \delta _{ij}   </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (4.a)
|}

<span id="eq-4.b"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>s _{ij}=2\mu \left(\dot \varepsilon _{ij}-\frac{1}{3}\delta  _{ij}\dot \varepsilon _{ii} \right)\quad ,\quad \dot \varepsilon _{ij}=\frac{1}{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;" | (4.b)
|}

In Eqs.([[#eq-4|4]]), <math display="inline">s _{ij}</math> is the deviatoric stresses, <math display="inline">p</math>  is the pressure (assumed to be positive in compression), <math display="inline">\dot \varepsilon _{ij}</math>  is the rate of deformation, <math display="inline">\mu </math> is the viscosity and <math display="inline">\delta _{ij}</math> is the Kronecker delta. In the following we will assume the viscosity <math display="inline">\mu </math> to be a known function of temperature, i.e <math display="inline">\mu =\mu (T)</math>.

Indexes in Eqs.([[#eq-1|1]])-([[#eq-4|4]]) range from <math display="inline">i,j=1,n_{d}</math>, where <math display="inline">n_d</math>  is the number of space dimensions of the problem (i.e. <math display="inline">n_{d} = 2</math> for two-dimensional problems). The standard sum notation for repeated indices is assumed, unless otherwise specified.

Eqs.([[#eq-1|1]])-([[#eq-4|4]]) are completed with the standard boundary conditions of prescribed velocities and surface tractions in the mechanical problem and prescribed temperature and prescribed normal heat flux in the thermal problem <span id='citeF-2'></span><span id='citeF-9'></span>[[#cite-2|[2]],[[#cite-9|9]]].

We note that Eqs.([[#eq-1|1]])-([[#eq-3|3]]) are the standard ones for modeling the deformation of viscoplastic materials using the so called “flow approach” <span id='citeF-49'></span><span id='citeF-50'></span><span id='citeF-51'></span>[[#cite-49|[49]],[[#cite-50|50]],[[#cite-51|51]]]. In our work the dependence of the viscosity with the strain typical of viscoplastic flows has been simplified to the Newtonian form of Eq.([[#eq-4.b|4.b]]).

===3.2 Discretization of the equations===

A key problem in the numerical solution of Eqs.([[#eq-1|1]])-([[#eq-4|4]]) is the satisfaction of the incompressibility condition (Eq.([[#eq-2|2]])). A number of procedures to solve his problem exist in the finite element literature <span id='citeF-9'></span><span id='citeF-49'></span>[[#cite-9|[9]],[[#cite-49|49]]]. In our approach we use a stabilized formulation based in the so-called finite calculus procedure <span id='citeF-27'></span>[[#cite-27|[27]]&#8211;<span id='citeF-29'></span>[[#cite-29|29]],<span id='citeF-36'></span><span id='citeF-38'></span><span id='citeF-40'></span>[[#cite-36|36]],[[#cite-38|38]],[[#cite-40|40]]]. The essence of this method is the solution of a ''modified mass balance'' equation which is written  as

<span id="eq-5"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>{\partial v_i \over \partial x_i}+ \sum _{i=1}^{3}\tau {\partial  \over \partial x_i}\left[{\partial p \over \partial x_i} +\pi _i\right]=0 </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (5)
|}

where <math display="inline">\tau </math> is a stabilization parameter given by <span id='citeF-12'></span>[[#cite-12|[12]]]

<span id="eq-6"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\tau = \left(\frac{2\rho \vert \mathbf{v}\vert }{h}+\frac{8\mu }{3h^2} \right)^{-1}  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (6)
|}

In the above, <math display="inline">h</math> is a characteristic length of each finite element (such as <math display="inline">[A^{(e)}]^{1/2}</math> for 2D elements) and <math display="inline">\vert \mathbf{v}\vert </math> is the modulus of the velocity vector. In Eq.([[#eq-5|5]]) <math display="inline">\pi _i</math> are auxiliary pressure-gradient projection variables chosen so as to ensure that the second term in Eq.([[#eq-5|5]]) can be interpreted as weighted sum of the residuals of the momentum equations and therefore it vanishes for the exact solution. The set of governing equations for the velocities, the pressure and the <math display="inline">\pi _i</math> variables is completed by adding the following constraint equation to the set of governing equation <span id='citeF-36'></span><span id='citeF-40'></span>[[#cite-36|[36]],[[#eq-40|40]]]

<span id="eq-7"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\int _V \tau w_i\left({\partial p \over \partial x_i} +\pi _i\right)dV=0 \quad , \quad i=1,n_d \hbox{ not sum in } i  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (7)
|}

where <math display="inline">w_i</math> are arbitrary weighting functions.

The rest of the integral equations are obtained by applying the standard Galerkin technique to the governing equations ([[#eq-1|1]]), ([[#eq-2|2]]), ([[#eq-3|3]]), ([[#eq-5|5]]) and ([[#eq-7|7]]) and the corresponding boundary conditions <span id='citeF-36'></span><span id='citeF-40'></span>[[#cite-36|[36]],[[#eq-40|40]]].

We interpolate next in the standard finite element fashion the set of problem variables. For 3D problems these are the three velocities <math display="inline">v_i</math>, the pressure <math display="inline">p</math>, the temperature <math display="inline">T</math> and the three pressure gradient projections <math display="inline">\pi _i</math>.  In our work we use equal order ''linear interpolation'' ''for all variables'' over meshes of 3-noded triangles (in 2D) and 4-noded tetrahedra (in 3D) <span id='citeF-36'></span><span id='citeF-40'></span><span id='citeF-53'></span>[[#cite-36|[36]],[[#eq-40|40]],[[#cite-53|53]]]. The resulting set of discretized equations has the following form

'''Momentum'''

<span id="eq-8"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\mathbf{M} \dot{\bar{\boldsymbol v}} + \mathbf{K} (\mu )\bar {\boldsymbol v} - \mathbf{G} \bar  {\boldsymbol p}= {\boldsymbol f}  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (8)
|}

'''Mass balance'''

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

'''Pressure gradient projection'''

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

'''Heat transport'''

<span id="eq-11"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\mathbf{C}\dot{\bar{\boldsymbol T}} + \mathbf{H} \bar {\boldsymbol T} = {\boldsymbol q}  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (11)
|}

In Eqs.([[#eq-8|8]])&#8211;([[#eq-11|11]]) <math display="inline">\bar{(\cdot )}</math> denotes nodal variables and <math display="inline">\dot{\bar{(\cdot )}}=  {\partial  \over \partial t}\bar{(\cdot )}</math>. The different matrices and vectors are given in the Appendix.

The solution in time of Eqs.([[#eq-8|8]])&#8211;([[#eq-11|11]]) can be performed using any time integration scheme typical of the updated Lagrangian finite element method. A basic algorithm following the conceptual process described in Section [[#2.1 Basic steps of the PFEM|2.1]] is presented in Box [[#box-1|I]]. <math display="inline">^{t+\Delta t}(\bar{\boldsymbol a})^{j+1}</math> denotes the values of the nodal variables <math display="inline">\bar{\boldsymbol a}</math> at time <math display="inline">t+\Delta t</math> and the <math display="inline">j+1</math> iterations. We note the coupling of the flow and thermal equations via the dependence of the viscosity <math display="inline">\mu </math> with the temperature.

<span id="box-1"></span>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Samper_357825070-BoxIcon301.png|389px|]]
|}

==4 Overview of the coupled FSI algoritm==

Figure [[#img-5|5]] 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-5'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[File:Draft_Samper_357825070_5717_img-5a.JPG|330px|]]
|-
|[[Image:Draft_Samper_357825070-Figure3b.png|330px|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="1" | '''Figure 5:''' 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.
<span id="eq-12"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>{}^t S  :=  [{}^t {\boldsymbol x}_s,{}^t {\boldsymbol u}_s, {}^t {\boldsymbol v}_s, {}^t {\boldsymbol a}_s, {}^t {\boldsymbol \varepsilon }_s, {}^t {\boldsymbol \sigma }_s,{}^t {T}_s ]^T</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (12)
|}

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

where <math display="inline">{\boldsymbol x}</math> is the nodal coordinate vector, <math display="inline">{\boldsymbol u}</math>, <math display="inline">{\boldsymbol v}</math> and <math display="inline">{\boldsymbol a}</math> 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, <math display="inline">T</math> is the temperature and  subscripts <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 denotes nodal values.

The coupled fluid-structure interaction (FSI) problem of Figure [[#img-4|4]] 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>

</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 <span id='citeF-52'></span>[[#cite-52|[52]]]
<span id="eq-14"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>{\boldsymbol M}_s \bar {\boldsymbol a}_s +{\boldsymbol g}_s -  {\boldsymbol f}_s = {\boldsymbol 0} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (14)
|}

where <math display="inline">\bar {\boldsymbol a}_s</math> is the vector of nodal accelerations and <math display="inline">{\boldsymbol M}_s, {\boldsymbol g}_s</math> and <math display="inline">{\boldsymbol f}_s</math> are the mass matrix, the internal node force vector and the external nodal force vector in the solid domain. Indeed, the solid model can include any type of material and geometrical non-linearity using standard non-linear solid mechanics procedures <span id='citeF-52'></span>[[#cite-52|[52]]]. The time integration of Eq.([[#eq-14|14]]) is performed using a standard Newmark method.

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.   An incremental iterative scheme is implemented within each time step to account for non linear geometrical and material effects.

Iterate between 1 and 2 until convergence.

The above FSI solution algorithm is shown schematically in Box [[#box-2|II]].
<span id="box-2"></span>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[File:Draft_Samper_357825070_4896_box2.JPG|400px|]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Box II:'''Staggered solution scheme for the FSI problem (Figure [[#img-5|5]]). <math display="inline">S</math>: variables in the solid domain. <math display="inline">F</math>: variables in the fluid domain
|}

==5 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. Any fast meshing algorithm can be used for this purpose. In our work the mesh is generated at each time step using the extended Delaunay tesselation (EDT) <span id='citeF-13'></span><span id='citeF-15'></span><span id='citeF-16'></span>[[#cite-13|[13]],[[#cite-15|15]],[[#cite-16|16]]]. The EDT allows one to generate non standard meshes combining elements of  arbitrary polyhedrical shapes (triangles, quadrilaterals and other polygons in 2D and tetrahedra, hexahedra and arbitrary polyhedra in 3D) in a computing time of order <math display="inline">n</math>, where <math display="inline">n</math> is the total number of nodes in the mesh (Figure [[#img-6|6]]). The <math display="inline">C^\circ </math> continuous shape functions of the elements can be simply obtained using the so called meshless finite element interpolation (MFEM). In our work the simpler linear <math display="inline">C^\circ </math> interpolation has been chosen <span id='citeF-13'></span><span id='citeF-15'></span><span id='citeF-16'></span>[[#cite-13|[13]],[[#cite-15|15]],[[#cite-16|16]]].

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

Figure [[#img-7|7]] 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 <span id='citeF-40'></span>[[#cite-40|[40]]]. The figure shows the CPU time in seconds for each time step of the algorithm of Section [[#3.2 Discretization of the equations|3.2.]] 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 <math display="inline">20%</math> of the total CPU time for each time step, while the solution of the equations and the assembly of the system consume approximately <math display="inline">50%</math> and <math display="inline">20%</math> 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.

<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_357825070-etiquetas.png|340px|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 7:''' 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
|}

==6 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-10'></span>[[#cite-10|[10]]]. Figure [[#img-8|8]] shows an example of the boundary recognition using the Alpha Shape technique.

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

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 [[#img-2|2]] and [[#img-8|8]]).

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 the next section.

We note that the main difference between the PFEM and the classical FEM is just the remeshing technique and the identification of the domain boundary at each time step. The rest of the steps in the computation are coincident with those of the classical FEM.

==7 Treatment of contact  conditions in the PFEM==

===7.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-9|9]]). This simple way to treat the fluid-wall contact at mesh generation level is a distinct and attractive feature of the PFEM.

<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_357825070-Figure9.png|350px|]]
|-
|[[Image:Draft_Samper_357825070-Figure9b.png|350px|Automatic treatment of contact conditions at the fluid-wall interface]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 9:''' Automatic treatment of contact conditions at the fluid-wall interface
|}

===7.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-10|10]]).

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-11|11]]&#8211;[[#img-14|14]] show examples of application of the contact algorithm to the  bumping of a ball falling in a container, 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-11|11]] and [[#img-14|14]] show explicitely  the  layer of contact elements which controls the accuracy of the contact algorithm.

<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_357825070-Figure10b.png|400px|Contact conditions at a solid-solid interface]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 10:''' Contact conditions at a solid-solid interface
|}

<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_357825070-bola.png|311px|Bumping of a ball within a container. The layer of contact   elements is shown ]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 11:''' Bumping of a ball within a container. The layer of contact   elements is shown 
|}

<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_357825070-arch.png|334px|Failure of an arch   formed by  stone blocks under seismic loading]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 12:''' Failure of an arch   formed by  stone blocks under seismic loading
|}

<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_357825070-tetrapods-six.png|338px|Motion of five tetrapods on an inclined plane]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 13:''' Motion of five tetrapods on an inclined plane
|}

==8 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-22'></span><span id='citeF-48'></span>[[#cite-22|[22]],[[#cite-48|48]]]. The effect of water velocity on soil erosion was studied in <span id='citeF-42'></span>[[#cite-42|[42]]]. In a recent work we have proposed an extension of the PFEM to model bed erosion <span id='citeF-39'></span>[[#cite-39|[39]]]. 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-32'></span><span id='citeF-43'></span>[[#cite-1|[1]],[[#cite-32|32]],[[#cite-43|43]]].

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">\hat \tau </math> induced by the fluid motion. In 3D problems <math display="inline">\hat \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">{\boldsymbol n}</math> at the bed   node. The value of <math display="inline">\hat \tau </math> for 2D problems can be estimated as follows:
<span id="eq-15"></span>
<span id="eq-15a"></span>
{| 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;" | (15.a)
|}</li>

with
<span id="eq-15b"></span>
{| 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;" | (15.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-15|15]]).

<li>Compute the frictional work originated by the tangential stresses at the   bed surface as
<span id="eq-16"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>

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;" | (16)
|}</li>

Eq.([[#eq-16|16]]) is integrated in time using a simple scheme as
<span id="eq-17"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>

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

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

<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_357825070-detail-tetrapods.png|300px|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 14:''' Detail of five tetrapods on an inclined plane. The layer of elements modeling the   frictional contact conditions is shown
|}

<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_357825070-Erosion_fluid_forces.png|400px|Modeling of bed erosion by dragging of bed material]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 15:''' Modeling of bed erosion by dragging of bed material
|}

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

==9 Modelling and simulation of excavation and wear of rock cutting tools==

The PFEM has been successfully applied for modelling excavation processes in civil and mining engineering. The method can also accurately predict the wear of the rock cutting tools during the excavation.

The process to model surface erosion and tool wear during excavation follows the lines explained for modelling soil erosion in river beds (Section [[#8 Modeling of bed erosion|8]]). Material is removed from the excavation front or the tool surface when the work of the frictional forces at the rock/soil-tool interface exceeds a prescribed value. A new boundary is defined with the volume that remains in the analysis domain using the alpha-shape approach as it is typical in the PFEM (Section [[#6 Identification of boundary surfaces|6]]). The surface properties control the wear occurring during the frictional contact.

Mass loss in a cutting tool and the amount of excavated material that is extracted by the machine is modeled via a wear rate function. When a steady state position in the wear mechanism is reached, wear rate is described by a linear Archard-type equation <span id='citeF-1'></span><span id='citeF-5'></span><span id='citeF-43'></span>[[#cite-1|[1]],[[#cite-5|5]],[[#cite-43|43]]] as:
<span id="eq-18"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>V_{w}=K\;\frac{\| \mathbf{{f}}_{n}\| }{H}\;s </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (18)
|}

where <math display="inline">V_{w}</math> is the volume loss of the material along the contact surface due to wear, <math display="inline">s\,(m)</math> is the sliding distance, <math display="inline">\mathbf{{f}}_{n}</math> is the normal force vector to the contact surface and <math display="inline">H</math> is the hardness of the material. Constant <math display="inline">K</math> is a non-dimensional wear coefficient which depends on the relative contribution of the body under abrasion, adhesion and wear processes <span id='citeF-5'></span><span id='citeF-43'></span>[[#cite-5|[5]],[[#cite-43|43]]].

In the PFEM each node on the contact surface has a mesh of elements associated to it. The volume of material wear is compared with the volume associated to each contact node. When both volumes coincide, the node is released and all the elements associated to it are eliminated. The incremental equation for updating the volume loss due to wear at a node is as follows:
<span id="eq-19"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>V_{w}^{t+\triangle t}= V_{w}^{t}+K\;\frac{\| \mathbf{{f}}_{n}\| }{H}\;(\| \mathbf{{{v}}}_{t}\| \cdot \triangle  t) </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (19)
|}

where all variables are nodal variables, <math display="inline">\mathbf{{{v}}}_{t}</math> is the relative tangent velocity between the contact surfaces and <math display="inline">\triangle t</math> is the time step.

When the volume of worn material associated to a node and the volume of material are the same, the node is released. Elements that contain the released node are eliminated in the next time step. Some particles are also eliminated and hence the global volume of the problem changes. The historical value of the variables in these particles is lost as these particles do not contribute to the system anymore. A scheme of the geometry updating process  is shown in Figure&nbsp;[[#img-16|16]].

<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_357825070-4_Excavationscheme.png|400px|Removing material and boundary update in an excavation process]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 16:''' Removing material and boundary update in an excavation process
|}

The remeshing process allows the boundary recognition and the update of the analysis domain due to  excavation. The geometry of the  domain is changed at each time step as  excavation moves forward.

The flowchart for solving an excavation problem with the PFEM using an updated Lagrangian approach and an implicit integration scheme is the following:

<ol>

<li>Read initial geometrical, mechanical and kinematic conditions from a reference mesh.   </li>
<li>Transfer the elemental variables to the particles (i.e. the nodes).   </li>
<li>For each time step and each Newton iteration:

</li>
* Compute internal forces and contact forces at nodes

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>{\boldsymbol r} := {\boldsymbol M} \bar {\boldsymbol a}_s + {\boldsymbol g}_s + {\boldsymbol f}_c -{\boldsymbol f}_s</math>
|}
|}

* Compute displacement increments and update displacement values

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math> \delta \bar {\boldsymbol u}={\boldsymbol A}^{-1} {\boldsymbol r} \longrightarrow {}^{t+\Delta t}\Delta \bar{\boldsymbol u}^{i+1} = {}^{t+\Delta t}\Delta \bar{\boldsymbol u}^i + \delta  \bar{\boldsymbol u}</math>
|}
|}

where <math display="inline">{\boldsymbol A}</math> is the Jacobian matrix. Typically

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>{\boldsymbol A} = \frac{1}{\beta \Delta t^2}{\boldsymbol M}+ {\boldsymbol K}_T + {\boldsymbol K}_c</math>
|}
|}

where <math display="inline">\beta </math> is a parameter of the Newmark scheme <span id='citeF-52'></span>[[#cite-52|[52]]], <math display="inline">{\boldsymbol K}_T</math> is the tangent stiffness matrix of the solid mechanics problem accounting for material and non-linear geometrical effects <span id='citeF-5'></span><span id='citeF-6'></span><span id='citeF-52'></span>[[#cite-5|[5]],[[#cite-6|6]],[[#cite-52|52]]].

<li>Compute internal variables, strains and stresses at integration points within each element.   </li>
<li>Check convergence of Newton iterations.   </li>
<li>Once the iterative solutions has converged

</li>
* Update particle positions: <math display="inline">{}^{t+\Delta t} {\boldsymbol x} = {}^t {\boldsymbol x} + {}^{t+\Delta t}\Delta \bar{\boldsymbol u}</math>
* Compute velocities (<math display="inline">{}^{t+\Delta t} \bar{\boldsymbol v}</math>) and accelerations at particles (<math display="inline"> {}^{t+\Delta t} \bar{\boldsymbol a}</math>).
* Transfer strains and stresses from elements to particles:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>{}^{t+\Delta t}  {\boldsymbol \sigma }_p = {}^t {\boldsymbol \sigma }_p + \Delta {\boldsymbol \sigma }_p</math>
|}
|}

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>{}^{t+\Delta t}  {\boldsymbol \varepsilon }_p = {}^t {\boldsymbol \varepsilon }_p + \Delta {\boldsymbol \varepsilon }_p</math>
|}
|}

where <math display="inline">(\cdot )_p</math> denotes values at each particle. Note that ''the strain and stress history is stored at the particles''.
* Update constitutive law parameters.

<li>Check damage and erosion (wear) on particles. Remove eroded particles from the excavation front and worn particles from cutting tools.  </li>
<li>Boundary recognition via the alpha shape method. Create new mesh. Update problem dimensions if the number of particles has changed.  </li>
<li>Identify interface elements for contact.  </li>
<li>Initiate solution for next time step. </li>

</ol>

A detailed description of above algorithm, together with many applications, can be found in <span id='citeF-5'></span><span id='citeF-6'></span>[[#cite-5|[5]],[[#cite-6|6]]].

==10 Examples==

===10.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-17|17]] 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-18|18]]).

<div id='img-17'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Samper_357825070-fig16.png|400px|2D simulation of the penetration and evolution of a cube and a cylinder in a water container. The colours denote the different sizes of the elements at several times]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 17:''' 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-18'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Samper_357825070-Figure18.png|400px|Detail of element sizes during the motion of a rigid cylinder within   a water container]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 18:''' Detail of element sizes during the motion of a rigid cylinder within   a water container
|}

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

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

<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_357825070-diap17.png|405px|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 20:''' Impact of a wave on a prismatic column on a slab sustained by four   pillars.
|}

===10.3 Dragging of objects by water streams===

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

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

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

Figure [[#img-23|23]] shows the simulation of the falling of two tetrapods in a water container. Figure [[#img-24|24]] shows the motion of a collection of ten tetrapods placed in the slope of a breakwaters under an incident wave.

Figure [[#img-25|25]] shows a detail of the complex three-dimensional interactions between  water particles and tetrapods and between the tetrapods themselves.

Figures [[#img-26|26]] and [[#img-27|27]] show the analysis of the effect of breaking waves on two different sites of a breakwater containing reinforced concrete blocks (each one of <math display="inline">4\times 4</math> mts). The figures correspond to the study of Langosteira harbour in A Coruña, Spain using PFEM.

Figure [[#img-28|28]] displays the effect of an overtopping wave on a truck circulating by the perimetral road of the harbour adjacent to the breakwater.

<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_357825070-Fig22.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 22:''' Interaction of a wave with a vertical pier formed by   reinforced concrete cylinders.
|}

<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_357825070-Figure23.png|400px|Motion of two tetrapods falling in a water container.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 23:''' Motion of two tetrapods falling in a water container.
|}

<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_357825070-Figure24.png|400px|Motion of ten tetrapods on a slope under an incident wave.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 24:''' Motion of ten tetrapods on a slope under an incident wave.
|}

<div id='img-25'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|
[[File:Draft_Samper_357825070_1441_img-25.JPG|400px|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="1" | '''Figure 25:''' 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.
|}

<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_357825070-bloques2.png|400px|]]
|-
|[[Image:Draft_Samper_357825070-bloques1.png|400px|Effect of breaking waves on a breakwater slope containing reinforced concrete blocks. Detail of the mesh of 4-noded tetrahedra near the slope at two different times]]
|- style="text-align: center; font-size: 75%;"
| colspan="2" | '''Figure 26:''' Effect of breaking waves on a breakwater slope containing reinforced concrete blocks. Detail of the mesh of 4-noded tetrahedra near the slope at two different times
|}

<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_357825070-Dique_invernada2.png|300px|Study of breaking waves on the edge of a breakwater structure formed by reinforced concrete blocks]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 27:''' Study of breaking waves on the edge of a breakwater structure formed by reinforced concrete blocks
|}

<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_357825070-camion.png|400px|Effect of an overtopping wave on a truck passing by the perimetral road of a harbour adjacent to the breakwater]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 28:''' Effect of an overtopping wave on a truck passing by the perimetral road of a harbour adjacent to the breakwater
|}

===10.5 Soil erosion===

Figure [[#img-29|29]] shows a very illustrative example of the potential of the PFEM to model soil 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   show the progressive erosion of the dam until the whole dam is dragged out by the fluid flow <span id='citeF-39'></span>[[#cite-39|[39]]].

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

<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_357825070-Erosion4.png|500px|Erosion, transport and deposition of particles at a river bed due to a jet stream.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 30:''' Erosion, transport and deposition of particles at a river bed due to a jet stream.
|}

Figure [[#img-30|30]] shows the capacity of the PFEM to modelling soil erosion, sediment transport and material deposition in a river bed. The soil particles are first detached from the bed surface under the action of the jet stream. Then they are transported by the flow and eventually  fall down due to gravity forces and are deposited on the bed surface at a downstream point.

Figure [[#img-31|31]] shows the progressive erosion of the unprotected part of a break water slope in the Langosteira harbour in A Coruña, Spain. Note that the upper shoulder zone not protected by the concrete blocks is progressively eroded under the action of the sea waves.

Figure [[#img-32|32]] displays the progressive erosion and dragging of soil particles in a river  bed adjacent to the foot of bridge pile due to a water stream (water is not shown in the figure). Note the disclosure of the bridge foundation due to the removal of the adjacent soil due to erosion.

<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_357825070-Erosion2D.png|357px|Erosion of unprotected part of a breakwater slope due to sea waves.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 31:''' Erosion of unprotected part of a breakwater slope due to sea waves.
|}

<div id='img-32'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Samper_357825070-Erosionpila.png|400px|Progressive erosion and dragging of soil particles in a river bed adjacent to the foot of a bridge pile due to a water stream. Water is not shown.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 32:''' Progressive erosion and dragging of soil particles in a river bed adjacent to the foot of a bridge pile due to a water stream. Water is not shown.
|}

===10.6 Melting, spread and burning of polymer objects in fire===

We show an application of the PFEM for simulating an experiment performed at the National Institute for Stanford and Technology (NIST) in which a slab of polymeric material is mounted vertically and exposed to uniform radiant heating on one face. It is assumed that the polymer melt flow is governed by the equations of an incompressible fluid with a temperature dependent viscosity. A quasi-rigid behaviour of the polymer object at room temperature is reproduced by using a very high value of the viscosity parameter. As temperature increases in the thermoplastic object due to heat exposure, the viscosity decreases in several orders of magnitude as a function of temperature and this induces the melt and flow of the particles in the heated zone. Polymer melt is captured by a pan below the sample.

A rectangular polymeric sample of dimensions 10 cm high by 10 cm wide by 2.5 cm thick is mounted upright and exposed to uniform heating on one face from a radiant cone heater placed on its side (Figure [[#img-33|33]]). The sample is insulated on its lateral and rear faces. The melt flows down the heated face of the sample and drips onto a surface below.  Measurements include the mass of polymer remaining in the sample, and the mass of polymer falling onto the catch surface  <span id='citeF-4'></span>[[#cite-4|[4]]].

Figure [[#img-33|33]] shows all three curves of viscosity vs. temperature for the polypropylene type PP702N, a low viscosity commercial injection molding resin formulation.  The relationship used in the model, as shown by the black line, connects the curve for the undegraded polymer to points A and B extrapolated from the viscosity curve for each melt sample to the temperature at which the sample was formed.  The result is an empirical viscosity-temperature curve that implicitly accounts for molecular weight changes.

<div id='img-33'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Samper_357825070-Viscosity_temperature.png|500px|Polymer melt experiment. Viscosity vs. temperature for PP702N   polypropylene in its initial undegraded form and after exposure to 30   kW/m² and 40 kW/m² heat fluxes.  The black curve follows the   extrapolation of viscosity to high temperatures.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 33:''' Polymer melt experiment. Viscosity vs. temperature for PP702N   polypropylene in its initial undegraded form and after exposure to 30   kW/m<math>^{2}</math> and 40 kW/m<math>^{2}</math> heat fluxes.  The black curve follows the   extrapolation of viscosity to high temperatures.
|}


The finite element mesh has 3098 nodes and 5832 triangular elements.  No nodes are added during the course of the run. The addition of a catch pan to capture the dripping polymer melt tests the ability of the PFEM model to recover mass when a particle or set of particles reaches the catch surface.  Heat flux is only applied to free surfaces above the midpoint between the catch pan and the base of the sample.  However, every free surface is subject to radiative and convective heat losses.  To keep the melt fluid, the catch pan is set to a temperature of 600 K. Figure [[#img-34|34]] shows four snapshots of the melt flow into the catch pan.

To test the ability of the PFEM to solve this type of problem in three dimensions, a 3D problem for flow from a heated sample was run. The same boundary conditions are used as in the 2D problem illustrated in Figure [[#img-33|33]], but the initial dimensions of the sample are reduced to <math display="inline">10\times 2.5 \times 2.5</math> cm.  The initial size of the model is 22475 nodes and 97600 four-noded tetrahedra. The shape of the surface and temperature field at different times after heating begins are shown in Figure [[#img-35|35]].

<div id='img-34'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[File:Draft_Samper_357825070_5873_img-34a.JPG|248px|]]
|[[Image:Draft_Samper_357825070-Evolution_melt_2.png|248px|]]
|-
|[[File:Draft_Samper_357825070_6866_img-34c.JPG|248px|]]
|[[Image:Draft_Samper_357825070-Evolution_melt_4.png|250px|Polymer melt experiment. Evolution of the melt flow into the catch pan at t = 400s, 550s, 700s and 1000s]]
|- style="text-align: center; font-size: 75%;"
| colspan="2" | '''Figure 34:''' Polymer melt experiment. Evolution of the melt flow into the catch pan at t = 400s, 550s, 700s and 1000s
|}

Although the resolution for this problem is not fine enough to achieve high accuracy, the qualitative agreement of the 3D model with 2D flow and the ability to carry out this problem in a reasonable amount of time suggest that the PFEM can be used to model melt flow and spread of complex 3D polymer geometry.

Figure [[#img-36|36]] shows results for the analysis of the melt flow of a triangular thermoplastic object into a catch pan. The material properties for the polymer are the same as for the previous example.  The PFEM succeeds to predicting in a very realistic manner the progressive melting and slip of the polymer particles along the vertical wall separating the triangular object and the catch pan. The analysis follows until the whole object has fully melt and its mass is transferred to the catch pan.

We note that the total mass was preserved with an accuracy of 0.5% in all these studies. Gasification, in-depth absorption or radiation were not taken into account in these analysis. More examples of application of the PFEM to the melting and dripping of polymers are reported in <span id='citeF-41'></span>[[#cite-41|[41]]].

<div id='img-35'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Samper_357825070-Fig7_polymers.png|400px|Simulation of a 3D polymer melt problem with the PFEM. Melt flow from a   heated prismatic sample at different times.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 35:''' Simulation of a 3D polymer melt problem with the PFEM. Melt flow from a   heated prismatic sample at different times.
|}

<div id='img-36'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Samper_357825070-Fig9_polymers.png|451px|Melt flow of a heated triangular object into a catch pan.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 36:''' Melt flow of a heated triangular object into a catch pan.
|}

<div id='img-37'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Samper_357825070-Silla.png|600px|Simulation of the burning, melting and dripping of a chair modelled as a 2D prismatic polymer object.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 37:''' Simulation of the burning, melting and dripping of a chair modelled as a 2D prismatic polymer object.
|}

The PFEM has been recently extended for modelling the combined melting and burning of polymer objects under fire. The equation governing the coupled thermal-flow problem are extended with a combustion model governing the burning of combustible and the heat interchanges between the object and  the air during  combustion <span id='citeF-21'></span>[[#cite-21|[21]],<span id='citeF-24'></span>[[#cite-24|24]],<span id='citeF-46'></span>[[#cite-46|46]]]. Figure [[#img-37|37]] shows a 2D application of the PFEM to the burning of a prismatic polymer object simulating a chair. The sequence of images shows the change of shape of the object as it burns, melts and drips on the floor surface and the intensity of the flame at different times.

===10.7 Simulation of excavation process and wear of rock cutting tools===

''Disc cutting of a ground section''

The first example is an elastic cutting disc in 2D acting against a solid wall. The disc has an imposed rotation in order to generate friction when contacting with the  solid wall. The material is modelled with a simple damage law.

The problem is solved first for the case of a soft wall material. Figure [[#img-38|38]] shows that contact is detected when the disc comes near the wall. An interface mesh of contact elements is generated and it anticipates the contact area. The contacting forces are transmitted thought the contact elements to each domain. This interaction damages the solid wall until it crashes. Contact forces are computed at the axis of the disc in order to yield force and momentum reactions.

The mesh is coarse so as to show better the process and the contact interface mesh. In a fine mesh contact elements are quite small and are difficult to visualize. It can be seen how as contact forces erode the wall, the excavated particles are taken away from the model. This generates a hollow in the surface while at the same time the material experiences large deformations. Figures [[#img-39|39]] and [[#img-40|40]] show a similar examples of excavation of a soft soil mass with rotating discs.

Figure [[#img-41|41]] displays the action of a rotating disc on a stiff wall. Note the change in the pattern of the excavation front and the progressive wear of the disc surface.


<div id='img-38'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Samper_357825070-5_Example2D.png|400px|Simulation of a disc excavating a soft wall with the PFEM]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 38:''' Simulation of a disc excavating a soft wall with the PFEM
|}

<div id='img-39'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Samper_357825070-Excavation_wear_modulus2.png|400px|Example of application of the PFEM to the excavation of a soft soil mass with a rotating disc]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 39:''' Example of application of the PFEM to the excavation of a soft soil mass with a rotating disc
|}

<div id='img-40'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Samper_357825070-Excavation_acceleration_two.png|400px|Simulation of the excavation of a soft soil mass with a rotating gear disc with the PFEM. Contour of the modulus of the acceleration vector in the soil at two instances ]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 40:''' Simulation of the excavation of a soft soil mass with a rotating gear disc with the PFEM. Contour of the modulus of the acceleration vector in the soil at two instances 
|}

<div id='img-41'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Samper_357825070-Excavation_disc_three.png|400px|Simulation of the excavation of a stiff rock wall with the PFEM. Note the change of the rotating disc edge due to wear]]
|- style="text-align: center; font-size: 75%;"
| colspan="2" | '''Figure 41''' Simulation of the excavation of a stiff rock wall with the PFEM. Note the change of the rotating disc edge due to wear 
|}

''Roadheader penetrating in the ground''

The next example is the simulation of a roadheader digging a portion of ground. This is an illustrative example of the capability of the PFEM for modeling ground excavation and wear of the cutting tools at the same time.

The results are shown in Figure [[#img-42|42]]. A rotation and a displacement have been imposed to the roadheader. Note that contact elements only appear in the contact zone. The cone that models the roadheader loses material at the tip due to wear. Ground geometry suffers big changes during the simulation. Remeshing and detection of the boundary via the alpha-shape technique are crucial for capturing the fast and drastic changes of the domain boundary.

<div id='img-42'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Samper_357825070-7_Roadheader3D.png|400px|Simulation of an excavation with a roadheader using the PFEM. Note the geometry change in the roadheader tip due the wear]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 42:'''Simulation of an excavation with a roadheader using the PFEM. Note the geometry change in the roadheader tip due the wear 
|}

''Simulation of an excavation with a TBM''

Figures [[#img-43|43]]&#8211;[[#img-45|45]] show a simulation of a tunneling process with a TMB (Tunnel Boring Machine) acting on a 3D soil/rock domain. This  example evidences the capability of the PFEM  to model complex excavation settings. The discretization of the TMB  and the soil/rock region is displayed in Figure [[#img-43|43]]. Figure [[#img-44|44]] shows an overview of the simulation as the tunneling process advance and the stress contour lines and Figure [[#img-45|45]] shows the wear of the rock cutting discs in the TBM induced by the excavation forces. Far away from the rotating axis the displacement is bigger for the same rotation velocity and it generates larger friction forces at the edges of the tunneling head.

<div id='img-43'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Samper_357825070-Fig1_particles.png|400px|Simulation of a tunneling process with a TBM using the PFEM. Discretization of soil mass and TBM geometry with 4-noded tetrahedra]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 43:'''Simulation of a tunneling process with a TBM using the PFEM. Discretization of soil mass and TBM geometry with 4-noded tetrahedra
|}

<div id='img-44'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Samper_357825070-8_TBM3D.png|400px|Simulation of a tunneling operation with a TBM using the PFEM]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 44:''' Simulation of a tunneling operation with a TBM using the PFEM
|}

<div id='img-45'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Samper_357825070-Fig2_particles.png|400px|Wear of the rock cutting discs in a TBM during the simulation of a tunneling operation using the PFEM. Circles denote worn cutting discs.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 45:''' Wear of the rock cutting discs in a TBM during the simulation of a tunneling operation using the PFEM. Circles denote worn cutting discs.
|}

The previous examples illustrate the good  capabilities of the PFEM for modelling ground excavation processes.

==11 Conclusions==

The particle finite element method (PFEM) is a powerful computational technique for solving coupled problems in engineering,  involving fluid-structure interaction, large motion of  fluid or solid particles, surface waves, water splashing, separation of water drops, frictional  contact situations, bed erosion, coupled thermal flows, melting, dripping and burning of objects, etc. The success of the PFEM lies in the accurate and efficient solution of the equations of fluid and of solid mechanics 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 identification of the domain boundaries via the Alpha Shape technique and the efficient regeneration of the finite element mesh at each time step, and the  algorithm to treat frictional contact conditions at fluid-solid and solid-solid interfaces via mesh generation. The examples presented have shown the  potential of the PFEM for solving a wide class of coupled problems in engineering. Examples of validation of the PFEM results with data from experimental tests are reported in <span id='citeF-23'></span>[[#cite-23|[23]]].

==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 (MEC) of Spain, project XPRES of the National I+D Programme of MEC (Spain) and projects REALTIME and SAFECON of the European Research Council (ERC). Thanks are also given to the Spanish construction company Dragados for financial support for the study of harbour engineering and tunneling problems.

==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>
[[#citeF-3|[3]]] K.M. Butler, T.J. Ohlemiller, G.T. Linteris, A Progress Report on Numerical Modeling of Experimental Polymer Melt Flow Behavior, Interflam (2004) 937&#8211;948.

<div id="cite-4"></div>
[[#citeF-4|[4]]]  K.M. Butler, E. Oñate, S.R. Idelsohn, R. Rossi, Modeling and simulation of the melting of polymers under fire conditions using the particle finite element method, 11th Int. Fire Science & Engineering Conference, University of London, Royal Halbway College, UK, (2007) 3-5 September.

<div id="cite-5"></div>
[[#citeF-5|[5]]]  J.M. Carbonell, Modeling of ground excavation with the Particle Finite Element method. Ph.D. Thesis, Technical University of Catalonia (UPC), Barcelona, (2009).

<div id="cite-6"></div>
[[#citeF-6|[6]]]  J.M. Carbonell, E. Oñate, B. Suárez, Modeling of ground excavation with the Particle Finite Element method. Journal of Engineering Mechanics (ASCE), April (2010).

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

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

<div id="cite-11"></div>
[[#citeF-11|[11]]]  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-12"></div>
[[#citeF-11|[11]]]  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-13"></div>
[[#citeF-13|[13]]]  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-14"></div>
[[#citeF-14|[14]]]  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-15"></div>
[[#citeF-15|[15]]]  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-16"></div>
[[#citeF-16|[16]]]  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-17"></div>
[[#citeF-17|[17]]]  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-18"></div>
[[#citeF-18|[18]]]  S.R. Idelsohn, J. Marti, A. Limache, E. Oñate, Unified Lagrangian formulation for elastic solids and incompressible fluids: Application to fluid-structure interaction problems via the PFEM. Comput Methods Appl Mech Engrg. (2008) 197 1762&#8211;1776.

<div id="cite-19"></div>
[[#citeF-19|[19]]]  S.R. Idelsohn, M. Mier-Torrecilla, E. Oñate, Multi-fluid flows with the Particle Finite Element Method. Comput Methods Appl Mech Engrg. 198 (2009) 2750&#8211;2767.

<div id="cite-20"></div>
[[#citeF-20|[20]]]  S.R. Idelsohn, M. Mier-Torrecilla, N. Nigro, E. Oñate, On the analysis of heterogeneous fluids with jumps in the viscosity using a discontinuous pressure field. Comput. Mech. (2010) 46 (1) 115&#8211;124.

<div id="cite-21"></div>
[[#citeF-21|[21]]]  S.R. Idelsohn, J. Marti, E. Oñate, R. Rossi, K. Butler, A flame model for melting and dripping of polymers. 12th International Interflam Fire Science and Engineering Conference, 5-7 July 2010, Nottingham, UK.

<div id="cite-22"></div>
[[#citeF-22|[22]]]  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-23"></div>
[[#citeF-23|[23]]]  A. Larese,  R. Rossi, E. Oñate, S.R. Idelsohn,  Validation of the Particle Finite Element Method (PFEM) for simulation of free surface flows, Engineering Computations 25 (4) (2008) 385&#8211;425.

<div id="cite-24"></div>
[[#citeF-24|[24]]]  J. Marti, P. Ryzhakov, S.R. Idelsohn, E. Oñate, V. Novozhilov, A new approach for simulation of the polymers in fire situations. International Congress on Combustion and Fire Dynamics, 20 -23 October 2010, Santander, Spain.

<div id="cite-25"></div>
[[#citeF-25|[25]]]  M. de Mier Torrecilla, Numerical Simulation of Multi-Fluid Flows with the Particle Finite Element Method. ''Ph.D. Thesis'', Technical  University of Catalonia (UPC), July 2010.

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

<div id="cite-30"></div>
[[#citeF-30|[30]]]  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-31"></div>
[[#citeF-31|[31]]]   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-32"></div>
[[#citeF-32|[32]]]  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-33"></div>
[[#citeF-33|[33]]]  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-34"></div>
[[#citeF-34|[34]]]  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-35"></div>
[[#citeF-35|[35]]]  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, J. Wiley, Vol 3,   (2004a) 579&#8211;610.

<div id="cite-36"></div>
[[#citeF-36|[36]]]  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-37"></div>
[[#citeF-37|[37]]]  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-38"></div>
[[#citeF-38|[38]]]  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-39"></div>
[[#citeF-39|[39]]]  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-40"></div>
[[#citeF-40|[40]]]  E. Oñate, S.R. Idelsohn, M.A. Celigueta, R. Rossi, Advances in the particle finite element method for the analysis of fluid-multibody interaction and bed erosion in free surface flows, Comput. Meth. Appl. Mech. Engng. 197 (19-20) (2008) 1777-1800.

<div id="cite-41"></div>
[[#citeF-41|[41]]]  E. Oñate, R. Rossi, S.R. Idelsohn, K. Butler, Melting and spread of polymers in fire with the particle finite element method. Int. J. Num. Meth. in Engng., 81 (8) (2010) 1046-1072.

<div id="cite-42"></div>
[[#citeF-42|[42]]]  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-43"></div>
[[#citeF-43|[43]]]  Rabinowicz, E., ''Friction and Wear of materials''. Wiley, (1995).

<div id="cite-44"></div>
[[#citeF-44|[44]]]  R. Rossi, P.B. Ryzhakov, E. Oñate, A monolithic FE formulation for the analysis of membranes in fluids. Journal of Spatial Structures 24 (4) (2009) 205&#8211;210.

<div id="cite-45"></div>
[[#citeF-45|[45]]]  P.B. Ryzhakov, R. Rossi, S. Idelsohn, E.Oñate, A monolithic Lagrangian approach for fluid-structure interaction problems. Journal of Computational Mechanics 46 (6) (2010) 883&#8211;899.

<div id="cite-46"></div>
[[#citeF-46|[46]]]  P.B. Ryzhakov, R. Rossi, E. Oñate, An algorithm for polymer melting simulation. Conference Proceedings METNUM-2009,  Barcelona, Spain, 29 June -02 July (2009).

<div id="cite-47"></div>
[[#citeF-47|[47]]]  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, 3, (J. Wiley, 2004) 545&#8211;578.

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

<div id="cite-50"></div>
[[#citeF-50|[50]]]  O.C. Zienkiewicz, P.C. Jain, E. Oñate, Flow of solids during forming and extrusion: Some aspects of numerical solutions. Int. Journal of Solids and Structures  14 (1978) 15&#8211;38.

<div id="cite-51"></div>
[[#citeF-51|[51]]]  O.C. Zienkiewicz, E. Oñate, J.C. Heinrich, A general formulation for the  coupled thermal flow of metals using finite elements. Int. Journal for Numerical Methods in  Engineering  17 (1981) 1497&#8211;1514.

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

<div id="cite-53"></div>
[[#citeF-53|[53]]]  O.C. Zienkiewicz, R.L. Taylor, J.Z. Zhu,  The finite element method. Its basis and fundamentals,  Elsevier, (2005).

==APPENDIX==

The  matrices and vectors in Eqs.([[#eq-8|8]])-([[#cite-11|11]]) for a 4-noded tetrahedron are:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\mathbf{M}_{ij} =\int _{V^{e}} \rho \mathbf{N}_{i}^{T} \mathbf{N}_{j} dV \quad , \quad \mathbf{K}_{ij}=\int _{V^{e}} \mathbf{B}_{i}^{T} \mathbf{D} \mathbf{B}_{j} dV </math>
|}
|}

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math> \mathbf{G}_{ij} = \int _{V^{e}} \mathbf{B}_{i}^{T}  \mathbf{m}{N}_{j} dV \quad , \quad \mathbf{f}_{i}=\int _{V^{e}} \mathbf{ N}_{i}^{T}\mathbf{b}dV+ \int _{\Gamma ^e} \mathbf{ N}_{i}^{T}\mathbf{t}d\Gamma \quad , \quad \hat {\boldsymbol M}_{ij} = \int _{V^{e}} \tau \mathbf{N}_{i}^{T} \mathbf{N}_{j} dV </math>
|}
|}

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>{L}_{ij} = \int _{V^e} {\boldsymbol \nabla }^T N_i \tau {\boldsymbol \nabla }N_j dV  \quad ,\quad {\boldsymbol \nabla } = \left[{\partial  \over \partial x_1},{\partial  \over \partial x_2},{\partial  \over \partial x_3}\right]^T</math>
|}
|}

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>{\boldsymbol Q}= [{\boldsymbol Q}_1,{\boldsymbol Q}_2,{\boldsymbol Q}_3]\quad ,\quad [\mathbf{Q}_{k}]_{ij} = \int _{V^e}\tau {\partial N_i \over \partial x_k} \mathbf{N}_j dV \quad ,\quad  {\boldsymbol m} =[1,1,1,0,0,0]^T</math>
|}
|}

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>{\boldsymbol B}=[{\boldsymbol B}_1,{\boldsymbol B}_2,{\boldsymbol B}_3,{\boldsymbol B}_4];\, {\boldsymbol B}_i =\left[\begin{matrix} \displaystyle {\partial N_i \over \partial x}&0&0\\ 0 & \displaystyle {\partial N_i \over \partial y} &0\\ 0&0& \displaystyle {\partial N_i \over \partial z} \\ \displaystyle {\partial N_i \over \partial y} & \displaystyle {\partial N_i \over \partial x}&0\\\\ \displaystyle {\partial N_i \over \partial z} &0&\displaystyle {\partial N_i \over \partial x}\\\\ 0& \displaystyle {\partial N_i \over \partial z} &\displaystyle {\partial N_i \over \partial y} \end{matrix} \right]\quad ,\quad {\boldsymbol D} =\mu \left[\begin{matrix}2{\boldsymbol I}_3 & {\boldsymbol 0}\\ {\boldsymbol 0} & {\boldsymbol I}_3 \end{matrix}\right]</math>
|}
|}

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>{\boldsymbol N}=[{\boldsymbol N}_1,{\boldsymbol N}_2,{\boldsymbol N}_3,{\boldsymbol N}_4] \quad , \quad  {\boldsymbol N}_i={N}_i {\boldsymbol I}_3\quad ,\quad  {\boldsymbol I}_3: \quad 3\times 3 \hbox{ unit matrix}</math>
|}
|}

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>{C}_{ij}=\int _{V^e} \rho {c}{N}_i {N}_j dV \quad ,\quad  {H}_{ij}= \int _{V^e}{\boldsymbol \nabla }^T {N}_i [k] {\boldsymbol \nabla }{N}_j dV </math>
|}
|}

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>[k]=  \left[\begin{matrix}k_1 &0 &0 \\ 0 & k_2&0\\ 0&0& k_3\end{matrix}\right]\quad , \quad   q_i =\int _{V^e}N_iQdV - \int _{\Gamma _q^{e}}  N_i q_n d\Gamma </math>
|}
|}

In the above equations indexes <math display="inline">i,j</math> run from 1 to the number of element nodes (4 for a tetrahedron), <math display="inline">q_n</math> is the heat flow prescribed at the external boundary <math display="inline">\Gamma _q</math>, '''t''' is the surface traction vector <math display="inline">\mathbf{t}=[t_x,t_y,t_z]^T</math> and <math display="inline">V^{e}</math> and <math display="inline">\Gamma ^e</math> are the element volume and the element boundary, respectively.