Simulation of Dynamic 3D Crack Propagation within the Material Point Method

Y.J. Guo and J.A. Nairn

1 Introduction

Numerical simulation of dynamic crack propagation is an extremely important research subject not only for academic interest but also for the establishment of a safety design methodology that prevents structures from catastrophic failures. Numerical simulation of dynamic crack propagation, however, remains a challenging problem due to various inherent difficulties. Even though a variety of methods are available for using fracture mechanics in two-dimensional dynamic crack propagation analysis, much less progress has been made on three-dimensional dynamic fracture simulation.

The major numerical techniques currently used in fracture mechanics are the finite element method (FEM) and the boundary element method (BEM). The suitability of those techniques for crack propagation simulation depends on the complexity and versatility of handling crack evolution, computational accuracy and efficiency. For mesh methods, such as FEM, additional techniques are usually needed to simulate crack propagation since the crack is part of the discretized material body. A commonly used technique is the nodal release method [Chen and Wilkins (1977); Kishimoto and Sakata (1987); Ivankovic and Williams (1995)]. The major disadvantage of nodal release is that the crack propagation path is limited to pre-existing mesh lines. It may therefore have severe limitations formodeling mixed-mode crack propagation (where cracks may propagate in unknown directions), and presents even more challenges for modeling three-dimensional propagation. Basically, the node-release method is only reliable for two-dimensional selfsimilar crack propagation. A second frequently used technique is the moving mesh method developed by [Nishioka, Murakami and Takemoto (1990); Nishioka, Okada and Nakatani (1997)], and extended by [Koh, Lee and Haber (1988); Gallego and Dominguez(1992); Koh, Lee and Jeong (1995)]. The method seems difficult for modeling arbitrarily curved crack growth since excessive mesh distortion might occur, and its versatility is still questionable.

FEM and BEM can be adapted to handle arbitrary crack propagation simulations by using remeshing techniques after each increment of crack propagation, which involves implementation of some automatic remeshing scheme [Swanson and Ingraffea (1988);Bittencourt, Wawrzynek and Ingraffea (1996)]. The disadvantages are obvious: remeshing is complicated and takes time. The situation only gets worse when considering arbitrary threedimensional, mixed-mode crack propagation problems.

The extended finite element (XFEM) was developed to avoid some FEM issues by representing the crack independently of the mesh [Belytschko and Black (1999); Mo? s,Dolbow and Belytschko (1999)] and adding enriched shape functions near crack tips and crack planes. The need to remesh is replaced by the need to alter shape functions and to make assumptions about crack tip displacements that make modeling of plastic materials a challenge. Furthermore, XFEM cannot handle interacting cracks, such as two cracks propagating through the same element.

The challenges faced by FEM crack methods might explain the limited progress made in numerical crack propagation, most importantly in three-dimensional crack propagation simulation. The question arises, can some meshless method, such as the smooth particle hydrodynamics (SPH) method [Monaghan (1992)], the element free Galerkin (EFG)method [Belytschko, Lu and Gu (1994)], or the material point method (MPM) [Sulsky,Chen and Schreyer (1994)] provide an alternate to FEM for problems involving crack propagation? In meshless methods, a crack is treated as an independent entity, which may provide greater potential in crack propagation simulations, especially for threedimensional problems.

This paper focuses on the MPM particle-based method for crack modeling with an emphasis on three-dimensional explicit cracks. MPM is a numerical method to solve dynamic solid mechanics problems Sulsky, Chen and Schreyer (1994); Sulsky, Zhou and Schreyer (1995); Sulsky and Schreyer (1996)]. The motivation for MPM development was to simulate problems with history-dependent internal state variables, contact, impact,penetration/perforation, and metal forming without needing master/slave nodes or global remeshing. For modeling cracks, MPM was extended to handle explicit cracks by a method called CRAMP (forCRAcks inMPM) [Nairn (2003)], which was shown to provide excellent results for crack-front fracture parameters such asJintegral [Guo and Nairn (2004, 2006)]. CRAMP can be described as an “extended'” MPM, but it is a more straightforward extension MPM than XFEM is of FEM. It does not require addition of"enriched shaped" functions that depend on crack tip stress state and can handle two interacting cracks in the same background grid cell [Nairn and Aimene (2016)]. As a consequence, MPM may have advantages for some fracture problems. Some examples are problems with crack contact, cracks in inelastic, large-deformation materials, and crack fronts propagating up to and intersecting with or passing through other cracks[Nairn and Aimene (2016)].

In brief, in 2D CRAMP, an explicit crack is described by a series of massless particles connected by line segments, with the first and last particle being the crack tips. Like MPM particles, crack particles are not tied to the background grid. Crack propagation is modeled by adding a new particle ahead of the current crack tip and that propagation can move and propagate in any direction. Each time step in MPM involves extrapolating particle state to the grid, solving the equations of motions on the grid, and then updating the particles. When explicit cracks are present, geometric relations between particles,crack planes, and nodes are used to subdivide the extrapolated grid velocity into two separate velocity fields describing motion “above” and “below” the crack plane. By tracking these two fields, the standard MPM algorithm can be extended to model the crack and implement contact laws on crack surfaces [Nairn (2003), Guo and Nairn(2004)]. In 3D CRAMP, a crack is discretized into a set of triangular elements instead of line segments, but otherwise the 3D algorithm is essential the same as the 2D one [Guo and Nairn (2006)]. A major difference in 3D CRAMP is for crack propagation problems.Rather then dealing only with two crack tips, 3D crack propagation requires analysis of the entire crack front and propagation of that front can be predicted by 3DJintegral methods.

This paper presents the principles and algorithms to simulate arbitrary dynamic crack propagation with the material point method (MPM). MPM has previously been used to model explicit, 2D crack propagation in wood [Nairn (2007)], crack growth with cohesive zones [Nairn (2009); Bardenhagen, Nairn and Lu (2011)] or fiber bridging[Matsumoto and J. A. Nairn (2009)], fracture during cutting [Nairn (2015, 2016)], and geomechanical modeling of hydraulic fracturing [Aimene and Nairn (2014, 2016)]. This paper summarizes CRAMP andJintegral methods as a 3D crack modeling tool. Two 2D examples are discussed and directly compared to FEM and experiments, but the emphasis of this paper is on the extra work needed for simulation of mixed-mode, 3D crack propagation. Two examples show 3D crack propagation in a bar and around the circumference of a hollow tube.

2 Numerical methods

2.1 Stress analysis of cracks in MPM

The numerical methods for stress analysis of cracks with the material point method(MPM) have been developed by the authors of the paper for two-dimensional [Nairn(2003)] and three-dimensional problems [Guo and Nairn (2004)]. Here we give a brief summary as a 3D problem for completeness of description.

In MPM, a material body is discretized into a set of particles, as shown in Figure 1. A particle carries all the information of the material body at that position, such as mass,volume, velocity, strains, and stresses during the simulation. There is no geometrical connection among the particles; they interact through a background grid, which is typically a regular mesh of orthogonal elements. When modeling cracks using CRAMP,the material point description of the object is supplemented by one or more explicit crack planes. In 2D, each crack plane is divided into a set of massless points connected by linesegments [Nairn (2003)] as shown in Fig. 1. In 3D, each crack surface is described by a collection of massless points connected by triangular elements [Guo and Nairn (2006)].

Figure 1: Discretization of a material body and a crack. The solid circles represent material points. The hollow circles are massless particles that define the path of the 2D crack. The background grid is used for solving MPM equations of motion.

During each MPM time step, particle information is extrapolated to the background grid,which is used to solve the momentum equation. In the absence of cracks, particle momenta are extrapolated to a single velocity field resulting in discretized momentum equation for time stepnof:

wherepi(n),fi(n)andfi,T(n)on momentum, internal force, and external traction force on nodei[Sulsky, et al.(1994); Bardenhagen and Kober (2004)]. But such an update enforces a continuous deformation field that cannot describe displacement discontinuities at crack planes. The extension of CRAMP is to introduce two velocity fields to represent motion on opposite sides of the crack plane [Nairn (2003)]. Let each node accommodate one or two velocity fields. When cracks are present, each CRAMP time step starts by assigning a set of crack velocity fields,v(p,i) = 0 or 1, to each particlep– nodeipair. These fields are determined by tracing a line from the particle to the node to see if it crosses any crack plane. Field 0 is when the line crosses no crack (particle and node on same side of the crack) or field 1 is when it crosses one crack (particle and node on opposite sides of the crack). This calculation is best done during the initialization phase of the MPM time step.Note that two velocity fields on each node can handle any number of cracks, but assumes non-interacting cracks such that no node “sees” a line crossing for more one crack. An extension to explicitly handle two cracks by allowing up to four velocity fields per node is give elsewhere [Nairn and Aimene (2016)]. This paper limits discussion to one crack or non-interacting crack problems.

Becausev(p,i) has to look at all possible particle–node pairs, it is the most time consuming task of the CRAMP algorithm. Although it appears to scale asN*n*cwhereNis number of particles,nis number of nodes, andcis number of cracks and therefore to be very expensive for large calculations, the scaling can be made much better. First, by exploiting the local nature of shape functions, only a small number of nodes (which does not increase with problem size) near each particle need to be checked. The CRAMP algorithm thus scales withN*c(or justNfor a constant number of cracks). Second, its efficiency can be greatly improved by screening out the vast-majority of crack-crossing calculations by describing the crack plane using a hierarchical binary tree. The methods are analogous to hidden line removal methods in 3D graphics [Kay and Kajiya (1986)].

Once allv(p,i) are found, the remaining MPM tasks are similar to standard MPM [Sulsky,Chen and Schreyer (1994); Bardenhagen and Kober (2004)] except that they must solve the equations in each crack velocity field, resolve crack contact situations at nodes with more that one crack velocity field [Nairn (2003, 2013)], and update positions of crack surfaces. The revised MPM explicit update for momentum, by solving Eq. (1) for each velocity field, becomes:

where subscriptjindicates velocity fieldjon nodei. The momenta and forces are given by slightly modified MPM extrapolations:

whereSip(n)andGip(n)are the GIMP shape functions and shape function gradients[Bardenhagen and Kober (2004)],mpis particle mass, τp(n)is particle Kirchoff stress, ρ0isundeformed particle density,bpis particle body force,Fp(n)is external force on the particle, andNi(x) is nodal shape functions on the grid, andTpis traction applied to a particle’s surface. These extrapolations are modified from standard MPM extrapolations by the Kronecker δj,v(p,i)that selects the appropriate crack velocity field on each node (i.e.,standard MPM extrapolations omit δj,v(p,i)and drop subscriptj).

Once the grid momenta are updated, the forces and velocities on the grid are used to update the particle velocities (vp) and positions (xp) using:

whereag→p(n)andvg→p(n)are accelerations and velocities extrapolated from the grid to the particle. But, these extrapolations must use the appropriate velocity fields:

wheremi,j(n)is nodal mass extrapolated from the particles:

This particle update can be extended to include various damping schemes as explained in[Nairn (2015)]. The reader is referred to standard MPM papers for more details on these basic equations, without cracks [Sulsky, Chen and Schreyer (1994); Bardenhagen and Kober (2004)] or with cracks [Nairn (2003); Guo and Nairn (2004)]

An additional task need for MPM simulations with cracks is to move the crack plane according to the velocities fields on the two sides of the crack. An improvement over the initial development of CRAMP, which moved crack particles in the center-of-mass velocity field [Nairn (2003)], is for each crack particle to track its’ position as well as locations for the crack surfaces above and below the crack plane. The locations on top and bottom crack surfaces are moved in velocity field associated with that surface. In other words, the crack surfaces move by

Here subscript “s” indicates side of the crack (0 or 1) being updated andvg→s(n)is velocity extrapolated from the grid to the location for that surface of the crack. It is similar to the velocity extrapolation for particle updates (see Eq. (9)) except that is uses velocity fieldv(s,i) appropriate for the surface locations-nodeipair, which is determined by same methods used to determinev(p,i).

Note that the GIMP shape function,Sip(n), in Eq. (9) is replaced by a grid shape functionNi*(xs(n)). This replacement occurs because GIMP shape functions integrate over particle domain. Because crack particles are massless with no domain, they should use GIMP shape functions with a Dirac delta function for domain, which converts them to ordinary grid shape functions [Bardenhagen and Kober (2004)]. The “*” on grid shape function arises because the extrapolation may need renormalization to preserve partition of unity caused by potentially inactive nodes near the crack surface (i.e.,Ni*(xs(n)) isNi(xs(n)) divided by Σi Ni(xs(n))). The crack particle location can update by center-of-mass velocity [Nairn (2003)],or, more efficiently, by moving to the midpoint between the two crack surfaces. Some advantages of tracking locations of crack surfaces are that it allows for improved contact modeling [Nairn (2013)] and implementation of imperfect interface [Nairn (2007)] or cohesive laws [Nairn (2009)] that depend on crack opening displacement. These displacements are also used to partitionJintegral into mode I and II stress intensity factors,as explained in the next section.

2.2 Dynamic j-integral and stress intensity factors

The CRAMP algorithm in the previous section explains how MPM can account for explicit discontinuities caused by cracks. To handle crack propagation, the next task is to calculate crack front fracture parameters that are needed to predict when and in which direction a crack will propagate. The general method is to use the crack frontJ-integral.MPM methods for findingJintegral and stress intensity factors in 2D and 3D simulations are in [Guo and Nairn (2004, 2006)]. The two results are summarized here together for completeness.

The dynamicJ-integral components in crack-front coordinates was formulated by[Nishioka and Atluri (1983) and Nishioka (1995)], as follows:

wherek= 1 or 2 for the two components.WandKdenote the stress-work density and kinetic energy density, respectively; σijis a component of Cauchy stress,uia component of displacement (accordingly, ?ui/?xkis an element of displacement gradient),nkis a component of the unit normal vector to the near-fieldJ-integral contour Γc, ρ is density,and repeated indicesiandjare summed. For 2D problems, the Γccontour is most conveniently drawn along background grid mesh lines near the crack tip [Guo and Nairn(2004)]. For 3D problems, theJ-integral is calculated at each crack-front node with near-field integral contour (Γc) being a circular section of a short cylinder on anx1-x2plane in crack-tip local coordinates with the origin located at the crack front node and a radius ofr,as shown in Figure 2. Herex2coincides with the direction of the normal to the crack plane at the node,x3(which is the thickness direction of the cylinder) is tangent to the crack front, andx1points in the direction of bi-normal at that position [Guo and Nairn(2006)]. The second term is an integral over the volume contained by Γc. It is zero for quasi-static problems but is needed to provide path-independent results in dynamic problems [Nishioka (1995)]. In MPM, the second term can be found by using the material points within the contour as numerical integration points [Guo and Nairn (2004)].

Figure 2: Dynamic J-integral calculation scheme for 3D crack problems: (a) crack plane discretization; (b) J-integral contour (Γc), which is a circle on x1-x2 plane with the origin located at crack front node with a radius of r.

Once the dynamicJ-integral is calculated, it can be converted into mode I and mode II crack tip stress intensity factors-KIandKII. (these first calculations were for problems with little or no mode III and therefore ignoresKIII). For linear elastic, isotropic materials,the formulae are given by [Nishioka (1995)]:

whereJ1is the first component of dynamicJ-integral vector in crack-front coordinates. It equals the total strain energy release rate for elastic materials.μis the shear modulus and δIand δIIdenote crack opening and shearing displacements near the crack tip. βI, βII,AI,andAIIare parameters related to crack propagating velocityC. They are given by:

whereCsandCdare the shear and dilatational wave speeds:

where κ = (3-ν)/(1+ν) for 3D or 2D plane stress, κ = 3-4ν for 2D plane strain, and ν is Poisson’s ratio. For slow crack growth, the calculations can useC= 0, βI= βII= 1 andAI=AII= (κ+1)/4. This slow limit applies well unless the crack speeds get very close to wave speed of the material.

2.3 Modeling of crack propagation

2.3.1 Crack propagation criteria

Because the effect of the tearing mode stress intensity (KIII) on crack propagation direction in general 3D mixed-mode loading is not well understood, we neglect the effect of mode III stress intensity, and assume that only modes I and II affect the direction of crack propagation (it could be added in the future, when needed). It is common to find thatKIIIis small compared toKIorKIIin many real-world problems, and examples chosen here all have negligible mode III component. As a result, a crack-front node will propagate in thex1-x2plane in the crack-front local coordinates at that node. In order to perform crack propagation analysis, we need to predict on what conditions the crack will propagate, and the propagation direction for the crack under mixed-mode loading, where mixed-mode ratio is α =KI/KII. Several criteria have been proposed, and the following two are widely used: (1) the maximum hoop (or principal) stress criterion proposed by[Erdogan and Sih (1963)]; (2) the minimum strain energy density criterion presented by[Sih (1972, 1974)].

The maximum hoop stress criterion [Erdogan and Sih (1963)] states that the crack will propagate in the direction normal to the maximum principal stress when the equivalent stress intensity factor reaches a critical value,KIc,i.e., when

where θcis angle for the direction with the maximum principal stress. It is given by:

where we take the positive sign forKII≤ 0, and the negative sign otherwise.

The minimum strain energy criterion [Sih (1972, 1974)] assumes that the crack propagates in the direction along which the strain energy density is minimum when the equivalent stress intensity reaches the critical valueKIc,i.e., when

where θcis the direction along which the stain energy density is minimum. It is found by numerically solving:

Note that the two propagation direction criteria (Eqs. (20) and (24)) give nearly identical functions for θcas a function of α, which makes it impossible to prefer one over the other by comparison to experimental observations of propagation direction with typical uncertainties in measuring θc.

2.3.2 Crack evolution algorithms

In dynamic fracture analysis within MPM, we calculate dynamic stress intensity factors(KIandKII) through dynamicJ-integral at crack-front nodes in crack-front local coordinates. These values are compared to a material critical value (using Eq. (19) or (21))to predict whether or not the crack should propagate. If it propagates, Eq. (20) or (24) are used to determine the direction of propagation. Propagation in 2D involves adding a new crack particle in the calculated direction at a distance of less than a background cell (we typically extend it half a grid cell length). The new particle becomes the new crack tip.The situation in 3D is more complex and described next.

After propagating 3D crack-front nodes, we need to build additional crack elements to represent the new crack growth. Figure 3 illustrates the scheme used for 3D MPM crack propagation analysis, where A through G are assumed to be the crack-front nodes before crack propagation, and B' through G' are the new positions of the nodes that result from propagation. If only one end of a crack segment propagates (such as segment AB), a single new crack surface element will be created (such as ABB'). If both ends of a segment propagate (such as segment BC), then four new crack surface elements are created (such as the four elements within BCC'B'). If a new crack-front segment becomes too long after propagation, say two times longer than the initial length, the segment will be broken into two, and five new crack elements will be created (such as five elements within CDD'C' with added point P1along the long segment in the new crack front). If a crack-front segment becomes too short after propagation, say less than 25% of the initial length, the segment will be shrunk into a point, and only one new crack element will be created (such as E'F' combining to P2). This adaptive crack-front discretization procedure is needed to simulate curved crack fronts and to assure computational accuracy.

Figure 3: Illustration of building new crack plane elements after crack propagation

2.3.3 Interaction between crack plane and material boundary

Special treatment is needed when crack propagation interacts with the boundary of an object. For example, Fig. 4 shows a crack front (solid blue lines) propagating within a cross section of a material where the dashed lines represent the real material boundaries.After increments of propagation (the solid blue lines are the crack fronts at increments of time that move from bottom to top and toward both sides). In 2D crack growth, when a crack tip reaches the boundary its’ propagation stops. But, in 3D, only part of the front may reach a boundary while the rest should continue to propagate. 3D crack propagation methods have to deal with these boundary interactions.

In MPM, a material body is discretized into a set of particles. Although this discretization does not explicitly track locations of boundaries, we can determine them approximately along the grid lines with the aid of the gradient of extrapolated nodal volume. We call these approximate boundaries the MPM material boundaries, resulting in boundaries such as the dotted red lines in Fig. 4. At each time step, if a newly generated crack-front reaches the MPM material boundaries, it will get trimmed at the MPM material boundaries, as shown in Figure 4. If both ends of a crack-front segment are at or outside the MPM material boundary, that segment is beyond the material, such as segments AB and BC. Such segments will be eliminated as a crack-front, which means they can no longer propagate. It is the 3D analogy of 2D model that stops all propagation when a crack tip reaches the edge, but now only a portion of the crack front stops propagating.

Figure 4: Illustration showing interaction between crack front and material boundaries during crack propagation, where dashed lines are real material boundaries, dotted lines are MPM material boundaries on the grid, and solid lines are crack fronts showing increments in crack propagation.

3 Numerical crack propagation examples

Prior papers have used MPM for 2D crack propagation [e.g., Nairn (2007, 2015, 2016);Bardenhagen, Nairn and Lu (2011); Aimene and Nairn (2014, 2016)]. This section supplements 2D results by direct comparison to FEM and experiments. It also adds, for the first time, two examples of 3D crack propagation using MPM.

3.1 Crack propagation length versus time

A three-point bending specimen subjected to central impact was modeled (see Fig. 5). This problem was experimentally investigated, and crack length versus time was measured by[Nishioka, Tokudom eand Kinoshita (2001)]. The specimen was made of PMMA, and the following material properties were used in the analysis: modulusE= 2.94 GPa, Poisson’s ratio ν = 0.3, density ρ = 1190 kg/m3, and fracture toughnessKIc= 1.2 MPa√m (adjusted to match experiments). The dimensions of the specimen were as follows: span 2l= 400 mm,heighth= 100 mm, thickness = 10 mm, and the initial crack lengtha0= 50 mm. A stiff, 5.05 kg drop weight (modeled withEw= 100 GPa and νw= 0.25) impacted the specimen’s middle with the velocity of 5 m/s. The geometry of the drop weight was estimated from images in[Nishioka, Tokudom eand Kinoshita (2001)] with its density adjusted to give proper final weight (the full drop weight details were not available). Because the drop weight mass was 10 times the beam mass, it moved at nearly constant velocity (<10% reduction during the impact event). The 2D calculations used plane strain conditions. TheJintegral calculations used both terms in Eq. (13). As expected, due to specimen symmetry, the crack propagated along the initial direction as self-similar crack growth.

Figure 5: A three-point bending specimen subjected to central impact.

The results for crack length versus time predicted by MPM using the minimum strain energy criterion with various cell sizes (c= 1, 2, and 4 mm) are plotted in Figure 6 and compared to experimentally measured results. The results were nearly identical (and hence converged) for all cell sizes. Thus, not only is remeshing not needed, but also MPM crack calculations converge well. The largest cell size (c = 4 mm), however,showed discrete jumps caused by each increment in crack length being half a cell length or 2 mm. The results for cell sizes of 1 and 2 mm were smoother. The calculations were very fast (from 2 sec for 4 mm cells to 90 sec for 1 mm cells on a desktop computer with a 2.7 GHz 12-Core Intel Xeon E5 processor). The MPM simulations and experiments agreed well. Perhaps the MPM results had slightly faster crack growth at later times,which could be due to dynamic changes in stress intensity factor. According to Eq. (14),KIdecreases at higher crack growth rates, which could cause propagation to slow, but these simulations did not implement that adjustment. Because the maximum simulated(and experimental) crack velocity was about 450 m/s andCd= 1240 m/s for PMMA, this refinement would only changeKIby 7% at most. The slight discrepancies are more likely caused by details about the impactor or the contact conditions on bottom edge of the beam; these effects could be modeled with more information on the experimental arrangement. The corresponding results using the maximum hoop stress criterion were identical (and therefore not shown).

Figure 6: Comparison of the results of crack length versus time predicted by MPM with various mesh sizes (cell size c=1, 2, and 4 mm) to experimental results.

3.2 Two-dimensional crack propagation path

Figure 7: A three-point bending specimen subjected to impact with the eccentricity of e/l=0.1.

The same three-point bending specimen was modeled to simulate mixed-mode crack propagation path caused by the drop weight impacting off the center of the specimen, as shown in Figure 7. The dimensions and material properties of the specimen were the same in section 3.1. The loading eccentricity,e/l, was 0.1. Due to asymmetry of the specimen, the crack tip is loaded under mixed-mode conditions, which is expected to lead to a change in crack direction during crack propagation. This specimen was studied experimentally and also simulated with the finite element method (FEM) by [Nishioka,Tokudom eand Kinoshita (2001)].Crack paths predicted by MPM simulations with 1 mm cells using either the maximum hoop stress (“mh”) or minimum strain energy (“mse”) criterion are shown in Fig. 8 where they are compared to experimental results and two alternate FEM predictions [Nishioka,Tokudom eand Kinoshita (2001)]. The FEM predictions used either the maximum hoop stress criterion (“mh”) or a local symmetry criterion based on the direction whereKII= 0.The experimental and numerical results are for total crack growth att= 210 μs after the impact. All methods, including both MPM and FEM simulations, agreed reasonably well with experimental results; the crack paths are within 2 mm of experimental results. The maximum hoop stress criterion worked best for these experiments, with MPM slightly improved over FEM (especially in the beginning). TheKII= 0 criterion by FEM had the worst agreement. The MPM results were fast (about 80 sec on desktop computer with a 2.7 GHz 12-Core Intel Xeon E5 processor) and did not need to use remeshing that was required for FEM results.

Figure 8: Comparison of crack growth paths simulated by the MPM and FEM using maximum hoop stress (“mh”), minimum strain energy (“mse”), and local symmetry condition (“KII=0”) criteria to the experimentally measured path. The experiments and simulations were for t = 210 μs after the impact. The inset shows experimental snap shot from [Nishioka, Tokudom eand Kinoshita (2001)]. The black circle is a caustic used to monitor crack tip stress state.

3.3 Three-dimensional dynamic crack propagation

Three-dimensional (3D) dynamic crack propagation simulation is essentially an unexplored area, with very few reports available [Nishioka and Stan (2003); Krysl and Belytschko (1999)]. The next two examples are used to demonstrate the capabilities of MPM to simulate full 3D, mixed-mode crack propagation. As mentioned previously, a crack plane in MPM is treated as an independent entity, not connected geometrically to the material body or background grid. Therefore, the initial mesh can be used throughout the whole process of crack propagation, and no remeshing is needed, which makes it possible to model full 3D, mixed-mode crack propagation in arbitrary directions.

The first 3D example simulated crack propagation process far an inclined corner crack (a quarter penny-shaped crack) in a square rod subjected to step loading (σ = 400 MPa), see the “Initial” condition in Fig. 9. The length of the specimen was 100 mm, and the cross section of the specimen was 50×50 mm. The angle between the load direction and the initial crack plane was 60o. The origin of the corner crack was located on a generator of the specimen, and the center of the crack plane crossed the middle cross-section of the specimen. The radius of the initial corner crack was 18.475 mm. The material properties used in the analysis were: modulusE= 200 GPa, Poisson’s ratio ν = 0.298, density ρ =7900 kg/m3, and fracture toughnessKIc= 20 MPa√m. TheJintegral calculations used only the first term in Eq. (13). This change was done for efficiency (i.e., to avoid searching for material points within the contour) when needing all the extraJintegral calculations compared to 2D problems. Previous results have shown that the second term is small when the contour is within a few cells of the crack front [Guo and Nairn (2004)].Crack patterns at different time instants simulated by MPM are shown in Figure 9 with the last frame being when the entire crack front has reached the specimen edge. It can be seen that the crack pattern changed suddenly from mixed-mode to mode I dominated type after it starts to propagate. The crack front contour length grew in length during the first stage of propagation, but then got shorter and shorter after it reached the diagonal of the cross section of the specimen. The proposed method for dealing with crack front propagation and interactions with material boundaries worked well throughout the propagation process. In particular, the method correctly interacted with edges and was able to propagate through to complete failure where the entire front had reached a specimen edge.

Figure 9: Crack patterns at different time increments simulated by MPM for an inclined corner crack in a square rod subjected to step tension.

The second 3D example was a thin, hollow tube with an inclined, through-the-wall crack subjected to step tension (σ = 400 MPa), as shown in the initial condition of Figure 10.The length of the tube was 300 mm and outer the radius of the tube wasR= 45 mm. The wall thickness wast= 10 mm. The angle between the crack plane and the loading direction was 450. The through-the-wall crack plane was centrally located along the length direction of the tube with a crack length of 14.142 mm. The material properties used in the analysis were the same as those for the square rod.

Figure 10: Crack patterns at different time increments simulated by MPM for an inclined,through-the-wall crack in a hollow tube subjected to step tension.

The initial crack and subsequent crack propagation patterns at different time increments are shown in Fig. 10. Both ends of the inclined crack turned quickly from the mixedmode to mode I propagation and then continued to propagate around the entire circumference of the tube. The two cracks did not meet because they started from offset crack tips. This example demonstrated the robust capability of this 3D MPM crack propagation method to interact with materials boundaries while still propagating within the material.

The above two examples of full 3D, mixed-mode dynamic crack propagation were not compared to experiments because none are available [Nishioka and Stan (2003)]. One challenge of experiments is how to map crack plane surfaces within opaque materials.The modeling results, however, are giving expected results (crack turning to mode I conditions) and robustly modeling crack propagation through to complete failure and separation (had the simulations continued) along with stable interactions with edges.

4 Conclusions

This paper presented the principles and algorithms of dynamic crack propagation simulations within the material point method (MPM). The two-dimensional crack propagation summarized and expanded on some previous results and included new comparisons to experiments and FEM. The full three-dimensional, mixed-mode crack propagation was described and demonstrated for the first time. MPM exhibits the flexibility of meshless methods in crack propagation simulation due to the fact that the crack plane is treated as an independent geometrical entity. The crack is not restrained by the background grid and is able to follow any arbitrary path dictated by the crack tip stress state. Furthermore, the initial grid can be used repeatedly throughout the whole process of crack propagation, and no remeshing is needed. These features are advantages of MPM compared to the typical and widely used mesh methods such as the finite element method (FEM) or the boundary element method (BEM). This paper also described an adaptive crack plane scheme for crack-front division, and the algorithms to interact and trim the crack plane at the material boundaries. The applications of the methodology to several typical problems have shown that MPM is a reliable and powerful approach to simulate dynamic crack propagation. The convergence in terms of background cell size is efficient. The methods were shown to work from relative simple two-dimensional problems to very complicated three-dimensional, mixed-mode problems.

Acknowledgement:This work, some of which was done several years ago, was supported by the University of Utah Center for Simulation of Accidental Fires and Explosions (C-SAFE), funded by the Department of Energy, Lawrence Livermore National Laboratory, under Subcontract B341493. We also thank J. Davison de St.Germain for help in 3D graphics of crack propagation.

Aimene, Y. E.; Nairn, J. A.(2014): Modeling multiple hydraulic fractures interacting with natural fractures using the material point method. InSPE/EAGE European Unconventional Conference and Exhibition, Vienna, Austria, Feb. pp.25-27.

Bardenhagen, S. G.; Kober, E. M.(2004): The generalized interpolation material point method.Computer Modeling in Engineering & Sciences, vol. 5, pp. 477-496.

Bardenhagen, S. G.; Nairn, J. A.; Lu, H.(2011): Simulation of dynamic fracture with the material point method using a mixed J-integral and cohesive law approach.Int. J.Fracture, vol. 170, pp. 49-66.

Bathe, K. J.(1996):Finite Element Procedures. Prentice-Hall, Inc., pp.768.

Belytschko, T.; Lu, Y. Y.; Gu, L.(1994): Element Free Galerkon Methods.International Journal for Numerical Methods in Engineering, vol. 37, pp. 229-256.

Belytschko, T.; Black, T.(1999): Elastic crack growth in finite elements with minimal remeshing.Int. J. for Numerical Methods in Engineering, vol. 45, pp. 601-620.

Bittencourt, T. N.; Wawrzynek, P. A.; Ingraffea, A. R.(1996): Quasi-Automatic Simulation of Crack Propagation for 2-D LEFM problems. Engineering Fracture Mechanics, vol. 55,pp.321-334.

Chen, Y. M.; Wilkins, M. L.(1977): Numerical Analysis of Dynamic Crack Problems. In:Elastodynamic Crack problems, Edited by G.C. Sih, Noordhoff, Leiden.

Erdogan, F.; Sih, G. C. (1963): On the Crack Extension in Plates under Plane Loading and Transverse shear.Journal of Basic Engineering, ASME, vol. 85, pp. 519-527.

Gallego, R.; Dominnguez, J.(1992): Dynamic Crack Propagation Analysis by Moving Singular Boundary Elements.J. Appl. Mech. Trans., ASME, vol. 59, pp. S158-S162.

Guo, Y. J.; Nairn, J. A.(2004): Calculation of J-Integral and Stress Intensity Factors Using the Material Point Method.Computer Modeling in Engineering & Science, in press.

Guo, Y. J.; Nairn, J. A.(2006): Three-Dimensional Dynamic Fracture Analysis Using the Material Point Method. Computer Modeling in Engineering & Science, vol. 16, pp.141-156.

Ivankovic, A.; Williams, J. G.(1995): The Finite Element Analysis of Linear Elastic Dynamic Fracture. In:Dynamic Fracture Mechanics, edited by M.H. Aliabadi, Southampton,UK, and Boston, USA,Computational Mechanics Publications, pp.101-136.

Kay, T. L.; Kajiya, J. T.(1986): Ray tracing complex scenes.ACM SIGGRAPH, vol. 20,no. 4, pp. 269.

Kishimoto, S.; Sakata, M.(1987): Finite Element Computation of Dynamic Stress Intensity Factor for A rapidly Propagating Crack Using J-Integral.Computational Mechanics, vol. 2, pp. 54-62.

Koh, H. M.; Lee, H. S.; Haber, R. B.(1988): Dynamic Crack Propagation Analysis Using Eulerian-Lagrangian Kinematic Description.Computational Mechanics, vol. 3, pp.141-155.

Koh, H. M.; Lee, H. S.; Jeong, U. Y.(1995): An Incremental Formulation of the Moving Grid Finite Element Method for the Prediction of Dynamic Crack Propagation.Nucl. Engng.Des., vol. 158, pp. 295-309.

Krysl, P.; Belytschko, T.(1999): The Element Free Galerkin Method for Dynamic Propagation of Arbitrary 3-D Cracks.International Journal for Numerical Methods in Engineering, vol. 44, pp. 767-800.

Matsumoto, N.; Nairn, J. A.(2009): The fracture toughness of medium density fiberboard (MDF) including the effects of fiber bridging and crack-plane interference.Eng. Fract. Mech., vol. 78, pp. 2748-2757.

Mo? s, N.; Dolbow, J.; Belytschko, T.(1999): A finite element method for crack growth without remeshing.Int. J. for Numerical Methods in Engineering, vol. 46, pp. 131-150.

Monaghan, J. J.(1992): Smooth Particle Hydrodynamics.Annual Review of Astronomy and Astrophysics, vol. 30, pp. 543-574.

Nairn, J. A.(2003): Material Point Method Calculations with Explicit Cracks.Computer Modeling in Engineering & Science, vol. 4, pp. 649-664.

Nairn, J. A.(2007): Material point method simulations of transverse fracture in wood with realistic morphologies.Holzforschung, vol. 61, pp. 375-281.

Nairn, J. A.(2009): Analytical and numerical modeling of R curves for cracks with bridging zones.Int. J. Fract., vol. 155, pp. 167-181.

Nairn, J. A.(2013): Modeling of imperfect interfaces in the material point method using multimaterial methods.Computer Modeling in Eng. & Sci., vol. 92, no. 3, pp. 271-299.

Nairn, J. A.(2015): Numerical simulation of orthogonal cutting using the material point method.Engineering Fracture Mechanics, vol. 149, pp. 262-275.

Nairn, J. A.(2016): Numerical modeling of orthogonal cutting: Application to woodworking with a bench plane.Interface Focus, vol. 6, no. 3, p. 20150110.

Nairn, J. A.; Aimene, Y.(2016): Modeling interacting and propagating cracks in the material point method: Application to simulation of hydraulic fractures interacting with natural fractures.in preparation.

Nishioka, T.(1995): Recent Developments in Computational Dynamic Fracture Mechanics.In:Dynamic Fracture Mechanics, edited by M.H. Aliabadi.Computational Mechanics Publications, Southampton, UK, pp.1-58.

Nishioka, T.; Atluri, S. N.(1983): “Path-Independent Integrals, Energy Release Rates,and General Solutions of Near-Tip Fields in Mixed-Mode Dynamic Fracture Mechanics.Engineering Fracture Mechanics, vol.18, no.1, pp.1-22.

Nishioka, T.; Murakami, R.; Takemoto, Y.(1990): The Use of Dynamic J-Integral (J’)in Finite Element Simulation of Mode I and Mixed-Mode Dynamic Crack Propagation.International Journal of Pressure Vessels and Piping, vol. 44, pp. 329-352.

Nishioka, T.; Okada, T.; Nakatani, M.(1997): Mixed-Phase Simulation with Fracture-Path Prediction Mode for Dynamically Curving Fracture.Advanced in Fracture Research,Pergamon Press, New York, vol. 4, pp. 2063-2070.

Nishioka, T.; Stan, F.(2003): A Hybrid Experimental-Numerical Study on the Mechanism of Three-Dimensional Dynamic Fracture.Computer Modeling in Engineering & Science, vol.4, pp.119-139.

Nishioka, T.; Tokudome, H.; Kinoshita, M.(2001): Dynamic Fracture-Path Prediction in Impact Fracture Phenomena Using Moving Finite Element Method Based on Delaunay Automatic Mesh Generation.International Journal of Solids and Structures, vol. 38, pp.5273-5301.

Sih, G. C.(1972): Introductory Chapter: A Special Theory of Crack Propagation.Mechanics of Fracture, Noordhoof, Hooland, vol. 1, pp. XXI-XLV.

Sih, G. C.(1974): Strain-Energy-Density Factor Applied to Mixed Mode Crack Problems.International Journal of Fracture, vol. 10, pp. 305-321.

Sulsky, D.; Chen, Z.; Schreyer, H. L.(1994): A Particle Method for History-Dependent Materials.Computer Methods in Applied Mechanics and Engineering, vol. 118, pp. 179-196.

Sulsky, D.; Schreyer, H. L.(1996): Axisymmetric Form of the Material Point Method with Applications to Upsetting and Taylor Impact Problems.Computer Methods in Applied Mechanics and Engineering, vol. 139, pp. 409-429;

Sulsky, D.; Zhou, S. J.; Schreyer, H. L.(1995): Application of A Particle-in-Cell Method to Solid Mechanics.Computer Physics Communications, vol. 87, pp. 236-252.

Swanson, D. V.; Ingraffea, A. R.(1988): Modeling Mixed-Mode Dynamic Crack Propagation Using Finite Elements: Theory and Applications.Computational Mechanics,vol. 3, pp. 381-397.