MPM vs. CEL: Numerical modelling of penetration processes

Abstract

The simulation of large deformation problems in geotechnical engineering has gained increasing interest in recent years. This is partly due to new technologies and industry requirements (e.g. the installation of new foundations in the offshore industry), the improved availability of numerical methods that can handle large deformations and the overall improvement in computational hardware.
 A frequently adopted numerical scheme is the Coupled-Eulerian-Lagrangian (CEL) method [1,2,3], which is available in the commercial software Abaqus. Originally developed for solid-fluid interactions, it is now also used in geotechnics, modelling the soil by Eulerian finite elements. These elements decouple mesh and material movement, allowing for arbitrary large deformations. However, since elements can be partially filled by material, special considerations for integration and reconstruction of material surfaces within elements are necessary. A one point integration rule is adopted and the volume-weighted fraction of the material is used in case of integration of partially filled elements.
 Another approach is adoped in the Material-Point-Method (MPM) [4], where the the continuum body is represented by a cloud of points, called material points (MPs). They carry all information of the continuum body, such as density, velocity, strain, stress, material parameters and external loads. They are not physical particles, like for example single solid grains described in Discrete Element Method, but they represent a portion of continuum body. At the beginning of each timestep, the nodal forces are computed based on the information stored at the material points by means of the shape functions. The governing equations are then solved at the nodes and afterwards used to compute the strain, stress and density, and to update the position of the material points. At the end of the timestep, the mesh is reset. For penetration problems, the “moving mesh” concept is adopted, where the mesh moves together with the structural element to keep the contact surface on the same computational nodes shared between the soil and the penetrating object, ensuring higher accuracy. For the MPM formulation used in this work, the mean stress and the state variables are averaged in each element, which can be interpreted as a “stress-smoothing” algorithm.
 Up to now, no comprehensive study comparing the two methods exists. In this study, both schemes are applied to the simulation of cone penetration testing (CPT). An explicit time integration scheme is adopted for both methods.
 The field of deviatoric stress (stress due to self-weight of the soil not accounted for) in a simulation of a CPT using the CEL method in conjunction with the Tresca model is depicted in Fig. 1a). Note that undrained conditions are assumed and the undrained shear strength su is constant in the model. Note in addition that the built-in Mohr-Coulomb model with a friction angle of zero is used for the simulation with Abaqus. Undrained conditions are realized by setting the Poisson’s ratio to 0.495. Figure 1b) shows the average vertical stress at the tip of the CPT during penetration for simulations using MPM and the CEL method for perfectly smooth surface conditions between pile and soil. A spatially constant effective vertical stress of 66 kPa is included in the value of qc. While the trend of the curves is similar for both methods, the CEL method tends to give more unstable results, manifesting in strong oscillations of the stresses. In comparison, the MPM provides a nearly constant stress curve for greater depths, which would be expected for constant values of su. To confirm these findings, additional simulations with variation of mesh and time increment size where performed in case of the CEL method. Previous research [5] has shown that both can have an influence on the results of CEL analyses. However, the observed oscillations could only be marginally reduced by using a finer mesh and a smaller time increment size.
 The comparison of the cone factor Nc, defined as the ratio of the cone resistance (excluding initial stress due to self-weight) to su, for different values of shear stiffness G to su for simulations using MPM and the CEL method is given in Fig. 1c). Perfectly undrained conditions and a smooth contact are assumed. For the evaluation, the cone resistance at the end of the simulation, i.e. at a penetration depth of 0.4 m, is used. For the entire range of G/su, good agreement is found between the two methods.
 Overall, it is concluded that both numerical schemes give similar results for typical benchmark problems, but the CEL method tends not to be as numerically stable as the MPM. This is somewhat surprising considering that the MPM code (a version of Anura3D developed at Deltares is used) is developed by a comparably small community of reseachers, while the CEL method is implemented in the framework of a massively commercial code. However, the superiority of the MPM is in parts due to the adoption of axialsymmetric elements, which is not possible using the CEL method implemented in Abaqus (even though one could model only a narrow slice instead of the quarter model used here). This is also why the MPM simulations are computationally more efficient in the present case. More comparisons of the two methods will follow in future work.