Development of an implicit material point method for geotechnical applications

An implicit material point method (MPM), a variant of the ﬁnite element method (FEM), is presented in this paper. The key feature of MPM is that the spatial discretisation uses a set of material points, which are allowed to move freely through the background mesh. All history-dependent variables are tracked on the material points and these material points are used as integration points similar to the Gaussian points. A mapping and re-mapping algorithm is employed, to allow the state variables and other information to be mapped back and forth between the material points and background mesh nodes during an analysis. In contrast to an explicit time integration scheme utilised by most researchers, an implicit time integration scheme has been utilised here. The advantages of such an approach are twofold: ﬁrstly, it addresses the limitation of the time step size inherent in explicit integration schemes, thereby potentially saving signiﬁcant computational costs for certain types of problems; secondly, it enables an improved algorithm accuracy, which is important for some constitutive behaviours, such as elasto-plasticity. The main purpose of this paper is to provide a uniﬁed MPM framework, in which both quasi-static and dynamic analyses can be solved, and to demonstrate the model behaviour. The implementation closely follows standard FEM approaches, where possible, to allow easy conversion of other FEM codes. Newton’s method is used to solve the equation of motion for both cases, while the formation of the mass matrix and the required updating of the kinematic variables are unique to the dynamic analysis. Comparisons with an Updated Lagrangian FEM and an explicit MPM code are made with respect to the algorithmic accuracy and time step size in a couple of representative examples, which helps to illustrate the relative performance and advantages of the implicit MPM. A geotechnical application is then considered, illustrating the capabilities of the proposed method when applied in the geotechnical ﬁeld. (cid:1) 2015 The Authors. Published by Elsevier Ltd. ThisisanopenaccessarticleundertheCCBYlicense(http:// creativecommons.org/licenses/by/4.0/).


Introduction
The material point method (MPM) has been shown to be a robust spatial discretisation method for simulating multi-phase interactions involving large deformations and failure evolution.During 1994-96, Sulsky et al. [1][2][3] first developed and applied the method for modelling solid materials.This led to researchers, from different fields, recognising the potential of the method and adapting it to various applications, e.g.silo discharge and plastic forming [4,5], explosion problems, exploiting its ability to represent an arbitrary geometry [6,7], large-scale response of cellular constructs in biomechanics [8], and, more recently, for geotechnical analysis, including the modelling of retaining wall failure [4], anchor pull-out [9], soil column collapse [10,11], landslides and debris flows [12], landslide induced interactions with structures [13], and quasi-static analyses of slope stability [14,15].
MPM uses two spatial discretisations, the first one that discretises a continuum body with a set of material points carrying all the state variables, and the second one that discretises the background grid (a computational mesh) to solve the equations of motion.The computational mesh may be maintained in its original position, or it can be adjusted in an appropriate way to avoid mesh distortion after each time/loading step, thereby removing the disadvantage of the finite element method (FEM) for which extreme mesh distortion may occur due to large deformations.As with FEM, the time integration of MPM can be either implicit or explicit, in which the latter has been employed for most of the MPM developments so far.This paper is concerned with the implementation of an implicit MPM framework, and the validation of the resulting implicit solver, as well as with comparisons with an explicit MPM with respect to the time step size and accuracy of the results.In this paper, the term implicit MPM refers to a framework where both dynamic and quasi-static problems (with inertial terms neglected) can be solved effectively.Although implicit dynamic MPM formulations [16,17] have been reported, this paper aims to provide a clear and straightforward description of all the necessary techniques for adapting an existing FEM implementation into one based on the implicit MPM.
In the remaining sections of the paper, the theoretical formulation is first presented.This closely follows the standard FEM procedure, thereby clearly demonstrating the similarities between MPM and FEM.Implementation details are then discussed, where a special treatment for MPM is needed.The subsequent section focuses on a series of representative examples to investigate and validate the presented framework for quasi-static and dynamic analyses, respectively, with the results being compared with those obtained from an explicit code, in order to gain a thorough understanding of how the implicit algorithm behaves.

General framework
To describe the implicit MPM, Fig. 1 demonstrates the standard mapping and remapping procedure between the material points and background computational mesh.In the first phase (Fig. 1(a)) the state variables are mapped from the material points to the nodes of the background mesh; in the second phase (Fig. 1(b)), the equation of motion is solved over the background mesh to find the current acceleration, with the element integration being based on the material points (rather than on the information mapped to Gauss points); and, in the third phase (Fig. 1(c)), the state variables on the material points are updated via remapping from the deformed background mesh, and the mesh is then reset, leaving the material points at their updated locations.These phases are repeated until the end of the time/loading steps.
Connectivity can be set up between the material points and background grid nodes, and thus information can be mapped back and forth between them.Due to the different ways that may be adopted for solving the equation of motion in time in the second phase, the formulation can yield either implicit or explicit MPM approaches.

Continuum equations
At the continuum scale, the governing differential equations under purely mechanical loading can be derived from the respective conservation equations for mass and momentum, supplemented with a suitable constitutive equation to describe the stress-strain relation.In Eqs. ( 1) and ( 2), q is the mass density, v is the velocity, r is the Cauchy stress, and b is the body force due to, for example, gravity.
The mass of a given material point is independent of time, and hence Eq. ( 1) is automatically satisfied.For Eq. ( 2), a derivation based on the static equilibrium between the internal force, represented by r, and external force, represented by b, is introduced first for simplicity, i.e.
By applying the principle of virtual displacement, followed by the use of the divergence theorem, the equilibrium equation expressed in the weak form [18] with respect to the current configuration, at time t + Dt, is given by Z where S is the second Piola-Kirchhoff stress tensor, de is the Green-Lagrange strain tensor, du represents the virtual displacement, s denotes the prescribed part of the traction on the surface S and the volume of the body is represented by V.
Using the last known configuration at time t as a reference, the stress can be decomposed into the incremental form, whereas the strain at time t + Dt, with respect to the time t, is actually the incremental strain e tþDt ¼ De.The incremental strain is then divided into two parts; a linear part as commonly used in small strain analysis, plus a high order term, i.e.
where u is the incremental displacement.By substituting Eqs. ( 5) and ( 6) into the weak formulation (4) and assuming, for the moment, that the loading is deformation independent, then, by expressing the right hand side of Eq. (4) as ext , which is the external loading accounting for the effects of both body loads and tractions, and neglecting the high order term R V t Dr Á dDg, the small strain equation of motion in the Updated Lagrangian (UL) formulation is obtained as, To handle a large strain increment in implicit MPM, an appropriate stress measure should be used.Generally, the incremental stress can be calculated either from the second Piola-Kirchhoff stress and Green-Lagrange strain tensors (see [19]), or, as in this paper, via a rate dependant formulation using the Jaumann stress rate and velocity strain tensors [19].For consistency, the stress and strain rates are here shown in incremental form, with the Jaumann stress rate being given by, where Dr J is the Jaumann stress increment and C is the incremental stress-strain tensor.
Following Bathe [19], the Cauchy stress increment can then be written as where Dx is the spin tensor (also called vorticity tensor) increment at time t, in which Substituting Eqs. ( 8) and ( 9) into the equation of motion for small strain, Eq. ( 7), the final equilibrium equation for large deformation analysis is obtained as

Spatial discretisation
To solve Eq. ( 11) it must be spatially discretised.MPM discretises a continuum body in the original configuration into a finite set of N p material points that are tracked throughout the deformation process.The points are selected to represent a material sub-domain and do not have a defined shape.The method solves the equations on the background mesh; hence this spatial discretisation is undertaken utilising typical finite element methodology, with the major difference being that the integration uses the material points directly.
Taking the first term in Eq. (11) as an example, which is an incrementally linear stiffness term multiplied by the unknown displacement, and using the same method as in FEM, the shape functions and nodal values of displacement are used to approximate the continuum field, i.e.
u. Using the strain-displacement transformation matrix, De ¼ B L u, and the method of weighted residuals, the stiffness part of the term can be easily transformed into matrix form, i.e.
where B L is the linear strain-displacement transformation matrix.The integrals of the weak form are then converted to the sums of quantities evaluated at the material points, which yields where, for a single element, p is the number of material points within the element and the shape function differential B L is a function of the material point positions x p , which are moved after each time step.Similarly, the second and third terms in Eq. ( 11), referring to the non-linear contributions due to the geometry change and Jaumann stress, are written as (with the displacement part omitted), where the nonlinear strain-displacement transformation matrix is expressed as: where @N p =@x and @N p =@y are the shape function differentials with respect to the Cartesian coordinates, x and y, respectively, for the configuration at time t.r _ p is the Cauchy stress matrix, while the matrix rp , used in the second term of the integrand of K NL , is defined (for plane strain problems) as: Note that the component K t NL becomes highly significant when the stresses are of the same order as the material stiffness moduli.
The internal force, i.e. the last term in the Eq. ( 11), at the reference time t and the external force in Eq. ( 4) are respectively expressed as, where Nðx t p Þ are the shape functions at location x t p at time t.After substituting Eqs. ( 13), ( 14), ( 17) and ( 18) into Eq.( 11), the equilibrium equation in matrix form can then be formulated as: where , and is comprised of both linear and nonlinear terms.

Dynamic form
A dynamic governing equation is obtained by adding an inertial term in Eq. ( 19), satisfying Eq. ( 2).Hence where M t is the mass matrix at time t, which, in lumped form, is given by, where M p is the material point mass matrix, of size N dim Â N dim (where N dim is the number of dimensions).Following Newmark's time integration scheme [20], where Dt is the time step, a tþDt ; v tþDt and u tþDt are the accelerations, velocities and displacements at time t + Dt, and a and d are time stepping parameters which influence the integration accuracy and stability, and are chosen as d = 0.5 and a = 0.25 in the following analyses.
The incremental displacement is given by u ¼ u tþDt À u t .Hence, substituting Eq. ( 23) into Eq.( 20), and rearranging the resulting equation in the form of Eq. ( 19), leads to, Let K t represent a modified stiffness matrix, which takes the , and assume the new external force contains the kinetic terms from the previous time step, i.e.
Hence the rewritten governing equation takes the form, so that both the dynamic and quasi-static formulations can be solved in the same manner using the Newton-Raphson method.
To increase the numerical stability, it is common to include a damping force in the governing equation, which is assumed to be proportional to the nodal velocity.''Rayleigh" damping, by assuming the damping matrix to be a linear combination of the mass and stiffness matrices, has been widely proven to be stable and useful [21]; however, it is frequency-dependent, and therefore prior knowledge about the frequency of the system is needed.Cundall [22] described a local non-viscous damping to overcome the issue of frequency dependence, in which the damping force on a node is proportional to the magnitude of the out-of-balance force, with a directional function that ensures that vibrational modes are damped, i.e.
where c is a dimensionless damping factor, here chosen to be 0.75 [23], f i are the nodal resultant forces, and signðv i Þ is the velocity direction.Therefore, by adding the damping force into Eq.( 25), the final governing equation for the analysis can be formulated as,

Update of kinematic variables
In a quasi-static analysis, after the incremental displacement is obtained by solving the equation of motion, Eq. ( 19), the next step is only to move the material point in accordance with the displacement of the computational grid.However, dynamic analysis involves the solution of Eq. ( 27) to update accelerations and velocities.Details of kinematic variable updating are therefore presented in this sub-section.
Corresponding to the three MPM phases illustrated in Fig. 1, the kinematic variables need to be updated three times within the MPM computational cycle.

Mapping phase
Since information, such as velocity and acceleration, is initially stored at the material points and the background mesh is reset regularly after each time step, it is necessary to map (i.e.transfer) the associated kinematic information from the material points to the grid nodes at the start time t of each time step by utilising the shape functions, i.e.
in which i refers to a grid node, p refers to material points surrounding the grid node, m p is the material point mass, v p is the material point velocity and m i is the node mass which is assembled from the material points within the adjacent elements, see Eq. ( 21).
The acceleration is updated in the same way, with the displacement being initialised to zero at the beginning of each time step, i.e. u t i ¼ 0.

Updated-Lagrangian phase
For this computational stage, Newton's method is used to calculate the incremental displacement u i (Eq.( 27)).Solving for the velocity and acceleration vectors, using Eqs.( 22) and ( 23) with d ¼ 0:5 and a ¼ 0:25, the following is obtained: 2.5.

Convection phase
The final stage is to map the information from the grid nodes back to the material points.The acceleration and displacement vectors are directly mapped from the grid nodes using the shape functions, i.e.
where N n is the number of grid nodes which provide support to the material point and in this case is the number of nodes in an element (local support), while, in cases of utilising generalised interpolation MPM [24], the nodes of surrounding elements can also be included to provide a non-local support so that the cell-crossing error could be reduced.The trapezoidal rule is used here to update the velocity,

Explicit MPM
To compare the results and illustrate the performance of the implicit MPM presented in this paper, a version of the code utilising explicit time integration has also been produced, following the approach by Sulsky et al. [1,3], where the discrete form of the momentum equation (Eq.( 2)) becomes Once the solution for the nodal acceleration is obtained, other kinematic variables, i.e., velocities and displacements, are advanced explicitly.Furthermore, to improve the algorithmic stability of the explicit MPM, the authors follow Sulsky et al. [3], who proposed a ''modified update stress last" (MUSL) approach, in which the updated material point momentum vector is used to calculate the nodal velocities and thereby update the material point strain and stress rates.The explicit MPM version is herein applied only as a comparison, and therefore detailed aspects are not discussed.

Numerical examples
In this section, numerical examples for both quasi-static and dynamic analyses are provided.Firstly, a cantilever beam is considered as a quasi-static case to illustrate how the algorithm behaves.In particular, it is shown that, by adding an extra stiffness to the background mesh, known as a ''soft stiffness", the numerical stability of the solution is improved when elements exist that are not fully filled by material points.Moreover, the influence of the num-ber of material points per element on the numerical results is discussed.Secondly, the collision between disks has been investigated.This is a typical dynamic test case in MPM, and has been investigated by other researchers [17,25].Here it is considered for a direct comparison between the implicit and explicit MPM codes.Finally, a geotechnical application, using the dynamic MPM, is considered, which investigates the effect of an excavation on slope instability and post-failure behaviour.

Cantilever beam
A 2-D (plane strain) linear elastic cantilever beam has been analysed, and three aspects investigated: (a) the algorithm to initialise the positions of the material points; (b) the influence of the initial number of material points per element; (c) the influence of the extra ''soft stiffness" used across the background mesh.
The beam, as shown in Fig. 2, has a length of 1 m, a depth of 0.3 m, a Young's modulus of 100 kPa and a Poisson's ratio of 0. It is built-in along its left edge and its initial configuration is shown in grey, representing the originally filled background mesh.The self-weight was increased from 0 kN/m 3 to 4 kN/m 3 in 20 loading steps, with each applied increment being equal to 0.2 kN/m 3 .The analysis used 8-noded quadrilateral elements, with 16 material points being initially located in each element.
In order to increase the numerical stability, a ''soft stiffness", as described by Lim [26], has been assigned to the background mesh, which has been assembled using conventional Gauss point integration of the mesh (i.e.independently of the material points).The ratio of the soft stiffness to the actual stiffness is represented by g, where g = 0.01 is applied in this case.
The final configuration of the beam is also shown in Fig. 2, with the tip displacement (taken from the tracked material point located in the bottom right corner) being À0.38 m and À0.61 m, respectively, in the horizontal and vertical directions.Accordingly, the final activated mesh/elements are shown as a wireframe, which was determined by tracing the material point positions.The colouration of the material points represents the longitudinal stress along the beam, with blue representing compression and red representing tension, ranging from À40.8 kPa (compression) to 27.6 kPa (tension) on the material points located nearest the fixed end.By using extrapolation, the stresses at the element nodes can be obtained, ranging from À47.5 kPa to 32.0 kPa.UL FEM solutions, which are considered to be capable of providing accurate results for large displacement, but small strain, problems [14], are provided here as a direct comparison, where, for the utilised 10 Â 2 element mesh, the stress range of À47.4 kPa to 32.6 kPa at the fixed end is in good agreement with the MPM solution.
Due to the use of high order elements in the quasi-static code, the potential problems associated with using linear elements, i.e. locking, cell crossing [24], etc., are avoided.However, by omitting the inertial terms, the code stability decreases, which leads to the need for more material points and a ''soft stiffness" being used.After a brief description of the method for initialising the material points, the influences of the number of material points and ''soft stiffness" are investigated in the following.

Subdivision algorithm
There are many potential ways to initialise the locations of the material points, such as evenly spaced over the area and using Gauss point positions.Here, the subdivision algorithm of Lim [26] is chosen.For example, an 8-noded quadrilateral parent element with 4 Gauss points is adopted initially, and then, for obtaining more material points, the parent element is subdivided into smaller cells; e.g.into 2 by 2 cells, and then, by placing the material points on the local Gauss point positions of the small cells, to give a total of 16 material points for the original finite element.

Soft stiffness influence
The addition of an extra small stiffness across the background mesh is to increase the numerical stability of the code.For example, as a single material point moves near a finite element boundary, the shape function (at the material point) corresponding to the farthest element node may approach zero, while the gradient of the shape function will not be equal to zero; hence it is very hard to get the resultant force on the farthest node converging to zero, thereby leading to a non-convergent analysis.For implicit MPM, the problem normally occurs during the stiffness matrix formation, or factorisation, and can lead to an extremely large/unrealistic solution being output, as illustrated by the near vertical dotted line in Fig. 3.In the second loading step, the code stops due to an extremely large displacement at the beam tip.In contrast, even with a very ''soft"/small background stiffness, the outputs are realistic.Results obtained with four different stiffness ratios, i.e. ratio of the soft stiffness to the actual stiffness, of 0.1, 0.01, 0.001 and 0.0001, have been included here, in order to gain a picture of the influence of the ''soft stiffness".16 material points are initially placed in each element.
For comparisons with the UL-FEM solution, the vertical tip displacement of the beam (which is taken to be the average vertical displacement of the right-most column of material points near  the beam tip) is shown in Fig. 3.It is seen that, the higher the stiffness ratio, the greater the error.When the ratio is increased to 0.1, at the final weight of 4 kN=m 3 the error relative to the UL solution reaches 15%, although, for ratios equal to or smaller than 0.01, there is only a small difference in the results.Meanwhile, as the stiffness ratio increases, the code gets more stable.Therefore, for the beam case considered here, a ratio of 0.01 is adopted for all other simulations.

Influence of the number of material points per element
By using the subdivision algorithm introduced above, three different numbers of material points per element (in the initial mesh) have been considered, i.e. 4, 16 and 36.The influence of the number of material points has been investigated by comparing the averaged material point displacements across the beam tip with the UL-FEM results.Fig. 4 shows the relationship between the self-weight and the recorded displacement ((a) horizontal and (b) vertical).It is seen that, with initially 4 material points per element, there is a large divergence from the UL-FEM result, whereas, as the number increases, the differences become negligible, i.e. for 16 and 36 material points per element.
By putting more material points within an element, it increases the possibility of material points being distributed over the mesh more evenly, and thus it increases the code stability.However, as the number increases, the computation cost goes up.Hence, a balance has to be made, and here, initially 16 material points per 8noded element are considered appropriate to analyse the beam.
In terms of generalising the results in this section, it is difficult to determine a priori the appropriate level of discretisation in terms of both the number of elements and number of material points per element.As with the number of finite elements, good numerical modelling practice is recommended to ensure spatial convergence of the results, i.e. by increasing the level of discretisation and checking the consistency of the results.In the following two examples the presented results are for levels of discretisation that have been checked in this way.

Collision between two elastic disks
For analysing the collision between two elastic disks, the implicit dynamic method is utilised.Two things are worthy of note: (a) due to the presence of the mass matrix, in which negative terms may arise for high-order elements, low-order elements are adopted in the dynamic analyses (i.e. a 4-node quadrilateral element in this paper); (b) as the inertial term is included, initially 4 material points per element are found to be enough for ensuring accuracy when carrying out dynamic analyses, and an extra ''soft stiffness" is not necessary.
Two disks with radii equal to 0.2 m are simulated, assuming plane strain conditions.Initially, the two disks are located in the lower left and upper right corners of a square background mesh, as shown in Fig. 5(a).The mesh comprises 400 equal-sized square elements of side length equal to 0.05 m.Initial velocities (1.0 m/s, 1.0 m/s) and (À1.0 m/s, À1.0 m/s) are assigned to the two disks, so that they move towards each other along the diagonal of the square.The Young's modulus is 100.0 kPa, the Poisson's ratio is 0.3 and the density is 1.0 kg/m 3 .The simulation is run to a final time of t = 37.5 s.
Fig. 5 shows the collision process of the two disks, in which the colouration represents the mean stress distribution within the disks.In Fig. 5(a), the initial positions of the two disks are shown.Before impact, due to the constant velocities assigned to the material points, no velocity gradient will be experienced during this phase, and hence no stress magnitudes are shown on the disks.Fig. 5(b) shows the simulation during the impact, revealing that the material points near the contact region have higher compressive stresses, as would be expected.Finally, Fig. 5(c) shows the disks after impact, where the velocities are now in the opposite direction to those before the impact.Small residual stresses can be observed on the disks, due to the free vibration of the disks after the impact.
The energy conservation errors of the system were tracked, as a function of time, throughout the analysis.At time t + Dt, the total energy of the system can be expressed as where E kin is the kinetic energy and E strain is the strain energy stored in the material points.These are defined as Fig. 6 shows the system energy conservation errors with time for different applied time steps, where the initial velocities were changed to 0.005 m/s in order to simulate a longer time during which the two disks remain in contact.Fig. 6(a) shows the results using the explicit MPM, where the time step was chosen to be Dt = 0.0025 s by applying the Courant-Friedrichs-Lewy (CFL) condition (CFL is a necessary condition for stability with explicit time integration schemes).The results obtained with the implicit MPM are provided in Fig. 6(b)-(d), where time steps were chosen to be the same, 10 times and 100 times that of the explicit code, respectively.
For the cases shown in Fig. 6(a)-(c), the energy is well preserved before and after the impact.Before impact, the material points comprising the disks are moving with uniform velocities, and no strains or stresses are generated; hence the energy is purely due to the uniform movement of the disks.During the impact, the kinetic energy decreases to zero, with the strain energy reaching its maximum value at the point of maximum deformation.After that, the disks start to bounce back each other and, at the same time, the deformation is recovered, leading to a decrease in the strain energy and a resurgence of kinetic energy.After the impact, a small amount of strain energy remains, which is due to the free vibration of the disks after separation.
Fig. 6(d) shows the energy conservation errors corresponding to a time step size of Dt = 0.25 s, i.e. 100 times the explicit Dt.It is seen that the energy not conserved is severe, although the material point trajectories were not much affected.As pointed out by others, (e.g.[21]), implicit integration schemes such as the trapezoidal rule are unconditionally stable in linear analyses, which enables the code to work at extremely large time steps.However, for the problem presented herein, there is a bound on the implicit time step imposed by the accurate resolution of the collision, known as the characteristic collision time, i.e. the time that it takes for a wave to traverse the disk [17].In this example, the wave speed is about 10.0 m/s, and the disk diameter is approximately 0.4 m, thus giving the characteristic collision time to be about 0.04 s, i.e. 16 times the explicit time step size.It is worth noting that a reduction in the required number of time steps would typically result in a reduction in the computational time.However, as the solution scheme is different in the two approaches, that is, using an implicit solver (of which a number of types are available) in implicit integration schemes compared to explicit solution schemes where the new result can be simply calculated, it means that a unique direct relationship cannot be established.Moreover, as demonstrated in the above example, the specific characteristics of the analysis (in this case the wave speed and disk diameter) affect possible reductions in the number of time steps required.

Slope collapse under the influence of excavation at the slope toe
Excavation at the toe of a slope can trigger the collapse of the slope.As well as changing the slope geometry, it may also expose significant geologic features such as shear zones, faults and folds in some circumstances, thereby causing slope instability.
To investigate the effect of excavation on a slope, a naturally stable slope was first considered, as shown in Fig. 7, and then a vertical cut of height 0.2 m was excavated at the slope toe, as indicated.The problem is assumed to be plane strain, and the (a) Initial locations (b) At impact (c) After impact computational grid is made up of 4-noded square elements of 0.05 m side length, with each element initially containing 4 material points.The height of the slope is 1 m and the slope angle is 45°.Roller boundaries are prescribed at the left side of the domain to allow only vertical displacements.The interaction between the slope base and the ground is modelled using a Coulomb frictional boundary, which allows the slope base to move horizontally when friction is overcome.
A linear elastic and cohesion softening Von Mises model was used to describe the soil behaviour, as shown in Fig. 8.The initial yield stress s p is defined by a peak shear strength of c p = 5.0 kPa, whereas the residual strength s r is given by c r = 1.0 kPa.The ing modulus is H ¼ dc=d e p ¼ À50:0 kPa, in which the plastic strain invariant e p , in its incremental form, is defined as d e p ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi Since a local softening model is used here for demonstrating the proposed solution procedure, the numerical results are meshdependent.However, as the proposed method is based on a FEM formulation, the same regularisation techniques as used in FEM [27] can be used to address this issue.
The dynamic implicit MPM was utilised in this case.The computation was divided into two stages.Initially, the self-weight of c ¼ 20 kN/m 3 is applied to each material point in the slope in one step to generate the initial total stresses.Quasi-static equilibrium is detected by using the criterion based on the energy and force ratio [23], with the out of balance force ratio expressed as, and a dimensionless energy ratio as, where E kin denotes the kinetic energy of the system, W ext is the work induced by the external force, and a tolerance of 0.01 was used for both criteria.To obtain faster convergence in the quasi-static equilibrium stage, a local damping of c = 0.75 was used.Secondly, the triangular soil block was ''excavated" instantly to trigger the slope collapse.The collapse process is illustrated in Fig. 9, with the colouration representing the accumulated plastic strain invariant e p .
Due to the toe excavation, the slope became unstable and a slope slide was triggered.In Fig. 9(a), a complete shear band starting from the new slope toe was formed.Fig. 9(b) shows that the soil above the shear band moves as a block along the failure surface, and the final quasi-static configuration is displayed.
This analysis illustrates the ability of the proposed method to simulate geotechnical behaviour at large strains, which generally occur in arbitrary directions, i.e. not aligned with the mesh.The consequence of a geotechnical instability can then be better quantified, as can further potential unstable situations be observed, for example, in retrogressive and progressive failures, which will be investigated in future publications.

Conclusions
An implicit material point method (IMPM) framework has been developed for both quasi-static and dynamic analyses.The improved characteristics in solving large-deformation elasto-   plastic problems as compared to the finite element method (FEM) are highlighted.The paper provides a detailed description of the IMPM formulation and numerical implementation, following, where possible, standard FEM methods, thereby enabling the adaptation of an FEM code to an MPM code in a relatively straightforward manner.
Three numerical examples of increasing complexity are provided so as to clarify and illustrate the implementation procedures, performance and behaviour in a step by step manner, including features particular to MPM, e.g. the number of material points per element and a background ''soft stiffness".The cantilever beam problem represents a quasi-static case, and is used to demonstrate various details of the MPM, including spatial discretisation, the influence of using a ''soft stiffness" for the background mesh, and the influence of the number of material points per element.The dynamic case is illustrated via the disk collision example, where comparisons with an explicit MPM are provided to identify the advantage of using IMPM over explicit MPM in terms of the chosen time step size.Finally, an analysis of slope collapse under the influence of toe excavation shows that IMPM is a useful tool to capture post-failure behaviour in geotechnical applications.For large deformation analyses, however, cell-crossing errors are an important issue to be addressed.Recent efforts have been made by other research teams, to reduce the cell-crossing errors and to better combine MPM with FEM for modern engineering applications [28][29][30].Future work will be conducted to improve IMPM for general applications.

Fig. 3 .
Fig. 3. Applied load versus tip vertical displacement for different ratios of soft stiffness.

q
kde p k, and de p is the plastic strain increment, normal to the plastic potential surface.The Young's modulus and Possion's ratio of the soil are E = 200 kPa and t = 0.33, respectively.
(a) Shear band formed due to toe excavation (b) Final configuration of the collapsed slope

Fig. 9 .
Fig. 9. Evolution of slope failure due to excavation at the toe.