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 [9,49] and the references there included. A survey of recent works in fluid-structure 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 ${\textstyle div}$-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 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 [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 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 ${\textstyle t}$. The new set of variables in both domains are sought for in the next or updated configuration at time ${\textstyle t+\Delta t}$ (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

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.

Figure 2: Sequence of steps to update a “cloud” of nodes representing a domain containing a fluid and a solid part from time $n$ ($t=t_{n}$) to time $n+2$ ($t=t_{n}+2\Delta t$)

A typical solution with the PFEM involves the following steps.

1. The starting point at each time step is the cloud of points in the fluid and solid domains. For instance ${\textstyle {}^{n}C}$ denotes the cloud at time ${\textstyle t=t_{n}}$ (Figure 2).
2. Identify the boundaries for both the fluid and solid domains defining the analysis domain ${\textstyle {}^{n}V}$ 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 [10] is used for the boundary definition (Section 5).
3. Discretize the fluid and solid domains with a finite element mesh ${\textstyle {}^{n}M}$. In our work we use an innovative mesh generation scheme based on the extended Delaunay tesselation (Section 4) [13,14,16].
4. 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 ${\textstyle t+\Delta t}$: velocities, pressure, viscous stresses and temperature in the fluid and displacements, stresses, strains and temperature in the solid.
5. Move the mesh nodes to a new position ${\textstyle {}^{n+1}C}$ where ${\textstyle n+1}$ denotes the time ${\textstyle t_{n}+\Delta t}$, in terms of the time increment size. This step is typically a consequence of the solution process of step 4.
6. Go back to step 1 and repeat the solution process for the next time step to obtain ${\textstyle {}^{n+2}C}$. The process is shown in Figure 2.

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 ${\textstyle q}$ acting at a boundary.

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 $q$. Crossed circles denote prescribed nodes at the boundary. The same figure explains the application of the PFEM to modelling a rock cutting problem

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 ${\textstyle q}$ 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 four-noded rectangles discretizing the fluid domain and the solid walls. Figures 4b and 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

$\rho {\partial v_{i} \over \partial t}={\partial \sigma _{ij} \over \partial x_{j}}+b_{i}\qquad {\hbox{in }}\Omega$

(1)

Mass balance

${\partial v_{i} \over \partial x_{i}}=0\qquad {\hbox{in }}\Omega$

(2)

Heat transport

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

(3)

In above equations ${\textstyle v_{i}}$ is the velocity along the ith global (cartesian) axis, ${\textstyle T}$ is the temperature, ${\textstyle \rho }$, ${\textstyle c}$ and ${\textstyle k_{i}}$ are the density (assumed constant), the specific heat and the conductivity of the material, respectively, ${\textstyle b_{i}}$ and ${\textstyle Q}$ are the body forces and the heat source per unit mass, respectively and ${\textstyle \sigma _{ij}}$ are the (Cauchy) stresses related to the velocities by the standard constitutive equation (for incompressible Newtonian material)

$\sigma _{ij}=s_{ij}-p\delta _{ij}$

(4.a)

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

(4.b)

In Eqs.(4), ${\textstyle s_{ij}}$ is the deviatoric stresses, ${\textstyle p}$ is the pressure (assumed to be positive in compression), ${\textstyle {\dot {\varepsilon }}_{ij}}$ is the rate of deformation, ${\textstyle \mu }$ is the viscosity and ${\textstyle \delta _{ij}}$ is the Kronecker delta. In the following we will assume the viscosity ${\textstyle \mu }$ to be a known function of temperature, i.e ${\textstyle \mu =\mu (T)}$.

Indexes in Eqs.(1)-(4) range from ${\textstyle i,j=1,n_{d}}$, where ${\textstyle n_{d}}$ is the number of space dimensions of the problem (i.e. ${\textstyle n_{d}=2}$ for two-dimensional 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).

3.2 Discretization of the equations

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

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

(5)

where ${\textstyle \tau }$ is a stabilization parameter given by [12]

$\tau =\left({\frac {2\rho \vert \mathbf {v} \vert }{h}}+{\frac {8\mu }{3h^{2}}}\right)^{-1}$

(6)

In the above, ${\textstyle h}$ is a characteristic length of each finite element (such as ${\textstyle [A^{(e)}]^{1/2}}$ for 2D elements) and ${\textstyle \vert \mathbf {v} \vert }$ is the modulus of the velocity vector. In Eq.(5) ${\textstyle \pi _{i}}$ are auxiliary pressure-gradient 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 ${\textstyle \pi _{i}}$ variables is completed by adding the following constraint equation to the set of governing equation [36,40]

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

(7)

where ${\textstyle w_{i}}$ 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 ${\textstyle v_{i}}$, the pressure ${\textstyle p}$, the temperature ${\textstyle T}$ and the three pressure gradient projections ${\textstyle \pi _{i}}$. 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) [36,40,53]. The resulting set of discretized equations has the following form

Momentum

$\mathbf {M} {\dot {\bar {\boldsymbol {v}}}}+\mathbf {K} (\mu ){\bar {\boldsymbol {v}}}-\mathbf {G} {\bar {\boldsymbol {p}}}={\boldsymbol {f}}$

(8)

Mass balance

$\mathbf {G} ^{T}{\bar {\boldsymbol {v}}}+\mathbf {L} {\bar {\boldsymbol {p}}}+\mathbf {Q} {\bar {\boldsymbol {\pi }}}=\mathbf {0}$

(9)

Pressure gradient projection

${\hat {\boldsymbol {M}}}{\bar {\boldsymbol {\pi }}}+\mathbf {Q} ^{T}{\bar {\boldsymbol {p}}}=\mathbf {0}$

(10)

Heat transport

$\mathbf {C} {\dot {\bar {\boldsymbol {T}}}}+\mathbf {H} {\bar {\boldsymbol {T}}}={\boldsymbol {q}}$

(11)

In Eqs.(8)–(11) ${\textstyle {\bar {(\cdot )}}}$ denotes nodal variables and ${\textstyle {\dot {\bar {(\cdot )}}}={\partial \over \partial t}{\bar {(\cdot )}}}$. 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. ${\textstyle ^{t+\Delta t}({\bar {\boldsymbol {a}}})^{j+1}}$ denotes the values of the nodal variables ${\textstyle {\bar {\boldsymbol {a}}}}$ at time ${\textstyle t+\Delta t}$ and the ${\textstyle j+1}$ iterations. We note the coupling of the flow and thermal equations via the dependence of the viscosity ${\textstyle \mu }$ with the temperature.

4 Overview of the coupled FSI algoritm

Figure 5 shows a typical domain ${\textstyle V}$ with external boundaries ${\textstyle \Gamma _{V}}$ and ${\textstyle \Gamma _{t}}$ where the velocity and the surface tractions are prescribed, respectively. The domain ${\textstyle V}$ is formed by fluid (${\textstyle V_{F}}$) and solid (${\textstyle V_{S}}$) subdomains (i.e. ${\textstyle V=V_{F}\cup V_{S}}$). Both subdomains interact at a common boundary ${\textstyle \Gamma _{FS}}$ 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 $V$ into fluid and solid subdomains. Equality of surface tractions and kinematic variables at the common interface

Let us define ${\textstyle {}^{t}S}$ and ${\textstyle {}^{t}F}$ the set of variables defining the kinematics and the stress-strain fields at the solid and fluid domains at time ${\textstyle t}$, respectively, i.e.

${}^{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}$

(12)

${}^{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}$

(13)

where ${\textstyle {\boldsymbol {x}}}$ is the nodal coordinate vector, ${\textstyle {\boldsymbol {u}}}$, ${\textstyle {\boldsymbol {v}}}$ and ${\textstyle {\boldsymbol {a}}}$ are the vector of displacements, velocities and accelerations, respectively, ${\textstyle {\boldsymbol {\varepsilon }},{\dot {\boldsymbol {\varepsilon }}}}$ and ${\textstyle {\boldsymbol {\sigma }}}$ are the strain vector, the strain-rate (or rate of deformation) vectors and the Cauchy stress vector, respectively, ${\textstyle T}$ is the temperature and subscripts ${\textstyle F}$ and ${\textstyle S}$ 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 4 is solved in this work using the following strongly coupled staggered scheme:

1. We assume that the variables in the solid and fluid domains at time ${\textstyle t}$ (${\textstyle {}^{t}S}$ and ${\textstyle {}^{t}F}$) are known.
2. Solve for the variables at the solid domain at time ${\textstyle t+\Delta t}$ (${\textstyle {}^{t+\Delta t}{S}}$) under prescribed surface tractions at the fluid-solid boundary ${\textstyle \Gamma _{FS}}$. The boundary conditions at the part of the external boundary intersecting the domain are the standard ones in solid mechanics.

The variables at the solid domain ${\textstyle {}^{t+\Delta t}{S}}$ are found via the integration of the equations of dynamic motion in the solid written as [52]

${\boldsymbol {M}}_{s}{\bar {\boldsymbol {a}}}_{s}+{\boldsymbol {g}}_{s}-{\boldsymbol {f}}_{s}={\boldsymbol {0}}$

(14)

where ${\textstyle {\bar {\boldsymbol {a}}}_{s}}$ is the vector of nodal accelerations and ${\textstyle {\boldsymbol {M}}_{s},{\boldsymbol {g}}_{s}}$ and ${\textstyle {\boldsymbol {f}}_{s}}$ 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 [52]. The time integration of Eq.(14) is performed using a standard Newmark method.

Solve for the variables at the fluid domain at time ${\textstyle t+\Delta t}$ (${\textstyle {}^{t+\Delta t}{F}}$) under prescribed surface tractions at the external boundary ${\textstyle \Gamma _{t}}$ and prescribed velocities at the external and internal boundaries ${\textstyle \Gamma _{V}}$ and ${\textstyle \Gamma _{FS}}$, 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). ${\textstyle S}$: variables in the solid domain. ${\textstyle F}$: 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) [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 ${\textstyle n}$, where ${\textstyle n}$ is the total number of nodes in the mesh (Figure 6). The ${\textstyle C^{\circ }}$ 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 ${\textstyle C^{\circ }}$ 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 ${\textstyle 20\%}$ of the total CPU time for each time step, while the solution of the equations and the assembly of the system consume approximately ${\textstyle 50\%}$ and ${\textstyle 20\%}$ of the CPU time for each time step, respectively. These figures prove that the generation of the mesh has an acceptable cost in the PFEM.

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 ${\textstyle h(x)}$ distribution, where ${\textstyle h(x)}$ 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 $\alpha h$, are considered as boundary nodes. In practice ${\textstyle \alpha }$ is a parameter close to, but greater than one. Values of ${\textstyle \alpha }$ 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.

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 ${\textstyle \Gamma _{FS}}$, 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 fluid-wall contact at mesh generation level is a distinct and attractive feature of the PFEM.

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 (${\textstyle h_{c}}$) between two solid boundaries. If the distance exceeds the minimum value (${\textstyle h_{c}}$) 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 solid-solid 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

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

1. Compute at every point of the bed surface the resultant tangential stress ${\textstyle {\hat {\tau }}}$ induced by the fluid motion. In 3D problems ${\textstyle {\hat {\tau }}=(\tau _{s}^{2}+\tau _{t})^{2}}$ where ${\textstyle \tau _{s}}$ and ${\textstyle \tau _{t}}$ are the tangential stresses in the plane defined by the normal direction ${\textstyle {\boldsymbol {n}}}$ at the bed node. The value of ${\textstyle {\hat {\tau }}}$ for 2D problems can be estimated as follows:

   $\tau _{t}=\mu \gamma _{t}$

   (15.a)

with

$\gamma _{t}={1 \over 2}{\partial v_{t} \over \partial n}={v_{t}^{k} \over 2h_{k}}$

(15.b)

where ${\textstyle v_{t}^{k}}$ is the modulus of the tangential velocity at the node ${\textstyle k}$ and ${\textstyle h_{k}}$ is a prescribed distance along the normal of the bed node ${\textstyle k}$. Typically ${\textstyle h_{k}}$ is of the order of magnitude of the smallest fluid element adjacent to node ${\textstyle k}$ (Figure 15).

2. Compute the frictional work originated by the tangential stresses at the bed surface as

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

   (16)

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

${}^{n}W_{f}={}^{n-1}W_{f}+\tau _{t}\gamma _{t}\,\Delta t$

(17)

3. The onset of erosion at a bed point occurs when ${\textstyle {}^{n}W_{f}}$ exceeds a critical threshold value ${\textstyle W_{c}}$ defined empirically according to the specific properties of the bed material.
4. If ${\textstyle {}^{n}W_{f}>W_{c}}$ 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.
5. 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.

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.

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). 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). 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 [1,5,43] as:

$V_{w}=K\;{\frac {\|\mathbf {f} _{n}\|}{H}}\;s$

(18)

where ${\textstyle V_{w}}$ is the volume loss of the material along the contact surface due to wear, ${\textstyle s\,(m)}$ is the sliding distance, ${\textstyle \mathbf {f} _{n}}$ is the normal force vector to the contact surface and ${\textstyle H}$ is the hardness of the material. Constant ${\textstyle K}$ is a non-dimensional 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:

$V_{w}^{t+\triangle t}=V_{w}^{t}+K\;{\frac {\|\mathbf {f} _{n}\|}{H}}\;(\|\mathbf {v} _{t}\|\cdot \triangle t)$

(19)

where all variables are nodal variables, ${\textstyle \mathbf {v} _{t}}$ is the relative tangent velocity between the contact surfaces and ${\textstyle \triangle t}$ 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:

1. Read initial geometrical, mechanical and kinematic conditions from a reference mesh.
2. Transfer the elemental variables to the particles (i.e. the nodes).
3. For each time step and each Newton iteration:
• Compute internal forces and contact forces at nodes

${\boldsymbol {r}}:={\boldsymbol {M}}{\bar {\boldsymbol {a}}}_{s}+{\boldsymbol {g}}_{s}+{\boldsymbol {f}}_{c}-{\boldsymbol {f}}_{s}$

• Compute displacement increments and update displacement values

$\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}}}$

where ${\textstyle {\boldsymbol {A}}}$ is the Jacobian matrix. Typically

${\boldsymbol {A}}={\frac {1}{\beta \Delta t^{2}}}{\boldsymbol {M}}+{\boldsymbol {K}}_{T}+{\boldsymbol {K}}_{c}$

where ${\textstyle \beta }$ is a parameter of the Newmark scheme [52], ${\textstyle {\boldsymbol {K}}_{T}}$ is the tangent stiffness matrix of the solid mechanics problem accounting for material and non-linear geometrical effects [5,6,52].

4. Compute internal variables, strains and stresses at integration points within each element.
5. Check convergence of Newton iterations.
6. Once the iterative solutions has converged
• Update particle positions: ${\textstyle {}^{t+\Delta t}{\boldsymbol {x}}={}^{t}{\boldsymbol {x}}+{}^{t+\Delta t}\Delta {\bar {\boldsymbol {u}}}}$
• Compute velocities (${\textstyle {}^{t+\Delta t}{\bar {\boldsymbol {v}}}}$) and accelerations at particles (${\textstyle {}^{t+\Delta t}{\bar {\boldsymbol {a}}}}$).
• Transfer strains and stresses from elements to particles:

${}^{t+\Delta t}{\boldsymbol {\sigma }}_{p}={}^{t}{\boldsymbol {\sigma }}_{p}+\Delta {\boldsymbol {\sigma }}_{p}$

${}^{t+\Delta t}{\boldsymbol {\varepsilon }}_{p}={}^{t}{\boldsymbol {\varepsilon }}_{p}+\Delta {\boldsymbol {\varepsilon }}_{p}$

where ${\textstyle (\cdot )_{p}}$ denotes values at each particle. Note that the strain and stress history is stored at the particles.

• Update constitutive law parameters.
7. Check damage and erosion (wear) on particles. Remove eroded particles from the excavation front and worn particles from cutting tools.
8. Boundary recognition via the alpha shape method. Create new mesh. Update problem dimensions if the number of particles has changed.
9. Identify interface elements for contact.
10. Initiate solution for next time step.

A detailed description of above algorithm, together with many applications, can be found in [5,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 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

10.2 Impact of water streams on rigid structures

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.

10.3 Dragging of objects by water streams

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.

10.4 Impact of sea waves on piers and breakwaters

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 three-dimensional 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 ${\textstyle 4\times 4}$ 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 4-noded 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

10.5 Soil erosion

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.

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 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 viscosity-temperature curve that implicitly accounts for molecular weight changes.

Figure 33: Polymer melt experiment. Viscosity vs. temperature for PP702N polypropylene in its initial undegraded form and after exposure to 30 kW/m$^{2}$ and 40 kW/m$^{2}$ 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 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 ${\textstyle 10\times 2.5\times 2.5}$ 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 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, 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 [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 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 [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.

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 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 alpha-shape 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 4-noded 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.

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

[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 convection-diffusion, Computer & Structures 83(17-18) (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) 3-5 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 Navier-Stokes 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 free-surface flows and fluid-object 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 free-surface 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 fluid-structure 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(22-24) (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 free-surfaces and breaking waves, Int. J. Num. Meth. Engng. 61 (2004) 964-989.

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

[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–1776.

[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–2767.

[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–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, 5-7 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 Multi-Fluid Flows with the Particle Finite Element Method. Ph.D. Thesis, Technical University of Catalonia (UPC), July 2010.

[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–694.

[27] E. Oñate, Derivation of stabilized equations for advective-diffusive 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 advective-diffusive transport and fluid flow problems, Computational Mechanics 21 (1998) 283–292.

[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–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) 267-307.

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

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

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

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

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

[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 fluid-structure 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 METNUM-2009, 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).

APPENDIX

The matrices and vectors in Eqs.(8)-(11) for a 4-noded tetrahedron are:

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

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

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

${\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}$

${\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]$

${\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}}$

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

$[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_{i}QdV-\int _{\Gamma _{q}^{e}}N_{i}q_{n}d\Gamma$

In the above equations indexes ${\textstyle i,j}$ run from 1 to the number of element nodes (4 for a tetrahedron), ${\textstyle q_{n}}$ is the heat flow prescribed at the external boundary ${\textstyle \Gamma _{q}}$, t is the surface traction vector ${\textstyle \mathbf {t} =[t_{x},t_{y},t_{z}]^{T}}$ and ${\textstyle V^{e}}$ and ${\textstyle \Gamma ^{e}}$ are the element volume and the element boundary, respectively.