A State of the Art Review of the Particle Finite Element Method (PFEM)

Abstract

The particle finite element method (PFEM) is a powerful and robust numerical tool for the simulation of multi-physics problems in evolving domains. The PFEM exploits the Lagrangian framework to automatically identify and follow interfaces between different materials (e.g. fluid–fluid, fluid–solid or free surfaces). The method solves the governing equations with the standard finite element method and overcomes mesh distortion issues using a fast and efficient remeshing procedure. The flexibility and robustness of the method together with its capability for dealing with large topological variations of the computational domains, explain its success for solving a wide range of industrial and engineering problems. This paper provides an extended overview of the theory and applications of the method, giving the tools required to understand the PFEM from its basic ideas to the more advanced applications. Moreover, this work aims to confirm the flexibility and robustness of the PFEM for a broad range of engineering applications. Furthermore, presenting the advantages and disadvantages of the method, this overview can be the starting point for improvements of PFEM technology and for widening its application fields.

1 Introduction

The last decades have seen a growing interest in the development of computational methods for the simulation of engineering problems. A robust and efficient numerical simulation is particularly complex in the presence of multi-physics phenomena and/or large deformations of the physical domains. Typical examples can be found in unsteady free-surface fluid dynamics problems, fluid–structure interaction applications with large motions of fluid–solid interfaces, non-linear solid mechanics with large changes of the topology and contact of solid bodies, and thermal-mechanical coupled analysis in the presence of phase-change phenomena.

To tackle these complex problems, the Finite Element Method (FEM) has been generally privileged. In order to solve a problem in mechanics with the FEM, the reference configuration¹ should be provided with a mesh. Depending on the framework considered, different FEM approaches arise.

For continuum mechanics problems, in a Eulerian approach, the finite element mesh is fixed and the material moves across the grid, being the mesh nodes dissociated from physical particles. Due to the relative motion between the material and the grid, convective terms appear in the definition of the time derivatives. Eulerian meshes are particularly suited for large deformation problems in enclosed domains, as those generally considered in standard Computational Fluid Dynamics (CFD). On the other hand, they do not provide a natural definition of evolving interfaces (like a free surface in fluid flows), calling for ad-hoc techniques such as level set [97] or volume of fluid [42] approaches.

On the contrary, in a Lagrangian approach, the finite element mesh moves along with the continuum body. Consequently, boundaries and interfaces are naturally tracked during the motion allowing for a simpler imposition of boundary conditions. As the material points coincide with the grid nodes, no convective terms appear in the governing equations and material derivatives reduce to time derivatives [23]. Also the integration points move with the material, so constitutive laws are evaluated at the same material points for all the duration of the analysis. This feature is particularly useful for materials with history-dependent behaviour. Nevertheless, in large-deformation problems, the mesh can undergo excessive distortion leading to accuracy loss or even to compromise the FEM calculation.

Arbitrary Lagrangian Eulerian (ALE) description gives a possibility to overcome the typical difficulties related to Lagrangian and Eulerian approaches. In the ALE description, the reference configuration is chosen ad-hoc to reduce mesh distortion and, in general, does not coincide either with material or spatial configurations. The mesh is defined in this reference configuration and moves independently from the material motion. ALE strategy exploits some of the advantages of both Lagrangian and Eulerian descriptions, however, the method has important limitations for very large and unpredictable domain deformations. In these cases, also the ALE mesh can become too distorted, reducing or compromising the accuracy of the computation.

In this work, we focus on the use of Lagrangian approaches in the presence of very large deformations of the domain. It is worth mentioning that there exist different examples of Lagrangian mesh-based solvers for the solution of mechanical problems also in the presence of significant deformations. For example, Hassager and Bisgaard [40] proposes a Lagrangian FEM for the solution of non-Newtonian fluid flows. Ramaswamy et al. [102] uses a Lagrangian FEM in conjunction with a fractional step method to solve small-amplitude sloshing problems. In [101], the same research group solves fluid dam-break problems and large-amplitude sloshing, and in [41], a solitary wave propagation. All these Lagrangian methods maintain the same discretization for all the duration of the analysis, leaving the mesh free to move and, eventually, deteriorate. Consequently, these approaches apply only to analyses in which the discretized domain does not undergo very large deformations, thus limiting the applications to real case problems.

Mesh deterioration constituted for a long time the intrinsic limit of Lagrangian mesh-based solvers. In the literature, there exist two different options to overcome this endemic feature of Lagrangian mesh-based solvers: to introduce a remeshing technique or to abandon completely the concept of mesh. The latter option gives rise to the so-called meshless methods. These techniques represent the behaviour of a physical problem by a collection of particles. Each particle has assigned all physical properties and moves according to its weight and the interaction forces with the neighbour particles. Over the last decades, several meshless methods have been proposed, based either on weak or strong forms of the conservation equations. Meshless methods fall outside the scope of this review. An interested reader can refer to e.g. [67].

The second possibility to avoid mesh distortion is a remeshing technique. Here, when the Lagrangian mesh becomes too distorted, a new mesh is created with an ad-hoc procedure. An example of the application of this technique to fluid flow problems can be found in [36], where the authors reconstruct locally triangular meshes to solve fluid dynamics problems with a finite difference approach. Alternatively, in [3] a new finite element mesh is built from scratch whenever the elements become too distorted. Muttin et al. [80] uses a Lagrangian finite element method to simulate metal casting problems with an automatic remeshing technique to avoid mesh distortion. Radovitzky and Ortiz [100] and Malcevic and Ghattas [68] propose a continuous and adaptive remeshing at each time step. All these Lagrangian approaches are based on some form of remeshing techniques. In all cases, when a new mesh is generated, the results need to be remapped from the old to the new mesh. Unavoidably, these operations introduce unwanted numerical diffusion into the numerical solution [36].

An innovative idea to exploit the Lagrangian framework and overcome mesh distortion issues has been proposed by Idelsohn et al. [52]. They introduced the so-called Particle Finite Element Method (PFEM), an innovative numerical tool able to solve complex non-linear problems in significantly evolving domains [52, 87].

The PFEM combines the accuracy and robustness of mesh-based techniques with the advantages of particle-based methods. The PFEM discretizes the physical domain with a mesh on which the differential governing equations are solved with a standard finite element approach. Following a Lagrangian description, the mesh nodes move according to the equations of motion, behaving like particles and transporting their momentum together with all their physical properties. In the PFEM, the mesh distortion issue, typical of Lagrangian mesh-based solvers, is overcome by generating a new mesh when the current one gets too distorted. However, unlike the previously mentioned methods, to avoid remapping from mesh to mesh, the PFEM keeps the nodes of the previous mesh fixed. The new connectivity is built using the Delaunay Tessellation and a specific technique is used to identify internal and external boundaries. The obtained mesh is then used as the support over which the differential equations are solved in a standard FEM fashion.

Although the PFEM was initially conceived for fluid dynamics and fluid–structure interaction problems, the method soon was been extended to the solution of non-linear solid mechanics problems and to different types of multi-physics problems interesting for varied fields of engineering and technology. This paper provides an extended overview of the theory and applications of the method and it contains all the ingredients to understand the PFEM from its basics to the more advanced features. A very extended literature review of the method is also given by highlighting the main advances of each contribution. The new techniques arisen from the original PFEM formulation are also presented.

The paper is structured as follows. Section 2 introduces the underlying and general concepts of the Particle Finite Element Method (PFEM). Section 3 gives an extended review of the application of the PFEM to fluid dynamics problems. In Sect. 4, the use of the PFEM in fluid–structure interaction problems is examined, while in Sect. 5, the application of the method to non-linear solid mechanics is considered. Section 6 is dedicated to the description of PFEM for other coupled problems, namely multi-fluids and thermal–mechanical coupled analysis. In Sect. 7, some significant applications of the method to several engineering and industrial problems are shown, while Sect. 8 shows some of the recent advances in the PFEM formulation. Finally, the concluding remarks of this work are given in Sect. 9.

2 The Particle Finite Element Method (PFEM)

The PFEM is a numerical technique developed for the solution of multi-physics problems involving large deformations of the domain. It was originally conceived to treat free-surface fluid flows [52] and fluid–structure interaction phenomena [53], but later it has been applied to many other physical problems (see next sections for details).

The key idea of the PFEM is to combine a Lagrangian Finite Element Method (FEM) with an efficient and fast remeshing procedure. In the PFEM, the domain is defined by a set of particles (coinciding with the mesh nodes) that move in a Lagrangian manner according to the calculated nodal variables (e.g. velocity or displacements) and bringing their physical properties (e.g. density, viscosity). Unlike meshless approaches, the interacting forces between particles are evaluated using a finite element mesh. In this sense, the PFEM can be seen as both a FEM-based and a particle method.

One of the most characteristic features of the PFEM is the mesh regeneration algorithm. Whenever the Lagrangian motion of the nodes leads to an excessively distorted mesh, such mesh is deleted and a new one is generated over the same set of nodes. To do so, after the elimination of the elements of the previous distorted mesh, a Delaunay triangulation algorithm is used to rebuild the element connectivity, and an alpha shape scheme [25] is used to define the internal and external boundaries (see the next section). In case of data stored at the element integration points (typically, historical variables in solid mechanics problems), Gauss points data must be transferred from the old mesh to the new one (Sect. 5.1).

To summarize, the fundamental features of the PFEM are the following:

• 1. Lagrangian framework for the description of motion.

• 2. Mesh nodes are treated as physical particles.

• 3. All information is stored at the mesh nodes.

• 4. The FEM is used to solve the governing equations.

• 5. Mesh connectivity is regenerated with a Delaunay Tessellation.

• 6. Boundaries are recovered through ad-hoc techniques (e.g. alpha-shape method).

2.1 The PFEM Steps

Excluding the FEM solution of the differential equations, which differ for each specific problem, the PFEM solution algorithm is independent on the physics of the problem. A general solution scheme of the PFEM can be summarized as follows.

• 1. Fill the domain with a set of points referred to as “particles” (Fig. 1a).

• 2. Generate a finite element mesh using the particles as nodes (Fig. 1b).

• 3. Identify the external and internal boundaries of the computational domain (Fig. 1c).

• 4. Solve the Lagrangian form of the governing equations with the FEM.

• 5. Update the positions of the nodes (Fig. 1d).

• 6. Proceed to the next time step. If remesh is needed go to step 2, otherwise, go directly to step 5. Figures 1e–g show the solution step for the following time instant.

Fig. 1

Steps of Particle Finite Element Method

In step 2, the mesh can be regenerated with any tessellation algorithm. Typically, the Delaunay triangulation is used in the PFEM (Sect. 2.2.1). The identification of boundaries in step 3, needed to compute the domain integrals and to impose correctly the boundary conditions, is performed using the alpha Shape method (Sect. 2.2.2).

It is also important to remark that equations of motions solved in step 4, can be non-linear and so they may require an iterative solution scheme. For those formulations using historical variables, e.g. in non-linear solid mechanics (Sect. 5), a remap of the historical variables stored at the Gauss points on the new mesh is needed before step 4.

At the end of each time step, the quality of the mesh is evaluated (step 6). The original works of the PFEM propose to perform the triangulation at every time step. Other papers suggest performing the remeshing only when the mesh is too distorted globally. This second strategy leads to a reduced computational cost, but also to a lower quality of the mesh. Examples of both alternatives can be found in implicit PFEM formulations (see e.g. [28]). On the contrary, for explicit strategies, due to the high computational cost, only the second strategy is suitable [73].

2.2 The Mesh

A key step in the PFEM solution is the generation of the finite element mesh. In large deformation problems, this operation is performed frequently, eventually at each time step for the most critical cases. Therefore, a very fast, efficient and robust algorithm is required. In the PFEM, this is performed using an enhanced Delaunay triangulation algorithm. It is important to highlight that this procedure should be considered a redefinition of the element connectivities rather than a real remeshing because the mesh nodes of the previous mesh are kept in the same position. Due to the crucial importance of this step in the PFEM solution scheme, a detailed description of the basic concepts, of the implementation details and of the implications of the remeshing step are given in this section.

2.2.1 The Delaunay Triangulation

Before defining the Delaunay triangulation, the Voro-noï diagram has to be introduced. Given a set of N points $(n_1,..., n_N)$, the Voronoï diagram is defined as the partition of ${\mathbb{R}}^3$ in convex regions $T_i$ where a node $n_i$ is associated to each region $T_i$, such that every point of $T_i$ is closer to $n_i$ than to any other nodes $n_j$ with $i\ne j$, i.e.

$$\begin{aligned} T_i=\left\{ {\mathbf{x}}\in {\mathbb{R}}^3 : d({\mathbf{x}},{\mathbf{x}}_i) \le d({\mathbf{x}},{\mathbf{x}}_j) \quad \forall j\ne i \right\} \end{aligned}$$

(1)

where $d({\mathbf{x}},{\mathbf{x}}_i):= \Vert {\mathbf{x}}-{\mathbf{x}}_i \Vert$ is the Euclidean norm. Every region $T_i$ is called Voronoï cell. Each Voronoï cell is convex and closed if internal, open if placed at the boundary.

The Delaunay triangulation can be constructed by joining the points whose Voronoï cells have a common boundary. The Delaunay triangulation can be considered the dual of the Voronoï diagram because two nodes of the Delaunay triangulation are joined by an edge, only if the respective Voronoï cells share a boundary. As a consequence, the Delaunay tessellation generates a mesh of tetrahedra (in 3D) and triangles (in 2D). Figure 2 shows a 2D example of points, Voronoï cells (dashed line) and Delaunay triangulation (solid line).

Fig. 2

Example of Voronoï diagram (dashed line) and Delaunay triangulation (solid line)

A fundamental property of the Delaunay triangulation is that none of its vertices lays inside any tetrahedron’s circumsphere (in 3D) or triangle’s circumcircle (in 2D), see the blue circle in Fig. 2. Moreover, the vertices of Voronoï cells represent the center of tetrahedron’s circumsphere (in 3D) or triangle’s circumcircle (in 2D) of the Delaunay triangulation (Fig. 2). Given a set of points in the space, the Voronoï diagram is unique, whereas different Delaunay triangulations may exist.

In 2D, the Delaunay triangulation has remarkable properties such as the minimization of the maximum radius of an element circumcircle and the maximization of the minimum angle among all the elements (max-min property). On the other hand, the 3D Delaunay algorithm loses some of the optimal properties of its 2D counterpart. Unfortunately, this has important consequences on the possible presence of bad quality tetrahedra in the mesh, such as zero-volume elements (slivers).

The mesh given by the union of the generated tetrahedra or triangles is the convex hull of the points, i.e. the convex figure with minimum volume (area in 2D) that encloses all points. Consequently, the Delaunay triangulation can generate only convex domains. In the next section, it is shown how in the PFEM, this geometrical limitation is overcome through the application of the alpha-shape technique. It is worth noting that non-convex domains can be created by Constrained Delaunay tessellation algorithms. However, these methods do not allow the reconnection of different parts of the computational domain and for this reason, cannot be used for several free-surface fluids applications, while they have potential for solid mechanics problems.

The overall good properties of the Delaunay triangulation together with the availability of several fast open-source algorithms, explain the popularity of this tessellation procedure in the PFEM framework. Nevertheless, we remark that using Delaunay meshes is not mandatory for the PFEM solution scheme. The key point is to obtain a mesh very rapidly starting from a points distribution. Hence, also non-Delaunay meshes could be used without altering the nature of the method.

2.2.2 Boundary Recovery through the Alpha-Shape Technique

In a Lagrangian framework, the current volume $\varOmega _t$ and its external boundary $\varGamma _t=\partial \varOmega _t$ are defined by the position of the material points. As introduced in the previous section, the Delaunay triangulation generates a convex figure which encloses all the nodes belonging to the set. As the convex hull may be not conformal with the actual external boundaries of the computational domain, the new contours have to be identified every time the Delaunay triangulation is done. To clarify this problem, a 2D example is presented in Fig. 3: a set of points is shown in Fig. 3a and its Delaunay triangulation in Fig. 3b. It is clear that the Delaunay triangulation does not match the real internal and external boundaries.

Fig. 3

a Distribution of points, b delaunay triangulation

As originally proposed in [54], a possible method to recover the real shape of the point distribution is the so-called alpha-shape method [25]. This technique is based on the observation that the unphysical elements which do not belong to the real domain, are in general the largest and most distorted ones because they connect nodes that are far from each other.

The basic idea is to remove these unnecessary elements from the mesh using a geometrical criterion based on the mesh distortion. For each element e of the mesh, an index of elemental distortion $\alpha _e$ is defined as:

$$\begin{aligned} \alpha _e=\frac{R_e}{h_{mean}} \end{aligned}$$

(2)

where $R_e$ is the radius of the circumsphere (or circumcircle in 2D) to the considered element and $h_{mean}$ is a characteristic mesh size. An example of $h_{mean}$ can be the average of the minimum element side among all the elements of the initial mesh [18].

A threshold value ${\bar{\alpha }}$ for the distortion of the mesh can be fixed and, consequently, all the elements that do not satisfy the condition:

$$\begin{aligned} \alpha _e \le {\bar{\alpha }} \end{aligned}$$

(3)

are removed from the mesh. Figure 4 shows the Delaunay triangulation coupled with the alpha-shape scheme for the example introduced previously.

Fig. 4

Delaunay triangulation with alpha-shape: a ${\bar{\alpha }}=1.9$, b ${\bar{\alpha }}=1.4$

A proper definition of the distortion threshold ${\bar{\alpha }}$ is mandatory for the correct boundary recognition. Different values of ${\bar{\alpha }}$ can lead to different configurations (see e.g. Fig. 4). If ${\bar{\alpha }}$ is too large, some overly distorted or too large elements can be incorporated into the mesh. On the contrary, if ${\bar{\alpha }}$ is too small, too many elements are removed creating unphysical holes within the analysis domain. Obviously, for ${\bar{\alpha }}\rightarrow \infty$, the original Delaunay tessellation is recovered. A critical review of the effect of the parameter ${\bar{\alpha }}$ can be found in [28].

An alternative definition of the alpha-shape method can be found in the first papers on the PFEM (see e.g. [52, 87]). Considering that the particles follow a variable distribution h(x), where h(x) is the minimum distance between two particles, the following criterion has been used:

All particles on an empty sphere with a radius r(x) bigger than ${\bar{\alpha }} h(x)$ are considered as boundary particles.

It is important to note that the two definitions give the same results.

More advanced versions of the PFEM, allow using different values of ${\bar{\alpha }}$ for distinct parts of the domain. A possibility is the use of a lower value of ${\bar{\alpha }}$ for the elements whose nodes belonged to the boundaries in the previous time step and larger values for the internal elements.

An alternative technique to efficiently manage the boundary definition in the PFEM is the Constrained Delaunay tessellation. However, it is difficult to apply this method to problems with significant changes of the material boundaries, like in many free-surface fluids applications. Moreover, the separation and reconnection of different parts of the computational domain are difficult to be handled with this technique. On the other hand, the Constrained Delaunay tessellation can be used for some solid mechanics problems where no significant changes in time of the boundary are expected.

It must be noted that the local changes of topology due to the alpha-shape technique may produce the lack of preservation of the volume. If the sum of the volume of the erased elements is different than the one of the new simplices, the overall volume of the analysis is not conserved. Nevertheless, it must be noted that this inconvenience of the PFEM remeshing procedure is proportional to the mesh size and can be reduced to the desired accuracy by refining the mesh. An accurate description of the effect of the parameter $\alpha$ on the accuracy of the results can be found in [28].

2.2.3 Separation and Reconnection of Nodes and Subdomains

One of the main strengths of the PFEM lays in its capability to model separation and reconnection of parts of the computational domain and also single isolated particles. This is of great relevance for the simulation of several engineering applications and natural phenomena, such as breaking waves and splashes formation in free-surface fluid flows.

The identification of the parts detaching from the rest of the domain is done automatically by the alpha-shape method. When a boundary node belongs to a too distorted element, the criterion of Eq. (3) removes the element and the node is separated from the domain. After the separation, the particle motion is governed by the body force and the initial velocity which it is subjected to. At each new Delaunay triangulation, the detached particle becomes again a vertice of the new tessellation and its connection with the rest of the domain is evaluated. If the particle has approached enough the boundary, the element is not eliminated by the alpha-shape check, and it is again incorporated in the main mass. Figure 5 shows an example of separation and incorporation of a particle.

Fig. 5

a Separation of a fluid particle from the bulk. b Incorporation of a fluid particle in the bulk

The same considerations described for isolated nodes apply for groups of nodes linked together by elements. Nevertheless, in this case, the motion of these subdomains is computed through the FEM equations.

2.2.4 Adding and Removing Nodes

In mesh-based Lagrangian approaches, the nodes move as a consequence of the equation of motion. Hence, some nodes may concentrate in a region of the domain and, on the contrary, in another region, the number of nodes becomes too low to obtain an accurate solution. Such a situation may affect the quality of the mesh and, consequently, the accuracy and effectiveness of the FEM solver, especially in 3D problems. To overcome this drawback, it is possible to remove and add nodes in the mesh zones that need it. The same idea was originally proposed in the very early work on Lagrangian FEM for fluid flow [36]. It is important to remark that in the PFEM, insertion and removal of mesh nodes can be safely done because the mass is not associated with the nodes but with the elements, as in standard FEM.

In the PFEM literature, different algorithms to add and remove particle have been proposed. However, the key idea is always based on the two following concepts:

• If a node comes too close to another node (or to a boundary) the node should be removed (or moved to another location), see Fig. 6;

• If an element becomes too large, a new node should be inserted, see Fig. 7.

The nodes insertion can be done inside the element (for example, in its center of mass, Fig. 7b) or along an edge (for example, in the middle of its largest edge, Fig. 7c).

Fig. 6

Example of node removal. a Original mesh; b mesh without the removed node

Fig. 7

Example of node addition. a Original mesh; b addition of a new node inside an element; c addition of a new node along an edge

Removing nodes is also important to avoid the undesired situation of artificial leakage. This pathological situation mainly refers to fluid dynamics problems and it occurs when a node comes so close to a boundary wall that the alpha-shape removes the element connecting it to the boundary nodes. A void is then produced and material could pass through the wall nodes. To prevent this, it is recommended to remove those nodes that come to close to the fixed domain boundaries.

It is important noticing that the previous operations can be performed without generating a new mesh with the Delaunay triangulation and by only changing the local elemental topology. However, in this case, the mesh cannot be longer considered as a Delaunay mesh.

Adding and removing nodes can also be done without altering the total number of nodes. In this case, a new node is added only if a node can be removed from another position. Doing so, the mesh size tends to remain constant during the analysis. Moreover, this is also useful from the implementation point of view, as it enables the use of simplified data structure with dimensions fixed in time.

2.2.5 Mesh Refinement

The Delaunay tessellation coupled with the alpha-shape leads to uniform meshes if a unique threshold value ${\bar{\alpha }}$ in Eq. (3) is used for the whole domain. However, there exist situations in which a non-uniform mesh is recommendable, for example in the regions close to boundary layers, in nozzles or at the interface between different materials.

In these cases, one may define different values of the parameter ${\bar{\alpha }}$ depending on the position of the nodes. Consequently, nodes can concentrate automatically in the desired zones and get rarefied elsewhere. A similar mesh refinement technique was applied in [15].

Another possibility is to modify the alpha-shape criterion according to error estimation considerations. For example, in manufacturing processes to simulate complex large deformation, the insertion of particles is based on the equidistribution of the plastic power and the removal were driven by a Zienkiewicz-Zhu error estimator [104, 106].

2.3 Implications of Remeshing on the FEM Solution Scheme

The use of remeshing has some important implications on the FEM model that can be used in a PFEM framework. As explained in previous sections, a key characteristic of PFEM remeshing is that the new mesh is not created from scratch, but is built keeping the nodes of the previous mesh. This feature, allows avoiding data interpolation from the old to the new mesh when only nodal variables are used (typically in fluid dynamics), on the other hand, it limits the applicability of the method to linear elements. Higher-order interpolations would require to remap from mesh to mesh the unknowns placed in the middle of the edges of the elements. Furthermore, the element curvature will be lost at each step of the Delaunay triangulation. As a consequence, only linear shape functions should be used in a PFEM framework.

In this respect, we note that low order elements can be unstable for incompressible or quasi-incompressible problems, also in case of using mixed formulations (e.g. linear shape functions for velocity and pressure). This issue explains why the PFEM formulations for fluid mechanics are generally provided with some stabilization methods. A description of the most commonly used stabilization procedures in PFEM can be found in Sect. 3.6.

Another important implication of remeshing on the FEM solution is related to the continuous elimination of the elements. In standard methods using Gaussian integration, this makes necessary the implementation of specific techniques for the recovery of historical variables [96]. Nevertheless, remapping can be avoided in case of using a nodal integration scheme, as the historical variables are stored at the nodes [126]. More details can be found in Sect. 5.1.

PFEM remeshing also affects the choice of the reference configuration used in the FEM solution. Typically, three possible reference configurations can be chosen:

• 1. The initial configuration;

• 2. The configuration at the beginning of each time step;

• 3. The configuration at the beginning of each non-linear iteration.

Due to the excessive distortion of the mesh, the first choice (Total Lagrangian method) is practically unfeasible and only the other two choices (Updated Lagrangian methods) are exploitable. Keeping the reference configuration constant in a time step, the shape functions, their derivatives and all the geometrical quantities can be computed only once in the time step, but the computation of the deformation gradient has to be performed at each non-linear iteration. Instead, if the configuration at the last non-linear iteration is used, the deformation gradient coincides with the identity matrix, but the geometrical quantities must be updated at every iteration. In both cases, the connectivity does not change within the same time step.

3 PFEM for Fluid Dynamics Problems

The PFEM was originally conceived for the solution of free-surface fluid flow problems [52]. The method was designed to deal with the large motion of the fluid domain and to track its highly deforming free boundaries, which could eventually fold, break and reconnect. The range of possible applications of such a technology is wide and interests several branches of engineering and applied sciences.

In fluid dynamics, the automatic capability of tracking the evolving free-surface is not the only benefit of the PFEM. Thanks to its Lagrangian description of the motion, the convective terms do not enter in the PFEM governing equations. This is a great advantage as the convective terms are responsible for non-linearity, non-symmetry and non-self-adjoin operators, thereby complicating significantly the solution of the governing equations in an Eulerian framework and typically requiring the introduction of stabilization terms to avoid numerical oscillations.

In contrast, in the Lagrangian framework, the non-linearity appears because the governing equations are written in the unknown current deforming configuration, which may differ by large displacements from the reference one.

The absence of convective terms in the governing equations represents an important advantage of Lagrangian methods, such as PFEM versus standard Eulerian formulations. However, in fluid dynamics, the price to pay is the need to continuously remesh the computational domain. This fact has important implications on several aspects of the numerical solver, such as time integration and spatial discretization, the imposition of boundary conditions, or mass conservation. The following sections aim to analyze in detail the most important aspects of using the PFEM for fluid dynamics problems.

3.1 Fluid Dynamics Problem Statement

Let consider a moving fluid domain $\varOmega _t$ in the time interval [0, T]. The motion of the fluid body is governed by the Navier–Stokes equations. Introducing the velocity $\mathbf{u}=\mathbf{u}(\mathbf{x},t)$ and the Cauchy stress tensor $\boldsymbol{\sigma }=\boldsymbol{\sigma }(\mathbf{x},t)$, momentum balance and mass conservation read:

$$\begin{aligned}&\rho \frac{D {\mathbf{u}}}{D t}= \text{ div } \boldsymbol{\sigma }+ {\mathbf{b}} \quad&\text{ in } \; \varOmega _t \times (0,T) \end{aligned}$$

(4)

$$\begin{aligned}&\frac{D \rho }{D t} + \rho \text{ div } \,{\mathbf{u}}=0 \quad&\text{ in } \; \varOmega _t \times (0,T) \end{aligned}$$

(5)

where $\rho (\mathbf{x})$ represents the fluid density, ${\mathbf{b}}(\mathbf{x},t)$ the external body forces per unit mass, and $D \cdot$/Dt denotes the material time derivative. It is important to remark again that, due to the Lagrangian nature of the method, the convective term does not appear in the governing equations and the total time derivatives reduce to a local time derivative.

Many real fluids tend to an incompressible behaviour. In this case, the dependence of the density on time disappears, and Eq. (5) simplifies as:

$$\begin{aligned}&\text{ div } \,{\mathbf{u}}=0 \quad&\text{ in } \; \varOmega _t \times (0,T) \end{aligned}$$

(6)

which represents the mass conservation for an incompressible fluid. However, it must be underlined that a small amount of compressibility still exists in all the cases idealized as incompressible ones. It is interesting to note that considering small fluid compressibility can be convenient from the numerical point of view also for problems that would be classified as fully incompressible. Such an assumption is referred to as weakly compressible fluid hypothesis. For suitable small Mach numbers (much smaller than one), the weakly compressible model well approximates the incompressible limit. However, the algebraic nature of the problems changes: the system of equations is hyperbolic-parabolic for the compressible case and elliptic-parabolic for the incompressible one. The actual benefit of weakly compressible solvers versus incompressible ones is the possibility to compute explicitly the pressure, thereby avoiding the need of using a Poisson solver. This aspect explains its popularity in PFEM formulations, see e.g. [32, 44, 51, 112]

It is important to remark that in a PFEM framework, also real compressible fluids can be considered [21, 57]. In this case, the system of governing equations (Eqs. (4)–(5)) must be complemented by adding an energy conservation equation and an equation of state.

Typically, in fluid dynamics, the Cauchy stress tensor $\boldsymbol{\sigma }$ is decomposed into isotropic and deviatoric parts:

$$\begin{aligned} \boldsymbol{\sigma }=-p {\mathbf{I}}+\boldsymbol{\tau }\end{aligned}$$

(7)

where $p=p(\mathbf{x},t)$ is the pressure field, ${\mathbf{I}}$ the second-order identity tensor, and $\boldsymbol{\tau }$ is the deviatoric stress tensor, which is generally related to the deviatoric strain rate $\varvec{\epsilon }$ through a rheological law as:

$$\begin{aligned} \boldsymbol{\tau }= 2 \mu (\varvec{\epsilon }) \varvec{\epsilon }\end{aligned}$$

(8)

where $\mu = \mu (\varvec{\epsilon })$ is the viscosity. The deviatoric strain rate is obtained from the velocity field as:

$$\begin{aligned} \varvec{\epsilon }({\mathbf{u}})=\frac{1}{2} \left( \text{ grad }{\mathbf{u}} + \text{ grad }{\mathbf{u}}^T \right) - \frac{1}{3} (\text{ div }{\mathbf{u}}) {\mathbf{I}} \end{aligned}$$

(9)

Once the pressure field has been introduced, it is important to note that Eq. (5) can also be expressed as:

$$\begin{aligned}&\frac{D p }{D t} + \kappa \text{ div } \,{\mathbf{u}}=0 \quad \text{ in } \; \varOmega _t \times (0,T) \end{aligned}$$

(10)

where $\kappa$ represents the bulk modulus of the fluid.

3.2 Space Discretization and Stabilization

In the PFEM, a standard Galerkin approach is used to discretize in space the Eqs. (4)–(5) (or alternatively, Eqs. (4)–(10)). As already underlined in Sect. 2.3, in a standard PFEM framework, only linear shape functions are used to approximate the unknown variables. Introducing an isoparametric finite element discretization, the velocity and the pressure can be expanded in terms of the nodal vectors $\mathbf{U}$ and $\mathbf{P}$, respectively. The semi-discretized equations of motions read [23]:

$$\begin{aligned}&\mathbf{M}_{\mathbf{u}} {\dot{\mathbf{U}}} + \mathbf{K}\mathbf{U}+ \mathbf{D}^T \mathbf{P}= \mathbf{F}\end{aligned}$$

(11)

$$\begin{aligned}&\mathbf{M}_{p} {\dot{\mathbf{P}}} + \kappa \mathbf{D}\mathbf{U}= {\mathbf{0}} \end{aligned}$$

(12)

where $\mathbf{M}_{\mathbf{u}}$ and $\mathbf{M}_{p}$ are the mass matrices for velocity and pressure unknowns, $\mathbf{K}$ is the fluid matrix emanating from the viscosity term, $\mathbf{D}$ is the discretized divergence operator, and $\mathbf{F}$ is the vector of body forces and boundary conditions. The symbol $\dot{(\;)}$ represents a time derivative.

Using equal order interpolation for both the velocity and pressure unknowns, the LBB inf-sup compatibility condition is not fulfilled [9]. Hence, the formulation must be stabilized. In the PFEM literature, different kinds of stabilization procedures have been applied. For instance, the Finite Increment Calculus (FIC) formulation has been frequently used to stabilize the PFEM equations (see e.g. [52, 53, 83]). Alternatively, the Pressure Stabilizing Petrov-Galerkin technique has been used (see e.g. [13, 18]). In [61], the Algebraic Sub-Grid Scale stabilization technique is introduced to stabilize mixed pressure velocity PFEM formulation. Examples of other stabilization techniques can be found in [111, 113]. In [2], the authors propose to make use of stable elements belonging to the bubble family.

3.3 Time Integration

Equations (11)–(12) must be integrated in time and thus an approximation for the time derivative of pressure and velocities should be provided. Note that for really incompressible materials, the time derivative of pressure (or density) does not appear. Both implicit and explicit approaches can be used in a PFEM framework, although, typically, an implicit time integration has been preferred. For example, in [18, 83], implicit first-order scheme were used while in [33] second-order schemes were proposed. To simplify the numerical treatment of the coupled equations, fractional step or similar partitioned schemes have been used in several PFEM applications (see e.g. [2, 52, 86, 109]) .

Only recently, explicit time integration schemes have been used in PFEM formulations (see e.g. [17, 72, 73]. Explicit solvers are very appealing for fast dynamics problems and also for non-linear problems that may suffer from numerical issues of convergence. Moreover, they have good parallelization skills as a system of fully decoupled equations can be obtained. On the other hand, the explicit time integration is conditionally stable and the choice of the time step size is governed by the CFL (Courant, Friedrichs, Lewy) stability condition. This feature gives to mesh quality a crucial role for the computation efficiency of explicit methods because the presence of excessively distorted elements may lead to vanishing stable time step size and compromise the computation. For this reason [72] proposes an efficient mesh smoothing approach based on an elastic analogy to improve the worst elements, with a computational cost compatible with the frequent remeshing procedure required by the PFEM.

3.4 Rheological Law

The first PFEM works focused on Newtonian flows only ( [52, 87]). However, soon also non-Newtonian fluids were analyzed [15, 19].

In Newtonian fluids, the deviatoric stress $\boldsymbol{\tau }$ is linearly related to the deviatoric strain rate $\varvec{\epsilon }$ as:

$$\begin{aligned} \boldsymbol{\tau }= 2 \mu \varvec{\epsilon }\end{aligned}$$

(13)

where the fluid viscosity $\mu$ is constant.

On the contrary, non-Newtonian fluids are characterized by a non-linear relationship between the deviatoric stress and deviatoric strain rate. The Bingham law is one of the most commonly used models. The law is well representative of some geophysical flows, such as mudflow, lahars or debris flows but also fresh concrete and other types of materials, like some fresh paints or alimentary sauces. A Bingham material undergoes shear deformations only if the shear stress overcomes a fixed limit, $\tau _0$, called yield stress. The Bingham law can be expressed with the following piecewise equations:

$$\begin{aligned}&\boldsymbol{\tau }(\mathbf{u}) = 2 \mu \varvec{\epsilon }(\mathbf{u}) + \tau _0 \frac{\varvec{\epsilon }(\mathbf{u})}{\Vert \varvec{\epsilon }(\mathbf{u}) \Vert }&\text{ if } \, \Vert \boldsymbol{\tau }\Vert > \tau _0 \end{aligned}$$

(14)

$$\begin{aligned}&\varvec{\epsilon }(\mathbf{u}) = 0&\text{ otherwise } \end{aligned}$$

(15)

being

$$\begin{aligned} \Vert \varvec{\epsilon }(\mathbf{u}) \Vert =\sqrt{\frac{1}{2} \varvec{\epsilon }: \varvec{\epsilon }} \qquad \Vert \boldsymbol{\tau }\Vert =\sqrt{\frac{1}{2} \boldsymbol{\tau }: \boldsymbol{\tau }} \end{aligned}$$

(16)

The incrementally discontinuous behaviour of Eqs. (14) and (15) and the unbounded value of viscosity, introduce numerical difficulties in the solution scheme which can be avoided using an exponential smoothing approximation [98]:

$$\begin{aligned} \boldsymbol{\tau }(\mathbf{u}) = 2 \, {\tilde{\mu }} \, \varvec{\epsilon }(\mathbf{u}) = \left[ 2 \mu + \frac{\tau _0}{\Vert \varvec{\epsilon }\Vert } \left( 1 - e^{-n \Vert \varvec{\epsilon }\Vert } \right) \right] \, \varvec{\epsilon }(\mathbf{u}) \end{aligned}$$

(17)

where the apparent viscosity ${\tilde{\mu }}$ has been introduced. When $n\rightarrow \infty$, the original Bingham model is recovered. Due to its good adaptability to fluid dynamics solvers, the Bingham model with exponential regularization has been largely used in PFEM works, e.g. [15, 64, 103, 117].

The Bingham model was conceived for materials with fixed yield stress. However, in some cases, the definition of the yield value can depend on the characteristic of the material itself or by its intrinsic multi-physics (e.g. presence of the water in partially drained soils). A more general model, usually called frictional model, introduces non-constant yield stress depending on an effective pressure $p'$ and a friction angle $\phi$:

$$\begin{aligned} \tau _0 = p' \, \text{ tan }(\phi ) \end{aligned}$$

(18)

This relationship, which can be interpreted as a Mohr–Coulomb failure criterion (see e.g. [16]), can be incorporated into Eq. (17) to represent the behaviour of frictional viscoplastic materials. This model has been extensively used in the PFEM framework to simulate the flow of soil or granular material in [22, 61, 117].

3.5 Boundary Conditions

To be well-posed, the problem (4)–(5) has to be supplemented with an appropriate set of initial and boundary conditions. As underlined in the previous sections, the Lagrangian nature of the PFEM is very useful for the definition of evolving free-surfaces, which are automatically detected by the position of the external nodes of the domain.

On the contrary, the Lagrangian approach leaves much less flexibility to the treatment of constrained contours, which are represented by nodes with prescribed velocity. In Lagrangian methods, the external boundaries are defined by the position of the material particles (mesh nodes in the PFEM). These methods have some limitations when the boundary condition assigns non-zero velocity components to the nodes belonging to a fixed boundary. For these cases, the PFEM solver must be complemented by appropriate techniques. This section aims to shed light on the treatment of the different boundary conditions in the PFEM.

3.5.1 Free Surfaces and Interfaces

The tracking of the free contours of a fluid is a complex task for many numerical approaches and, in some cases, it requires the additional implementation of ad-hoc techniques. This is the case of all Eulerian finite elements or finite volumes formulations that must be complemented by techniques like Volume of Fluid [42] or Level set [97] methods.

On the contrary, the PFEM can automatically define the free-surface position without any specific technique, because, each node of the domain is tracked during the motion, and the free surface can be easily identified directly by the position of the nodes of free contours.

The same considerations apply for the detection of the interface between different fluids or between a fluid and a structure. The Lagrangian nature of the PFEM allows for the automatic identification of these interfaces and to follow them on time.

Depending on the physics of the problem, different conditions can be imposed on the evolving interfaces. For free-surface contours, it must be ensured that the normal component of the stress tensor should vanish. This is much preferable to the imposition in a strong form of null pressure, which can lead to the violation of mass conservation [51, 83, 112]. In the multi-fluid analysis, surface tension has significant importance and should be included at the interface between the different fluids. Examples of PFEM formulations for multi-fluid flows accounting for surface tension effects can be found in [74, 110]. The accurate imposition of boundary conditions at the interface is particularly important for fluid–structure interaction problems. The PFEM treatment of these boundary conditions is explained in Sect. 4.3.

3.5.2 Slip Boundary Conditions

No-slip boundary conditions between fluid and confining walls are generally considered in fluid dynamics. However, there exist cases in which a relative slip between the fluid and the surface is observed. A classical example is the so-called Navier slip boundary condition, which defines a linear correlation between the slip velocity and the tangential stress at the slipping surface. The amount of slip is defined through a parameter called slip length, which can range from the no-slip condition to free-slip. Free-slip conditions are particularly interesting as they may help to avoid using very fine boundary layer meshes.

Slip boundary conditions have been applied in a PFEM framework using different approaches. Cerquaglia et al. [13] proposes to use a layer of contact elements between the fluid and the wall to account for the slip effects (both in the weak or in strong forms). In [26], the implications in terms of mass conservation of using free-slip conditions in the PFEM are discussed. Cremonesi et al. [20] suggests the use of a Lagrangian–Eulerian approach to describe the slip condition. In this method, the fluid domain is solved in a standard Lagrangian way, while the slip interface is assumed to be composed by Eulerian nodes and solved accordingly.

Note that a symmetry plane can be considered a special case of slip boundary conditions where no normal velocity is present and free-slip in allowed in the tangent direction. Consequently, the same techniques used for the slip can be extended also to the symmetry boundary conditions.

3.5.3 Inflow and Outflow

The same complexity associated with slip boundary conditions arises also for inflow and outflow boundaries. The imposition of this kind of boundary conditions is of great importance in the simulation of real engineering problems but can be critical in a Lagrangian FEM-based method. When a velocity profile (or a pressure profile) is imposed on a boundary, the nodes belonging to that boundary move following the fluid velocity and consequently the definition of the boundary is lost.

Ryzhakov et al. [110] presents an inlet technique in which the inflow region is treated in a standard Lagrangian form: the nodes belonging to the boundary move with the prescribed velocity creating empty space which is then replaced with a new set of nodes. In [119], the same technique was used to simulate hydraulic channel conditions. Cremonesi et al. [20] suggests describing the inflow as an Eulerian boundary and the rest of the domain as Lagrangian. This technique overcomes the mesh issues of the purely Lagrangian strategy, but it requires the modification of the solver for the inlet nodes.

For the outflow conditions, in many relevant situations where no particular condition is imposed at the outlet surface, it is sufficient to remove from the mesh the Lagrangian nodes which have crossed the outflow limit, i.e. nodes that are outside the computational domain. On the contrary, in the case of prescribing the pressure or velocity at the outlet, the same techniques used for the inflow should be applied.

3.6 Mass Conservation

One of the critical points of free-surface fluid dynamics analysis is the unfulfillment of mass conservation. For the PFEM, like for other FEM-based methods, mass conservation may be violated due to two distinct reasons:

• 1. Due to the numerical solution of the governing equations;

• 2. Due to the free-surface tracking technique.

The first source of mass variation is due to the Galerkin finite element solution of mass conservation (Eq. (6)) because the inaccuracy of the numerical solver affects directly the total mass conservation. Moreover, the typical stabilization techniques used to circumvent the unfulfillment of the LBB condition, have the effect of relaxing the incompressibility constrain leading to mass variations. Therefore, it is crucial to rely on a fluid dynamics solver with good mass preservation properties [51, 83, 112]. It is also important to highlight that this type of mass conservation violation can be experienced in all the standard Eulerian and ALE approaches, not only in the PFEM.

The second mass violation source depends exclusively on the technique used to track the evolving fluid contours. In the PFEM, the fluid free surface is detected by the position of the nodes. However, during the remeshing, the connectivity may change (see Sect. 2.2.2), and this may lead to local topological variations of the fluid domain. These changes, besides perturbating locally the equilibrium configuration obtained at the previous step, can lead to a global variation in the volume of the computational domain, and so a global violation of mass conservation. It is important to remark that this volume variation reduces progressively by refining the mesh and with a proper definition of the parameter $\alpha$ of the alpha Shape method [18, 28]. A detailed analysis of the effects of the remeshing procedure on the volume conservation for free-surface incompressible fluid problems can be found in [28].

3.7 An Illustrative Example: The Dam Break Test

The application of different PFEM formulations to the dam break test, a classical benchmark for free-surface dynamics analysis, is here considered. The problem consists of a column of water initially located at the left part of a tank and sustained by a removable wall (dam) on its right side. At the initial time, the vertical wall is suddenly lifted and the water flows under the effect of gravity on a rectangular channel until it collides with the right vertical wall of the tank. The geometry of the problem is depicted in Fig. 8. This test was originally proposed in [59], where both experimental observations and numerical simulation results were provided.

Fig. 8

Dam break test. Geometry of the problem. $L=0.146$ m, $h=0.175$ m

This problem has been also solved many time with PFEM, both in 2D (see e.g. [13, 52, 53, 138]) and in 3D (see e.g [35, 72]). In Fig. 9, snapshots of a PFEM simulation have been synchronized with the experimental test observations of [59], showing the good capability of the method to reproduce the complex phenomena occuring in this test, such as breaking waves, splashes and strong impacts with solid boundaries.

Fig. 9

Dam break test. Snapshots of the simulation at different time steps compared with the corresponding experimental results. Pictures taken from [72]

Fig. 10

Dam break test. Time evolution of the front wave position

Figure 10 shows the time evolution of the front wave. Different PFEM formulations are compared showing a good agreement among them and with experimental data. However, all the numerical simulations (PFEM and not) show a slightly faster front advancement of the front versus laboratory observations. The reason for this small discrepancy lays in the modelling of the raising wall motion. In all the numerical simulations, the wall is removed instantaneously while it takes inevitably a finite amount of time in the experiment (see the comments in [43, 120]).

It is worth to note that this test has also been solved many times with the PFEM including a rigid obstacle in the middle of the container (see e.g. [32, 63, 72, 83, 121]).

4 PFEM for Fluid–Structure Interaction Problems

Since the very first works, the PFEM has been applied to the simulation of free-surface fluids interacting with moving solids. In the early applications, the solid objects were treated as rigid bodies [53, 87]. However, soon the method was extended also to elastic [44, 45] and elastoplastic [34, 137] bodies.

The main reason that explains the success of PFEM in the framework of fluid–structure interaction (FSI) analysis, lays undoubtedly in its accurate tracking of moving interfaces. This feature is extremely useful for the solution of a wide range of engineering and industrial problems with fluid–solid interfaces undergoing large motions.

Furthermore, as it will be detailed later, the body-fitted mesh arisen from the PFEM interface detection algorithm, enables an easy transfer of boundary conditions between the different bodies, creating the good basis for an accurate solution of the FSI problem.

Finally, it is important to remark that the PFEM allows for a very natural coupling with all types of Lagrangian FEM because it does not impose any restrictions on their solution schemes, enabling the re-utilization of existing solvers.

In this section, after a general classification of FSI methods, an extended review of theory and applications of the PFEM for fluid–structure interaction problems is given.

4.1 General Classification of FSI Methods

Let us consider a continuum domain $\varOmega ^t$ evolving in the time interval $\left[ 0, T \right]$. The domain is constituted by two non-overlapping subdomains: a fluid one, $\varOmega ^t_F$ and a structural one, $\varOmega ^t_S$.² Let us define the subdomains boundaries: $\varGamma _F^t=\partial \varOmega ^t_F$ for the fluid and $\varGamma _S^t=\partial \varOmega ^t_S$ for the solid. The fluid–structure interface is given by $\varGamma _{FSI}^t= \partial \varOmega ^t_F \cap \partial \varOmega ^t_S$ (Fig. 11).

Fig. 11

Domain of the FSI problem

At the fluid–solid interface $\varGamma _{FSI}^t$, dynamic and kinematic coupling conditions should be enforced:

$$\begin{aligned}&\mathbf{u}_{F} = \mathbf{u}_{S}&\text{ on } \, \varGamma _{FSI}^t \end{aligned}$$

(19)

$$\begin{aligned}&\boldsymbol{\sigma }_{F} \cdot \mathbf{n}= \boldsymbol{\sigma }_{S} \cdot \mathbf{n}&\text{ on } \, \varGamma _{FSI}^t \end{aligned}$$

(20)

where $\mathbf{u}_F$ and $\mathbf{u}_S$ are the fluid and solid velocities, $\boldsymbol{\sigma }_{F}$ and $\boldsymbol{\sigma }_{S}$ the Cauchy stress tensors in the fluid and solid domains and $\mathbf{n}$ is the normal to the interface $\varGamma _{FSI}^t$.

A large variety of schemes has been proposed in the literature to solve FSI problems. A first classification of these depends on the degree of coupling of the method. An algorithm that enforces both the kinematic and dynamic transmission conditions across the fluid–structure interface is defined as strongly coupled. If instead, an algorithm does not satisfy exactly the coupling conditions, is defined as weakly coupled or loosely coupled

Another classification of FSI algorithms takes into account the structure of the solver. A method is said to be monolithic if a unique solver is used to solve the fluid–structure interaction problem. In this approach, the coupling conditions are enforced exactly, leading to a strongly coupled scheme and preserving accuracy and stability. This represents the main advantage of monolithic strategies together with the fact that they do not suffer from added mass effect problems. However, monolithic methods do not allow reusing existing fluid and solid solvers, and they potentially lead to large ill-conditioned linear systems, as both fluid and solid contributions are contained therein.

Alternatively, partitioned (or staggered) schemes solve fluid and solid sub-problems independently, and then they couple the solutions via transmission conditions. Partitioned methods can provide strong coupling, for example through sub-iterations or predictor/corrector techniques, or weak coupling otherwise. Staggered schemes allow the reuse of existing codes and to solve smaller and better conditioned linear systems. On the other hand, in some cases, partitioned schemes can suffer convergence issues.

4.2 Detecting Fluid–Structure Interface in PFEM Framework

One of the key features of the application of the PFEM to FSI problems is the possibility to exploit a fully Lagrangian description of the fluid and solid subdomains. Therefore, the mesh evolves in time with the motion of the different bodies, enabling accurate tracking of the evolving fluid–solid interface.

In the most general version of the PFEM, the interface detection is done by superposing a set of fictitious fluid particles to the structural boundary surfaces: these particles have the fluid physical properties, but, in the beginning, they are just used for the FSI boundaries identification (Fig. 12a or d). These particles are involved in the remeshing step of the fluid analysis, as the real fluid nodes. The resulting Delaunay mesh discretizes the whole fluid domain and creates the interface elements between the fluid and the solid domain (Fig. 12b or e). After that, the alpha-shape technique makes a selection of the interface elements eliminating those too large or too distorted. At this point, if all the interface elements are removed, the fluid and solid domains are not interacting any longer (Fig. 12c). On the contrary, if some interface elements connecting solid and fluid parts are not removed (Fig. 12e), a coupled analysis is performed (Fig. 12f).

Fig. 12

FSI for conforming fluid and solid meshes: a Fluid as a set of particles, structure as a mesh of quadrilateral elements with fictitious fluid particles at its boundary. b Delaunay triangulation. c Application of alpha-shape method. d same as a but subdomains in contact. e Delaunay triangulation. f Application of the alpha-shape method

It is important to note that, depending on the location of the fictitious nodes, different body-fitted strategies may arise. If the fictitious nodes coincide with the solid nodes, a conforming-mesh interface is obtained. If, on the contrary, the fictitious nodes do not coincide with the solid nodes, a non-conforming mesh strategy is obtained. This latter method enables the use of different mesh sizes for the fluid and the solid elements at the interface, but, it requires mapping the variables at the interface to transfer accurately the boundary conditions.

4.3 Literature Review of PFEM for Fluid–Structure Problems

FSI problems have been solved with the PFEM using different strategies. In all cases, the PFEM is used to solve the fluid domain and to detect the interface, while the FEM is used for the solid solution. Practically speaking, the remeshing step affects only the fluid parts, while the solid discretization remains unchanged during the entire analysis.

The full monolithic scheme has been used on different occasions. The first contribution in this category is [44], where the authors solve the FSI problem in a general and unified Lagrangian framework, exploiting the similarities between the Newtonian model (used for the fluid) and the hypoelastic one (used for the solid). The method was extended to quasi-incompressible solids in [33] and to thermal-coupled problems in [34]. Recently, in [27] the method has been adapted to a nodal-integration framework. [111, 113] presented a monolithic approach with a global pressure condensation which enables the definition of a purely displacement-based linear system of equations. A matrix-free technique is used for the solution of such a linear system. In [138,139,140], ill-conditioning issues of the monolithic formulation had been overcome by applying a fractional step approach to the linear system, segregating pressures and velocities unknowns into smaller systems of equations.

There also exist several partitioned PFEM approaches for FSI problems. In [53], a Gauss-Seidel technique was used to treat the FSI problem, although only rigid bodies were considered. In [18], a partitioned Dirichlet-Neumann algorithm solving the fluid with the PFEM and the elastic solid with a displacement-based FEM was presented. A similar algorithm was applied in [114] to the simulation of sea-landing of aerial vehicles. Meduri et al. [73] proposed a partitioned approach allowing for a non-conforming mesh at the interfaces and different time steps for the fluid and solid solutions. This work also showed the successful coupling of a PFEM fluid solver with commercial software used for the solid parts. Recently [14] presented a fully partitioned Lagrangian framework using an Interface Quasi-Newton Inverse Least Squares strategy to avoid added mass effects.

4.4 An Illustrative Example: The Dam Break with an Elastic Obstacle

The collapse of a water column against a deformable membrane is a typical benchmark problem for FSI analysis with free-surface fluid flows [123]. The geometry of the problem is depicted in Fig. 13. Following [123], the following geometrical parameters are used:

$$\begin{aligned} \text{ L }=14.6 \; \text{ cm } \qquad \text{ h } = 1.2 \; \text{ cm } \qquad \text{ d } = 0.8 \; \text{ cm } \end{aligned}$$

(21)

The fluid is water and its physical properties are: viscosity $\mu =0.001\, \text{ kg/ms }$, density $\rho _f=1000\, \text{ kg/m}^3$. The solid is elastic having Young modulus $E=1000\, \text{ kPa }$, Poisson ratio $\nu =0$ and density $\rho _s=2500 \, \text{ kg/m}^3$.

Fig. 13

Collapse of a water column on a elastic object. Geometry of the problem

Figure 14 shows some representative snapshots of the simulation. The water column collapses after the removal of the vertical wall and spreads on the rectangular channel until hits the vertical elastic membrane (Fig. 14a). The obstacle blends under the effect of the fluid impact (Fig. 14b). Water slips on the deflected membrane and collides with the terminal vertical rigid wall (Fig. 14c). After that, the water volume fills the right part of the container (Fig. 14d) impacting again on the right part of the solid object (Fig. 14e, f).

Fig. 14

Collapse of a water column on a elastic object. Picture from [71]

Figure 15 plots the time evolution of the horizontal displacement of the left upper corner of the solid object obtained by different PFEM strategies. The curves show an overall good mutual agreement.

Fig. 15

Collapse of a water column on a elastic object. Comparisons of the PFEM results of Idelsohn et al. [44], Zhu and Scott [113], Meduri et al. [73], and Cerquaglia et al. [14]

5 PFEM for Non-linear Solid and Contact Mechanics

Although the PFEM was originally designed for fluid dynamics and FSI applications, the method was soon extended to non-linear solid mechanics. Indeed, a Lagrangian approach for large-deformation problems was appealing for several industrial and natural processes, where solid bodies undergo so large motions to behave like a fluid. This is the case, for example, of several manufacturing processes and many geotechnical applications. In this section, we will refer to PFEM formulations for non-linear solid mechanics, in short PFEM-solid, as those PFEM strategies that make use of historical variables defined inside the elements.

The first PFEM-solid formulation [95] was applied to complex industrial processes, such as metal forging, machining, or powder filling, and showed the capability of the method to track accurately the deforming shape of the material and to deal with the complex interactions between the different solid bodies. This first work opened up the way to many other PFEM-solid formulations. References [84, 104, 106, 107] applied their PFEM-solid methods to different types of manufacturing processes, [4, 11, 12] analyzed tunneling and excavation applications, while bed erosion in river dynamics was tackled in [81, 82, 86]. Several PFEM-solid formulations have been proposed in the field of soil mechanics and geotechnal engineering, especially for the modelling of frictional materials and granular flows [10, 24, 56, 66, 127, 128] and for different types of geomechanics problems [75,76,77,78]

5.1 Historical Variable Recovery

One of the most critical points of PFEM-solid formulations is the management and conservation of the historical variables during the remeshing step. The solution accuracy depends on how well the historical information is preserved along with the transition from the previous mesh to the new one.

The PFEM remeshing procedure consists of erasing all the elements of the distorted mesh and creating the new tessellation over the cloud of points composed by the mesh nodes (Sect. 2). For PFEM-solid formulations using Gaussian integration, this remeshing strategy implies that all elemental information must be transferred from the elements of the old mesh to the nodes, and then, from these to the new mesh elements.

This remapping procedure leads inevitably to smooth the historical solution, also if the element connectivity is not changed. To limit this drawback, Oliver et al. [95] proposed to remap not the whole historical variable but just its time step increment. Hence, at the end of a generic computation step $\left[ t^n;t^{n+1}\right]$, the nodal stresses $\varvec{{\bar{\sigma }}}^{n+1}$ are obtained from the incremental elemental stresses $\varDelta \varvec{\sigma }^{n+1}$ as follows

$$\begin{aligned} \varvec{{\bar{\sigma }}}^{n+1}=\varvec{{\bar{\sigma }}}^{n}+ \left[ \varvec{M_{\sigma }} \right] ^{-1} \int _{\varOmega _n}\varvec{N_{\sigma }} \varDelta \varvec{\sigma }^{n+1} d\varOmega \end{aligned}$$

(22)

where $\varvec{M_{\sigma }}$ is a standard mass-type matrix and $\varvec{N_{\sigma }}$ is a transfer matrix using the shape functions [95].

At the beginning of the new time step, the nodal stresses are mapped on the new Gauss points with standard FEM interpolation procedures as:

$$\begin{aligned} \varvec{\sigma }^{n+1}\left( {\textit{\textbf{x}}}_{gp}\right) =\sum ^{nn}_{i=1} N_i\left( {\textit{\textbf{x}}}_{gp}\right) \varvec{\sigma }^{n+1}_i \end{aligned}$$

(23)

where ${x}_{gp}$ is the position of the Gauss point, N is the shape function and nn is the number of nodes of the element.

In order to avoid excessive smoothing of historical variables [104, 106] proposed to transfer the information directly from the Gauss points of the previous mesh to the closest ones of the new mesh. This way, all information on elements that did not modify their connectivity during the remeshing, is perfectly preserved. However, in the parts of the mesh affected by the change of connectivity, this technique introduces inevitably an error, whose magnitude depends on the distance between the original and the new positions of the Gauss points.

Recently [56, 124, 126] proposed to use nodal integration rather than the Gaussian one to improve the preservation of historical variables. In nodal integration schemes, all variables (including stresses and strains) are stored at the nodes. Hence, this information is not erased during the remeshing and neither is affected by remapping operations. This solution appears to be a good way to preserve the historical variables information in a PFEM-solid formulation. Nevertheless, we note that PFEM-solid strategies with nodal integration, analogously to elemental integration methods, are not completely insensitive to mesh variations, such as change of connectivity or the elimination/creation of elements. These mesh modifications, perturbating the equilibrium configuration obtained at the previous computation step, may affect the convergence and accuracy of the solution, also when nodal integration is used.

5.2 Contact Problems

Since its first approach to non-linear solid mechanics, the potential of the PFEM for dealing with contact interaction between solids was explored. The first study that showed the suitability of the method for contact problems was [81]. In this work, the authors studied FSI problems where the solid bodies, dragged by the fluid motion, could eventually hit the walls of the computational domain.

The PFEM contact algorithm uses a mesh of interlayer elements between the boundaries of the interacting solid bodies. This auxiliary mesh is created following the same steps of the standard PFEM remeshing procedure (Sect. 2). Nevertheless, not all the elements fulfilling the alpha-shape criterion are used to compute the contact forces, but only those having a size smaller than a fixed critical value ($h_c$). Over these active contact elements, the elastic and frictional forces are computed either with penalty methods or using Lagrangian multipliers.

This methodology was also used in [95] to solve mutual contacts in manufacturing processes. The paper emphasized the property of the auxiliary PFEM mesh, called anticipating interface mesh, to recognize in advance the solid parts getting in contact and to impose the contact constraints in a diffusive manner.

Inevitably, the accuracy of this method depends on the size of the mesh and, in particular, on the critical distance $h_c$, whose value affects the timing of the contact and the size of the gap between the solid interfaces. To improve this aspect [95] proposed to apply an artificial contraction to the solid boundaries to capture more accurately the contact time and to reduce the distance between the interacting solid bodies. The same methodology was also used in the so-called Contact Domain Method [96] and formalized in [115].

The first three-dimensional (3D) application of the PFEM contact algorithm was presented in [86], where the capabilities of the method were proved against complex multi-body interactions, either in the presence of water or not.

A different contact algorithm was used in the PFEM-solid formulations [127, 129] for the simulation of granular flows and landslides. In these works, the auxiliary mesh was used only to detect the colliding solid boundaries, but not to calculate the contact forces, which were computed with a strategy originally conceived for discrete element methods (DEM) [60].

6 The PFEM for Other Coupled Problems

6.1 Multi-fluid Problems

Heterogeneous fluid flows are involved in several natural phenomena and engineering applications. The numerical simulation of mixing processes of immiscible fluids is particularly complex in the presence of very different physical properties (density and viscosity) and multiple and articulated interfaces. In this context, the success of the numerical simulation mainly depends on the ability of the method to track accurately these interfaces and to model the phenomena taking place there.

In the Eulerian framework, standard interface-capturing methods, such as Level Set [97] and Volume of Fluid [42] should be introduced. However, these front-capturing methods have some difficulties to avoid the smoothing of the interface, in particular in unsteady flows. On the contrary, the PFEM can automatically track the evolution of many sharp interfaces, thanks to its Lagrangian nature and to the definition of the material properties at the mesh nodes.

The PFEM multi-fluid interfaces can naturally fold, break and merge, analogously to the fluid free surface (Sect. 3). This represents a key capability of the method in the framework of heterogeneous fluids simulation. Furthermore, in the PFEM, the finite elements located at the interface can be properly enhanced to deal with the numerical issues arisen from the abrupt jump of material properties. As in a standard FEM, to deal with localized pressures jumps and instabilities, ad-hoc stabilization and pressure enrichment at the interface should be introduced [46, 48]. In [46], fluids with different density are considered, in [48], also the viscosity jump is introduced. Moreover [48] shows that discontinuous pressure fields can avoid errors in the incompressibility condition. Mier-Torrecilla et al. [74] presents the PFEM for large jumps in the physical properties, including also surface tension. Figure 16, taken from [46], shows the mixing process of two fluids with different density and initial temperature.

Fig. 16

Examples of multi-fluid problem. Mixing of two fluids with different densities and initial temperature conditions. Pictures from [46] plot the evolution of nodal density

In the PFEM, the interface separating two materials can be considered inside an element (elemental interface) or along its edges (nodal interface) [46]. Elemental interfaces are more stable as they do not change much when remeshing is performed. On the other hand, nodal interfaces are more accurate because they allow representing exactly the gradient pressure jump that normally occurs when with a jump in the density.

One of the most crucial aspects of the PFEM applied to multi-fluid flows, is the conservation of fluid volumes during the remeshing step. In fact, on the one hand, the change of connectivity may modify locally the path of the interface elements, and, on the other hand, the insertion and removal of mesh nodes may benefit or disadvantage one fluid versus the others. For this reason, it is very important to control the remeshing operations at the interface to avoid or limit the artificial volume variations of the involved fluids [122].

6.2 Thermal Coupled Problems

Thermally coupled flows are of great relevance for many fields of engineering and technology as well as for many natural phenomena. Since its origin, the PFEM has paid great attention to these kinds of problems. The first work in the field [1] showed that the Lagrangian nature of the PFEM can be very useful to model accurately thermal convection (Fig. 17). The work was then extended to 3D analysis in [2]. The potential of the method for dealing with thermal problems was further proved in a multi-fluid framework [46, 47].

Fig. 17

Examples of thermal coupled problem. Pictures from [1]

Taking advantage of the capability of the method to deal with large changes of topologies, the PFEM has been also applied to the simulation of melting problems accounting for phase changes. The first approach to phase change analysis was presented in [90] for the simulation of the melting and spreading of polymers. The solid objects were modelled as highly viscous fluids with temperature-dependent viscosity. The same approach was extended to 3D simulations in [88, 91] and used in [58, 69] to reproduce a small-scale fire test used to assess the flammability of polymers.

An immersed approach was proposed in [70], where a burning polymer was modelled with the PFEM and the surrounding air with an Eulerian formulation. This hybrid method allowed for the solution of the energy equation for both subdomains on the Eulerian mesh.

The first extension to thermally coupled fluid–structure analysis was presented in [85]. In this work, the fluid was modelled with the PFEM, while an elastoplastic FEM model was used for the structure. Melting phenomena were modelled by transferring the external solid elements to the fluid domain when they fulfilled a melting criterion. The same approach was extended to 3D analysis in [34] and applied to the simulation of hypothetical scenarios of nuclear core melting during a severe accident in a nuclear plant.

Obviously, thermal effects are also of paramount importance for a wide range of manufacturing processes whose overview can be found in Sect. 7.4.

7 Advanced Applications of the PFEM

This section aims to analyze four of the main application fields of the PFEM to engineering and environmental problems.

7.1 Hydraulic Engineering

Fluid dynamics has been the first field of application of the PFEM. The accurate simulation of the fluid free surface, also in the presence of breaking waves and splashes, together with the automatic modelling of fluid convection and the good energy conservation properties, make the PFEM an ideal tool for the analysis of different hydraulic engineering problems.

The first PFEM simulation of hydraulic laboratory tests was presented in [63]. In this work, a thorough comparison of the PFEM results against experimental observations is carried on four different free-surface flow configurations.

The suitability of PFEM in modelling wave propagation was clearly demonstrated in [92,93,94] by reproducing accurately different types of propagating waves in laboratory channels. Zhu et al. [136] extended the application field by also considering the presence of solid structures, reproducing the situation of tsunami waves loading on a bridge. Recently [79] showed the successful application of the PFEM to the simulation of water dam break impacting water at rest and creating impulse waves.

Dam engineering is another main field of application of the PFEM. The first works in this area [62, 64, 65] used a hybrid FEM-Eulerian approach to simulate overtopping and failure of rockfill dam and the related seepage phenomena. Salazar et al. [118] studied a real dam geometry and modelled the 3D air-water interaction to estimate the air demand at the bottom outlets. On the other hand [119] focused on the water shock-waves that form at the exit of dam spillways. Figure 18, taken from [119], shows a view of the real dam spillway and the 3D simulation with PFEM.

Fig. 18

View from downstream of dam spillway. Comparison between numerical and experimental results. Picture from [119]

7.2 Granular Flows

Granular flows are involved in many fields ranging from geophysics and geotechnical engineering to mining, pharmaceutical, and alimentary industries. Their numerical simulation is challenging because, depending on the flow regime, granular materials can behave as solids, fluids or gases, and, in general, multiple phases may appear simultaneously. Granular flows have been approached with PFEM formulations both in a solid and a fluid mechanics framework.

Zhang et al. [127] modelled the granular material as a rigid-plastic body and used PFEM-solid formulation to reproduce quasi-static and dynamic granular flows, while in [128, 131], an axisymmetric PFEM with rate-independent plasticity was employed.

A PFEM-solid formulation was also used in [24] considering large strains plasticity and a Drucker Prager model provided yield surface. The method was applied in [10] to complex industrial applications, such as silo discharge and tumbling mills.

Recently, a PFEM-solid formulation has been validated against several experimental tests of the collapse of granular columns [66], while in [56] a granular flow simulation has been carried out using a nodal integration method.

There exist also examples of PFEM-fluid formulations for granular flow simulation. A frictional viscoplastic model was used in [116] to simulate a sandy flow on a slope and impacting water at rest. Alternatively [29] used a regularized $\mu$(I)-rheology to model the 2D and 3D flow of dense granular material. Figure 19 shows a series of consecutive snapshots of the collapse of a cylindrical granular column taken from [29].

Fig. 19

Example of collapse of an cylindrical granular column solved with the $\mu (I)$-rheology. Pictures from [29]

7.3 Landslides

Landslides are one of the most destructive and dangerous natural hazards. Each year they cause billions of euros in damages and thousands of casualties worldwide. Forecasting the effects of these catastrophic events on the civil and natural environment is a complex task. In particular, it is hard to predict the landslide runout dynamics, due to the complex characterization of the landslide bulk material, the highly deforming shape of the sliding bodies, and the large size of the events. In this scenario, numerical methods can help to reduce the uncertainties in landslide events prediction and to better evaluate the associated risk. In the last decade, the PFEM has shown to have some potential in this field thanks to its capacity of capturing evolving free surfaces and the possibility of using accurate constitutive models for the landslide material.

The first attempt to simulate landslides events with the PFEM was reported in [19]. The work focused on the impulse water waves induced by a landslide, a dangerous multi-hazard scenario affecting above all mountains’ natural and artificial reservoirs. In this first contribution, the landslide material was modelled as a non-Newtonian Bingham fluid. Similar landslides impulse wave events were analyzed in [116] using a frictional viscoplastic model. Preliminary 3D results were also presented in the same work. 3D simulations of sliding material on an unstable slope have been also presented in [16] focusing on the importance in the landslide runout of considering the slip velocity between the flowing mass and the basal surface.

In [129], a mathematical programming framework is introduced to simulate landslide with a plasticity model using a Mohr–Coulomb yield criterion. In [132], an elastic-viscoplastic model for progressive failure analysis of sensitive clays is presented while in [133], its application to landslide is shown. The same ideas have been extended to simulate submarine landslides in [130] and used to analyze the Saint Jude landslide case study in [134]. Very recently, a large scale PFEM model has been used to reproduce the Vajont landslide and the consequent impulse wave in the hydroelectric reservoir [30]. Figure 20 shows some snapshots of this 3D analysis.

Fig. 20

Example of landslide-water interaction modelled with the PFEM (the Vajont landslide). Pictures from [30]

7.4 Manufacturing

The simulation of manufacturing processes represents one of the relevant areas of application of the PFEM. The capability of the method to deal with large deformations, complex contact interaction, and constitutive models, explains the large number of PFEM works on manufacturing processes. Furthermore, typically, these problems also include coupled thermal effects that can be easily handled with the PFEM.

Traditionally, in the PFEM, manufacturing processes have been tackled via a solid mechanics framework. The first PFEM-solid formulation applied to this highly non-linear analysis is reported in [95]. In this work, a PFEM-solid formulation was used to reproduce industrial metal forging, machining, or powder filling problems. PFEM can be also efficiently used to simulate cutting processes in which phenomena of friction, adiabatic shear bands, excessive heating, large strains, and high strain rates are involved. Oñate et al. [84] showed examples of extrusion of steel plates, forging of metal pieces and cutting of metals. Figure 21 shows some representative results of a cutting process modelled with the PFEM. Other examples of manufacturing processes can be found in [104, 105, 107] which are focused on the simulation of the segmental chips generated during metal cutting processes.

Fig. 21

Examples of a cutting problem. Pictures from [84]

Some manufacturing processes involve so large deformation of the contours that are more conveniently approached in a fluid dynamics framework than in a solid one. For example, mould filling and casting problems are modelled with a PFEM fluid dynamics formulation in [89], while the application to forming problems is presented in [84]. These latter processes were also considered in [109], where the heat equation is coupled to the mechanical model through a temperature-dependent viscosity to simulate glass forming processes. An axisymmetric PFEM formulation was proposed in [108] to simulate the forming of glass bottles.

8 Recent Advancement on the PFEM

8.1 The PFEM of Second Generation (PFEM-2)

An alternative PFEM technique to solve the incompressible Navier–Stokes problem with large time steps has been presented in [50]. The key idea of the method lays in the X-IVAS (eXplicit Integration following the Velocity and Acceleration Streamlines) method, which consists of integrating the convective terms following the streamlines rather than the particle trajectories. The X-IVAS technique has been coupled with the standard PFEM giving rise to the so-called Particle Finite Element Method—Second Generation (PFEM-2).

Two different versions of PFEM-2 have been proposed. The first one, the PFEM-2 with moving mesh, is based on the standard PFEM scheme and it creates a new mesh using the position of the nodes. Conversely, in the PFEM-2 with fixed mesh, the initial background mesh is kept unchanged during the analysis and the information is mapped on this fixed mesh. The first technique needs to re-build the mesh when is too distorted, as it happens in standard PFEM, and it is very efficient for free-surface flows. The approach with fixed mesh, on the contrary, maps the variables on the mesh avoiding the generation of a new one, and it is particularly suited for fluid flow problems in closed domains.

An interesting feature of PFEM-2 is the possibility to use an explicit time integration independently on the Courant number. The method remains explicit and stable independently on the mesh size. The time step is established following only accuracy considerations, besides the limits given by the Fourier number.

The PFEM-2 has been applied successfully to different engineering problems. In [55], the ideas presented in [49, 50] were generalized for multifluid flows with large time steps. In [38], an extended validation of the method for academic problems is presented. Gimenez et al. [37] shows the potential of PFEM-2 to simulate industrial problems of large time duration. An application of the method to jet atomization simulation can be found in [39]. Becker and Idelsohn [5] shows the potential of PFEM-2 to simulate large scale landslides events. FSI problems are tackled with a monolithic PFEM-2 approach in [6]. Finally [8] shows an application of the method to the simulation of sediment transport phenomena in rivers.

8.2 PFEM with Nodal Integration

Traditionally, the PFEM has been formulated for standard elemental integration, storing stresses and strains at the Gauss points. However, in PFEM with Gaussian integration, due to the continuous elimination of the elements done during the remeshing steps, it may be required to transfer the elemental information from the old mesh to the new one. This is avoided in fluid dynamics problems, where the measures of stresses and strains are computed from the scratch at each time step, but is mandatory for non-linear solid mechanics methods that require the storage of historical variables. Remapping procedures, besides having a certain computational cost, introduce interpolation errors into the numerical scheme and cause the smoothing of the historical variables (Sect. 2).

On the other hand, in nodal integration methods, stresses and material historical variables are computed and stored at the mesh nodes [99]. Consequently, a PFEM strategy with nodal integration does not require variable remapping procedures along the remeshing step.

This feature motivated recent research on the use of nodal integration in a PFEM framework [56, 124, 126]. The method, called by the authors Smoothed Particle Finite Element Method, took inspiration from the Smoothed Finite Element Method [135] and was successfully applied to 2D geomechanics problems with large deformations.

The use of nodal integration is also appealing in fluid dynamics analysis. FEM models with nodal integration are expected to suffer less from the low quality of the mesh [125]. This feature is particularly important for PFEM models because it would allow reducing the number of remeshing events and the associated drawbacks. Good resilience to mesh distortion was also found in the first application of PFEM with nodal integration to free-surface fluid dynamics problems [31]. The same work showed that the scheme gives a faster stress convergence than the standard element-based method. On the other hand, using stiffness matrices with larger bandwidth for the same mesh, nodal methods have a higher computational cost for the solution of the linear systems.

Very recently [27] showed the accuracy of PFEM with nodal integration for the solution of FSI problems in presence of free-surface fluid flows.

From a broader perspective, the demonstrated suitability of PFEM with nodal integration for fluid dynamics [31], solid mechanics [56, 124, 126], and FSI analysis [27] opens the field to a unified treatment of a general multi-physics continuum within a unique PFEM framework.

9 Conclusions

The PFEM is a powerful and well-assessed numerical technique that has been extensively used to simulate complex engineering problems. This work aimed to give an extended overview of the method describing its basic ideas, the main strengths and weaknesses, and the spectrum of applications.

In the first part of the work, the method has been described in its general form without focusing on specific physics or application fields. The remeshing procedure of the PFEM has been described in detail. The properties and implications of the Delaunay Tesselation, used to re-build the elemental connectivities, are analyzed and the alpha-shape method used to identify internal and external boundaries, is presented in details. Particular attention has been devoted to describing the techniques to improve the mesh quality, as the insertion and removal of mesh nodes, and to highlight the capability of the method to reproduce separation and reconnection of parts of the computational domain.

After this general description, different physical problems have been analyzed separately focusing on the particularities of their respective PFEM solution scheme. Free-surface fluid dynamics was historically the first field of application of the PFEM and, for this reason, has been analyzed first. Crucial aspects of fluid dynamics analysis and their treatment with the PFEM have been presented. Particular care has been devoted to the modelling of the different boundary conditions and to discuss mass conservation issues. A benchmark problem for free-surface fluid dynamics solved with different PFEM formulations has also been presented, showing the capability of the method to deal with complex phenomena such as breaking waves, strong impacts and water splashes formation. Then, the application of the PFEM to fluid–structure interaction (FSI) problems has been shown. After a general classification of the FSI methods, the PFEM strategy to track the evolving fluid–solid interfaces has been accurately described. An extended literature review of the different PFEM formulations for FSI has also been provided. Finally, different PFEM solutions of the collapse of a water column against a deformable membrane have been presented, confirming the suitability of PFEM for solving complex FSI problems with large solid–fluid interface motions. Finally, the application of PFEM to non-linear solid mechanics (PFEM-solid) has been described. The advantages and disadvantages of PFEM-solid formulations have been highlighted. Particular attention has been devoted to describing the re-mapping methodologies used to recover the historical information during the remeshing step, and the strategy used to model solid–solid contact interaction. An extended literature review of PFEM-solid formulations has been also provided.

Other main multi-physics problems tackled with the PFEM have been analyzed. First, the suitability of the method to multi-fluid flows has been highlighted and then thermally coupled problems have been analyzed. The PFEM has shown to have a high potential for both multi-physics problems, proving to be able to handle multiple, articulate and sharp fluid–fluid interfaces, and to model naturally and accurately thermal convection in different application fields.

Interesting engineering and industrial applications of the method have been also presented. Obliviously, many other applications of the method have been reported in the literature. In this review, we limited our analysis only to the applications of PFEM to hydraulic engineering, granular flows, manufacturing and landslides problems. For all the mentioned fields, the potential of the PFEM has been demonstrated and the main works have been referenced.

Finally, the recent advances in the method have been presented. First, the second generation of the PFEM (PFEM-2) has been introduced. Then the use of nodal integration in the PFEM has been described. PFEM-2 allows reducing computational cost enlarging the time steps of the analysis. Nodal integration algorithms give the possibility to avoid remapping and reducing remeshing events. The potential of both methods has been generally presented together with the significant literature.

In conclusion, this work aimed to show the potential of the PFEM for solving a broad range of applications in engineering and applied sciences. By highlighting the advantages and disadvantages of the method, this work wants also to stimulate future improvements of the PFEM technology and to widen its field of application.

Appendix: The Lagrangian Framework

Consider a continuum body occupying an evolving domain $\varOmega _t$ in the time interval [0: T]. The domain occupied by the body a time $t=0$ is defined as $\varOmega _0$ and called initial configuration (see Fig. 22). The domain of the body at a generic time instant is denoted as $\varOmega _t$ and defined as current configuration. Moreover, we define reference configuration the configuration at which the equation are referred [7].

Fig. 22

Description of the motion

In the reference configuration, the position vector of a material point is defined as material coordinates (or equivalently Lagrangian coordinates) and labelled as $\mathbf{X}$. Similarly, in the current configuration, the position of a point is defined as $\mathbf{x}$ and called spatial or Eulerian coordinates.

The motion of the continuum body is described by:

$$\begin{aligned} \mathbf{x}= \mathbf{\varPhi }(\mathbf{X},t) \end{aligned}$$

(24)

where the function $\mathbf{\varPhi }$ maps the reference configuration into the current configuration.

Two different approaches are typically used to describe the response of a continuum: the Lagrangian description and the Eulerian description. In the Lagrangian framework, the material particles are followed in their motion and, in this description, the independent variables are the material coordinates $\mathbf{X}$ at time t. On the contrary, in the Eulerian framework, a fixed region of space where the material particles pass is considered, in this case, the independent variables are the spatial coordinates $\mathbf{x}$ at time t. Typically, the Lagrangian framework is preferred in solid mechanics, while the Eulerian framework is typically used in fluid mechanics. However, the focus of this work is the application of the Lagrangian description to fluid, solid and multi-physical problems.

According to the definition of the reference configuration (i.e. where the equation are written) different Lagrangian approaches can be defined. In the total Lagrangian formulation the reference configuration is the initial configuration and it is fixed in time. In the updated Lagrangian formulation, the reference configuration changes in time and typically corresponds to the last known configuration (i.e. the configuration at the previous time instant or at the previous load increment) [23]. In both formulations, the variables are a function of material coordinates, but derivatives and integration are defined differently. In the total Lagrangian, derivatives are taken with respect to the material coordinates and the integrals (e.g. for the weak form of the equations) are defined over the initial configuration. In the updated Lagrangian method, the derivatives are written with respect to spatial coordinates and the integration is performed over the current configuration.