Overview of Material Point Method

Notations

Name of the methods

MPM; Material Point Method (Sulsky et al 1994) [3]

GIMP; Generalized Interpolation Material Point Method (Bardenhagen and Kober 2004) [4]

Bspline MPM: Bspline MPM (Steffen et al. 2008) [9]

WLS: Weighted Least Square (Wallstedt and Guilkey 2009) [14]

CPDI; Convected Particle Domain Interpolation (Sadeghiad et al. 2011) [12]

DDMP: Dual Domain Material Point Method (Zhang et al. 2011) [10]

CPDI2: Convected Particle Domain Interpolation 2 (Sadeghirad et al. 2013) [13]

iMPM: Improving Material Point Method (Sulsky and Gong 2016) [15]

CMPM. Composite Material Point Method (Acosta et al. 2017) [11]

iMLS: Improved Moving Least Squares MPM (Tran et al. 2019) [17]

TLS: Taylor Least Squares MPM (Wobbes et al. 2019) [16]

CPLS: Convected Particle Least Squares Interpolation MPM

The first category is to consider the material point as a discrete point. In this path, the MPM development trend is to employ higher order grid basis function rather than linear basis function such as B-spline functions [9], higher order gradient of the shape function [10] or using Lagrange basis function [11].

The second category is to consider the material point as a domain. In this path, the MPM development trend is to improve the evolvement of the material point domain while still using the linear grid basis function. For example, Bardenhagen and Kober [4] introduced the undeformed GIMP, the contiguous particle GIMP which the material point domain is undeformed or deformed only in axial directions. Convected Particle Domain Interpolation [12] [13], on the other hand, allows the material point domain to capture the shear deformation.

In addition, the function reconstruction can be used to improve the accuracy of the MPM. In this trend, Wallstedt and Guilkey [14] firstly introduced the function reconstruction by combining GIMP and weight least squares method. Latter, Sulsky and Gong [15] proposed the moving least squares for MPM; Wobbes et al. [16] presented the Taylor least squares for DDMP and B-spline MPM; Tran et al. [17] introduced the improved moving least squares for MPM and combined CPDI and the improved moving least squares to improve the accuracy and stability of the function reconstruction method, refer as Convected Particle Least Squares Interpolation MPM.

Application of Material Point Method

Granular materials

Material Point Method has been applied to study the granular materials since 1999 by Wieckowski et al. [18] and Bardenhagen et al. [19]. Other researches investigated the granular behaviour under silo discharge [20], [21]; bucket filling [22]; impact force of granular material into the structure [23] [24] [25] and granular column collapses [26] [27].

Geomechanics and GeoEngineering

Material point method has been adopted in large deformation in geomechanics and geo-engineering including landslides [28] [29] and for back analysis of the case study including Aznalcollar dam [30], Hongshiyan landslide [31], Oso landslide [32], Sainte-Monique landslide [33] or to simulate the slab failure of snow avalanches [34], [35].

Penetration problem is a typical problem in geoengineering. Material Point Method has been used to study the penetration of foundation [26], fall cone tests [33] [36], quickness tests [37] and cone penetration [38] [39].

Membrane

In 1997, York II [40] modified the material point method to simulate the thin membranes and studied the fluid-membrane interaction [41]. Ionescu et al. [42] apply membrane MPM to study the soft tissues during penetrating trauma. MPM was also used to model the behaviour of geomembranes [43] [44] and the thin-walled tubes [45].

Implicit MPM

First implementation of the implicit time integration for MPM was presented by Cummins and Brackbill [46] in 2002. In this research, the matrix-free Newton Krylov was used for the implicit integration, similar to the work of Sulsky and Kaul [47]. Later, Love and Sulsky [48] [49] used consistent mass matrix with implicit MPM to achieve the conservation of energy for MPM.

Guilkey and Weiss [50] applied the Newton’s method to solve the increment nodal displacement in a linear system with a tangent stiffness matrix. This approach allows larger time step (limit by a strain of 50% per time step) than the matrix-free Newton Krylov (limit by a strain of 1% per time step). Similar formulation was developed by Wang et al. [51] for geotechnical applications. Other works used the quasi-static formulation with the implicit solver such as Beuth et al. [52] [53] and implicit GIMP by Charlton et al. [54].

Fracture, crack and material failure in MPM

Material Point Method was exploited to model the fracture and crack propagation by Tan and Nairn [55], Gilabert et al. [56], Nairn [57], Wang et al. [58] or using the cohesive zone models by Daphalapurkar et al. [59]; Guiamatsia and Nguyen [60]; Bardenhagen et al. [61]. Other approaches can be used to model cracking such as using kernel‐based damage field model by Homel and Herbold [62], [63], [64]; using enrichment by Homel et al. [65], Liang et al. [66] or using phase-field model by Kakouris and Triantafyllou [67], [68].

Contact algorithm in MPM

Although MPM algorithm can handle automatically the non-slip contact, many engineering problems often require different physical contact interactions between bodies. Bardenhagen et al. [19] introduced the friction contact for MPM in 2000 and this formulation is commonly used in literature. It was further improved and applied to stress propagation of granular materials [69]. The friction contact is extended to implicit MPM [70] or friction contact considering the heating [71]. Other contact algorithms were proposed based on multi-mesh [72] [73] and multi-material method [74]; or different numerical treatment [75]; or for the adhesive contact [76].

Multi-phase coupling

MPM can be coupled with hydro-mechanical formulation to solve the solid-fluid interaction problems. The one-layer material point formulation uses only a single set of material points to represent both the solid and liquid phases. For example, Zhang et al. [77] proposed a velocity - pressure (v-p) formulation [78] in the one-layer MPM framework and used it to solve a problem involving impact of solid bodies into a saturated porous media. Similarly, Zabala and Alonso [31] applied a hydro-mechanical MPM formulation with strain softening Mohr-Coulomb to model the progressive failure of Aznalcollar dam. Later, Jassim et al. [79] applied a velocity - velocity (v-w) hydro-mechanically coupled MPM formulation [80] and modelled the wave attack on a sea dike. Compared with the v-p scheme, the velocity v-w scheme brings several advantages such as (i) simpler imposing of impervious boundary condition, (ii) no need for zero pressure detection and (iii) being suitable for large deformation at failure when the separations between liquid and solid phases may occur. Other works presents different strategies for hydro-mechanical coupling for MPM such as using dragging interactions by Mackenzie-Helnwein et al. [81]; combining with CPDI by Zheng et al. [82]; combining with GIMP by Liu et al. [83] and Homel et al. [84]; combining with DDMP by Tran and Solowski [85], [86].

Latter, the two-phase formulation was extended into three-phase formulation for modelling the unsaturated soil by using three balances equations (soil-liquid-air) by Yerro et al. [87] or using two balances equations (solid-liquid) by Bandara et al. [88] and Wang et al. [89].

Other works studied the thermal coupling using GIMP such as by Tao et al. [90]; or using MPM by Pinyol et al. [91] or chemical coupling using CPDI by Grittion et al. [92].

Apart from the one-layer formulation for multi-phase coupling, another alternative way to model the soil-water interaction is to use a two-layer Material Point Method formulation. The two-layer formulation provides opportunities to solve the dynamic flow problems in geomechanics such as internal and external erosion as well as the fluid flow failure (seepage failure). The two-layer approach was adopted by Abe et al. [93] and Bandara and Soga [94] to model the seepage failure of a river bank embankment or for solid-fluid interaction in the animation [95].

Multi-scale coupling

Different strategies for multi-scale coupling have been investigated including coupling MPM with molecular dynamics [96], [97], [98], [99], [100]; coupling MPM with dissipative particle dynamics [101] or coupling MPM with Discrete Element method in a hierarchy approach [102] or a surface approach [103].

Accuracy and Stability of MPM

Cell-crossing errors

Cell-crossing errors of original MPM was pointed out by Zhou et al. [43] in 1999. That error is caused by the discontinuous gradient of the shape functions. The discontinuity in the shape function gradient results in a sudden change of the stress when a material point crosses to a new cell. It can be mitigated by using a smooth gradient of the shape function such as in GIMP, DDMP, B-spline.

Quadrature errors

Steffen et al. [104] indicates that quadrature errors in MPM are due to the first-order integration of the field quantities. Therefore, MPM integration is less accurate than FEM which use the Gauss integration. The quadrature errors reduce the accuracy, leading to the reduction of the spatial convergence rate for the numerical solution. Steffen et al. [104] also demonstrate that this error can be reduced by using higher order shape function such as GIMP and B-spline.

Conservation errors

Ideally in a physic-based simulation, three quantities should be conserved including linear momentum, angular momentum and energy. In MPM algorithm, the linear momentum is conserved in the transfer between material points and grid nodes as well as in the balance equations. The former is satisfied because of the partition of unity of the shape function and the latter is enforced in governing equation. The other quantities (angular momentum and energy) are not always conserved for MPM.

Bardenhagen [105] investigated the energy conservation errors with different stress update schemes for MPM including Update-Stress-Last (USL) and Update-Stress-First (USF). The former updates the stress in the end of the time step while the latter updates the stress at the beginning of the time step. Bardenhagen showed that the USL is energy-dissipative while USF is energy-conservative for the vibration example. However, he suggested the USL and it is commonly used in the MPM community. Other schemes including modified USL (MUSL) presented by Sulsky [7]. This method mapped particle velocity to the grid nodes before computing the particle velocity gradient for stress computation. Nairn [57] proposed the Update-Stress-Average (USAVG) as the stress updating is weighted from USF and USL.

To conserve the angular momentum, Jiang et al. [106], [107] proposed a locally affine transfer and an implicit solver in the framework of PIC. The new method, refer as affine PIC, showed a conservation of linear and angular momentum but no energy conservation due to the dissipative properties of PIC.

Love and Sulsky [48], [49] demonstrated a MPM formulation to achieve an exact conservation of linear momentum, angular momentum and energy by using a implicit solver and a consistent mass matrix. Although a full conservation is obtained, this approach is not preferable in practice as the consistent mass matrix may be singular for certain material point configurations.

Null space errors (Ringing instability)

Null space errors (ringing instability) was firstly introduced by Brackbill [108] for PIC in 1988. The null space errors occur because of a mismatch between the material points and grid nodes as the Fourier modes in material points with wavelengths shorter than the grid spacing are indistinguishable at the grid nodes. This error was further examined and eliminated by a null space filter proposed by Gritton et al. [109], [110] for MPM. Other works also suggested to reduce the ringing instability and null space errors by a PIC interpolation (Jiang et al. [107, 106]), or using a Moore-Penrose inverse mapping in MPM [111]. In Chapter 4, the null space errors are discussed, and a new null space filter is proposed.

Non-linear instability

In explicit MPM formulation, the Courant-Friedrichs-Lewy (CFL) condition is commonly used to determine the critical time step for numerical solutions. However, CFL condition, which was used for the linear stability analysis, does not always ensure a stable non-linear solution for MPM. Some detailed investigations were performed by Brannon [112] who demonstrate the unstable solution in MPM, refer it as kinematic anomaly. Indeed, Jiang et al. [107] selected the critical time step considering the maximum velocity of material points. Berzins [113] proposed a critical time step using the non-linear stability analysis for MPM and GIMP.

Volumetric locking and shear locking

Some instabilities, which inherent in other numerical methods such as in Finite Element Method, also occurs in MPM including volumetric locking for incompressible materials and shear locking. Different anti-locking strategies can be used such as by Mast et al. [114] for MPM, by Yang et al. [115] for water-solid interaction in MPM, by Coombs et al. [116] using average deformation gradient technique.

QUOC ANH

APRIL - 2019