Phase‐field material point method for brittle fracture

INTERNATIONAL JOURNAL FOR NUMERICAL METHODS IN ENGINEERING Int. J. Numer. Meth. Engng 2016; 00:1–47 Published online in Wiley Online Library (wileyonlinelibrary.com). DOI: 10.1002/nme Phase-Field Material Point Method for Brittle Fracture E. G. Kakouris ∗ and S.P. Triantafyllou ∗† Centre for Structural Engineering and Informatics, The University of Nottingham, Coates Building, Nottingham, NG7 2RD, UK SUMMARY The Material Point Method for the analysis of deformable bodies is revisited and originally upgraded to simulate crack propagation in brittle media. In this setting, phase field modelling is introduced to resolve the crack path geometry. Following a particle in cell approach, the coupled continuum/ phase-field governing equations are defined at a set of material points and interpolated at the nodal points of an Eulerian, i.e. non-evolving, mesh. The accuracy of the simulated crack path is thus de-coupled from the quality of the underlying finite element mesh and relieved from corresponding mesh-distortion errors. A staggered incremental procedure is implemented for the solution of the discrete coupled governing equations of the phase field brittle fracture problem. The proposed method is verified through a series of benchmark tests while comparisons are made between the proposed scheme, the corresponding finite element implementation as well as experimental results. Copyright c 2016 John Wiley & Sons, Ltd. Received . . . KEY WORDS: Material Point Method; Phase-Field Model; Fracture Mechanics ∗ Correspondence to: S.P. Triantafyllou, Centre for Structural Engineering and Informatics, The University of Nottingham, Coates Building, Nottingham, NG7 2RD, UK. ∗ E-mail: evxek3@nottingham.ac.uk † E-mail: savvas.triantafyllou@nottingham.ac.uk Copyright c 2016 John Wiley & Sons, Ltd. Prepared using nmeauth.cls [Version: 2016/03/02 v3.01] 2 E. G. KAKOURIS, S. P. TRIANTAFYLLOU 1. INTRODUCTION Simulation of damage pertinent to crack initiation and crack growth is an intriguing and challenging aspect of Computational Mechanics. Damage modelling has received considerable attention during the past thirty years as it is relevant to a number of natural and industrial processes, e.g., composite material behaviour [1], concrete fracture [2] and ice mechanics [3] amongst many. Within this setting, damage is being treated either within a continuum (or smeared) phenomenological framework [4] or through discrete methods where the geometry of the crack is explicitly approximated, see, e.g., [5]. Thus, discrete methods can provide a better insight on the actual cracked configuration of a deformable body and form the basis for the study of related phenomena, e.g., corrosion [6]. Initial efforts in discrete crack approaches include element deletion method and re-meshing strategies whereas more sophisticated techniques involve the eXtended Finite Element Method (XFEM) [7, 8, 9], cohesive element methods [10, 11] and cohesive segments methods [12]. Cohesive element methods also based on the notion of configurational force [13] are being used to address the crack initiation and propagation problem with the accuracy of the solution depending on the quality of the underlying finite element mesh. In these methods, the evolution of complex crack paths, including merging cracks, needs to be tracked algorithmically. This increases the complexity of the underlying computational scheme and also the required computational resources. Variational methods for fracture emerged in an effort to address such computational issues. Within this set of methods, Bourdin et al. [14] utilized the mathematical framework of phase-field theory [15] to provide a consistent theoretical framework of the analysis of crack propagation problems. Phase field models represent cracks by means of an additional continuous field (phase field) that smoothly varies from zero (on the crack) to one (away from the crack) (see, e.g., [16]). The evolution of the additional field is defined on the basis of additional governing equations pertaining to the mathematics of phase-field theory [17] linked however to a phenomenological framework such as Griffith’s theory for brittle fracture [18]. The phase field evolution equations are weakly coupled Copyright c 2016 John Wiley & Sons, Ltd. Prepared using nmeauth.cls Int. J. Numer. Meth. Engng (2016) DOI: 10.1002/nme PHASE-FIELD MATERIAL POINT METHOD FOR BRITTLE FRACTURE 3 to the standard governing equations (i.e., equilibrium, compatibility and constitutive equations) of the continuum, effectively introducing a coupled-field problem. This is solved using any standard discretization procedure such as the Galerkin method [19]. In this setting, the crack path emerges from the direct solution of the coupled field framework. This renders phase field methods a promising computational tool to tackle fracture mechanics, at the cost however of introducing additional unknowns, i.e., the phase field. Phase field modelling has been successfully applied within grid based methods i.e. the finite element method (see, e.g., [20, 21, 22]) and its isogeometric variant [23] for the case of quasi-static fracture. Furthermore, Phase field fracture modelling has been effectively applied to treat dynamic fracture propagation problems [16, 24]. However, treating crack propagation using a grid based method introduces further challenges as robustness and accuracy directly depend on mesh quality and corresponding mesh distortion errors. Avoiding numerical errors due to mesh distortion is not a trivial task in grid based Lagrangian methods (see, e.g., [13, 25]). Failure to bound such mesh-dependent errors may result in considerable loss of accuracy especially if large displacements and/or large deformations are taken into account. Discrete element methods [26], smooth particle hydrodynamics [27] and peridynamics [28] can also efficiently deal with problems of fracture mechanics where large deformation take place. Recently, phase field modelling has been introduced within the context of local maximumentropy meshfree approximants to address the problem of fracture in thin-shells [29]. Although robust, especially when dealing with complex geometrical domains, purely meshless methods are computationally taxing as a set of additional procedures is required to achieve convergence, i.e., high-order integration schemes and neighbour searching [30]. To mitigate such issues, Material Point Method (MPM) [31] has been introduced as an extension of Particle-In-Cell methods that efficiently treats history-dependent variables. Combining concepts pertinent to both the Eulerian and Lagrangian description of classical mechanics [32], MPM has been proven advantageous in the analysis of large scale problems involving material and geometric non-linearities, especially within Copyright c 2016 John Wiley & Sons, Ltd. Prepared using nmeauth.cls Int. J. Numer. Meth. Engng (2016) DOI: 10.1002/nme 4 E. G. KAKOURIS, S. P. TRIANTAFYLLOU the context of coupled field problems, e.g., fluid-structure interaction [33] and poro-mechanics [34], also within a large deformation hydrodynamic setting [35, 36]. Very recently, the method has been further extended to account for axisymmetric problems further generalizing its remit. In MPM, the continuum is represented by a set of material points that are moving within a fixed (Eulerian) computational grid where solution of the governing equations is performed. The grid is used to evaluate the gradient and divergence terms of each material point. MPM has been found to offer significant computational advantages when compared to purely meshless methods since it does not require time-consuming neighbour searching. With regards to fracture, the fact that material behaviour is monitored at material points that move within a fixed Eulerian grid implies that the transition from continuous to discontinuous displacement field can be modelled without the need for remeshing the computational grid and without the requirement to account for and mitigate mesh-distortion due to crack propagation. Despite this, few research has been conducted to model the problem of damage modelling and in particular crack growth utilizing the MPM. In [37] and also [38] decohesion was treated by introducing a cohesive material constitutive framework at the material point level. Brittle fracture within an MPM setting was examined for the first time in [39] although considering only the case of pre-existing, i.e., explicit, crack geometries by allowing multiple velocity fields to be defined at the background grid. More recently, cohesive modelling approaches have been introduced in an effort to further generalize the applicability of the MPM for problems pertinent to arbitrary crack paths [40, 41, 42, 43]. A continuum damage based approach has been introduced in [44] also demonstrating the advantages of utilizing domain decomposition methods to accelerate MPM. The aforementioned approaches demonstrated the merits of MPM in simulating damage in terms of computational simplicity in particular when considering the case of large deformations and contact-fracture related problems. Further to the current state-of-the-art, phase field modelling for brittle fracture is introduced in this work within a material point method setting to address the general problem of quasi-static crack propagation in brittle materials using the MPM. By Copyright c 2016 John Wiley & Sons, Ltd. Prepared using nmeauth.cls Int. J. Numer. Meth. Engng (2016) DOI: 10.1002/nme 5 PHASE-FIELD MATERIAL POINT METHOD FOR BRITTLE FRACTURE introducing phase fields at the material point level rather than the nodal points of a fixed Lagrangian grid, the proposed method succeeds in monitoring crack initiation and growth in an efficient and robust manner. Numerical investigations demonstrate that compared to the standard phase-field finite element implementation, the proposed method is advantageous in terms of accuracy. A staggered strategy is utilized for the solution of the governing coupled equations of the problem. An implicit rather than explicit phase-field MPM implementation is established due to the quasistatic nature of the problem under investigation. However, the proposed coupled scheme renders itself naturally to an explicit solution approach. This paper is organized as follows. In Section 2 phase field modelling for brittle fracture is briefly described to facilitate presentation of subsequent derivations. MPM implementation for brittle fracture which constitutes the core contribution of this work is presented in Section 3. The numerical procedure implemented for the solution of the resulting coupled field system of equation is described in Section 4. In Section 5 a set of benchmark problems are examined to verify the proposed formulation compared to the standard phase field Finite Element implementation. Validation of the method is also performed with respect to published experimental results. 2. PRELIMINARIES 2.1. Brittle Fracture The purpose of this work is to introduce a phase-field approximation for brittle fracture within MPM. Thus, derivations presented herein pertain to Griffith’s theory for brittle fracture [18], although generalization to the case of ductile fracture can be also considered (see, e.g., [23, 45, 46]). In Fig. 1(a) an arbitrary deformable medium is shown. The initial configuration C0 ⊂ Rd of the body at time t0 = 0 has a volume Ω0 . The superscript d corresponds to the dimension of the problem, i.e., d ∈ 1, 2, 3. The boundary of the initial configuration is denoted as ∂Ω0 . The medium is subjected to body forces b = {b1 , b2 , b3 } and tractions t̄ applied on ∂Ωt̄ ∈ ∂Ω. Γ0 corresponds to an initial crack within the medium at time t0 . Copyright c 2016 John Wiley & Sons, Ltd. Prepared using nmeauth.cls Int. J. Numer. Meth. Engng (2016) DOI: 10.1002/nme 6 E. G. KAKOURIS, S. P. TRIANTAFYLLOU Under the action of the applied loads, the body undergoes a motion φ, that maps the initial configuration to the current configuration Ct ⊂ Rd at t > t0 with a volume Ωt . Furthermore, the initial crack Γ0 evolves to Γt (Γ0 ⊆ Γt ), which corresponds to the crack path at time t. According to Griffith’s theory, the total potential energy Ψpot of an elastic deformable body with an evolving crack along the path Γ is defined as Ψpot = Ψel + Ψf = Z ψel dΩ + Ω Z Gc dΓ (1) Γ where Ψel is the elastic strain energy, Ψf is the fracture energy, ψel is the elastic energy density and Gc is the critical fracture energy density. The critical fracture energy density is a material parameter corresponding to the energy required to create a unit area of fracture surface [18]. The elastic energy density is readily expressed as a function of the strain field ε according to equation (2) below ψel =   1 λTr2 [ε] + µTr ε2 2 (2) where λ and µ are the Lamé constants. Small strains are considered in this work, with the corresponding strain tensor ε defined as ε=  1 ∇u + ∇uT 2 (3) where the (∇) stands for the gradient operator. Considering a spectral decomposition of the strain tensor as (see, e.g., [20, 16]), the elastic energy density can be established in the following convenient form where the total elastic potential energy is additively decomposed in parts of purely tensile and purely compressive origin, i.e., + − ψel = ψel + ψel Copyright c 2016 John Wiley & Sons, Ltd. Prepared using nmeauth.cls (4) Int. J. Numer. Meth. Engng (2016) DOI: 10.1002/nme PHASE-FIELD MATERIAL POINT METHOD FOR BRITTLE FRACTURE 7 + is the elastic energy density due to tension where ψel + = ψel h
2 i 1 2 λhTr [ε]i + µ Tr ε+ 2 (5) − and ψel is the elastic energy density due to compression − = ψel h
2 i 1 2 λ(Tr [ε] − hTr [ε]i) + µT r ε− 2 (6) respectively. The h.i symbol in relations (5) and (6) denotes the Macaulay brackets hX i = X for
X > 0 and hX i = 0 for X ≤ 0 . The positive part of the strain tensor ε+ in equation (5) is defined through the following spectral decomposition ε + = P Λ+ P T (7) where P is a matrix whose columns comprise the eigen-vectors of the strain tensor ε and Λ+ is a diagonal matrix defined as Λ+ = diag (hλ1 i, hλ2 i, hλ3 i) (8) where λi , i = 1 . . . 3 are the eigen-values of the strain tensor. The negative part of strain tensor is evaluated as ε− = ε − ε+ (9) Substituting equation (4) in relation (1), the expression for the brittle fracture potential energy assumes the following form Ψpot = Z Ω Copyright c 2016 John Wiley & Sons, Ltd. Prepared using nmeauth.cls + dΩ + ψel Z Ω − dΩ + ψel Z Gc dΓ (10) Γ Int. J. Numer. Meth. Engng (2016) DOI: 10.1002/nme 8 E. G. KAKOURIS, S. P. TRIANTAFYLLOU Equation (10) effectively decomposes the total potential energy into purely tensile, compressive and fracture energy parts, thus constituting an efficient platform for the phase field derivations described in the following Section. 2.2. Phase-field fracture As evaluation of the surface crack energy in equation (10) requires prior knowledge of the crack path Γ, computational fracture mechanics revert to crack tracking algorithms in order to identify the crack-path during the solution procedure [47]. To avoid such procedures, the phase field method approximates the path integral of the fracture energy with a volume integral defined over the entire domain of the deformable medium according to the following expression [14] Z Gc dΓ ≈ Γ Z Gc Zc dΩ (11) Ω where Zc is a crack density functional. Several expressions are provided in the literature for the definition of Zc involving the case of second-order [14, 48] and fourth-order functionals. The latter have been found to allow for increased convergence rates [49]. Very recently, an anisotropic definition for the crack surface functional addressing the case of anisotropic fracture has been introduced [50]. In this work, the second order definition for crack density functional Zc of equation (12) [14] is adopted to facilitate verification of the proposed method. However, utilization of higherorder functionals is a straightforward procedure. " 2 (c − 1) + l0 |∇c|2 Zc = 4l0 # (12) In equation (12) c(x, t) ∈ [0, 1] is a phase-field defined over the domain Ω. By considering minimization of the functional with respect to c it can be shown that a value of c = 1 corresponds to un-cracked regions of the domain Ω, i.e., regions away from the crack Γ. Similarly, values of c = 0 are retrieved on regions coinciding with the crack-surface Γ. Involving the gradients of the phase-field on the functional definition introduces a smooth variation of the phase-field from 0 to 1. Copyright c 2016 John Wiley & Sons, Ltd. Prepared using nmeauth.cls Int. J. Numer. Meth. Engng (2016) DOI: 10.1002/nme PHASE-FIELD MATERIAL POINT METHOD FOR BRITTLE FRACTURE 9 The width of the region over which this smooth transition takes place is controlled by the length scale parameter l0 ∈ R (see Fig. 1(b)). The length scale parameter l0 can be considered to correspond to a domain of degrading material parameters in the vicinity of the crack surface. From a purely mathematical standpoint, l0 is a regularization parameter with values of l0 → 0 allowing for the phase field theory to practically converge to Griffith’s theory. In practise, convergence is achieved by using a finite value for l0 . In view of relation (11), the potential energy introduced in equation (10) assumes the following form Ψpot = Z Ω + dΩ ψel + Z − dΩ ψel + Ω Z Gc Zc dΩ (13) Ω Having established through relation (11) that as a crack propagates within the domain Ω the value of the crack surface energy integral will be increasing, the corresponding decrease in the elastic energy due to the degradation of the material properties needs also to be considered in the vicinity of the crack Γ. This is achieved by introducing a degradation function g(c) that is superimposed on the positive part of the elastic strain energy density. The degradation function should be continuously differentiable, monotonically decreasing with properties g(0) = 0, g(1) = 1 and g ′ (0) = 0 as explained in [51]. To facilitate verification of the proposed procedure, the degradation function introduced in [48] is utilized herein, i.e., g = (1 − k)c2 + k (14) where 0 ≤ k ≪ 1 is a model parameter introduced in [52] to avoid ill-posedness. According to the arguments provided in [53] as well as the numerical investigations presented in [16] this parameter can be considered redundant. Results derived from our set of numerical experiments also seem to agree with the aforementioned. Therefore in this work also k = 0. In view of the aforementioned, Copyright c 2016 John Wiley & Sons, Ltd. Prepared using nmeauth.cls Int. J. Numer. Meth. Engng (2016) DOI: 10.1002/nme 10 E. G. KAKOURIS, S. P. TRIANTAFYLLOU the expression for the elastic potential energy finally assumes the following form Ψpot = Z + dΩ gψel + Z − dΩ ψel + Gc Zc dΩ (15) Ω Ω Ω Z The elastic stress field on the medium is readily derived from the elastic potential [32] through the following relation σ=g ∂ψe− ∂ψe+ + = gσ + + σ − ∂ε ∂ε (16) where σ + = λhTr [ε]iI + 2µ  ε+
 (17) and σ − = λ(Tr [ε] − hTr [ε]i)I + 2µ  ε−
 (18) respectively, whereas I denotes the 3x3 identity matrix. The damage elastic tangent constitutive matrix can be analytically derived as D= ∂σ ∂ε (19) using any symbolic programming language. By superimposing the effect of the degradation parameter on the positive part of the energy density only, crack-propagation under compressive stresses is a priori avoided. Finally, considering the Euler-Lagrange equations of both the displacement u and phase field c, the coupled strong form of the brittle-fracture phase field formulation is established as     ∇ · σ + b = ρü,       4l0 (1−k)H + 1 c − 4l02 ∆c = 1, Gc on [Ωt0 , Ωt ] (20) on [Ωt0 , Ωt ] + obtained in time space [t0 , t]. The history where H is the history field defined as the maximum ψel field H (see, e.g., [48] ) essentially enforces the necessary irreversibility condition pertinent to Copyright c 2016 John Wiley & Sons, Ltd. Prepared using nmeauth.cls Int. J. Numer. Meth. Engng (2016) DOI: 10.1002/nme 11 PHASE-FIELD MATERIAL POINT METHOD FOR BRITTLE FRACTURE the crack propagation problem, i.e., Γ (t) ⊆ Γ (t + ∆t) and satisfies the following Kuhn-Tucker conditions for elastic loading and unloading, i.e., ψe+ − H ≤ 0 Ḣ ≥ 0 Ḣ (ψe+ − H) = 0 The coupled field equations (20) are subject to the following set of boundary and initial conditions     σn = t̄,           u = ū,          u = u0 ,      u̇ = u̇0 ,          ü = ü0 ,          ∇c · n = 0,         c = c0 , on [∂Ωt̄0 , ∂Ωt̄ ] on [∂Ωū0 , ∂Ωū ] on Ω0 (21) on Ω0 on Ω0 on [∂Ωt0 , ∂Ωt ] on Ω0 where n is the outward unit normal vector of the boundary, ū is the prescribed displacement field on ∂Ωū boundary, u̇ is the velocity field, ü is the acceleration field and ρ is the mass density. 3. MATERIAL POINT METHOD FOR BRITTLE FRACTURE 3.1. Material Point Method approximation In MPM, a deformable body is approximated with a set of material points p = 1, 2, . . . , Np , where Np ∈ Z is the total number of material points. The material point discretization can be defined by any appropriate tessellation of Ω. Under the action of φ, the initial position vector xp0 of a material point is mapped to the current position vector xpt at t > 0 (see Fig. 1(b)). Consequently, the current position of a material point always depends on the initial position and time t. The displacement vector of the material point is defined as utp = xpt − xp0 ∗ . ∗ Subscript p refers to the material point value. Copyright c 2016 John Wiley & Sons, Ltd. Prepared using nmeauth.cls Int. J. Numer. Meth. Engng (2016) DOI: 10.1002/nme 12 E. G. KAKOURIS, S. P. TRIANTAFYLLOU (a) (b) Figure 1. (a) Deformation process (b) Material point discretization Through this discretization, the mass density distribution of the deformable body is readily defined as ρ (x, t) = Np X ρp δ (x − xpt ) (22) p=1 where ρp = Mp /Ωp is the mass density of the material point, Mp is the material point mass, Ωp is the material point volume and δ is the Dirac function. Similarly, the domain volume Ω (x, t) is additively decomposed into the corresponding material point domain contributions according to the following expression Ω (x, t) = Np X Ωp δ (x − xpt ) (23) p=1 In this work, the tributary volumes Ωp of each material point are defined according to the following methodology. An isoparametric descritization of the material domain is first performed using quadrilateral elements. Material points are then defined at the positions of the Gauss points of each individual element at their natural coordinate system. The corresponding volumes are then mapped back to the Cartesian system by means of the isoparametric transformation. Defining material points at the Gauss points of the corresponding finite element mesh has been chosen to facilitate comparison against standard finite element implementation. Copyright c 2016 John Wiley & Sons, Ltd. Prepared using nmeauth.cls Int. J. Numer. Meth. Engng (2016) DOI: 10.1002/nme 13 PHASE-FIELD MATERIAL POINT METHOD FOR BRITTLE FRACTURE 3.2. Eulerian Mesh Integral to the computational scheme of MPM is the definition of a computational grid where the solution of the governing equations of motion is performed. This grid is termed the Eulerian Mesh (or Eulerian Grid). The Eulerian Grid is a non-deforming mesh corresponding to the space that the material points move through (see. Fig. 2). Referring to Fig. 2, the Eulerian Mesh is divided to active cells, i.e., cells where one or more material points exist at a certain time t and inactive cells where no material point exists. In this work, the Eulerian grid is constantly updated according to the topology of the material points, thus reducing the solution space at any time instant. 3.3. Equilibrium discrete equations The discrete form of MPM is derived through a Galerkin approximation of the strong from defined in equation (20). In this work only quasi-static problems are considered. The weak form is readily derived as [32] by weighting the equilibrium equation (20) with an arbitrary set of weighting functions w and partially integrating over the domain Ω. Z (σ : ∇w) dΩ = Ω Z (t̄ · w) d∂Ωt̄ + ∂Ωt̄ Z (b · w) dΩ (24) Ω The weighting functions w satisfy the essential boundary conditions of the problem. Substituting equation (23) into the weak form (24), the following relation is established Np X (σp : ∇wp )Ωp = p=1 Z (t̄ · w) d∂Ωt̄ + ∂Ωt̄ Np X (bp · wp )Ωp (25) p=1 where wp , σp and bp are the test function, stress field and body forces evaluated at material point xp . Relation (25) which essentially collocates the weak equilibrium of the continuum into the material points derived from the tessellation of the deformable body. Considering the following Galerking interpolation scheme for the test functions and their spatial derivatives Copyright c 2016 John Wiley & Sons, Ltd. Prepared using nmeauth.cls Int. J. Numer. Meth. Engng (2016) DOI: 10.1002/nme 14 E. G. KAKOURIS, S. P. TRIANTAFYLLOU wp = Nn X NI (xp )wI (26) ∇NI (xp )wI (27) I=1 and ∇wp = Nn X I=1 respectively, where Nn is the number of grid nodes, NI (xp ) are the interpolation functions evaluated at the material points and wI are the test function nodal values † . Both the interpolation functions and the test function nodal values are defined with respect to the underlying Eulerian mesh as described in Section 3.2. Thus, standard finite element interpolation functions can be utilized to interpolate material point defined quantities at the nodal points of the corresponding parent cell. In this work, bi-linear shape functions are used although implementation of higher order shape functions can also be made in a straight-forward manner, see, e.g., [54]. Substituting (26) and (27) in relation (25) and performing the necessary algebraic manipulation, the following expression is derived Nn X I=1   wI · FIint − FIext = 0 (28) where int FI,i = Np X (σpjk · BIijk (xp ))Ωp (29) p=1 FIext = Z (t̄NI (xp )) d∂Ωt̄ + ∂Ωt̄ Np X bp NI (xp )Ωp (30) p=1 In relation (29), σpjk denotes the stress components σp = {σpjk } whereas BI (xp ) = {BIijk (xp )} is defined as BIijk (xp ) = † Subscript 1 2  ∂NI (xp ) ∂NI (xp ) δik + δij ∂xj ∂xk  (31) I refers to grid node value. Copyright c 2016 John Wiley & Sons, Ltd. Prepared using nmeauth.cls Int. J. Numer. Meth. Engng (2016) DOI: 10.1002/nme PHASE-FIELD MATERIAL POINT METHOD FOR BRITTLE FRACTURE 15 As the test functions in equation (25) are chosen arbitrarily, subject to the essential conditions of the problem, equation (28) should hold for every set of nodal values wI . Thus, the following equilibrium equation is finally established RIu (uI ) = FIint − FIext = 0, I = 1 . . . , Nn (32) where FIint is the vector of corresponding internal forces, FIext corresponds for the equivalent vector of external forces evaluated at grid node I and RIu (uI ) are the residual nodal values for the displacement field and uI is the displacement at the grid node I . Equation (32) corresponds to nodal force equilibrium established at the nodes of the background mesh with the material point to background node mapping performed through relations (29) and (30) for the internal and external forces respectively. Further considering the strain-displacement relation defined in equation (3), the strain components εp = {εpjk } in each material point can be expressed as εpjk = Nn X BIijk (xp )uI,i (33) I=1 Substituting Eqs. (19) and (33) into Eq. (29) and using Eq. (30), the following compact form is eventually derived K u u = F ext (34) u where K u is the global stiffness matrix of the structure whose KI,J,i,j component is expressed as u KI,J,i,j = Np  X p=1 Dplkmn BJjmn (xp )  · BIilk (xp ) ! Ωp (35) The term Dplkmn , accounts for the components of the constitutive matrix Dp evaluated at pth material point. Copyright c 2016 John Wiley & Sons, Ltd. Prepared using nmeauth.cls Int. J. Numer. Meth. Engng (2016) DOI: 10.1002/nme 16 E. G. KAKOURIS, S. P. TRIANTAFYLLOU 3.4. Phase field discrete equations The weak form of the phase-field governing equation assumes the following form Z Z Z   4l0 (1 − k)H + 1 cq dΩ + 4l02 (∇c : ∇q) dΩ = q dΩ Gc Ω Ω Ω (36) where c is the phase field and q are the corresponding weighting functions for the phase-field. In this work, the continuous phase field and the corresponding weighting functions introduced in equation (36) are collocated at material points, resulting in the following discrete form (37) Np X Fp cp qp Ωp + p=1 Np X 4l02p (∇cp : ∇qp )Ωp = p=1 Np X q p Ωp (37) p=1 where cp and qp are values of the the phase-field and weighting functions respectively at the material point p. Fp is defined as Fp = 4l0p (1 − kp )Hp +1 G cp (38) where l0p , kp , Hp and Gcp are the length scale parameter, model parameter, history field and critical fracture energy density of material point xp . Next, both cp and qp and interpolated at the nodal points of the background mesh. The value of the test function and its spatial derivatives at the pth material point are expressed as qp = Nn X NI (xp )qI (39) ∇NI (xp )qI (40) I=1 and ∇qp = Nn X I=1 respectively, where NI (xp ) are the background mesh shape functions pertinent to the phase field interpolation and qI are nodal values of the corresponding test functions. Copyright c 2016 John Wiley & Sons, Ltd. Prepared using nmeauth.cls Int. J. Numer. Meth. Engng (2016) DOI: 10.1002/nme 17 PHASE-FIELD MATERIAL POINT METHOD FOR BRITTLE FRACTURE The value of the phase field at the pth material point is written as cp = Nn X NI (xp )cI (41) I=1 where cI are the phase field nodal values. Attention is drawn to the fact that the phase-field is a scalar quantity. It follows from relation (41) that the gradients of the the phase field at the pth material point are defined accordingly as ∇cp = Nn X ∇NI (xp )cI (42) I=1 Similar shape functions are considered for both the phase-field and the corresponding weighting functions according to the Galerkin approximation [32]. Furthermore, in this work, the same family interpolation functions is considered for both the displacement and the phase field for brevity (see also [48]). Substituting equations (39) and (40) in relation (37) and re-arranging terms, the following expression is derived Nn X qI · I=1 h (I) i S1 +(I) S2 −(I) S3 = 0 (43) where in equation (43) above (I) S1 = Np X Fp cp NI (xp )Ωp (44) 4l02p (∇cp · ∇NI (xp ))Ωp (45) p=1 while (I) S2 = Np X p=1 and (I) S3 = Np X NI (xp )Ωp (46) p=1 respectively. Since the choice of the weighting functions is arbitrary, it must hold that RIc (cI ) =(I) S1 +(I) S2 −(I) S3 = 0, Copyright c 2016 John Wiley & Sons, Ltd. Prepared using nmeauth.cls I = 1 . . . , Nn (47) Int. J. Numer. Meth. Engng (2016) DOI: 10.1002/nme 18 E. G. KAKOURIS, S. P. TRIANTAFYLLOU where RIc (cI ) are the residual nodal values for the phase field. By substituting relations (44) to (46) into equation (47) results in Np X Fp cp NI (xp )Ωp + Np X 4l02p (∇cp · ∇NI (xp ))Ωp = NI (xp )Ωp , I = 1 . . . , Nn (48) p=1 p=1 p=1 Np X Further considering the phase-field interpolation schemes defined in equations (41) and (42), and substituting in equation (48), the following relation is established Np X p=1 Fp Nn X ! NJ (xp )cJ NI (xp )Ωp + J=1 Np X Nn X 4l02p p=1 ∇NJ (xp )cJ J=1 ! ! · ∇NI (xp ) Ωp = Np X NI (xp )Ωp p=1 (49) Re-arranging and collecting terms, equation (49) gives rise to the following convenient form K cc = F c (50) c component is defined as where K c is an (Nn × Nn ) coefficient matrix whose KI,J c KI,J = Np X  Fp NJ (xp )NI (xp ) + 4l02p ∇NJ (xp ) · ∇NI (xp ) p=1  ! Ωp (51) while c is the (Nn × 1) vector of unknown nodal phase fields and F c is the (Nn × 1) vector whose FIc component is defined as FIc = Np X NI (xp )Ωp (52) p=1 The vector quantity F c will be termed herein as the phase-field forcing term. Similarly to the MPM displacement based equilibrium equations defined in Section 3.3, equation (50) is established and solved at the nodal points of the background mesh with the corresponding material point to background node mapping performed in equations (44) to (46). Copyright c 2016 John Wiley & Sons, Ltd. Prepared using nmeauth.cls Int. J. Numer. Meth. Engng (2016) DOI: 10.1002/nme PHASE-FIELD MATERIAL POINT METHOD FOR BRITTLE FRACTURE 19 4. NUMERICAL IMPLEMENTATION 4.1. Material Point Method Computational Cycle The computational cycle of MPM is divided into six main steps (Fig. 3). These are: (i) Initialisation phase (ii) Mapping from material points to grid nodes (Fig. 3(a)) (iii) Solving governing equations in grid nodes (Fig. 3(b)) (iv) Mapping from grid nodes to material points (Lagrangian Phase - Fig. 3(c)) (v) Update material point properties (Fig. 3(d)) (vi) Reset the computational grid (Convective Phase) In the first step the computational grid and material point configuration and properties are defined. Then the solution phase begins. In this, the displacements, strains and stresses defined at the material points are interpolated at the nodes of their corresponding background parent cell. Following, the governing equations are formulated and solved at the background grid in an updated Lagrangian fashion. Next, the solution is mapped back to the material points. This is termed the Lagrangian phase of the MPM. Finally, the material point properties are updated and the computational grid is reset. The latter is termed the MPM Convective Phase. In the convective phase of the method, the displacements evaluated on the background mesh are disregarded and the initial background mesh configuration is reused. In the implementation used in this work, background cells that are found at this stage to contain no material points are considered inactive, thus reducing the order of the system to be solved at the next cycle. The computational cycle of material point method is illustrated in Fig. 3. Further information on the numerical scheme of MPM can be found in [31, 38]. It should be noted that compared to the standard finite element implementation, the implicit material point method does introduce additional computational costs as a re-factorization of the stiffness matrix introduced in equation (34) is required when material points move across background cells. However, this results in a high-fidelity procedure for simulating complex Copyright c 2016 John Wiley & Sons, Ltd. Prepared using nmeauth.cls Int. J. Numer. Meth. Engng (2016) DOI: 10.1002/nme 20 E. G. KAKOURIS, S. P. TRIANTAFYLLOU (a) (b) Figure 2. (a) Initial Configuration (b) Deformed Configuration Figure 3. Computational Cycle (MPM) phenomena, (see, e.g., [34, 41, 44]) and also provides very accurate estimates for the crack geometry as will shown in Section 5. In this work, both the displacement field and the phase field which are defined at material points are mapped at the corresponding parent cell nodes resulting in an updated MPM Lagrangian phase. 4.2. Phase field Material Point Method solution scheme The coupled equilibrium and phase field evolution equations can be solved in a so called monolithic fashion, i.e., simultaneously within each incremental step. However, it has been demonstrated that a staggered solution approach (see also, [55] for the case of thermo-mechanical coupling) can be utilized where the phase field equations are solved independently and the resulting phase field prediction is then used to iteratively solve for the equilibrium equations [48]. In this work, a staggered solution procedure has been implemented and the corresponding computational scheme is presented in Algorithm 1. Copyright c 2016 John Wiley & Sons, Ltd. Prepared using nmeauth.cls Int. J. Numer. Meth. Engng (2016) DOI: 10.1002/nme 21 PHASE-FIELD MATERIAL POINT METHOD FOR BRITTLE FRACTURE With regards to Algorithm 1, Ep , vp and gp refer to the Young’s modulus, Poisson’s ratio and degradation function of the pth material point. Pre-existing cracks can be modelled by defining an initial history field (Hp0 ) in all material points around the crack similar to [16]. Alternately, preexisting cracks can also be introduced as discrete cracks in the geometry of the structure. A displacement control incremental analysis procedure is implemented in this work for the solution of the quasi-static brittle fracture problem, considering a set of Nsteps incremental steps. In the beginning of each time step h, the active cells of the Eulerian Grid are identified according to the material point positions and the inactive cells are discarded (Fig. 2). Next, the total number of grid nodes as well as grid degree of freedom are redefined according to total active grid nodes (Nn ) and total active unconstrained grid degree of freedom (Ndof s ). Furthermore, the basis functions (N (xp )) as well as their derivatives (∇N (xp )) at all material points need to be evaluated at each time step h. This is one of the main differences between Finite Element Method and Material Point Method as in the former the number of nodes, degree of freedom, cells (Ncells ) as well as basis functions and their derivatives remain constant during the analysis. Following, the staggered iterative scheme (k = 1, 2, . . . , Nstaggs ) initiates within the current incremental step. The phase field equations are solved for the current value of the history field H and the phase field nodal values cI are derived. Using this phase field prediction, updated values for the degradation function at each material point gp are derived and the displacement field equations are iteratively solved in the inner iterative loop (j = 1, 2, . . . , Niters ). From this, the incremental displacement field nodal values ∆uI are obtained. The displacement field equations (32) are solved by incrementally applying the external forces ∆FIext to obtain the increments of the displacement field ∆uI . and the following equations are solved using a Newton-Raphson method (inner iterations j = 1, .., Niters ). δRIu (∆uI ) = ∆FIint − ∆FIext = 0, Copyright c 2016 John Wiley & Sons, Ltd. Prepared using nmeauth.cls I = 1 . . . , Nn (53) Int. J. Numer. Meth. Engng (2016) DOI: 10.1002/nme 22 E. G. KAKOURIS, S. P. TRIANTAFYLLOU where the symbol ∆ denotes incremental quantities e.g. ∆X = X(h) − X(h−1) , whereas the symbol δ denotes iterative quantities e.g. δX = X (j) − X (j−1) . Convergence of the equilibrium equation iterative procedure is achieved when the Euclidean norm of the residual force vector introduced in equation (53) assumes a sufficiently small value, i.e., when kδRu(j) k ≤ tolu . Upon convergence, updated values for the history field H are evaluated and the residual of the phase field equation is established as the difference between the initial phase field forcing term estimate and the updated one. Outer, phase field iterations terminate c(k) when kR(h) k ≤ tolc where tolc is a predefined tolerance. Although robust, the staggered until convergence scheme is prone to low convergence rates and practically bounds the maximum allowable incremental displacement step in a displacement controlled analysis. Very recently, a line-search assisted iterative scheme has been developed to treat such issues and further improve the convergence speed of the method [56]. 5. NUMERICAL EXAMPLES In this Section the proposed method is compared against the finite element phase field implementation through a set of representative tests both in terms of accuracy and computational efficiency. In all cases external loads are directly applied at material points. Kinematical constraints are imposed by means of the Penalty Method [32]. As these constraints are imposed on the material points rather than the background grid, the corresponding numerical implementation is presented in Appendix A. This is contrary to the finite element implementation where essential boundary conditions are imposed directly on the domain boundary. However, as shown from the actual verification results provided this does not affect the accuracy of the method. The density of material points utilized as necessitated by the fracture propagation problem ensures that material points are sufficiently close to the actual domain boundary where displacement variations can be considered negligible. Several methods have been examined for direct and accurate implementation of essential boundary conditions in MPM, either directly at material points, see, e.g., [57] or at the background mesh , see, e.g., [58]. Copyright c 2016 John Wiley & Sons, Ltd. Prepared using nmeauth.cls Int. J. Numer. Meth. Engng (2016) DOI: 10.1002/nme 23 PHASE-FIELD MATERIAL POINT METHOD FOR BRITTLE FRACTURE Data: Define computational grid, material point properties (xp(0) , Ωp(0) , Ep , vp , l0p , kp Gcp , Hp(0) , σp(0) , εp(0) ) for each time step h = 1, 2, .., Nsteps do Reset the computational grid: Find active part of Eulerian Grid, Nn , Ndof s , Ncells ; Compute: N (xp(h) ), ∇N (xp(h) ), BI (xp(h) ), for all material points. ; ext Define: δRu(1) = ∆F(h) ; for each stagger iteration k = 1, 2, .., Nstaggs do c(k) Compute: F(h) (see Eq. (52)). ; Compute: K c (see Eq. (51)). ; (k) c(k) Solve: K c c(h) = F(h) ; (k) (k) (k) (k) Map phase field (c(h) ) from grid nodes to material points. Evaluate: cp(h) , ∇cp(h) , gp(h) , for all material points (see Eq. (41), (42) and (14)). ; Initialize ∆u(0) = 0 ; for each inner iteration j = 1, 2, .., Niters do Compute: K u (see Eq. (35), for constitutive matrix see Eq. (19)) ; Solve: K u δu(j) = δRu(j) , with displacement contol. ; Compute: ∆u(j) = ∆u(j−1) + δu(j) ; (j) Compute: ∆εp , for all material points (see Eq. (33)). ; (j) (j) Compute: εp(h) = εp(h−1) + ∆εp , for all material points. ; (j) Compute: σp(h) , for all material points (see Eq. (16)) ; Compute: ∆F int(j) = {∆FIint }, PNp (j) (j) Ωp(h) (σp(h) − σp(h−1) ) · BI (xp(h) ) ; ∆FIint = p=1 ext Compute Residual (Displacement-Field): δRu(j) = ∆F(h) − ∆F int(j) ; Convergence Check (Displacement Field): If kδRu(j) k ≤ tolu or j ≥ Niters then ”exit” from loop else j = j + 1 go to next inner iteration. ; end + , for all material points (see Eq. (5)) Compute: ψel p(h)   ψ + , for ψ + elp(h) elp(h) > Hp(h−1) → Hp(h) = ;   Hp , otherwise (h−1) (k) c(k) (k) Compute Residual (Phase-Field): R(h) (see Eq. (47)) according to cp(h) , ∇cp(h) , Hp(h) ; c(k) Convergence Check (Phase Field): If kR(h) k ≤ tolc or k ≥ Nstaggs then ”exit” from loop else k = k + 1 go to next stagger iteration. ; end Compute: ∆up(h) = PNn I=1 (j) NI (xp(h) )∆uI , for all material points. ; Compute: up(h) = up(h−1) + ∆up(h) , for all material points. ; Compute: xp(h) = xp(h−1) + ∆up(h) , for all material points. ; end Algorithm 1: Phase-Field Material Point Method pseudo-code (Stagger Solution Algorithm). Copyright c 2016 John Wiley & Sons, Ltd. Prepared using nmeauth.cls Int. J. Numer. Meth. Engng (2016) DOI: 10.1002/nme 24 E. G. KAKOURIS, S. P. TRIANTAFYLLOU For the purpose of verification both the Phase Field Material Point Method (PF-MPM) and Phase Field Finite Element Method (PF-FEM) have been implemented in Fortran code. In all examples considered herein, the staggered until convergence solution strategy (see also Algorithm 1) was adopted. The phase field residual tolerance was set at tolc = 1.0e − 6. Simulation parameters, i.e., number of incremental steps, convergence tolerance and maximum number of iterations are similar for both schemes as defined in the corresponding Sections below. All tests were performed in a PC fitted with an Intel Xeon E5-1620 CPU and 32 GBs of RAM. 5.1. Single-edge notched tension test In this example a square plate under pure tension is examined and results are compared to the standard PF-FEM. The geometric configuration, boundary conditions and material parameters considered are presented in Fig. 4(a). The square plate consists of 249000 material points. The Eulerian grid is formulated by 67600 (260 x 260) 4-node isoparametric quadrilateral elements with a uniform mesh size equal to h = 0.004 mm. Element size of the background mesh is defined such that h < l0 . The overall dimensions of the Eulerian grid are 1.04 mm x 1.04 mm (xmin = ymin = −0.02 mm , xmax = ymax = 1.02 mm). Material points are initially located at the Gauss point position of their corresponding parent cells and plane strain conditions are assumed. For the PF-FEM case, the corresponding finite Element mesh comprises 62250 four node quadrilateral plane strain elements with bilinear basis functions. Full integration is considered in each element with 4 Gauss points. The material parameters considered are E = 210 kN/mm2 , v = 0.30, l0 = 0.0075 mm and Gc = 0.0027 kN/mm for the Young’s modulus, Poisson’s ratio, length scale and fracture energy density respectively. Zero displacement boundary conditions, i.e., upx = upy = 0, are imposed in all material points (nodes in the PF-FEM case) in the bottom edge of the specimen. Both in the PF-MPM and PF-FEM implementations, a displacement control nonlinear static analysis scheme is utilized with a constant displacement increment ∆u = 10−6 mm. Displacement is monitored Copyright c 2016 John Wiley & Sons, Ltd. Prepared using nmeauth.cls Int. J. Numer. Meth. Engng (2016) DOI: 10.1002/nme 25 PHASE-FIELD MATERIAL POINT METHOD FOR BRITTLE FRACTURE and controlled in the upper edge of the specimen where the vertical displacements of all material points (nodes in the PF-FEM cases) are kinematically constrained. The kinematic constraint penalty parameter (see Appendix A) was chosen to be a = 10000000. The solution is implemented within a stagger solution algorithm with a single prediction step (Nstaggs = 1) and tolu = 10−5 . The load paths derived from both PF-MPM and PF-FEM are presented in Fig. 4(b). The load paths are practically identical. Results obtained by both solution approaches also agree with the results provided in [48]. In particular, the critical vertical displacement and critical load obtained by PF-FEM are ucrP F −F EM = 0.005626 mm and FcrP F −F EM = 0.7051 kN respectively. The critical vertical displacement and critical load obtained by PF-MPM are ucrP F −M P M = 0.005627 mm and FcrP F −M P M = 0.7052 kN respectively. The results derived and the agreement between the two different approaches is justified by the fact that due to small displacements (≈ 0.60% of the total length of the plate), material points only marginally move from the Gauss points of the corresponding Finite Element mesh. The phase field distribution over the plate domain for both the Finite Element and Material Point Method are presented in Figs. 5. The observed crack paths derived from both methods are identical. Analysis time for PF-FEM was approximately 98 hrs whereas for PF-MPM 111 hrs. The increase in computational time due to the MPM implementation was of the order of 13% corresponding to the re-factorization of the stiffness matrix when material points move across background cells. 5.2. Single-edge notched shear test In this case, the response of the square plate considered in Section 5.1 is investigated under pure shear conditions. The same example has been previously examined in [48] and [16] considering a standard Finite Element scheme and its isogeometric formulation respectively. The geometry, boundary conditions and material properties are shown in Fig. 6(a). The discretization of the Eulerian Grid as well as the number of material points are the same as in the tension experiment of Section 5.1. Both in the PF-MPM and PF-FEM implementation the simulation is performed with a constant horizontal displacement increment ∆u = 10−5 , mm monitored at the upper edge of the plate. Copyright c 2016 John Wiley & Sons, Ltd. Prepared using nmeauth.cls Int. J. Numer. Meth. Engng (2016) DOI: 10.1002/nme 26 E. G. KAKOURIS, S. P. TRIANTAFYLLOU (a) (b) Figure 4. (a) Tension test. Geometry and boundary conditions. (b) Load-Displacement curve. Comparison between Material Point Method and Finite Element Method (FEM) The load displacement curve is presented in Fig. 6(b). The results obtained by PF-MPM are compared to the results from PF-FEM. The latter have been derived considering a 40658 constant strain triangle finite element mesh and are in perfect agreement with the results reported in [48]. The critical horizontal displacement and critical load obtained by Finite Element Method are ucrP F −F EM = 0.0087 mm and FcrP F −F EM = 0.5310 kN respectively; whereas the critical vertical displacement and critical load obtained by Material Point Method are ucrP F −M P M = 0.0089 mm and FcrP F −M P M = 0.5416 kN respectively. Figs. 7 illustrate the phase field of both Finite Element Method and Material Point Method. The two methods illustrate good agreement in regards of the crack path with minor differences observed in the post-peak regime. These are attributed to the severe distortion of the triangular finite elements observed in the FEM-PF case which is however by definition avoided in the MP-PF solution. In particular, the distortion of the elements along the crack path is presented in Figs. 8(a)-(c). This is avoided in the material point method as shown in Figs. 8(d)-(f) as the material points naturally follow the geometry of the crack. The evolution of the hydrostatic stress for the case of the Material Point implementation is shown in Figs. 9 for several timesteps. Comparing Figs. 7(d)-(f) to Figs. 9 one is able to verify that the crack is propagates only due to tension as a result of additive decomposition of the elastic energy introduced in equation (15). The computational times for PF-FEM and PF-MPM were approximately 23 hrs and 26 hrs respectively resulting in an overhead of approximately 13%. Copyright c 2016 John Wiley & Sons, Ltd. Prepared using nmeauth.cls Int. J. Numer. Meth. Engng (2016) DOI: 10.1002/nme 27 PHASE-FIELD MATERIAL POINT METHOD FOR BRITTLE FRACTURE (a) (b) (c) (d) (e) (f) Figure 5. Phase-Field in single-edge notched tension test with Finite Element Method for (a) u = 0.0056 mm, (b) u = 0.0058 mm and (c) u = 0.0059 mm respectively. Phase-Field in single-edge notched tension test with Material Point Method for (d) u = 0.0056 mm, (e) u = 0.0058 mm and (f) u = 0.0059 mm respectively. 5.3. L-Shaped panel test In this example an L-Shaped concrete panel is examined under cyclic loading. This example has also been considered in [22] utilizing the phase field finite element scheme. The geometry, boundary conditions and material properties are presented in Fig. 10(a). Herein, a series of simulations is carried out to investigate the effect of the length scale parameter (l0 ), cell spacing of the underlying Eulerian grid (h) and cell density on the accuracy of the PF-MPM scheme. In all simulations (both in the PFMPM and PF-FEM implementation) a constant displacement increment ∆u = 10−3 mm is considered for 2000 time increments. The load history is shown in Fig. 10(b). The solution is implemented within Copyright c 2016 John Wiley & Sons, Ltd. Prepared using nmeauth.cls Int. J. Numer. Meth. Engng (2016) DOI: 10.1002/nme 28 E. G. KAKOURIS, S. P. TRIANTAFYLLOU (a) (b) Figure 6. (a) Shear test. Geometry and boundary conditions. (b) Load-Displacement curve. Comparison between Material Point Method and Finite Element Method (FEM) (a) (d) (b) (e) (c) (f) Figure 7. Phase-Field in single-edge notched shear test with Finite Element Method for (a) u = 0.009 mm, (b) u = 0.011 mm and (c) u = 0.0134 mm respectively. Phase-Field in single-edge notched shear test with Material Point Method for (d) u = 0.0092 mm, (e) u = 0.0117 mm and (f) u = 0.0136 mm respectively. Copyright c 2016 John Wiley & Sons, Ltd. Prepared using nmeauth.cls Int. J. Numer. Meth. Engng (2016) DOI: 10.1002/nme 29 PHASE-FIELD MATERIAL POINT METHOD FOR BRITTLE FRACTURE (a) (b) (c) (d) (e) (f) Figure 8. Distorted Elements in Finite Element Method (a) u = 0.009 mm, (b) u = 0.011 mm and (c) u = 0.0134 mm respectively. No distortion in Material Point Method for (d) u = 0.0092 mm, (e) u = 0.0117 mm and (f) u = 0.0136 mm respectively. (a) (b) (c) Figure 9. Hydrostatic Stress in single-edge notched shear test for (a) u = 0.0092 mm, (b) u = 0.0177 mm and (c) u = 0.0136 mm respectively. Material Points with cp < 0.08 have been removed. Copyright c 2016 John Wiley & Sons, Ltd. Prepared using nmeauth.cls Int. J. Numer. Meth. Engng (2016) DOI: 10.1002/nme 30 E. G. KAKOURIS, S. P. TRIANTAFYLLOU (a) (b) Figure 10. (a) L-Shape panel test. Geometry and boundary conditions. (b) Load History a stagger solution algorithm for a single prediction step (Nstaggs = 1) and tolu = 10−4 . In all cases, the background grid is formulated using four node quadrilateral elements with bilinear basis functions. Plane stress conditions are assumed with a thickness th = 100 mm. Initially, the PF-MPM implementation is compared to a PF-FEM solution considering a length scale parameter l0 = 2.5 mm while the mesh size is h = 2.5 mm. Each background cell in the Material Point implementation is populated with 2 x 2 material points while full integration (4 Gauss points per element) is considered for the Finite Element Method formulation . The corresponding load paths are shown in Fig. 11(a) where the two methods demonstrate very good agreement. The relative divergence of the two methods is presented in Fig. 11(b). Overall, the nominal value of the relative divergence is smaller than 0.1% except from the regions of zero imposed displacement where discrepancies are observed which are attributed to cut-off and round-off errors. The relative divergence significantly increases in the final stages of the loading scenario, i.e., on softening regime of the member response. As in the case of the Shear Test examined in Section 5.2, this pertains to the different kinematics between the two solution procedures, with the PF-MPM implementation providing a more accurate representation of the actual crack path. The L2 norm of the divergence is ǫL2 = 0.3%. The crack pattern between the experimental results and the proposed method demonstrate good agreement as Copyright c 2016 John Wiley & Sons, Ltd. Prepared using nmeauth.cls Int. J. Numer. Meth. Engng (2016) DOI: 10.1002/nme 31 PHASE-FIELD MATERIAL POINT METHOD FOR BRITTLE FRACTURE (a) (b) Figure 11. MPM vs FEM phase field implementation comparison (a) Load paths (b) Relative Error vs applied displacement. The cell density for the MPM method is 2 x 2, whereas in both cases l0 = 2.5mm and h = 2.5mm. (a) (b) Figure 12. Comparison of crack path between experimental ([59], see also [22]) (a) and simulation data (b). The red line represents the crack path obtained from simulation. Material Points with cp < 0.08 have been removed. presented in Figs. 12(a) and 12(b) respectively. The analysis times for the PF-FEM was approximately 15 hrs whereas for PF-MPM 18 hrs resulting in an increase of 20%. Next, the sensitivity of the PF-MPM implementation on the cell density is investigated. Two cell densities are examined, namely 2 x 2 and 4 x 4. In the first case, each active cell contains approximately 2 x 2 = 4 Material Points whereas in the second approximately 4 x 4 = 16 Material Points are used to represent the deformable domain. The total number of material points utilized in each case is 30204 and 120404 respectively. The grid is formulated by 110 x 110 = 12100 cells with a cell spacing h = 5 mm for all runs. Four cases are considered for the length scale, namely l0 = {10, 5, 2.5, 1} mm. Copyright c 2016 John Wiley & Sons, Ltd. Prepared using nmeauth.cls Int. J. Numer. Meth. Engng (2016) DOI: 10.1002/nme 32 E. G. KAKOURIS, S. P. TRIANTAFYLLOU The corresponding load paths are summarized in Figs. 13(a)-(b). Indicative results are presented in Table I where it can be seen that for the same length scale parameter, the influence of the cell density on the both the peak load and the corresponding critical displacement is marginal (less than 1%). To further examine the robustness of the proposed scheme, the sensitivity of the analysis results on the cell spacing (h) is also investigated. Four cases are considered for the size h of the background grid, i.e., h = 10 mm, h = 5 mm, h = 2.5 mm and h = 1 mm whereas the cell density is kept constant at 2 x 2. The corresponding number of cells for cell spacing h = 10 mm, h = 5 mm, h = 2.5 mm and h = 1 mm are 55 x 55 = 3025, 110 x 110 = 12100, 220 x 220 = 48400 and 510 x 510 = 260100 respectively. The derived results are presented in Table II whereas the corresponding load paths are shown in Figs. 14. Contrary to the behaviour identified when varying the cell density, the cell spacing seems to bear a stronger impact on the structural response. Differences in the measured peak load and corresponding displacement presented in Table II significantly increase with decreasing values of the length scale parameter l0 . This is expected as the smaller the value of the length scale parameter for a given cell size, the less the diffusion due to the phase field evolution, thus the more mesh-dependent the crack path is. When h ≤ l0 , the results derived from both the different cell spacings converge; this is in accordance with the remarks presented in [20] regarding the effect of the the length scale to mesh size ratio on the accuracy of the results. The crack paths corresponding to the different length scale parameters are shown in Figs. 15(a)-(d) and 15(e)-(h) for cell spacing h = 5 mm and h = 2.5 mm respectively. The experimentally observed crack path [59], see also [22], is shown in Fig. 12(a). The geometry of the crack path is only marginally affected when h ≤ l0 . However, when h > l0 (see. Figs 15(c),(d) and (h)) the crack pattern diverges from the experimental observation. In Figs. 16 the evolution of the phase field is shown for l0 = 2.5 mm and h = 2.5 mm. Fig. 16(a) represents the degradation of structure when the critical load is observed. A degradation is also observed from time steps u = 0.30 mm to u = 1 mm in the region around of the load due to the cyclic load. In particularly, this region is on tension for load steps 300 until 800. Copyright c 2016 John Wiley & Sons, Ltd. Prepared using nmeauth.cls Int. J. Numer. Meth. Engng (2016) DOI: 10.1002/nme 33 PHASE-FIELD MATERIAL POINT METHOD FOR BRITTLE FRACTURE l0 [mm] 10 5 2.5 1 Cell spacing 5 mm Critical vertical displacement [mm] Critical load [kN] Cell Density Diff Cell Density 2 x2 4 x4 [%] 2x2 4x4 0.255 0.254 0.392 12.5692 12.6099 0.264 0.262 0.757 14.0292 14.0343 0.286 0.284 0.699 15.9206 15.8939 0.345 0.345 0.000 19.2705 19.2039 Table I. L-Shape panel test for different cell density Diff [%] 0.323 0.036 0.167 0.345 5.3.1. Distortion error: PF-FEM vs PF-MPM To assess the advantages of PF-MPM as compated to PF-FEM with regards to mesh distortion errors, the L-shaped panel benchmark is re-run considering a significantly larger critical fracture energy density Gc = 8.9 · 10−1 kN/mm. All the other material parameters, boundary conditions and solution procedure parameters remain unchanged. The cell spacing is chosen to be h = 10 mm with cell density 3 x 3 while the length scale parameter is equal to l0 = 10 mm. A constant displacement increment ∆u = 10−1 mm is considered until complete failure. The resulting load-displacement diagrams of both schemes are represented in Fig. 17. The critical vertical displacement and critical load obtained by Finite Element Method are ucrP F −F EM = 25.6 mm and FcrP F −F EM = 1324.12 kN respectively; whereas the critical vertical displacement and critical load obtained by Material Point Method are ucrP F −M P M = 25.6 mm and FcrP F −M P M = 1316.16 kN respectively. The differences in critical values are less than 0.60%. However, after the critical load (crack initiation) the PF-FEM equilibrium path significantly diverges from the corresponding PF-MPM solution. Both the evolution of phase-field and the deformed configuration of the specimens are shown in Figs. 18 and 19 for the PF-FEM and PF-MPM respectively. From Figs. 18, it is obvious that in PF-FEM the elements are highly distorted especially after the critical load (see Figs. 18(b) and 18(c)). PF-MPM is free of mesh-distortion errors thus allowing for a better representation of the actual crack path with the specimen being able to rotate until complete failure (see Fig. 19(d)). Contrary to PF-MPM, PFFEM fails to converge for displacements u > 278 mm. Figs. 19 also represent the active cells in the corresponding time-steps. Copyright c 2016 John Wiley & Sons, Ltd. Prepared using nmeauth.cls Int. J. Numer. Meth. Engng (2016) DOI: 10.1002/nme 34 E. G. KAKOURIS, S. P. TRIANTAFYLLOU (a) (b) Figure 13. Influence of cell density for length scale parameter (a) l0 = 1 mm and l0 = 2.5 mm (b) l0 = 5 mm and l0 = 10 mm. The cell spacing is h = 5 mm. l0 Cell density 2x2 mm Critical vertical displacement [mm] Cell spacing [mm] Critical load [kN] Cell spacing [mm] [mm] 10 10 5 2.5 1 0.255 0.270 0.310 0.424 5 2.5 1 10 5 2.5 0.254 0.252 0.252 13.323 12.617 12.569 0.264 0.256 0.256 15.213 14.029 13.783 0.286 0.271 0.266 18.046 15.921 15.216 0.345 0.320 0.286 25.156 19.271 17.933 Table II. L-Shape panel test for different cell spacing (a) 1 12.455 13.631 14.816 16.915 (b) Figure 14. Influence of cell spacing h for length scale parameter (a) l0 = 1 mm and l0 = 2.5 mm (a) l0 = 5 mm and l0 = 10 mm. The cell density is 2 x 2. 5.4. Notched Plate with Hole In the final example, a notched plate with hole is examined and compared to the experimental crack path obtained by [22]. The geometry and material parameters are presented in Fig. 20(a). The Copyright c 2016 John Wiley & Sons, Ltd. Prepared using nmeauth.cls Int. J. Numer. Meth. Engng (2016) DOI: 10.1002/nme 35 PHASE-FIELD MATERIAL POINT METHOD FOR BRITTLE FRACTURE (a) (b) (c) (d) (e) (f) (g) (h) Figure 15. Phase-Field in L Panel Test for cell spacing h = 5 mm ((a)-(d)) and h = 2.5 mm ((e)-(h)) respectively. Fig.s (a) and (e) are for l0 = 10 mm, (b) and (f) are for l0 = 5 mm,(c) and (g) are for l0 = 2.5 mm and (d) and (h) are for l0 = 1 mm. (a) (b) (c) (d) Figure 16. Phase-Field in L-Panel test for l0 = 2.5 mm and h = 2.5 mm. (a) u = 0.27 mm, (b) u = 0.30 mm, (c) u = 0.45 mm and (d) u = 1.00 mm. rectangular plate consists of 178236 material points. The Eulerian grid is formulated by 73080 (203 x 360) cells with cell spacing hx ≈ 0.3497 mm and hy = 0.35 mm in x and y direction respectively. Copyright c 2016 John Wiley & Sons, Ltd. Prepared using nmeauth.cls Int. J. Numer. Meth. Engng (2016) DOI: 10.1002/nme 36 E. G. KAKOURIS, S. P. TRIANTAFYLLOU Figure 17. MPM vs FEM phase field implementation comparison for h = 10 mm, l0 = 10 mm and Gc = 8.9 · 10−1 kN/mm. (a) (b) (c) Figure 18. Phase-Field Finite Element Method in L-Panel test for h = 10 mm, l0 = 10 mm and Gc = 8.9 · 10−1 kN/mm. (a) u = 25.60 mm, (b) u = 90 mm, (c) u = 278 mm. The dimensions of the grid are 71 mm x 126 mm (xmin = ymin = −3.00 mm , xmax = 68 mm ymax = 123 mm). The cell spacing is chosen in order to be less or equal than the length scale parameter l0 . Four node cells with bilinear basis functions are used for the grid. The active cells at the beginning of the analysis are shown in Fig. 20(b). Plane stress conditions are assumed. Material points are randomly distributed in cells with the cell density varies from 1 to 4 material points per cell. Zero displacement boundary conditions, i.e., upx = upy = 0, are imposed in all material points in the boundary of the lower pin. Next, the vertical displacements of all material points in the boundary of upper pin are kinematically constrained to have the same vertical displacement. The penalty parameter Copyright c 2016 John Wiley & Sons, Ltd. Prepared using nmeauth.cls Int. J. Numer. Meth. Engng (2016) DOI: 10.1002/nme PHASE-FIELD MATERIAL POINT METHOD FOR BRITTLE FRACTURE (a) (b) (c) (d) 37 Figure 19. Phase-Field Material Point Method in L-Panel test for h = 10 mm, l0 = 10 mm and Gc = 8.9 · 10−1 kN/mm. (a) u = 25.60 mm, (b) u = 90 mm, (c) u = 278 mm and (d) u = 392 mm. was chosen to be a = 1000000. Finally, the displacement is monitored and controlled in the boundary of upper pin. To investigate the influence of stagger solution algorithm four cases are considered with constant displacement increment ∆u = 10−2 , ∆u = 5 · 10−3 mm, ∆u = 10−3 mm and ∆u = 5 · 10−4 mm. Table III presents the influence of stagger iterations Nstaggs and displacement increment ∆u on critical load and its corresponding displacement. In all cases tolu = 10−5 . Increasing the number of stagger iterations, the fidelity of the solution is improved. After the third stagger iteration the algorithm is converging to a value and additional iterations marginally affect the results at the cost of increased number of evaluation. As shown in Table III, i.e., ∆u = 10−3 mm and ∆u = 5 · 10−4 mm, when the displacement increment ∆u is sufficiently small then there is no need for additional stagger iterations. Copyright c 2016 John Wiley & Sons, Ltd. Prepared using nmeauth.cls Int. J. Numer. Meth. Engng (2016) DOI: 10.1002/nme 38 E. G. KAKOURIS, S. P. TRIANTAFYLLOU (a) (b) Figure 20. (a) Notched Plate with Hole. Geometry and boundary conditions. (b) Active cells in the beginning of the analysis. ∆u [mm] 10−2 5 · 10−3 10−3 5 · 10−4 Critical vertical displacement [mm] Critical load [kN] Nstaggs = 1 Nstaggs = 2 Nstaggs = 3 Nstaggs = 1 Nstaggs = 2 Nstaggs = 3 0.400 0.350 0.330 10.233 9.421 9.103 0.355 0.325 0.315 9.501 9.006 8.823 0.309 0.301 0.298 8.722 8.585 8.534 0.302 0.297 0.296 8.595 8.516 8.487 Table III. Influence of Stagger Solution Algorithm Fig. 21(a) presents the influence of displacement increment ∆u, for Nstaggs = 1. For large values of displacement increment the results are overestimated. However, as long as the displacement increment is decreased the results are converging. Whereas, Fig. 21(b) represents the influence of stagger iterations on the results for ∆u = 5 · 10−3 mm. The load-displacement curve is stabilized in third stagger iteration. Figs. 22 represent the evolution of phase field in four timesteps u = 0.28 mm, u = 0.35 mm, u = 0.96 mm and u = 1.20 mm respectively. Whereas, Figs. 23 present the evolution of hydrostatic stress for the same timesteps. Both Figs. 22 and 23 are referred to the solution obtained by stagger Nstaggs = 1 iteration with constant displacement increment ∆u = 10−3 mm. the analysis time for this simulation was approximately 37hrs. The crack paths of both the proposed method and the experimental data from [22] have good agreement and they are presented in Figs. 24. Copyright c 2016 John Wiley & Sons, Ltd. Prepared using nmeauth.cls Int. J. Numer. Meth. Engng (2016) DOI: 10.1002/nme 39 PHASE-FIELD MATERIAL POINT METHOD FOR BRITTLE FRACTURE (a) (b) Figure 21. Load displacement response dependence on (a) Displacement increment size for 1 stagger iteration (Nstaggs = 1) (b) Number of stager iterations (∆u = 5 · 10−3 mm) (a) (b) (c) (d) Figure 22. Phase-Field in Notched Plate with Hole test for (a) u = 0.28 mm, (b) u = 0.35 mm, (c) u = 0.96 mm and (c) u = 1.20 mm respectively. Displacement increment ∆u = 10−3 and Stagger with 1 iteration. 6. CONCLUSIONS In this work, a Material Point Method (MPM) for the simulation of crack propagation pertinent to brittle fracture is formulated. The crack geometry and underlying brittle fracture mechanics have been considered on the basis of a phase-field formulation, appropriately adapted to be introduced within the MPM framework. Using this approach, the need for algorithmically tracking the crack path is alleviated thus reducing the underlying computational complexity. The deformable domain is approximated using a set of material points that are allowed to move within a fixed Eulerian mesh. Copyright c 2016 John Wiley & Sons, Ltd. Prepared using nmeauth.cls Int. J. Numer. Meth. Engng (2016) DOI: 10.1002/nme 40 E. G. KAKOURIS, S. P. TRIANTAFYLLOU (a) (b) (c) (d) Figure 23. Hydrostatic stress in Notched Plate with Hole test for (a) u = 0.28 mm, (b) u = 0.35 mm, (c) u = 0.96 mm and (d) u = 1.20 mm respectively. Displacement increment ∆u = 10−3 and Stagger with 1 iteration. Material Points with cp < 0.08 have been removed. (a) (b) Figure 24. Comparison of crack path between experimental [22] (a) and simulation data (b). Material Points with cp < 0.08 have been removed. Fusing MPM with phase field modelling results in a coupled system of governing equations, namely the equilibrium and phase field evolution equations. Coupling is achieved through the definition of a history field that tracks the evolution of the tensile part of the elastic energy density. The latter is derived on the basis of a spectral decomposition of the elastic strain tensor. Copyright c 2016 John Wiley & Sons, Ltd. Prepared using nmeauth.cls Int. J. Numer. Meth. Engng (2016) DOI: 10.1002/nme 41 PHASE-FIELD MATERIAL POINT METHOD FOR BRITTLE FRACTURE The resulting coupled field equations is numerically treated using a staggered solution scheme. In this, the phase field evolution equations are solved for constant displacement fields and a prediction for the phase field distribution is derived which is then used to iterate for the updated displacement field through the equilibrium equations. A set of benchmark applications is presented to verify the proposed scheme whereas validation is also performed through experimental data available in the literature. From these, the accuracy of the method is established both in predicting the equilibrium paths as well as representing the actual geometry of the crack paths. The method succeeds in providing realistic crack path geometries as the background mesh is reused within each computational cycle, thus avoiding mesh distortion errors pertinent to the standard finite element method. However, this comes at an increased computational cost for the implicit MPM implementation considered in this work, as re-factorization of the underlying stiffness matrix is required in each computational cycle. In the proposed scheme discontinuities are naturally created, since material points are naturally separated. This is not the case in the classical Phase-Field Finite Element Method where Gauss points are always located in Gauss positions. From the preceding analysis, Phase-Field Material Point Method is established as a promising computational tool, able to address very challenging tasks in computational mechanics. Due to its significant advantages in treating dynamics and large displacement problems, the introduced PF-MPM can be further extended to treat impact-fracture and ductile fracture problems. Furthermore, due to its particle in cell formulation, the method can be coupled with both mesh-based and particle based methods within a multiscale setting. APPENDIX A MATERIAL POINT PENALTY METHOD FOR IMPOSING KINEMATICAL CONSTRAINTS In the following, we consider the case of the two dimensional problem shown in Fig. 25. This involves 2 active cells with 6 active grid nodes totalling 12 degrees of freedom for the background Eulerian grid. Eight material points points are considered. Kinematical constraints need to be imposed on both displacement components of the material point p, namely up = {upx , upy }. As described in Section Copyright c 2016 John Wiley & Sons, Ltd. Prepared using nmeauth.cls Int. J. Numer. Meth. Engng (2016) DOI: 10.1002/nme 42 E. G. KAKOURIS, S. P. TRIANTAFYLLOU Figure 25. Imposition of constraints with Penalty Method. 3.3, in the material point method the displacement components of a material point are interpolated to the corresponding nodes of the active cell the material point resides. For the case of material point p shown in Fig. 25, the interpolation expressions for both the horizontal and vertical components of the displacement vector are expressed in the following form     upx = N1 (xp )u3 + N2 (xp )u5 + N3 (xp )u11 + N4 (xp )u9 (54)    upy = N1 (xp )u4 + N2 (xp )u6 + N3 (xp )u12 + N4 (xp )u10 Following the standard procedure for the Penalty Method [32], the following matrices are formulated  0 0 B=  0 0 N1 (xp ) 0 N2 (xp ) 0 0 0 N4 (xp ) 0 N3 (xp ) 0 0 N1 (xp ) 0 N2 (xp ) 0 0 0 N4 (xp ) 0 N3 (xp )     (55) and V = {upx , upy } (56) where B is an (Nconstr × Ndof s ) coefficient matrix whereas V is an (1 × Nconstr ) coefficient vector. Nconstr and Ndof s denotes to the total number of imposed constraints and total number of active unconstrained degree of freedom of the whole structure. Therefore, the incremental external forces as well as the incremental internal forces are modified according to equations (57) and (58) below ∆Fext := ∆Fext + αBT VT Copyright c 2016 John Wiley & Sons, Ltd. Prepared using nmeauth.cls (57) Int. J. Numer. Meth. Engng (2016) DOI: 10.1002/nme 43 PHASE-FIELD MATERIAL POINT METHOD FOR BRITTLE FRACTURE ∆Fint := ∆Fint + αBT B∆u (58) where α is the Penalty Method parameter. If a direct solver is utilized then the global stiffness matrix of the structure is redefined according to equation (59) below Ku := Ku + αBT B (59) ACKNOWLEDGEMENT The research described in this paper has been financed by the University of Nottingham through the Dean of Engineering Prize, a scheme for pump priming support for early career academic staff. The authors are grateful to the University of Nottingham for access to its high performance computing facility. REFERENCES 1. Wu J, McAuliffe C, Waisman H, Deodatis G. Stochastic analysis of polymer composites rupture at large deformations modeled by a phase field method. Computer Methods in Applied Mechanics and Engineering 2016; 312(1):596–634. 2. Ferté G, Massin P, Moës N. 3D crack propagation with cohesive elements in the extended finite element method. Computer Methods in Applied Mechanics and Engineering 2016; 300(1):347–374. 3. Konuk I, Gürtner A, Yu S. A Cohesive Element Framework for Dynamic Ice-Structure Interaction Problems Part II: Implementation. ASME 2009 28th International Conference on Ocean, Offshore and Arctic Engineering, American Society of Mechanical Engineers, 2009; 185–193. 4. Murakami S. Continuum Damage Mechanics: A Continuum Mechanics Approach to the Analysis of Damage and Fracture, vol. 185. Springer Science & Business Media, 2012. 5. Sukumar N, Dolbow JE, Moës N. Extended finite element method in computational fracture mechanics: a retrospective examination. International Journal of Fracture 2015; 196(1):189–206. 6. Duddu R. Numerical modeling of corrosion pit propagation using the combined extended finite element and level set method. Computational Mechanics 2014; 54(3):613–627. 7. Moës N, Belytschko T. Extended finite element method for cohesive crack growth. Engineering Fracture Mechanics 2002; 69(7):813–833. Copyright c 2016 John Wiley & Sons, Ltd. Prepared using nmeauth.cls Int. J. Numer. Meth. Engng (2016) DOI: 10.1002/nme 44 E. G. KAKOURIS, S. P. TRIANTAFYLLOU 8. Fries TP, Belytschko T. The intrinsic XFEM: a method for arbitrary discontinuities without additional unknowns. International Journal for Numerical Methods in Engineering 2006; 68(13):1358–1385. 9. Agathos K, Chatzi E, Bordas SP. Stable 3D extended finite elements with higher order enrichment for accurate non planar fracture. Computer Methods in Applied Mechanics and Engineering 2016; 306(1):19–46. 10. Radovitzky R, Seagraves A, Tupek M, Noels L. A scalable 3D fracture and fragmentation algorithm based on a hybrid, discontinuous Galerkin, cohesive element method. Computer Methods in Applied Mechanics and Engineering 2011; 200(1):326–344. 11. Snozzi L, Molinari JF. A cohesive element model for mixed mode loading with frictional contact capability. International Journal for Numerical Methods in Engineering 2013; 93(5):510–526. 12. Remmers J, de Borst R, Needleman A. A cohesive segments method for the simulation of crack growth. Computational mechanics 2003; 31(1-2):69–77. 13. Kaczmarczyk Ł, Nezhad MM, Pearce C. Three-dimensional brittle fracture: configurational-force-driven crack propagation. International Journal for Numerical Methods in Engineering 2014; 97(7):531–550. 14. Bourdin B, Francfort GA, Marigo JJ. The Variational Approach to Fracture. Journal of Elasticity 2008; 91(1-3):5– 148. 15. Fix GJ. Phase field methods for free boundary problems 1982; . 16. Borden MJ, Verhoosel CV, Scott MA, Hughes TJ, Landis CM. A phase-field description of dynamic brittle fracture. Computer Methods in Applied Mechanics and Engineering 2012; 217(1):77–95. 17. Francfort GA, Marigo JJ. Revisiting brittle fracture as an energy minimization problem. Journal of the Mechanics and Physics of Solids 1998; 46(8):1319–1342. 18. Sun CT, Jin ZH. Fracture Mechanics. Academic Press, 2012. 19. Zienkiewicz OC, Taylor RL. The Finite Element Method: Solid Mechanics, vol. 2. Butterworth-Heinemann, 2000. 20. Miehe C, Welschinger F, Hofacker M. Thermodynamically consistent phase-field models of fracture: Variational principles and multi-field FE implementations. International Journal for Numerical Methods in Engineering 2010; 83(10):1273–1311. 21. Verhoosel CV, Borst R. A phase-field model for cohesive fracture. International Journal for Numerical Methods in Engineering 2013; 96(1):43–62. 22. Ambati M, Gerasimov T, De Lorenzis L. A review on phase-field models of brittle fracture and a new fast hybrid formulation. Computational Mechanics 2015; 55(2):383–405. 23. Borden MJ, Hughes TJ, Landis CM, Anvari A, Lee IJ. A phase-field formulation for fracture in ductile materials: Finite deformation balance law derivation, plastic degradation, and stress triaxiality effects. Computer Methods in Applied Mechanics and Engineering 2016; 312(1):130–166. 24. Li T, Marigo JJ, Guilbaud D, Potapov S. Gradient damage modeling of brittle fracture in an explicit dynamics context. International Journal for Numerical Methods in Engineering 2016; 108(11):1381–1405. Copyright c 2016 John Wiley & Sons, Ltd. Prepared using nmeauth.cls Int. J. Numer. Meth. Engng (2016) DOI: 10.1002/nme 45 PHASE-FIELD MATERIAL POINT METHOD FOR BRITTLE FRACTURE 25. Rangarajan R, Chiaramonte MM, Hunsweck MJ, Shen Y, Lew AJ. Simulating curvilinear crack propagation in two dimensions with universal meshes. International Journal for Numerical Methods in Engineering 2015; 102(34):632–670. 26. Scholtès L, Donzé FV. Modelling progressive failure in fractured rock masses using a 3D discrete element method. International Journal of Rock Mechanics and Mining Sciences 2012; 52:18–30. 27. Batra R, Zhang G. Search algorithm, and simulation of elastodynamic crack propagation by modified smoothed particle hydrodynamics (MSPH) method. Computational Mechanics 2007; 40(3):531–546. 28. Bobaru F, Hu W. The Meaning, Selection, and Use of the Peridynamic Horizon and its Relation to Crack Branching in Brittle Materials. International Journal of Fracture 2012; 176(2):215–222. 29. Amiri F, Milln D, Shen Y, Rabczuk T, Arroyo M. Phase-field modeling of fracture in linear thin shells. Theoretical and Applied Fracture Mechanics 2014; 69:102–109. 30. Nguyen VP, Rabczuk T, Bordas S, Duflot M. Meshless methods: A review and computer implementation aspects. Mathematics and Computers in Simulation 2008; 79(3):763–813. 31. Sulsky D, Chen Z, Schreyer HL. A particle method for history-dependent materials. Computer Methods in Applied Mechanics and Engineering 1994; 118(1-2):179–196. 32. Bathe KJ. Finite element procedures. Klaus-Jurgen Bathe, 2006. 33. Li J, Hamamoto Y, Liu Y, Zhang X. Sloshing impact simulation with material point method and its experimental validations. Computers & Fluids 2014; 103:86–99. 34. Jassim I, Stolle D, Vermeer P. Two-phase dynamic analysis by material point method. International Journal for Numerical and Analytical Methods in Geomechanics 2013; 37(15):2502–2522. 35. Abe K, Soga K, Bandara S. Material Point Method for Coupled Hydromechanical Problems. Journal of Geotechnical and Geoenvironmental Engineering 2013; 140(3). 36. Bandara S, Soga K. Coupling of soil deformation and pore fluid flow using material point method. Computers and Geotechnics 2015; 63:199–214. 37. Schreyer H, Sulsky D, Zhou SJ. Modeling delamination as a strong discontinuity with the material point method. Computer Methods in Applied Mechanics and Engineering 2002; 191(23-24):2483–2507. 38. Sulsky D, Kaul A. Implicit dynamics in the material-point method. Computer Methods in Applied Mechanics and Engineering 2004; 193(1214):1137–1170. 39. Nairn JA. Material Point Method Calculations with Explicit Cracks. Computer Modeling in Engineering and Sciences 2003; 4(6):649–664. 40. Daphalapurkar NP, Lu H, Coker D, Komanduri R. Simulation of dynamic crack growth using the generalized interpolation material point (GIMP) method. International Journal of Fracture 2007; 143(1):79–102. 41. Bardenhagen SG, Nairn JA, Lu H. Simulation of dynamic fracture with the Material Point Method using a mixed J-integral and cohesive law approach. International Journal of Fracture 2011; 170(1):49–66. 42. Yang P, Gan Y, Zhang X, Chen Z, Qi W, Liu P. Improved decohesion modeling with the material point method for simulating crack evolution. International Journal of Fracture 2014; 186(1):177–184. Copyright c 2016 John Wiley & Sons, Ltd. Prepared using nmeauth.cls Int. J. Numer. Meth. Engng (2016) DOI: 10.1002/nme 46 E. G. KAKOURIS, S. P. TRIANTAFYLLOU 43. Sanchez J, Schreyer H, Sulsky D, Wallstedt P. Solving quasi-static equations with the material-point method. International Journal for Numerical Methods in Engineering 2015; 103(1):60–78. 44. Homel MA, Herbold EB. Field-gradient partitioning for fracture and frictional contact in the material point method. International Journal for Numerical Methods in Engineering 2016; doi:10.1002/nme.5317. 45. Ambati M, Gerasimov T, De Lorenzis L. Phase-field modeling of ductile fracture. Computational Mechanics 2015; 55(5):1017–1040. 46. Miehe C, Aldakheel F, Raina A. Phase field modeling of ductile fracture at finite strains: A variational gradientextended plasticity-damage theory. International Journal of Plasticity 2016; 84:1–32. 47. Trädegard A, Nilsson F, Östlund S. FEM-remeshing technique applied to crack growth problems. Computer Methods in Applied Mechanics and Engineering 1998; 160(1-2):115–131. 48. Miehe C, Hofacker M, Welschinger F. A phase field model for rate-independent crack propagation: Robust algorithmic implementation based on operator splits. Computer Methods in Applied Mechanics and Engineering 2010; 199(45-48):2765–2778. 49. Borden MJ, Hughes TJ, Landis CM, Verhoosel CV. A higher-order phase-field model for brittle fracture: Formulation and analysis within the isogeometric analysis framework. Computer Methods in Applied Mechanics and Engineering 2014; 273:100–118. 50. Li B, Peco C, Milln D, Arias I, Arroyo M. Phase-field modeling and simulation of fracture in brittle materials with strongly anisotropic surface energy. International Journal for Numerical Methods in Engineering 2015; 102(34):711–727. 51. Kuhn C, Schlter A, Mller R. On degradation functions in phase field fracture models. Computational Materials Science 2015; 108, Part B:374–384. 52. Ambrosio L, Tortorelli VM. On the approximation of free discontinuity problems 1992; 6(1):105–123. 53. Braides A. Approximation of Free-Discontinuity Problems. Springer, Berlin, 1998. 54. Steffen M, Kirby RM, Berzins M. Analysis and reduction of quadrature errors in the material point method (MPM). International journal for numerical methods in engineering 2008; 76(6):922. 55. Simo J, Miehe C. Associative coupled thermoplasticity at finite strains: Formulation, numerical analysis and implementation. Computer Methods in Applied Mechanics and Engineering 1992; 98(1):41–104. 56. Gerasimov T, Lorenzis LD. A line search assisted monolithic approach for phase-field computing of brittle fracture. Computer Methods in Applied Mechanics and Engineering 2016; 312:276 – 303, doi:http://dx.doi.org/10.1016/j. cma.2015.12.017. 57. Guilkey JE, Weiss JA. Implicit time integration for the material point method: Quantitative and algorithmic comparisons with the finite element method. International Journal for Numerical Methods in Engineering 2003; 57(9):1323–1338. 58. Cortis M, Coombs W, Augarde C. Implicit essential boundaries in the material point method. Cardiff University, 2016. Copyright c 2016 John Wiley & Sons, Ltd. Prepared using nmeauth.cls Int. J. Numer. Meth. Engng (2016) DOI: 10.1002/nme PHASE-FIELD MATERIAL POINT METHOD FOR BRITTLE FRACTURE 47 59. Winkler B. Traglastuntersuchungen von unbewehrten und bewehrten Betonstrukturen auf der Grundlage eines objektiven Werkstoffgesetzes für Beton. Dissertation. University of Innsbruck, Austria, 2001. Copyright c 2016 John Wiley & Sons, Ltd. Prepared using nmeauth.cls Int. J. Numer. Meth. Engng (2016) DOI: 10.1002/nme