Review of the modified finite particle method and application to incompressible solids

This paper focuses on the application of the Modified Finite Particle Method (MFPM) on incompressibile elasticity problems. MFPM belongs to the class of meshless methods, nowadays widely investigated due to their characteristics of being totally free of any kind of grid or mesh. This characteristic makes meshless methods potentially useful for the study of large deformations problems and fluid dynamics. In particular, the aim of the work is to compare the results obtained with a simple displacement-based formulation, in the limit of incompressibility, and some formulations proposed in the literature for full incompressibility, where the typical divergence-free constraint is replaced by a different equation, the so-called Pressure Poisson Equation. The obtained results show that the MFPM achieves the expected second-order accuracy on formulation where the equations imposed as constraint satisfies also the original incompressibility equation. Other formulations, differently, do not satisfy the incompressibility constraint, and thus, they are not successfully applicable with the Modified Finite Particle Method.


INTRODUCTION
Meshless methods are a recent methodology in the field of numerical methods whose main characteristic is to be free of any kind of grid or mesh. For this reason, meshless methods can be particularly useful when attacking problems of large deformations and fluid-dynamics, where classical element-based methods may suffer from element distorsion, spurious numerical errors, and mesh sensitivity.
The first meshless method introduced in the literature is the Smoothed Particle Hydrodynamics (SPH) [1,2], where the domain is represented by a set of particles, each one characterized by mass, velocity, and energy. In the SPH the approximation of functions can be interpreted as the projection on the function itself on a bell-shaped kernel function, which simulates the Kronecker Delta distribution. The original SPH formulation, however, presents some inaccuracies, especially at the boundary of the domain, and therefore other SPHderived formulations have been introduced, such as the Reproducing Kernel Particle Method [3,4], the Corrective Smoothing Particle Method [5,6], and the Modified Smoothing Particle Hydrodynamics [7].
In addition to the group of the SPH-derived methods (that have both properties of a meshless and a particle method), some other meshless methods have been introduced in the literature. An important group is the one based on shape functions, that is, unknown fields are reproduced through linear combination of known shape functions. Among shape functions, an important role in the literature has been played by the Radial Basis Functions (RBF), investigated from a theoretical point of view by Kansa [8,9] and applied to elastodynamics [10], wave propagation [11], and fluid dynamics [12].
A further group of meshless methods is the one based on the Taylor series expansion, and they generalize the idea of the Finite Difference Method. This group includes the Generalized Finite Difference Method, the Meshless Least-Square based Difference Method, and the Modified Finite Particle Method.
The Generalized Finite Difference Method (GFDM) aprroximates derivatives starting from the minimization of the error between the value of a function and its Taylor series expansion. GFDM has been applied to the heat equation, Poisson equation, hyperbolic equation [13], advection-diffusion equation [14], and to the dynamics of plates [15,16].
The Meshless Least-Square based Finite Difference Method (LSFDM) [17,18] is similar to the GFDM, with some technical differences based essentially on the selection algorithm of the particles which contribute to the approximation of derivatives in a point.
Finally, the Modified Finite Different Method (MFPM) bases the derivative approximation on the projection of the Taylor series expansion on some known projection functions. MFPM is applied to 1D elasticity problems [19], to 2D Poisson problems [20] and to elasticity and plasticity problems [21]. In Asprone et al. [22] a novel MFPM formulation is introduced, based on the projection of vectors and not of functions, resulting in a reduction of the method computational cost and also in a reduction of the approximation error. The novel formulation has been applied to 2D and 3D elasto-statics and to 2D explicit elasto-dynamics.
In this paper, we apply the Modified Finite Particle Method on incompressible and quasiincompressible elasticity problems. In particular, the displacement-based formulation is investigated in the limit of incompressibility (v → 0.5), and then the Stokes equations for full incompressible solids are investigated. In the field of Finite Difference Method it is well known [23] that the classical discretization of the Stokes Equation on non-staggered grids leads to spurious numerical errors, known as checkerboard instability of pressure. These oscillations are due to the non satisfaction of the inf-sup condition, first studied by Brezzi [24] in the field of Finite Element Method. For this reason, in order to discretize the Stokes problem on non-staggered grids (and then on meshless methods, where staggered grids are not permitted), some different formulations have to be introduced. In particular, the incompressibility constraint equation is replaced by a derived equation, called Pressure Poisson Equation, in which the respect of the inf-sup condition is not requested. However, on this formulation it is not evident which set of boundary condition is needed. A significant contribution to this discussion has been given in [25,26], where the problem is faced using a weak formulation.
The paper is organized as follows: in Section 2, we recall the derivative approximation procedure in the Modified Finite Particle Method; in Section 3 we introduce the equations which describe the statics of solids, first in the compressible form, and then in the limit of incompressibility and, finally, we introduce the Stokes Equations for full incompressibility. In Section 4 we introduce the Poisson Pressure Equation formulation, and the problem of the correct coice of boundary conditions, and in Section 5 we apply the Modified Finite Particle Method on two benchmark problems, using the formulations discussed in Section 4. Finally, in Section 6, we discuss the obtained results and draw some conclusions.

NOVEL MODIFIED FINITE PARTICLE METHOD: MULTIDIMENSIONAL FORMULATION
In the present section we recall the Modified Finite Particle Method (MFPM), as presented in Asprone et al. [22], method that we use for derivative approximation of a function u(x), with x = [x y z] T ∈ Ω ⊂ . The domain Ω is discretized into a set of points X.
For each point x i ∈ X, the approximation procedure considers the Taylor series expansion of u(x) up to the second order, centered in x i : Then, for each x i we select a node subset X i ⊂ X, which serves as support for the derivative approximation in x i . Conceptually X i could coincide with the whole set of nodes X, but the choice of a limited number N i of "supporting nodes" has a beneficial effect on the final computational cost of the method.
Eqn. (1) is then evaluated in the points x j ∈ X i , obtaining in which we neglect higher order terms. The evaluations of the left and right hand sides of Eqn. (2) in the points x j ∈ X i are collected in two vectors p i and q i , respectively. It is important to highlight that at this stage we consider to know the nodal values of u (i.e., u(x i ) and u(x j )), and, therefore, in Eqn. (2) the unknown terms are the derivative evaluations at the point x i . In order to compute the derivative values, we introduce 9 arbitrary .., 9, and evaluate each of them in the points x j ∈ X i . Such evaluations are collected in 9 vectors W i α , α = 1, ..., 9, defined as: Finally, by projecting both p i and q i on W i α , we obtain (4) that can be rewritten in the form Eqns. (5), repeated for α = 1, ..., 9, can be rearranged in matrix form as Now, inverting Eqn. (6) it is possible to compute derivative approximations in x i .

LINEAR ELASTICITY
In the following, we introduce the equations that discribe the equilibrium of an elastic, incompressible body. We first introduce the equations in the classical displacement-based formulations, then we switch to a mixed, displacement-pressure based formulation, in order to fully enforce the incompressibility constraint. In the applications, we will show that the limit to incompressibility of the displacement-based formulation leads to numerical problems, that is the main cause for which the displacement-pressure formulation is introduced.
Review of the modified finite particle method and application to incompressible solids We consider an elastic body on a domain Ω, subjected to internal forces where ρ is the mass density of the material, a is the material acceleration, n is the outward normal vector at the boundary, u = u(x, t) is the vectorial displacement field; σ = (∇u) S is the symmetric Cauchy stress tensor. We use the notation (•) S to denote the symmetric part of a tensor The fourth order linear elastic isotropic tensor  is expressed in index notation as follows where λ and μ are the Lamé constants, which can be expressed in terms of the Young modulus E and the Poisson ratio v: The condition of incompressibility is imposed when the Poisson ratio v is set to 0.5. Unfortunately, when v approaches 0.5, the parameter λ tends to infinity, leading to an ill conditioned discrete system of equations, with consequent degradation of the solution [27]. Therefore, a different formulation is needed where the incompressibility constraint is enforced in a different way.
For an incompressible body, the constitutive relation is modified in the form where p is the pressure, considered, as usual in fluid-dynamics, positive in case of compression. I is the identity tensor, μ is the second Lamé constant and ε is the symmetric part of the gradient of the vector u = u(x). By replacing (10) into the first equation of system (7), we obtain where the incompressibility constraint is introduced. Eqns. (11) and (12) are known as the Stokes equations in primitive variables (u, p), and describe the dynamics of fully incompressible bodies. They have to be completed with suitable boundary conditions, that can be Dirichlet boundary conditions (when the boundary displacement is known), or Neumann boundary conditions (when the boundary traction is known).

SOLUTION OF THE STOKES PROBLEM WITH COLLOCATION METHODS
The discretization of Eqns. (11) and (12), performed using the same spatial discretization for u and p, leads to a well known instability of the pressure field, due to the non satisfaction of the so-called inf-sup condition [28]. This means that alternative formulations have to be introduced in order to overcome this numerical difficulty.
In the Finite Element Method, the classical way to overcome pressure instability is the use of different interpolations for the velocity and pressure fields, the first being discretized using quadratic elements (i.e. six-nodes triangles), while the pressure is discretized using linear interpolation. In this way, the respect of the LBB condition is ensured, and spurious oscillations of the pressure are avoided.
In the field of collocation methods, in particular in the Finite Difference Method, the standard method to satisfy the LBB condition is the use of staggered grids, called also MAC grids [29] (see Figure 1). This kind of grids, however, require rectangular domains and regular node distribuutions, and therefore they are not suitable for meshless methods, where, in general, non regular distributions of points are permitted.
In order to solve the Stokes problem on non-staggered grids, many different formulations have been introduced in the literature [25,26,11,30]. In particular, the 240 Review of the modified finite particle method and application to incompressible solids  previous works are concentrated on whether boundary conditions are required or not, at a discrete level, for the incompressibility equations. In fact the constraint equation holds both on the interior and on the boundary of the domain, and then, no additional boundary condition is required. A reference work regarding this discussion is the one by Sani et al. [26], in which a deep mathematical analysis is done, in the context of the weak formulation. In particular, the analysis is done on the so called Stokes problem with the Poisson Pressure Equation, where the constraint equation of incompressibility is replaced by an equation on the pressure obtained applying the divergence operator on the equations of equilibrium (11).
Restricting to the static case (that is, a = 0) we obtain that is, separating the different components at the left-hand side Eqn. (14) is referred, in [26], as the Consistent Poisson Pressure Equation (CPPE). Changing the order between the Laplacian and divergence operators (that are commutative differential operators) in the term μ∇ ⋅ (Δu) we obtain μΔ(∇ ⋅ u) that is evidently zero due to the incompressibility equation. This permits to simplify Eqn. (14), obtaining the so called Simplified Pressure Poisson Equation.
Δp = ∇ ⋅ b (15) In [26] the discussion is performed in particular on which boundary conditions are required for the solution of the incompressibility problem using the Consistent Pressure Poisson Equation and the Simplified Pressure Poisson Equation. In particular, using the CPPE, no boundary conditions are required for the constraint equation; on the constrary, using the SPPE, a Neumann boundary condition for the pressure is required, obtained projecting the equilibrium equation on the outward normal at the boundary, that is (16) In this paper we solve the Stokes problem using both the Consistent and the Simplified Pressure Poisson Equation, and using the Modified Finite Particle Method to discretize spatial derivatives.
We also solve the incompressibility problem using, instead of boundary condition (16), the discretization of the divergence constrain at the boundary. We refer to this possibility as the SPPE-div formulation.

APPLICATIONS
In this section we apply the Modified Finite Particle Method to discretize the spatial derivatives of the formulations presented in the previous section: the Consistent Pressure Poisson Equation, the Simplified Pressure Poisson Equation (with the boundary condition proposed by [26]), and the SPPE-div. In particular, we test the effectiveness of these formulation on two test cases: an incompressible annulus clamped on all its boundary under a polynomial body load, and an incompressible square under a vertical body load, clamped on two edges. For both problems, we test the MFPM on a displacement-based formulation in the limit of incompressibility (v → 0.5) and then on the mentioned ( ) incompressible formulations. We remark that on the incompressible formulations, the incompressibility constraint is not enforced strongly, but through a derived equation. For this reason we investigate, on both problems, how the incompressibility is respected.

QUARTER OF ANNULUS UNDER POLYNOMIAL BODY LOAD
In this section we study a quarter of annulus, clamped on all its boundary, under a polynomial body load. The geometry of the problem is depicted in Figure 2(a), with R = 4 and r = 1. The problem has been studied in Auricchio et al. [31] using a streamline formulation and isogeometric analysis for the spatial discretization, exploiting the high continuity of the isogeometric shape-functions, and also the possibility of reproducing exactly the geometry of the domain. The analytical solution of this problem is (17) The internal body loads are obtained using the manufactured solution (17). The discretization of this problem using the Consistent Pressure Poisson Equation as constrain (and no boundary condition, as stated in Sani et al. [26]), leads to the non-respect of the incompressibility, as it is shown in Figure 3, where the values of ∇ ⋅ u is shown inside the domain, for a discretization of 10201 nodes. We remark, however, that the demonstration published in Sani et al. [26] et al. had been done using a weak formulation, and then it is possible that, for a collocation method, it is not true.
In Figure 4 we show the performance of MFPM on formulations SPPE and SPPE-div, together with the displacement-based formulations in the limit to incompressibility. All formulations exhibit second-order convergence, with exception of the displacementbased formulation with v = 0.4999, where the solution is affected by the numerical pathology of locking, and the SPPE formulation, where there is convergence, but with a lower order. 242 Review of the modified finite particle method and application to incompressible solids This problem has been solved in Auricchio et al. [31] using the stream-function formulation and an isogeometric approach for the spatial discretization. The second Lamé constant is μ = 40.
Here we solve this problem using a displacement-based formulation in the limit of incompressibility (v = 0.49, v = 0.499, v = 0.4999, v = 0.49999) and using the CPPE, SPPE and SPPE-div.
Using CPPE formulation, similarly to the previous case, the incompressibility is not verified, as it can be seen from Figure 6, where the values assum by the divergence of displacements is shown, using 10201 particles. For what concerns other formulations, the 244 Review of the modified finite particle method and application to incompressible solids  Figure 7 for the vertical displacement of the point B. We notice that the displacement based formulations have results compatible with the numerical problems of the locking; SPPE formulation also shows no convergence, while the SPPE-div formulation shows convergence even faster than the expected second-order. We remark also that in this case is not possible to compute a 2nd norm of the error, since we do not have an analytical solution. We only can compute the relative error in some sampling points, as reported in Auricchio et al. [31].
In Figure 8 we show the deformed configuration obtained with MFPM and a displacementbased formulation (58081 nodes) and v = 0.4999. A comparison with Figure 9, in which an overkilled deformed structure is shown, highlights that the displacement-based methods, in the limit of incompressibility, suffer from volumetric locking.

CONCLUSIONS
In the present paper we applied the Modified Finite Particle Method on some incompressible elasticity problems. In particular, some different formulations have been investigated: a displacement-based formulation, in the limit of incompressibility, with v → 0.5, and three different formulations of the Stokes problem. For these formulations, in particular, the incompressibility constraint (∇ ⋅ u) is not imposed strongly, but it is replaced by a derived one, in which the Lapacian operator is applied to pressure. This choice is done to overcome the difficulties related to the non-respect of the inf-sup condition, which results in unphysical oscillations of the pressure field.
Unfortunately, these derived formulation may need some boundary conditions for the constrain equation, that are not needed by the original Stokes problem in the divergence form. In this paper we investigate three different fully-incompressible formulations: the so Consistent Pressure Poisson Equation, the Simplified Pressure Poisson Equation and the Simplified Pressure Poisson Equation with the divergence constraint at the boundary. These formulations differ among each other for the boundary conditions imposed on the constraint equation: in the first case, according to Sani et al. [26], no boundary conditions are required; in the second case, where the normal component of the displacement is know, a boundary condition for the pressure is obtained projeting the equalibrium equation on the outward normal; in the SPPE-div, instead, the divergence-free constrain is applied as boundary condition for the Pressure Poisson equation on the whole boundary.
Here we see that for the CPPE formulation, the incompressibility constrain is not respected for both problems under investigation; the SPPE exhibits lower convergence of the error with respect to the expected second order; and finally the SPPE-div formulation exhibits correct second-order accuracy, even if with a high constant of the error. The displacement-based formulations, even if correctly discretized with MFPM, exhibits the numerical pathology of locking.