Published in ParticleBased Methods, Computational Methods in Applied Sciences, Springer, E. Oñate and R. Owen (Eds.), Vol. 25, pp. 148, 2011
DOI: 10.1007/9789400707351_1
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 fluidstructure 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 fluidstructure problem is described. The procedure for modelling frictional contact conditions at fluidsolid and solidsolid 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.
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, offshore and harbour structures, spillways in dams, free surface channel flows, environmental flows, liquid containers, stirring reactors, mould filling processes, etc.
Typical difficulties of fluidmultibody 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 [9,49] and the references there included. A survey of recent works in fluidstructure interaction (FSI) analysis can be found in [26,35,47,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 ([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 stability condition [9,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 [13,15]. The theory and applications of the PFEM are reported in [2,8,13,14,16,17,34,35,36,38,40,44,45,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 [40], the simulation of excavation and tunneling problems and the study of wear in rock cutting tools [5,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 thermalflow 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 [24,41,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 fluidwall 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.
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 . The new set of variables in both domains are sought for in the next or updated configuration at time (Figure 1). The finite element method (FEM) is used to solve the continuum equations in both domains. Hence a mesh discretizing these domains must be generated in order to solve the governing equations for both the fluid and solid problems in the standard FEM fashion. Recall that the nodes discretizing the fluid and solid domains are treated as material particles which motion is tracked during the transient solution. This is useful to model the separation of fluid particles from the main fluid domain in a splashing wave, or soil particles in a bed erosion problem, and to follow their subsequent motion as individual particles with a known density, an initial acceleration and velocity and subject to gravity forces. The mass of a given domain is obtained by integrating the density at the different material points over the domain.
The quality of the numerical solution depends on the discretization chosen as in the standard FEM. Adaptive mesh refinement techniques can be used to improve the solution in zones where large motions of the fluid or the structure occur.
Figure 1: Updated lagrangian description for a continuum containing a fluid and a solid domain 
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.
Figure 2: Sequence of steps to update a “cloud” of nodes representing a domain containing a fluid and a solid part from time () to time () 
A typical solution with the PFEM involves the following steps.
Figure 3 shows another conceptual example of application of the PFEM to modelling the melting and dripping of a polymer object under a heat source acting at a boundary.
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 3 can be also used to explain the application of the PFEM to rock cutting problems. In those cases 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 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 [16,36]. Figure 4a shows the initial grid of fournoded rectangles discretizing the fluid domain and the solid walls. Figures 4b and 4c show the deformed mesh at two later times.
The key equations to be solved in the incompressible thermal flow problem, written in the Lagrangian frame of reference, are the following:
Momentum

(1) 
Mass balance

(2) 
Heat transport

(3) 
In above equations is the velocity along the ith global (cartesian) axis, is the temperature, , and are the density (assumed constant), the specific heat and the conductivity of the material, respectively, and are the body forces and the heat source per unit mass, respectively and are the (Cauchy) stresses related to the velocities by the standard constitutive equation (for incompressible Newtonian material)

(4.a) 

(4.b) 
In Eqs.(4), is the deviatoric stresses, is the pressure (assumed to be positive in compression), is the rate of deformation, is the viscosity and is the Kronecker delta. In the following we will assume the viscosity to be a known function of temperature, i.e .
Indexes in Eqs.(1)(4) range from , where is the number of space dimensions of the problem (i.e. for twodimensional problems). The standard sum notation for repeated indices is assumed, unless otherwise specified.
Eqs.(1)(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 [2,9].
We note that Eqs.(1)(3) are the standard ones for modeling the deformation of viscoplastic materials using the so called “flow approach” [49,50,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.(4.b).
A key problem in the numerical solution of Eqs.(1)(4) is the satisfaction of the incompressibility condition (Eq.(2)). A number of procedures to solve his problem exist in the finite element literature [9,49]. In our approach we use a stabilized formulation based in the socalled finite calculus procedure [27–29,36,38,40]. The essence of this method is the solution of a modified mass balance equation which is written as

(5) 
where is a stabilization parameter given by [12]

(6) 
In the above, is a characteristic length of each finite element (such as for 2D elements) and is the modulus of the velocity vector. In Eq.(5) are auxiliary pressuregradient projection variables chosen so as to ensure that the second term in Eq.(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 variables is completed by adding the following constraint equation to the set of governing equation [36,40]

(7) 
where are arbitrary weighting functions.
The rest of the integral equations are obtained by applying the standard Galerkin technique to the governing equations (1), (2), (3), (5) and (7) and the corresponding boundary conditions [36,40].
We interpolate next in the standard finite element fashion the set of problem variables. For 3D problems these are the three velocities , the pressure , the temperature and the three pressure gradient projections . In our work we use equal order linear interpolation for all variables over meshes of 3noded triangles (in 2D) and 4noded tetrahedra (in 3D) [36,40,53]. The resulting set of discretized equations has the following form
Momentum

(8) 
Mass balance

(9) 
Pressure gradient projection

(10) 
Heat transport

(11) 
In Eqs.(8)–(11) denotes nodal variables and . The different matrices and vectors are given in the Appendix.
The solution in time of Eqs.(8)–(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 is presented in Box I. denotes the values of the nodal variables at time and the iterations. We note the coupling of the flow and thermal equations via the dependence of the viscosity with the temperature.
Figure 5 shows a typical domain with external boundaries and where the velocity and the surface tractions are prescribed, respectively. The domain is formed by fluid () and solid () subdomains (i.e. ). Both subdomains interact at a common boundary 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.
Figure 5: Split of the analysis domain into fluid and solid subdomains. Equality of surface tractions and kinematic variables at the common interface 
Let us define and the set of variables defining the kinematics and the stressstrain fields at the solid and fluid domains at time , respectively, i.e.

(12) 

(13) 
where is the nodal coordinate vector, , and are the vector of displacements, velocities and accelerations, respectively, and are the strain vector, the strainrate (or rate of deformation) vectors and the Cauchy stress vector, respectively, is the temperature and subscripts and denote the variables in the fluid and solid domains, respectively. In the discretized problem, a bar over these variables denotes nodal values.
The coupled fluidstructure interaction (FSI) problem of Figure 4 is solved in this work using the following strongly coupled staggered scheme:
The variables at the solid domain are found via the integration of the equations of dynamic motion in the solid written as [52]

(14) 
where is the vector of nodal accelerations and and 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 nonlinearity using standard nonlinear solid mechanics procedures [52]. The time integration of Eq.(14) is performed using a standard Newmark method.
Solve for the variables at the fluid domain at time () under prescribed surface tractions at the external boundary and prescribed velocities at the external and internal boundaries and , 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 II.
Box II:Staggered solution scheme for the FSI problem (Figure 5). : variables in the solid domain. : variables in the fluid domain 
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) [13,15,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 , where is the total number of nodes in the mesh (Figure 6). The 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 interpolation has been chosen [13,15,16].
Figure 6: Generation of non standard meshes combining different polygons (in 2D) and polyhedra (in 3D) using the extended Delaunay technique. 
Figure 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 [40]. The figure shows the CPU time in seconds for each time step of the algorithm of Section 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 of the total CPU time for each time step, while the solution of the equations and the assembly of the system consume approximately and 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.
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 
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 distribution, where 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 , are considered as boundary nodes. In practice is a parameter close to, but greater than one. Values of 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 [10]. Figure 8 shows an example of the boundary recognition using the Alpha Shape technique.
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 2 and 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.
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 , 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 9). This simple way to treat the fluidwall contact at mesh generation level is a distinct and attractive feature of the PFEM.
Figure 9: Automatic treatment of contact conditions at the fluidwall interface 
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 () between two solid boundaries. If the distance exceeds the minimum value () then the generated elements are treated as fluid elements. Otherwise the elements are treated as contact elements where a relationship between the tangential and normal forces and the corresponding displacement is introduced so as to model elastic and frictional contact effects in the normal and tangential directions, respectively (Figure 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 11–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 11 and 14 show explicitely the layer of contact elements which controls the accuracy of the contact algorithm.
Figure 10: Contact conditions at a solidsolid interface 
Figure 11: Bumping of a ball within a container. The layer of contact elements is shown 
Figure 12: Failure of an arch formed by stone blocks under seismic loading 
Figure 13: Motion of five tetrapods on an inclined plane 
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 [22,48]. The effect of water velocity on soil erosion was studied in [42]. In a recent work we have proposed an extension of the PFEM to model bed erosion [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 [1,32,43].
The algorithm for modeling the erosion of soil/rock particles at the fluid bed is the following:

(15.a) 
with

(15.b) 
where is the modulus of the tangential velocity at the node and is a prescribed distance along the normal of the bed node . Typically is of the order of magnitude of the smallest fluid element adjacent to node (Figure 15).

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

(17) 
Figure 14: Detail of five tetrapods on an inclined plane. The layer of elements modeling the frictional contact conditions is shown 
Figure 15: Modeling of bed erosion by dragging of bed material 
Figure 15 shows an schematic view of the bed erosion algorithm proposed.
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). Material is removed from the excavation front or the tool surface when the work of the frictional forces at the rock/soiltool interface exceeds a prescribed value. A new boundary is defined with the volume that remains in the analysis domain using the alphashape approach as it is typical in the PFEM (Section 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 Archardtype equation [1,5,43] as:

(18) 
where is the volume loss of the material along the contact surface due to wear, is the sliding distance, is the normal force vector to the contact surface and is the hardness of the material. Constant is a nondimensional wear coefficient which depends on the relative contribution of the body under abrasion, adhesion and wear processes [5,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:

(19) 
where all variables are nodal variables, is the relative tangent velocity between the contact surfaces and 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 16.
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:


where is the Jacobian matrix. Typically

where is a parameter of the Newmark scheme [52], is the tangent stiffness matrix of the solid mechanics problem accounting for material and nonlinear geometrical effects [5,6,52].


where denotes values at each particle. Note that the strain and stress history is stored at the particles.
A detailed description of above algorithm, together with many applications, can be found in [5,6].
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 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 18).
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 
Figure 18: Detail of element sizes during the motion of a rigid cylinder within a water container 
Figure 19 shows an example of a wave breaking within a prismatic container including a vertical cylinder. Figure 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.
Figure 19: Evolution of a water column within a prismatic container including a vertical cylinder 
Figure 20: Impact of a wave on a prismatic column on a slab sustained by four pillars. 
Figure 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.
Figure 21: Dragging of a cubic object by a water stream. 
Figure 22 shows the 3D simulation of the interaction of a wave with a vertical pier formed by a collection of reinforced concrete cylinders.
Figure 23 shows the simulation of the falling of two tetrapods in a water container. Figure 24 shows the motion of a collection of ten tetrapods placed in the slope of a breakwaters under an incident wave.
Figure 25 shows a detail of the complex threedimensional interactions between water particles and tetrapods and between the tetrapods themselves.
Figures 26 and 27 show the analysis of the effect of breaking waves on two different sites of a breakwater containing reinforced concrete blocks (each one of mts). The figures correspond to the study of Langosteira harbour in A Coruña, Spain using PFEM.
Figure 28 displays the effect of an overtopping wave on a truck circulating by the perimetral road of the harbour adjacent to the breakwater.
Figure 22: Interaction of a wave with a vertical pier formed by reinforced concrete cylinders. 
Figure 23: Motion of two tetrapods falling in a water container. 
Figure 24: Motion of ten tetrapods on a slope under an incident wave. 
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. 
Figure 26: Effect of breaking waves on a breakwater slope containing reinforced concrete blocks. Detail of the mesh of 4noded tetrahedra near the slope at two different times 
Figure 27: Study of breaking waves on the edge of a breakwater structure formed by reinforced concrete blocks 
Figure 28: Effect of an overtopping wave on a truck passing by the perimetral road of a harbour adjacent to the breakwater 
Figure 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 [39].
Figure 29: Erosion of a 3D earth dam due to an overspill stream. 
Figure 30: Erosion, transport and deposition of particles at a river bed due to a jet stream. 
Figure 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 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 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.
Figure 31: Erosion of unprotected part of a breakwater slope due to sea waves. 
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. 
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 quasirigid 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 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 [4].
Figure 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 viscositytemperature curve that implicitly accounts for molecular weight changes.
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 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 33, but the initial dimensions of the sample are reduced to cm. The initial size of the model is 22475 nodes and 97600 fournoded tetrahedra. The shape of the surface and temperature field at different times after heating begins are shown in Figure 35.
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 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, indepth 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 [41].
Figure 35: Simulation of a 3D polymer melt problem with the PFEM. Melt flow from a heated prismatic sample at different times. 
Figure 36: Melt flow of a heated triangular object into a catch pan. 
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 thermalflow problem are extended with a combustion model governing the burning of combustible and the heat interchanges between the object and the air during combustion [21,24,46]. Figure 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.
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 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 39 and 40 show a similar examples of excavation of a soft soil mass with rotating discs.
Figure 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.
Figure 38: Simulation of a disc excavating a soft wall with the PFEM 
Figure 39: Example of application of the PFEM to the excavation of a soft soil mass with a rotating disc 
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 
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 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 alphashape technique are crucial for capturing the fast and drastic changes of the domain boundary.
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 43–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 43. Figure 44 shows an overview of the simulation as the tunneling process advance and the stress contour lines and Figure 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.
Figure 43:Simulation of a tunneling process with a TBM using the PFEM. Discretization of soil mass and TBM geometry with 4noded tetrahedra 
Figure 44: Simulation of a tunneling operation with a TBM using the PFEM 
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.
The particle finite element method (PFEM) is a powerful computational technique for solving coupled problems in engineering, involving fluidstructure 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 fluidsolid and solidsolid 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 [23].
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.
[1] J.F. Archard, Contact and rubbing of flat surfaces, J. Appl. Phys. 24(8) (1953) 981–988.
[2] R. Aubry, S.R. Idelsohn, E. Oñate, Particle finite element method in fluid mechanics including thermal convectiondiffusion, Computer & Structures 83(1718) (2005) 1459–1475.
[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–948.
[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) 35 September.
[5] J.M. Carbonell, Modeling of ground excavation with the Particle Finite Element method. Ph.D. Thesis, Technical University of Catalonia (UPC), Barcelona, (2009).
[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).
[7] R. Codina, O.C. Zienkiewicz, CBS versus GLS stabilization of the incompressible NavierStokes equations and the role of the time step as stabilization parameter, Communications in Numerical Methods in Engineering (2002) 18(2) (2002) 99–112.
[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 freesurface flows and fluidobject interactions. Computers & Fluids 36 (2007) 27–38.
[9] J. Donea, A. Huerta, Finite element method for flow problems, J. Wiley, (2003).
[10] H. Edelsbrunner, E.P. Mucke, Three dimensional alpha shapes, ACM Trans. Graphics 13 (1999) 43–72.
[11] J. García, E. Oñate, An unstructured finite element solver for ship hydrodynamic problems, J. Appl. Mech. 70 (2003) 18–26.
[11] S.R. Idelsohn, E. Oñate, F. Del Pin, N. Calvo, Lagrangian formulation: the only way to solve some freesurface fluid mechanics problems, Fith World Congress on Computational Mechanics, H.A. Mang, F.G. Rammerstorfer, J. Eberhardsteiner (eds), July 7–12, Viena, Austria, (2002).
[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–912.
[14] S.R. Idelsohn, E. Oñate, F. Del Pin, A lagrangian meshless finite element method applied to fluidstructure interaction problems, Computer and Structures 81 (2003b) 655–671.
[15] S.R. Idelsohn, N. Calvo, E. Oñate, Polyhedrization of an arbitrary point set, Comput. Method Appl. Mech. Engng. 192(2224) (2003c) 2649–2668.
[16] S.R. Idelsohn, E. Oñate, F. Del Pin, The particle finite element method: a powerful tool to solve incompressible flows with freesurfaces and breaking waves, Int. J. Num. Meth. Engng. 61 (2004) 964989.
[17] S.R. Idelsohn, E. Oñate, F. Del Pin, N. Calvo, Fluidstructure interaction using the particle finite element method, Comput. Meth. Appl. Mech. Engng. 195 (2006) 21002113.
[18] S.R. Idelsohn, J. Marti, A. Limache, E. Oñate, Unified Lagrangian formulation for elastic solids and incompressible fluids: Application to fluidstructure interaction problems via the PFEM. Comput Methods Appl Mech Engrg. (2008) 197 1762–1776.
[19] S.R. Idelsohn, M. MierTorrecilla, E. Oñate, Multifluid flows with the Particle Finite Element Method. Comput Methods Appl Mech Engrg. 198 (2009) 2750–2767.
[20] S.R. Idelsohn, M. MierTorrecilla, 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–124.
[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, 57 July 2010, Nottingham, UK.
[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–183.
[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–425.
[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.
[25] M. de Mier Torrecilla, Numerical Simulation of MultiFluid Flows with the Particle Finite Element Method. Ph.D. Thesis, Technical University of Catalonia (UPC), July 2010.
[26] R. Ohayon, Fluidstructure interaction problem, in: E. Stein, R. de Borst, T.J.R. Hugues (Eds.), Enciclopedia of Computatinal Mechanics, Vol. 2, (J. Wiley, 2004) 683–694.
[27] E. Oñate, Derivation of stabilized equations for advectivediffusive transport and fluid flow problems, Comput. Meth. Appl. Mech. Engng. 151 (1998) 233–267.
[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–2) (2000) 355–370.
[29] E. Oñate, Possibilities of finite calculus in computational mechanics Int. J. Num. Meth. Engng. 60(1) (2004) 255–281.
[30] E. Oñate, S.R. Idelsohn, A mesh free finite point method for advectivediffusive transport and fluid flow problems, Computational Mechanics 21 (1998) 283–292.
[31] E. Oñate, J. García, A finite element method for fluidstructure interaction with surface waves using a finite calculus formulation, Comput. Meth. Appl. Mech. Engrg. 191 (2001) 635–660.
[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–3128.
[33] E. Oñate, C. Sacco, S.R. Idelsohn, A finite point method for incompressible flow problems, Comput. Visual. in Science 2 (2000) 67–75.
[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).
[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–610.
[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) 267307.
[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 (45) (2006a) 440455.
[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 (2324) (2006b) 30013037.
[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) 237252.
[40] E. Oñate, S.R. Idelsohn, M.A. Celigueta, R. Rossi, Advances in the particle finite element method for the analysis of fluidmultibody interaction and bed erosion in free surface flows, Comput. Meth. Appl. Mech. Engng. 197 (1920) (2008) 17771800.
[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) 10461072.
[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–178.
[43] Rabinowicz, E., Friction and Wear of materials. Wiley, (1995).
[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–210.
[45] P.B. Ryzhakov, R. Rossi, S. Idelsohn, E.Oñate, A monolithic Lagrangian approach for fluidstructure interaction problems. Journal of Computational Mechanics 46 (6) (2010) 883–899.
[46] P.B. Ryzhakov, R. Rossi, E. Oñate, An algorithm for polymer melting simulation. Conference Proceedings METNUM2009, Barcelona, Spain, 29 June 02 July (2009).
[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–578.
[48] C.F. Wan, R. Fell, Investigation of erosion of soils in embankment dams, J. Geotechnical and Geoenvironmental Engineering 130 (2004) 373–380.
[49] O.C. Zienkiewicz, R.L. Taylor, P. Nithiarasu, The finite element method for fluid dynamics, Elsevier, (2006).
[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–38.
[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–1514.
[52] O.C. Zienkiewicz, R.L. Taylor, The finite element method for solid and structural mechanics, Elsevier, (2005).
[53] O.C. Zienkiewicz, R.L. Taylor, J.Z. Zhu, The finite element method. Its basis and fundamentals, Elsevier, (2005).
The matrices and vectors in Eqs.(8)(11) for a 4noded tetrahedron are:








In the above equations indexes run from 1 to the number of element nodes (4 for a tetrahedron), is the heat flow prescribed at the external boundary , t is the surface traction vector and and are the element volume and the element boundary, respectively.
Published on 22/05/19
DOI: 10.1007/9789400707351_1
Licence: CC BYNCSA license
Are you one of the authors of this document?