Dynamic analysis of interfacial crack problems in anisotropic bi-materials by a time-domain BEM

Engineering Fracture Mechanics 76 (2009) 1996–2010

Contents lists available at ScienceDirect

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

Dynamic analysis of interfacial crack problems in anisotropic bi-materials by a time-domain BEM Jun Lei a, Felipe Garcia-Sanchez b, Michael Wünsche c, Chuanzeng Zhang c,*, Yue-Sheng Wang d, Andres Saez e a

Department of Engineering Mechanics, Beijing University of Technology, Beijing 100124, PR China Departamento de Ingenieria Civil, de Materiales y Fabricacion, Universidad de Malaga, 29013 Malaga, Spain c Department of Civil Engineering, University of Siegen, D-57068 Siegen, Germany d Institute of Engineering Mechanics, School of Civil Engineering, Beijing Jiaotong University, Beijing 100044, PR China e Departamento de Mecanica de Medios Continuos, Escuela Superior de Ingenieros, Universidad de Sevilla, 41092 Sevilla, Spain b

a r t i c l e

i n f o

Article history: Received 2 October 2008 Received in revised form 30 April 2009 Accepted 9 May 2009 Available online 18 May 2009 Keywords: Dynamic fracture mechanics Dynamic stress intensity factors Anisotropic bi-materials Interface cracks Time-domain boundary element method

a b s t r a c t A time-domain boundary element method (BEM) together with the sub-domain technique is applied to study dynamic interfacial crack problems in two-dimensional (2D), piecewise homogeneous, anisotropic and linear elastic bi-materials. The bi-material system is divided into two homogeneous sub-domains along the interface and the traditional displacement boundary integral equations (BIEs) are applied on the boundary of each sub-domain. The present time-domain BEM uses a quadrature formula for the temporal discretization to approximate the convolution integrals and a collocation method for the spatial discretization. Quadratic quarter-point elements are implemented at the tips of the interface cracks. A displacement extrapolation technique is used to determine the complex dynamic stress intensity factors (SIFs). Numerical examples for computing the complex dynamic SIFs are presented and discussed to demonstrate the accuracy and the efficiency of the present time-domain BEM. Ó 2009 Elsevier Ltd. All rights reserved.

1. Introduction Dissimilar or layered materials such as ceramic–metal, polymer-composite and composite-metal have been widely used in engineering structures for the purpose of strengthening or lightening. In these materials, interfacial debonding or cracking often becomes the principal failure mode in the process of the manufacture or the in-service applications. Additionally, composite materials are either intrinsically anisotropic or anisotropic on some length scale of the deformation and frequently subjected to dynamic loading conditions. As a result, dynamic interfacial failure in dissimilar anisotropic materials is an important and promising research subject in material science, fracture and damage mechanics and non-destructive material evaluation. In contrast to cracks in homogeneous, isotropic/anisotropic and linear elastic materials, the local asymptotic displacement and stress fields near the tips of interfacial cracks in dissimilar, piecewise homogeneous, isotropic/anisotropic and linear elastic materials are much more complicated. For instance, the asymptotic stress field near the tips of an interfacial crack possesses an oscillatory singularity, which is not present for a crack in homogeneous, isotropic/anisotropic and linear elastic materials (e.g., [28,25]). The relationship between the stress intensity factors and the energy release rates for interfacial cracks between dissimilar isotropic elastic solids was presented and discussed by Mantic and Paris [23]. Some fundamental features of dynamic interfacial cracks were discussed by Willis [28], Yang et al. [32] and Wu [29] among others. The steady * Corresponding author. Tel.: +49 271 7402173; fax: +49 271 7404074. E-mail address: [email protected] (C. Zhang). URL: http://www.statics.uni-siegen.de (C. Zhang). 0013-7944/$ - see front matter Ó 2009 Elsevier Ltd. All rights reserved. doi:10.1016/j.engfracmech.2009.05.006

J. Lei et al. / Engineering Fracture Mechanics 76 (2009) 1996–2010

Nomenclature Latin symbols a characteristic crack-length compliance matrix aij ð1Þ

ð2Þ

ð3Þ

Ai ; Ai ; Ai Anj kl

Cabcd Cij cab cm cL Chi d, da, ea E1, E2 Ei g^ðÞ G12 Gnj ; Hnj kl kl ðDÞ

ðDÞ

Gkl ; Hkl

H(t) K |K| K1, iK2 K 1 ; iK 2 L Lbm, Pbm, na N pGab r r Rm ab s Shi sm t ta ta ua a u uGab Va m x xnl ; ynl y zm ; z0m

parameters in the asymptotic crack-tip displacement field

system matrix at (n-j)th time-step elasticity tensor elasticity matrix free-term matrix phase wave velocity wave velocity along the principal material axis ‘‘2” cosine integral parameters to compute the amplitude of the complex stress intensity factor Young’s modulus exponential integral Laplace transform of a function g() shear modulus time-domain system matrices at (n  j)th time-step Laplace-domain system matrices Heaviside function complex stress intensity factor amplitude of the complex stress intensity factor real part and complex part of the complex stress intensity factor normalized real part and complex part of the complex stress intensity factor length of the quarter-point element Qma matrices in the static fundamental solutions outward unit normal vector number of time-steps traction fundamental solutions polar coordinate parameter in Lubich’s quadrature formula matrix in the dynamic part of fundamental solutions Laplace-transform parameter sine integral complex variable in Lubich’s quadrature formula time traction vector prescribed traction vector displacement vector prescribed displacement vector displacement fundamental solutions mth eigenvector of the Christoffel-tensor source or field point unknown and known vectors at nth time-step observation point complex counterparts of the field and collocation points

Greek symbols Dirac-delta dab Dua crack-opening-displacements (CODs) Dt time-step e bi-material constant e error parameter in Lubich’s quadrature formula u phase angle of complex stress intensity factor quadratic shape-functions /i cm projection operator in the fundamental solutions ab C boundary of a domain CE external boundary Cu part of the external boundary with prescribed displacements

1997

1998

Ct CC C12 Cab

gi lm m12 q qc2m rab r0 xn-j X W, WR fm

J. Lei et al. / Engineering Fracture Mechanics 76 (2009) 1996–2010

part of the external boundary with prescribed tractions crack-faces interface between domain 1 and 2 Christoffel-tensor outward unit normal vector of the unit-circle complex roots of the characteristic equation Poisson’s ratio mass density eigenvalues of the Christoffel-tensor stress tensor amplitude of the impact loading weights in Lubich’s quadrature formula domain functions in the fundamental solutions complex variable in Lubich’s quadrature formula

Other symbols Riemann convolution * det [] determinant of the matrix [] o ()/ ox partial derivative of the function () boundary of the ith domain Xi oXi 0 () = d()/dz derivative of the function () f_ time derivative of the function f Re() real part sign signum function

state singular stress fields close to the tip of a crack growing along the interface between two anisotropic elastic solids were presented by Yang et al. [32]. However, for general transient dynamic interfacial crack problems in dissimilar, piecewise homogeneous, anisotropic and linear elastic materials, numerical methods such as the finite element method (FEM) or the boundary element method (BEM) are more feasible due to the mathematical complexity of the corresponding initial boundary value problems. BEM has certain advantages in fracture analysis and was successfully applied to static and dynamic crack analysis in homogeneous, isotropic and linear elastic materials since many years (e.g., [1,3,4,5,9,11,15]). Its extension to static analysis of interfacial cracks in dissimilar, piecewise homogeneous, anisotropic and linear elastic materials was done by Yuuki and Cho [33], Cho et al. [10], Tan et al. [26], Ang et al. [2], and Pan and Amadei [24]. In contrast, its extension to dynamic interfacial crack problems in dissimilar, piecewise homogeneous, anisotropic and linear elastic materials is not straightforward, since the required elastodynamic fundamental solutions for homogeneous, anisotropic and linear elastic materials have a quite complex mathematical structure and cannot be reduced to closed-form expressions. For this reason, only few investigations on dynamic interfacial crack problems in dissimilar, anisotropic and linear elastic materials by BEM have been reported in literature till now. Recently, a conventional time-domain BEM was developed for this purpose by Beyer [6] and Beyer et al. [7], who used a collocation method for both temporal and spatial discretizations. Another time-domain BEM was recently implemented by Wünsche et al. [30] and Wünsche [31], who applied a collocation method for the temporal discretization and a Galerkin-method for the spatial discretization. In both methods, time-domain elastodynamic fundamental solutions were required and the classical explicit time-stepping scheme was utilized. An important disadvantage of both methods is the fact that both methods are only conditionally stable and the accuracy of the numerical results could be strongly dependent on the choice of the time-steps. To circumvent this difficulty, a novel time-domain BEM for dynamic crack analysis in homogeneous, anisotropic and linear elastic materials was implemented by García-Sánchez and Zhang [12] and García-Sánchez et al. [14], who used the convolution quadrature formula of Lubich [21,22] for the temporal discretization and a collocation method for the spatial discretization. Unlike the conventional time-domain BEM, this method requires Laplace-domain instead of time-domain fundamental solutions. An essential advantage of this method is that it is more stable and less sensitive to the choice of the time-steps than the conventional time-domain BEM [13]. In this paper, the time-domain BEM developed by García-Sánchez and Zhang [12], García-Sánchez et al. [14] in conjunction with the sub-domain technique is applied to analyze the dynamic response of interfacial cracks in 2D, dissimilar, piecewise homogenous, anisotropic and linear elastic bi-materials. The bi-material system with an interfacial crack is divided into two homogeneous sub-domains along the interface. The displacement boundary integral equations (BIEs) are applied to each sub-domain. For the temporal discretization, the convolution quadrature formula of Lubich [21,22] is adopted, while a collocation method with quadratic shape functions is implemented for the temporal discretization. The present BEM is formulated in the time-domain, but requires only Laplace-domain instead of time-domain dynamic fundamental solutions. Time-dependent solutions can be obtained directly without an inverse Laplace-transform. Quarter-point elements are

1999

J. Lei et al. / Engineering Fracture Mechanics 76 (2009) 1996–2010

applied at the tips of the interfacial cracks. A displacement extrapolation method proposed by Cho et al. [10] is used to determine the complex dynamic stress intensity factors (SIFs). Numerical results are presented and discussed to verify the accuracy and the efficiency of the present time-domain BEM, and to show the effects of the material anisotropy on the dynamic overshoot and other relevant characteristics of the complex dynamic SIFs. 2. Problem statement and time-domain BIEs Let us consider a pre-existing interfacial crack in a 2D, piecewise homogeneous, anisotropic and linear elastic bi-material system as shown in Fig. 1. The deformation of the cracked system is either in a state of generalized plane strain or generalized plane stress. Both materials are assumed to be homogeneous, anisotropic and linear elastic. The ith domain Xi (i = 1, 2) is surrounded by @ Xi ¼ CiE þ C12 þ CC . Here C12 and CC denote the bonded interface and the interfacial crack-faces, respectively; CiE denotes the external boundary which is constituted by the part Cit with prescribed tractions ti and the remaining external  i , i.e., CiE ¼ Ciu þ Cit . The superscript i = 1 and 2 represents material 1 and 2, boundary Ciu with prescribed displacements u respectively. Both anisotropic materials satisfy the equations of motion

riab;b ¼ qi u€ia ;

ð1Þ

the Hooke’s law

riab ¼ C iabcd uic;d ;

ð2Þ

the initial conditions

uia ðx; 0Þ ¼ u_ ia ðx; 0Þ ¼ 0;

ð3Þ

the boundary conditions

tia ðx; tÞ ¼ t ia ðx; tÞ; x 2 Cit ;

ð4Þ

 ia ðx; tÞ; x 2 Ciu ; uia ðx; tÞ ¼ u

ð5Þ

ta ðx; tÞ ¼ 0; x 2 CC ;

ð6Þ

and the continuity conditions on the interface

t1a ðx; tÞ ¼ t2a ðx; tÞ u1a ðx; tÞ ¼ u2a ðx; tÞ;

x 2 C12 ;

ð7Þ i

i

i

where rab , ua and ta denote the stress, the displacement and the traction components; q is the mass density; C abcd is the fourth order elasticity tensor; a comma after a quantity designates spatial derivatives; while the superscript dots stand for temporal derivatives. Unless otherwise stated, the conventional summation rule over double indices is applied with Greek indices a, b, c, d = 1, 2 for the present 2D problem. Time-domain displacement boundary integral equations (BIEs) for homogeneous, anisotropic and linear elastic materials in conjunction with the sub-domain technique as used by Lei et al. [17–20] for piecewise homogeneous, isotropic and linear elastic bi-materials are adopted in the present analysis. To apply the time-domain BIEs for a homogeneous, anisotropic and linear elastic domain, the bi-material system is split from the interface into two homogeneous sub-domains X1 and X2 as depicted in Fig. 2. So the following displacement BIEs i

i

Z ciab ðxÞuib ðx; tÞ þ --

@ Xi

i pGi ab ðx; y; tÞ  ub ðy; tÞdsðyÞ ¼

Z @ Xi

i uGi ab ðx; y; tÞ  t b ðy; tÞdsðyÞ

ð8Þ

R can be applied to the boundary oXi of the ith sub-domain Xi (i = 1, 2). In Eq. (8), -- stands for the Cauchy-principal value inteGi gral, x = (x1, x2) and y = (y1, y2) represent the source and the observation points, pGi ab and uab are the 2D elastodynamic fundamental solutions for the tractions and the displacements; and an asterisk * denotes Riemann convolution which is defined by

Γ12

Γ1E x2

x3

Γ 2E

Ω

1

Γ1u

x1

ΓC

Ω2

Γ t2

Γ12

u Fig. 1. An interfacial crack in a 2D bi-material system.

t

2000

J. Lei et al. / Engineering Fracture Mechanics 76 (2009) 1996–2010

perfectly bonded interface Γ12 subregion 2

subregion 1

Γ1E

Ω1

traction free crack-faces Γ C

Ω2

Γ 2E

perfectly bonded interface Γ 12

Fig. 2. The sub-domain technique for a 2D bi-material system.

gðx; tÞ  hðx; tÞ ¼

Z

t

gðx; t  sÞhðx; sÞds:

ð9Þ

0

The term ciab ðxÞ in Eq. (8) is a constant matrix which depends on the smoothness of the boundary and reduces to 0.5dab for a smooth boundary, where dab is the Kronecker delta. The displacement fundamental solutions uGab and the traction fundamental solutions pGab are related by

pGab ¼ C bcfd

@uGaf nc ðxÞ; @yd

ð10Þ

where nc(x) represents the outward unit normal vector to the boundary at the collocation point. 3. 2D Laplace-domain fundamental solutions Due to the application of the quadrature formula of Lubich [21,22] for approximating the convolution integrals in the time-domain BIEs (8), the present analysis requires only the Laplace-domain dynamic fundamental solutions, which has been introduced by García-Sánchez and Zhang [12]. The Laplace-transform f ðsÞ of a time-dependent function f(t) is defined by

f ðsÞ ¼

Z

1

f ðtÞest dt;

ReðsÞ P 0;

ð11Þ

0

where s is the Laplace-transform parameter. 2D Laplace-domain displacement fundamental solutions for homogeneous, anisotropic and linear elastic solids can be written as [27]

uGab ðx; y; sÞ ¼

1 8p2

Z jgj¼1

2 X cmab Wðsjg  ðx  yÞj=cm ÞdSg ; qc2m m¼1

ð12Þ

where g = (g1, g2) is the wave propagation vector, and cm ab ¼ V am  V bm is the projection operator, with Vam and qcm being the eigenvectors and the eigenvalues of the Christoffel tensor Cab

Cab ðg1 ; g2 Þ ¼ C adbc gd gc :

ð13Þ

The function W(z) is given by

WðzÞ ¼ ½ez EiðzÞ þ ez EiðzÞ;

ð14Þ

where z is a complex variable and Ei(z) is the exponential integral [16]. Using the hyperbolic sine and cosine integrals, the function W(z) can also be expressed as

WðzÞ ¼ 2½sinhðzÞShiðzÞ  coshðzÞChiðzÞ:

ð15Þ

It can be easily shown that the function W(z) has a logarithmic singularity at |z| = 0. To deal with this singularity it is advantageous to split the displacement fundamental solutions into a singular static part and a regular dynamic part as GðRÞ

GðSÞ

uGab ðx; y; sÞ ¼ uab ðx; y; sÞ þ uab ðx; yÞ;

ð16Þ

where the regular dynamic part is given by GðRÞ

uab ðx; y; sÞ ¼

1 8p2

Z jgj¼1

2 X cmab R W ðs=cm ; jg  ðx  yÞjÞdSg ; qc2m m¼1

WR ðx; yÞ ¼ Wðx; yÞ þ 2 logðyÞ;

ð17Þ

ð18Þ

J. Lei et al. / Engineering Fracture Mechanics 76 (2009) 1996–2010

2001

while the singular static part has the following form: GðSÞ

uab ðx; yÞ ¼ 

Z

1 4p2

jgj¼1

2 X cmab log jg  ðx  yÞjdSg : qc2m m¼1

ð19Þ

The static displacement fundamental solution (19) can be reduced to an explicit form as GðSÞ

uab ðx; yÞ ¼ 

1

p

Re

" 2 X

# Pbm Q ma lnðzm  z0m Þ  lnði þ lm Þ ;

ð20Þ

m¼1

where

z0m ¼ x1 þ lm x2

zm ¼ y1 þ lm y2 ;

ð21Þ

are the complex counterparts of the collocation and the integration points, and Pbm, Qma and lm are determined by the anisotropic material constants which are given in Appendix A. by substituting Eqs. (17) and (19) into (10), the corresponding traction fundamental solutions pGab can be obtained as GðRÞ

GðSÞ

pGab ðx; y; sÞ ¼ pab ðx; y; sÞ þ pab ðx; yÞ;

ð22Þ

where GðSÞ

pab ðx; yÞ ¼

GðRÞ

1

p

pab ðx; y; sÞ ¼

Re

" 2 X

Lbm Q ma

m¼1

1 8p2

Z jgj¼1

lm n1 ðyÞ  n2 ðyÞ zm  z0m

# ;

2 X Rm ab s  w0 ðsjg  ðx  yÞj=cm Þ  sign½g  ðx  yÞdSg 2 c q c m m m¼1

ð23Þ

ð24Þ

with w0 = dw(z)/dz. The auxiliary functions Lbm and Rm ab are also given in Appendix A. 4. Numerical solution procedure To solve the strongly singular time-domain displacement BIEs (8), the quadrature formula of Lubich [21,22] is applied for the temporal discretization to approximate the Riemann convolution integrals in the time-domain BIEs, while a collocation method is used for the spatial discretization by using quadratic elements. 4.1. Time discretization In the convolution quadrature formula of Lubich [21,22], the Riemann convolution integral is approximated by

f ðtÞ ¼ gðtÞ  hðtÞ ¼

Z

t

gðt  sÞhðsÞds ffi

0

n X

xnj ðDtÞhðjDtÞ;

ð25Þ

j¼0

where the time t is divided into N equal time-steps Dt, and the weights xnj(Dt) are determined by

xnj ðDtÞ ¼

N1 r ðnjÞ X g^½dðfm Þ=Dte2piðnjÞm=N ; N m¼0

ð26Þ

in which g^ðÞ stands for the Laplace-transform of the function g(t), and

dðfm Þ ¼

2 X ð1  fm Þj =j;

fm ¼ re2pim=N ;

r ¼ e1=ð2NÞ

ð27Þ

j¼1

with e being the numerical error in computing the Laplace-transform g^ðÞ. More details on the convolution quadrature formula can be found in the original works of Lubich [21,22]. 4.2. Spatial discretization For the spatial discretization of the boundary oXi of each homogeneous sub-domain Xi, standard continuous quadratic elements are applied to the external boundary CiE , the interface C12 and the crack-faces CC except crack-tips, where quadratic quarter-point elements are adopted (see Fig. 3). A detailed description of the standard continuous quadratic shape functions can be found in recent works of García-Sánchez and Zhang [12] and García-Sánchez et al. [14]. Here we just outline the quarter-point element. According to Blandford et al. [8], the displacement components at any point inside the quarter-point element with a distance r from the crack-tip can be written as

2002

J. Lei et al. / Engineering Fracture Mechanics 76 (2009) 1996–2010

Crack- tip 3 2

1

2′ L

3′

4

2

3

L 4

L

L

Fig. 3. The quarter-point element near the tip of an interfacial crack.

ð1Þ

ð2Þ

ui ðrÞ ¼ Ai þ Ai

rffiffiffi   r ð3Þ r þ Ai L L

ð28Þ ð1Þ

ð1Þ

where L denotes the length of the quarter-point element and r is the distance from the crack-tip, Ai ¼ ui , ð2Þ ð1Þ ð2Þ ð3Þ ð3Þ ð1Þ ð2Þ ð3Þ ðkÞ Ai ¼ ð3ui þ 4ui  ui Þ and Ai ¼ ð2ui  4ui þ 2ui Þ, with ui being the displacement vector of the kth node of the quarter-point element. It is worthwhile to note that the quarter-point element is able to capture correctly the local asymptotic behavior near the tips of a crack in homogeneous, isotropic/anisotropic and linear elastic materials. However, it cannot pffiffiffi describe properly the local oscillatory behavior of the displacement field ua  r cosðe ln rÞ near the tips of an interfacial crack, where e is the bi-material constant. In the present analysis, the SIFs are computed by using a displacement extrapolation technique proposed by Cho et al. [10]. Indeed, the present approach can be refined by constructing a special crack-tip shape function to take the asymptotic oscillatory crack-tip field into account. An efficient way to achieve this was presented by Tan et al. [26] and Ang et al. [2] for interfacial cracks between dissimilar, piecewise homogeneous, anisotropic and linear elastic materials subjected to static loading. 4.3. Time-stepping scheme By applying the displacement BIEs (8) to each sub-domain, performing the temporal and spatial discretizations, and considering the interfacial continuity conditions (7), a system of linear algebraic equations for all discrete boundary quantities in both sub-domains X1 and X2 can be obtained as n X

j Hnj kl  ul ¼

j¼0

n X

j Gnj kl  t l ; ðn ¼ 0; 1; . . . ; N; k; l ¼ 1; 2; . . . ; EÞ;

ð29Þ

j¼0

nj nj where Hkl and Gkl are the time-domain system matrices, ujl is the vector containing the discrete boundary displacements, t jl is the corresponding vector containing the discrete boundary tractions, and E denotes the total number of the collocation points. According to Eq. (26), the system matrix at the (n  j)th time-step can be obtained as

Hnj kl ¼

N1 r ðnjÞ X ðDÞ H ðsm Þe2piðnjÞm=N þ dkl cl ; N m¼0 kl

ð30Þ

Gnj kl ¼

N1 r ðnjÞ X ðDÞ G ðsm Þe2piðnjÞm=N : N m¼0 kl

ð31Þ ðÞ

ðÞ

In Eqs. (30) and (31), sm ¼ dð1m Þ=Dt, cl = [cab(xl)]. The Laplace-domain system matrices Hkl and Gkl can be obtained by ðDÞ

Hkl ðsm Þ ¼

XZ el

 G  UdC; p

Ce

ðDÞ

Gkl ðsm Þ ¼

XZ el

 G  UdC; u

ð32Þ

Ce

where the summation applies to all elements containing the collocation point l, and the matrices are defined by

G ¼ u  U¼

 G11 u  G21 u

!  G12 u ;  G22 u

! G12 p ; G22 p  0 ; /3

G11 p G21 p

G ¼ p

/1

0

/2

0

/3

0

/1

0

/2

0

ð33Þ ð34Þ ðÞ

ðÞ

with ui (i = 1, 2, 3) being the spatial shape functions. It should be remarked here that Hkl is strongly singular while Gkl is weakly singular. Strongly singular integrals are computed analytically and numerically by a special regularization technique. Weakly singular integrals are computed numerically by using a logarithmic quadrature formula. Line-integrals arising in the

J. Lei et al. / Engineering Fracture Mechanics 76 (2009) 1996–2010

2003

dynamic fundamental solutions are computed numerically. More details on the numerical implementation of the time-domain BEM can be found in the recent works of García-Sánchez and Zhang [12] and García-Sánchez et al. [14]. By considering the boundary conditions (4)–(7), leaving all the unknown boundary quantities on the left-hand side, Eq. (29) can be rearranged as n X

j j Anj kl  xl ¼ yk ;

ð35Þ

j¼0 nj where Akl is the rearranged system matrix, xjl is the vector containing the unknown boundary data, and yjk is the vector containing the prescribed or known boundary quantities. By invoking the zero initial conditions (3), the following explicit timestepping scheme

xnk

¼

ðA0kl Þ1



ynl



n1 X

! Anj kl



xjl

ð36Þ

j¼1

can be obtained for computing the unknown boundary quantities at the nth time-step. Here ðA0kl Þ1 is the inverse of the system matrix A0kl at the time-step n = 0. The stability and the mesh-sensitivity of the present time-domain BEM have been investigated by García-Sánchez et al. [13]. 5. Computation of dynamic stress intensity factors For interface cracks between dissimilar, anisotropic and linear elastic materials, the local asymptotic crack-tip field can be described conveniently by a complex stress intensity factor, which can be written in the following form (e.g., [28,25])

KðtÞ ¼ jKðtÞjei/ðtÞ r ie ¼ K 1 þ iK 2 ;

ð37Þ

where e is the bi-material constant. In the case of two dissimilar, general anisotropic and linear elastic materials, the absolute value |K| and the phase angle u are given by [10]

rffiffiffiffiffiffiffi qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 þ 4e2 2p jKðtÞj ¼ K 21 þ K 22 ¼ 4 cosh pe r

rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi h i ðe2 Du1  e1 Du2 Þ2 þ ðd1 Du2  d2 Du1 Þ2 d1 e2  d2 e1

;

ð38Þ

2

tan u ¼

K 2 d  d1 ðDu2 =Du1 Þ ¼ ; K 1 e1 ðDu2 =Du1 Þ  e2

ð39Þ

where the parameters d, da, ea and e can be determined by the material constants and they are given in Appendix B. As mentioned in Section 4.2, the local oscillatory behavior of the displacement field near the tips of the interfacial crack, i.e., pffiffiffi ua  r cosðe ln rÞ, is not considered in the present numerical method. To compute the complex dynamic SIFs accurately in spite of this fact, the displacement extrapolation technique as suggested by Cho et al. [10] is utilized in the analysis. In the displacement extrapolation technique, numerical results of the displacements at element-nodes away from the cracktips are used, which avoids both the numerical errors and the oscillatory behavior of the displacements at the crack-tips. As shown by Cho et al. [10] for the static case, this simple technique can provide sufficiently accurate results for the complex SIFs. An application of special crack-tip elements by considering the local oscillatory behavior of the displacements at the crack-tip as suggested by Tan et al. [26] and Ang et al. [2] can certainly improve the numerical accuracy of the present time-domain BEM. However, it may not be very practically useful because of the complexity of the crack-tip elements and the fact, that the dominant zone of the oscillatory crack-tip field is anyway confined to an extremely small region of the crack-tip. 6. Numerical results and discussion To show the accuracy and the efficiency of the present time-domain BEM, several numerical examples are presented in this section. All examples presented here are performed for two bonded, rectangular, anisotropic and linear elastic bi-material plates subjected to an impact tensile loading, although the present time-domain BEM described above has no limitations on the plate geometry, crack configuration and loading conditions. Plane stress condition is assumed and the dynamic loading is taken as uniformly distributed impact loading of the form r(t) = r0H(t), where r0 is the loading amplitude and H(t) is the Heaviside step function. For convenience, the real part K1 and the imaginary part iK2 of the complex dynamic stress intensity factor for interfacial cracks are normalized by

K 1 ðtÞ ¼ K 1 ðtÞ=ðr0

pffiffiffiffiffiffi paÞ;



iK 2 ðtÞ ¼ iK 2 ðtÞ=ðr0

pffiffiffiffiffiffi paÞ;

where a is the length or the half-length of the crack.

ð40Þ

2004

J. Lei et al. / Engineering Fracture Mechanics 76 (2009) 1996–2010

6.1. An edge interfacial crack As the first example, we consider a rectangular bi-material plate with an edge interface crack subjected to an impact tensile loading r0 HðtÞ as shown in Fig. 4. A special case that the two materials have the same orthotropic material properties given by

E1 ¼ 118:30 GPa; E2 ¼ 54:80 GPa; G12 ¼ 8:79 GPa;

m12 ¼ 0:083; q ¼ 1900 kg=m3

is firstly considered to check the present time-domain BEM. This system is theoretically identical to the homogeneous case which has been studied by García-Sánchez et al. [14]. The geometric pffiffiffiffiffiffiffiffiffiffiffiffiffi parameters are taken as a = 12 mm, W = 52 mm and H = 20 mm. A time-step cLDt = H/50 is selected, where cL ¼ C 22 =q is the wave velocity along the material-axis with E2. The normalized real part of the complex dynamic SIF is plotted in Fig. 5. The corresponding numerical result for an edge crack of length a in a homogeneous, orthotropic and linear elastic material obtained by García-Sánchez et al. [14] is also given for the comparison purpose. A comparison of both numerical results shows a very good agreement, which partially verifies the correctness and the accuracy of the present time-domain BEM. It should be noted here that the crack-tip field in a 2D, homogeneous, anisotropic and linear elastic material is characterized by the mode-I and the mode-II SIFs. In the present example, the mode-II dynamic SIF is identically zero, because the considered material is orthotropic and the crack-line coincides with one of the material axes. In this special case, the mode-I dynamic SIF is equal to the real part of the complex dynamic SIF. Fig. 5 shows a dynamic overshoot of the SIF over its corresponding static value. The maximum dynamic K 1 is 3.26 and the corresponding static value obtained by ANSYS is 1.72. So the dynamic overshoot of K 1 is about 89.5%. Next, a different material combination is considered in the same bi-material plate as shown in Fig. 4. In this calculation, the geometric parameters are taken as a = 5 mm, W = 30 mm and H = 20 mm. A graphite–epoxy composite with a composition of 65% graphite and 35% epoxy is investigated, which has a mass density of q1 = q2 = 1600 kg/m3 and the following elastic constants:

2

155:43

6 C ij ¼ 4 3:72 0:0

3:72 16:34 0:0

0:0

3

7 0:0 5 GPa: 7:48

ð41Þ

Fully anisotropic material properties in plane stress are investigated by using different values of the inclination angle h between the principal material-axis 1 and the crack-line. The elastic constants in Eq. (41) are transformed by the rotation of the angles h1 = 30° and h2 = 60°. Each sub-domain is discretized into 39 elements, where 18 elements are used for the external boundary and 21 elements for the bi-material interface including 10 for the crack-face. A time-step of Dt = 0.2 ls is chosen. The numerical results of the present time-domain BEM and the FEM using ANSYS are presented in Fig. 6 as a function of the  time. The real part K 1 and the imaginary part iK 2 of the complex dynamic SIF show a very good agreement between the numerical results of the present time-domain BEM and that of the FEM. Fig. 6 implies that the amplitudes of the real and the imaginary parts of the complex dynamic SIF have typical peaks and dips, which are related to wave propagation and  reflections on the boundaries and interfaces of the cracked plate. The static K 1 and iK 2 obtained by ANSYS are given as   1.635 and 0.691, while the maximum values of the corresponding dynamic K 1 and iK 2 provided by the time-domain BEM are 2.956 and 1.37. Therefore, their dynamic overshoots are 80.82% and 98.15%, respectively, which shows the significant effects of the interface and the material anisotropy on the complex dynamic SIF.

σ (t )

W

H

a

Cij1 ,θ1 Cij2 ,θ2

H

y

x σ (t ) Fig. 4. An edge interfacial crack in a layered rectangular plate.

2005

J. Lei et al. / Engineering Fracture Mechanics 76 (2009) 1996–2010

3.5 3.0

2.0

*

K I (t)

2.5

1.5 1.0 Garcia-Sanchez et al. (2008) This work

0.5 0.0

0

2

4

6

tcL/H

8

10

12

Fig. 5. Normalized dynamic SIF for an edge crack in a homogeneous and anisotropic plate.

3.2

K1

2.4

*

*

K (t)

1.6 0.8

FEM This work

0.0 -0.8 -1.6

iK2 0

4

8

*

12

16

20

24

28

t [μs] Fig. 6. Normalized complex dynamic SIF for an edge interfacial crack.

6.2. A central interfacial crack Fig. 7 shows a rectangular, anisotropic and linear elastic bi-material plate with a central interfacial crack. The geometric parameters are fixed as 2a = 4.8 mm and W = H = 20 mm. The bi-material plate is composed of two different graphite–epoxy composites and they have the following elastic constants

2

122:77

6 C 1ij ¼ 4 3:88 0:0

3:88 16:34 0:0

0:0

3

7 3 0:0 5 GPa; q1 ¼ 1600 kg=m 6:94

ð42Þ

and

2

C 2ij

3 65:41 4:29 0:0 3 6 7 ¼ 4 4:29 16:34 0:0 5 GPa; q2 ¼ 1600 kg=m 0:0 0:0 5:58

ð43Þ

respectively. In our computation, each sub-domain is discretized into 54 elements, where 30 elements are used for the external boundary and 24 elements for the bi-material interface including eight for the crack-face. A time-step of Dt = 0.11 ls is chosen. The present numerical results and other results using a conventional time-domain BEM by Beyer [6] and Beyer et al.

2006

J. Lei et al. / Engineering Fracture Mechanics 76 (2009) 1996–2010

σ (t )

W

H 2a

Cij1 ,θ1 Cij2 ,θ2

y H

x σ (t ) Fig. 7. A central interfacial crack in a layered rectangular plate.

4

K1

2

*

K (t)

3

Beyer et al. (2008) This work

1

iK2

0 0

4

8

12

16

20

24

t [μs] Fig. 8. Normalized complex dynamic SIF for a central interfacial crack.



[7] are shown in Fig. 8. The normalized real part K 1 and the normalized imaginary part iK 2 of the complex dynamic stress intensity factor show a very good agreement between our results and the results of Beyer [6] and Beyer et al. [7]. In contrast to the edge interfacial crack, the real part of the complex dynamic SIF exhibits two distinct peaks. On the other hand, the imaginary part of the complex dynamic SIF is very small, because the considered materials are orthotropic and the loading is a pure tensile loading perpendicular to the interfacial crack. Here again, a dynamic overshoot of the SIFs over their corre sponding static values is observed in Fig. 8. The maximum dynamic K 1 and iK 2 are 3.63 and 0.172, while the corresponding static values obtained by ANSYS are 1.473 and 0.0385. This corresponds to a dynamic overshoot of 146.44% for K 1 and  346.75% for iK 2 . 6.3. A pair of interacting collinear interfacial cracks In the last numerical example, we consider a pair of interacting collinear interfacial cracks in a rectangular, anisotropic and linear elastic bi-material plate as shown in Fig. 9. The cracked bi-material plate has the following geometrical data: H = 15 mm, W = 20 mm, 2a = 4 mm and d = 12 mm. The temporal discretization is performed by using a time-step of Dt = 0.11 ls. The external boundary and the interface are approximated by elements of a length 2.0 mm and both interfacial cracks are divided into 10 elements. The mass densities of the materials are taken as q1 = q2 = 1600 kg/m3 and the elastic constants of the upper and the lower sub-domains are given by Eq. (41) with the rotation angles h1 = 30° and h2 = 60°.

J. Lei et al. / Engineering Fracture Mechanics 76 (2009) 1996–2010

2007

σ (t )

W

W

H 2a

d

2a Cij1 ,θ1 Cij2 ,θ2

y H

x σ (t ) Fig. 9. A pair of collinear interfacial cracks in a layered rectangular plate.

2.5 2.0

K1

-

1.0

*

K (t)

1.5

FEM This work Wünsche et al. (2008)

0.5 0.0 -0.5

iK2

-1.0 -1.5

0

5

-

10

15

20

t [μs] 2.5 2.0 1.5

K1

+

*

K (t)

1.0 FEM This work Wünsche et al. (2008)

0.5 0.0 -0.5

+

iK2

-1.0 -1.5

0

5

10

15

20

t [μs] Fig. 10. Normalized complex dynamic SIFs for the left interfacial crack.

The numerical results for the left and the right interfacial cracks are presented in Figs. 10 and 11 versus the time. The normalized real and imaginary parts of the complex dynamic SIF for the left () and the right (+) crack-tips obtained by the

2008

J. Lei et al. / Engineering Fracture Mechanics 76 (2009) 1996–2010

2.5 2.0 1.5

K1

*

K (t)

1.0 0.5

FEM This work Wünsche et al. (2008)

0.0 -0.5

iK2

-1.0 -1.5

-

0

-

5

10

t [μs]

15

20

2.5 2.0 +

K1

1.5

*

K (t)

1.0 0.5

FEM This work Wünsche et al. (2008)

0.0 -0.5

iK2

-1.0 -1.5

0

5

+

10

15

20

t [μs] Fig. 11. Normalized complex dynamic SIFs for the right interfacial crack.

Table 1 Dynamic overshoots of the SIFs for collinear interfacial cracks. Crack

Crack-tip

K 1 (%)

iK 2 (%)

Left crack

Left tip Right tip

77.03 86.43

109.36 104.34

Right crack

Left tip Right tip

98.21 106.03

97.17 99.43

present time-domain BEM, the time-domain BEM of Wünsche [31] and the FEM results using ANSYS show a very good agreement. Figs. 10 and 11 imply that the material anisotropy has significant effects on the complex dynamic SIFs with respect to their maximum/minimum amplitudes and their time-dependence. The left and the right interfacial cracks possess quite different dynamic SIFs, although the geometry of the bi-material plate and the external loading are symmetric with respect to the interfacial cracks. The material anisotropy of the bi-material plate leads to an asymmetry in the dynamic response of the left and the right interfacial cracks. For the left interfacial crack, the real part of the complex dynamic SIFs for both cracktips shows two distinct peaks in the considered time interval, while the corresponding real part of the complex dynamic SIFs for the right interfacial crack presents several picks. In contrast, the differences between the imaginary parts of the complex dynamic SIFs for the left and the right interfacial cracks are relatively small. Furthermore, the material anisotropy and the crack interaction result in a notable imaginary part of the complex dynamic SIFs, which is not negligible and may play an important role in the propagation and kinking of the interfacial cracks. The dynamic overshoots of the SIFs are summarized in Table 1, which vary from 77.03% to 109.36%.

J. Lei et al. / Engineering Fracture Mechanics 76 (2009) 1996–2010

2009

7. Conclusions This paper presents a transient dynamic analysis of interfacial cracks in 2D, piecewise homogeneous, generally anisotropic and linear elastic bi-materials. A time-domain BEM together with the sub-domain technique is applied for this purpose. For spatial discretization, a collocation method with quadratic elements is implemented, while the convolution quadrature formula of Lubich [21,22] is applied for the temporal discretization. The present time-domain BEM requires only the Laplace-domain fundamental solutions and yields directly time-dependent solutions without an inverse Laplace-transform. It is very general and has no limitations on the material anisotropy, the geometry of the bi-material systems, the configuration of the interfacial cracks and the dynamic loading conditions. For simplicity, straight quarter-point elements are implemented at the crack-tips, which can only capture the square-root but not the oscillatory behavior of the displacement field near the tips of an interfacial crack. Complex dynamic stress intensity factors are computed by using a displacement extrapolation technique. The accuracy and the efficiency of the present time-domain BEM are verified by reference solutions obtained by FEM and other time-domain BEM. Numerical results imply that the material anisotropy has significant influences on the dynamic overshoots and other characteristics of the complex dynamic SIFs. Acknowledgements This work is supported by the German Research Foundation (DFG) under the Project Nos. ZH 15/5-1 and ZH 15/5-2, which is gratefully acknowledged. Appendix A. Auxiliary functions In Eqs. (20) and (21), lm are the complex roots of the following characteristic equation

det½C 1ij1 þ ðC 1ij2 þ C 2ij1 Þlm þ C 2ij2 l2m  ¼ 0;

ðA:1Þ

and the matrices Pjm and Qmi are determined by

Cij ð1; lm ÞP jm ¼ 0;

ðA:2Þ

Q ¼ P1 ðB1 þ B1 Þ1 ;

ðA:3Þ

B ¼ iPL1 ;

ðA:4Þ

where

L ¼ ½Lim ;

Lim ¼ ðC 2ij1 þ C 2ij2 lm ÞPjm ¼ 

1

lm

ðC 1ij1 þ C 1ij2 lm ÞPjm ;

ðno sum over mÞ:

ðA:5Þ

The matrix Rm ab in Eq. (24) is determined by m Rm ab ¼ C bdkj caj nd ðxÞgk :

ðA:6Þ

Appendix B. Parameters for computing the dynamic SIFs In Eqs. (37)–(39), the parameters related to material constants are given by

da ¼ Ga ðcos h þ 2e sin hÞ þ Pa ðsin h  2e cos hÞ

ðB:1Þ

ea ¼ Ga ð2e cos h  sin hÞ þ P a ðcos h þ 2e sin hÞ; d¼



H11 H22 ð1  12 Þ  H212 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ; 2 H11 H22  H212

G1  iP 1



G2  iP 2

2

H12 6 H22

¼ d4

i

ðB:2Þ

H11 H22 þH212 H22

pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 H11 H22 H12

2H12 1  i pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2

3 7 5;

ðB:3Þ

H11 H22 H12

H11 ¼ ½a11 ðn1 þ n2 ÞI þ ½a11 ðn1 þ n2 ÞII ; " H22 ¼ a22

n1

a21 þ n21

þ

n2

a22 þ n22

!#

" þ a22 I

ðB:4Þ n1

a21 þ n21

þ

n2

a22 þ n22

!# ; II

ðB:5Þ

2010

J. Lei et al. / Engineering Fracture Mechanics 76 (2009) 1996–2010

H12 ¼ ½a11 ða1 n2 þ a2 n1 ÞI þ ½a11 ða1 n2 þ a2 n1 ÞII ;

ðB:6Þ

1¼

½a11 ða1 a2  n1 n2 Þ  a12 I  ½a11 ða1 a2  n1 n2 Þ  a12 II pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ; H11 H22

ðB:7Þ

e¼

1 1b ln ; 2p 1 þ b

pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi H11 H22 b ¼ 1 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ; H11 H22  H212

ðB:8Þ

r h ¼ e ln : a

In Eqs. (B.1)–(B.8), a is the length of the crack, aij is the compliance matrix, ak and nk are the real part and the imagine part of the complex root lk, and the subscripts I and II in Eqs. (B.4)–(B.7) indicate material I and material II, respectively. References [1] Aliabadi MH. The boundary element method, Applications in solids and structures. vol. 2. Computational Mechanics Publications, John Wiley & Sons; 2002. [2] Ang HE, Torrance JE, Tan CL. Boundary element analysis of orthotropic delamination specimens with interface cracks. Engng Fract Mech 1996;54:601–15. [3] Balas J, Sladek J, Sladek V. Stress analysis by boundary element methods. Amsterdam: Elsevier; 1989. [4] Beskos DE. Boundary element methods in dynamic analysis. Appl Mech Rev 1987;40:1–23. [5] Beskos DE. Boundary element methods in dynamic analysis: part II (1986–1996). Appl Mech Rev 1997;50:149–97. [6] Beyer S. A time-domain BEM for transient crack analysis in anisotropic materials. Ph.D. Thesis, TU Freiberg, Germany; 2008 [in German]. [7] Beyer S, Zhang Ch, Hirose S, Sladek J, Sladek V. Transient dynamic analysis of interface cracks in 2-D anisotropic elastic solids by a time-domain BEM. In: Abascal R, Aliabadi MH, editors. Advances in boundary element techniques IX. UK: EC Ltd.; 2008. p. 419–26. [8] Blandford GE, Ingraffea AR, Liggett JA. Two-dimensional stress intensity factors computations using the boundary element method. Int J Numer Methods Engng 1981;17:387–404. [9] Bonnet M. Boundary integral equation methods for solids and fluids. John Wiley & Sons Ltd.; 1999. [10] Cho SB, Lee KB, Choy YS, Yuuki R. Determination of stress intensity factors and boundary element analysis for interface cracks in dissimilar anisotropic materials. Engng Fract Mech 1992;43:603–14. [11] Dominguez J. Boundary elements in dynamics. Southampton, London: Computational Mechanics Publications; 1993. [12] García-Sánchez F, Zhang Ch. A comparative study of three BEM for transient dynamic crack analysis of 2-D anisotropic solids. Comput Mech 2007;40:753–69. [13] García-Sánchez F, Zhang Ch, Sáez A. Mesh-sensitivity analysis of dynamic BEM for cracked anisotropic solids. In: Minutolo V, Aliabadi MH, editors. Advances in boundary element techniques VIII. UK: EC Ltd.; 2007. p. 31–6. [14] García-Sánchez F, Zhang Ch, Sáez A. A 2-D time-domain BEM for dynamic crack problems in anisotropic solids. Engng Fract Mech 2008;75:1412–30. [15] Gaul L, Kögl MK, Wagner M. Boundary element methods for engineers and scientists. Berlin, Heidelberg: Springer-Verlag; 2003. [16] Gradshteyn IS, Ryzhik IM. Table of integrals, series, and products. 6th ed. Academic Press; 2000. [17] Lei J, Wang YS, Gross D. Dynamic interaction between a sub-interface crack and the interface in a bi-material: time-domain BEM analysis. Arch Appl Mech 2003;73:225–40. [18] Lei J, Wang YS, Gross D. Time-domain BEM analysis of a rapidly growing crack in a bi-material. Int J Fract 2004;126:103–21. [19] Lei J, Wang Y-S, Gross D. Numerical simulation of crack kinking from an interface under dynamic loading by time domain boundary element method. Int J Solids Struct 2007;44:996–1012. [20] Lei J, Wang Y-S, Gross D. Numerical simulation of crack deflection and penetration at an interface in a bi-material under dynamic loading by timedomain boundary element method. Int J Fract 2008;149:11–30. [21] Lubich C. Convolution quadrature and discretized operational calculus. Part I. Numer Math 1988;52:129–45. [22] Lubich C. Convolution quadrature and discretized operational calculus. Part II. Numer Math 1988;52:413–25. [23] Mantic V, Paris F. Relation between SIF and ERR based measures of fracture mode mixity in interface cracks. Int J Fract 2004;130:557–69. [24] Pan E, Amadei B. Boundary element analysis of fracture mechanics in anisotropic bimaterials. Engng Anal Bound Elem 1999;23:683–91. [25] Suo Z. Singularities, interfaces and cracks in dissimilar anisotropic media. Proc Roy Soc Lond, Ser A, Math Phys Sci 1990;427:331–58. [26] Tan CL, Gao YL, Afagh FF. Boundary element analysis of interface cracks between dissimilar anisotropic materials. Int J Solids Struct 1992;29:3201–20. [27] Wang C-Y, Zhang Ch. 3-D and 2-D dynamic Green’s functions and time-domain BIEs for piezoelectric solids. Engng Anal Bound Elem 2005;29:454–65. [28] Willis JR. Fracture mechanics of interfacial crack. J Mech Phys Solids 1971;19:353–68. [29] Wu KC. Explicit crack-tip fields of an extending interface crack in an anisotropic bimaterial. Int J Solids Struct 1991;27:455–6. [30] Wünsche M, Zhang Ch, Sladek J, Sladek V, Hirose S. Interface crack in anisotropic solids under impact loading. Key Engng Mater 2007;348–349:73–6. [31] Wünsche M. A Galerkin boundary element method for the solution of dynamic crack problems in anisotropic materials. Ph.D. Thesis, TU Freiberg, Germany; 2008 [in German]. [32] Yang W, Suo Z, Shih CH. Mechanics of dynamic debonding. Proc Roy Soc Lond A 1991;433:679–97. [33] Yuuki R, Cho SB. Efficient boundary element analysis of stress intensity factors for interface cracks in dissimilar materials. Engng Fract Mech 1989;34:179–88.