A homogenized XFEM approach to simulate fatigue crack growth problems

Computers and Structures 150 (2015) 1–22

Contents lists available at ScienceDirect

Computers and Structures journal homepage: www.elsevier.com/locate/compstruc

A homogenized XFEM approach to simulate fatigue crack growth problems Sachin Kumar, I.V. Singh ⇑, B.K. Mishra Department of Mechanical and Industrial Engineering, Indian Institute of Technology Roorkee, Uttarakhand, India

a r t i c l e

i n f o

Article history: Received 13 March 2014 Accepted 11 December 2014

Keywords: Edge crack plate Homogenized XFEM Local mesh refinement Discontinuities Fatigue life CPU time

a b s t r a c t In this work, a homogenized XFEM approach has been proposed to evaluate the fatigue life of an edge crack plate in the presence of discontinuities (holes, inclusions and minor cracks). In this approach, the central 20% region of the plate is discretized by fine-mesh whereas remaining 80% region is discretized by coarse mesh. In case of holes and inclusion, the 80% continuum region of the plate has been modeled with homogeneous material properties. Special transition elements are developed to maintain the displacement continuity at the junction nodes. Several problems containing discontinuities in 20% region are solved by the proposed approach. Ó 2014 Elsevier Ltd. All rights reserved.

1. Introduction In order to verify the tolerance of the structures to fatigue damage, one needs to investigate the effects of the existing flaws on the structures. In general, all natural or manmade materials are heterogeneous at microscopic level. Typical examples are metal alloys, porous and cracked media, polycrystalline materials and composites. From the time and cost point of view, macroscopic analysis is performed to determine the overall characteristics of the heterogeneous materials. However, the heterogeneous nature of the material attributes to the crack initiation, propagation and coalescence of the micro-cracks within the material. Thus, the microscopic analysis of static and growing cracks becomes quite important for the structures/components subjected to fatigue loading. Many components such as aircraft engine, turbine blades and rotating shafts subjected to cyclic loads may sometimes lead to the fatigue failure. The fatigue life of a component is highly influenced by the crack closure during cyclic loading. Thus, the design analysis of the structures subjected to fatigue loading needs to be performed effectively and accurately. To consider the effect of the heterogeneity of the material, the macroscopic analysis can be performed by replacing the actual heterogeneous material by an equivalent homogenous material. This process is called effective medium approximation (EMA) or ⇑ Corresponding author. Tel.: +91 1332 285888; fax: +91 1332 285665. E-mail address: [email protected] (I.V. Singh). http://dx.doi.org/10.1016/j.compstruc.2014.12.008 0045-7949/Ó 2014 Elsevier Ltd. All rights reserved.

homogenization which involves finding out the properties of the equivalent homogenous material. So far, several homogenization approaches have been proposed such as Mori–Tanaka model [39], self-consistent approach [23,15], differential scheme [35]. Some other numerical methods such as asymptotic method [64], multipole expansion method [29] and fast Fourier transform (FFT) based approach [41,37,6]. Recently many researchers [73,55,1] have used EMA approach to approximate the composite behavior. This approach is very effective from time and cost point of view but the effect of micro-cracks and flaws cannot be obtained exactly on fracture and fatigue behavior. To model the crack and crack growth behavior, several methods are available such as finite element method [13], boundary element method [77,78], meshfree methods [3,4,17] and extended finite element method [5,38,67]. The extended finite element method (XFEM), based on the partition of unity [33] has been extensively used for the simulation of crack growth problems. In this method, a crack can be modeled independent of finite element mesh. In XFEM, standard finite element approximation is enriched locally with some additional functions [5], which are obtained from the theoretical background of the problem. The level set method proposed by Osher and Sethian [47] was first applied in XFEM by Stolarska et al. [66] and Sukumar et al. [67] to track a moving discontinuity. To reproduce the partition of unity in the blending elements, a ramp function concept was introduced by Fries [19] for two-dimensional problems, and then extended for three-dimensional problems by Loehnert et al. [32].

2

S. Kumar et al. / Computers and Structures 150 (2015) 1–22

The advantage of XFEM is quite clear in crack growth problems because the level set method can be easily modified to determine the new crack location and geometry. Daux et al. [16] applied the extended finite element method to model the voids, cracks with multiple branches, and crack emerging from the holes. Zi and Belytschko [79] modeled the elements intersected by the crack (split and tip elements) with the sign function only. Sukumar and Prevost [68] carried out the quasi-static crack growth analysis without remeshing in homogeneous and bi-material media using XFEM. Loehnert and Belytschko [30] performed the numerical simulations of interaction between micro cracks with a macro crack using XFEM. They concluded that for practical problems, XFEM leads to better results than analytical methods. In 2006, Bordas and Moran [7] extended XFEM to simulate the out-of-plane propagation of two cracks in an aeroengine component under combined mechanical and thermal loading. Further, Wyart et al. [75] implemented XFEM in the finite element commercial codes to simulate the three-dimensional crack growth problems. In this approach, the cracked and un-cracked domains were modeled by XFEM and FEM respectively. In 2012, Singh et al. [65] used XFEM to evaluate the fatigue life of homogenous plate in the presence of multiple discontinuities under cyclic load. In 2013, various bi-material interface crack problems were modeled by XFEM [11,52]. A direct numerical simulation of the micro-scale problems with uniform mesh is not viable choice even with the modern supercomputers. The major difficulty in direct solution is the requirement of very fine mesh to model the discontinuities, which results in the increase of degree of freedoms, and hence requires a tremendous amount of computer memory and CPU time. The situation can be alleviated to some extends by parallel computing, but the size of the problem cannot be reduced. To avoid this problem, a multiscale approach is proposed, which reduces the numerical effort without compromising on the accuracy. Several different concepts of the multiscale approach have been proposed in literature i.e. homogenized Dirichlet projection method [80,46], variational multiscale methods [25], domain decomposition methods [18] and homogenization based multiscale methods [72]. Guidault et al. [20] proposed a multiscale strategy, based on domain decomposition method, for the modeling of micro cracks. In 2007, Loehnert and Belytschko [31] proposed a multiscale method based on two scale decomposition of the displacements. They observed the shielding effect of micro cracks over the macro cracks in the vicinity of macro cracks. A generalized finite element method based on two scale decomposition approaches is used for the simulation of fracture problems by Kim et al. [26] and Pereira et al. [53]. Holl et al. [24] proposed an adaptive multiscale technique for the crack propagation and crack coalescence of macro cracks and micro cracks. In 2009, Rannou et al. [60] proposed a local multigrid XFEM strategy for the simulation of crack growth in three-dimensional objects. This strategy provides the higher computational efficiency due to local mesh refinement near the cracks, inclusions and steep gradients, etc. Pierres et al. [54] simulated three-dimensional fatigue crack growth problems of interfacial contact using two-scale XFEM. This concept was further extended to evaluate the stress intensity factors using a three-scale concurrent multigrid XFEM [50]. Further, an alternative algorithm, based on a localized multigrid algorithm [50,51] was proposed to simulate the mixed-mode crack propagation. In 2008, Nguyen et al. [45] presented a review on meshless methods, where a detailed explanation of a simple element-free Galerkin MATLAB code is provided. In this article, a comparison between local (intrinsic or extrinsic) and global enrichment is also provided for simple one-dimensional problem. An h-adaptivity with a posteriori error estimation is applied in meshfree particle method by Rabczuk and Belytschko [57]. In this work, the initial particle arrangement is structured along with a background mesh,

and interior interfaces and outside boundaries are described by the implicit functions. Gravouil et al. [22] is the first who implemented XFEM with level set for three-dimensional linear elastic fracture mechanics problems. Further, some articles [58,59,9,61] describe the application of meshfree methods to three-dimensional crack growth problems. The strain smoothing technique proposed by Chen and Wang [12] for meshless methods developed the smoothed finite element method (SFEM). Further, a numerical analysis was performed by Bordas et al. [10] to obtain the performance of strain smoothing for quadratic enriched FE approximations. They observed that for higher order elements it fails the patch test but improves the results of polynomial enrichment functions. In enriched elements (split and tip elements), special care was required for numerical integration. To avoid this, a technique based on replacing nonpolynomial functions by equivalent polynomials was proposed by Ventura [74]. Thereafter, a recent polygonal integration technique was also proved to be an effective alternative [42,43]. In 2011, Menk and Bordas [36] observed that the stiffness matrix may become ill-conditioned by using additional enrichment functions in XFEM. Therefore, to avoid this problem, a robust preconditioning technique was described to obtain the stiffness matrices whose condition number was close to the FE. This technique provides well-conditioned matrices for any type of enrichment. Natarajan and Song [44] extended SBFEM in XFEM and replaced the asymptotic enrichments around the crack tip with the semianalytical solution obtained by SBFEM. Hence, special numerical integration technique was not required in XFEM, which improves the capability of the XFEM to model the cracks in homogeneous/ heterogeneous materials. The main purpose of the adaptive refinement strategy is to adjust the mesh resolution in order to get the better accuracy. In 1978, Gupta [21] developed a set of conforming interpolation functions that maintain the inter-element compatibility for twodimensional analysis. Later on the concept of transition elements is extended for the simulation of three-dimensional problems by McDill et al. [34], Choi and Lee [14] and Morton et al. [40]. Today, the most advanced refinement technique is quadtree data structure, which is based on the principle of recursive decomposition [63]. This technique provides a simple and efficient way for the h-refinement and p-refinement. In this technique, hanging nodes are generated and temporary triangular or rectangular elements are added to those elements which have hanging nodes to get the compatible meshes [48,56]. In 2005, Tabarraei and Sukumar [70] developed an h-adaptive finite element method based on quadtree data structure with conforming polygon interpolants. Wu et al. [76] developed new enhanced assumed strain and hybrid stress transition element families for the refinement of two and three-dimensional elastic problems. Sze and Wu [69] also introduced the transition elements which comprise five to seven node quadrilateral elements for adaptive analysis of axisymmetric elasticity problems. This paper is organized as follows: first the governing equations are given for the linear elastic fracture mechanics. A review of crack growth criterion and fatigue life estimation is presented in Section 2. In Section 3, the parametric study is performed to define the region for discontinuities (holes, inclusions and minor cracks). The homogenization process and the evaluation of the homogenized material properties are described in Section 4. The transition elements shape functions are derived in Section 5 to combine the coarse and fine mesh regions. In Section 6, the proposed homogenized XFEM approach is applied to several crack interaction and fatigue crack growth problems. The results obtained by proposed method are compared with the results in available literature and standard XFEM. Section 7 provides the conclusions of the present work.

3

S. Kumar et al. / Computers and Structures 150 (2015) 1–22

2. Numerical formulation

2.2. XFEM formulation

2.1. Governing equations

In 2-D, at a particular node of interest xi, the displacement approximation with discontinuities such as cracks, holes and inclusions can be written as [5,67,36,27],

Consider a domain (X) shown in Fig. 1. The boundary of domain is partitioned into displacement (Cu), traction (Ct) and crack surface (Cc) boundaries. The strong form of the equilibrium equations and boundary conditions can be written as [67,38,65],

2 n X 6 uh ðxÞ ¼ Ni ðxÞ4u i þ ½HðxÞ  Hðxi Þai þ ½uðxÞ  uðxi Þbi |fflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflffl} |fflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflffl} i¼1 i2nr

r  r þ b ¼ 0 in X

ð1aÞ

r  n^ ¼ t on Ct : external traction r  n^ ¼ 0 on Cc : traction free crack

ð1bÞ

 on Cu : prescribed displacement u¼u

ð1dÞ

ð1cÞ

where r is the Cauchy stress tensor, u is the displacement field vec^ is the unit outward normal and b and t are the body force and tor, n external traction vector. For small strains and displacements, strain–displacement relation can be written as,

e ¼ eðuÞ ¼ rs u

ð2Þ

In the above equation rs is the symmetric part of the gradient operator and e is the linear strain tensor. The constitutive relation for linear elastic material is given by Hook’s law,

r¼De

ð3Þ

where D is the Hook’s tensor. The variational form of the equilibrium equation can be written as,

Z

rðuÞ : eðv ÞdX ¼

X

Z X

b  v dX þ

Z

t  v dC

ð4Þ

Ct

After substituting the trial and test functions and using the arbitrariness of nodal variations, the following discrete system of equations are obtained,

½Kfdg ¼ ffg

ð5Þ

i2ni

8 > >
9 3 > > = 7 a þ ½wðxÞ  wðxi Þci þ ½b ðxÞ  ba ðxi Þdi RðxÞ7 5 a¼1 a |fflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflffl} > > > :|fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl}> ; i2n

ð6Þ

i2nA

h

 i is a nodal displacement vector of the standard finite elewhere u ment part at node i; n is the set of all nodes in the mesh; nr is the set of nodes associated with those elements which are completely cut by the crack; ni is the set of nodes associated with those elements which are cut by the inclusions and material interfaces; nh is the set of nodes associated with those elements which are cut by the hole; nA is the set of nodes associated with those elements which are partially cut by the crack. R(x) is a ramp function [19], which is used to model the blending elements and is defined as, P RðxÞ ¼ i2nenr Ni ðxÞ. nenr is the set of crack tip enriched nodes. A typical discretized geometry as illustrated in Fig. 2 shows all types of nodes, where, H(x) is the discontinuous enrichment function or Heaviside function, defined for those elements which are completely cut by the cracks, and takes the value +1 on one side of crack and 1 on other side of the crack; ai is the nodal enriched degree of freedom associated with H(x). bi is the nodal enriched degree of freedom associated with u(x), where u(x) is level set function. ci is the additional nodal degrees of freedom associated with w(x), where w(x) takes a value +1 outside the hole and zero inside the a hole. di is the nodal enriched degrees of freedom vector associated with crack tip enrichment function, ba(x). Considering a local polar coordinates system r and h at the crack tip, the asymptotic crack tip enrichment functions can be written as [8],

ba ðxÞ ¼

 pffiffiffi  pffiffiffi pffiffiffi pffiffiffi r sin 2h ; r cos 2h ; r sin 2h cos h; r cos 2h cos h ð7Þ

Using the approximation function defined in Eq. (6), the elemental matrices, k and f can be expressed as,

where K is the global stiffness matrix, d is the vector of nodal unknowns (both standard and enriched) and f is the external force vector.

Split nodes of inclusion/hole

Inclusion/hole Blending element Split nodes of crack Tip nodes

Tip element

Fig. 1. Domain with a crack and prescribed boundary conditions.

Crack

Split elements

Fig. 2. Typical discretization of the 2-D domain with crack and inclusions/hole

4

S. Kumar et al. / Computers and Structures 150 (2015) 1–22

2

uu

ua

k 6 ij 6 au 6 kij 6 6 e kij ¼ 6 kbu 6 ij 6 cu 6 kij 4 du kij

ub

3 ud

uc

kij

kij

kij

kij

aa kij

ab kij

ac kij

ad kij

ba

bb

bc

kij

kij

kij

ca kij

cb kij

cc kij

da

db

kij

kij

7 7 7 7 bd 7 kij 7 7 cd 7 kij 7 5 dd kij

dc

kij

ð8Þ

The sub-matrices and vectors that comes in the foregoing equations are defined as,

n

e

fi ¼ rs kij

¼

u

fi ¼ a

fi ¼ b fi

¼

c

fi ¼ da

u

b

c

fi

fi

d1

d2

fi

d3

fi

fi

d4

fi

oT

ð9Þ

Xe

Xe

Z Xe

Z Xe

Z Xe

T ðBri Þ D Bsj hd

Z

Ni b dX þ

X where; r; s ¼ u; a; b; c; d

ð10Þ

Nit dC

ð11Þ

Ct

Ni ðHðxÞ  Hðxi ÞÞb dX þ

Z

Ni ðHðxÞ  Hðxi ÞÞ t dC

ð12Þ

Ct

Ni ðuðxÞ  uðxi ÞÞb dX þ

Z

Ni ðuðxÞ  uðxi ÞÞt dC

Ni ðwðxÞ  wðxi ÞÞb dX þ

Ni ðwðxÞ  wðxi ÞÞt dC

ð14Þ

Ct

Z Xe

Ni ba ððxÞ  ðxi ÞÞb dX þ

Z Ct

2

Ni;x

6 ¼4 0 Ni;y 2

6 Bai ¼ 4 2 6 Bbi ¼ 4 2 6 Bci ¼ 4

0

ð15Þ

3 ð16aÞ

ðNi ðHðxÞ  Hðxi ÞÞÞ;x

ðN i ðHðxÞ  Hðxi ÞÞÞ;y 7 5

0 ðNi ðHðxÞ  Hðxi ÞÞÞ;y

Ni ðuðxÞ  uðxi ÞÞ;y 7 5 ðNi ðuðxÞ  uðxi ÞÞÞ;x

ðNi ðuðxÞ  uðxi ÞÞÞ;y

ðNi ðwðxÞ  wðxi ÞÞÞ;y

2 6 ¼4

Bd4 i

ð16dÞ

ðNi ðwðxÞ  wðxi ÞÞÞ;x

i

ðNi ðba ðxÞ  ba ðxi ÞÞÞ;x 0 ðNi ðba ðxÞ  ba ðxi ÞÞÞ;y

ð16eÞ 0

3

ðNi ðba ðxÞ  ba ðxi ÞÞÞ;y 7 5 ðNi ðba ðxÞ  ba ðxi ÞÞÞ;x

 1  ð1Þ ð2Þ ð1Þ ð1Þ ð2Þ ð2Þ ð1Þ rij eij þ rð2Þ ¼ rij eij ¼ rij eij ij eij 2

ð17Þ

ð18Þ

where rij and eij are the stress and strain fields respectively, 1 and 2 signify the actual state and auxiliary state respectively. The SIF values are calculated from following equation by using the Eq. (17),

Ið1;2Þ ¼

 2  ð1Þ ð2Þ ð1Þ ð2Þ 0 K I K I þ K II K II E

ð19Þ

In this section, the criterion for the crack growth direction is briefly described along with the fatigue life computation. The fatigue crack growth simulations are performed by XFEM under constant amplitude cyclic loading. The discrete equations are solved to obtain the displacements, and then values of stress intensity factors are extracted from Eq. (19). The range of SIF for constant amplitude cyclic loading is defined as,

DK ¼ K max  K min

ðNi ðwðxÞ  wðxi ÞÞÞ;y 7 5

0

Bd3 i

ð16cÞ

3

0

ðNi ðwðxÞ  wðxi ÞÞÞ;x

Bd2 i

3

0

0

h Bdi ¼ Bd1 i

ð16bÞ

ðNi ðHðxÞ  Hðxi ÞÞÞ;x

ðNi ðuðxÞ  uðxi ÞÞÞ;x

# ð2Þ ð1Þ @ui @q ð2Þ @ui þ rij  W ð1;2Þ d1j dA @xj @x1 @x1

2.5. Crack growth direction and fatigue life

3

0

rijð1Þ

where q is a smooth weighting function, used to convert the contour integral form into equivalent domain form, which has value one at the crack tip and zero on the contour. W(1,2) is the interaction strain energy and is given as,

W ð1;2Þ ¼

Ni ðba ðxÞ  ðxi ÞÞt dC;

7 Ni;y 5 N i;x

Z "

ð13Þ

Ct

Z

In this section, the domain form of the interaction integral is described for the extraction of individual stress intensity factors KI and KII. The interaction integral is an effective technique for extracting the mixed-mode stress intensity factors [16,62,65]. For two independent equilibrium states of a cracked body, the domain form of interaction integral can be written as,

A

where Ni are finite element shape functions, Bui ; Bai ; Bbi ; Bci and Bdi are the matrices of shape function derivatives and are given below,

Bdi a

Tip element: 7 Gauss points in each sub-triangle. Split element: 5 Gauss points in each sub-triangle. Tip blending element: 16 Gauss points. Standard element: 4 Gauss points.

Ið1;2Þ ¼

where; a ¼ 1; 2; 3; 4

Bui

   

2.4. Evaluation of stress intensity factors

Z

Z

fi ¼

a

fi

fi

effectively. Thus, these partitioned elements are divided into sub-triangles [65,52] for the purpose of integration. In these subtriangles, the integrand remains continuous. As the displacement approximation is different from element to element therefore, the following integration orders are adopted in the various elements,

ð16fÞ

where; a ¼ 1; 2; 3; 4

where Kmax and Kmin are the stress intensity factors corresponding to maximum and minimum applied loads respectively. In real life problems, a crack path is curved in nature. Therefore, it is modeled through many small line segments. To determine the direction of crack growth, several methods have been developed. Maximum principal stress criterion is employed to determine the crack growth direction, which employs that the crack propagates in the direction of maximum circumferential stress. The equivalent mode-I SIF and the direction of crack growth hc at each crack increment is obtained by the following expression [52,28],

   hc hc hc  3DK II cos2 sin 2 2 2 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi1 0 2 2 K I  K I þ 8K II A hc ¼ 2 tan1 @ 4K II

DK Ieq ¼ DK I cos3 2.3. Integration schemes for discontinuous elements For the elements partitioned by the discontinuities, the standard Gauss quadrature rule cannot execute the integration

ð20Þ

ð21Þ

ð22Þ

5

S. Kumar et al. / Computers and Structures 150 (2015) 1–22

a

H

H

H

L

L

L

(a)

(b)

(c)

H

H

H

L

L

L

(d)

(e)

(f)

Fig. 3. Plate with edge crack and holes for different cases: (a) without holes (b) holes in 20% region, (c) holes in 30% region, (d) holes in 40% region, (e) holes in 60% region, (f) holes in 100% region.

For quasi-static analysis, the fatigue life is obtained using the generalized Paris law [49,17], which can be written as,

da ¼ CðDK Ieq Þm dN

ð23Þ

where a is the crack length, N is number of loading cycles, C and m are material properties. The failure takes place whenever (KIeq)max > KIC, where (KIeq)max is the equivalent SIF corresponding to maximum load and KIC is the fracture toughness of the material. At the time of failure, the number of cycles elapsed gives the fatigue life of the material. 3. Region selection for modeling discontinuities In order to compute the effect of small size discontinuities over the major crack, different cases of edge crack rectangular plate in

the presence of multiple discontinuities are simulated under the prescribed traction and boundary conditions. To capture the effect of discontinuities over the major crack using single scale FEM or XFEM is either numerically inaccurate with a coarse mesh or computationally expensive with a fine mesh. Therefore, to avoid this problem, a small region near the main crack is modeled with discontinuities (holes, inclusions and minor cracks). These discontinuities are taken in such a way that they do not very close to the main crack tip. To fix the region for discontinuities, several problems of edge crack having discontinuities in 20%, 30%, 40%, 60% and 100% of the domain are simulated by XFEM. A rectangular plate of size 100 mm  200 mm along with an edge crack of length a = 20 mm is taken for the simulation as shown in Fig. 3a. Several cases of randomly distributed multiple discontinuities are simulated under plane stress condition. In single scale XFEM, a very fine mesh is required to model the problem accurately, hence the domain is discretized with a uniform mesh

6

S. Kumar et al. / Computers and Structures 150 (2015) 1–22

1100 Holes in 100% Region Holes in 60% Region Holes in 40% Region Holes in 30% Region Holes in 20% Region Without Holes

SIF (MPa.mm1/2)

1000

900

800

700

600

500 18

20

22

24

26

28

30

32

34

36

Crack length (mm) Fig. 4. SIF variation with crack length for an edge crack plate in the presence of multiple holes.

Table 1 Variation of SIF and CPU time with holes modeling region for a crack length of 20 mm. S. no.

Holes modeling region (%)

Stress intensity factor SIF (MPa mm 2 )

% Difference

1 2 3 4 5 6

Without holes 20 30 40 60 100

539.7 565.4 580.7 593.8 606.3 612.7

11.9 7.7 5.2 3.1 1.0 –

1/

CPU time (s)

215 274 304 334 379 491

size of 88 by 176 nodes in x and y directions respectively. An increment of 2 mm is given at each step of crack growth until the final failure i.e. when (KIeq)max > KIC. The plate is subjected to a tensile load of rmin = 0 MPa and rmax = 50 MPa at the top edge. The bottom edge of the plate is constrained in the y-direction. The aluminum alloy 6061-T6, is used for all simulation and properties are taken from [2].

Elastic modulus E ¼ 68:9 GPa Elastic modulus for inclusions EI ¼ 20 GPa Poisson ratio for material Yield strength

m ¼ 0:33

0 y

r ¼ 276 MPa mI ¼ 0:33

Poisson ratio for inclusions

pffiffiffiffiffi Fracture toughness K IC ¼ 29 MPa m

Paris exponent m ¼ 3:17 Paris constant C ¼ 5:88  108 In the first problem, multiple holes of radii varying from 2.0 mm to 3.5 mm are randomly distributed in the plate. Five different cases of randomly distributed holes are used to define the region of discontinuities near the major crack i.e. 20, 30, 40, 60 and 100 holes are randomly distributed in 20% (Fig. 3b), 30% (Fig. 3c), 40% (Fig. 3d), 60% (Fig. 3e) and 100% (Fig. 3f) portions of the plate respectively. In all five cases, the volume fraction (10%) of the holes with respect to the region in which they are distributed is kept constant. The results are compiled 10 times to obtain the average effect of randomly distributed holes on the major crack. The variation of stress intensity factor with crack length for different cases is shown in

Fig. 4. For an initial crack length of 20 mm, the values of stress intensity factor without holes and holes in 100% region of the total domain are found as 539.7 MPa mm1/2 and 612.7 MPa mm1/2 respectively. The value of SIF becomes 565.4 MPa mm1/2 for holes in 20% region, 580.7 MPa mm1/2 for holes in 30% region, 593.8 MPa mm1/2 for holes in 40% region and 606.3 MPa mm1/2 for holes in 60% region. From the results presented in Fig. 4, it is observed that the percentage difference in SIF (with respect to holes in 100% region) is found as 7.7%, 5.2%, 3.1% and 1.0% for holes in 20%, 30%, 40% and 60% region respectively. Table 1 presents the CPU time required to solve these problems having crack length of 20 mm. The CPU time required without holes is found as 215 s, whereas it is increased to 274 s, 304 s, 334 s, 379 s and 491 s for holes in 20%, 30%, 40%, 60% and 100% regions respectively. Although, the CPU time significantly reduces if the holes are modeled in 60%, 40%, 30% and 20% domain, but their effect on SIF is not be properly captured in these cases. In the second problem, multiple inclusions of different radii varying from 2.0 mm to 3.5 mm are randomly distributed in the plate. Five different cases of randomly distributed inclusions are used to fix the region of discontinuity near the major crack i.e. 20, 30, 40, 60 and 100 inclusions are randomly distributed in 20% (Fig. 5a), 30% (Fig. 5b), 40% (Fig. 5c), 60% (Fig. 5d) and 100% (Fig. 5e) region of the plate respectively. In all five cases, the volume fraction (10%) of the inclusions is kept constant. The results are compiled ten times to get an average effect of randomly distributed inclusions on the major crack. The variation of stress intensity factor with crack length for different cases is shown in Fig. 6. For an initial crack length (a) of 20 mm, the values of stress intensity factor without inclusions and inclusions in 100% region of the total domain are found as 539.7 MPa mm1/2 and 572.4 MPa mm1/2 respectively. The SIF becomes 550.5 MPa mm1/2 for inclusions in 20% region, 557.6 MPa mm1/2 for inclusions in 30% region, 563.7 MPa mm1/2 for inclusions in 40% region and 569.9 MPa mm1/2 for inclusions in 60% region. From the results presented in Fig. 6, it is observed that the percentage difference in SIF (with respect to inclusions in 100% region) is found as 3.8%, 2.6%, 1.5% and 0.4% for inclusions in 20%, 30%, 40% and 60% region respectively. Table 2 presents CPU time required to solve these problems with an initial crack length of 20 mm. The CPU time without inclusions is found as 215 s, whereas it is increased to 333 s, 421 s, 532 s, 821 s and 1450 s for inclusions in 20%, 30%, 40%, 60% and 100% regions respectively. Even though, the CPU time significantly reduces if the inclusions are modeled in 60%, 40%, 30% and 20% region, but their effect on SIF cannot be properly captured in these cases. In the third problem, multiple minor cracks of different sizes varying from 3.0 mm to 4.5 mm are randomly distributed in the plate. Again, to define the region for modeling discontinuity near the major crack, five different cases of randomly distributed minor cracks are used i.e. 16, 24, 32, 48 and 80 minor cracks are randomly distributed in 20% (Fig. 7a), 30% (Fig. 7b), 40% (Fig. 7c), 60% (Fig. 7d) and 100% (Fig. 7e) region of the plate respectively. The results are compiled ten times to obtain the average effect of randomly distributed minor cracks on the major crack. The variation of stress intensity factor with crack length for different cases is shown in Fig. 8. For an initial crack length of 20 mm, the values of SIF without minor cracks and minor cracks in 100% region of the total domain are found as 539.7 MPa mm1/2 and 548.9 MPa mm1/2 respectively. Furthermore, the value of stress intensity factor becomes 544.9 MPa mm1/2 for minor cracks in 20% region, 545.8 MPa mm1/2 for minor cracks in 30% region, 547.4 MPa mm1/2 for minor cracks in 40% region and 548.4 MPa mm1/2 for minor cracks in 60% region of the domain. From the results presented in Fig. 8 it can be observed that the percentage difference in SIF (with respect to the case of minor cracks

7

S. Kumar et al. / Computers and Structures 150 (2015) 1–22

H

H

H

L

L

L

(a)

(b)

(c)

H

H

L

L

(d)

(e)

Fig. 5. Plate with edge crack and inclusions for different cases: (a) inclusions in 20% region, (b) inclusions in 30% region, (c) inclusions in 40% region, (d) inclusions in 60% region, (e) inclusions in 100% region.

in 100% region) is found as 0.7%, 0.6%, 0.3% and 0.1% for minor cracks in 20%, 30%, 40% and 60% regions respectively. Table 3 presents CPU time required to solve these problems for a crack length of 20 mm. The CPU time without minor cracks is found as 215 s, whereas it is increased to 224 s, 229 s, 234 s, 243 s and 271 s for minor cracks in the 20%, 30%, 40%, 60% and 100% region respectively. Thus, the CPU time decreases with the decrease in minor cracks modeling region (from 60% to 20%) without much difference in the SIF values. From the simulations of above three problems, it can be seen that with the decrease in the region of discontinuities from 100% to 20%, there is certainly a substantial reduction in the CPU time but the difference in the results is also quite significant for holes and inclusions. Therefore, to improve the accuracy of the simulated results, a homogenization approach based on strain energy density is used in the present work for modeling holes and inclusions. As per this approach, first the effective properties of the plate in the

presence of the discontinuities will be evaluated using representative volume element (RVE) approach. The properties thus evaluated will be used as input for the continuum region (where holes and inclusions are not modeled) of the plate, and the region near the main crack will be modeled with discontinuities (holes and inclusions). The results obtained by this homogenization approach will be compared with those obtained by modeling discontinuities in 100% region of the domain. Thus, the next section presents a scheme for evaluating the effective properties of the continuum region. 4. Evaluation of effective properties All real materials used for engineering applications are essentially heterogeneous in nature. Their structure/property relations are complex and time consuming, hence difficult to analyze them analytically. To analyze the heterogeneous material, first

8

S. Kumar et al. / Computers and Structures 150 (2015) 1–22

 Pure in-plane shear load (with specified boundary conditions as shown in Fig. 10b) is used to evaluate the effective shear modulus (Geff).  Poisson’s ratio can be expressed in terms of shear and Young’s moduli.

1050 Inclusions in 100% Region Inclusions in 60% Region Inclusions in 40% Region Inclusions in 30% Region Inclusions in 20% Region Without Inclusions

1000 950

SIF (MPa.mm1/2)

900 850 800 750 700 650 600 550 500 18

20

22

24

26

28

30

32

34

36

Crack length (mm) Fig. 6. SIF variation with crack length for an edge crack plate in the presence of multiple inclusions.

Table 2 Variation of SIF and CPU time with inclusions modeling region for a crack length of 20 mm. S. no.

Inclusions modeling region (%)

Stress intensity factor SIF (MPa mm1/2)

% Difference

1 2 3 4 5 6

Without inclusions 20 30 40 60 100

539.7 550.5 557.6 563.7 569.9 572.4

5.7 3.8 2.6 1.5 0.4 –

Now as per the requirement, identify the boundary conditions from the above explained cases and obtain the corresponding material constant. To evaluate the other constant, choose other set of boundary conditions. An algorithm to calculate the properties of the equivalent homogenous material is given in the Flowchart. In the present work, Young’s modulus is taken different for matrix and inclusions while the Poisson’s ratio is taken same, hence one set of loading and boundary conditions is sufficient to obtain the Young’s modulus of the equivalent homogenous material. A two-dimensional plane stress condition is assumed for the simulations. The plate is subjected to uniform tensile strain as shown in Fig. 10a. Since, the strain energy density is a function of Young’s modulus only, thus Young’s modulus of equivalent homogenous material can be obtained by solving the equations. The effective properties for holes, inclusions and combined (holes and inclusions) case are obtained using strain energy density approach. The effective Young’s modulus in the presence of 10% volume fraction of discontinuities is found as 52.4 GPa for holes, 60.8 GPa for inclusions and 56.5 GPa for holes and inclusions.

CPU time (s)

215 333 421 532 821 1450

homogenization of the material is performed. In this process, a heterogeneous material is replaced by an equivalent homogenous material. It involves finding out the mechanical properties such as Young’s modulus and Poisson’s ratio of equivalent homogenous material. A schematic diagram for the homogenization process is shown in Fig. 9. So far, several homogenization processes have been proposed by many investigators. In this section, the principle of energy equivalence is used to obtain the mechanical properties of the equivalent materials. The strain energy density (SED) of the actual heterogeneous material under prescribed loading and boundary conditions must be equal to strain energy density of the equivalent homogenous material under the same loading and boundary conditions. Extended finite element method (XFEM) is used to obtain the strain energy density of the heterogeneous materials while, the strain energy density of the equivalent homogenous material can be obtained theoretically in terms of unknown materials properties. The strain energy density of heterogeneous material is equated with the strain energy density of the equivalent homogenous material to obtain the unknown properties of equivalent homogenous material. The equivalent homogenous material is assumed to be isotropic in nature, therefore can be defined by two material constants, e.g. Young’s modulus and Poisson’s ratio. To obtain these constants, the simulations are performed for two types of loading and boundary conditions as explained below:  Pure uniaxial tension load (with specified boundary conditions as illustrated in Fig. 10a) is used to evaluate the effective Young’s modulus (Eeff).

Using effective properties, the variation of stress intensity factor with crack length for different cases of multiple holes are shown in Fig. 11 and the values of SIF for initial crack length are provided in Table 4. The value of stress intensity factor for a crack length of 20 mm is found as 597.4 MPa mm1/2 for holes in 20% region, 602.0 MPa mm1/2 for holes in 30% region, 607.5 MPa mm1/2 for holes in 40% region, 611.8 MPa mm1/2 for holes in 60% region and 612.7 MPa mm1/2 for holes in 100% region of the total domain. From the results presented in Fig. 11, it is observed that the percentage difference in SIF (with respect to holes in 100% region) is

9

S. Kumar et al. / Computers and Structures 150 (2015) 1–22

H

H

H

L

L

L

(a)

(b)

(c)

H

H

L

L

(d)

(e)

Fig. 7. Plate with edge crack and minor cracks for different cases: (a) minor cracks in 20% region, (b) minor cracks in 30% region, (c) minor cracks in 40% region, (d) minor cracks in 60% region, (e) minor cracks in 100% region.

2.5%, 1.7%, 0.8% and 0.1% for holes in 20%, 30%, 40% and 60% region respectively. The simulations show that the use of effective properties in the continuum region improves the results drastically. Now using effective properties the variation of stress intensity factor with crack length for different cases of multiple inclusions are shown in Fig. 12 and the values of SIF for a crack length of 20 mm are provided in Table 5. The value of stress intensity factor is found as 564.7 MPa mm1/2, 567.1 MPa mm1/2, 569.9 MPa mm1/2, 572.4 MPa mm1/2 and 572.5 MPa mm1/2 for inclusions in 20%, 30%, 40%, 60% and 100% region respectively. From the results presented in Fig. 12, it is observed that the percentage difference in SIF (with respect to inclusions in 100% region) is 1.4%, 1.0%, 0.5% and 0.02% for inclusions in 20%, 30%, 40% and 60% region respectively. Thus, the use of effective properties in the region away from the major crack not only improves the results but also reduces the CPU time.

The homogenization is not performed for minor cracks since the minor cracks modeling region does not have significant effect on the SIF. Apart from this, the presence of minor cracks does not affect the material properties. The results show that the use of effective properties in case of holes and inclusions improves the results. Hence, from the computational time point of view, the discontinuities are modeled only in the 20% region of the domain near the major crack without much compromising on the accuracy of the results. Since, discontinuities are modeled in 20% region near the major crack then the use of fine mesh in the entire domain is not a viable option. Therefore, to further reduce the computational time, only 20% region of the domain containing discontinuities is discretized by fine mesh and the rest of the domain is discretized by coarse mesh with effective properties. To maintain the displacement continuity at the junctions, special types of transition elements (explained in next section) are used.

10

S. Kumar et al. / Computers and Structures 150 (2015) 1–22

uy=0

1050 Minor cracks in 100% Region Minor cracks in 60% Region Minor cracks in 40% Region Minor cracks in 30% Region Minor cracks in 20% Region Without Minor cracks

1000 950

SIF (MPa.mm1/2)

900

ux=0

ux=0

uy=0

ux=C

RVE

uy=C

RVE

850 800

y

750

y x

700

uy=0

650

(a)

x

ux=0

(b)

Fig. 10. Schematic representation of boundary conditions to evaluate (a) Young’s modulus; (b) shear modulus.

600 550 500 18

20

22

24

26

28

30

32

34

36

1000

Crack length (mm)

Holes in 100% Region Holes in 60% Region Holes in 40% Region Holes in 30% Region Holes in 20% Region Without Holes

950 Fig. 8. SIF variation with crack length for an edge crack plate in the presence of multiple minor cracks.

900

SIF (MPa.mm1/2)

850 Table 3 Variation of SIF and CPU time with minor-cracks modeling region for a crack length of 20 mm. S. no.

Minor-cracks modeling region (%)

Stress intensity factor SIF (MPa mm1/2)

% Difference

1 2 3 4 5 6

Without minor-cracks 20 30 40 60 100

539.7 544.9 545.8 547.4 548.4 548.9

1.7 0.7 0.6 0.3 0.1 –

CPU time (s)

800 750 700 650

215 224 229 234 243 271

600 550 500 18

20

22

24

26

28

30

32

34

36

Crack length (mm) 5. Shape functions for transition elements

Fig. 11. SIF variation with crack length for an edge crack plate in the presence of multiple holes using effective properties.

In case of uniform mesh, only four node isoparametric quadrilateral elements (varies linearly along all the edges) are used for the simulation whereas, in case of non-uniform mesh, few transition elements are used as shown in Fig. 13. These elements possess 1, 2 or more mid-side nodes in addition to the standard corner nodes. The variation of displacement field along the edges in transition elements depends on the order. To ensure the compatibility, piecewise displacement interpolation functions are introduced for the transition elements, which are compatible with the adjacent regular four node quadrilateral elements. To obtain the shape functions for 6 node transition element, consider a quadrilateral transition element as shown in Fig. 14. The shape functions for the

Inclusions

Table 4 Variation of SIF and CPU time with holes modeling region for a crack length of 20 mm (effective properties). S. no.

1 2 3 4 5

Holes modeling region (%)

20 30 40 60 100

Voids

Equivalent Homogeneous Material with properties: Eeff , Geff and eff

Fig. 9. Schematic representation of homogenization process.

Stress intensity factor SIF (MPa mm1/2)

% Difference

597.4 602.0 607.5 611.8 612.7

2.5 1.7 0.8 0.1 –

S. Kumar et al. / Computers and Structures 150 (2015) 1–22

1050 Inclusions in 100% Region Inclusions in 60% Region Inclusions in 40% Region Inclusions in 30% Region Inclusions in 20% Region Without Inclusions

1000 950

SIF (MPa.mm1/2)

900

11

 27 2 N7 ðn; gÞ ¼ ð1  n2 Þð1 þ gÞ þn ; 32 6  27 2 2 ð1  n Þð1 þ gÞ n N8 ðn; gÞ ¼ 32 6 The derivatives of the shape functions in above equation are discontinuous at n = 0 and/or g = 0. Therefore, the shape functions for the corner nodes can be modified as,

850 800 750

N1 ðn; gÞ ¼

1 1 1 ð1  nÞð1  gÞ  fac5  N5 ðn; gÞ  fac6  N 6 ðn; gÞ 4

N2 ðn; gÞ ¼

1 2 2 ð1 þ nÞð1  gÞ  fac5  N5 ðn; gÞ  fac6  N 6 ðn; gÞ 4

N3 ðn; gÞ ¼

1 3 3 ð1 þ nÞð1 þ gÞ  fac7  N7 ðn; gÞ  fac8  N 8 ðn; gÞ 4

N4 ðn; gÞ ¼

1 4 4 ð1  nÞð1 þ gÞ  fac7  N7 ðn; gÞ  fac8  N 8 ðn; gÞ 4

700 650 600 550 500 18

20

22

24

26

28

30

32

34

36

Crack length (mm) Fig. 12. SIF variation with crack length for an edge crack plate in the presence of multiple inclusions using effective properties.

Table 5 Variation of SIF and CPU time with% of inclusions modeling region for initial crack length of 20 mm using effective properties. S. no.

1 2 3 4 5

Inclusions modeling region (%)

20 30 40 60 100

Stress intensity factor SIF

% Difference

564.7 567.1 569.9 572.4 572.5

1.4 1.0 0.5 0.02 –

Transition Element

j

where faci stands for factor of node i associated with node j and the 1 1 2 2 values are given as; fac5 ¼ 2=3, fac6 ¼ 1=3, fac5 ¼ 1=3, fac6 =2/3, 3 3 4 4 fac7 ¼ 2=3, fac8 ¼ 1=3, fac7 ¼ 1=3, fac8 ¼ 2=3. 6. Numerical simulations and discussion In this work, the fatigue life of an edge crack plate is simulated in the presence of multiple discontinuities (holes, inclusions and minor cracks). As explained in Section 4, 20% region of the plate is selected for modeling the discontinuities near the major crack to reduce the computational time without compromising on the accuracy. To further reduce the computational time, a local mesh refinement scheme is proposed in this work. In this approach, only 20% region of the domain containing the discontinuities is discretized by fine the mesh and the remaining 80% region of the domain is discretized by the coarse mesh. In the coarse mesh, effective properties are used for the problems of multiple holes and inclusions, while original material properties are used in the case of minor-cracks. In order to check the convergence, accuracy and effectiveness of the proposed approach, several fracture problems are simulated. 6.1. Finite size edge crack plate

Vary cubically along this edge

Vary linearly along this edge

Transition Element

Fig. 13. Schematic representation of transition element between coarse and fine mesh.

In this case, a rectangular plate with an initial edge crack of 20 mm is taken for the simulation. The loading, boundary conditions and material properties used for the simulations are taken same as per Section 3. A typical non-uniform discretization of the edge crack plate is illustrated in Fig. 15. The standard and transition elements with enriched nodes are also shown in Fig. 15. The convergence for stress intensity factor with number of elements is given in Table 6. From the table, it is observed that the results are converged for a non-uniform mesh of 4379 elements. Thus, 4379 elements mesh is used further to discretize the domain for the solution of the problems. To further validate the accuracy of the proposed XFEM, the SIF values under mode-I loading are computed analytically as given by Tada et al. [71],

KI ¼ r quadrilateral transition element in terms of element’s natural coordinates ðn; g 2 ½1; þ1Þ can be expressed as,  27 2 ð1  n2 Þð1  gÞ n ; N 5 ðn; gÞ ¼ 32 6  27 2 2 N 6 ðn; gÞ ¼ ð1  n Þð1  gÞ þn 32 6

pffiffiffiffiffiffi pa Fða=LÞ

where r is the uniform tensile load applied at the top edge of the plate and, F(a/L) = 1.122  0.231(a/L) + 10.55(a/L)2  21.71(a/ L)3 + 30.382(a/L)4. For an edge crack plate, the variations of SIF and fatigue life with crack length are shown in Figs. 16a and 16b respectively. The failure fatigue life obtained by analytical and proposed XFEM is found as 15886 cycles and 16,223 cycles corresponding to the crack length

12

S. Kumar et al. / Computers and Structures 150 (2015) 1–22

3

4

3

4

standard element

4

transition element

2

1

1

5

6

7

8

3

transition element

2

1

2

Fig. 14. Standard element and transition elements with different edge nodes.

Coarse Mesh Upper Region

Transition Elements

Fine Mesh

H

Middle Region

Tip Nodes

Coarse Mesh Lower Region

Split Nodes

L Fig. 15. A typical discretization of an edge crack plate with non-uniform mesh.

Table 6 Convergence study for an edge crack plate (20% region is discretized with fine mesh) for proposed XFEM. No of elements SIF (MPa mm1/2)

1881 537.6

3000 539.2

4379 539.6

1000 Analytical Solution Proposed XFEM

950

6018 539.7

900

of 33.4 mm and 33.5 mm respectively. The computation time required to simulate the edge crack problem, is found as 681 s and 67 s for standard and proposed XFEM respectively. It has been observed that the proposed XFEM results are in good agreement with the analytical results and is also computationally quite effective as compared to standard XFEM. In order to further check the effectiveness of the proposed approach, several problems of crack plate with holes, inclusions and minor cracks are simulated. In all these cases, the effect of discontinuities on the SIF is observed.

SIF (MPa.mm1/2)

850 800 750 700 650 600 550 500 18

6.2. Finite size edge crack interaction with other discontinuities In this section, four edge crack problems in the presence of hole, inclusion and minor crack have been solved by proposed XFEM and the results are compared with the available literature results. First, a finite size edge crack plate with an aligned minor crack is considered under mode-I loading. In the second problem, a finite size edge crack plate with two symmetric minor cracks is considered

20

22

24

26

28

30

32

34

Crack length (mm) Fig. 16a. SIF variation with crack length for an edge crack plate.

for the simulation. In third problem, again a finite size edge crack plate with a circular hole is solved by the proposed XFEM. Finally, an edge crack plate with a circular inclusion is taken for the simulation.

13

S. Kumar et al. / Computers and Structures 150 (2015) 1–22

34 Analytical Solution Proposed XFEM

32

Crack Length (mm)

30 28 26 24 22 20 18 0.2

0.4

0.6

0.8

1

1.2

1.4

1.6

1.8

x 104

No of cycles

Fig. 16b. Fatigue life variation with crack length for an edge crack plate.

6.2.2. Finite cize edge crack plate with two symmetric minor cracks In this case, an edge crack plate with two symmetric minor cracks as shown in Fig. 19 is taken for the simulation. The plate dimensions, minor crack length, macro crack length, loading and boundary conditions are kept similar as explained in article by Loehnert and Belytschko [31]. The mesh is also kept similar as used above. In this case, the shielding effect of minor cracks on the macro crack is observed. The normal stress contour plot in the vicinity of the crack tip is depicted in Fig. 20. From the figure, it is observed that the stresses due to minor cracks reduce the stresses in the vicinity of the macro crack tip; hence it provides the crack shielding effect. The value of mode-I SIF obtained by proposed approach is found as 0.699 whereas the SIF value in literature is found as 0.685 [31]. Thus, the result obtained by proposed XFEM shows a good agreement with the available literature result.

b

a0

a 1

l 2

H

3

6.2.1. Finite size edge crack plate with an aligned minor crack In this problem, a minor crack aligned with a macro crack as shown in Fig. 17 is taken for the simulation. The plate dimensions, macro and minor crack lengths, loading and boundary conditions are kept similar as explained in the article [30]. A non-uniform mesh of 37,247 elements which contains a very fine mesh in the discontinuous region and a coarse mesh in the continuous region is taken for the simulation. The numerical study performed by Loehnert and Belytschko [30] shows that the SIF for the three crack tips depend on the ratio of a/b as shown in Fig. 17. For the validation purpose, the length of minor crack is kept constant while the distance of the minor crack varies from the macro crack tip. In this case, the variation of mode-I SIF at the three crack tips is shown in Fig. 18 and the results are compared with the literature [30]. From these results, it has been observed that the minor crack causes crack amplification rather than shielding ahead of the macro crack tip and its effect decreases as the ratio of a/b increases. Further, the results obtained by proposed approach are found in good agreement with the literature results.

L Fig. 17. An edge crack plate with a minor crack of length l = b  a aligned with an edge crack.

6.2.3. Finite size edge crack plate with a hole In this case, a square plate with an edge crack (as explained in Section 6.2.1) and a circular hole of radius r0 = 0.1 is taken for the simulation as shown in Fig. 21. Initially, the center of hole is placed at x0 = 0.95 and y0 = 1.15. The loading and boundary conditions are also kept similar as explained in Section 6.2.1. The purpose of the numerical simulation is to analyze the effect of hole on the stress intensity factor value at the crack tip. The variation of mode-I SIF with x0/L is shown in Fig. 22. The results show that the hole just

1.4 Ref: (Loehnert and Belytschko, 2007) Proposed XFEM

1.2 1

at tip 1 at tip 2

b

KI

0.8

a0

0.6

H

at tip 3

0.4

H

0.2 0

b

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

2

1

a/b Fig. 18. The variation of mode-I stress intensity factor values with ratio of a/b at the macro and minor crack tips for constant minor crack length.

L Fig. 19. An edge crack plate with two symmetric minor cracks.

14

S. Kumar et al. / Computers and Structures 150 (2015) 1–22

5

center point: x0,y0

4

radius: r0 3

a0

H

2

1

0

L Fig. 20. Normal stress contour plot in the vicinity of the macro crack due to two symmetric minor cracks.

Fig. 23. An edge crack plate with circular inclusion.

1.25

center point: x0,y0

1.2

radius: r0

1.15

a0

1.1

KI

H

1.05 1 0.95 0.9 0.45

L

0.5

0.55

0.6

0.65

0.7

0.75

0.8

0.85

0.9

x0/L

Fig. 21. An edge crack plate with circular hole. Fig. 24. Variation of mode-I stress intensity factor at the crack tip due to interaction of inclusion.

1.3 1.25 1.2 1.15

KI

1.1 1.05 1 0.95 0.9 0.85 0.8 0.45

0.5

0.55

0.6

0.65

0.7

0.75

0.8

0.85

0.9

x0/L Fig. 22. Variation of mode-I stress intensity factor at the crack tip due to interaction of hole.

behind the crack tip provides the shielding effect and ahead of the crack tip intensify the stress intensity factor. Further, the effect of hole on the SIF decreases, as we move away from the crack tip, and finally it is negligible.

6.2.4. Finite size edge crack plate with a inclusion In this case, a square plate with an edge crack (as explained above) and a circular inclusion of radius r0 = 0.1 mm is taken for the simulation as shown in Fig. 23. Initially, the center of inclusion is placed at x0 = 0.95 mm and y0 = 1.15 mm. The loading and boundary conditions are kept similar as explained in Section 6.2.1. The numerical simulations are performed to examine the effect of inclusion on the SIF. The variation of mode-I SIF with x0/L is shown in Fig. 24. The results show that the inclusion just behind the crack tip reduces the SIF whereas it amplifies the SIF ahead of the crack tip. The effect of inclusion on SIF decreases as it moves away from the macro crack tip and finally becomes negligible. To further check the accuracy and computational efficiency of the proposed XFEM, the plate with multiple holes, inclusions and minor cracks (as explained in Section 4) lying in 20% region of the domain are simulated by XFEM. The SIF and CPU time obtained by proposed and standard XFEM using effective properties for a crack length of 20 mm are given in Table 7. From the results presented in Table 7, it is seen that the standard XFEM with uniform mesh and proposed XFEM with non-uniform mesh provides almost same values of SIF for all the cases considered. From Table 7, it is also observed that the main advantage of the proposed XFEM is that it reduces the computation time considerably without

15

S. Kumar et al. / Computers and Structures 150 (2015) 1–22 Table 7 SIF and CPU time for holes, inclusions and minor-cracks lying in 20% region with a crack length of 20 mm. Type of discontinuity

Holes Inclusions Minor-cracks Without discontinuities

SIF (MPa mm1/2)

CPU Time (s)

Standard XFEM

Proposed XFEM

% Difference

Standard XFEM

Proposed XFEM

597.4 564.7 544.9 539.7

595.6 563.9 545.9 539.6

0.3 0.1 0.2 0.3

274 422 224 215

34 87 15 14

1000 950

Standard XFEM (Holes in 100% region) Proposed XFEM (Holes in 20% region with eff. prop)

900

SIF (MPa.mm1/2)

Effective Properties Eeff , Geff and eff

850 800 750 700 650

H

600 550 500 18

20

22

24

26

28

30

32

Crack length (mm)

Effective Properties Eeff , Geff and

eff

Fig. 26a. SIF variation with crack length for an edge crack plate in the presence of holes.

32 30

Fig. 25. Edge crack plate along with multiple holes.

compromising on the accuracy. Therefore, for the further analysis, the discontinuities are modeled in 20% region by a non-uniform mesh. 6.3. Finite size edge crack plate with multiple holes In this case, an edge crack plate of size 100 mm  200 mm having multiple holes in the central 20% region of the total domain is taken for the simulation as shown in Fig. 25. The modeling conditions i.e. loading, boundary conditions, mesh size and crack increment are kept similar as explained in Section 6.1. In actual problem, 100 holes of 10% volume fraction are distributed in the entire domain but as per proposed approach, only 20 holes (10% by volume) are modeled in the central 20% region of the plate. The remaining 80% region (lower and upper region) of the plate is given effective properties. The effective Young’s modulus obtained by homogenization approach is found as 52.4 GPa. The radius of the holes is varying from 2.0 mm to 3.5 mm. The results are compiled ten times to obtain the average effect of randomly distributed holes on the major crack. To check the accuracy of the proposed XFEM, edge crack plate with randomly distributed 100 holes (10% volume fraction) is also solved by standard XFEM. In standard XFEM, a uniform mesh size of 88 nodes  176 nodes is used for the simulations. The average stress intensity factor and fatigue life with crack length are shown in Figs. 26a and 26b respectively. For this case, the final fatigue life obtained by proposed and standard XFEM is found as 10,575 cycles and 10,014

Crack Length (mm)

L

Standard XFEM (Holes in 100% region) Proposed XFEM (Holes in 20% region with eff. prop)

28 26 24 22 20 18 2000 3000 4000 5000 6000 7000 8000 9000 10000 11000

No of cycles Fig. 26b. Fatigue life variation with crack length for an edge crack plate in the presence of holes.

cycles with a corresponding failure crack length of 30.7 mm and 30.4 mm respectively. These results also show that the fatigue life of the plate is reduced by 32.3% due to the presence of 10% holes. 6.4. Finite size edge crack plate with multiple inclusions In this case, a rectangular plate of size 100 mm  200 mm along with an edge crack of initial length, 20 mm is taken for simulation as shown in Fig. 27. The multiple inclusions (10% by volume) are randomly distributed in 20% region of the domain. The loading, boundary conditions, mesh size and crack increments are kept

16

S. Kumar et al. / Computers and Structures 150 (2015) 1–22

34 32

Standard XFEM (Inclusions in 100% region) Proposed XFEM (Inclusions in 20% region with eff. prop)

Crack Length (mm)

30

Effective Properties Eeff , Geff and eff

28 26 24 22

H

20 18 2000

4000

6000

8000

10000

12000

14000

No of cycles

Effective Properties Eeff , Geff and

Fig. 28b. Fatigue life variation with crack length for an edge crack plate in the presence of inclusions.

eff

L Fig. 27. Edge crack plate along with multiple inclusions.

similar as above. In actual problem, 100 inclusions (10% by volume) are randomly distributed in the entire domain but from the computational time point of view, only 20 inclusions (10% by volume) in central 20% region of the total domain are modeled. The remaining 80% region (lower and upper region) of the plate is provided effective properties. The effective Young’s modulus obtained by homogenization approach is found as 60.8 GPa. The radii of the inclusions vary from 2.0 mm to 3.5 mm. The results are compiled ten times to get the average effect of randomly distributed inclusions on the major crack. To verify the accuracy of the proposed approach, the problem is also solved by standard XFEM with 10% inclusion in the entire plate. A uniform nodal distribution of 88  176 is used in standard XFEM simulations. The variation of average stress intensity factor and fatigue life with crack length

1000 950

L

Standard XFEM (Inclusions in 100% region) Proposed XFEM (Inclusions in 20% region with eff. prop)

Fig. 29. Edge crack plate along with multiple minor cracks.

900

SIF (MPa.mm1/2)

H

850

is shown in Figs. 28a and 28b respectively. The final fatigue life evaluated by proposed and standard XFEM is found as 13,544 cycles and 12,966 cycles with corresponding failure crack length of 32.2 mm and 32.0 mm respectively. From these results, it is observed that the fatigue life of the plate is reduced by 16.5% due to the presence of 10% inclusions.

800 750 700 650 600

6.5. Finite size edge crack plate with multiple minor cracks

550 500 18

20

22

24

26

28

30

32

34

Crack length (mm) Fig. 28a. SIF variation with crack length for an edge crack plate in the presence of inclusions.

In this case, a rectangular plate of size 100 mm  200 mm with a major edge crack of initial length a = 20 mm is taken for simulation as shown in Fig. 29. The minor cracks are randomly distributed in 20% region of the total domain. The loading, boundary conditions, mesh size and crack increments are kept similar as explained

17

S. Kumar et al. / Computers and Structures 150 (2015) 1–22

1000 Standard XFEM (Minor cracks in 100% region) Proposed XFEM (Minor cracks in 20% region)

950

SIF (MPa.mm1/2)

900 850

Effective Properties Eeff , Geff and eff

800 750 700 650

H

600 550 500 18

20

22

24

26

28

30

32

34

Crack length (mm)

Effective Properties Eeff , Geff and eff

Fig. 30a. SIF variation with crack length for an edge crack plate in the presence of minor cracks.

34 32

Standard XFEM (Minor cracks in 100% region) Proposed XFEM (Minor cracks in 20% region)

L Fig. 31. Edge crack plate along with multiple holes and inclusions.

28 26 24 22 20 18 2000

4000

6000

8000

10000

12000

14000

16000

No of cycles Fig. 30b. Fatigue life variation with crack length for an edge crack plate in the presence of minor cracks.

above. In actual problem, 80 minor cracks are distributed in the entire domain but in the present case, only 20% region of the total domain is modeled with 16 minor cracks of size varying from 3.0 mm to 4.5 mm. The results are compiled ten times to obtain the average values of SIF. To further check the accuracy and effectiveness of the proposed XFEM, plate with 80 randomly distributed minor cracks in the entire domain is simulated by standard XFEM using a uniform mesh size of 88  176 nodes. The variation of SIF with crack length is plotted in Fig. 30a. A plot of fatigue life with crack length is shown in Fig. 30b. The fatigue life obtained by proposed and standard XFEM for this case is found as 15,573 cycles and 15,077 cycles with corresponding failure crack length of 33.0 mm 32.8 mm respectively. Thus, the life of cracked plate is reduced by 4.0% due to the presence of 10% minor cracks.

holes and inclusions are distributed in 20% region of the total domain with effective properties in remaining 80% of the domain. The effective Young’s modulus obtained by homogenization approach is found as 56.5 GPa. The radius of both holes and inclusions vary from 2.0 mm to 3.5 mm. The loading, boundary conditions, mesh size and crack increments are kept same as above. The results are compiled ten times to obtain the average effect of randomly distributed discontinuities (holes and inclusions). To ensure the accuracy of the proposed XFEM, the same problem with randomly distributed discontinuities in the entire domain is also solved by standard XFEM using a uniform mesh size of 88  176 nodes. The variation of SIF and fatigue life with crack length is plotted in Figs. 32a and 32b respectively. For this case, the failure fatigue life obtained by proposed XFEM and standard

1000 950

Standard XFEM (Holes & inclusions in 100% region) Proposed XFEM (Holes & inclusions in 20% region with eff. prop)

900

SIF (MPa.mm1/2)

Crack Length (mm)

30

850 800 750 700 650 600

6.6. Finite size edge crack plate with multiple holes and inclusions In this case, fatigue crack growth is modeled in the presence of multiple holes and inclusions. A rectangular plate of size 100 mm  200 mm with an edge crack of initial length a = 20 mm containing 10% discontinuities i.e. 10 holes and 10 inclusions are taken for the simulation as shown in Fig. 31. These

550 500 18

20

22

24

26

28

30

32

Crack length (mm) Fig. 32a. SIF variation with crack length for an edge crack plate in the presence of holes and inclusions.

18

S. Kumar et al. / Computers and Structures 150 (2015) 1–22

32

1000 950

Standard XFEM (Holes & Mics in 100% region) Proposed XFEM (Holes & Mics in 20% region with eff. prop)

30

SIF (MPa.mm1/2)

Crack Length (mm)

900 28 26 24

850 800 750 700 650

22

600 20

Standard XFEM (Holes & inclusions in 100% region) Proposed XFEM (Holes & inclusions in 20% region with eff. prop)

550

18 3000 4000 5000 6000 7000 8000 9000 10000 11000 12000 13000

500 18

No of cycles

20

22

24

26

28

30

32

Crack length (mm)

Fig. 32b. Fatigue life variation with crack length for an edge crack plate in the presence of holes and inclusions.

Fig. 34a. SIF variation with crack length for an edge crack plate along with holes and minor cracks.

32

Standard XFEM (Holes & Mics in 100% region) Proposed XFEM (Holes & Mics in 20% region with eff. prop)

Crack Length (mm)

30

Effective Properties Eeff , Geff and eff

H

28 26 24 22 20 18 2000 3000 4000 5000 6000 7000 8000 9000 10000 11000 12000

No of cycles

Effective Properties Eeff , Geff and eff

L Fig. 33. Edge crack plate along with multiple holes and minor cracks.

XFEM is found as 12,271 cycles and 11,736 cycles with corresponding crack length of 31.4 mm and 31.2 mm respectively. In this case, the fatigue life of the plate is reduced by 24.4% due to the presence of 10% holes and inclusions.

Fig. 34b. Fatigue life variation with crack length for an edge crack plate along with holes and minor cracks.

same as explained earlier. The results are compiled ten times to obtain an average effect of randomly distributed holes and minor cracks. To verify the accuracy of the proposed XFEM, the same problem with randomly distributed holes and minor cracks in the entire domain is also solved by standard XFEM using a uniform mesh size of 88  176 nodes. The variation of SIF and fatigue life with crack length is plotted in Figs. 34a and 34b respectively. From these results, it can be noticed that the final fatigue life obtained by proposed and standard XFEM is found as 10,281 cycles and 9809 cycles corresponding to 30.4 mm and 30.3 mm crack length respectively. The life of the plate is reduced by 36.6% due to the presence of 10% holes and minor cracks.

6.7. Finite size edge crack plate with multiple holes and minor cracks In this case, a major crack of length 20 mm is taken at the edge of the plate (100 mm  200 mm) as shown in Fig. 33. Twenty holes/voids (10% by volume) and sixteen minor cracks are randomly distributed in 20% region of the total domain. The holes radii vary from 2.0 mm to 3.5 mm whereas the size of minor cracks varies from 3.0 mm to 4.5 mm. The remaining 80% region (lower and upper region) of the plate is given effective properties. The loading, boundary conditions, mesh size and crack increments are taken

6.8. Finite size edge crack plate with multiple inclusions and minor cracks In this case, a rectangular plate of size 100 mm  200 mm with an edge crack of length 20 mm is taken for simulation in the presence of randomly distributed 20 inclusions and 16 minor cracks as shown in Fig. 35. The inclusions radii vary from 2.0 mm to 3.5 mm and the size of minor cracks varies from 3.0 mm to 4.5 mm. Only 20% central region of the plate is modeled with discontinuities,

19

S. Kumar et al. / Computers and Structures 150 (2015) 1–22

32

Crack Length (mm)

30

Effective Properties Eeff , Geff and eff

H

28 26 24 22 20 Standard XFEM (Inclusions & Mics in 100% region) Proposed XFEM (Inclusions & Mics in 20% region with eff. prop)

18 2000

4000

6000

8000

10000

12000

14000

No of cycles

Effective Properties Eeff , Geff and eff

Fig. 36b. Fatigue life variation with crack length for an edge crack plate in the presence of inclusions and minor cracks.

6.9. Finite size edge crack plate with multiple holes, inclusions and minor cracks

L Fig. 35. Edge crack plate along with multiple inclusions and minor cracks.

1000 Standard XFEM (Inclusions & Mics in 100% region)

950

Proposed XFEM (Inclusions & Mics in 20% region with eff. prop)

900

SIF (MPa.mm1/2)

850 800

In the final case, the fatigue life of the edge crack plate (100 mm  200 mm) is obtained in the presence of holes, inclusions and minor cracks. In this study, 10% discontinuities (10 holes, 10 inclusions and 16 minor cracks) are randomly distributed in the 20% region of the total domain as shown in Fig. 37. The holes and inclusions radii vary from 2.0 mm to 3.5 mm while the size of minor cracks varies from 3.0 mm to 4.5 mm. The modeling conditions i.e. loading, boundary conditions, mesh size and crack increments are kept same as above. The results are complied ten times to obtain the average effect of randomly distributed discontinuities on the major crack. To ensure the accuracy of the proposed

750 700 650 600 550 500 18

20

22

24

26

28

30

32

Effective Properties Eeff , Geff and eff

Crack length (mm) Fig. 36a. SIF variation with crack length for an edge crack plate in the presence of inclusions and minor cracks.

H while the remaining 80% continuum region (lower and upper region) of the plate is given effective properties. The programme is compiled ten times to get the average effect of randomly distributed inclusions and minor cracks on the major crack. To check the accuracy of the proposed XFEM, the problem with inclusions and minor cracks in the entire domain is solved by standard XFEM. A uniform nodal distribution of 88  176 nodes is used for standard XFEM. In this case, the stress intensity factor and fatigue life with crack length is plotted in Figs. 36a and 36b respectively. The fatigue life obtained by proposed and standard XFEM is found as 12,759 cycles and 12,203 cycles corresponding to 31.4 mm and 31.2 mm crack lengths respectively. From these results, it is observed that the fatigue life of the plate is reduced by 21.3% due to the presence of 10% inclusions and minor cracks.

Effective Properties Eeff , Geff and eff

L Fig. 37. Edge crack plate along with multiple holes, inclusions and minor cracks.

20

S. Kumar et al. / Computers and Structures 150 (2015) 1–22

SIF with crack length. A plot of the fatigue life with crack length is shown in Fig. 38b. These simulations show that the fatigue life obtained by proposed and standard XFEM is found as 11,475 cycles and 10,935 cycles corresponding to a crack length of 30.7 mm and 30.4 mm respectively. Hence, the fatigue life of the plate is reduced by 29.3% due to the presence of 10% discontinuities i.e. holes, inclusions and minor cracks. To examine the computational effectiveness of the proposed homogenized XFEM, all above problems are also solved by standard XFEM with uniform mesh size of 88  176 nodes. In standard XFEM, the multiple discontinuities (holes, inclusions and minor cracks) are considered in the entire domain. The CPU time required to simulate one step in case of proposed and standard XFEM is tabulated in Table 8 for all the problems. From the results, it has been observed that the proposed XFEM is computationally accurate and quite efficient as compared to standard XFEM.

950 900

SIF (MPa.mm1/2)

850 800 750 700 650 600 550

Standard XFEM (Holes, Incs & Mics in 100% region) Proposed XFEM (Holes, Incs & Mics in 20% region with eff. prop)

500 18

20

22

24

26

28

30

32

Crack length (mm) 7. Conclusions Fig. 38a. SIF variation with crack length for an edge crack plate in the presence of holes, inclusions and minor cracks.

32

Crack Length (mm)

30 28 26 24

In the present study, a homogenized XFEM has been proposed for the evaluation of fatigue life of an edge crack plate in the presence of multiple discontinuities i.e. holes, inclusions and minor cracks. A detailed numerical study is performed to decide the modeling region for the discontinuities. To improve the results, a homogenization scheme based on strain energy density approach is used to determine the effective properties of heterogeneous region containing the holes, inclusions and both. A non-uniform meshing is used to discretize the entire domain. The central 20% region with multiple discontinuities is discretized by very fine mesh while the remaining 80% region is discretized by the coarse mesh. The transition elements are used to maintain the displacement continuity at the junction nodes. On the basis of present simulations, the following conclusions can be drawn,

22 20 Standard XFEM (Holes, Incs & Mics in 100% region) Proposed XFEM (Holes, Incs & Mics in 20% region with eff. prop)

18 2000 3000 4000 5000 6000 7000 8000 9000 10000 11000 12000

No of cycles Fig. 38b. Fatigue life variation with crack length for an edge crack plate in the presence of holes, inclusions and minor cracks.

Table 8 CPU time for standard and proposed XFEM in case of multiple defects for a crack length of 20 mm. Case

CPU time (s) Standard XFEM

Plate with an edge Plate with an edge Plate with an edge inclusions Plate with an edge Plate with an edge Plate with an edge cracks Plate with an edge cracks Plate with an edge minor cracks

crack crack and multiple holes crack and multiple

Proposed XFEM

215 491 1450

14 34 87

crack and minor cracks crack, holes and inclusions crack, holes and minor

271 588 553

15 57 41

crack, inclusions and minor

1565

104

crack, holes, inclusions and

621

62

method, the same problem with multiple discontinuities in the entire domain is solved by standard XFEM using a uniform nodal distribution of 88  176 nodes. Fig. 38a shows the variation of

 A strain energy based homogenization scheme is used to evaluate the effective properties of the heterogeneous plate.  Special transition elements are developed to maintain the displacement continuity at the junction nodes.  Only 20% central region of the plate is found suitable for modeling the discontinuities.  The proposed homogenized XFEM approach drastically reduces the CPU time without compromising on the accuracy.  The presence of discontinuities significantly affects the fatigue life.  The effect of holes is found more severe as compared to the minor crack and inclusions.  10% holes (by volume) reduce the fatigue life of the plate nearly by 32%.  The proposed approach shows a good convergence, and saves huge computational time.

References [1] Agarwal A, Singh IV, Mishra BK. Evaluation of elastic properties of interpenetrating phase composites by mesh-free method. J Compos Mater 2012;47:1407–23. [2] Aerospace Specification Metals (ASM) material data sheet, Aluminum alloy 6061-T6. . [3] Belytschko T, Gu L, Lu YY. Fracture and crack growth by element-free Galerkin methods. Model Simul Mater Sci Eng 1994;2:519–34. [4] Belytschko T, Lu YY, Gu L. Crack propagation by element-free Galerkin methods. Eng Fract Mech 1995;51:295–315. [5] Belytschko T, Black T. Elastic crack growth in finite elements with minimal remeshing. Int J Numer Meth Eng 1999;45:601–20. [6] Bilger N, Auslender F, Bornert M, Michel JC, Moulinec H, Suquet P, et al. Effect of a nonuniform distribution of voids on the plastic response of voided

S. Kumar et al. / Computers and Structures 150 (2015) 1–22

[7]

[8] [9]

[10]

[11]

[12]

[13] [14] [15] [16]

[17] [18] [19] [20]

[21] [22]

[23] [24] [25]

[26]

[27] [28]

[29] [30] [31] [32]

[33] [34]

[35] [36] [37]

[38] [39] [40]

materials: a computational and statistical analysis. Int J Solids Struct 2005;42:517–38. Bordas S, Moran B. Enriched finite elements and level sets for damage tolerance assessment of complex structures. Eng Fract Mech 2006;73:1176–201. Bordas SP, Nguyen PV, Dunant C, Guidoum A, Dang HN. An extended finite element library. Int J Numer Meth Eng 2007;71:703–32. Bordas S, Rabczuk T, Zi G. Three-dimensional crack initiation, propagation, branching and junction in non-linear materials by an extended meshfree method without asymptotic enrichment. Eng Fract Mech 2008;75:943–60. Bordas SPA, Natarajan S, Kerfriden P, Augarde CE, Mahapatra DR, Rabczuk T, et al. On the performance of strain smoothing for quadratic and enriched finite element approximations (XFEM/GFEM/PUFEM). Int J Numer Meth Eng 2011;86:637–66. Bouhala L, Shao Q, Koutsawa Y, Younes A, Nunez P, Makradi A, et al. An XFEM crack-tip enrichment for a crack terminating at a bi-material interface. Eng Fract Mech 2013;102:51–64. Chen JS, Wang HP. Some recent improvements in meshfree methods for incompressible finite elasticity boundary value problems with contact. Comput Mech 2000;25:137–56. Cheung S, Luxmoore AR. A finite element analysis of stable crack growth in an aluminum alloy. Eng Fract Mech 2003;70:1153–69. Choi CK, Lee NH. Three dimensional transition solid elements for adaptive mesh gradation. Struct Eng Mech 1993;1:61–74. Christensen RM, Lo KH. Solutions for effective shear properties in three phase sphere and cylinder models. J Mech Phys Sol 1979;27:315–30. Daux C, Moes N, Dolbow J, Sukumar N, Belytschko T. Arbitrary branched and intersecting cracks with the extended finite element method. Int J Numer Meth Eng 2000;48:1741–60. Duflot M, Dang HN. Fatigue crack growth analysis by an enriched meshless method. J Comput Appl Math 2004;168:155–64. Farhat C. A Lagrange multiplier based divide and conquer finite element algorithm. J Comput Syst Eng 1991;2:149–56. Fries TP. A corrected XFEM approximation without problems in blending elements. Int J Numer Meth Eng 2008;75:503–32. Guidault PA, Allix O, Champaney L, Cornuault C. A multiscale extended finite element method for crack propagation. Comput Meth Appl Mech Eng 2008;197:381–99. Gupta AK. A finite element for transition from a fine to a coarse grid. Int J Numer Meth Eng 1978;12:35–45. Gravouil A, Moes N, Belytschko T. Non-planar 3D crack growth by the extended finite element and level sets – Part II: Level set update. Int J Numer Meth Eng 2002;53:2569–86. Hill R. A self-consistent mechanics of composite materials. J Mech Phys Sol 1965;13:213–22. Holl M, Loehnert S, Wriggers P. An adaptive multiscale method for crack propagation and crack coalescence. Int J Numer Meth Eng 2013;93:23–51. Hughes TJR. Multiscale phenomena: Green’s functions, the Dirichlet-toNeumann formulation, subgrid scale models, bubbles and the origins of stabilized methods. Comput Methods Appl Mech Eng 1995;127:387–401. Kim DJ, Pereira JP, Duarte CA. Analysis of three-dimensional fracture mechanics problems: a two-scale approach using coarse-generalized FEM meshes. Int J Numer Meth Eng 2010;81:335–65. Kumar S, Singh IV, Mishra BK. XFEM simulation of stable crack growth using J– R curve under finite strain plasticity. Int J Mech Mater Des 2014;10:165–77. Kumar S, Singh IV, Mishra BK. A multigrid coupled (FE-EFG) approach to simulate fatigue crack growth in heterogeneous materials. Theoret Appl Fract Mech 2014. Kushch VI. Interacting cracks and inclusions in a solid by multipole expansion method. Int J Solids Struct 1998;35:1751–62. Loehnert S, Belytschko T. Crack shielding and amplification due to multiple microcracks interacting with a macrocrack. Int J Fract 2007;145:1–8. Loehnert S, Belytschko T. A multiscale projection method for macro/ microcrack simulations. Int J Numer Meth Eng 2007;71:1466–82. Loehnert S, Mueller-Hoeppe DS, Wriggers P. 3D corrected XFEM approach and extension to finite deformation theory. Int J Numer Meth Eng 2011;86:431–52. Melenk JM, Babuska I. The partition of unity finite element method: basic theory and applications. Comput Methods Appl Mech Eng 1996;139:289–314. McDill JM, Goldak JA, Oddy AS, Bibby MJ. Isoparametric quadrilaterals and hexahedrons for mesh-grading-algorithms. Commun Appl Numer Meth 1987;3:155–63. McLaughlin R. A study of the differential scheme for composite materials. Int J Eng Sci 1977;15:237–44. Menk A, Bordas SPA. A robust preconditioning technique for the extended finite element method. Int J Numer Meth Eng 2011;85:1609–32. Michel JC, Moulinec H, Suquet R. Effective properties of composite materials with periodic microstructure: a computational approach. Comput Methods Appl Mech Eng 1999;172:109–43. Moes N, Dolbow J, Belytschko T. A finite element method for crack growth without remeshing. Int J Numer Meth Eng 1999;46:131–50. Mori T, Tanaka K. Average stress in the matrix and average elastic energy of materials with misfitting inclusions. Acta Metall 1973;21:571–4. Morton DJ, Tyler JM, Dorroh JR. New 3D finite element for adaptive hrefinement in 1-irregular meshes. Int J Numer Meth Eng 1995;38:3989–4008.

21

[41] Muller WH. Mathematical versus experimental stress analysis of inhomogeneities in solids. J Phys IV 1996;6. pp. C1-139–C1-148. [42] Natarajan S, Bordas S, Mahaptra DR. Numerical integration over arbitrary polygonal domains based on Schwarz-Christoffel conformal mapping. Int J Numer Meth Eng 2009;80:103–34. [43] Natarajan S, Mahaptra DR, Bordas S. Integrating strong and weak discontinuities without integration subcells and example applications in an XFEM/GFEM framework. Int J Numer Meth Eng 2010;83:269–94. [44] Natarajan S, Song C. Representation of singular fields without asymptotic enrichment in the extended finite element method. Int J Numer Meth Eng 2013;96:813–41. [45] Nguyen VP, Rabczuk T, Bordas S, Duflot M. Meshless methods: a review and computer implementation aspects. Math Comput Simul 2008;79:763–813. [46] Oden JT, Zohdi TI. Analysis and adaptive modeling of highly heterogeneous elastic structures. Comput Methods Appl Mech Eng 1997;148:367–91. [47] Osher S, Sethian JA. Fronts propagating with curvature-dependent speed: algorithms based on Hamilton-Jacobi formulations. J Comput Phys 1988;79:12–49. [48] Palle N, Dantzig JA. An adaptive mesh refinement scheme for solidification problems. Metall Mater Trans A 1996;27A:707–18. [49] Paris PC, Gomez MP, Anderson WE. A rational analytic theory of fatigue. Trend Eng 1961;13:9–14. [50] Passieux JC, Gravouil A, Rethore J, Baietto MC. Direct estimation of generalized stress intensity factors using a three-scale concurrent multigrid X-FEM. Int J Numer Meth Eng 2011;85:1648–66. [51] Passieux JC, Rethore J, Gravouil A, Baietto MC. Local/global non-intrusive crack propagation simulation using a multigrid X-FEM solver. Comput Mech 2011;52:1381–93. [52] Pathak H, Singh A, Singh IV. Fatigue crack growth simulations of bi-material interfacial cracks under thermo-elastic loading by extended finite element method. Euro J Comput Mech 2013;22:79–104. [53] Pereira JP, Kim DJ, Duarte CA. A two-scale approach for the analysis of propagating three-dimensional fractures. Comput Mech 2012;49:99–121. [54] Pierres E, Baietto MC, Gravouil A. A two-scale extended finite element method for modelling 3D crack growth with interfacial contact. Comput Methods Appl Mech Eng 2010;199:1165–77. [55] Poniznik Z, Salit V, Basista M, Gross D. Effective elastic properties of interpenetrating phase composites. Comput Mater Sci 2008;44:813–20. [56] Provatas N, Goldenfeld N, Dantzig J. Adaptive mesh refinement computation of solidification microstructures using dynamic data structures. J Comput Phys 1999;148:265–90. [57] Rabczuk T, Belytschko T. Adaptivity for structured meshfree particle methods in 2D and 3D. Int J Numer Meth Eng 2005;63:1559–82. [58] Rabczuk T, Belytschko T. A three-dimensional large deformation meshfree method for arbitrary evolving cracks. Comput Methods Appl Mech Eng 2007;196:2777–99. [59] Rabczuk T, Bordas S, Zi G. A three-dimensional meshfree method for continuous multiple-crack initiation, propagation and junction in statics and dynamics. Comput Mech 2007;40:473–95. [60] Rannou J, Gravouil A, Baietto-Dubourg MC. A local multigrid X-FEM strategy for 3-D crack propagation. Int J Numer Meth Eng 2009;77:581–600. [61] Rabczuk T, Bordas S, Zi G. On three-dimensional modelling of crack growth using partition of unity methods. Comput Struct 2010;88:1391–411. [62] Rao BN, Rahman S. An interaction integral method for analysis of cracks in orthotropic functionally graded materials. Comput Mech 2003;32:40–51. [63] Samet H. Application of spatial data structure. New York: Addison-Wesley; 1990. [64] Sanchez-Palencia E, Zaoui A. Homogenization Techniques for composite media, lecture notes in physics, vol. 272. Berlin: Springer; 1987. [65] Singh IV, Mishra BK, Bhattacharya S, Patil RU. The numerical simulation of fatigue crack growth using extended finite element method. Int J Fatigue 2012;36:109–19. [66] Stolarska M, Chopp DL, Moes N, Belytschko T. Modeling crack growth by level sets in the extended finite element method. Int J Numer Meth Eng 2001;51:943–60. [67] Sukumar N, Chopp DL, Moes N, Belytschko T. Modeling of holes and inclusions by level sets in the extended finite element method. Comput Methods Appl Mech Eng 2001;190:6183–200. [68] Sukumar N, Prevost JH. Modeling quasi-static crack growth with extended finite element method Part I: Computer implementation. Int J Solids Struct 2003;40:7513–37. [69] Sze KY, Wu D. Transition finite element families for adaptive analysis of axisymmetric elasticity problems. Finite Elem Anal Des 2011;47: 360–72. [70] Tabarraei A, Sukumar N. Adaptive computations on conforming quadtree meshes. Finite Elem Anal Des 2005;41:686–702. [71] Tada H, Paris PC, Irwin GR. The stress analysis of cracks handbook. 3rd ed. ASME Press; 2000. [72] Terada K, Kikuchi N. A class of general algorithms for multi-scale analyses of heterogeneous media. Comput Methods Appl Mech Eng 2001;190:5427–64. [73] Torquato S. Modeling of physical properties of composite materials. Int J Solids Struct 2000;37:411–22. [74] Ventura G. On the elimination of quadrature subcells for discontinuous functions in the extended finite-element method. Int J Numer Meth Eng 2006;66:767–95.

22

S. Kumar et al. / Computers and Structures 150 (2015) 1–22

[75] Wyart E, Duflot M, Coulon D, Martiny P, Pardoen T, Remacle JF, et al. Substructuring FE-XFE approaches applied to three-dimensional crack propagation. J Comput Appl Math 2008;215:626–38. [76] Wu D, Sze KY, Lo SH. Two- and three-dimensional transition element families for adaptive refinement analysis of elasticity problems. Int J Numer Meth Eng 2009;78:587–630. [77] Yan AM, Dang N. Multiple-cracked fatigue crack growth by BEM. Comput Mech 1995;16:273–80.

[78] Yan X. A boundary element modeling of fatigue crack growth in a plane elastic plate. Mech Res Commun 2006;33:470–81. [79] Zi G, Belytschko T. New crack-tip elements for XFEM and applications to cohesive cracks. Int J Numer Meth Eng 2003;57:2221–40. [80] Zohdi TI, Oden JT, Rodin GJ. Hierarchical modeling of heterogeneous bodies. Comput Methods Appl Mech Eng 1996;138:273–98.