On the implementation of a material point‐based arc‐length method

The material point method is a versatile technique which can be used to solve various types of solid mechanics problems, especially those involving large deformations. However, the capability of the material point method to track a load‐displacement response can deteriorate once a limit point, such as snap‐through or snap‐back, in the response is encountered. One way of overcoming this is to use path following techniques, such as an arc‐length method. This technique is well established in finite element analysis but not within any material point method formulation. This paper provides for the first time an arc‐length controlled implicit, quasi‐static material point method. The modifications to the standard arc‐length scheme to allow for the stable execution of an arc‐length solver within the material point method are detailed. The capability of the material point‐based arc‐length method is demonstrated through a number of problems, which include linear elastic, non‐linear elastic, linear elastic‐perfectly plastic and linear elastic‐plastic softening material behaviour under large deformations. The techniques presented in this paper are essential for arc‐length techniques to be applied effectively to the material point method and the combination of these techniques makes the method suitable for new problems that cannot be solved with existing implicit material point approaches.


INTRODUCTION
In the finite element method (FEM), for a quasi-static, implicit solution of a non-linear problem, the applied load (or applied displacement) is incremented at the start of each load step and a Newton-Raphson scheme is often used to determine the corresponding displacement (or load) such that the out-of-balance forces of the system are below a given tolerance.However, this procedure only works if the incremented variable is strictly increasing throughout the simulation such that the equilibrium path can be traced successfully.Figure 1A shows an increasing equilibrium path in which a single value of force corresponds to a single value of displacement.It is possible to continue increasing either variable and determine a solution on the equilibrium path.On the other hand, problems involving large deformations 1 and some types of non-linear material behaviour 2 will often have critical points along the equilibrium path where, after the critical point, it is not possible to increment the variable further and obtain a solution.A snap-through response is when an increase in displacement for a decrease in load is seen, as shown in Figure 1B.In the case of a snap-through response, a single force value can be associated with multiple displacement states and it is may not be possible to determine the correct state, meaning that it is not possible to model with a load controlled scheme.It is still possible to model a snap-through response with a displacement controlled scheme as the equilibrium path is strictly increasing in terms of displacement, but this can also change the physical problem being modelled depending on the nature of the imposed boundary conditions.A snap-back response is when a decrease in displacement is seen after a snap-through response, as shown in Figure 1C.As with a snap-through response, it is not possible to use a load controlled scheme due to the multiple possible displacement states for a given force, a snap-back response also has multiple load states for a given displacement.As the equilibrium path is not strictly increasing in terms of displacement for a snap-back response, it is not possible to model this response using a displacement controlled scheme.][5] The arc-length method was independently developed by Riks 3,4 and Wempner 5 in the 1970s and is one of the main techniques used in finite element analysis (FEA) to overcome the issues associated with the solving of snap-through and snap-back equilibrium paths. 6Crisfield 7 adapted the method proposed by Riks such that the arc-length method is more suited for use in the FEM when used with the modified Newton-Raphson method.Since its conception, many variations of the arc-length method have been proposed to solve complex equilibrium paths including (but not limited to) Bergan et al., 8 Ramm, 9 Powell and Simons, 10 Lam and Morley, 11 Fafard and Massicotte, 12 Krenk, 13 Nicholson, 14 May 15 and Pretti et al. 16 The essence of the arc-length method comes from the constraint equation, or arc-length equation, which is included in the non-linear equations of the problem and solved to give an incremental load factor and corresponding incremental displacements.By solving for both load and displacement throughout the simulation, the issues which may arise in traditional load/displacement controlled schemes due to snap-through or snap-back responses are mitigated.The arc-length method has been used to solve a variety of non-linear problems including buckling of cylindrical shells, [17][18][19][20] structural analysis, 21 modelling of crack propagation 22,23 and failure of brittle materials. 24,25hile the arc-length method is well established in the FEM, when modelling problems involving large deformation mechanics the FEM can run into issues involving volumetric locking and large mesh distortion 26,27 leading to inaccurate results and may necessitate the expensive process of re-meshing.The material point method (MPM) is an alternative to standard Lagrangian schemes and was developed by Sulsky and co-workers 28,29 as a solid mechanics extension of the FLuid Implicit Particle (FLIP) method 30 which is itself an extension of the Particle-In-Cell (PIC) method. 31As a mixture of Eulerian and Lagrangian schemes, the physical domain is discretised by Material Points (MPs) which are immersed in a background mesh of finite elements.The MPs act as integration points and deform through the background mesh which remains unchanged throughout the analysis, avoiding the large mesh distortion seen in the FEM.The MPM is a tool that has been widely used to model many engineering problems including manufacturing processes, 32 large deformation problems in geomechanics, [33][34][35] fracture mechanics 36 and contact. 37This paper will focus on an implicit quasi-static implementation of the MPM, adapting the AMPLE (A Material Point Learning Environment) framework 38 in order to incorporate the arc-length solver.
To the authors' knowledge, there is no current literature detailing the implementation of the arc-length method in a MPM framework and this paper describes the first Material Point-Based Arc-Length Method (MP-BALM).The combination of an arc-length method with the MPM expands the applicability of the MPM to load-controlled problems involving snap-through, which are currently not possible with techniques in the published literature.Section 2 covers the basic principles of the arc-length method as it is implemented in the FEM, followed by an explanation of the MPM and the adaptations needed for the arc-length method to be used successfully with the MPM in Section 3. Section 4 then presents numerical examples using the MP-BALM, including problems involving linear elastic, non-linear elastic, linear elastic-perfectly plastic and linear elastic-plastic softening material behaviour under large deformations.While the arc-length method and the concepts presented in this paper are able to analyse problems involving both snap-through and snap-back responses, the MPM is a continuum method for which, to the authors' knowledge, there are no known problems involving global snap-back responses.Therefore, only problems involving snap-through responses are analysed in Section 4.

THE ARC-LENGTH METHOD
In the standard FEM, an arc-length solver (or Riks solver) solves for both load and displacement based on the converged state of the previous load step and a specified arc length.The equilibrium path is formed of load factor-displacement pairs, {, u}, where  is the load factor and u are the nodal displacements.The load factor acts as a scalar multiplier of the external force vector applied to the material domain such that where f ext 0 is the initial applied force vector and f ext is the force vector corresponding to the current load factor.A key feature of the arc-length method is that the distance from the current converged state to the next state along the equilibrium path is limited by a specified arc length, Δl.The most general form of the arc-length constraint is the spherical arc-length constraint 2 which is given by where Δu and Δ are the converged values of the change in nodal displacements and the change in load factor respectively for the current load step, and  is a prescribed scaling parameter.This paper utilises the more widely used version of the arc-length constraint, the cylindrical arc-length method, which sets the scaling parameter  to be 0, thus only ensuring that the Euclidean norm of the change in nodal displacements is equal to the prescribed arc length. 2 The cylindrical arc-length constraint is given by Figure 2 shows the use of the cylindrical arc-length method to track through the equilibrium path for a simple two degree of freedom system.Given a converged state, ( i , u i ), when visualised in terms of two displacement components, u 1 and u 2 , the constraint equation in (3) creates a cylindrical zone around the converged point.The equilibrium path will intersect the cylinder at two points, A and B. These are the two points on the equilibrium path which satisfy the arc-length constraint based on the change in nodal displacements from the converged state.Based on the direction of travel along the equilibrium path taken in the previous load step, either point A or point B is chosen to be the next converged state and the process is repeated for the next load step.The arc-length scheme is often paired with a Newton-Raphson solution scheme such that equilibrium is achieved as well as satisfying the specified arc-length constraint.In each iteration of a load step, the load factor increment, , is found by solving the quadratic equation 7 The cylindrical arc-length method (reproduced from de Souza Neto et al. 2 ).
with the coefficients of the quadratic given by where Δu (k−1) is the change in nodal displacements calculated so far in the load step.The tangential displacement, u, and the Newton-Raphson displacement * , u * , are given by where K is the global stiffness matrix and f oobf is the out-of-balance force vector.As the form of ( 4) is a quadratic equation, two possible solutions for  will be found; one root will progress along the equilibrium path while the other will track back along the path.Therefore, the selection of the appropriate root is required to ensure that the the simulation progresses in the desired direction.One option to find the most appropriate load factor increment is the "angle criteria" which aims to minimise the angle between the change in nodal displacements so far in the load step and the newly calculated change in nodal displacements.This is done by maximising the dot product of the two vectors in question such that the value of  is chosen such that A second option for root selection is presented in Hellweg and Crisfield 39 which was proposed to counter the fact that the angle criterion encounters issues for problems involving sharp snap-back responses.This approach involves calculating the two roots of  (k) as normal, and the total displacements, stresses and internal forces are computed for each root and the root which gives the smallest residual is taken.As stated by Hellweg and Crisfield, this method need only be used in the presence of sharp snap-backs, however, this paper only focuses on snap-through responses thus no further consideration of this method is taken.
The incremental change in nodal displacements, u, is then calculated as a combination of the Newton-Raphson displacement and the tangential displacement (which is scaled by the load factor increment determined by ( 4) and ( 7)) such that At the start of the simulation, an initial arc length is required.It is possible to simply prescribe this value, however, this is not easy for the user to determine a reasonable initial arc length as it is dependent on the mesh size.Instead, the initial arc length can be calculated as the Euclidean norm of the iterative increment of the tangential displacement at the start of the problem, where where the u used is the initial tangential displacement based on the reference load where an initial iterative load factor increment is prescribed as 1 to produce an initial arc length value.A choice can be made to set this initial arc length value from (9) to be the maximum allowed arc length throughout a simulation, however, this choice is not a condition to be satisfied in the arc-length method.In principle, the arc length can take any value but there may be cases where the arc length becomes too large and causes issues in the analysis, so imposing an upper limit for this value is often the best action.If the number of Newton-Raphson iterations needed in a load step exceeds a prescribed desired number, the arc length for the following load step is adjusted accordingly, allowing for an increase or decrease in the distance away from the current converged state to the next converged state.The arc length for the following load step is calculated at the end of each load step such that where k des and k n are the desired and actual number of Newton-Raphson iterations needed in the load step respectively and  and  are positive scalars.In the literature,  is taken to be 1 7 or 0.5, 11 however, Ramm 9 explains that when  = 1, oscillations in the number of Newton-Raphson iterations needed in each load step are seen and therefore taking  = 0.5 is a more stable approach.The value of  controls the minimum allowed arc length through the entire simulation such that Δl i is a value in the interval , this ensures that if the required Newton-Raphson process continues to exceed the desired number of iterations, the arc length will not continue to decrease so much that the tracing of the equilibrium path begins to halt as the allowed distance becomes tiny.There is no theory on the choice of value of , therefore  = 4 is chosen here based on Pretti et al. 16 In the first iteration of each load step, there is clearly no information as to the change in nodal displacements, Δu (0) = 0.In this case, solving (4) produces the solutions for the initial load factor increment of giving a magnitude of the initial increment.However, it is not possible to use (7) to determine which solution to choose and a predictor is required.It is important to choose the correct value of the predictor solution as one will progress the simulation along the equilibrium path and the other will track back along the path. 1 Multiple options of predictor solution criteria have been presented in the literature to choose the correct sign for  (1) to progress along the path, including: 1. Stiffness determinant 40 -sign determined by the determinant of the current stiffness matrix: 2. Incremental work 10 -sign determined by predictor work increment: 3. Secant path 41,42 -sign determined by deformation of the previous load step (Δu n ): Option 1 is widely used and is easily implemented as the stiffness matrix will be known, however, issues may arise if the stiffness matrix is close to singular 7 or in the presence of bifurcations where the solution will oscillate around the bifurcation point. 2 Option 2 is able to trace an equilibrium path through a bifurcation point but is not able to effectively trace snap-back responses. 2Both options 1 and 2 use information based on the current equilibrium state and nothing about the history of the equilibrium path.Option 3 involves the use of Δu n which takes the deformation of the mesh in the previous load step into account when determining the deformation of the current load step.The secant path predictor is able to overcome the issues seen with the stiffness determinant and incremental work predictor schemes.An alternative predictor solution was proposed by Kadapa et al. 20 where the nodal displacements and load factor for the first iteration of a load step are calculated as a linear combination of the respective nodal displacements and load factor of the previous two load steps such that change in nodal displacements and load factor increment in the predictor step are given by u (1) = Δu n  (1) This is a simpler and less expensive predictor solution which does not require any additional computation to determine the direction of the progression along the equilibrium path in the predictor step.
Figure 3 shows the process of the arc-length solver through the iterations of a single load step.From the previously converged point, ( i , u i ), the predictor iteration determines the initial load factor increment (with direction determined by the secant predictor) which is used to calculate an initial change in nodal displacement in order to satisfy the arc-length constraint.In subsequent iterations, the load factor increment is calculated using (4) and then the change in nodal displacements are calculated using (8), these steps are repeated until the Newton-Raphson error is below an allowed tolerance to give the next converged point, ( i+1 , u i+1 ).
F I G U R E 3 Arc-length solution scheme using Newton-Raphson iterations.
The arc-length method detailed in this section is presented in the form of a FEM-based method to solve problems with the path following technique.However, there are some key differences that must be considered when implementing the arc-length method into the MPM formulation, and these will be explained in the next section.

THE MATERIAL POINT-BASED ARC-LENGTH METHOD
The key steps of the standard MPM algorithm for solid mechanics are depicted in Figure 4, starting from the upper left image and moving left to right, top to bottom, the procedure for a load step is: 38 (a) the physical domain is discretised by MPs which carry data such as mass, stress, deformation, and so forth, of the material and the positions of these points within the background mesh are determined; (b) the MP data for stress, volume and stiffness is mapped onto the element nodes associated with that MP using the basis functions to generate the internal force and stiffness quantities; (c) the governing equations are assembled on the background mesh nodes; (d) the governing equations are solved for the unknown nodal displacements; (e) the mesh information is then mapped back to the MPs to update the deformation, stress, volume, and so forth, and the MP displacements are calculated; (f) the MP positions are updated and the background mesh is reset.
For more information on the formulation of the MPM, the reader is directed towards the original work of Sulsky and co-workers 28,29 as well as the many papers 34,36,43,44 and textbooks [45][46][47] that cover the MPM.
The MPM approach used in this paper is based on the open-source AMPLE code, 38 a quasi-static implicit finite deformation MPM based on the updated Lagrangian formulation with the weak form of equilibrium given by assuming a field of admissible virtual displacements,  i , where Ω is the domain of the body with boundary Ω.The domain is subjected to a traction, t i , on the boundary of the domain (with surface, s) and body forces, b i , applied to the volume, v, of the domain.These tractions and body forces result in a Cauchy stress field,  ij , through the body.Discretising (16)  with background mesh elements, E, and imposing Galerkin's weak form produces where A is the standard assembly operator (in this case acting over all elements in the mesh), S vp is a matrix containing the basis functions to provide the mapping between the MPs and the background mesh, ∇ x S vp is the strain-displacement matrix containing the derivatives of the basis functions with respect to the global coordinates.In this paper, the boundary of the MP domain is generated using a B-spline representation 48 which allows Neumann boundary conditions (such as the tractions in ( 17)) to be applied directly to the MP domain.The shape functions used to integrate the traction forces over the domain boundary, S v , are the standard Finite Element shape functions.The weak form in ( 17) can be discretised by the MPs such that where the standard assembly operator is applied over each MP, p, and V p is the MP volume.
There are multiple options for the choice of basis functions used to map data between the MPs and the mesh nodes.The original approach was to use standard FEM shape functions 28 in which the MPs were assumed to be concentrated point volumes and the shape functions of the finite element mesh are used.An alternative to standard FEM shape functions are Generalised Interpolation MPM (GIMPM) basis functions 43 which are generated as the convolution of the finite element shape functions with a point characteristic function (taken to be a hat unity function).This generates a set of C 1 continuous basis functions using a domain centred on the MP position over which the volume of the MP is evenly distributed.Another option is to use higher order B-splines 49 instead of the piecewise-linear basic functions, which is seen to significantly reduce spatial quadrature errors. 50The higher order B-spline basis functions are generated based on a knot vector, with values associated with positions of the mesh nodes, and the location of the point in question relative to the knots within the knot vector.The value of the basis function is calculated recursively based on the order of the B-spline.Using a cubic B-spline produces C 2 continuous basis functions centred at the MP location and spanning four mesh elements in each dimension.This paper however, utilises the GIMPM basis functions throughout.
It is widely documented that the MPM can be unstable, mostly due to cell crossing errors, a problem mitigated by the use of GIMPM basis functions, however there are other issues linked to boundary representation.A key feature of the MPM is that the physical domain of the problem is made up of MPs, therefore, there is no physical boundary representation.For illustrative purposes, the "boundary" of the physical domain is depicted with a thick black line in Figure 4, however, this boundary is not explicitly represented in the method and the edge of the physical domain is implied by the extent of the MPs and their associated domains.If desired, the boundary of the MP domain can be represented using B-splines 48 where points with zero volume are placed along the boundary of the domain which are used as sample points for the creation of a B-spline curve in 2D, or surface in 3D.These boundary points positions are updated at the end of each load step in the same fashion as the MPs, allowing for the tracking of the boundary throughout the analysis.This B-spline boundary also allows Dirichlet and Neumann boundary conditions to be applied directly to the MP domain by integrating over the B-spline curve/surface.The reader should refer to Bing et al. 48for more information on B-spline boundary representation in the MPM.However, irrespective if the boundary is explicitly represented or not, the fact that the extent of the physical domain will not necessarily coincide with the background mesh means that there is the potential for background elements to only overlap with a small portion of the physical domain.This "small cut" problem can lead to poor conditioning of the global system of equations and impact on the stability of the MPM.In more detail, the ill-conditioning of the system of equations due to the presence of partially filled elements can be quantified by the system's condition number (the ratio of the largest and smallest eigenvalues of the global system matrix).The smallest eigenvalue of the system is related to the smallest intersection between the physical domain and the background mesh, as this intersection can be much smaller than the size of the elements in the background mesh, the condition number is unbounded, which is detrimental to the solution stage of the analysis.The ghost stabilisation technique proposed by Burman 51 penalises jumps in the gradient of the solution field across elements cut by the physical domain. 52This is especially useful for the MPM due to the fact that it is often the case that MPs around the edge of the physical domain are located in partially filled elements.Coombs 53 implemented the ghost stabilisation method within the MPM where boundary and faces are determined over which the ghost penalty was integrated, resulting in a well-conditioned system and improved stability in analysis.
When modelling large deformations of nearly incompressible solids, the issue of volumetric locking can occur due to low order elements being used and the low order interpolation polynomials inadequacy to represent the volume preserving displacement fields. 54One common method of avoiding this volumetric locking is the F-bar approach developed by de Souza Neto and co-workers 54,55 where the deformation gradient is decomposed into deviatoric and volumetric components such that where F, F d and F v are the deformation gradient, its deviatoric and volumetric components respectively with the two components being defined as A deformation gradient, F 0 , is computed which results from the displacement at centroid of the element.The modified deformation gradient, F is calculated as the composition of the deviatoric component of F and the volumetric component of F 0 such that The modified deformation gradient is then used to compute the stiffness matrix and internal force vector.Coombs et al. 26 implemented the F-bar approach in the MPM, which was chosen as it: does not introduce any additional unknowns; is straightforward to implement; does not restrict the choice of constitutive model; does not require any tuning parameters; and does not require additional MPs to track volumetric behaviour.Numerical examples presented in this paper utilise the B-spline boundary representation, ghost stabilisation and the F-bar formulation.

Reformulation of the arc-length method into the MPM
As the MPM allows the MPs to move through the background mesh, it is possible that the MPs will populate small portions of an element around the boundary of the physical domain.This can be seen in Figure 4 where the domain represented by the MPs has a boundary which cuts through elements of the background mesh.This causes issues with elements that participate in the analysis having very small stiffnesses, resulting in extremely large nodal displacements relative to the applied external load.This means that if an arc-length solver is implemented into the MPM in its FEM form, the large nodal displacements will account for the majority of the allowed nodal displacements of the arc-length constraint, manifesting in very small increments of the load factor and an inefficient tracing of the equilibrium path.We reformulate the arc-length constraint such that it is based on the change in MP displacements to produce the Material Point-Based Arc-Length Method (MP-BALM).The change in MP displacement, Δu p , is calculated by mapping the change in nodal displacements to the location of the MP such that where n is the number of nodes in the background mesh, S vp (x p ) is the basis function of node (or vertex) v based on the position of the MP, x p , Δu v is the change in displacement of node v and N p is a 3 × 3n matrix containing the basis functions for all nodes in the mesh.A diagonal matrix, N, of size 3n × 3n can be formed which contains the sum of all MP basis functions for every node in the mesh such that The values within N are scaled by the maximum diagonal value such that the matrix contains values between 0 and 1.The values in N are representative of the influence of the change in nodal displacements on the physical body, where nodes surrounded by fully populated elements have a value equal or close to 1 and the nodes associated with less physical material have a value of less than 1.Due to the nature of the MPM, the displacements of the MPs have more of a physical meaning compared to the displacements of the nodes of the mesh as it is the MPs that represent the deformation of the physical domain in a MPM analysis.Therefore, in a similar fashion to the MP displacement calculation in (22), a vector of scaled nodal displacements, Δu s , is calculated as therefore, the arc-length constraint in (3) can then be modified using scaled nodal displacements giving where Ñ = N T N is a diagonal matrix with values between 0 and 1.The result of the modification in (25) is that the nodes associated with partially filled elements will not take up as much of the allowed change in nodal displacements, improving the performance of the arc-length solver in tracing the equilibrium path with more uniform increments.This will be demonstrated with numerical examples in Section 4. The values in N are scaled in order to ensure that changing the number of MPs used in the analysis does not have a significant effect on the progression of tracing along the equilibrium path.For example, if no scaling is present, increasing the number of MPs will increase the values on the diagonal of N thus increasing the left hand side of (25) and resulting in smaller displacements per load step and a slower progression through the analysis.Changing the number of MPs should not result in significant changes in how an analysis progresses.The N matrix scales the nodal displacements in the arc-length condition based on a node's influence on the physical domain, however, as the values in the matrix are rescaled such that they are between 0 and 1, N has no physical meaning in terms of the problem.The reason for the use of N is to dampen the large nodal displacements of the partially filled elements around the boundary of the domain (as the corresponding values in N will be close to zero) while maintaining the displacements of the nodes which have more influence on the physical domain.
The addition of the Ñ matrix in ( 25) is then incorporated throughout the MP-BALM formulation to calculate the load factor increment where the new coefficients of (4) become The initial arc-length calculation in (9) also includes Ñ such that as well as for the predictor solution where the initial load factor increment is calculated as It is still important to include the scaling matrix in (28), as u is calculated from the assembled stiffness matrix for the MPs in the current state meaning that the partially filled elements at the boundary of the MP domain may still see large nodal displacements.The large nodal displacements may also be seen within Δu n as these are simply the calculated nodal displacements from the previous load step where the large nodal displacements may have also occurred.

Nodal displacement reconstruction
This paper utilises the secant path predictor in (28) (scaled form of ( 14)) for the determination of the direction of the equilibrium path in the first iteration of a load step where the converged nodal displacement of the previous load step is required.However, since a key feature of the MPM is that the background mesh is reset at the end of every load step, obtaining these nodal displacements is not trivial and instead using the nodal displacement vector at the end of the previous load step is not feasible as the body has deformed and the MPs have changed position, meaning that the domain is not influenced by the mesh nodes in the same way as in the previous load step.A method of reconstructing the nodal displacement field based on the MP displacements of the previous load step and the updated MP positions is therefore needed.This is done using a least squares weighted residual technique 48 where Δu a p are the known MP displacements from the previous load step and Δu n are the unknown nodal displacements which are to be solved for.Rearranging for the unknown nodal displacements we obtain It is worth noting that Δu n is evaluated globally in order to minimise the residual error across the whole domain meaning that the integrals in Equation ( 30) must be assembled into a global matrix and vector respectively.This means that Δu n is calculated as As the nodal displacements of the previous load step are reconstructed based on the MPs' position in the current state, there may still be cases where large nodal displacements are calculated for partially filled elements around the boundary of the MP domain, therefore, the scaling matrix is needed in the predictor meaning (28) should be used over (14)  determine  (1) using predictor solution -( 28) else calculate Newton-Raphson displacement, u * -(6b) determine  (k) by solving quadratic equation -( 4) with ( 26) end if calculate incremental change in nodal displacement, u (k) -( 8) else zero incremental displacement in zeroth iteration, u (0) = 0 end if end if increment change in nodal displacements, Δu (k) = Δu (k−1) + u (k) assemble stiffness matrix, K, and internal force vector, f int calculate out-of-balance forces, f oobf calculate Newton-Raphson error end while update MP positions update Δl -(10) end for to the AMPLE framework, the MP-BALM algorithm is reasonably similar.The nodal displacement reconstruction takes place at the start of each load step before entering the Newton-Raphson iteration loop.The arc-length solver itself requires multiple if statements in order to determine how to find the load factor increment and the arc length needs to be updated at the end of each load step based on the iterative performance of said load step.This algorithm is used to solve the numerical examples presented in the next section.

NUMERICAL EXAMPLES
A series of numerical examples are presented to demonstrate the capability of the proposed arc-length solver with a MPM formulation.The first example demonstrates that the solver can be implemented into a MPM code and produce correct results for simple problems with analytical solutions.Further examples are given to demonstrate the capabilities of an arc-length solver within the MPM with problems which encounter snap-through responses.

Compression of an elastic column under self weight
To demonstrate that the implemented solver is functional within a MPM code, a simple benchmark test is utilised, namely the 1D compression of an elastic column with an initial height of l 0 = 50 m under self weight.Although this problem is in essence 1D, a 2D plane strain implementation is used.The material is modelled as a Henky material with a Young's modulus of 10 kPa and a Poisson's ratio of zero.The background mesh is generated using square elements with roller boundary conditions applied over the base and sides of the mesh.The column is discretised using a distribution of 2 × 2 evenly distributed generalised interpolation material points in each populated element of the background mesh.A target body load of 800 N m −2 (g = 10 m s −2 and a density of  0 = 80 kg m −3 ) is applied.As the load at any point in the analysis is determined by the initial applied load and the load factor (as per (1)), an initial gravitational acceleration is prescribed and the problem is run until the load at the end of a load step exceeds the target body load.Selecting a smaller initial load means that more load steps are required to reach the target load.
To study convergence, the number of elements along the height of the column (y-direction) is varied between 2 2 and 2 13 (between 4 and 8192) while a single element is used through the width of the column (x-direction) throughout.The initial MP to mesh element ratio is also maintained throughout.The problem is run with the arc-length solver for three different initial body loads, 8, 4 and 2 N m −2 which results in the simulations running over 17, 33 and 67 load steps respectively.
The analytical solution for the normal stress in the y-direction at the end of the analysis is where Y is the original position of the point in the body, l 0 is the original height of the column and g n is the gravitational acceleration due to the load factor at the end of the simulation.Figure 5 shows the numerical results for the column under self weight using a mesh of 2 6 elements and an initial body force of 2 N m −2 , the calculated vertical stress is plotted against the original vertical position of each MP.The solid line in Figure 5 is the the analytical solution in (32).The error between the calculated stress in the column and the analytical solution is calculated as where  p yy is the normal stress in the y-direction at each MP position and V 0 = ∑ V 0 p is the initial volume of the column.As the arc-length method will not perfectly achieve the target body load and the simulation is simply run until the body load first exceeds the target load, each run will have a different value for g n .Running over 17, 33 and 67 load steps gives final g n values of 12.9, 11.7 and 10.2 m s −2 respectively.Figure 6 shows the convergence of the solution as the mesh is refined.It is clear that the rate of convergence matches the rates of convergence when running the same problem with a linear solver.It is also clear to see that the error value of the solver compared to the linear solver depends on total number of load steps F I G U R E 5 Numerical solution of a column under self weight plotted against the analytical solution.

F I G U R E 6 Compression of column under self weight convergence.
required to reach the target body force, this is due to the fact that the linear solver is run over 40 evenly spaced load steps whereas in the arc-length solver, the load increment is not evenly distributed throughout the simulation.As the number of load steps used to run the simulation is increased, the convergence behaviour of the MP-BALM tends towards that of the linear solver.

Hyperelastic cylinder under internal pressure
The next example is that of a 2D incompressible, hyperelastic cylinder under internal pressure.The cylinder has an inner radius R i = 1 m with the outer radius, R o , set such that the ratio = 0.85.The cylinder is modelled as a hyperelastic material using the regularised, isotropic Ogden model 2,56 with the strain-energy function, W, given in terms of the principal stretches,  i (i = 1, 2, 3) † , such that where  r and  r are pairs of material parameters, K is the material bulk modulus, J is the volume ratio.It is worth noting that due to the fact the cylinder is incompressible, it is assumed that J ≈ 1, however, this is not strictly enforced, therefore, incompressibility is ensured through a suitably large value for the bulk modulus.The material parameters are defined as 57 where and  is the ground-state shear modulus.The ground-state shear modulus is chosen arbitrarily as 0.375 MPa, which is the equivalent shear modulus based on a Young's modulus of 1 MPa and a Poisson's ratio of 0.3.The bulk modulus is set such that it is sufficiently large to enforce near incompressibility and takes the value of 99.8 MPa.
Due to symmetry, one quarter of the cylinder cross section is modelled with a plane strain assumption.A diagram of the problem is shown in Figure 7.The background mesh is generated with an overall width and height of 3.5 m in order to provide enough space for the domain to deform up to a final inner stretch of 3, in other words, until the deformed inner radius is three times the initial inner radius.Each element in the background mesh is a square with size h = 0.05m.In order to generate the initial domain of the cylinder, the entire mesh is filled with MPs where each element contains 8 × 8 evenly spaced generalised interpolation material points.It is important to ensure that there is a large number of MPs initially populating the mesh for optimal integration over the domain.The MPs which do not lie within the boundaries of the cylinder are then removed, leaving the desired physical domain made up of a total of 7722 MPs.This check is performed at the centre of the MP domain, determining whether the centre of the domain is within the boundaries of the cylinder meaning that some portions of the GIMP domains may exist outside of the cylinder boundary.Roller boundary conditions are applied over each edge of the background mesh.The problem is run using the ghost stabilisation technique 53 in order to reduce errors and improve the stability of the simulation due to partially filled background elements.The ghost stabilisation parameter is set to  = 10 6 .An initial normalised pressure, P * = P  , of 5 × 10 −3 is applied over the inner boundary of the cylinder using a B-spline boundary representation. 48No pressure is applied over the outer radius of the cylinder.An analytical solution for this problem is given as 57 where Ŵ(,  z ) = W(,  z , ( z ) −1 , 1) and Ŵ is the derivative of Ŵ with respect to  and the integration is performed between the radial stretch at the inner radius,  a , and the radial stretch at the outer radius,  b .As the problem in question is run using the plane strain assumption, there is zero stretch in the axial direction meaning  z = 1 throughout the analysis.Figure 7 shows the calculated inner stretch against applied internal pressure (normalised with the ground-state shear modulus) for the analysis using the new (scaled) approach.There is a clear snap-through behaviour displayed with the maximum in pressure seen at an inner stretch of 1.9.Comparing the calculated load-displacement curve with the analytical solution calculated using (35), the MP-BALM matches the analytical solution very closely throughout the simulation.At large displacements, small fluctuations can be seen in the calculated load-displacement curve, this may be due to the fact that the physical domain does not have a defined boundary in the MPM with the use of the B-spline boundary being used as an approximation for the boundary.This means that the pressure will not be applied perfectly to the true boundary of the domain.Also, at large deformations, the incompressible cylinder is stretched by a large distance and the thickness of the cylinder is reduced, potentially meaning that a small number background elements are active through the thickness of the cylinder thus the integration over the physical domain is sub-optimal, especially over the boundaries.Figure 8 shows the undeformed and deformed configurations of the cylinder where it can be seen that at least two background elements span the thickness of the cylinder in the undeformed configuration whereas only a single element contains the full thickness of the cylinder in the deformed configuration for this problem setup.
The problem is run using the traditional (unscaled) and new arc-length approaches for the same problem setup.The traditional approach takes 299 load steps to reach the desired final inner stretch, whilst the new approach takes only 233 load steps in order to complete the analysis.As mentioned above, the results from the new approach fit the analytical solution very well, a similar set of results are also seen with the traditional approach, however, more erratic fluctuations in the curve are seen for the highly deformed state ( a > 2.6).In order to compare the two approaches, a study is performed over the first 100 load steps of the analysis for mesh sizes of 0.1, 0.05, 0.025 and 0.0125 m, each for initial MP distributions of 6 2 and 8 2 MPs per element.Figure 9 shows the inner stretch of the cylinder after 100 load steps F I G U R E 8 Two-dimensional hyperelastic cylinder under internal pressure in the undeformed (blue) and deformed (red) configurations.

F I G U R E 9
Two-dimensional hyperelastic cylinder under internal pressure comparison inner stretch after 100 load steps for both traditional and new arc-length approaches for varying mesh sizes and initial MP distributions.
for varying mesh size and initial distribution.The results shown in Figure 9 does not represent a plot of accuracy and is not representative of the convergence of the analyses to the exact solution.By showing the inner stretch of the cylinder after 100 load steps, it is possible to see the sensitivity of the arc-length approaches (both the traditional and new approaches) to the discretisation of the problem with the scaled (new) approach displaying much less sensitivity to mesh refinement and changes in MP distribution compared to the unscaled (traditional) approach.It is clear to see that the traditional approach is consistently slower in terms of progression along the equilibrium path due to the fact that the displacement of the cylinder domain is limited by the large nodal displacements of partially filled elements accounting for the majority of the allowed nodal displacements in each load step.Using the new arc-length approach also produces more consistent results when varying the initial MP displacements as seen in Figure 9 where for the case of 1∕h = 10 1 , there is a significant difference between the calculated inner stretches with the traditional approach for the two MP distributions whereas the inner stretches for the new approach are much closer.As the mesh is refined, the difference in inner stretches for the two MP distributions reduce, producing virtually the same results while the traditional approach still shows a disparity in results for the two MP distributions.While both the traditional and new approaches produce the same equilibrium curve, it is clear that scaling nodal displacements results in a more efficient analysis which requires less load steps overall and allows for the use of coarser meshes and a reduction in the number of MPs.The benefits of scaling the nodal displacement in the arc-length constraint will also be shown in Section 4.3, including a load-displacement curve for the traditional and new arc-length approaches to demonstrate their respective progression along the equilibrium path.

Double notched plate under tension
The third numerical example is the plane strain, elasto-plastic analysis of a deep double-edge-notched plate under tensile load to demonstrate the ability of the MP-BALM implementation to track a snap-through response when modelling an elasto-plastic material.The specimen has an overall width, 2w, of 10 m and length, 2l, of 30 m. Two deep notches are made along the horizontal midline of the plate creating a ligament with a width, 2b, of 2 m.The geometry of the plate is shown in Figure 10A, which also shows a pressure, P, applied along the top and bottom edges of the plate.Due to the symmetry of the problem, one quarter of the geometry is modelled as shown in Figure 10A.The computational setup of the problem is shown in Figure 10B with roller boundary conditions applied along the side of the mesh where x = 0 m and only part of the section along the bottom of the mesh which corresponds to the ligament of the plate, y = 0 m and x ≤ 1 m.The physical domain is discretised using 2 × 2 evenly spaced generalised interpolation material points within each background mesh element.The problem is analysed using an F-bar implementation 26 in order to reduce the effect of volumetric locking.The plate material is modelled using a Henky model with a von Mises plastic model using with Young's modulus, E, of 206.9 a Poisson's ratio, , of 0.29, and a uniaxial yield stress,  y , of 0.45 GPa.An initial pressure of 50 MPa is applied by integrating the traction over a B-spline boundary generated along the top edge of the domain and the simulation is run until a top edge displacement of 0.2 m is achieved.
The simulation is run four times with multiple stages of mesh refinement, starting at a mesh size of 2 −1 m and down to a mesh size of 2 −4 m, with the mesh size reducing by a factor of 2 each time.Throughout each simulation, the reaction force along the nodes which are situated along the ligament (i.e., where the roller boundary conditions are applied along the bottom of the mesh) are calculated as the top edge of the plate is displaced due to the applied force.Letting R denote the total reaction force of the restrained edge of the ligament, the net vertical stress through the ligament is given by The normalised load-displacement response for this problem is shown in Figure 11.The snap-through behaviour is clearly shown for all mesh sizes with varying intensities, with only a slight snap-through behaviour seen for the mesh size of 2 −1 m.Also shown in Figure 11 is the Prandtl limit load 58 which is the analytical small strain limit load for this problem which is controlled by the stress at the notch in which  lim ≈ 2.97 y .It is not until the mesh is refined to a size of 2 −3 m that the load does not exceed the theoretical small strain load limit.Due to the geometry of the problem, a singularity is formed at the end of the notch which is difficult to model accurately.As the mesh is refined, this singularity is amplified and causes the issue of MPs separating from the base.These points of separation can be seen in Figure 11 by the sudden changes in load-displacement response for mesh sizes 2 −3 m (point A) and 2 −4 m (points B and C).This is especially prominent for the mesh size of 2 −4 m as this singularity causes the simulation to fail at a displacement of 0.15 m. Figure 12 shows the MP positions for the final load step of the simulation for each mesh size, the stress distribution is also shown.It is clear from Figure 12A,B that the singularity is not prominent as there is a smooth change in the MPs' vertical position moving away from the ligament.In Figure 12C the separation of the MPs away from the base can be seen as well as the increased separation between MPs around the singularity.Figure 12D shows the MP distribution for the end displacement of 0.15 m as this is the point at which the simulation fails, it is clear to see that large stress concentrations of the MPs at the end of the ligament as well as the disordering of the MPs around the singularity.
In order to demonstrate the need for special adaptation of the arc-length method for the MPM, the double notched plate problem is run with a mesh size of 2 −1 m and an initial pressure of 10 MPa over the top edge of the domain.This problem is run twice, using the traditional arc-length constraint in (3) where there is no scaling of the nodal displacements and with the arc-length constraint in (25) in which the nodal displacements are scaled based on the basis function contributions of associated MPs. Figure 13 shows the normalised load-displacement curve for the simulation run over 100 load steps.It is clear that the case with the traditional arc-length constraint begins to stall at a vertical displacement of approximately 0.02 m.After this point, the progression along the equilibrium path slows such that the physical material of the plate shows only slight deformation over each load step.This is due to the fact that there are partially filled elements and extremely large nodal displacements are calculated for nodes that are weakly influenced by the MPs.Referring to (3), these large nodal displacements mean that the load factor increment must be very small in order to satisfy the arc-length constraint.On the other hand, by weighting the calculated nodal displacements by the strength of the node's association with the MP domain, using (25), the effect of these large nodal displacements are mitigated, allowing the simulation to progress along the equilibrium path much further than the unscaled case.In this case, no ghost stabilisation is utilised in order to solely demonstrate the effect of the scaling matrix with the new arc-length approach.
Figure 14 shows the normalised top edge displacement of the double notched plate after 100 load steps for varying mesh sizes and MP distributions with the same problem setup as above to compare the traditional and new arc-length approaches.In order to improve the stability of the analyses, a ghost stabilisation parameter of  = 2.069 × 10 6 is used (five orders of magnitude less than the Young's modulus of the material).As ghost stabilisation essentially incorporates extra stiffness around the boundary of the domain, the solution of the unknown nodal displacements is improved and a reduction in large nodal displacements is seen, therefore, the stalling seen in Figure 13 is less frequent.However, as seen in Figure 14, the new arc-length approach still results in more efficient analysis with further

F I G U R E 13
Normalised load-displacement response over 100 load steps for double notched plate problem using the traditional arc-length approach and the new arc-length approach for a mesh size of 2 −1 m.
progression along the equilibrium path for all mesh sizes and MP distributions compared to the traditional arc-length approach.It is not until a mesh size of 2 −4 m that the traditional and new arc-length approaches produce the same displacement after 100 load steps.When using the new arc-length approach, the top edge displacement after 100 load steps is consistent for all initial MP distributions meaning that it is reasonable to reduce the overall number of MPs needed for the analysis compared to the traditional approach as well as having the option to increase the mesh size to reduce the overall runtime if desired.As with the hyperelastic cylinder example, the results shown in Figure 13 are not representative of the convergence to the exact solution with mesh refinement but a means to demonstrate the fact that the traditional arc-length approach is more sensitive to the problem discretisation compared to the new approach.
F I G U R E 14 Double notched plate comparison of top edge displacement after 100 load steps for both traditional and new arc-length approaches for varying mesh sizes and MP distribution.

F I G U R E 15
Problem geometry and basic mesh schematic.

Arch with applied point load
The next numerical example is a 2D arch with a plane strain assumption subjected to a point load to demonstrate the ability to track the equilibrium path of snap-through behaviour using a force controlled scheme.This problem is akin to the Crisfield arch 1 also presented by Hrinda. 59Here we adapt the arch to be modelled with continuum elements.A schematic of the problem geometry and mesh setup is given in Figure 15 in which the arch is formed with an inner length, L i , and outer length, L , from the origin of 9.5 and 10 m respectively.The height of the inner arch wall, h i , and the outer wall, h o , are set as 4.09 and 4.45 m respectively.The inner and outer radius of the curves formed by the arch are then 13.09 and 13.45 m respectively.The mesh is generated such that the nodes of the mesh align with the base of the arch in order to apply fixed boundary conditions to the nodes coincident with the base of the arch.The mesh is also extended below the line of the base of the arch to allow for further displacement past the initial height of the arch.The arch is modelled as a Henky material with a Young's modulus of 100 kPa and a Poisson's ratio of 0.25.The mesh is fully populated with each element containing 4 × 4 evenly spaced generalised interpolation material points.The points which are not contained within the domain of the arch are then discarded in a similar fashion to the cylinder problem in Section 4.2.A mesh of square elements each with size 2 −3 m is used for the simulation.An initial point load of 100 N is applied at the centre of the inner edge of the arch (indicated by the red arrow in Figure 15), this load is applied evenly between the bottom-most MPs in the arch, one either side of the centre line.The simulation is run until the MPs which have the load applied to them have displaced by 8.2 m, taking a total of 139 load steps to complete.Figure 16 shows the load-displacement response at the MPs at which the load is applied.Four stages are highlighted along the equilibrium path, the unscaled MP deformations at each of these four points are shown in Figure 17.
From Figure 16 it is clear to see that there is a snap-through response as the arch is pulled down and through the baseline.A peak load is encountered at point (a) (at load step 19), where the centre of the arch has been pulled down and

Slope failure
The final numerical example examines the deformation of a 45 • elasto-plastic, plane strain slope collapse due to gravitational loading.The geometry of the slope is shown in Figure 18 where roller boundary conditions are applied over the The difference in the load-displacement response between the arc-length and load controlled cases are stark and can be explained that at the point at which the paths begin to diverge that the material points will continue to displace with smaller increases in gravity than imposed by the load-controlled case.Implicit MPMs satisfy equilibrium using the state of the material points at the end of the current load step, where the stress carried by the material at their deformed positions are in equilibrium with the externally applied load.Once this equilibrium point has been found the points are then advected (displaced) through the mesh based on the incremental nodal displacements associated with the current step.However, once the material point positions are updated there is no guarantee that the stresses at the material points are still in equilibrium with the externally applied loads; this is the case for all MPMs, not just implicit MPMs or the implementation considered in this paper.This can be clearly seen by considering Equation (18) where the body forces are transferred to the nodes using the basis functions whereas the stress in the material is transferred using the derivatives of the basis functions and once the material points have changed position within the element it is likely, apart from for trivial cases, that there will be a non-zero residual at the start of the next step without increasing the imposed load.This non-zero residual will lead to additional material displacement in order to satisfy (18) without any increase in load.Increasing the load will lead to a different equilibrium path, as seen in Figure 19.The use of the MP-BALM ensures that the unbounded displacements seen by the load controlled schemes are contained due to the fact that the displacements of the domain are limited within each load step, this means that the simulation is stabilised around the critical point where the unbounded displacements begin.
Figure 20 shows the deformed configurations of the slope collapse analysis using a load controlled scheme (left column) and the MP-BALM (right column).For each analysis the final shear stress (top) and yield stress (bottom) states are shown.It is clear from Figure 20 that the MP-BALM results in far greater overall displacement with a larger failure F I G U R E 20 Final deformed configurations of the slope collapse problem for the load controlled analysis (left column) and MP-BALM analysis (right column).For each analysis, the shear stress (top row) and yield stress (bottom row) is shown.band being seen in the material and importantly at a lower load compared to the load-controlled case, identifying a more critical failure mode for the slope.

CONCLUSION
This paper presents an arc-length method for the MPM framework to expand the physical problems and loading scenarios that the MPM can solve, notably those involving snap-through and snap-back responses.Simply implementing the arc-length method into the MPM as it is given in the FEM approach is not suitable due to the fact that partially filled elements are a common feature in a MPM simulation.These partially filled elements mean that elements around the boundary of the domain will be poorly integrated, understiff and have the potential to experience extremely large nodal displacements.While the effect of the large nodal displacements are reduced when mapping nodal displacements to the MP positions, the effects are seen within the arc-length constraint as the nodal displacements account for the majority of the allowed displacement based on the prescribed arc length.Therefore, the arc-length constraint is reformulated such that it is based on MP displacements, essentially scaling the nodal displacements based on the node's influence on the physical domain.As the MPM involves the resetting of the background mesh at the end of each load step, nodal displacement information is lost due to the fact that the MPs have changed positions and the nodal displacements cannot simply be carried over from the previous load step.This means that in order to use the secant path predictor, an estimate for the nodal displacements for the previous load step based on the current MP positions are needed.This nodal displacement reconstruction is carried out at the start of each load step and uses a least squares weighted residual technique.
A simple benchmark test of a 2D column under self weight is given to show that the MP-BALM is able to model standard problems that the MPM with a linear solver has been proven to solve.The expansion of a hyperelastic cylinder which demonstrates a snap-through response is shown and a comparison was made between the traditional arc-length scheme and the MP-BALM scheme presented in this paper.It is clear that through scaling the nodal displacement increments, a more robust and efficient analysis is achieved.More complex problems are shown to demonstrate that the MP-BALM can model problems with more complex deformation and model snap-through responses with problems of a double notched plate under tension, the deflection of a shallow arch and collapse of a slope including material softening.The final case has demonstrated the ability of the method to control the displacement response of a gravity-driven problem with unbounded deformation, typically of geotechnical slope collapse, which is not possible with conventional load controlled quasi-static methods.
As is the case for many arc-length papers in the current literature, the arc-length method has only been presented for the problems with load control.The case of displacement controlled problems is much more involved with an alternative arc-length constraint equation needed, please see Pretti et al. 16 for more details on the formulation of a displacement controlled arc-length method.

1
Equilibrium paths of (A) increasing response, (B) snap-through response and (C) snap-back response.

F I G U R E 7
Two-dimensional hyperelastic cylinder under internal pressure.Plot: inner stretch against normalised applied internal pressure comparing calculated results (solid black line) and analytical solution (red crosses).Inset: diagram of cylinder problem setup showing initial inner radius and outer radius and deformed inner radius.

F
I G U R E 10 Double notched plate: (A) Problem geometry and (B) computational mesh and MP distribution.
U R E 12 MP positions and stress distribution around the notch singularity at the end of the simulation for mesh sizes of (A) 2 −1 m, (B) 2 −2 m, (C) 2 −3 m and (D) 2 −4 m.

F
I G U R E 16Load-displacement response of elastic arch with applied point load.Unscaled MP deformation for shallow arch problem at stages (A) load step 19, (B) load step 48, (C) load step 84 and (D) load step 139 (see Figure16for position along equilibrium path).bulging of the arch can be seen on the left and right sides of the arch.After this point, the snap-through response can be observed as the material softens and deforms through point (b) (at load step 48), is pulled through the base of the arch until the next stationary point at point (c) (at load step 84).After point (c), the material begins to stiffen again as the load required to further deform the arch begins to increase until the centre of the arch has displaced by 8.2 m at point (d) (at load step 139).

F I G U R E 19
Load-displacement responses at specified positions around edge of slope domain throughout elasto-plastic collapse.a significant increase from the maximum displacement seen by the load controlled case, which fails at a position C displacement magnitude of approximately 0.19 m.The MP-BALM analysis took a total of 477 load steps to achieve a position C displacement magnitude of 0.5 m.The load-displacement curves of the MP-BALM analysis are shown in Figure19and are depicted as solid black lines in each plot.It is worth reinforcing the point that the reported displacements in Figure19are the magnitudes of the displacements of the positions in question.In terms of local deformation, it is clear that positions A and B display local snap-back responses and positions C and D undergo snap-through responses.
Double plate load-displacement curves for various mesh sizes.