Multiple discrete crack initiation and propagation in Material Point Method

versatility

On the other hand, the Material Point Method (MPM) [13] and its derivatives (e.g.[14][15][16][17]) allow for solving large deformation problems without mesh entangling while bearing substantial resemblance to FEM (see e.g.[18,19]).MPM discretizes a domain into particles arranged on top of a static background grid.tensile stress at an infinite distance  J−Integral energy release rate calculated via numerical J-integral  elastic modulus  ,J−integral mode I stress intensity factor calculated from  J−Integral [𝐷] standard 2D elastic constitutive tensor presented as matrix in Voigt notation All state variables are stored in the particles, while the problem solution is obtained on the grid nodes.Shape functions map the variables from the particles to the grid and, after the solution, back to the particles, to which particles map state variables to and from via shape functions.
In MPM, numerous contributions successfully implemented fracture using, among others, phase field [20][21][22][23][24][25][26], damage field gradient [27], XFEM-like extension of MPM [28,29], plasticity with softening [30,31], and even the use of irregular mesh (grid) [32].Nonetheless, an MPM-native fracture simulation method called CRAMP [33] introduced a versatile and intuitive implementation of discontinuity due to fracture.CRAMP utilizes the concept of multiple velocity fields via parallel grids such that the adjacent particles may move and deform independently to simulate discontinuity.The concept of multiple grids introduced in CRAMP also influenced other MPM algorithms beyond the simulation of cracks, e.g.Damage Field Gradient [34] and Domain of Interaction [35,36].
In this paper, we introduce significant improvements to CRAMP so that it is capable of multiple crack initiation, propagation, and merging while retaining the benefit of the Material Point Method, i.e. the lack of fixed elements allowing for straightforward simulation of large deformation problems.The novel changes to the algorithm include: 1. Additional grid to accommodate arbitrary number of crack paths 2. Crack initiation via Rankine criterion 3. Crack propagation evaluation adapted to also work with Rankine criterion 4. Minimization of the number of background grids via a modified grid mapping method 5. Detection and termination of crack path reaching the end of the material domain 6.Capability for simulating crack merging The proposed algorithm is implemented shown to work with Convected Particle Domain Interpolation MPM (CPDI) [15].. Nonetheless, implementation to GIMP [14] or the original MPM [13] would only require simplifying the interpolation functions.
The contribution is novel, as it efficiently supports an arbitrary number of cracks that can merge and initiate at any point of the domain during the simulation, while e.g.[37] managed to simulate only two crack paths.Other studies [35,36,38] require the evaluation of all particle pairs relative to all existing crack paths, therefore the computation effort grows quadratically with respect to the number of particles.We instead retain the original grid assignment method of CRAMP with an additional algorithm that reduces the density of the parallel grids, therefore requiring a similar amount of calculations compared to the original CRAMP algorithm.
A total of 8 verification test cases show good convergences of the displacement field, stress field, energy release rate, energy dissipation, and load-displacement curves, along with demonstrations of successful crack initiation, propagation, and merging.

CRAMP and multiple grids
CRAMP allows adjacent particles (which normally interact as a continuum) to separate, resulting in a discontinuous displacement field.This is achieved by associating particles on the different faces of the crack with different velocity fields and different grids, leading to independent velocity fields for the separated particles.
Let us first consider only one grid (i.e.without CRAMP) and the original MPM mapping of momentum and mass between particles and nodes: where  is the particle number, Node is the node number,  is the time,   is the number of particles,  Node is the number of nodes, ⃗   Node is the nodal momentum, ⃗    is the particle velocity,   is the particle mass,   Node is the nodal mass, and   Node, is the MPM shape function.
In MPM algorithm the particles are connected through the mapping into common nodes.CRAMP breaks some of these connections by associating relevant particle-grid mappings to different grids, therefore avoiding node sharing and ultimately allowing displacement-field discontinuity.Fig. 1 compares the original MPM (a) and CRAMP (b) mapping of particles.CRAMP models the discontinuity between particles 3 and 4 by rerouting the particle-node mappings that cross the discontinuity into grids #2 and #3.The mapping paths are represented by the blue and yellow lines.
CRAMP modifies Eq. ( 1), with an introduction of a Kronecker delta to map particles to the correct independent grids: where Grid can either be 1, 2, or 3; the nodal mass   Node,Grid is arranged as diagonal matrix; and  =  (, Node) returns which background grid number  is involved with mapping between particle number  and node number Node, whose possible values are shown in Eq. (3) (see [33]):  After computing the nodal momenta and masses, particle position update -as an example -follows Eq. ( 4)

CRAMP with CPDI
The simulations in this paper utilize the Convected Particle Domain Interpolation (CPDI) [15], a version of MPM with the particle domain shaped as a parallelogram, as shown in Fig. 2. To allow for fracture modelling, the concept of multiple grids needs an analogue implementation for CPDI.A top hat function defines the particle domain, described in Eq. (5).The shape function is acquired by integrating this top hat function with the standard ''tent'' grid shape function (  ( ⃗  ) ) from the original MPM [13], as shown in Eq. (6).
The integration of Eq. ( 6) is computationally impractical, thus CPDI utilizes an alternative shape function.The alternative shape function approximates the exact solution with significantly-reduced computational costs.For 2D application, let us first approximate the integrand of Eq. ( 6) with an approximation: While where ⃗    is the position of the  -th corner of particle )) where  Node is the standard MPM ''tent'' shape function as shown in Eq. ( 8), and  is the particle domain corner number.Therefore, mappings of e.g.mass and momentum from particle to grid are as shown in Eq. ( 14): where  is the particle number, Node is the node number,  is the timestep number,   is the number of particles, ⃗   Node is the nodal momentum, ⃗    is the particle velocity,   is the particle mass,   Node is the nodal mass (arranged as diagonal matrix), and   Node, is the CPDI shape function.
The determination of  =  (, Node) only requires the positions of nodes, crack path(s), and the centres of particle domains, i.e. particle positions.Therefore, the procedure remains the same as proposed in the original CRAMP.This approach is similar to [34] in terms of separating the particles and their domains to different sides of the discontinuity, as illustrated in Fig. 3.
The introduction of variable  to particle-grid mapping and additional dimension  of nodal values now allow for multiple background grids to work with CPDI.For example, Eq. ( 14) now becomes Eq. ( 15 Fig. 3.The separation of crack particles and their domains (CPDI); whether a particle domain is above or below a crack is determined by the position of the particle relative to the crack path.
where Grid is the grid number, Node is the node number, and  (, Node) determines the relevant grid number for mapping between particle  to node number Node.Note that particles associate with exactly one grid, while its particle-node mapping may associate with more than one.With a similar principle to Eq. ( 4), the mapping from the parallel grids to a particle, e.g. for updating the particle position, follows Eq. ( 16) where Δt is the time-step size,   is the number of nodes, and   is the number of grids.To satisfy the Courant stability condition, our Forward-Euler explicit time integration scheme requires time-step to follow  <  cr = (∕), where  is the time-step size,  cr is the critical time-step size,  is the cell size, and  is the speed of elastic wave propagation.In practice, we take half of the critical value, i.e.  = 0.5 (∕).

Crack initiation and propagation mechanisms
The initiation and propagation algorithm we propose utilizes the Rankine criterion [39].A new crack path of length equal to the grid cell size initiates or an existing crack path propagates when the maximum principal stress is higher than the defined material strength (see Eq. ( 18)).Previous cracking-particle method implementations evaluate and initiate new crack surfaces at particles (see e.g.[40][41][42][43]).However, in e.g.CPDI, the particle domain is a parallelogram that cannot be split easily.Therefore, this paper proposes crack initiation to occur at nodes (Fig. 4), while propagation simply extends the existing crack path (Fig. 5).The additional benefit of monitoring stress at the nodes instead of particles is better accuracy, as equations of motion are solved at the nodes, with nodal values free from fluctuations such as checkerboarding (see e.g.[44]) Nodal stress is calculated via the same method as for other variables (Eq.( 15)): where  is the particle number, Node is the node number,  is the timestep number,   is the number of particles,  , is the stress tensor components at particle ,   Node,Grid is the stress tensor components at node number Node and grid number Grid.Through continuous evaluation of the maximum principal stress at every node, a new crack initiates of length equal to the grid cell size, when the maximum principal stress reaches a certain value, 2 i.e. no longer satisfies inequality (18). where 2 Unless node is within the prevention zone and satisfies the prevention rule, explained in the next subsection.

Fig. 4.
A crack initiates when the maximum principal stress at a node exceeds the material strength.The length of the crack segment is equal to the cell size.
While   is the tensile strength,   ,   , and   are the components of 2-dimensional stress tensor.
The direction of the crack surface is perpendicular to the direction of maximum principal stress, i.e.Eq. (20).
where  crack is the direction of the new crack surface, and   1 is the direction of maximum principal stress, i.e.: Crack propagation utilizes a different mechanism.By taking a cell as a single quadrilateral finite element and identifying its nodal displacements, the stress field can be calculated at any arbitrary location inside the cell.In particular, it allows for calculating the value of maximum principal stress at the crack tip (see Appendix D for more details).Propagation occurs when the maximum principal stress reaches the material strength, where the direction is perpendicular to the direction of the maximum principal stress.Appendices B and C give details of the crack initiation and propagation algorithms.

Cohesion
The use of cohesive forces in CRAMP was introduced by [45], therefore allowing energy dissipation during crack propagation.Cohesion is applied as traction on crack surfaces near the crack tip.Since CRAMP utilizes zero-mass particles to represent crack surfaces, the traction field on the crack faces is discretized into forces acting on the crack particles, as illustrated in Fig. 6 The test cases presented in this paper utilizes the bilinear traction law (see Fig. 7), where the area under the curve is equal to the prescribed fracture energy   .

Prevention zone
Even though the maximum principal stress at a node exceeds the material strength, crack initiation may not proceed if the said node is close and similarly-angled to an existing crack path.This concept is similar to the prevention criteria implemented in [41,46] originally for the Element-Free Galerkin (EFG) method, with the goal of reducing spurious crack formation.
In the proposed method, the prevention zone surrounds the existing crack paths.Crack initiation outside of any prevention zones is always realized, as Fig. 8(a) shows.Crack initiation within this zone is prevented if the initiated crack would have similar orientation to the nearest existing crack path, illustrated in Fig. 8(b).However, a crack can still initiate within a prevention zone if the initiation angles the crack significantly compared to the nearest crack path segment, as shown in Fig. 8(c).For the verification cases presented in this paper, the angle is considered significantly different if the difference is more than ∕4.The prevention zone itself spans 1.5 times grid spacing distance.
In the unlikely event that the maximum principal stress at 4 nodes of a cell reaches the defined tensile strength at the same time, a crack will initiate at one at random.Spurious crack initiation at the other nodes is avoided due to the aforementioned prevention rule.Nonetheless, the positional uncertainty from this randomness can be lessened through mesh refinement.

Multiple cracks
An arbitrary crack path requires definition of 3 parallel grids.To guarantee the same level of discontinuity representation for arbitrary multiple crack paths as that of single crack cases, each additional crack path requires the split of each grid into another 3. Therefore, the presence of a second crack path requires 9 grids, a third 27 and so on (see Fig. 9).This exponential growth of grid number would result in an extremely-high computational cost for even a moderate number of crack paths.
To avoid this exponential growth and maintain a reasonable number of background grids, this paper proposes a Crossing Profile method.Instead of splitting each grid into three with the presence of a new crack path, the algorithm identifies unique particlecrack-node interactions.Each unique interaction reserves a grid for use by all entities with calculations featuring such interaction.Appendix A explains the method in detail.As an example, the verification case of a plate with multiple initial cracks simulates 10 crack paths with 21 grids at the initial configuration and 51 at the end of the simulation.

Merging and edges
For every active crack tip, let us define a point ahead of it called a ghost tip.The ghost tip is positioned at a distance equal to twice the defined propagation length in the direction of the crack tip trajectory.If the ghost tip detects crack path crossing, the real crack tip extends towards the crack path to simulate merging.The extended tip is then defined as inactive to prevent further propagation.If the ghost tip detects zero nodal mass (i.e. when the crack tip leaves the domain), the real crack tip extends towards the ghost tip plus an additional one-cell-sized extension.In total, this extension ensures complete separation between crack surfaces.The extended tip is then also marked inactive.Figs. 10 and 11 show the mechanisms of crack merging and crack propagation near an edge.
There is no prescribed behaviour regarding the change in grid layers and particle-grid mapping when merging occurs.They are governed by the proposed crossing profile method (see Appendix A).The resulting number of grid layers and how they map to the particles are emergent behaviours from the implementation of said method.

Results
We verified the crack-capable CPDI solver with the vortex case [47] and simulations from [47,48], with the results and convergence being very similar, see Figs. 14 and 15.
After the first non-fracture-related verification test, this paper presents six more cases to further verify the proposed method.The first one is a convergence study on the stress field in the simulation of an infinite plane with a horizontal crack.The second case shows crack initiation and propagation, while the third case shows crack propagation and merging.In the fourth case, we simulate a plane with multiple initial cracks under biaxial tension.The final fifth case shows the capability of a discrete crack in MPM for large crack opening and closing through the simulation of a ring.The ring has initial cracks arranged radially at the outer surface that propagate during simulation.Table 1 shows the material properties used in the simulations.

Vortex problem
The vortex problem [49] is a manufactured solution of displacement field throughout a ring-shaped problem domain.From the manufactured solution, the initial conditions, body forces, and tractions are derived and applied to numerical simulations to replicate a given displacement field.In the simulation, the ring is sheared in a way that zero displacements occur at the inner and outer surfaces of the ring while tangential displacement occurs in between.Fig. 12 illustrates the intended displacement field of the vortex problem.This paper simulated the vortex problem to verify that the CPDI formulation remains intact and correct.The simulation uses the average of 4 particles per cell arranged radially on a cartesian grid.Figs. 13 show the deformed shape at maximum displacement ( = 0.5 second) with a grid size of 0.03125 metre using the CPDI proposed in this paper (a) and that used by Tran, et al. [48] (b).Both plots appear to be qualitatively similar.To confirm their quantitative agreement, the convergence curves of three CPDI codes are compared: the code used in this paper, Tran, et al. [48], and Kamojjala et al. [47].The convergence study takes the root mean square (RMS) of the displacement field error at maximum displacement ( = 0.5 second) with five different grid sizes: 0.25, 0.125, 0.0625, 0.03125, and 0.015625 metre (deformed shapes in Fig. 14).Eq. ( 22) describes the calculation of root-mean-squared displacement error, a comparable quantity across simulations with a different number of particles.Fig. 15 shows the resulting convergence curve and comparison to those presented in [47,48].confirming the correctness of the CPDI solver in the code used in   this contribution.

RMS
where RMS is the root mean square, ⃗  analytical  is the displacement of particle  from analytical solution, ⃗  MPM  is the displacement of particle  from MPM simulation, and  particles is the number of particles in the MPM simulation

Horizontal crack in an infinite plane
The analytical solutions for both the stress field and the energy-release rate for a horizontal crack in an infinite plate under uniaxial tension are well-established for small deformations (see e.g.[50]), therefore the MPM simulation employs small deformations as well.
The MPM simulates the infinite plane by modelling a finite plane with traction applied around the edges calculated from the analytical solution (see Fig. 16).The applied edge traction and maximum-principal-stress field at equilibrium are shown in Fig. 17.
Qualitatively, the results show stress concentrations around the crack tip with the well-known ''butterfly'' pattern.To provide a quantitative comparison with the exact solution we evaluate the error in strain energy, defined in (23) for various mesh densities.
where  numerical is strain energy calculated from MPM solution; and  analytical is strain energy calculated from the analytical solution.
The error vs grid cell size on a log-log plot shows linear convergence with slope at maximum refinement levels of 1.01 (0.981 using all data points), as shown in Fig. 18(a), with the slope being very close to the theoretical FEM convergence rate of [51].22), RMS calculated the same way as [47]) convergence curves of vortex problem from literatures are closely-matched to that from this paper.
Additionally, this convergence study compares the stress intensity factor via numerical J-integral to the analytical value.Eq. (24) shows the analytical value of Mode I stress intensity factor for an infinite plane: where  ,analytical is the analytical mode I stress intensity factor,  is the Poisson's ratio, while  ∞ and  are defined in Fig. 16.
The numerical energy release rate is calculated via the path-independent dynamic J-integral equation for CRAMP established in [52].The stress intensity factor is then calculated following Eq.( 25):     where  J−Integral is the energy release rate calculated via numerical J-integral,  is the Poisson's ratio,  is the elastic modulus, and  ,J−integral is the mode I stress intensity factor calculated from  J−Integral .A normalized error for the numerical J-integral is calculated by comparing to the analytical value, as Eq. ( 26) shows: Plotted on a log-log scale, the error shows a linear convergence curve with the slope of 1.04 at the maximum refinement levels and 1.11 via the least square of all the data points, as shown in Fig. 18(b).These values exceed slightly the convergence rate of 1 given in [51].

Crack propagation simulation convergence
The propagation convergence test involves simulations of a 30 cm X 30 cm plate with a 7.5 cm initial crack on one edge at different refinement levels.Displacement-controlled top and bottom boundaries induce Mode I crack propagation.The convergence test uses is repeated with different grid cell sizes: 100/3 mm, 100/5 mm, 100/7 mm, 100/9 mm, 100/13 mm, and 100/17 mm.Figs.19 and 20 show the schematic and loading scenario of the problem respectively.
Energy dissipation in the simulations is calculated by keeping track of boundary condition work, strain energy, and kinetic energy, which are shown in Table Table

Crack initiation simulation
The following simulation shows the capability of the method to simulate crack initiation.For this case, a rectangular specimen with a notch is loaded in a displacement-controlled tension from both the top and bottom surfaces.Fig. 22 shows schematically the problem.
The presence of a notch via particle removal emulates that of an initial crack.The initial configuration features no crack paths.Additionally, in the context of the notch tip itself, the right-angled corners do not act as particular points of singularity due to a velocity-field-smearing effect (see e.g.[53]).Although not infinitely sharp, the resemblance of the notch to a crack results in a significant stress concentration build-up at the tip.Fig. 23(a) shows the tensile stress builds up within the notched area where a discrete crack initiates, as shown in Fig. 23(b), then initiated a discrete crack, shown in Fig. 23(b).As loading progresses, the crack continues to propagate until the other end of the domain (Figs.23(c) and 23(d)).After reaching the edge, the crack tip detects zero nodal mass (Fig. 23(d)), thus the crack propagation terminates.Fig. 24 shows another visualization of the crack path development.
Note that the resulting crack path is not completely straight due to the use of dynamic analysis with explicit time integration, therefore calculations involving crack propagation are interfered with by stress waves passing through the immediate area of the crack tip.Nonetheless, the crack propagation averages into the correct propagation direction.This issue could be significantly reduced by employing non-local stress evaluation around the crack tip.

Crack initiation simulation (large deformation)
The following simulation is the same as that in Section 3.4, except that we aim to simulate crack initiation and propagation with large deformations.This is achieved by significantly lowering the elastic modulus of the material.For this simulation, we define the elastic modulus of 80 MPa, 100 times less than the previous crack initiation simulation.
Figs. 25 and 26 shows the simulation result.Similar to the previous one, our proposed algorithm predicted the correct crack initiation location at the tip of the notch, albeit with the large deformation already occurring.Propagation then continues with a large crack opening.Note that the crack surfaces are quite jagged following the large crack opening.This fact has no physical consequences, since cohesion, which is determined by the crack surface positions, is already zero at this point.

Simulation of crack merging
The crack merging case involves two initial cracks.Crack path 1 goes towards the edge, while crack path 2 spans across the inside of the specimen vertically.Both the top and bottom edges are loaded via displacement control to induce tension.Fig. 27 shows the initial configuration and loading scheme of the crack merging case, while Figs.[28][29][30] show the simulation results.The displacement loading scheme follows that defined in Fig. 20.
In the simulation, stress concentrations initially build up around all three crack tips.Initial crack 1 then starts to propagate, eventually merging into crack path 2. Moments before merging, both tips 2 and 3 started to propagate together.The simulation ends with the specimen completely split into three pieces.Nonetheless, this result is unstable.A slight rotation to the vertical initial crack causes one of its tips to be more likely to propagate while the other propagates in the opposite direction.Figs.31-35 shows the simulation results with a slightly-angled vertical crack and the comparison to those with perfectly-vertical initial crack.

Simulation of a plate with multiple initial cracks
The plate with multiple initial cracks under biaxial tension introduced by Budyn et al. [54] has 10 initial cracks of equal lengths that propagate.Shown simulation has displacements prescribed on the four edges of the plate, while the Rankine criterion determines crack propagations and their directions.Figs.36 and 37 show the initial cracks and loading conditions, while Fig. 38 shows the propagation result by Sutula et al. [55].Figs.39 and 40 show the MPM result that qualitatively replicate most of the crack patterns.We do not expect an exact replication of the crack pattern as the MPM simulation employs dynamic explicit analysis, while the results presented in [55] are an outcome of an implicit quasi-static simulation with an energy minimization technique.
Throughout the simulation, the crack propagation reaches several bifurcation points (see e.g.[56]).As expected, the outcomes of said bifurcations are affected by both the number of material points (see Fig. 41) as well as the dynamic effects present in the explicit MPM simulations.The introduction of numerical damping lessens the impact of elastic waves, therefore comparison to non-damped results can illustrate the influence of the aforementioned dynamic effects, see Fig. 42.Damping implementation follows that in [33], with medium damping halving particle velocity every 10 s, while high damping does every 1 s.
Additionally, this simulation shows the effectiveness of the proposed Crossing Profile method in minimizing the required number of grids.The discontinuity caused by 10 cracks at the initial configuration requires 21 grids to simulate.The more complex crack pattern at the end of the simulation involves 51 grids.Fig. 36.Initial cracks and loading scheme of the plate with multiple crack problem by Budyn et al. [54] with results from Sutula et al. [55].

Cracked ring
To show the capability of MPM performing simulations with multiple cracks in the large displacement regime and large crack openings, this paper presents a simulation of a ring with eight radially-arranged initial cracks at its outer surface, see Fig. 43.The ring expands and contracts as a result of a displacement-controlled traction in a push-and-release cycles along the inner surface, see Fig. 44.The deformations lead to the initial crack propagation, opening and closing.The cracks propagate only in the initial material (Figs.46(a) and (b)).At maximum load, Fig. 46(c) shows the crack opening.The blue and red lines show the cracked surfaces (see [33] on the method to track cracked surfaces).Afterwards, traction is gradually reduced to zero, closing the crack surfaces.Apart from the propagated initial cracks, the particle configuration in Fig. 46(d) resembles that of the initial configuration.
Simulations with two mesh densities are carried out, the so-called coarse mesh and fine mesh (double and half the cell size, respectively), Figs.45(a)-(f) and Figs.46(a)-(f).The coarse mesh successfully simulated the large crack opening for both load cycles, and the final configuration shows a similar state to the initial configuration.The fine mesh also managed to simulate the large crack opening in the 1st load cycle, however, a significant distortion can be observed after the second load cycle.
Both the coarse and fine mesh undergo another load cycle for the 3rd time, which results are shown in Figs.47 and 48.At the end of the cycle, the coarse mesh shows mild distortion, while for the fine mesh the distortion is significant enough for the crack surface tracking to degrade.From these results, the reason for such instability is hypothesized to be stress concentration around the crack tip crossing through the cells.This hypothesis explains the higher distortion severity in the simulation with a fine mesh, since a finer mesh captures the stress concentration more prominently whilst at the same time the particles move through more cells.
Confirming the hypothesis requires further study.

Conclusions
This paper extends the capability of CRAMP to simulate an arbitrary number of crack paths by introducing additional grids while maintaining the original CRAMP algorithm of particle-grid assignment.A modification to the assignment algorithm (called Crossing Profile method) reduces the number of necessary grids from exponential to linear with respect to the number of crack paths without significant addition of computational cost.This paper also introduces a method of evaluating crack initiation and propagation via the Rankine Criterion.As crack initiation, propagation, and multiple cracks are possible, we also proposed an algorithm for handling crack merging and termination at the edge of the material domain.All of the proposed methods and modifications are implemented in the CPDI version of MPM.Additionally, as an inherent advantage of MPM with large deformations, this paper also shows the ability of the method to model large crack openings.
A wide range of test cases was used to verify the proposed method.The first two cases with analytical solutions show the expected convergence rate of quantities related to the stress field, displacement field, load-displacement curves, and energy dissipation.The next three cases demonstrate the method capabilities related to crack initiation, propagation, and merging, with simulations resulting in the correct crack patterns.The next case shows a simulation of a multiple crack propagation test case from [55] with obtained results being in reasonable agreement with the original result.Additionally, the method only requires a total of 21 and 51 grids to   [55] with energy minimization and 1200 × 1200 elements and the proposed method with dynamic MPM and 65536 particles (128 × 128 cells), which show similarity and some differences in propagation pattern.simulate 10 crack paths at the start and end of the simulation respectively.The last case exhibits the capability of discrete crack modelling and MPM to simulate large deformation problems with large crack openings.This case would be difficult to replicate in FEM, as it involves large deformations and dynamics.
Due to the use of Rankine criterion, the proposed method in its current form is restricted to elastic materials with tensile failure mode only.An additional consequence is the limit of prescribing only one fracture energy value.Nonetheless, the method can be easily expanded by employing more complex failure criterion, potentially with several failure modes.Elasto-plastic materials can also be implemented, with material failure predicted by e.g. the evaluation of acoustic tensor.With additional failure modes, several fracture energy values can be defined.
In conclusion, the proposed method introduces the capability of crack simulation, including multiple cracks, crack merging and crack closure, into the CPDI Material Point Method.Such implementation boasts an additional capability of fully-dynamic crack simulations with large displacements and finite strains due to the inherent capability of the parent method.Furthermore, the proposed method can easily be implemented into other MPM versions, as CPDI can straightforwardly be simplified into e.g.GIMP or even the original MPM through a special choice of functions and simplification of the algorithm.Profiling: The method evaluates and assigns a profile for each particle-grid mapping regarding whether the path crosses each existing crack paths from the top, from the bottom, or not crossing at all.This requires minimal modification to the original CRAMP algorithm itself, more specifically to the evaluation of particle-node and crack path crossing (see Appendix A, Task 1 of [33]).
The profile of a particle-grid mapping is written as integers, each digit either 0, 1, or 2, with the number of digits being equal to the number of crack paths.The first digit describes how the particle-grid mapping crosses crack path number 1: the number 0 represents no crossing, number 1 for crossing from above, and number 2 for crossing from below (see Appendix A, Line Crossing Algorithm of [33] for the definitions of above and below crack paths).The second digit then describes how the particle-grid mapping crosses crack path number 2, the third digit for crack path 3, and so on.For example, a particle-grid mapping which crosses crack path 1 from above, crack path 2 from below, and not crossing crack path 3 or 4 would be: 1200.
Grouping: In this process, the unique profiles of particle-grid mappings are assigned a unique grid.Any algorithm to achieve this goal is acceptable, although this paper recommends grouping to occur simultaneously with profiling, therefore eliminating the necessity of storing all the profiles for each particle-grid mapping.After identifying the profile, After profiling the particle-grid mappings, .The method also aims to group the particle-grid mappings based on the identified profiles.The grouping process occurs right after each particle-grid mappings are profiled, eliminating the need to store the profiles of each particle-grid mapping.We achieve this by keeping track and numbering unique profiles in a table as profiling takes place.Each unique profile corresponds to a unique grid, i.e. the profile number is equal to the grid number.If the recently-profiled particle-grid mapping matches one on the table, a grid with the profile number (and therefore the grid number) is assigned to the mapping.Otherwise, a new entry and numbering is added to the table.
Let us take the case shown in Figure Fig.A.1 as an example.We may start with the mapping of particle 1 to node 1, and call it mapping 1→1.The path between the two crosses no crack paths, i.e. profile 00, therefore the default grid number 1 (i.e.grid 1) is assigned for the mapping.Mapping from particle 1 to node 2 crosses crack path 1 from above, therefore profile 10 is assigned to mapping 1→2.Such a profile has never been discovered before, thus a new grid number 2 (i.e.grid 2) is defined and assigned for mapping 1→2 and any other mapping with the same profile.Mapping from particle 2 to node 1 crosses no crack paths, i.e. profile 00, therefore the default grid 1 is assigned for mapping 2→1.Mapping from particle 2 to node 2 crosses crack path 1 from above, thus having the same profile 10.This profile has grid 2 assigned to it already.Therefore, mapping 2→2 goes to grid 2. Mapping from particle 3 to node 1 crosses crack path 1 from below, which means the profile for mapping 3→1 is 20.The profile is unique, so a new grid -grid 3 -is assigned for this particle-grid mapping.The process then continues with the rest of the particle-grid mappings.
Table A.1 shows all of the assigned grid numbers to the particle-grid mappings and their profiles, while the reason regarding the choice of grid number is shown in Table A   Mapping from particle 3 to node 1 crosses crack path 1 from below, a profile never defined before (note that crossing a path from above and below are treated differently).Therefore, a new layer number 3 is defined for mapping 3→1.Mapping from particle 3 to node 2 crosses no crack paths, therefore the default grid 1 is assigned for the mapping.Mapping from particle 4 to node 1 crosses crack path 1 from below and crack path 2 from below.Such profile a is not defined yet, therefore a new grid 4 is defined and assigned for mapping 4→1.Mapping from particle 4 to node 2 crosses crack path 2 from below.The profile is not defined yet, therefore a new grid 5 is defined and assigned for mapping 4→2 (see Fig. Tables A.1 and A.2 summarize the result of the crossing profile method for this example case.
Note that at the end of every timestep, the generated crossing profile is discarded, and a new one is generated for the next timestep.Nonetheless, due to the minimum disruption of the crossing profile method to the original working principle, information is retained at the same level as CRAMP.Even though the proposed method discards grid mapping of every timestep, essentially all the variables are carried by the particles.
For a typical case with the crossing profile method, the number of grids grows linearly with additional crack paths.Additionally, the crossing profile method is relatively light and can easily be parallelized.
One possible drawback of this method is the lack of physical meaning of the grid number, as -unlike in the original CRAMPthe number no longer reflects particle or node position with respect to any nearby crack paths.Only grid 1 retains a distinct physical meaning, which is the default where particle-node mapping crosses no crack paths.Example of crossing profile method with 4 particles, 2 nodes, and 2 crack paths; arrow indicates crack path direction (see [33] on the definition and purpose of crack direction).

Fig. 1 .
Fig. 1.Demonstration of CRAMP implementation in 1D, where two 3-particle clusters are separated by a discontinuity.CRAMP allows independent displacement and velocity fields of the two clusters by defining two extra grids and rerouting mapping of particles around the discontinuity such that no particle pair from different sides share common nodes.

Fig. 2 .
Fig. 2. CPDI parallelogram particle domain of particle p and the ability to deform from initial configuration (0) to current configuration (t).

Fig. 5 .
Fig. 5. Crack propagation by treating one grid cell as one quadrilateral Finite Element.

Fig. 6 .
Fig. 6.Traction on crack surfaces (a) and equivalent forces acting on the crack particles (b).

Fig. 8 .
Fig.8.Considering prevention zone around an existing crack, initiation is either: (a) allowed, as it happens outside of the prevention zone; (b) prevented, as it happens inside the prevention zone with similar angle to nearby existing crack path; or (c) allowed, as it is significantly angled from nearby existing crack, even though it happens inside the prevention zone.

Fig. 9 .
Fig.9.Particle-grid mapping with two crack paths necessitates a total of 9 grids;  number of crack paths would require 3  grids, unless the proposed Crossing Profile algorithm is implemented to reduce the total.

Fig. 10 .
Fig. 10.Crack merging occurs when a ghost tip detects another crack path.

Fig. 11 .
Fig. 11.Crack termination at the edge occurs when a ghost tip detects a node with zero mass.

Fig. 12 .
Fig. 12. Manufactured solution of the vortex problem, where body forces are derived from the prescribed displacement field [47].

Fig. 13 .
Fig. 13.Deformed shape of the Vortex problem with Convected-Particle-Domain-Interpolation code in this paper (a) matches with the result with Convected-Particle-Domain-Interpolation code from [48] (b).

Fig. 14 .
Fig. 14.Deformed shape with grid sizes of (a) 0.25 m, (b) 0.125 m, (c) 0.0625 m, (d) 0.03125 m, and (e) 0.015625 m; simulation with our CPDI code shows qualitatively-improving approximation of the vortex problem with further refinement (e.g.reduced shear distortion in the middle and near the surfaces of the ring), while the quantitative convergence is shown in Fig. 15.

Fig. 15 .
Fig. 15.Displacement-field error RMS (see Eq. (22), RMS calculated the same way as[47]) convergence curves of vortex problem from literatures are closely-matched to that from this paper.

Fig. 16 .
Fig. 16.Schematic of the crack in infinite plane simulation by implementing symmetric boundary condition and traction according to the analytical solution, acquired via conformal mapping of an elliptical hole in infinite plane problem.

Fig. 17 .
Fig. 17.Traction on the edges from analytical solution and the maximum principal stress contour at equilibrium.

Fig. 18 .
Fig. 18. (a)Convergence of strain energy error of the crack in the infinite plane test case (slope at maximum refinements of 1.01 and 0.981 via least-square of all data points) and b convergence of stress intensity factor error of the crack in infinite plane test case; slopes of 1.04 at maximum refinement levels and 1.11 via least-square of all data points are close to the convergence rate of 1 given in[51].

Fig. 19 .
Fig.19.Propagation convergence test case consists of a rectangular domain with initial crack, while tension is applied through the displacement-controlled fixed boundaries at the top and bottom surfaces.

2
Fig. 21(a)  shows the convergence test result for the energy dissipation error, calculated by comparing actual and theoretical dissipation values.The dissipation error shows an overall linear convergence against cell size on a log-log scale, with an asymptotic rate of 0.964.Additionally, Fig.21(b)shows similar load-displacement curves throughout the mesh refinements.The load is calculated from the internal forces at the nodes where displacement-controlled boundaries are applied prior to imposing the Dirichlet boundary condition.

Fig. 20 .
Fig. 20.Displacement-controlled loading scheme of the propagation convergence test case.The loading starts gradually via a quarter cosine function and continued linearly.

Fig. 21 .
Fig. 21.Energy dissipation shows convergence towards lower energy dissipation error with further refinement (a), while load-displacement curve exhibits monotonous trend of lower and earlier peak as cell size decreases (b) (with failure marked when crack propagation reaches the other end of the domain).

Fig. 22 .
Fig. 22. Crack initiation problem consists of a rectangular domain with initial crack emulated by a notch, while tension is applied through the displacement-controlled fixed boundaries at the top and bottom surfaces.The algorithm calculates where the discrete crack should initiate and propagate.

Fig. 23 .
Fig. 23.Maximum principal stress plots of the initiation test case: (a) just before crack initiation stress concentration is visible at the tip of the notch; (b) crack initiated at the blunt tip of the notch due to stress at the nearest node exceeding the material strength; (c) propagation after the initiation of the discrete crack; (d) complete separation of the specimen.

Fig. 24 .
Fig. 24.Crack path development of the initiation problem: (a) propagation after the initiation of the discrete crack; (b) complete separation of the specimen.

Fig. 25 .
Fig. 25.Maximum principal stress plots of the initiation test case with large deformation: (a) stress concentration is visible at the tip of the notch, while large deformation is already visible; (b) crack initiated at the blunt tip of the notch due to stress at the nearest node exceeding the material strength; (c) propagation after the initiation of the discrete crack with large crack opening; (d) complete separation of the specimen.

Fig. 26 .
Fig. 26.Crack path development of the initiation problem with low elastic modulus: (a) A new crack initiates already with a large deformation; (b) propagation occurs with large crack opening; (c) the two sides of the specimen are completely separated.

Fig. 27 .
Fig.27.The merging test case consists of a rectangular domain with a horizontal initial crack at one free edge with tip T 1 and a vertical initial crack in the middle with tips T 2 and T 3 .Tension is applied through the displacement-controlled fixed boundaries at the top and bottom surfaces.

Fig. 28 .
Fig. 28.Crack propagation and merging case with perfectly-vertical middle initial crack: (a) stress concentration at the tip of the horizontal crack; (b) merging as the of the horizontal crack meets the vertical crack; and (c) both tips of the vertical crack (tips 2 and 3) propagate at the same time, which at the end separates the domain completely.

Fig. 29 .
Fig. 29.Crack path development of propagation case with perfectly-vertical middle initial crack: (b) after merging; and (b) complete separation of the specimen.

Fig. 30 .
Fig. 30.Crack propagation sequence and load-displacement curve of the merging case with perfectly-vertical middle initial crack: (a) stress concentration; (b) merging; and (c) complete separation.Propagation starts from tip 1, while both tips 2 and 3 propagates at the same time after merging.

Fig. 31 .
Fig. 31.Crack propagation and merging case with slightly-angled middle initial crack: (a) stress concentration; (b) merging; and (c) complete separation; crack tip 2 propagates first while tip 3 propagates later in the opposite direction.

Fig. 32 .
Fig. 32.Crack path development of propagation case with slightly-angled middle initial crack: (a) after merging; and (b) complete separation of the specimen.

Fig. 33 .Fig. 34 .
Fig. 33.Crack propagation sequence and load-displacement curve of the merging case with slightly-angled middle initial crack.Propagation starts from tip 1, continued by tip 3 after merging, while tip 2 propagates slightly.

Fig. 35 .
Fig. 35.Crack propagation and merging cases after 2x refinement for (a) perfectly-vertical and (b) slightly-angled middle crack; both results are similar to their coarser counterparts.

Fig. 37 .
Fig.37.Displacement-controlled loading scheme of the plate with multiple initial cracks problem.The displacement load applies gradually via a quarter cosine function and continues with a linear function.

Fig. 39 .
Fig. 39.Propagation pattern with the proposed model in MPM: (a) before any propagation where stress concentration starts to appear; (b) small propagations start to occur; (c) propagation advanced further, with some tips already merged; (d) propagation stopped as the four boundaries are already separated.

Fig. 40 .
Fig. 40.Comparison between propagation result of[55] with energy minimization and 1200 × 1200 elements and the proposed method with dynamic MPM and 65536 particles (128 × 128 cells), which show similarity and some differences in propagation pattern.

Fig. 41 .
Fig. 41.MPM simulations of the plate with initial cracks problem with 4096 (a), 16384 (b), and 65536 particles (c).The difference in refinement level affects the produced crack pattern.

Fig. 42 .
Fig. 42.MPM simulations of the plate with initial cracks problem with no damping (a), medium damping (b), and high damping (c).The presence of damping at various levels to minimize dynamic effects significantly influence the final crack pattern.

Fig. A. 1 .
Fig. A.1.Example of crossing profile method with 4 particles, 2 nodes, and 2 crack paths; arrow indicates crack path direction (see[33] on the definition and purpose of crack direction).

Table 1
Mechanical properties of simulation cases.

Table 2
Calculation of various energy sums at complete separation of specimen.

Table A . 1
Grid assignments by the crossing profile method for the example case.