Block structured adaptive mesh refinement and strong form elasticity approach to phase field fracture with applications to delamination, crack branching and crack deflection
AA survey of fracture, delamination, crack branching and crackdeflection in composites using a novel framework for phase fieldfracture
Vinamra Agrawal a, ∗ , Brandon Runnels b a Department of Aerospace Engineering, Auburn University, Auburn, AL, USA b Department of Mechanical and Aerospace Engineering, University of Colorado, Colorado Springs, CO USA
Abstract
Fracture is a ubiquitous phenomenon in most composite engineering structures, and is oftenthe responsible mechanism for catastrophic failure. Over the past several decades, manyapproaches have emerged to model and predict crack failure. The phase field method forfracture uses a surrogate damage field to model crack propagation, eliminating the arduousneed for explicit crack meshing. In this work a novel numerical framework is proposedfor implementing hybrid phase field fracture in heterogeneous composite materials. Theproposed method is based on the “reflux-free” method for solving, in strong form, theequations of linear elasticity on a block-structured adaptive mesh refinement (BSAMR) mesh.The use of BSAMR enables highly efficient and scalable regridding, facilitates the use oftemporal subcycling for explicit time integration, and allows for ultra-high refinement atcrack boundaries with minimal computational cost. The method is applied to a variety ofsimple composite structures: laminates, wavy interfaces, and circular inclusions. In eachcase a non-dimensionalized parameter study is performed to identify regions of behavior,varying both the geometry of the problem and the relative fracture energy release rate. Inthe laminate and wavy interface cases, regions of delamination and fracture correspond tosimple analytical predictions. For the circular inclusions, the modulus ratio of the inclusionis varied as well as the delamination energy release rate and the problem geometry. In thiscase, a wide variety of behaviors was observed, including deflection, splitting, delamination,and pure fracture.
Keywords:
Fracture mechanics, phase field fracture, delamination, composites,block-structured adaptive mesh refinement, strong form elasticity
1. Introduction
Crack propagation is known to be the prominent cause of failure of composite engineeringstructures ranging from catastrophic structural failure to aging and/or catastrophic failure ofsolid rocket propellants [1–3]. Composite engineering materials - including structural (suchas nanolayered Cu-Nb or carbon-reinforced plastic) and energetic (such as AP+HTPB orAluminum+PVDF solid propellant) - can exhibit an even wider range of fracture behavior ∗ Corresponding author: [email protected] a r X i v : . [ c ond - m a t . m e s - h a ll ] F e b ue to the combination of delamination, fracture, and crack deflection resulting from materialheterogeneity. Enhanced numerical methods are required in order to predict fracture behaviorin such materials, and to create materials that are failure-resistant by design.Understanding and modeling crack propagation in materials has been a topic of significantstudy since the early 70s. A variational approach to fracture [4] has enabled modeling cracksthrough phase field methods for ease of computational implementation and efficiency [5–7].This approach introduces a length scale parameter ξ for the crack, which is then used toregularize the underlying energy functional [4]. It has been shown that the variational methodis consistent with linear elastic fracture mechanics (LEFM). The variational approach has beenextended to study brittle fracture [6, 8, 9], dynamic fracture [10–13], anisotropic materials[14, 15], interface strength [16], crack nucleation [17–19], ductile formulation [20–22] andfatigue induced crack propagation [23]. Efforts have also been directed towards understandingthe role of damage degradation function and the damage energy penalty term on the solutions[24–29].While crack propagation in homogeneous materials have been extensively studied, hetero-geneous materials have received relatively less attention. In heterogeneous materials, crackcan undergo deflection, repeated arrest and nucleation depending on the material architectureand constituent properties. This is especially important when considering (micro-) crackpropagation in the micro and mesoscales. Micro-cracks can interact with material interfaces,pores and grain boundaries leading to complex crack paths [30]. Recently phase field methodshave been used to study crack propagation in heterogeneous material from the context ofeffective fracture toughness [31, 32]. Other works have focused on predicting crack paths andfailure in composites [33–35] and bio-inspired materials [36–38]. Insights into crack pathsand failure mechanisms can lead to designing failure resistant materials.The interaction of crack with material interfaces can lead to either crack arrest, punchingthrough or redirection along the interface leading to delamination. This behavior is theresult of energetically competing mechanisms involving mechanical work and the cost ofcreating new surfaces. Recently, phase field formulation of fracture mechanics has beenmodified to incorporate cohesive zone elements to account for interfacial strength [39–41].While these methods have successfully modeled interface failure [39, 41–44], there is a needfor a systematic understanding of the competing mechanisms driving crack behavior. Thisis especially important for designing conventional fiber-reinforced composites and novelmaterials such as bio-inspired sutured interface materials [45–47].From a numerical perspective, phase field approach allows simulating complicated fractureprocesses such as nucleation, branching and merging without requiring explicit trackingof crack field. This offers significant advantages over discrete fracture formulations wherediscontinuities have to be explicitly (using classic finite element method) or implicitly (usingextended finite element method) tracked. In general, phase field fracture problems are solvedusing a staggered approach where displacement and crack fields are solved alternately [48], orin a monolithic fashion [49, 50] where both fields are solved simultaneously. While monolithicsolvers are more efficient, staggered approaches have been shown to be more robust [48] andare therefore widely adopted. These strategies typically use finite element methods (FEM)with an adaptive mesh refinement (AMR) strategies [51]. Recently, phase field fractureproblems have also been studied through fast Fourier transforms (FFT) [52, 53] for improvedmemory efficiency. 2n this work, we present a reflux-free finite difference (FD) based block structured AMR(BSAMR) methodology to solve the phase field fracture problem in strong form. We use thereflux-free solver to solve the elastic equilibrium equation implicitly, coupled with explicit-time integration for solving crack evolution equation. FEM is ubiquitous in solid mechanicsproblems due to a) their ability to handle complex geometry, and b) the weak formulationthat removes the need for higher order derivatives. However, in this work, we focus oncrack propagation at a microstructural scale with simplified domain geometries. We treatthe internal material boundaries as diffused interfaces, making the advantages of FEM lessrelevant for our case. FD offers a simple, fast and robust way to solve equilibrium equationsin this case.We use BASMR to refine the mesh near the crack and material boundaries. BSAMRorganizes the mesh into levels, each comprising of grid cells of the same size [54]. It treatseach level independently avoiding excess communication overhead between levels. This allowsBSAMR to be highly scalable and well suited for massively parallel implementations. Weuse a reflux-free method to treat the coarse-fine mesh boundaries. Finally, we use temporalsub-stepping in explicit time integration of the crack field within different levels to avoidnumerical instabilities.This paper is organized as follows. In section 2, we provide a brief overview of the hybridphase field fracture including equilibrium equations. We describe the FD based reflux-freeBSAMR solver in section 3, followed by three examples in section 4: a crack impinging ona planar surface, delamination/fracture of a wavy interface, and crack interaction with acircular inclusion.
2. Hybrid model for phase field fracture
We use the hybrid formulation of the phase field fracture proposed by Ambati et al . [51]is a general reference for this section. The smooth crack field is represented by c ( x ) , suchthat c = 0 inside the crack and 1 outside. The hybrid formulation combines the isotropic [5]and anisotropic formulations [9] for improved computational efficiency. The isotropic energyfunctional [5] is given by F i ( u , c ) = (cid:90) Ω (cid:0) c + η (cid:1) Ψ ( ε ) dx + (cid:90) Ω G c (cid:104) (1 − c ) ξ + ξ |∇ c | (cid:105) dx, (1)where η (cid:28) G c is the material fracture toughness, u is thedisplacement field, ε = sym grad u is the strain field and Ψ is the elastic strain energydensity. For a linear elastic isotropic material with Lam´e constants λ and µ , the elastic strainenergy density is given by Ψ ( ε ) = 12 λ (tr ε ) + µ tr (cid:0) ε (cid:1) . (2)The Euler Lagrange equilibrium equations result indiv σ = , σ = (cid:0) c + η (cid:1) ∂ Ψ ∂ ε , c H − G c (cid:20) ξ ∆ c + 1 − s ξ (cid:21) , (3)3here H := max τ ∈ [0 ,t ] Ψ ( ε ( x, t )).The anisotropic model [9] assumes an additive decomposition of the elastic strain energyΨ = Ψ +0 + Ψ − , which uses the spectral decomposition of the strain tensor ε = (cid:80) di =1 ε i ˆ n i ⊗ ˆ n i ,where ε i are principal strains, ˆ n i are principal directions and d is the dimension. Finally,Ψ ± = 12 λ (tr ε ± ) + µ tr (cid:0) ε ± (cid:1) , ε ± = d (cid:88) i =1 ( ε i ) ± ˆ n i ⊗ ˆ n i (4)The regularized anisotropic energy functional is given by F a ( u , c ) = (cid:90) Ω (cid:2)(cid:0) c + η (cid:1) Ψ +0 + Ψ − (cid:3) dx + (cid:90) Ω G c (cid:34) (1 − c ) ξ + ξ |∇ c | (cid:35) dx. (5)The Euler Lagrange equilibrium equation for anisotropic formulation are given bydiv σ = , σ = (cid:0) c + η (cid:1) ∂ Ψ +0 ∂ ε + ∂ Ψ − ∂ ε , c H + − G c (cid:20) ξ ∆ c + 1 − s ξ (cid:21) , (6)where H + := max τ ∈ [0 ,t ] Ψ +0 ( ε ( x, t )).In this work, the hybrid model is used to leverage the computational efficiency of theisotropic model, while still retaining a crack driven by tensile component of the elastic energy.The equilibrium equations of the hybrid model are given bydiv σ = , σ = (cid:0) c + η (cid:1) ∂ Ψ ∂ ε , (7a)0 = 2 c H + − G c (cid:20) ξ ∆ c + 1 − s ξ (cid:21) , (7b) ∀ x : Ψ +0 < Ψ − ⇒ c = 1 , (7c)where (7c) is added to prevent crack interpenetration. In this work, we replace the crackequilibrium equation (7b) with a Ginzberg-Landau type evolution law [55]:˙ c = − M (cid:20) c H + − G c (cid:18) ξ ∆ c + 1 − s ξ (cid:19)(cid:21) , (8)where M ≥
3. Reflux-free Finite Difference Solver
Fracture mechanics problems, by definition, contain highly localized features (cracks)in an otherwise relatively uninteresting ambient medium. As a result, strategic meshing isrequired to provide adequate resolution at the features of interest without incurring excessivecomputational cost due to the large problem size. The majority of previous phase fieldfracture work uses the finite element method to solve the equations of elasticity on a meshthat is either refined a priori [10] or using adaptive mesh refinement [51, 56].4 igure 1:
Block structured adaptive mesh refinement (with seven total levels) of crack propagation into andaround a circular inclusion Mesh refinement is used to track the crack front as well as to resolve the materialinterface.
The phase field method for fracture, unlike many other methods such as the cohesivezone model, does not require any explicit meshing scheme. Indeed, phase field methods oftenperform optimally when implemented on a regular grid, rather than a specifically tailoredunstructured mesh. However, they are often implemented using such a mesh anyway, in orderto accommodate the requisite FEM elasticity solver. As a result, such numerical schemescarry the computational burden of FEM based adaptive mesh refinement, while offering fewbenefits to the implementation of the phase field fracture model itself.In this work, we employ an alternative numerical implementation that is optimal for boththe implicit and explicit calculations in phase field fracture. Block-structured AMR (BSAMR)is a specific variety of adaptive mesh refinement that partitions the mesh into individuallevels, where each level is associated with a specific degree of refinement (Figure 1). Eachlevel is broken up into individual cubic regions called patches, where each patch contains datastored as a regular, structured grid. Each patch, in turn, is stored using a space-filling Z-curvescheme in order to maximize data locality. Over the course of temporal integration, patchesare evolved independently from one another, then eventually synced up using averagingbetween levels. The multi-level structure enables efficient temporal subcycling, allowingfiner levels to evolve with a different timestep than the coarser levels. Though not exploredspecifically in this work, temporal subcycling is highly advantageous for fourth-order fractureevolution models, allowing for the timestep to decrease commensurate with the increase ofspatial discretization, per the CFL condition. In the BSAMR method, mesh refinement andcoarsening occur at regular intervals (see section below) and trigger the regeneration of gridsusing the Berger-Rigoutsos algorithm [57]. This scheme offers myriad advantages, includingan ability to scale easily to very large platforms, the ready portability to heterogeneousarchitectures such as GPUs, and the localization of data for the most numerically intensiveoperations.Despite its many advantages, however, BSAMR has seen very limited application to solidmechanics and elasticity problems. This is primarily due to the inherent difficulties in solvingthe equations of static elasticity on a BSAMR mesh. In this work, we build on previouscontributions by the authors toward the development of a strong form elasticity solver [58]for numerical methods using BSAMR. In this method, the equations of elasticity along withboundary conditions are posed as a single partitioned operator D : Ω × R n → R n on problem5omain Ω, defined as D ( x )( u ) = div (cid:2) C ( x ) (cid:0) grad u − ε ( x ) (cid:1)(cid:3) x ∈ int(Ω) u x ∈ ∂ Ω (Dirichlet boundary) C ( x ) (cid:0) grad u − ε ( x ) (cid:1) n ( x ) x ∈ ∂ Ω (Traction boundary)(grad u ) n ( x ) x ∈ ∂ Ω (Neumann boundary) (9)where C is the spatially varying modulus tensor, ε is a spatially varying eigenstrain tensor,and n is the outward-facing normal vector. The standard equations of elasticity, along withboundary conditions, are thus expressible as D ( x )( u ) = r ( x ) , (10)where the right-hand side term r is defined piecewise as r ( x ) = b ( x ) x ∈ int(Ω) u ( x ) x ∈ ∂ Ω t ( x ) x ∈ ∂ Ω q ( x ) x ∈ ∂ Ω , (11)with b the body force, u the prescribed displacement, t the prescribed traction, and q theprescribed gradient of u normal to the boundary.Operator (9) is implemented numerically in strong form, where C ( x ) describes thespatial variation of C in accordance with (7a). The traditional difficulties of a strong formimplementation are avoided by the use of C++ templating and operator overloading, allowingthe highly efficient evaluation of fourth-order matrix multiplication operations. Crucially, theuse of strong form allows for the definition of level-consistent numerical operators D n , suchthat the constants associated with D n may be interpolated and restricted between levels in amanner consistent with the numerical discretization. That is, if I is the injective interpolationoperator and R the surjective restriction operator, the implementation of D n , D n +1 must besuch that the following are maintained: D n = R ◦ D n +1 ◦ I (12)for all n . (Note that the reciprocal relation is generally not true.) It is worthwhile tomention, at this point, a pitfall that can arise when solving near-singular elasticity problems,such as phase field fracture, with highly localized boundaries: It is sometimes the case thatoperators . . . , D − , D , D , . . . are initialized separately, such that (12) does not hold. Such aninconsistency can cause substantial convergence issues, even though the separate initializationmay actually be more accurate. To avoid such problems, coarse operators should always beinitialized by D n − = R D in all the areas of overlap, to ensure consistency.BSAMR is highly conducive to the use of multigrid for solving (10), as the refinementlevels posess the same structure as coarsening levels in the geometric multigrid method. D , D , . . . correspond to the discretized operator on the unrefined level, the first level ofrefinement, and so on. D − corresponds to the first coarsened multigrid level, D − to the6econd coarsening level, and so on. The operators R, I for AMR levels are the same as thoseused for the multigrid interpolation and restriction operations.The method of successive over-relaxation is used for the multigrid smoothing operation,defined here to be S = ( D D + ωD L ) − (cid:16) ω r [ ωD U + ( ω − D D ] x (cid:17) , (13)where D D , D L , D U are the diagonal, lower, and uper components of D . It was determinedthat a relaxation factor of ω = 0 .
08 produced optimal convergence for the problems discussedhere. For the bottom-level linear solve, the smoothing operator S was used instead of the morecommonly used stabilized biconjugate gradient method, due to the near-singular behavior ofthe operator.The implementation of the restriction operator across the coarse/fine boundary is oneof the key difficulties in solving elasticity using BSAMR. This work builds on the previousdevelopment of the “reflux-free” algorithm, which eliminates the need for an explicit coarse/finetreatment by using additional ghost nodes at the boundary. For a full discussion of thereflux-free method, the reader is referred to [58].
4. Examples
We use an in-house code (Alamo) to solve crack equations associated with the hybridformulation of phase field fracture to systematically study the interaction of crack withmaterial interfaces in two dimensions. We study three cases under Mode-I loading: interactionof crack with angled planar interface, interaction of crack with wavy interface, and crackpropagation through a maze of spherical inclusions. To focus on the interplay betweendriving forces governing delamination and fracture, we consider simplified problems that arerepresentative examples of crack/interface interaction. Each material is represented by aphase variable φ i , while the crack field is represented by c . For all cases, we choose a squaredomain x, y ∈ [ − . , . × [ − . , .
01] with 64 ×
64 base mesh, with six levels of BSAMRrefinement.Mesh regridding occurred every 10 Level-0 timesteps. Refinement occurred until thefollowing criterion was met:max( |∇ c | , |∇ φ | ) | ∆ x n | > . ∀ n ∈ (0 , N ) , ∀ x n ∈ Ω n . (14)In the above, n indicates mesh level with a maximum of N possible levels, x n the positionvector within each level’s domain Ω n , and ∆ x n the mesh resolution vector for level n . Allmaterials (matrix and inclusion) are linear elastic isotropic with same Lam´e parameters λ = 121 . × , µ = 80 . × . Materials have same fracture energy release rate G c = G cf = 2 . × , crack length scale ξ = 10 − and mobility M = 10 − . We choose η = 5 × − for numerical stability. In all the examples outlined below, the followingboundary conditions are used. u = 0 , y = − . , ∀ x ∈ [ − . , . u = u , y = 0 . , ∀ x ∈ [ − . , . u · ˆ e = 0 , x = ± . , ∀ y ∈ [ − . , . σ ˆ e ) · ˆ e = 0 , x = ± . , ∀ y ∈ [ − . , . . (15)7inally, the material interface in the examples presented below is diffused with a length scale b = 5 × − to ensure a large enough b/ξ ratio to avoid re-scaling the interfacial delaminationenergy release rate [59]. To compute the overall fracture energy release rate, we computedthe total crack volume in the domain and the total elastic energy ( c + η ) Ψ ( ε ) at every step.We then smoothed the data using a Savitzky-Golay filter, before differentiating it to computethe incremental energy release rate. G cd / G cf ( d e g ) A BFEC DAnalyticFractureDelaminationTransition
Figure 2:
Laminate: (
Left ) Mode I loading of a material with twin interfaces angled at θ from thevertical axis. Colored lines indicate sample crack patterns: the blue line indicates pure fracture, green indicates maximum delamination, and yellow indicates the transition region. ( Right ) Paramater surveyof delamination and fracture for an angled planar interface. Callouts A-F correspond to visualizationsof crack/delamination patterns in Figure 3, and colors correspond to the convention in the Right figure.Contours are a visual aid to assist in identifying regions of behavior.
In this set of simulations we consider a crack impinging an angled laminate inclusion. Theinterface has a delamination energy release rate G cd , which is varied from 0 . G cf to G cf . Wevary the orientation of the laminate inclusion from 0 ◦ (vertical) to 90 ◦ (horizontal) (Figure 2Left). We investigate the relationship between the non-dimensional parameters G cd /G cf andthe orientation θ to determine the regions in which delamination and fracture are dominant.In the case of ideal crack propagation, it is possible to estimate these two regions using asimple analytic calculation. If the crack tip has just encountered the interface, the transitionoccurs when the differential energy released by fracture dE f equals that of delamination dE d .The energies are related to the crack length as dE f ∼ G cf da f dE d ∼ G cd da d (16)where da f and da d are the differential area increases corresponding to fracture and delami-nation, respectively. Since relationship between the areas is da f = da d sin θ, the transitionbetween delamination and fracture clearly occurs when G cd G cf = sin θ. (17)8 a) (b) (c)(d) (e) (f) Figure 3:
Laminate: Visualization of crack field corresponding to cases A-F highlighted in Figure 2 right.The six representative cases demonstrate pure fracture , delamination & transition behavior of the crack.The laminate interface is overlaid in black for visualization purposes. We emphasize that this relation does not account for the additional energy required to redirectthe crack in the event of delamination, which is apparently nontrivial. Therefore, we expectdelamination to be somewhat less likely than predicted by (17).A total of 110 simulations were performed with varying angle θ and G cd /G cf ratio (Fig-ure 2). Each simulation is colored based on the observed behavior, with green correspondingto delamination and blue to fracture. Delamination generally correspond to a clean crackpropagation along the interface (Figure 3a). However, it was observed when the laminate tiltangle exceeded 45 ◦ , that the delamination cracks generally redirected into the material as afracture (Figure 3c). This is a result of edge effects, and these events are still categorizedas pure delamination. Pure fracture events were identified as cracks that were generallyunaffected by the presence of the interface (Figure 3d) or where the effect of the interfacewas minimal (Figure 3f). The pure fracture behavior was generally identical to classic ModeI fracture.There were a substantial number of “transition” results, which are all indicated in yellow.These are cracks that were neither purely delamination nor purely fracture, but containedsome combination of both. Some cracks delaminated for a very short period but then splitoff to fracture the material (Figure 3b). Crack branching was very common, in which themode was predominantly fracture combined with small delaminations at the crack-interfaceintersection. There was a continuum between transition and fracture, and the categorizationof a pattern as transition vs fracture is somewhat subjective. The primary bifurcation occursat the line between the transition and the delamination points. The match between theobserved fracture/delamination bifurcation point and the highly simplified, idealized analyticsolution appears to be reasonable.Figure 4 shows the evolution of the fracture energy release rate G as the crack moves9 .4 1.6 1.8 2.0 2.2 2.4 2.6 2.8 Crack volume
1e 50.0260.0280.0300.0320.0340.0360.0380.0400.042 E n e r g y r e l e a s e r a t e G cd / G cf = 0.70, = 18.00 (Delamination - Fig 3A) G cd / G cf = 0.90, = 18.00 (Branch - Fig 3B) G cd / G cf = 0.55, = 54.00 (Delamination - Fig 3C) G cd / G cf = 0.95, = 36.00 (Fracture - Fig 3 D) G cd / G cf = 0.55, = 72.00 (Branch - Fig 3E) G cd / G cf = 0.80, = 81.00 (Fracture - Fig 3F) Figure 4:
Laminate: Plot depicting the evolution of fracture energy release rate as crack propagates throughthe domain. through the domain for different cases. In the case of delamination ( G cd /G cf = 0 . θ = 54 ◦ ,Figure 3c), the energy release rate first shows a sharp dip followed by an increase and thenfinally a dip. The first dip corresponds to when the crack first encounters an interfaceresulting in crack redirection and delamination. As the crack progresses traveling in anenergetically low trajectory, the stress increases and so does G . The crack finally reflectsoff the boundary into a fracture mode of failure resulting in a decline in G . In the case offracture ( G cd /G cf = 0 . θ = 54 ◦ , Figure 3f), the crack stays mostly constant except fortwo dips corresponding to the point is encounters the two interfaces. Finally, for the case ofbranching ( G cd /G cf = 0 . θ = 72 ◦ , Figure 3e), the energy release rate shows bigger andwider dips. After encountering the first interface, the energy release rate dips and the cracksplits off. Due to higher energies associated with crack deflection, the delaminated sectionof the crack eventually stops. This causes the energy release rate to stabilize. The processrepeats itself when the second interface is encountered. Delamination vs fracture behavior in non-planar interfaces is of great importance in anyapplication involving joined materials. Classic examples include material systems producedby thermal spraying [60], textile ceramic composites [61], and bio-inspired sutured interfaces[45–47] In this section we consider the behavior of a crack under Mode I loading along asinusoidally varying interface with G cd ≤ G cf . Once again, the interface was diffused with alength scale b = 5 × − . In such a configuration, there is an obvious competition betweenthe lower energetic cost of delamination vs the lower amount of surface area resulting fromfracture. As with the previous example there are three anticipated possibilities: pure fracture,pure delamination, and transitionary behavior comprising elements of both (Figure 5 Left).Following the above, we can estimate the transition point between delamination andfracture via a simple by-hand calculation. The relationship between the fracture and delami-nation areas is slightly more complex. Assuming a homogenized interface and that the crack10 .5 0.6 0.7 0.8 0.9 G cd / G cf A / L FA E DB CAnalyticFractureDelaminationTransition
Figure 5:
Wavy interface: (
Left ) Schematic of the sinusoidal interface problem. Colored lines indicatesample crack paths for ease of visualization. Colors correspond to delamination (green), transition (orange)and fracture (blue). (
Right ) Paramater survey of delamination/fracture for a sinusoidal interface identifyingdelamination, fracture and transition zones. Callouts A-F correspond to visualization of crack paths presentedin Figure 6. Simplified analytic prediction is overlaid for comparison. (Compare: [44] Figure 10) propagates uniformly, then the area elements are related by the arc length formula as da f = (cid:32) π (cid:90) π/ (cid:114) π A L sin ( x ) dx (cid:33) da d = 2 π E (cid:16) − π A L (cid:17) da d , (18)where E is the compete elliptic integral of the second kind, A is the amplitude of the sinudoid,and L is the period. Therefore, we find the relationship between the Gc ratio and the ratio A/L to be G cd G cf = π E ( − π A /L ) . (19)We emphasize that this highly simplified calculation does not account for curvature effects,for mode II effects, or for the energetic cost associated with crack kinking. Therefore weexpect this transition point to be, in many cases, artificially high.We considered 110 configurations with A/L ranging from 0 to 0 . G cd /G cf ratiosfrom 0.5 to 1.0. As before, we identify three zones: pure fracture , delamination , and transition (Figure 5 right). As before, there were significant transitionary cases where thecrack starts to delaminate, but then splits off to fracture the material. As expected theactual transition happens much lower than the transition point obtained from the analyticalpredictions. Figure 6 shows actual crack path for six different cases corresponding to separatezones marked in Figure 5 demonstrating transition from branching to delamination to fracturefor different amplitudes and G cd /G cf ratios.Figure 7 shows the evolution of fracture energy release rate G as the crack propagatesleading to delamination, fracture and branching. As expected, for pure fracture, G stabilizesand remains constant for most of the crack path, and finally dropping off. The final drop off in G happens because the crack reaches the edge of the domain. For delamination, the G shows11 a) (b) (c)(d) (e) (f) Figure 6:
Wavy interface: Visualization of a selection of crack behavior. The wavy interface is marked ingray. Each figure corresponds to a point on Figure 5. The outline boxes denote branching , delamination ,and fracture wave-like behavior as the crack follows the interface leading to delamination. For branching,the G curve exhibits fluctuations upon encountering the interface, and then proceeds with aconstant value. These results are consistent with [44] who studied kinking vs delaminationbehavior for cracks using cohesive zone modeling with different amplitude G cd /G cf ratios. Crack deflection and redirection in heterogeneous materials has been a topic studiedthrough experimental and numerical methods. In this section, we consider an idealizedproblem of a linear elastic isotropic circular inclusion of radius R and modulus E i embeddedin a linear elastic isotropic matrix of modulus E m . We choose the same fracture energy G cf of both materials but vary the delamination fracture energy of the interface G cd . Finally, wevary the offset h of the center of the inclusion from the geometric center of the square domain.We study the role of non-dimensional parameters E i /E m , G cd /G cf and h/R on potential crackpaths. As with the previous two cases, crack behavior may be categorized as:) penetrationof the inclusion ( inclusion fracture ), delamination of the interface ( delamination ) or acombination ( branching ). Here, an additional possibility is that the crack not interact withthe inclusion at all, perhaps as a result of deflection ( matrix fracture ).Figure 8 presents a summary plot for different crack behavior for varying offset ratios h/R , fracture energy ratios G cd /G cf , and modulus ratios E i /E m . As expected, the observedbehavior is highly complex, and often accompanied by transition modes. In general, for low G cd /G cf values, delamination (depicted by green) is the predominant mode of failure evenfor high offset ratios. For high offset ratios, the crack path transitions from delamination topure matrix fracture as G cd /G cf increases. For zero offset ratio, the crack path transitionsfrom a branching behavior to pure inclusion fracture as G cd /G cf increases. This is to be12 .4 1.6 1.8 2.0 2.2 2.4 2.6 2.8 Crack volume
1e 50.0260.0280.0300.0320.0340.0360.0380.0400.042 E n e r g y r e l e a s e r a t e G cd / G cf = 0.55, A / L = 0.35 (Delamination - Fig6B) G cd / G cf = 0.85, A / L = 0.45 (Fracture - Fig6C) G cd / G cf = 0.75, A / L = 0.25 (Branch - Fig6D) Figure 7:
Wavy Interface: Plot depicting the evolution of fracture energy release rate as crack propagatesthrough the domain. G cd / G cf h / R A BC DA BA B
Figure 8:
Inclusion: (
Left ) Schematic of the problem and observed crack paths. (
Right ) Plot depictingdifferent outcomes for different scenarios. Each point corresponds to a specific G cd /G cf ratio and offset. Theinner circles correspond to a modulus ratio E i /E m of 0 .
5, the middle circles to a modulus ratio of 1, and theouter circles to a modulus ratio of 2 .
0. Colors correspond to matrix fracture , pure inclusion fracture , delamination , and branching . expected since the energetic cost of delamination increases as G cd /G cf ratio increases. Fora fixed value of G cd /G cf ratio, the crack path changes from pure inclusion fracture to purematrix fracture as offset ratio increases. Once again, this is to be expected because as theoffset ratio increases, the cost of crack deflection also increases.The modulus ratio E i /E m plays a critical role in governing crack behavior. Changingmodulus ratio introduces complex stress distribution around the inclusion causing cracksto deflect. Increasing the modulus ratio from 0 . . . G cd /G cf and offset ratios. Figure 9 shows actual crack paths for13 i / E m = . E i / E m = . E i / E m = . (a) G cd G cf = 0 . hR = 0 . (b) G cd G cf = 1 . hR = 0 . (c) G cd G cf = 0 . hR = 0 (d) G cd G cf = 0 . hR = 1 . Figure 9:
Inclusion: Visualization of a selection of crack behavior. In the grid of images, each column ofthree images corresponds to a point on Figure 8. The top image corresponds to a modulus ratio of 0 . . . inclusion fracture , branching , delamination , matrix fracture ).
12 cases with different modulus ratios. Column (a) of Figure 9 shows one such case where thecrack deflects “up” to cause inclusion fracture for E i /E m = 0 .
5, and deflects “down” to causematrix fracture for E i /E m = 2 .
0. In such cases, the energetic cost associated with deflectionis lower than the overall cost of propagation. The effect of modulus ratio is also evident incolumn (d) of Figure 9 where the crack deflects “up” for E i /E m = 0 . . E i /E m = 1 . E i /E m = 1 .
0, there is no stress concentration around the inclusion. This, in conjunction withhigh offset and higher G cd /G cf , results in a classic Mode-I type matrix fracture.Figure 10 shows the evolution of the fracture energy release rate G as the crack propagatesthrough the domain. A baseline response of flat G vs crack curve is obtained (orange) for14 .4 1.6 1.8 2.0 2.2 2.4 Crack volume
1e 50.0240.0260.0280.0300.0320.0340.0360.0380.040 E n e r g y r e l e a s e r a t e G cd / G cf = 0.5, h / R = 0.4 E i / E m = 1.0 (Delamination) G cd / G cf = 1.0, h / R = 0.4 E i / E m = 1.0 (Inclusion Fracture) G cd / G cf = 0.75, h / R = 0.8 E i / E m = 0.5 (Inclusion Fracture) G cd / G cf = 0.75, h / R = 0.8 E i / E m = 1.0 (Delamination) G cd / G cf = 0.75, h / R = 0.8 E i / E m = 2.0 (Matrix Fracture) Figure 10:
Inclusion: Evolution of crack energy release rate as it propagates through the domain. G cd /G cf = 1 . h/R = 0 . E i /E m = 1 . G curves are outcome of complex stress states and energy distributionaround the crack and the inclusion. For cases other than the baseline, the G curve exhibits asharp dip corresponding to the crack being in the vicinity of the interface. For delaminationcases (blue and red), G first shows a dip and then shows an increase above the baselineMode-I type response. For the matrix fracture case (purple), the crack deflects down to avoidthe inclusion resulting in a G curve that mostly stays below the baseline response.
5. Conclusions
The contributions of this work are twofold: First, we have presented a novel frameworkfor efficiently solving phase field fracture problems using block-structured adaptive meshrefinement. Second, we have demonstrated the efficacy of this framework by investigating alarge set of problems in composite fracture mechanics.The proposed method has demonstrated efficiency and accuracy, but is not without anumber of limitations. Like most phase field fracture codes, solving the equations of elasticityrequires a high number of solver iterations due to the poor condition number of the elasticoperator for the degraded material. It should also be pointed out that the performance ofthe method generally decreases as the crack grows in length, as the entire crack length mustbe fully resolved at all times during the simulation. The use of BSAMR with a regular gridalso means that all non-rectangular geometry must be included in a diffuse-interface manner.Finally, the use of a rectangular grid limits the applicability of the method. While it has awide range of application to representative volume element (RVE) type simulations, it haslimited applicability to non-rectangular domains.Three types of composite structures were examined using the proposed method. A largesuite of simulations was completed for each type of structure. (We note that not all simulationsran to completion; however, in each case, those that did not had already progressed farenough to determine the crack propagation behavior. ) First, laminate-type structures wereinvestigated to determine the regions of crack behavior with respect to the delamination vs15racture energy release rate, and the angle of the interface. The results were compared withan analytic (albeit simplified) model, and reasonable agreement was determined. Second,the case of a sinusoidal interface was considered, with analysis similar to the above butwith respect to the amplitude of the interface rather than the angle. As with the first case,reasonable agreement was observed between the analytic estimation and the results of theparameter study. Finally, the case of a circular inclusion was considered, where the stiffnessof the conclusion varied from 1 / Acknowledgements
B Runnels was funded by the Office of Naval Research, grant number N00014-21-1-2113.The authors are grateful for the support of the Auburn University Easley (2020) Cluster forassistance with this work.
References [1] R. Hu, C. Prakash, V. Tomar, M. Harr, I. E. Gunduz, and C. Oskay, “Experimentally-validated mesoscale modeling of the coupled mechanical–thermal response of ap–htpbenergetic material under dynamic loading,”
International Journal of Fracture , vol. 203,no. 1-2, pp. 277–298, 2017.[2] P. Kadiresh and B. Sridhar, “Experimental evaluation and simulation on aging char-acteristics of aluminised ap-htpb composite solid propellant,”
Materials Science andTechnology , vol. 24, no. 4, pp. 406–412, 2008.[3] Z. Wang, H. Qiang, T. Wang, and G. Wang, “Tensile behaviors of thermal aged htpbpropellant at low temperatures under dynamic loading,”
Mechanics of Time-DependentMaterials , vol. 24, no. 2, pp. 141–159, 2020.[4] G. A. Francfort and J.-J. Marigo, “Revisiting brittle fracture as an energy minimizationproblem,”
Journal of the Mechanics and Physics of Solids , vol. 46, no. 8, pp. 1319–1342,1998.[5] B. Bourdin, G. A. Francfort, and J.-J. Marigo, “Numerical experiments in revisited brittlefracture,”
Journal of the Mechanics and Physics of Solids , vol. 48, no. 4, pp. 797–826,2000.[6] B. Bourdin, “Numerical implementation of the variational formulation for quasi-staticbrittle fracture,”
Interfaces and free boundaries , vol. 9, no. 3, pp. 411–430, 2007.[7] B. Bourdin, G. A. Francfort, and J.-J. Marigo, “The variational approach to fracture,”
Journal of elasticity , vol. 91, no. 1-3, pp. 5–148, 2008.[8] C. Kuhn and R. M¨uller, “A continuum phase field model for fracture,”
EngineeringFracture Mechanics , vol. 77, no. 18, pp. 3625–3634, 2010.169] C. Miehe, L.-M. Schaenzel, and H. Ulmer, “Phase field modeling of fracture in multi-physics problems. part i. balance of crack surface and failure criteria for brittle crackpropagation in thermo-elastic solids,”
Computer Methods in Applied Mechanics andEngineering , vol. 294, pp. 449–485, 2015.[10] M. J. Borden, C. V. Verhoosel, M. A. Scott, T. J. Hughes, and C. M. Landis, “A phase-field description of dynamic brittle fracture,”
Computer Methods in Applied Mechanicsand Engineering , vol. 217, pp. 77–95, 2012.[11] A. Schl¨uter, A. Willenb¨ucher, C. Kuhn, and R. M¨uller, “Phase field approximation ofdynamic brittle fracture,”
Computational Mechanics , vol. 54, no. 5, pp. 1141–1161, 2014.[12] A. Karma, D. A. Kessler, and H. Levine, “Phase-field model of mode iii dynamic fracture,”
Physical Review Letters , vol. 87, no. 4, p. 045501, 2001.[13] M. Hofacker and C. Miehe, “Continuum phase field modeling of dynamic fracture:variational principles and staggered fe implementation,”
International Journal of Fracture ,vol. 178, no. 1-2, pp. 113–129, 2012.[14] B. Li, C. Peco, D. Mill´an, I. Arias, and M. Arroyo, “Phase-field modeling and simulationof fracture in brittle materials with strongly anisotropic surface energy,”
InternationalJournal for Numerical Methods in Engineering , vol. 102, no. 3-4, pp. 711–727, 2015.[15] S. Teichtmeister, D. Kienle, F. Aldakheel, and M.-A. Keip, “Phase field modeling offracture in anisotropic brittle solids,”
International Journal of Non-Linear Mechanics ,vol. 97, pp. 1–21, 2017.[16] D. Schneider, E. Schoof, Y. Huang, M. Selzer, and B. Nestler, “Phase-field modelingof crack propagation in multiphase systems,”
Computer Methods in Applied Mechanicsand Engineering , vol. 312, pp. 186–195, 2016.[17] E. Tann´e, T. Li, B. Bourdin, J.-J. Marigo, and C. Maurini, “Crack nucleation invariational phase-field models of brittle fracture,”
Journal of the Mechanics and Physicsof Solids , vol. 110, pp. 80–99, 2018.[18] S. Brach, E. Tann´e, B. Bourdin, and K. Bhattacharya, “Phase-field study of cracknucleation and propagation in elastic–perfectly plastic bodies,”
Computer Methods inApplied Mechanics and Engineering , vol. 353, pp. 44–65, 2019.[19] A. Kumar, B. Bourdin, G. A. Francfort, and O. Lopez-Pamies, “Revisiting nucleation inthe phase-field approach to brittle fracture,”
Journal of the Mechanics and Physics ofSolids , p. 104027, 2020.[20] M. Ambati, T. Gerasimov, and L. De Lorenzis, “Phase-field modeling of ductile fracture,”
Computational Mechanics , vol. 55, no. 5, pp. 1017–1040, 2015.[21] C. Miehe, M. Hofacker, L.-M. Sch¨anzel, and F. Aldakheel, “Phase field modeling offracture in multi-physics problems. part ii. coupled brittle-to-ductile failure criteriaand crack propagation in thermo-elastic–plastic solids,”
Computer Methods in AppliedMechanics and Engineering , vol. 294, pp. 486–522, 2015.1722] C. Kuhn, T. Noll, and R. M¨uller, “On phase field modeling of ductile fracture,”
GAMM-Mitteilungen , vol. 39, no. 1, pp. 35–54, 2016.[23] A. Mesgarnejad, A. Imanian, and A. Karma, “Phase-field models for fatigue crackgrowth,”
Theoretical and Applied Fracture Mechanics , vol. 103, no. December 2018,p. 102282, 2019.[24] E. Lorentz and V. Godard, “Gradient damage models: Toward full-scale computations,”
Computer Methods in Applied Mechanics and Engineering , vol. 200, no. 21-22, pp. 1927–1944, 2011.[25] K. Pham, H. Amor, J.-J. Marigo, and C. Maurini, “Gradient damage models and theiruse to approximate brittle fracture,”
International Journal of Damage Mechanics , vol. 20,no. 4, pp. 618–652, 2011.[26] K. Pham, J.-J. Marigo, and C. Maurini, “The issues of the uniqueness and the stabilityof the homogeneous response in uniaxial tests with gradient damage models,”
Journalof the Mechanics and Physics of Solids , vol. 59, no. 6, pp. 1163–1190, 2011.[27] K. Pham and J.-J. Marigo, “From the onset of damage to rupture: construction ofresponses with damage localization for a general class of gradient damage models,”
Continuum Mechanics and Thermodynamics , vol. 25, no. 2-4, pp. 147–171, 2013.[28] C. Kuhn, A. Schl¨uter, and R. M¨uller, “On degradation functions in phase field fracturemodels,”
Computational Materials Science , vol. 108, pp. 374–384, 2015.[29] J.-Y. Wu and V. P. Nguyen, “A length scale insensitive phase-field damage model forbrittle fracture,”
Journal of the Mechanics and Physics of Solids , vol. 119, pp. 20–42,2018.[30] S. Kumar and W. A. Curtin, “Crack interaction with microstructure,”
Materials today ,vol. 10, no. 9, pp. 34–44, 2007.[31] M. Hossain, C.-J. Hsueh, B. Bourdin, and K. Bhattacharya, “Effective toughness ofheterogeneous media,”
Journal of the Mechanics and Physics of Solids , vol. 71, pp. 15–32,2014.[32] C. Hsueh, L. Avellar, B. Bourdin, G. Ravichandran, and K. Bhattacharya, “Stressfluctuation, crack renucleation and toughening in layered materials,”
Journal of theMechanics and Physics of Solids , vol. 120, pp. 68–78, 2018.[33] J. J. Espadas-Escalante, N. P. van Dijk, and P. Isaksson, “A phase-field model forstrength and fracture analyses of fiber-reinforced composites,”
Composites Science andTechnology , vol. 174, pp. 58–67, 2019.[34] P. Zhang, X. Hu, T. Q. Bui, and W. Yao, “Phase field modeling of fracture in fiberreinforced composite laminate,”
International Journal of Mechanical Sciences , vol. 161,p. 105008, 2019. 1835] W. Tan and E. Mart´ınez-Pa˜neda, “Phase field predictions of microscopic fracture andr-curve behaviour of fibre-reinforced composites,”
Composites Science and Technology ,p. 108539, 2020.[36] P. Murali, T. K. Bhandakkar, W. L. Cheah, M. H. Jhon, H. Gao, and R. Ahluwalia, “Roleof modulus mismatch on crack propagation and toughness enhancement in bioinspiredcomposites,”
Physical Review E , vol. 84, no. 1, p. 015102, 2011.[37] S. Khaderi, P. Murali, and R. Ahluwalia, “Failure and toughness of bio-inspired com-posites: Insights from phase field modelling,”
Computational materials science , vol. 95,pp. 1–7, 2014.[38] A. Singh, T. S. Sandhu, and S. Pal, “Interplay of various fracture mechanisms inbio-inspired staggered structure,”
Mechanics of Materials , vol. 139, p. 103215, 2019.[39] V. Carollo, J. Reinoso, and M. Paggi, “Modeling complex crack paths in ceramiclaminates: A novel variational framework combining the phase field method of fractureand the cohesive zone model,”
Journal of the European Ceramic Society , vol. 38, no. 8,pp. 2994–3003, 2018.[40] P. Tarafder, S. Dan, and S. Ghosh, “Finite deformation cohesive zone phase field modelfor crack propagation in multi-phase microstructures,”
Computational Mechanics , vol. 66,no. 3, pp. 723–743, 2020.[41] A. Quintanas-Corominas, A. Turon, J. Reinoso, E. Casoni, M. Paggi, and J. Mayugo,“A phase field approach enhanced with a cohesive zone model for modeling delaminationinduced by matrix cracking,”
Computer Methods in Applied Mechanics and Engineering ,vol. 358, p. 112618, 2020.[42] P. Roy, S. Deepu, A. Pathrikar, D. Roy, and J. Reddy, “Phase field based peridynamicsdamage model for delamination of composite structures,”
Composite Structures , vol. 180,pp. 972–993, 2017.[43] B. Dhas, M. Rahaman, K. Akella, D. Roy, J. Reddy, et al. , “A phase-field damagemodel for orthotropic materials and delamination in composites,”
Journal of AppliedMechanics , vol. 85, no. 1, 2018.[44] S. Sehr, S. Amidi, and M. R. Begley, “Interface delamination vs. bulk cracking alongwavy interfaces,”
Engineering Fracture Mechanics , vol. 206, pp. 64–74, 2019.[45] Y. Li, C. Ortiz, and M. C. Boyce, “A generalized mechanical model for suture interfacesof arbitrary geometry,”
Journal of the Mechanics and Physics of Solids , vol. 61, no. 4,pp. 1144–1167, 2013.[46] Y. Cao, W. Wang, J. Wang, and C. Zhang, “Experimental and numerical study ontensile failure behavior of bionic suture joints,”
Journal of the Mechanical Behavior ofBiomedical Materials , vol. 92, pp. 40–49, 2019.1947] Z. Liu, Z. Zhang, and R. O. Ritchie, “Interfacial toughening effect of suture structures,”
Acta Biomaterialia , vol. 102, pp. 75–82, 2020.[48] C. Miehe, M. Hofacker, and F. Welschinger, “A phase field model for rate-independentcrack propagation: Robust algorithmic implementation based on operator splits,”
Com-puter Methods in Applied Mechanics and Engineering , vol. 199, no. 45-48, pp. 2765–2778,2010.[49] P. Farrell and C. Maurini, “Linear and nonlinear solvers for variational phase-fieldmodels of brittle fracture,”
International Journal for Numerical Methods in Engineering ,vol. 109, no. 5, pp. 648–667, 2017.[50] D. Jodlbauer, U. Langer, and T. Wick, “Matrix-free multigrid solvers for phase-fieldfracture problems,”
Computer Methods in Applied Mechanics and Engineering , vol. 372,p. 113431, 2020.[51] M. Ambati, T. Gerasimov, and L. De Lorenzis, “A review on phase-field models of brittlefracture and a new fast hybrid formulation,”
Computational Mechanics , vol. 55, no. 2,pp. 383–405, 2015.[52] Y. Chen, D. Vasiukov, L. G´el´ebart, and C. H. Park, “A fft solver for variational phase-fieldmodeling of brittle fracture,”
Computer Methods in Applied Mechanics and Engineering ,vol. 349, pp. 167–190, 2019.[53] F. Ernesti, M. Schneider, and T. B¨ohlke, “Fast implicit solvers for phase-field fractureproblems on heterogeneous microstructures,”
Computer Methods in Applied Mechanicsand Engineering , vol. 363, p. 112793, 2020.[54] A. Dubey, A. Almgren, J. Bell, M. Berzins, S. Brandt, G. Bryan, P. Colella, D. Graves,M. Lijewski, F. L¨offler, et al. , “A survey of high level frameworks in block-structuredadaptive mesh refinement packages,”
Journal of Parallel and Distributed Computing ,vol. 74, no. 12, pp. 3217–3227, 2014.[55] C. Kuhn and R. M¨uller, “A phase field model for fracture,” in
PAMM: Proceedings inApplied Mathematics and Mechanics , vol. 8, pp. 10223–10224, Wiley Online Library,2008.[56] T. Heister, M. F. Wheeler, and T. Wick, “A primal-dual active set method and predictor-corrector mesh adaptivity for computing fracture propagation using a phase-field ap-proach,”
Computer Methods in Applied Mechanics and Engineering , vol. 290, pp. 466–495,2015.[57] M. Berger and I. Rigoutsos, “An algorithm for point clustering and grid generation,”
IEEE Transactions on Systems, Man, and Cybernetics , vol. 21, no. 5, pp. 1278–1286,1991.[58] B. Runnels, V. Agrawal, W. Zhang, and A. Almgren, “Massively parallel finite differenceelasticity using block-structured adaptive mesh refinement with a geometric multigridsolver,”
Journal of Computational Physics , p. 110065, 2020.2059] A. C. Hansen-D¨orr, R. de Borst, P. Hennig, and M. K¨astner, “Phase-field modellingof interface failure in brittle materials,”
Computer Methods in Applied Mechanics andEngineering , vol. 346, pp. 25–42, 2019.[60] C.-H. Hsueh, J. A. Haynes, M. J. Lance, P. F. Becher, M. K. Ferber, E. R. Fuller, S. A.Langer, W. C. Carter, and W. R. Cannon, “Effects of interface roughness on residualstresses in thermal barrier coatings,”
Journal of the American Ceramic Society , vol. 82,no. 4, pp. 1073–1075, 1999.[61] M. Blacklock, J. H. Shaw, F. W. Zok, and B. N. Cox, “Virtual specimens for analyzingstrain distributions in textile ceramic composites,”