A two-scale generalized finite element method for interaction and coalescence of multiple crack surfaces

Engineering Fracture Mechanics 163 (2016) 274–302

Contents lists available at ScienceDirect

Engineering Fracture Mechanics journal homepage: www.elsevier.com/locate/engfracmech

A two-scale generalized finite element method for interaction and coalescence of multiple crack surfaces P. O’Hara a,⇑, C.A. Duarte b, T. Eason c a

Universal Technology Corporation, 1270 North Fairfield Rd, Dayton, OH 45432, USA Dept. of Civil and Environmental Engineering, University of Illinois, Newmark Laboratory, 205 North Mathews Avenue, Urbana, IL 61801, USA c Structural Sciences Center, Air Force Research Laboratory, WPAFB, OH 45433, USA b

a r t i c l e

i n f o

Article history: Received 2 March 2016 Received in revised form 13 June 2016 Accepted 15 June 2016 Available online 17 July 2016 Keywords: G/XFEM Multi-scale methods Enriched finite element methods Crack interaction Crack coalescence

a b s t r a c t This paper presents the application of a two-scale generalized finite element method (GFEM) which allows for static fracture analyses as well as fatigue crack propagation simulations involving the interaction of multiple crack surfaces on fixed, coarse finite element (FE) meshes. The approach is based on the use of numerically-generated enrichment functions computed on-the-fly through the use of locally-defined boundary value problems (BVPs) in the regions of existing mechanically-short cracks. The two-scale GFEM approach is verified against analytical reference solutions as well as alternative numerical approaches for crack interaction problems, including the coalescence of multiple crack surfaces. The numerical examples demonstrate the ability of the proposed approach to deliver accurate results even in scenarios involving multiple, interacting discontinuities contained within a single computational element. The proposed approach is also applied to a crack shielding/crack arrest problem involving two propagating crack surfaces in a representative panel model similar in complexity to that which may be of interest to the aerospace community. Ó 2016 Elsevier Ltd. All rights reserved.

1. Introduction The pursuit of accurate fatigue life predictions for structural components has been an active area of research in the civil, mechanical and aerospace communities for decades. From a numerical modeling standpoint, fatigue propagation simulations have been broadly performed for many years using standard FE approaches [1–5]. In a standard FE approach, the crack surface itself is modeled via meshes with double nodes to define the discontinuity and degenerate quarter-point elements in the crack tip region to model the strong gradients and singularity developing in the stress field in this region. Conventional FE approaches to fracture problems often lead to large problem sizes requiring significant amounts of computational resources [6,7]. There has been much work in the development of adaptive meshing techniques [8,9,7,10–14] which are aimed at the generation of highly graded meshes in crack tip regions, often employing crack front template meshes to help resolve the singularities. The use of highly graded meshes in general may lead to conditioning issues in the resulting system of equations. Additionally, crack front template meshes which utilize brick elements can become very difficult to construct in the case of evolving three-dimensional crack surfaces [11] and user intervention may be required to generate a suitable mesh [15]. With adaptive meshing approaches, care must be taken to ensure a well-structured mesh in the crack tip regions ⇑ Corresponding author. E-mail address: [email protected] (P. O’Hara). http://dx.doi.org/10.1016/j.engfracmech.2016.06.009 0013-7944/Ó 2016 Elsevier Ltd. All rights reserved.

P. O’Hara et al. / Engineering Fracture Mechanics 163 (2016) 274–302

275

Nomenclature a BVP C

crack length boundary value problem parameter in the Paris-Erdogan equation da crack growth rate in length/cycle dN E Young’s Modulus GFEM generalized finite element method GFEMgl generalized finite element method with global-local enrichment functions H1 ðxa Þ first order Hilbert space defined on xa K I ; K II ; K III mode I; II, and III stress intensity factors K eq equivalent stress intensity factor Lai enrichment function defined on xa m exponent in the Paris-Erdogan equation N number of load cycles R load ratio r initial crack radius SIFs stress intensity factors u displacement field uhp finite element approximation of u x point in space XFEM extended finite element method a index of node in the finite element mesh Damax maximum crack front increment DK stress intensity factor range DK eqmax maximum equivalent stress intensity factor range m Poisson’s ratio h crack front kinking angle r01 maximum principal stress rh ; rz ; shz components of the stress tensor w crack front twisting angle ua finite element partition of unity function generalized finite element shape function /ai va local approximation space on xa X solid, finite element domain xa set of elements connected to node a

throughout the course of a simulation so as to avoid, for instance, elements with poor aspect ratios which may adversely impact the resulting solution quality [13]. There are many existing numerical techniques which are aimed at avoiding the aforementioned meshing requirements and computational costs associated with conventional FE approaches to fracture analyses, while still maintaining sufficient levels of accuracy. The generalized/extended finite element method (G/XFEM) [16–21] is one such technique, which has been actively developed over the past ten years or so. As opposed to standard FE approaches which require the generation of a volume mesh which ‘fits’ the crack surface at every crack step [7,3,12,10,11], the G/XFEM is based upon the notion of using enrichment functions to accurately model the crack surface as it propagates throughout the course of a simulation. In this manner the G/XFEM is able to alleviate some of the potentially cumbersome modeling considerations, and associated computational costs often encountered with a standard FE approach. To further complicate the matter, in a high cycle fatigue framework, much of the fatigue life is spent in the crack initiation phase. In this phase multiple, diffuse discontinuity surfaces propagate and link up with one another in the formation of a single, dominant flaw, which eventually cause the ultimate fatigue failure of a component. A key element of a numerical tool appropriate for use in such a framework is the ability to simulate the interactions and eventual coalescence of multiple crack surfaces. The problem of crack coalescence has previously been investigated in the standard G/XFEM context in [22,23]. In the latter it was noted that there are issues with the enrichment strategy used in the XFEM when crack fronts approach one another due to the use of an implicit crack surface representation. Both approaches relied upon available analytical enrichment functions, defined separately for each crack under consideration, up until the point of coalescence. A potential limitation of the standard G/XFEM approach is the a priori requirement of enrichment functions which can accurately model the physics of the problem at hand. Unfortunately, accurate closed-form enrichment functions

276

P. O’Hara et al. / Engineering Fracture Mechanics 163 (2016) 274–302

are often unavailable, except for in a few specific scenarios. This is the case for the problem under consideration, with multiple discontinuities in close enough proximity to interact with one another, closed-form enrichment functions are in general not available. To address the lack of a priori available enrichment bases, a more general form of GFEM, that which utilizes global-local enrichment functions (GFEM gl ) [24,25] has been proposed. In the GFEMgl framework, the a priori requirement of a known enrichment basis is removed, and an additional locally-defined BVP is solved in a region of interest so as to numerically generate an appropriate enrichment function to model the relevant physics of a particular problem. The GFEMgl is not the only multi-scale approach which relies upon a ‘fine scale’ BVP whose solution is used to improve coarse-scale solution accuracy. Similar approaches include the multi-scale FEM put forth by Hou and Wu [26]; as well as the variational multi-scale enrichment technique of Oskay [27]. While not an exhaustive list, other multi-scale techniques which don’t take the same approach as the GFEMgl , but could be, or in some instances have been applied to the class of problems considered in this paper include the multi-scale projection methods of Loehnert et al. [28–30]; the spectral overlay method of Belytschko et al. [31]; the superposition FEM (s-fem) of Fish et al. [32,33]; the multi-scale FE approaches of Krause et al. [34,35]; concurrent multi-level FE approaches proposed by Ghosh et al. [36,37] and the combined s-fem/XFEM of Lee et al. [38]. The current work focuses on an application of the GFEM utilizing numerically generated enrichment functions to scenarios involving multiple, interacting crack surfaces. The proposed approach is able to explicitly resolve the interaction between multiple crack surfaces on fixed, coarse FE meshes solely through the use of scale-bridging enrichment functions. Static fracture scenarios involving shielding/amplification effects of micro-cracks in the regions of macro-crack tips, as well as interaction and coalescence of multiple propagating crack surfaces are investigated. The remainder of the paper is structured as follows. A problem description along with the governing equations and a brief presentation of the crack propagation criteria is presented in the next section. Brief overviews of the GFEM and GFEMgl are presented in Sections 3 and 4, respectively. A series of numerical examples are provided in Section 5 to illustrate the effectiveness of the proposed approach. Finally, the main conclusions from the work are summarized.

2. Problem formulation 2.1. Governing equations Consider a 3D domain, X, with external boundary @ X. Portions of @ X are subjected to tractions (cyclic in the case of fatigue crack propagation) and there are potentially multiple existing discontinuity surfaces in X as illustrated in Fig. 1. The material comprising X is assumed to be homogeneous, isotropic and linear elastic. The external boundary @ X is subdivided into Cu and Cr in a manner that Cu [ Cr ¼ @ X and Cu \ Cr ¼ ;. Dirichlet and Neumann boundary conditions are assumed to be prescribed on Cu and Cr , respectively.

Fig. 1. Typical 3D LEFM boundary value problem.

277

P. O’Hara et al. / Engineering Fracture Mechanics 163 (2016) 274–302

Equilibrium for the associated problem is defined by

$  r þ f ðxÞ ¼ 0

ð1Þ

where r is the Cauchy stress tensor and f is a prescribed body force vector. The assumption of a linear elastic material yields the following constitutive equations

rðxÞ ¼ ktrfeðxÞgI þ 2leðxÞ

ð2Þ

where k and l are two independent Lame elastic constants, tr feg is the trace of the small-strain tensor and I is the identity tensor. Additionally assuming small deformations, the small-strain tensor is given by the symmetric component of the deformation gradient as

eðxÞ ¼

 1
$u þ $uT : 2

ð3Þ

The following boundary conditions are enforced on the appropriate portions of @ X:

rðxÞnðxÞ ¼
tðx; tÞ on Cr

ð4Þ


ðxÞ on Cu uðxÞ ¼ u

ð5Þ

rðxÞnðxÞ ¼ 0 on C

c

ð6Þ


ðxÞ are prescribed displacements. The boundary conditions prescribed on C indiwhere
tðx; tÞ are prescribed tractions and u cate traction-free crack surfaces. In the context of finite elements we require the weak form of (1), which can be posed as: c

e 1 ðXÞ such that Find u 2 H Bðu; v Þ ¼ Lðv Þ 8

ð7Þ

v 2 H10

in which Bð; Þ and LðÞ are the bilinear and linear forms, respectively, and are given by

Z Bðu; v Þ ¼ rðuÞ : eðv ÞdX Z X Z
t  v dC: Lðv Þ ¼ f  v dX þ

ð8Þ

Cr

X

We define the set of kinematically admissible displacement fields as

n o e 1 ðXÞ ¼ ujuðxÞ 2 H1 ðXÞ; uðxÞ ¼ u
ðxÞ on Cu H

ð9Þ

where H1 is the first order Hilbert space. The space of kinematically admissible virtual displacement fields is defined as the subset of functions in H1 ðXÞ which satisfy a homogeneous Dirichlet boundary condition on Cu , i.e.

H10 ðXÞ ¼

n

v jv ðxÞ 2 H1 ðXÞ; v ðxÞ ¼ 0

o on Cu :

ð10Þ

The above is the standard weak formulation of a small-strain 3D linear elasticity problem suitable for a typical finite element formulation. Subsequently, (7) will be discretized in the standard FE manner using generalized finite element shape functions (c.f. Section 3). 2.2. Crack propagation considerations 2.2.1. Crack trajectory The discussion on the crack propagation criteria used in this work is left intentionally brief. For additional details the reader is referred to [39–42]. In a fully 3D setting, prediction of crack growth trajectory requires the computation of both a kinking angle ðhÞ and twisting angle ðwÞ at the crack front. In this work Schöllmann’s criterion [39], based on the assumption that crack growth develops perpendicularly to the maximum principal stress r01 , is used for computation of both deflection angles, as a function of the stress intensity factors (SIFs) at the crack front. The value of r01 is computed from the near field stresses as

r01 ¼

rh þ rz 2

þ

1 2

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðrh  rz Þ2 þ 4s2hz

ð11Þ

where rh ; rz and shz are the components of the stress tensor. The stress components take into account contributions of all three fracture Modes, and are computed in terms of cylindrical coordinates as

K



  h

  3h

I ffi 3 cos þ cos rh ¼ pffiffiffiffiffiffiffiffi 2 2 4 2pr

     K II h 3h  pffiffiffiffiffiffiffiffiffi 3 sin þ 3 sin 2 2 4 2pr

278

P. O’Hara et al. / Engineering Fracture Mechanics 163 (2016) 274–302

  h

K

III ffi cos shz ¼ pffiffiffiffiffiffiffiffi 2 2pr



 

 

8m h h  K II sin : rz ¼ pffiffiffiffiffiffiffiffiffi K I cos 2 2 4 2pr

ð12Þ

The mixed-modality is then accounted for in terms of an equivalent SIF, K eq , computed at each crack front vertex as

K eq

2 sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi3           2 1 h 4 h 3 h h 3 h K I cos2  K II sin þ K I cos2  K II sin ¼ cos þ 4K 2III 5: 2 2 2 2 2 2 2 2

ð13Þ

The kinking angle, h, is computed so as to satisfy

@ r01 ¼ 0 and @h

@ 2 r01 @h2

<0

ð14Þ

whereas the twisting angle, w, is computed as

w¼

1 arctan 2



 2shz ðhÞ : rh ðhÞ  rz ðhÞ

ð15Þ

2.2.2. Crack growth rate For fatigue crack growth simulations aimed at modeling the stable growth of a crack subjected to cyclic loads of constant amplitude, there exists significant amounts of experimental data suggesting that the crack growth rate may be related to the SIF range at the crack front. In this work the Paris-Erdogan equation is used to model the crack growth rate [43] as a function of the variation in the local SIF values as

da ¼ C ðDK Þm dN

ð16Þ

da where dN is the fatigue crack growth rate per loading cycle and DK is the stress intensity factor range during one cycle. C and m are material parameters. In this work, the SIF values at the crack front are computed from the FE solution using the Cut-off Function Method as may be found in the literature [41,42]. Due to the fully 3D nature of the problems considered, DK eqmax is used in (16) to compute the number of loading cycles corresponding to a given crack advancement. For a given crack step i, an incremental equation for fatigue crack growth in a linear elastic material is adopted. With the maximum crack front increment Damax set at the beginning of the given crack step, the current number of loading cycles is estimated in an incremental form of (16) as

Damax Ni ¼ N i1 þ  m C DK eqmax

ð17Þ

where N i and N i1 are the number of loading cycles at the current and previous crack steps, respectively. In a general 3D mixed-mode crack simulation, DK varies along a crack front. The crack growth increments at a given crack front vertex must be scaled so as to subject each vertex to the same number of loading cycles, and allow (16) to retain its validity. For a particular vertex j, the advancement is scaled as

 

m DK eqj m Damax Daj ¼ C DK eqj  m ¼ Damax DK eqmax C DK eqmax

ð18Þ

where DK eqj is the SIF range for vertex j and DK eqmax is the maximum SIF range seen at a crack front. In the case of multiple crack surfaces DK eqmax is taken as the overall maximum value taking into account all of the crack surfaces present in the analysis. 2.2.3. Crack coalescence considerations Coalescence of adjacent crack surfaces is handled in the same manner as that described in detail in [22]. The basic approach is that at a crack propagation step the minimum distance between two adjacent crack surfaces is computed and compared to a defined tolerance distance. If the minimum distance is less than the tolerance distance, the crack surfaces are joined. A set of control points are passed to the crack surface re-mesher corresponding now to the crack front vertices for both surfaces, aside from those that fall within the coalesced portion of the newly joined crack surface. The use of a tolerance distance amounts to the assumption of a zone ahead of the crack fronts with a reduced loadcarrying capacity. This assumption is the same as that made in [44,45], which assert that the plastic zones developing ahead of the crack fronts exhibit a substantial reduction in load-carrying capacity. In the current work the tolerance distance is prescribed as in [22,46], but the computation of such an effective plastic zone size based upon local SIFs is the focus of future

P. O’Hara et al. / Engineering Fracture Mechanics 163 (2016) 274–302

279

min work. In the current work the tolerance size is taken as D2a  0:05 in which a is the crack length (radius) and Dmin is the minimum separation distance between crack front vertices of adjacent crack surfaces. Just after the coalescence point, the joined crack surfaces often exhibit a highly non-convex region along the front. The use of relatively small crack front increments along with a moving least squares smoothing of crack front vertex locations has proven sufficient to avoid self-intersection of the propagating crack surface. Further details can be found in [22].

3. Generalized finite element method: an overview The GFEM [16] is a partition-of-unity (PoU) based FE method [47–51]. The GFEM is similar to the standard, Galerkin FE method [52], but offers increased flexibility in terms of shape function construction. Considering a particular node, a, in the FE mesh, a GFEM shape function is defined as

/ai ðxÞ ¼ Na ðxÞLai ðxÞ

ð19Þ

where there is no implied summation on a. The first component, N a ðxÞ is the PoU function, simply implying that the function sums to one at any point in the domain, i.e. N X Na ðxÞ ¼ 1;

8 x 2 X:

ð20Þ

a¼1

Within the GFEM context, the PoU functions are provided by standard, Langrangian FE shape functions, which inherently possess this property. The second component of the GFEM shape function, Lai ðxÞ, are the enrichment functions. The enrichment functions form a basis of local approximation spaces, va , which are a priori known to approximate well the local physics of the given problem. A standard application for the GFEM, and that which is the focus of this work, is linear elastic fracture mechanics (LEFM) [17], in which enrichment functions in the form of discontinuous Heaviside functions as well as singular crack tip expansion functions are used to accurately and efficiently model the displacement/stress contours associated with 3D cracks in an elastic medium. The use of GFEM for LEFM applications will be further discussed in subsequent sections. The resulting GFEM shape function obtains C 0 -continuity from the compactly supported PoU functions and also retains the approximation potential of the enrichment basis due to the PoU functions ability to exactly represent the enrichment functions in the enriched solution space. The GFEM sits atop a robust mathematical foundation [49,48,16,53,54]. As a result, an error bound has been derived for the global approximation error within the GFEM. The error bound is formally posed as Theorem. Let u 2 H1 ðXÞ. Then there is uhp 2 SGFEM such that

b ku  uhp k2H1 ðXÞ 6 C

N X

a¼1

inf ku  ua k2H1 ðxa Þ

ua 2va

b is a mesh-dependent constant, k  k Proof of the theorem can be found in the literature [49]. In the above theorem, C denotes an H1 -norm, and xa is a local patch on which

va is defined.

Essentially what the error bound amounts to is that the global error is controlled by the quality of the local errors. In other words, if high quality enrichment functions (and small local errors) can be utilized, then we can expect good global accuracy as well. Unfortunately, in most cases, high quality, closed-form, analytical enrichment functions are not available a priori. This is particularly true for problems involving nonlinearity, structural dynamics/transiency, multiple interacting crack surfaces, as well as many other general multi-scale problems. Due to the GFEM’s reliance upon high quality enrichment functions in order to deliver good global errors, and the limited availability of closed-form analytical enrichment functions, a more general version of enrichment function construction is required. A significantly more general extension of the GFEM, referred to as the Generalized finite element method with global-local enrichment functions (GFEMgl ) [24] has been proposed. 4. Generalized finite element method with global-local enrichments The GFEMgl [24,25] is a more general version of the standard GFEM in which enrichment functions are computed numerically during the course of a simulation, so as to customize the enrichment function to the specific problem under consideration. The overall goal is to alleviate the a priori requirement of high quality local approximation spaces. The GFEMgl requires the use of two separate discretizations: (1) fixed, coarse global discretization and (2) highly-adapted, potentially evolving local discretization defined locally where needed. Each discretization constitues a separate BVP to be analyzed, but there is a two-way coupling between the models, as described in subsequent sections. The following is a brief formulation for GFEMgl appropriate for LEFM applications. Formulations for GFEMgl approaches appropriate for material nonlineariy as well as propagation of cohesive fractures can be found in [55–57], respectively. Detailed formulations for time-dependent heat conduction can be found in [58].

280

P. O’Hara et al. / Engineering Fracture Mechanics 163 (2016) 274–302

4.1. Initial global problem The initial phase of the GFEMgl involves the solution of the 3D elasticity problem posed in Section 2.1 on a coarse FE mesh. The global model can be meshed using either tetrahedron or hexahedron elements. The crack itself is not modeled in the global BVP, similar to the approach taken in [25]. This essentially means that Cc ¼ ; in Fig. 1 for the initial global problem. The initial global BVP is formally posed as:

Find uG 2 vG ðXÞ  H1 ðXÞ such that 8v G 2 vG ðXÞ Z Z Z   r uG : e v G dX þ gp uG  v G dC ¼ Cu

X

X

f  v G dX þ

Z Cr


t  v G dC þ g p

Z Cu


 v G dC: u

ð21Þ

We may formally define the initial global FE solution space as

(

vG ðXÞ ¼ uG ¼

NG X ^a; u ^ a 2 R3 ; a ¼ 1; . . . ; NG Na ðxÞu

) ð22Þ

a¼1

in which N G is the number of nodes in the global mesh. The use of the gp parameter in (21) indicates that a penalty formulation is used to enforce Dirichlet boundary conditions. Other methods are certainly possible for the application of this boundary condition type, but the penalty method is used here due to its ease of implementation. 4.2. Local problem The global-local enrichment function is subsequently computed in a local BVP, which is driven by boundary conditions (BCs) from the previously solved initial global problem. Four-node tetrahedral elements are used in the local problem, regardless of the global element type used. If a global hexahedral mesh is used, coarse global elements are converted from hexahedral elements to local tetrahedral elements using the smallest vertex algorithm proposed in [59] in order to generate the local mesh. The local BVP is formally posed as:

 Find uL 2 vL ðXL Þ  H1 XL such that 8v L 2 vL ðXL Þ Z Z Z   r uL : e v L dX þ gp uL  v L dC ¼ f  v L dX T XL CL nðCL Cr Þ XL Z Z Z L L G L

þ v d C þ g v d C þ g u  t  T T T u  v dC: p p CL

Cr

CL

Cu

CL nðCL

ð23Þ

@ XÞ

An important aspect of the local problem is the use of the global solution, uG , as boundary conditions applied at the local boundary which are interior to the global domain. In this manner, information is transferred from the global BVP level down to the local BVP level, and thus constitutes the global-to-local coupling of the algorithm. Portions of the local boundary which intersect the global boundary are subjected to boundary conditions which come directly from the global BVP. We formally define the local GFEM solution space as

(

vL ðXL Þ ¼ uL ¼

) NL X
 ^ a ðxÞ þ u ~ a ðxÞ : NLa ðxÞ upa ðxÞ þ Hu

ð24Þ

a¼1

In the local BVP the crack surface is modeled using the approach presented in [17]. The discontinuity surface itself is modeled ^ a ðxÞÞ applied to nodes whose support is fully cut by the crack surface. Nodes through the use of a Heaviside enrichment ðHu whose supports contain the crack front are enriched with the asymptotic expansions for a sharp crack in a 2D elastic medium ~ a ðxÞÞ. With this enrichment strategy, the crack surface may be modeled completely independent of the underlying FE ðu mesh. Enrichment with continuous monomials ðupa ðxÞÞ can also be used in order to generate a higher polynomial order of approximation in the local problem. A locally graded mesh in the near crack tip region is still required along with the proper enrichment strategy in order to obtain high levels of accuracy. The reason for the necessity of local mesh refinement lies in the fact that the enrichment functions come from the 2D elasticity solution and not a 3D solution. While 3D enrichment functions may be possible in certain cases, it is more practical to use 2D expansions of the elastic solution as enrichment function, along with a properly refined mesh. It is the requirement of a highly adapted local mesh which necessitates the use of tetrahedral elements in the local problem. Tetrahedral elements allow for a straightforward local refinement scheme based on the marked edge algorithm [60,61], which yields a conforming, highly graded mesh. It should again be noted that the local solution is used only in the construction of the global approximation space, and not directly as part of the resulting displacement field.

P. O’Hara et al. / Engineering Fracture Mechanics 163 (2016) 274–302

281

4.3. Enriched global problem The enriched global BVP is then solved using the local solution, uL as a scale-bridging enrichment function. The use of uL as an enrichment function is the local-to-global aspect of the two-way coupling provided by the GFEMgl framework. The twoway coupling is illustrated in Fig. 2 (1) global BCs drive the local BVP and (2) local solution is used as enrichment in the enriched global BVP. It should be noted that multiple global-local iterations can be used so as to improve the quality of the globally-provided BCs in order to generate more accurate enrichment functions, and ultimately more accurate enriched global solutions. It should be stressed that these iterations are performed only to improve the quality of the enrichment function. There are no iterative requirements to satisfy equilibrium, as would be the case in the iterative global-local approaches found in the literature [62,63]. In the present work, equilibrium is satisfied in both the global and local domains at all times. The enriched global BVP is formally posed as:

Find uE 2 vE ðXÞ  H1 ðXÞ such that 8v E 2 vE ðXÞ Z Z Z Z   r uE : e v E dX þ gp uE  v E dC ¼ f  v E dX þ Cu

X

X

Cr


t  v E dC þ g p

Z Cu


 v E dC: u

ð25Þ

The enriched global GFEM solution space is defined as

(

vE ðXÞ ¼ uE ¼

NG X

a¼1

^a þ Na ðxÞu

X

) N b ðxÞugl b ðxÞ

ð26Þ

b2IE

in which IE is the set of global nodes enriched with uL ðxÞ; N a ðxÞ and N b ðxÞ are PoU functions provided by the global mesh, the n oT wðxÞ u v in which the terms uuL ðxÞ; uvL ðxÞ; uw term ugl b ðxÞ ¼ ub uL ðxÞ; v b uL ðxÞ; wb uL L ðxÞ are the Cartesian components of the local solution vector, uL ðxÞ and ub ; v b ; wb are global degrees-of-freedom. The GFEMgl is now applied to LEFM simulations – both static fracture as well as stable fatigue crack propagation - specifically for problems involving multiple, interacting crack surfaces. 5. Numerical examples 5.1. Crack shielding and amplification effects of micro-cracks: static crack analyses 5.1.1. Mode I crack configuration In this section we analyze the problem of crack shielding/amplification when there are one or two micro-cracks in the vicinity of a macro-crack front, in a Mode I configuration. The problem domain is subjected to tractions corresponding to a Mode I displacement field, as shown in Fig. 3. Point Dirichlet boundary conditions are applied only to prevent rigid body motions. Material properties are taken as Young’s Modulus, E ¼ 1:0 ksi and Poisson’s ratio, m ¼ 0:0. The value of m ¼ 0:0 is selected so as to allow for direct comparison against 2D numerical results available in the literature [29]. The GFEMgl

Fig. 2. Illustration of the GFEM methodology used in this paper.

282

P. O’Hara et al. / Engineering Fracture Mechanics 163 (2016) 274–302

Fig. 3. Model problem to evaluate effect of micro-cracks on macro-crack response, in Mode I configuration. Domain is subjected to tractions corresponding to Mode I displacement field, as indicated by red arrows (a). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)



 or tetrahedral GFEMgl elements. approach is investigated with global models meshed with either hexahedral GFEMgl HEX TET The domain itself is 2  2  0:5 in: and discretized with 9  9  1 hexahedron elements in the x-, y- and z-directions, respecgl tively. This is the global model used in the GFEMgl HEX approach. For hp-GFEM and GFEM TET approaches, the hexahedral elements are subsequently divided into 6, nested tetrahedral elements, thus generating models with equivalent global mesh gl densities for purpose of comparison. In the GFEMgl TET and GFEM HEX approaches, only scale-bridging, global-local enrichment functions are used to enrich the coarse-scale mesh in order to resolve the fine-scale features. In all numerical examples,  the global BVP uses a quadratic approximation space pG ¼ 2 , and the local BVP uses a cubic approximation space  loc p ¼ 3 . Discontinuous and singular enrichment functions are used to model the crack surfaces only in the local BVPs. It should be stressed that the discontinuity surfaces, both macro- and micro-cracks are then explicitly modeled on the coarse mesh with only the use of the scale-bridging enrichment functions, provided by the solution of the fine-scale BVP. Numerical reference solutions are generated with hp-GFEM [17], in which localized refinement along with discontinuous/singular enrichment functions are used directly on the initially coarse mesh to resolve the discontinuities in the displacement field.

5.1.1.1. Crack shielding. The first case considered is the crack shielding configuration, illustrated in Fig. 4. The relative dimen^ ¼ 2a^ . The computed Mode I SIFs for the macro-crack are shown in Table 1. In the sions, shown in Fig. 3(a), are taken as b 135 table, superscripts a; b indicate the absence and presence of micro-cracks, respectively. The reference solution for case a is somewhat trivial for the BCs applied here, and as noted in [29], K aI ¼ 1:0. The reference solution for case b is taken from [29], in which the same problem was solved using the multi-scale projection method. This type of problem has also been investigated analytically in [64]. As can be seen, the results compare very well with the reference values, and the crack shielding effect of the micro-cracks is captured accurately. It may be noted that in all numerical examples analyzed in this pffiffiffiffiffiffi paper, relevant SIFs are extracted using the cut-off function method (CFM) [41,42,65], and all SIFs are reported as ksi in:. gl Enriched global displacement contours are provided in Fig. 5 for both the GFEM gl HEX and GFEM TET approaches. As can be

seen from the contours, as well as the SIF values in Table 1, the GFEMgl can resolve fine-scale features on the fine-scale model, and directly upscale their effects into the coarse-scale model. The GFEMgl does so in an explicit manner with the

Fig. 4. Zoomed-in view of crack shielding configuration.

P. O’Hara et al. / Engineering Fracture Mechanics 163 (2016) 274–302

283

Table 1 Comparison of extracted SIF values for Mode I through-thickness macro-crack in a crack shielding configuration. In the table, superscripts a; b indicate the absence and presence of the micro-cracks, respectively.

Reference hp-GFEM GFEMgl TET GFEMgl HEX

K aI

K bI

1.0 0.9985 0.9983

 0:683 0.6816 0.6804

0.9965

0.6808

Fig. 5. Enriched global displacement contours in the neighborhood of the macro-crack tip.

incorporation of the actual micro-scale discontinuities into the coarse-scale solution field, and not just the implicit effects of the local micro-crack stress fields. The GFEMgl HEX approach is again used to analyze the same model problem, with the relative crack dimensions taken as

^ ¼ 2a^ . In this instance the micro-cracks are in closer proximity to the macro-crack, but they are also smaller cracks. The b 243 ref  0:6867, again Mode I SIF computed with the GFEM gl HEX is K I ¼ 0:6864, which compares well with the reference solution, K I computed using the multi-scale projection method [29]. It should be noted in this type of scenario, with crack surfaces in such close proximity, that in order to obtain accurate results the extraction radii used to compute K I must be kept sufficiently small so as to avoid the potential overlapping of an extraction domain with a neighboring discontinuity surface. ^ are provided in Fig. 6. As The coarse-scale displacement contours near the macro-crack front region for both values of b

can be seen from the figures, the global-local enrichment functions are able to accurately represent the multiple discontinuity surfaces in close proximity to each other in the coarse-scale mesh. In this instance all three discontinuity surfaces are explicitly represented within a single computational element. 5.1.1.2. Crack amplification. The second case considered is Mode I crack amplification, illustrated in Fig. 7. The domain itself as well as the macro-crack have the same dimensions as those used in the previous example. The distances between the macro
¼ 0:10, both as illustrated in Fig. 7, yielding a ratio of a
¼ 0:2. This
¼ 0:02 and b crack and micro-crack fronts are taken as a
b numerical example has also been investigated numerically with the XFEM by Loehnert and Belytschko [66]. Analytical solutions have also been derived, which can be found in [64,67]. The analytical solutions provide K I values at all three crack

284

P. O’Hara et al. / Engineering Fracture Mechanics 163 (2016) 274–302

Fig. 6. GFEM gl HEX enriched global displacement contours in the neighborhood of the macro-crack tip for different relative crack sizes and spacing.

Fig. 7. Zoomed-in view of crack configuration for aligned micro-crack.


For completeness, the analytical values for K I at each of these points is included below, as presented in
, and b. fronts, x ¼ 0; a [67] ðx¼0Þ KI

 sffiffiffi
E 1  ab

b   K Iðx¼1Þ  ¼
K 1  a
a

ð27Þ


b


Þ ðx¼a KI

 0 1,0sffiffiffiffiffiffiffiffiffiffiffiffi1
E 1  ab


b b ðx¼1Þ   1A @  1A  K I ¼ @


a a a K 1

ð28Þ


b

 1, rffiffiffiffiffiffiffiffiffiffiffiffi! E 1  ba


a ðx¼1Þ A 1 
 KI ¼ @1 
a b K 1  b
0


ðx¼bÞ KI

ð29Þ

in which KðÞ and EðÞ are complete elliptic integrals of the first and second kind [68,69], respectively, and the far-field stress ðx¼1Þ

¼ K ref ¼ 1:0 in this case. intensity, K I I The extracted K I values at both the macro-crack and micro-crack fronts are provided in Table 2. As can be seen, the extracted K I values generated with the GFEMgl compare very accurately with the analytical values, falling within 2% of ¼ 1:0 in the absence of a micro-crack, so the K I values the reference values. It may also be pointed out that the value of K ref I are indeed amplified in this scenario. Enriched global displacement contours in the macro-crack front region are presented in gl Fig. 8 for both the GFEMgl TET and GFEM HEX approaches. It can once again be clearly seen in the plots that the micro-crack is accurately represented in an explicit manner in the coarse-scale model.

5.1.2. Mixed-mode crack configuration In this section, the effect of micro-crack shielding/amplification on a mixed-mode macro-crack is analyzed. As shown in  Fig. 9(a), a domain with unit thickness ðt ¼ 1:0Þ is subjected to unit normal tractions ty ¼ 1:0 in the y-direction. Material properties are taken as Young’s Modulus, E ¼ 1:0 ksi and Poisson’s ratio, m ¼ 0:0. Poisson’s ratio is taken as zero so as to allow for an initial comparison against a reference solution [41]. There is a macro-crack in a mixed-mode configuration, at an angle of 450 counterclockwise with respect to the global y-axis, and a length of amacro ¼ w2 ¼ 0:5.

P. O’Hara et al. / Engineering Fracture Mechanics 163 (2016) 274–302

285

Table 2 Comparison of extracted SIF values for Mode I through-thickness cracks in a crack amplification configuration. ðx¼0Þ

Reference hp-GFEM GFEMgl TET GFEMgl HEX


Þ ðx¼a


ðx¼bÞ

KI

KI

KI

1.1675 1.1874 1.1865

0.8053 0.8162 0.8117

0.5343 0.5457 0.5358

1.1914

0.8188

0.5453

(a) Global hexahedral elements.

(b) Global tetrahedral elements.

Fig. 8. Enriched global displacement contours in the neighborhood of the macro-crack tip. Again it can be seen the GFEMgl can explicitly model the effects of the micro-crack on the macro-crack in the coarse-scale domain.

5.1.2.1. Crack amplification. The first case analyzed is that of mixed-mode crack amplification, illustrated in Fig. 9(b). The dis
¼ 0:10, as shown in the bot
¼ 0:02 and b tances between the macro-crack tip and the micro-crack tips are again taken as a tom portion of Fig. 9(b), with the micro-crack parallel to and co-planar with the macro-crack. In this instance, both Mode I and Mode II SIFs are relevant and of the same basic magnitude. For the purpose of comparison with results from the literature, the SIFs are normalized as aðbÞ

aðbÞ

kIðIIÞ ¼

K IðIIÞ pffiffiffiffiffiffiffiffiffiffi t y 2p w

ð30Þ

in which IðIIÞ indicate the fracture mode under consideration, and aðbÞ indicate the absence and presence of micro-cracks, respectively. Due to the mixed-mode nature of the problem, the energy release rate is also computed at the crack front, to provide a scalar value that can be easily compared to illustrate the amplification/shielding effect of the micro-cracks. With the choice of E; m and the pseudo-2D nature of this problem, the energy release rate is computed based on the normalized SIF 2

2

values simply as, G ¼ kI þ kII . The results in Table 3 (as well as Table 4) were generated using 2 global-local iterations to gl improve the global boundary conditions and thus improve the extracted SIF values for the GFEMgl TET and GFEM HEX approaches. There is no numerical reference value for case b, as only the problem with a mixed-mode through-thickness macro-crack is analyzed in [41]. As such, case a is used for the initial comparison to establish confidence in the extraction of mixed-mode SIFs. As can be seen from the table, the results for case a compare very well with the reference values, which was computed with high-order FE approximation of p ¼ 8. The amplification effect of the micro-crack is also clearly seen in the overall

increase in Gb as compared to Ga , of approximately 47%. It should be noted that both GFEMgl approaches provide results which compare very favorably against the results generated with the hp-GFEM. Enriched global displacement contours are provided in Fig. 10, in which the effect of the micro-cracks are again accurately represented in the global model. 5.1.2.2. Crack shielding. The same micro-crack surface as was utilized in the mixed-mode amplification simulation is utilized again in this scenario, but there are 2 micro-crack surfaces, translated into a crack shielding configuration. In this case the ^ as illustrated in the top of Fig. 9(b), is taken as b ^ ¼ 0:08. The extracted SIF values are presented in relative dimension, b, Table 4, in which the shielding effect of the micro-cracks is clearly illustrated. The energy release rate decreases by approximately 29% due to the presence of the 2 micro-cracks in the vicinity of the macro-crack front. The SIF values for the case with no micro-cracks are included in Table 4 for ease of comparison because in this instance the values are not simply, 1:0, as

286

P. O’Hara et al. / Engineering Fracture Mechanics 163 (2016) 274–302

Fig. 9. Numerical model for micro-crack effects on macro-crack response for a mixed-mode crack configuration. Domain is subjected to unit normal  tractions t y ¼ 1:0 , as indicated by red arrows in (a). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Table 3 Comparison of extracted SIF values for mixed-mode through-thickness macro-crack in a crack amplification configuration. In the table superscripts a; b indicate the absence and presence of the micro-crack, respectively. a

Reference hp-GFEM GFEMgl TET GFEMgl HEX

kI

kII

a

Ga

kI

kII

Gb

0.6013 0.6008 0.6011

0.2910 0.2909 0.2901

0.4462 0.4455 0.4455

– 0.7258 0.7199

– 0.3549 0.3537

– 0.6527 0.6434

0.6009

0.2911

0.4458

0.7270

0.3560

0.6553

b

b

Table 4 Comparison of extracted SIF values for mixed-mode through-thickness macro-crack in a crack shielding configuration. In the table superscripts a; b indicate the absence and presence of the micro-cracks, respectively. a

Reference hp-GFEM GFEMgl TET GFEMgl HEX

kI

kII

a

Ga

kI

kII

Gb

0.6013 0.6008 0.6011

0.2910 0.2909 0.2901

0.4462 0.4455 0.4455

– 0.4591 0.4593

– 0.3240 0.3242

– 0.3158 0.3160

0.6009

0.2911

0.4458

0.4592

0.3242

0.3159

b

b

was the case in the Mode I configuration. It is again noted that results produced with both GFEMgl approaches compare favorably against those generated with hp-GFEM. Enriched global displacement contours are presented in Fig. 11, where the gl micro-cracks are accurately modeled in the global model using both the GFEMgl HEX and GFEM TET approaches.

P. O’Hara et al. / Engineering Fracture Mechanics 163 (2016) 274–302

(a) Global hexahedral elements.

287

(b) Global tetrahedral elements.

gl Fig. 10. Enriched global displacement contours in the neighborhood of the macro-crack tip for both GFEMgl HEX (a) and GFEM TET (b) approaches.

(a) Global hexahedral elements.

(b) Global tetrahedral elements.

gl Fig. 11. Enriched global displacement contours in the neighborhood of the macro-crack tip for both GFEMgl HEX (a) and GFEM TET (b) approaches.

5.2. Crack shielding effect on fatigue crack growth in aerospace structure In the previous numerical examples, with small geometric size (on the order of a material coupon), it could be possible to attempt to model sub-continuum level material features in order to give the proposed method a more ’multi-scale feel’. But as presented, the target application area for the proposed method is that of structural scale models. In such scenarios, the relevant structural scales are potentially on the order of feet, whereas the life of the structure is dependent upon the accurate resolution of geometric features (in this case cracks) which are in general, on the order of <1 in. The quantities of greatest interest are in fact the strong strain/stress gradients developing in the regions of localized interest (in this case crack fronts) which themselves develop over regions on the order of  1 in:. In such scenarios it is not a requirement, nor even practical to model micro/grain-level material features, and the use of a purely homogeneous material model throughout the entire domain still yields truly multi-scale application problem. In this section we simulate fatigue cracks propagating in a representative aerospace structure. The structure is a three-bay stiffened panel (Fig. 12) subjected to cyclic normal tractions with magnitude 12 ksi applied in the z-direction. Material

y z

x

δ

Fig. 12. Crack shielding in representative aerospace structure.

288

P. O’Hara et al. / Engineering Fracture Mechanics 163 (2016) 274–302

properties and Paris law parameters for a Ti6242s alloy were obtained from [70], and taken as in: E ¼ 17; 100 ksi;m ¼ 0:0; C ¼ 1:17e15 cycle and m ¼ 5:47. The value of m ¼ 0:0 was selected because it serves to simplify the structural response of the panel itself, and yields the Mode I type behavior which would be expected based on the crack configuration and loading conditions. The load ratio was taken as R ¼ 0:1. All load ratios, R, in this work are taken to be nonnegative to avoid the complications which arise due to crack face contact which becomes relevant when the load ratio is negative. All dimensions and locations reported in this section are in inches and the origin is taken as the geometric center of the top surface of the panel. The panel skin has dimensions ð19:6  0:065  34Þ in the x-, y- and z-directions, respectively. There are two flanges located at x ¼ 9:8 with dimensions ð0:065  2:5  34Þ in the x-, y- and z-directions, respectively. There are four stiffeners located at z ¼ 5; 10 with dimensions ð17:6  1:25  0:065Þ in the x-, y- and z-directions, respectively. The stiffeners taper linearly from a depth of 1:25 to 0:75 at each end. There are 2 through-thickness cracks in the flange located at x ¼ 9:8. The first crack, Crack 1 has initial length a1 ¼ 0:5 with the crack mouth located at ð9:8  2:5  2:67Þ. The second crack, Crack 2 has initial length a2 ¼ 0:2 with the crack mouth located at ð9:8  2:5  3:12Þ. The cracks are separated by d ¼ 0:45. 5.2.1. Static fracture analysis The global models for the GFEMgl analyses are provided by standard linear finite element models, 4-node tetrahedron elegl gl ments for the GFEM gl TET , and 8-node brick elements for the GFEM HEX . In both GFEM approaches, the global model is enriched G with monomials to a resulting p-order of p ¼ 2. Crack surfaces are not modeled explicitly in the global model, as they are inserted only via the scale-bridging global-local enrichment functions. Two global local iterations are preformed for the static fracture analyses. In the subsequent fatigue propagation analyses, two global-local iterations are only performed for the initial cracked configuration, but only one global-local iteration is used for the subsequent propagation steps. Simulations are first run for each crack separately, so as to provide baseline data for comparison with the subsequent

crack interaction analyses. For GFEM gl simulations containing both cracks, i.e. crack interaction analyses, there are two different options for selection of the local domain, as illustrated in Fig. 13. The first option is to select one local domain such that it contains both crack surfaces, as in Fig. 13(a). It is also possible to create two separate local domains, each containing one crack surface, as shown in Fig. 13(b). The second option serves to illustrate quite clearly the impact of the two-way coupling of the GFEMgl , as opposed to the one-way coupling used in a more traditional global-local FEM or sub-modeling strategy. Results for both option are presented. In Fig. 13 global ‘seed nodes’, as well as the local domain boundaries are indicated by circular glyphs and solid lines, respectively. It should be noted that the crack surfaces themselves are not modeled in the hexahedral mesh, and the crack surfaces are only shown in the hexahedral model to provide a location of reference between the local and global domains. Results in the form of the crack mouth opening displacement (CMOD) as well as K I are presented in Table 5. Results presented corresponding to an hp-GFEM approach are taken as the reference values for the purpose of computing relative differences. In the table, the superscript a indicates the presence of only one crack surface in the analysis, b indicates the

(a) One local domain, containing both cracks.

(b) Two, separate local domains.

Fig. 13. Two different options for the GFEMgl analysis of two interacting flange cracks.

289

P. O’Hara et al. / Engineering Fracture Mechanics 163 (2016) 274–302 Table 5 Comparison of extracted K I and CMOD values for Crack 1 and Crack 2. Crack 1

Rel. diff.

Crack 2

Rel. diff.

19.6339 19.4371

– –

10.8304 3.1434

– –

2.4303e3 2.3573e3

– –

8.4095e4 1.3928e4

– –

bð1Þ

19.3385 19.4992

0.0150 0.0032

10.7829 3.0351

0.0044 0.0344

bð2Þ

19.1768

0.0134

3.2377

0.0300

2.4218e3 2.3624e3 2.3390e3

0.0035 0.0022 0.0078

8.3530e4 1.2796e4 1.5845e4

0.0067 0.0813 0.1376

bð1Þ

20.4976 20.2343

0.0440 0.0410

11.3017 3.2037

0.0435 0.0192

bð2Þ

20.2441

0.0415

3.3667

0.0711

2.4286e3 2.3622e3 2.3474e3

0.0007 0.0021 0.0042

8.3837e4 1.3244e4 1.5345e4

0.0031 0.0491 0.1017

hp-GFEM K aI K bI CMODa CMODb GFEMgl TET K aI KI

KI CMODa CMODb(1) CMODb(2) GFEMgl HEX K aI KI

KI CMODa CMODb(1) CMODb(2)

presence of both crack surfaces, i.e. a crack interaction analysis. For the two GFEMgl approaches, the superscript ð#Þ, indicates bð1Þ

the number of local domains utilized. For instance, K I

refers to the scenario when 1 local domain containing both crack

bð2Þ KI

indicates the use of two, separate local domains, each containing one crack surface. surfaces is used, and The shielding effect of Crack 1 on Crack 2 is clearly illustrated in the table. The K I value for Crack 2 drops from approximately 10 to about 3, due to the presence of Crack 1. The CMOD value drops by over a factor of 4 in the presence of the larger crack. The values computed for Crack 1 are not significantly impacted due to the presence of the smaller crack surface. Good agreement is seen between the GFEMgl approaches and the reference values, with relative difference values below 5% in most cases. There are some larger errors seen for some CMOD values for Crack 2 in the interaction analyses. In these instances the reference values themselves are very small, as such even small differences in the computed values can yield relatively large errors. As mentioned previously, the use of two separate local domains to model the interacting cracks serves to clearly illustrate the benefit of the two-way coupling between the global and local domains in allowing for the crack interaction to be modeled accurately in the global model, even though each local problem has only one crack surface. The CMOD and K I values computed from the local and enriched global solutions corresponding to the first global-local iteration are provided in gl Table 6. Results are provided for both the GFEMgl TET and GFEM HEX . The values computed from the local solution show very little difference between Crack 1 and Crack 2, which is indicative of not only the poor boundary conditions provided by the global solution (no crack surfaces were modeled in the global domain), but also show that there is no detection of interaction between the two crack surfaces. This is the expected behavior at this point, and the best that one could expect to do with a traditional sub-modeling strategy in which the interaction between the crack surfaces is not detectable in the submodels. When looking at the enriched global solution, it is clear that even at this juncture in the analysis, the interaction between the two crack surfaces is accounted for, as the K I value for Crack 1 is now significantly larger than that for Crack 2, indicating the shielding effect. Results for the second global-local iteration are provided in Table 7. In the second iteration, the local values are now indicative of significantly more accurate globally-provided boundary conditions, which now also account for the interaction of the two crack surfaces, even though each local domain contains only one crack surface. This behavior is again attributed to the two-way coupling between the global and local domains. The enriched global values are now in good agreement with the

Table 6 Comparison of K I and CMOD values from local and enriched global domains, corresponding to first global-local iteration. Crack 1

Local (global tets) Enriched global (tets) Local (global hex) Enriched global (hex)

Crack 2

KI

CMOD

KI

CMOD

5.2064 18.8740 5.1178 20.3097

2.6986e4 2.2436e3 2.8562e4 2.3176e3

5.1744 3.5579 5.2083 3.6260

2.9098e4 1.9578e4 2.9188e4 1.7985e4

290

P. O’Hara et al. / Engineering Fracture Mechanics 163 (2016) 274–302

Table 7 Comparison of K I and CMOD values from local and enriched global domains, corresponding to second global-local iteration. Crack 1

Local (global tets) Enriched global (tets) Local (global hex) Enriched global (hex)

Crack 2

KI

CMOD

KI

CMOD

19.2207 19.1768 20.1337 20.2441

2.2625e3 2.3390e3 2.2976e3 2.3474e3

3.6843 3.2377 3.6630 3.3667

1.9339e4 1.5845e4 1.7770e4 1.5345e4

reference values due to the improved quality of the local enrichment functions, which themselves are a result of the improved local boundary conditions provided for the first enriched global solution. 5.2.2. Fatigue crack propagation 5.2.2.1. Case a: each crack analyzed separately. As in the static fracture case, baseline data is first generated through analyses run with each crack separately. In each analysis, the crack of interest is propagated to twice its initial length in 10 uniform crack steps, meaning Damax is set to 0:1Dainitial . Due to the thin through thickness dimension of the stiffener, K I is extracted only at the center vertex of the crack front, and the entire crack front is advanced uniformly. Load cycles corresponding to a particular crack advancement is then calculated incrementally according to (17). The crack advance as a function of loading cycles is plotted for Crack 1 in Fig. 14(a), and in Fig. 14(b) for Crack 2. Reference curves are also provided, and the reference data is again generated numerically through the use of hp-GFEM. The curves show fairly large differences for the three approaches analyzed. The large error values are in some sense a result of the cycle data being a derived quantity, which depends nonlinearly on a quantity computed from the solution itself, K I . The nonlinear relationship is clearly evident from (17), in which the cycle data is related to the K I value, raised to the m power. As such, errors in the K I values are amplified quite significantly for materials with large values of m, such as the Tialloy considered here. The values of K I are plotted as a function crack step, and thus increasing crack length, for both cracks in Fig. 15. As can be seen from the plot, there is relatively good agreement seen in the K I values computed using the three approaches considered. The value of m ¼ 5:47 for the selected Ti alloy is very high for a structural metal. The fatigue life can be computed again using the same extracted SIF ranges, but with a lower value of m ¼ 2:0. With an m value of m ¼ 2:0, the dependency of the cycle data on the computed K I values is not as strong. Crack advance versus load cycle data for both cracks is provided in Fig. 16 corresponding to m ¼ 2:0. From this plot, it is seen that there is much better agreement between the curves, simply due to the change in the material parameter, and not due to any sort of improvement in the computed GFEMgl solutions. To quantify the error levels more precisely, a discrete L2 -norm of the error in a quantity ðÞ, between any two data sets is computed as

vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi o2 uP n u ðÞref  ðÞGFEM u i i Lerror ðÞ ¼ u o2 t 2 Pn ð  Þ i ref i

gl Fig. 14. Crack length versus loading cycles computed with both GFEM gl TET and GFEM HEX approaches.

ð31Þ

P. O’Hara et al. / Engineering Fracture Mechanics 163 (2016) 274–302

291

Fig. 15. Extracted K I value as a function of increasing crack length for Crack 1 and Crack 2.

Fig. 16. Crack length versus loading cycles for Crack 1 and Crack 2 corresponding to m ¼ 2:0.

in which i is the number of crack steps considered. The Lerror values are provided in Table 8 for both GFEMgl approaches. As 2 can be clearly seen, the error values for cycle data corresponding to m ¼ 2:0 is significantly lower that those for m ¼ 5:47, as would be expected. The values for m ¼ 2:0 are still larger than those computed for K I , because again the dependency, while not as strong as in the case for m ¼ 5:47, is still greater than linear. This speaks to the need for accurate prediction of K I in a fatigue crack propagation simulation. 5.2.2.2. Case b: both cracks present. The interaction of Crack 1 and Crack 2 in a fatigue propagation simulation manifests itself in the form of crack arrest for the smaller crack, in this case Crack 2. The basis for this behavior is in the nonlinear scaling, (18), of vertex advancement so as to ensure that each crack vertex in the simulation is subjected to the same number of

Table 8 Computed Lerror values corresponding to K I as well as crack length 2 versus load cycles data for Crack 1 and Crack 2. Crack 1

Crack 2

0.0185 0.0727 0.0248

0.0245 0.0719 0.0297

0.0339 0.1708 0.0663

0.0322 0.1944 0.0725

GFEMgl TET K aI Cycle data (m ¼ 5:47) Cycle data (m ¼ 2:0) GFEMgl HEX K aI Cycle data (m ¼ 5:47) Cycle data (m ¼ 2:0)

292

P. O’Hara et al. / Engineering Fracture Mechanics 163 (2016) 274–302

loading cycles, and thus the cracked configuration is physically valid at any point in the simulation. Plugging in the values of K I from the static analyses into Eq. (18) with m ¼ 5:47, the crack advance of Crack 2 for a given number of cycles is  m  3:1434 5:47 K Crack2 I DaCrack 2 ¼ K Crack DaCrack 1 ¼ 19:4371 DaCrack 1 ¼ ð4:6985e  5ÞDaCrack 1 . Even for a value of m ¼ 2:0; DaCrack 2 ¼ 1 I

ð0:0261ÞDaCrack 1 , which is still essentially indiscernible, even for the first crack step. As illustrated in Fig. 15, K I increases for increasing crack length (and only Crack 1 is significantly increasing in length), so the discrepancy in the K I for Crack 1 and Crack 2 becomes larger with each crack step, and thus the advancement of Crack 2 becomes completely negligible, and Crack 2 ultimately arrests. This behavior is illustrated in Fig. 17(a) for m ¼ 5:47 and in Fig. 17(b) for m ¼ 2:0. Data is provided corresponding for both gl the GFEMgl TET and GFEM HEX approaches. In both instances it can be seen that Crack 1 propagates from 0:5 in: to 1:0 in:, whereas the length of Crack 2 remains unchanged throughout the course of the simulation. Reference curves are also provided for Crack 1 corresponding to Case a, i.e. Crack 2 not present in the analysis. As can be seen, there is only a small difference between the crack advancement versus load cycles for Case a and Case b due to the small difference in the K I values in both instances. As would be expected, this manifests itself as a larger discrepancy between the curves for m ¼ 5:47 than for m ¼ 2:0. This is expected due to the stronger dependence of the loading cycles on the SIF range seen at the crack front for a larger value of m. Zoomed in views of the crack advance versus load cycles for Crack 2 are provided in Fig. 18(a) for m ¼ 5:47 and in Fig. 18 (b) for m ¼ 2:0, over the load cycle range from Case b. The Reference curves in the figures show the crack advancement versus load cycles data from Case a. The data in these plots clearly illustrates the arresting of Crack 2 in both instances. The final data of interest are the extracted SIFs at each crack front throughout the course of the simulation. The values of K I as a function of crack step are provided in Fig. 19. As can be seen, K I is monotonically increasing with increased crack length for Crack 1, as would be expected. The K I values for Crack 2 however, show a decrease throughout the course of

Fig. 17. Crack length versus loading cycles computed for crack shielding scenario.

Fig. 18. Illustrates the crack arrest of Crack 2 due the presence of Crack 1.

P. O’Hara et al. / Engineering Fracture Mechanics 163 (2016) 274–302

Fig. 19. Illustrates the crack shielding on the values of K I computed for Crack 2.

Fig. 20. Illustration of crack interactions for various, small values of d.

293

294

P. O’Hara et al. / Engineering Fracture Mechanics 163 (2016) 274–302

the simulation (even though the length of Crack 2 is technically increasing, albeit by negligible amounts) until K I ¼ 0:0. This is consistent with what is expected from the exhibited crack arrest behavior. 5.3. Coalescence of multiple crack surfaces In this section we consider a domain with two initially parallel, but vertically offset through-thickness cracks propagating towards each other in a rectangular domain with dimensions 2  6  1 in:. The domain is subjected to unit, normal tractions applied to the top face, with a built-in support on the bottom face, as shown in Fig. 20(a). The material properties are taken as Young’s Modulus, E ¼ 2e5 ksi, and Poisson ratio, m ¼ 0:0, to yield essentially two-dimensional behavior. For each simulation, the initial crack surfaces have dimensions 0:6  1:0 in:, and are separated by a vertical distance, d ðin:Þ, again as shown in Fig. 20(a). The problem of multiple cracks propagating towards, and interacting with one another is a very complicated problem. There is experimental evidence which indicates that two crack fronts approaching each other in an initially Mode I manner will tend to deflect away from each other, and then turn back towards each other as shown in for instance [71,72]. This issue of crack path stability in such scenarios has also been investigated analytically in [71,73]. The anticipated interaction behavior is captured and illustrated in Fig. 20 for various values of d. As expected, the crack fronts avoid one another upon approach, and then propagate back towards each other some distance behind the crack fronts. As would be expected, the interaction behavior dies out as d becomes larger, and the two crack paths are no longer significantly effected by the presence of the other crack. Fig. 21 shows that when d becomes large enough, the crack trajectories show very little influence on each other. This is particularly evident in the case with d ¼ 2:0, in which the cracks show virtually no influence on each others trajectories. Of particular interest is then the case for which d ¼ 0:0, and the two crack surfaces are initially co-planar. In this instance, the crack fronts do not change trajectory upon approach, as shown in Fig. 22. The crack surfaces coalesce and form one large crack surface, which is in accordance with experimental data [74,75], as well as analytical derivation [71] in the absence of an initial disturbance from a Mode I path. In this simple example, there is little interest beyond the point of coalescence, since the domain is completely severed. In practice, however; crack coalescence may play a very important role in the life prediction of engineering components. Upon coalescence a larger, potentially dominant crack surface is formed which can significantly impact the expected service life of a component as discussed in [76]. With this in mind, the subsequent numerical examples focus on crack coalescence.

Fig. 21. Illustration of crack interactions for large values of d.

(a) 3-DView

(b) SideView

Fig. 22. Final propagation step for d ¼ 0:0.

P. O’Hara et al. / Engineering Fracture Mechanics 163 (2016) 274–302

(b) SideView

(a) 3-DView Fig. 23. Initial configuration for two, co-planar circular crack surfaces.

(a) hp-GFEM

(c) GFEM gl TET

(e) GFEM gl HEX

(b) hp-GFEM

(d) GFEM gl TET

(f) GFEM gl HEX

gl Fig. 24. Diagonal view of global models for crack step 9 and crack step 37 for hp-GFEM, GFEMgl TET and GFEM HEX approaches.

295

296

P. O’Hara et al. / Engineering Fracture Mechanics 163 (2016) 274–302

gl Fig. 25. Top view of global models for crack step 9 and crack step 37 for hp-GFEM, GFEMgl TET and GFEM HEX approaches.

P. O’Hara et al. / Engineering Fracture Mechanics 163 (2016) 274–302

297

5.3.1. Circular cracks In this section we analyze a problem with two initially parallel, co-planar circular cracks in a domain with dimensions 6  2  6. The material properties are taken as Young’s Modulus, E ¼ 2e5 ksi, and Poisson ratio, m ¼ 0:3. The domain is subjected to cyclic unit, normal tractions (red arrows) applied to the top and bottom surfaces with load ratio R ¼ 0:0, and point boundary conditions (green arrows) are applied only to prevent rigid body motions, as shown in Fig. 23. The values of C; m are not provided as no cycle data is presented in this section. The crack surfaces propagate towards one another, and coalesce. The stress intensities are then greatest near the point of coalescence, and these vertices propagate faster for a given number of load cycles. As such, the coalesced surface tends to grow into a large circular crack surface, which is consistent with other results from the relevant literature [22,23]. Crack propagation results for crack steps 9 and 37 are presented for three different approaches: hp-GFEM, GFEMgl TET and GFEMgl HEX in Fig. 24. As can be seen, each of the three approaches yield similar results in terms of the final propagated crack surface. The significant difference in the mesh densities required for the hp-GFEM approach versus the GFEMgl approaches is also evident in the figures. A top view of the same propagation results is presented in Fig. 25, in which the fixed, coarse mesh utilized in the GFEMgl approach is more clearly illustrated, due to the vantage point. As can be seen, the global discretizations remain unchanged, and thus regions of high mesh density do not adaptively follow the propagated crack front, which is in stark contrast to the results presented for the hp-GFEM approach (Fig. 25(a) and (b)). The crack surface re-meshing algorithm [22] for the GFEM gl HEX case utilizes a different characteristic element size away from the crack front than that used in the GFEMgl TET case. The impact of varying the parameter in the two cases is clearly seen in Figs. 24 and 25 in the the significantly coarser crack surface mesh created behind the crack front in the GFEMgl HEX results. 5.3.2. Inclined penny cracks In this section we analyze two inclined penny-shaped cracks embedded in a domain with dimensions 4:0  2:0  4:0 in: in the x-, y- and z-directions respectively, as illustrated in Fig. 26. The domain is subjected to cyclic, unit normal tractions with load ratio R ¼ 0:0 applied to the top and bottom faces of the domain in the y-direction, and point Dirichlet boundary

Fig. 26. Solid model used for GFEM gl TET approach to simulate coalescence of inclined penny-shaped cracks.

Fig. 27. Side view of initial cracked configuration.

298

P. O’Hara et al. / Engineering Fracture Mechanics 163 (2016) 274–302

conditions are applied only to prevent rigid body motions. The material properties are taken as Young’s Modulus, E ¼ 1:0e3 ksi and Poisson’s ratio, m ¼ 0:3. The values of C; m are again not provided as no cycle data is presented in this section. The initial crack surfaces are penny-shaped with radius, r ¼ 0:1 in:, and they are inclined at an angle of a ¼ 45 with respect to the global z-axis. The cracks are separated by a distance of dz ¼ 0:2828 in:, and they are offset in the y-direction by a distance of, dy ¼ 0:0808 in: Fig. 27 illustrates the initial crack configuration in detail. In this problem, which is analyzed in great detail within the hp-GFEM context [22], the cracks grow towards one another in a mixed-mode manner. At the point of coalescence, the two crack surfaces are not co-planar. Upon coalescence there is again a highly non-convex region of the crack front, which is successfully propagated out into a circular crack front, at an elevation near the centroid of the initially cracked configuration. Fig. 28 shows the final propagated crack configuration. gl Fig. 29 shows GFEMgl TET and GFEM HEX results for crack steps 0; 14; 28 and 65 from a top view. In both cases, it can be seen that the crack is able to propagate out from the initial highly non-convex peanut shape seen in Fig. 29(a) and (b) into the ultimately circular shape seen in Fig. 29(e) and (f). Off-diagonal views of the same crack propagation steps are provided gl for the GFEMgl TET and GFEM HEX models from slightly different vantage points so as to illustrate the mixed-modality of the initial cracked configurations as well as the resulting propagated crack surfaces. In addition to the vantage points being from gl different sides of the global model, there is also a different angle of inclination used for the GFEMgl TET and GFEM HEX results so as to minimize the clutter which arises from the internal element edges present in the global meshes. It may be noted that the propagated crack surfaces in this numerical example compare very favorably with the results generated using hp-GFEM in [22] (see Fig. 30).

(a) TopView

(b) SideView Fig. 28. Final propagated crack surface, generated with GFEM gl TET approach.

P. O’Hara et al. / Engineering Fracture Mechanics 163 (2016) 274–302

299

gl Fig. 29. Top view of global models for inclined penny-shaped crack coalescence results for GFEMgl TET and GFEM HEX approaches. Results provided corresponding to crack steps 14; 28 and 65.

300

P. O’Hara et al. / Engineering Fracture Mechanics 163 (2016) 274–302

gl Fig. 30. Diagonal view of global models for inclined penny-shaped crack coalescence results for GFEM gl TET and GFEM HEX approaches. Results provided gl corresponding to crack steps 0; 14; 28 and 65. Different angles of inclination are used in the GFEMgl TET and GFEM HEX results in order to minimize clutter in the images due to the internal element edges visible in the global models.

6. Conclusions This paper presents an extension of the GFEMgl appropriate for fracture mechanics analyses involving multiple, interacting crack surfaces on fixed, coarse FE meshes. In the approach taken, the cracks themselves are not modeled directly in the coarse scale GFEM model. The crack surfaces are modeled in the global mesh exclusively through the use of numericallygenerated scale-bridging enrichment functions. The enrichment functions are provided in the form of the solution of local

P. O’Hara et al. / Engineering Fracture Mechanics 163 (2016) 274–302

301

BVPs defined only in the regions of existing mechanically-short cracks. The hp-GFEM as presented in [17] is used to model the crack only in the fine-scale BVPs. The proposed approach was successfully verified against hp-GFEM as well as analytical reference solutions for crack amplification and shielding scenarios with micro-cracks in the vicinity of a macro-crack tip. Both Mode I as well as mixed-mode crack configurations were analyzed. In the GFEMgl HEX approach, multiple discrete discontinuity surfaces were accurately represented entirely within one computational element. This type of application, i.e., multiple discrete discontinuities in close proximity to one another, is a prime example of a relevant engineering application in which the GFEMgl methodology is able to circumvent the lack of a priori available analytical enrichment functions. The approach was also successfully applied to a crack shielding scenario involving two propagating crack surfaces in the flange of a representative aerospace panel. The GFEMgl approach produced SIFs as a function of crack advance length which compared very favorably with the reference results generated with hp-GFEM. In this instance each crack surface was modeled in its own local domain, but the interaction between the two surfaces - arresting of the smaller crack - was still captured accurately thanks to the twoway coupling of the GFEMgl methodology. The final application area investigated was that of crack coalescence. A co-planar Mode I as well as a fully mixed-mode (and thus not co-planar) numerical example were presented. In both instances the coalescence, as well as the propagation out of the highly non-convex regions of the crack front was able to be accurately modeled on the fixed, coarse FE meshes. While beyond the scope of the current work, improved methods of computing the coalescence tolerance from the SIFs at the crack front is an area of focus for future work - for use within both the hp-GFEM as well as GFEMgl approaches. Acknowledgments The financial support of the Air Force Office of Scientific Research through Task No. 09RB07COR (Dr. David Stargel, Program Officer), is gratefully acknowledged. References [1] NASGRO/FLAGRO. Fatigue crack growth computer program NASGRO/FLAGRO version 2.0. Structures and mechanics division, NASA/Lyndon B. Johnson Space Center, Houston (TX); 1992. [2] FEACrack version 2.7. Quest reliability, LLC. Boulder, Colorado. . [3] Wawrzynek P, Martha L, Ingraffea A. A computational environment for the simulation of fracture processes in three-dimensions. In: Rosakis AJ, editor. Analytical, numerical and experimental aspects of three dimensional fracture process, vol. 91. New York: ASME AMD, ASME; 1988. p. 321–7. [4] Zentech International Limited. . ZENCRACK; 2008. [5] Schöllmann M, Fulland M, Richard H. Development of a new software for adaptive crack growth simulations in 3-D structures. Eng Fract Mech 2003;70:249–68. [6] Ural A, Wawrzynek P, Ingraffea A. Simulating fatigue crack growth in spiral bevel pinion Tech rep ARL-CR-0531. NASA/Army Research Laboratory; August 2003. [7] Ural A, Heber G, Wawrzynek P, Ingraffea A, Lewicki D, Neto J. Three-dimensional, parallel, finite element simulation of fatigue crack growth in a spiral bevel pinion gear. Eng Fract Mech 2005;72:1148–70. [8] Bouchard P, Bay F, Chastel Y, Tovena I. Crack propagation modelling using an advanced remeshing technique. Comput Methods Appl Mech Eng 2000;189:723–42. [9] Bouchard P, Bay F, Chastel Y. Numerical modelling of crack propagation: automatic remeshnig and comparison of different criteria. Comput Methods Appl Mech Eng 2003;192:3887–908. [10] Maligno A, Rajaratnam S, Leen S, Williams E. A three-dimensional (3D) numerical study of fatigue crack growth using remeshing techniques. Eng Fract Mech 2010;77:94–111. [11] Paluszny A, Zimmerman R. Numerical simulation of multiple 3D fracture propagation using arbitrary meshes. Comput Methods Appl Mech Eng 2011;200:953–66. [12] Khoei A, Eghbalian M, Moslemi H, Azadi H. Crack growth modeling via 3D automatic adaptive mesh refinement based on modified-SPR technique. Appl Math Modell 2013;37:357–83. [13] Davis B, Wawrzynek P, Ingraffea A. 3-D simulation of arbitrary crack growth using an energy-based formulation – Part I: Planar growth. Eng Fract Mech 2014;115:204–20. [14] Davis B, Wawrzynek P, Carter B, Ingraffea A. 3-D simulation of arbitrary crack growth using an energy-based formulation – Part II: Non-planar growth. Eng Fract Mech 2016;154:111–27. [15] Fracture Analysis Consultants, Inc., Ithaca, NY (USA). FRANC3D, version 7.0. ; 2016. [16] Duarte C, Babuška I, Oden J. Generalized finite element methods for three dimensional structural mechanics problems. Comput Struct 2000;77:215–32. http://dx.doi.org/10.1016/S0045-7949(99)00211-4. [17] Pereira J, Duarte C, Guoy D, Jiao X. Hp-generalized FEM and crack surface representation for non-planar 3-D cracks. Int J Numer Methods Eng 2009;77 (5):601–33. http://dx.doi.org/10.1002/nme.2419. [18] Belytschko T, Black T. Elastic crack growth in finite elements with minimal remeshing. Int J Numer Methods Eng 1999;45:601–20. [19] Sukumar N, Moës N, Moran B, Belytschko T. Extended finite element method for three-dimensional crack modelling. Int J Numer Methods Eng 2000;48 (11):1549–70. [20] Belytschko T, Gracie R, Ventura G. A review of extended/generalized finite element methods for material modeling. Modell Simulat Mater Sci Eng 2009;17:24. http://dx.doi.org/10.1088/0965-0393/17/4/043001. [21] Fries T-P, Belytschko T. The generalized/extended finite element method: an overview of the method and its applications. Int J Numer Methods Eng 2010:253–304. [22] Garzon J, O’Hara P, Duarte C, Buttlar W. Improvements of explicit crack surface representation and update within the generalized finite element method with application to three-dimensional crack coalescence. Int J Numer Methods Eng 2014;97:231–73. http://dx.doi.org/10.1002/nme.4573. [23] Chopp D, Sukumar N. Fatigue crack propagation of multiple coplanar cracks with the coupled extended finite element/fast marching method. Int J Eng Sci 2003;41:845–69. [24] Duarte C, Kim D-J. Analysis and applications of a generalized finite element method with global-local enrichment functions. Comput Methods Appl Mech Eng 2008;197(6–8):487–504. http://dx.doi.org/10.1016/j.cma.2007.08.017.

302

P. O’Hara et al. / Engineering Fracture Mechanics 163 (2016) 274–302

[25] Kim D-J, Pereira J, Duarte C. Analysis of three-dimensional fracture mechanics problems: a two-scale approach using coarse generalized FEM meshes. Int J Numer Methods Eng 2010;81(3):335–65. http://dx.doi.org/10.1002/nme.2690. [26] Hou T, Wu X-H. A multiscale finite element method for elliptic problems in composite materials and porous media. J Comput Phys 1997;134:169–89. [27] Oskay C. A variational multiscale enrichment for modeling coupled mechano-diffusion problems. Int J Numer Methods Eng http://dx.doi.org/10.1002/ nme.3258. [28] Holl M, Loehnert S, Wriggers P. An adaptive multiscale method for crack propagation and crack coalescence. Int J Numer Methods Eng 2013;93 (1):23–51. [29] Loehnert S, Belytschko T. A multiscale projection method for macro/microcrack simulations. Int J Numer Methods Eng 2007;71(12):1466–82. [30] Holl M, Rogge T, Loehnert S, Wriggers P, Rolfes R. 3D multiscale crack propagation using the XFEM applied to a gas turbine blade. Comput Mech 2014;53:5–6. [31] Belytschko T, Fish J, Bayliss A. The spectral overlay on finite elements for problems with high gradients. Comput Methods Appl Mech Eng 1990;81:71–89. [32] Fish J. The s-version of the finite element method. Comput Struct 1992;43:539–47. [33] Fish J, Nath A. Adaptive and hierarchical modelling of fatigue crack propagation. Int J Numer Methods Eng 1993;36:2825–36. [34] Rank E, Krause R. A multi-scale finite element method. Comput Struct 1997;64(1–4):139–44. [35] Krause R, Rank E. Multiscale computations with a combination of the h- and p-versions of the finite-element method. Comput Methods Appl Mech Eng 2003;192:3959–83. [36] Ghosh S, Lee K, Raghavan P. A multi-level computational model for multi-scale analysis in composite and porous materials. Int J Solids Struct 2001;38:2335–85. [37] Ghosh S, Raghavan P. Multi-scale model for damage analysis in fiber reinforced composites with interfacial debonding materials. Int J Multiscale Comput Eng 2004;2(4):621–45. [38] Lee S-H, Song J-H, Yoon Y-C, Zi G, Belytschko T. Combined extended and superimposed finite element method for cracks. Int J Numer Methods Eng 59 1119–36. http://dx.doi.org/10.1002/nme.908. [39] Schöllmann M, Richard H, Kullmer G, Fulland M. A new criterion for the prediction of crack development in multiaxially loaded structures. Int J Fatigue 2002;117:129–41. [40] Pereira J, Duarte C, Jiao X. Three-dimensional crack growth with hp-generalized finite element and face offsetting methods. Comput Mech 2010;46 (3):431–53. http://dx.doi.org/10.1007/s00466-010-0491-3. [41] Szabo BA, Babuška I. Computation of the amplitude of stress singular terms for cracks and reentrant corners. In: Cruse T, editor. Fracture mechanics: nineteenth symposium, ASTM STP 969. San Antonio (TX): Southwest Research Institute; 1988. p. 101–24. [42] Pereira J, Duarte C. Extraction of stress intensity factors from generalized finite element solutions. Eng Anal Bound Elem 2005;29:397–413. http://dx. doi.org/10.1016/j.enganabound.2004.09.007. [43] Anderson T. Fracture mechanics: fundamentals and applications. Boca Raton (FL): CRC Press; 2005. [44] Swift T. Damage tolerance capability. In: Fatigue of aircraft materials. Delft University Press; 1992. [45] Babuska I, Andersson B. The splitting method as a tool for multiple damage analysis. SIAM J Sci Comput 2005;26(4):1114–45. [46] Legrand L, Lazarus V. Front shape and loading evolution during cracks coalescence using and incremental perturbation method. Eng Fract Mech 2015;133:40–51. [47] Babuška I, Caloz G, Osborn J. Special finite element methods for a class of second order elliptic problems with rough coefficients. SIAM J Numer Anal 1994;31(4):945–81. [48] Babuška I, Melenk J. The partition of unity method. Int J Numer Methods Eng 1997;40:727–58. [49] Melenk J, Babuška I. The partition of unity finite element method: basic theory and applications. Comput Methods Appl Mech Eng 1996;139:289–314. [50] Duarte C, Oden J. Hp clouds – an hp meshless method. Numer Methods Part Diff Eq 1996;12:673–705. http://dx.doi.org/10.1002/(SICI)1098-2426 (199611)12:6<673::AID-NUM3>3.0.CO;2-P. [51] Duarte C, Oden J. An hp adaptive method using clouds. Comput Methods Appl Mech Eng 1996;139:237–62. http://dx.doi.org/10.1016/S0045-7825(96) 01085-7. [52] Becker EB, Carey GF, Oden JT. Finite elements: an introduction. Texas finite element series, vol. I. Englewood Cliffs: Prentice-Hall; 1981. [53] Oden J, Duarte C, Zienkiewicz O. A new cloud-based hp finite element method. Comput Methods Appl Mech Eng 1998;153:117–26. http://dx.doi.org/ 10.1016/S0045-7825(97)00039-X. [54] Strouboulis T, Copps K, Babuška I. The generalized finite element method. Comput Methods Appl Mech Eng 2001;190:4081–193. [55] Kim D-J, Duarte C, Proenca S. A generalized finite element method with global-local enrichment functions for confined plasticity problems. Comput Mech 2012;50(5):563–78. [56] Gupta V, Kim D-J, Duarte C. Extensions of the two-scale generalized finite element method to nonlinear fracture problems. Int J Multiscale Comput Eng 2013;11(6):581–96. [57] Kim J, Duarte C. A new generalized finite element method for two-scale simulations of propagating cohesive fractures in 3-D. Int J Numer Methods Eng 2015;104(13):1139–72. [58] O’Hara P, Duarte C, Eason T, Garzon J. Efficient analysis of transient heat transfer problems exhibiting sharp thermal gradients. Comput Mech 2013;51:743–64. http://dx.doi.org/10.1007/s00466-012-0750-6. [59] Dompierre J, Labb P, Vallet M-G, Camarero R. How to subdivide pyramids, prisms and hexahedra into tetrahedra, rapport CERCA R99-78; 1999. [60] Arnold DN, Mukherjee A, Pouly L. Locally adapted tetrahedral meshes using bisection. SIAM J Sci Comput 2000;22(2):431–48. [61] Bansch E. Local mesh refinement in 2 and 3 dimensions. Impact Comput Sci Eng 1991;3:181–91. [62] Whitcomb J. Iterative global/local finite element analysis. Comput Struct 1991;40:1027–31. [63] Gendre L, Allix O, Gosselet P, Comte F. Non-intrusive and exact global/local techniques for structural problems with local plasticity. Comput Mech 2009;44:233–45. [64] Rose L. Microcrack interaction with main crack. Int J Fract 1986;31:233–42. [65] Garzon J, Duarte C, Pereira J. Extraction of stress intensity factors for simulations of 3-D crack growth with the generalized finite element method. Key Eng Mater 2013;560:1–36. special issue on Advances in Crack Growth Modelling. [66] Loehnert S, Belytschko T. Crack shielding and amplification due to microcracks interacting with a macrocrack. Int J Fract 2007;145:1–8. [67] Rubinstein A. Macrocrack interaction with semi-infinite microcrack array. Int J Fract 1985;27:113–9. [68] Byrd P, Friedman M. Handbook of elliptic integrals. 2nd ed. Berlin: Springer Verlag; 1971. [69] Abramowitz M, Stegun I. Handbook of mathematical functions. New York: Dover Publications, Inc.; 1972. [70] O’Hara P, Hollkamp J. Modeling vibratory damage with reduced-order models and the generalized finite element method. J Sound Vib 2014;333:6637–50. [71] Melin S. Why do cracks avoid each other? Int J Fract 1983;23:37–45. [72] Swain M, Hagan J. Some observations of overlapping interacting cracks. Eng Fract Mech 1978;10:299–304. [73] Cotterell B, Rice J. Slightly curved or kinked cracks. Int J Fract 1980;16(2):155–69. [74] McComb T, Pope J, Grant A. Growth and coalescence of multiple fatigue cracks in polycarbonate test specimens. Eng Fract Mech 1986;24(4):601–8. [75] Kamaya M. Growth evaluation of multiple interacting surface cracks. Part I: Experiments and simulation of coalesced crack. Eng Fract Mech 2008;75:1336–49. [76] Chang R. On crack-crack interaction and coalescence in fatigue. Eng Fract Mech 1982;16(5):683–93.