Comparison of Material Point Method and Finite Element Method for Post-Failure Large Deformation Geotechnical Analysis

Finite Element Method (FEM) has been the state-of-the-art method in the geotechnical analysis since it was first formulated in the 40s. This method can handle Multiphysics simulation, soil-structure, soil-water interaction, and time history analysis. Though powerful, the standard Lagrangian FEM suffers from mesh distortion when handling large strain deformation problems. This mesh entan-glement problem makes modelling post-failure analysis considerably challenging, if not impossible, using FEM. Therefore, the Material Point Method (MPM) was introduced to solve these large strain deformation problems. Adapted from the Particle in Cell (PIC) method, MPM is a hybrid method that combines the Eulerian and Lagrangian approaches by utilizing moving material points that are moving over spatially fixed computational mesh. The method’s features enables it to calculate not only fluid mechanics such as in PIC but also solid mechanics and their intermediatory states. To demonstrate the capability of this method and its consistency with FEM in geotechnical analysis, this article presents a comparison of FEM and MPM analysis on a hypothetical slope using the Mohr-Coulomb constitutive model. The results obtained from the simulation showed that both FEM and MPM are consistent with each other, specifically in a small strain scheme. However, in large strain deformation, MPM was able to obtain convergent results while FEM could not. The MPM simulation was also able to animate post-failure behavior, calculate post-failure strains and stress distribution, and present the model’s final geometry.


INTRODUCTION
In geotechnical engineering, problems encountered are mostly associated with soil deformation. One common problem that holds catastrophic danger is slope instability which occurs in both natural and man-made environments. To address this imminent danger, engineers developed various analytical methods to predict the possibility of such an event occurring. Currently, most slope stability analyses involve failure prediction and safety factor determination. Although such analysis could help minimize the possibility of the event occuring, it does not prevent it entirely. It is thus, also very important to understand the postfailure soil behavior, indicating further damage could be prevented.
One popular method that has been used for this purpose over a long period is the limit equilibrium method (LEM). The notion consists of equi-librium analysis of the forces around the assumed failure plane (i.e., circular and planar failure surfaces). Other notable methods that emerged include Bishop, Janbu, and Fellenius among others.
As a force-based method, LEM is simple and userfriendly, but it cannot calculate the basic physics of the stress-strain relationship. Therefore, any deformation that occurs in the model is practically ignored and any problems associated with stress concentration and deformation compatibility cannot be solved using this method.
Another method that has been used extensively in the past decades is the Finite Element Method (FEM). While it is difficult to pinpoint the first time this method was formulated, it has been in use since the early 1940s. Presently, FEM can be considered the state-of-the-art method in geotechnical engineering analysis. This is because of the method's ability to work in a continuum domain and its capability to model history-dependent material behavior such as soil. In FEM, the continuous domain is discretized into smaller elements before they are solved. Also, this method incorporates a weak formulation of momentum balance laws, reducing the order of the approximation function needed to calculate element behavior.
The Conventional FEM algorithm uses Lagrangian mesh to keep track of the material domain. The mesh has well-defined free surfaces and can accurately model deformation and this feature is most favorable in modeling solid mechanics. However, soil behaviors can differ from normal solid materials due to their granular nature. Granular materials have peculiar properties. These materials strength is determined by their ability to withstand shear stresses. When these shear stresses exceed the material's capacity, it will deform extensively along with large shear strain. Such conditions limits the effectiveness of FEM considering the fact that under large deformation, the Lagrangian mesh could be so distorted that accuracy is lost. Figure 1 shows a mesh distortion in slope stability analysis. The mesh geometry, as shown in the figure, has been extremely deformed to a point that it practically resembles a line, whereas its original shape was a triangle. Typically, proprietary FEM software will stop calculating when such condition is met since the results will just diverge.
The way FEM behaves under large strain deformation is a problem in post-failure analysis because conclusions cannot be taken from diverging results. Attempts have been made to solve this problem using various re-meshing techniques, but there is still no optimal procedure currently available for such function. Forcing to re-mesh only creates an exponential increase in computational cost, which is undesirable.
Furthermore, to deal with large deformations, variations of FEM such as Eulerian FEM were formulated. The Eulerian FEM uses spatially fixed mesh to track material deformation and this solution directly solves the mesh distortion problem. However, because the computational mesh is now separated from the material, convective terms appear in the equations and create more complexity due to their matrix non-symmetrical properties.
In practice, the Eulerian FEM is more often coupled with standard Lagrangian FEM. It is called a Coupled Eulerian-Lagrangian (CEL) FEM. This approach presents the benefit of Eulerian Mesh strictly to the high strain region of the model, theteby making the computation efficient. Although its success in modeling large deformation, CEL is mostly effective in modeling solid and fluid interaction in which both material behavior is presumably consistent in the simulation. With soil, things are different since the behavior could be continually changing in the simulation. For details of CEL implementations in the geotechnical engineering field, readers are referred to (Henke et al., 2010).
Following this, other than mesh-based methods such as FEM, meshless or mesh-free methods have also been developed to accommodate large strain deformation problems. These methods are referred to as meshless because the space is discretized into several material points/particles instead of meshed elements as in FEM. These material points, or particles, interacts with each other in a relatively flexible manner.
One of the oldest meshless methods developed is the Smoothed Particle Hydrodynamics (SPH). This method originated from astrophysicists who were studying dust clouds and star explosions (Lucy, 1977;Pastor et al., 2009), in their study, presented a depth-integrated, coupled model discretized using SPH. This model was then applicated in a landslide run-out analysis which was conducted by Pastor et al. (2014). Furthermore, another notable meshless method is the Material Point Method (MPM). According to Sulsky et al. (1994), this method was originally an extension of the Particle-in-Cell method which was developed by Harlow (1964) to facilitate solid mechanics. In this study, MPM and conventional FEM were compared, specifically for post-failure large deformation analysis.

Finite Element Method (FEM)
Similar with any other numerical method, the principle of FEM is rooted in the idea of discretizing a continuum body into smaller elements. These elements, which are depicted in Figure 2, are assumed to behave linearly. Furthermore, in computational mechanics, the Equation 1 that governs all material behavior is the conservation law of mass, momentum, and thermodynamics. When neither mass nor heat is leaving and/or entering the system, however, these laws could be narrowed down to the conservation law of momentum alone. This is convenient because in terms of FEM analysis, the total mass in the system is constant and heat flux is ignored. The governing Equation 1 in a Lagrangian representation can thus be expressed as: Where ( ⃗ dv)/dt is the vector of acceleration,σ is the Cauchy's stress tensor, and ⃗ g is the gravity vector. Notice that in the standard momentum conservation there exist a second order term on the lefthand side. The presence of this term causes problems because the approximation function needed to solve the Equation 1 must be differentiable, at Subsequently, since the element used in FEM consists of several nodes, a shape function can be implemented to interpolate values from a discrete nodal value to any point within the elements. In order to obtain a precise value, the values calculated with the shape function were located on the gaussian points of each element, as shown in Figure 2. Subsequently, there was a divergence in the implementations as various people try to optimize the solution and improve the FEM performance. For further implementation details, readers are referred to PLAXIS (2021).

MPM is in a nutshell is an extension of FEM.
Their only difference is that in MPM, discretization comes in the form of material points (MPs). Figure 3 shows the material point discretization scheme in MPM. These MPs carry all the information, including stresses, displacements, and velocities, regarding the Lagrangian description of motion. While the points were carrying all the information, the momentum balance equations calculation was carried out in the background of the Eulerian computational mesh. Since the information is stored in the material points, the calculated value in the nodes can be discarded after each calculation step. To simplify understanding, these MPs can be compared to the gaussian integration points in FEM. Unlike gaussian points, these MPs are able to move in the computational mesh and cross element boundaries freely. Thus, leveraging the benefits of both the mesh-based and mesh-free methods. The complete MPM for- mulation and its implementation are discussed by Kafaji (2013) rigorously.
The MPM calculation procedure for explicit onestep time increment is depicted in Figure 4. The cycle starts with mapping information from MPs to the nodes of the computational mesh using shape functions. Similar to FEM, the mass in each MP is strictly constant, hence, the mass balance Equation 1 is automatically satisfied. Therefore, in the next step, only momentum balance equations are solved in each node to obtain primary unknown variables (i.e., nodal accelerations). These nodal values are then remapped into each MPs to further update their information. At this point, the nodal values can be discarded and prepared for the next cycle. Finally, the updated MPs information is then used to determine the new locations, stresses, strains, etc of the points.
Subsequently, the usage of Eulerian computational mesh in MPM automatically eliminates the problem of mesh entanglement. It also makes treating fictional contacts between bodies straightforward and efficient (Hamad et al., 2017). Because most MPM algorithms were derived from FEM, transitioning to MPM seems more natural than other meshless methods. However, considering the fact that no method is by any means perfect, some drawbacks need to be addressed. They include cell-crossing instability, relatively higher computational cost than FEM, and enforcement of boundary conditions (de Vaucorbeil et al., 2020).

Model Configuration
To perform the comparison study, a homogenous 2D clay slope with geometry, as shown in Figure  5, was selected and simulated using both FEM and MPM. The slope and the vertical rolling bountry rested on a fixed boundary condition on the bottom and on each side of the model, respectively. Furthermore, to facilitate the possibility of deepseated landslides, a 5 m deep and 25 m wide layer of base soil was provided. The slope had a 4.5-m height and 26.6°inclination.
The FEM simulation in this study was performed using PLAXIS and that of the MPM was performed using codes provided by the CB-Geo computational geomechanics research group (Kumar et al., 2019). In PLAXIS , the domain was discretized using triangular mesh automatically with a "Fine" element distribution. The mesh produced by the computer program had 15 nodes by default and this implies smoother interpolation capabilities.
Meanwhile, in CB-Geo MPM, the domain meshed with uniform quadrilateral elements with 4-nodes, and size similar to that of PLAXIS's mesh. This produced 0.5 cm square elements in the discretized domain. After which the MPs were injected into the domain in such a way that each element/cell can have at most 16 MP. The configurations for the FEM and MPM analysis are listed in Tables 1 and 2, respectively.
Following this, to successfully create larger deformation and comparable results between the two methods, an initial load of 15 kPa was added to the top of the slope after the initial phase was calculated. The load was then incremented by 15 kPa after each step until the calculation results diverge. This load was selected in accordance with the study conducted by Zhou and Sun (2020), where it was assumed that 15 kPa is equivalent to a load of a story building, hence, the incremental load was applied to generate large deformations. In addition, stressed points from 3 areas (top, middle, and bottom) were selected to monitor the deformation which had occurred due to the loading scheme. The areas selected were consecutively located at 0 m, -3 m, and -6 m from the center of the load into the soil.  The examples shown in this study used Mohr-Coulomb constitutive model which is widely used in geotechnical engineering practice. The model is simple and relatively easy to use. There are 2 parameters used in the model to calculate the soil shear strength: the cohesion c' and the friction angle φ'. The relationship between the 2 parameters can be expressed as: Where τ is the shear stress and σ is the normal stress. Furthermore, the input parameters for the Mohr-Coulomb soil model are listed in Table  3. These parameters were used to simulate loose sand with N SPT of 8. The cohesion, however, was determined at 5 kPa to maintain the numerical calculation consistency and prevent failure prematurely.
In terms of the model parameter, there was no significant difference between FEM and MPM. Both models used the Mohr-Coulomb model to simulate soil behavior. However, unlike FEM whose basic algorithm calculates problems implicitly, the MPM algorithm calculates problems explicitly. It is required that the user should thus provide a magnitude of time step that MPM can work with. Although this could be any arbitrary small number, it is recommended to use a critical time step such as de Béjar and Danielson (2016) proposed for explicit integration of dynamic higher-order FEM to produce an efficient running time and an economic analysis.
Another thing to consider is the number of particles per element. This could be opted based on the ideal number of gaussian points for integration. In a 2D element, for example, the number could be from 4 MPs/cell up to 16 MPs/cell or even more depending on the fineness expected in the model.

Simulation of Small-Strain Deformation
The results of the finite element and material point simulations are shown in Figure 6. As expected, the finite element analysis had no problems calculating displacements and maintaining a convergent result in small strain phases. The first row of the snapshots was taken at stage 4, in which the slope was given a load of 60 kPa. This was where the figure started to showing noticeable deformation on the top surface of the slope. The next row shows the slope's condition at stage 7, in which another 45 kPa was added to the initial load. At this stage, some areas of the slope had reached a critical state and had become plastic. Finally, the last row shows the slope's condition at stage 10, which was the stage before PLAXIS started feedbacking calculation errors due to the soil collapsing. The plastic failure distribution in this stage is shown in Figure 7 which indicated that there was a most likely punching shear failure occurring.
MPM, as shown in the right column of Figure 6, provided a seemingly consistent result compared to FEM analysis. This result is further supported by Figure 8 which presents the comparison between volumetric strains in all 3 control points (i.e., top, middle, and bottom) tracked in each stage of loading, with the top point located di- Subsequently, the strain value discrepancies calculated with FEM and MPM could be addressed according to the shape of the mesh used (i.e., quadrilateral and triangular). The mesh shale could cause differences in the failure planes predicted which explains the location in which both methods' results diverge. Accordingly, the problem with mesh shapes needs to be further investigated to properly understand how its impact differs quantifiably. The difference could also be addressed with regard to the fact that the material points, which act as gaussian points in MPM, are generally not in the optimal position for numerical integration. Moreover, in some locations, MPM shows sudden spikes of strains or stresses which usually occur due to cell-crossing noise. Such errors are common in quasi-static problems, in which the rate of loading is very slow, hence, the effect of inertia becomes insignificant. In this sit- uation, Bardenhagen and Kober (2004) introduced the Generalized Interpolation MPM (GIMPM) to mitigate this effect. This interpolation method was then further improved by Sadeghirad et al. (2011) into what is now known as Convected Particle Domain Interpolation (CPDI) MPM.

Simulation of Large-Strain Deformation
The load applied to the slope only slightly deformed it. However, the advantages of MPM are more apparent in large deformation settings. To achieve this condition, the slope received an exceptional load of 390 kPa at stage 10 as opposed to the prior 150 kPa.
Furthermore, to create a bigger picture of the deformation, the model was first to reset to its initial condition. The 390 kPa load was added all together to the model. In order to also prevent material oscillation due to sudden loading, the load was applied gradually in a linear manner with a time interval of 3s, and the time step used to calculate this simulation was reduced to 0.0001s to prevent the occurrence of numerical errors and create smoother results. Figure 9 shows the simulation result in the form of snapshots. The color contour in these snapshots represents the strain magnitude which was calculated by Equation 3.
Where,ε is the total strain tensor calculated from the constitutive Equation 3. The first snapshot was taken at t = 5s in which the load was just fully applied to the model. The second snapshot was taken at t = 25s. At this time, some areas near the foot of the load had already deformed significantly and the compression stress from the load created a wedge zone which was also apparent in Terzaghi's general failure plane assumption. Accordingly, the third snapshot was taken at t = 50s and this was when the failure plane became apparent. The fourth was taken at t = 75s which shows the sliding of the landmass over the failure plane. Lastly, the fifth shot was taken at t = 100s. At this stage, the final deposition profile of the failure was depicted.
This extreme load simulation indicated the occurrence of a sliding failure. The material's final deposition shows that the soil will settle at an elevation of 8.3 m from the bottom boundary. In a real case scenario, these simulation results could be useful, for example Zhou and Sun (2020) used a combination of MPM and Monte Carlo simulation to determine landslide risk quantitatively.