- Email: [email protected]
BIEM for 2D steady-state problems in cracked anisotropic materials
Engineering Analysis with Boundary Elements 29 (2005) 689–698 www.elsevier.com/locate/enganabound
BIEM for 2D steady-state problems in cracked anisotropic materials P. Dinevaa, T. Rangelovb,*, D. Grossc b
a Institute of Mechanics, Bulgarian Academy of Sciences, 1113 Sofia, Bulgaria Institute of Mathematics and Informatics, Bulgarian Academy of Sciences, 1113 Sofia, Bulgaria c Institute of Mechanics, Darmstadt University of Technology, 64289 Darmstadt, Germany
Received 28 October 2004; revised 17 February 2005; accepted 17 February 2005 Available online 21 April 2005
Abstract The two-dimensional ‘in-plane’ time-harmonic elasto-dynamic problem for anisotropic cracked solid is studied. The non-hypersingular traction boundary integral equation method (BIEM) is used in conjunction with closed form frequency dependent fundamental solution, obtained by Radon transform. Accuracy and convergence of the numerical solution for stress intensity factor (SIF) is studied by comparison with existing solutions in isotropic, transversely-isotropic and orthotropic cases. In addition a parametric study for the wave field sensitivity on wave, crack and anisotropic material parameters is presented. q 2005 Elsevier Ltd. All rights reserved. Keywords: Steady-state elastodynamics; Cracked anisotropic material; BIEM; SIF computation AMS Subject Classification: 74E10; 74S15; 74R10; 74H35
1. Introduction In many cases, engineering and construction materials such as modern fibre-reinforced composite, many crystals, piezoelectric materials, ice, wood, diamond, etc. possess both anisotropy and crack-like defects induced due to manufacturing or in-service conditions. Elastodynamic anisotropic crack analysis has attracted the attention of scientists in such fields as ultrasonic non-destructive evaluation for material quality inspection or seismology and exploration geophysics where observation of elastic waves generated by earthquakes or explosions may reveal the presence of structural discontinuities in the earth’s interior. Studies on crack problems in anisotropic materials have been done since many years. The initial groundwork was proposed by Sih et al. [1] and Sih and Liebowitz [2]. Lekhnitski’s [3,4] and Savin’s [5] systematic formulations of anisotropy are regarded as a foundation for the anisotropic elasticity. The methods used for the solution
* Corresponding author. Tel.: C359 2 979 2845; fax: C359 2 971 3649. E-mail address: [email protected] (T. Rangelov).
0955-7997/$ - see front matter q 2005 Elsevier Ltd. All rights reserved. doi:10.1016/j.enganabound.2005.02.006
of anisotropic crack problems are semi-analytical in Ohyoshi [6,7], Dhawan [8,9], Karim and Kundu [10,11], Kundu and Bostrom [12], Sarkar et al. [13]; numerical as finite-difference schemes in Hua et al. [14]; finite-element method in Brebbia and Venturini [15], Pageau and Biggers [16], Pageau et al. [17], BIEM in Snyder and Cruse [18], Chan and Cruse [19], Zhao et al. [20], Tan and Gao [21], Aliabadi and Brebbia [22], Chandra et al. [23], Pan and Amadei [24], Saez et al. [25], Ang et al. [26], Zhang [27], Saez and Dominguez [28–30], Abuquerque et al. [31], Wang and Achenbach [32–34] and hybrid methods in Song and Wolf [35]. In most of the semi-analytical papers various types of integral equation techniques are used. Fourier or Laplace transforms reduce the boundary-value problem to dual integral equations that are expressed in terms of Fredholm integral equations. Using this approach Ohyoshi [6,7] solved the anti-plane and the in-plane diffraction problems of waves in a cracked orthotropic plane. Dhawan [8,9] treated the same problem but for the transversely-isotropic case. Karin and Kundu [10,11] considered the dynamic response of a layered anisotropic half-space with anti-plane interfacecracks and an orthotropic half-space with a subsurface inplane crack. Kundu and Bostrom [12] studied wave scattering by a circular crack in a transversely-isotropic
690
P. Dineva et al. / Engineering Analysis with Boundary Elements 29 (2005) 689–698
solid and the dynamic response of three co-planar cracks in an infinite orthotropic medium in [13]. Although BIEM is a very well suited approach for the type of problems considered here, only a few works can be found for dynamic anisotropic crack problems. The main difficulty when applying BIEM to anisotropic dynamic crack problems is the availability of a fundamental solution that can be efficiently evaluated and implemented in a BEM code. The fundamental solution for 2D elastostatics obtained by Cruse and Swedlow [36] is applied by Abuquerque et al. [31] to transient anisotropic problems. Using multi-domain BIEM formulation, dual reciprocity method and traction singular quarter-point boundary elements (TSQP-BE) in [31] is determinated SIFs for a crack in a quasi-isotropic infinity long strip and for a rectangular orthotropic plate with a central crack and with a slanted edge crack. Wang and Achenbach [32–34] derived fundamental solutions by Radon transforms for general anisotropic solids. These fundamental solutions are related to transient and time-harmonic cases of 3D and 2D out of plane deformation states and they have been used also in Zhang [27] and in Saez and Dominguez [28–30]. The main aim of the present paper is the development of a non-hypersingular traction-based BIEM for solving 2D steady-state elastic in-plane wave scattering problems in plane anisotropic solids. The procedure is based on the non-hypersingular traction BIE proposed by Zhang and Gross [37] for the isotropic case, in conjunction with a closed form 2D fundamental solution, obtained by Radon transform and with crack-tip quarter-point boundary elements (QP-BE) implemented in the used parabolic approximation. The novelty of the present work lies in the generality of the proposed methodology and in the parametric study results revealing how the anisotropy of the material, the existence of cracks, the wave parameters as frequency and incident angle and mutual disposition of the crack line, axis of the material symmetry and wave propagation direction influence on the dynamic stress concentration fields. The advantages of the proposed method as being BIEM include all its important features as reduction of the dimension of the problem, implicit satisfaction of radiation condition at infinite, semi-analytical character of the method, high accuracy at solution of problems with stress gradients, in addition: (a) The method is valid for general type of anisotropy. It allows considering the cases when the coordinate axes do not coincide with the axis of material symmetry. (b) The proposed numerical procedure works for arbitrary incident wave angle, while in the most cases the used for the purpose methods as singular integral method has been applied only for the normal incident wave. (c) ‘In-plane’ crack is used in the proposed mechanical model, while in the most cases in the literature the simpler ‘anti-plane’ crack case is considered.
(d) The method has a potential to be applied for transient dynamic loads, after applying the Laplace or Fourier transformation. The paper is organized as follows. The setting of the boundary-value problem and steady-state non-hypersingular traction BIE formulation is given in Section 2. Fundamental solution in a closed form is presented in Section 3. Paragraph 4 focuses on the numerical solution procedure. Some results concerning the accuracy and the convergence of the proposed numerical scheme together with a parametric study on wave-crack interaction in anisotropic media are discussed in Section 5. The paper is finishing by some concluding remarks in Section 6.
2. Problem statement and traction BIM formulation We consider an infinite, homogeneous, anisotropic and linearly elastic solid containing a finite line crack G with a length 2a, see Fig. 1. The solid is subjected to an incident time-harmonic plane L-wave, and the deformation of the solid is as in a state of a plane strain. The non-zero quantities are the displacements components ui and the stresses components sij, i,jZ1,2. In cases of general anisotropy, the in-plane and the out-of-plane displacement are usually coupled. In contrast, we will assume, see Xu [38], that the material possesses a mid-plane symmetry such as the inplane (x1, x2 coordinate plane) and the out-of-plane deformations are decoupled. In the following we consider the in-plane problem, described by the equation of motion sij;j þ ru2 ui ¼ 0 in R2 nG;
(1)
Hooke’s law sij Z Cijkl uk;l
(2)
and the traction-free boundary condition on the crack-faces tj Z 0 on G:
Fig. 1. Incident L-wave in a cracked anisotropic plane.
(3)
P. Dineva et al. / Engineering Analysis with Boundary Elements 29 (2005) 689–698
Here r is the mass density, u is the frequency, tiZsijnj is the traction on the crack boundary GZGCgGK, ni are the components of the unit normal vector on GC and Cijkl are the elastic constants that satisfy the conditions
691
and the traction is obtained by using the Hooke’s law (2) ð tpsc ðx; uÞ ZKCpqkl nq sijk;l ðx; x; uÞDui ðx; uÞnj dGx ; x;GC GC
Cijkl Z Cjikl Z Cklij ; Cijkl gij gkl O 0 for every non zero real symmetric tensor gij :
(9) (4)
All subscripts vary 1, 2 and the conventional summation rules over double indices is implied. Finally, it is assumed that the displacement ui satisfies Sommerfeld type radiation condition at infinity. Substituting Eq. (2) into Eq. (1) the equation of motion can be written as c11 u1;11 C c66 u1;22 C 2c16 u1;12 C c16 u2;11 C c26 u2;22 Cðc12 C c66 Þu2;12 C ru2 u1 Z 0 c u C c u C ðc C c Þu C c u 16 1;11 26 1;22 12 66 1;21 66 2;11 Cc22 u2;22 C 2c26 u2;12 C ru2 u2 Z 0:
GC
(5)
sc ui ðx; uÞ Z uin i ðx; uÞ C ui ðx; uÞ
ð6Þ
where xZ(x1, x2) and u is the frequency. The functions in uin i ; ti are the displacement and traction components of the incident wave field in the absence of the crack, while sc functions usc i ; ti are the corresponding displacement and traction of the scattered field due to the interaction of the incident wave with the crack. It is assumed the incident wave to be known, while the scattered wave field is unknown and has to be determined. The scattered wave would satisfy the equation of motion (5) and the boundary condition (3), which can be rewritten as tisc Z Ktiin on G
(7)
The scattered displacement usc i ðx; uÞ can be represented as in Zhang and Gross [37] by the boundary integral ð ðx; uÞ Z K skij ðx; x; uÞDuk ðx; uÞni dGx ; x ;GC (8) usc j GC
Kru2 Cijkl nj ðxÞ KCijkl nj ðxÞ
Constants cmn (m,nZ1,2,6) denote the anisotropic material parameters in abbreviated notation, obtained from Cijkl following the rules (11)41, (22)42, (12Z21)46. There are altogether six independent constants. In the orthotropic case with the elasticity axes parallel to the coordinate axes c16Zc26Z0 is fulfilled and the independent material constants are four: c11, c12, c22, c66, see Su and Sun [39]. When the coordinate axis and the elasticity axis do not coincide there are 6 material constants. The interaction of an incident wave with the crack induces scattered waves. The total wave field can be presented as a sum of the incident wave field and the scattered wave field:
and ti ðx; uÞ Z tiin ðx; uÞ C tisc ðx; uÞ
Here sijk ðx; x; uÞ is the stress Green’s function and Dui Z ui jGC K ui jGK is the crack-opening displacement (COD). Following the approach in Zhang and Gross [37], for the isotropic case, the BIE for the anisotropic case is obtained by using the boundary condition (7) and the radiation condition at infinity: ð in ti ðx;uÞZCijkl nj ðxÞ sphk ðx;x;uÞDup;h ðx;uÞdll nl dGx
ð
ð
Udk ðx;x;uÞDud ðx;uÞdll nl dGx
GC
smlk ðx;x;uÞDum;l ðx;uÞnl dGx ; x2GC
GC
(10) Eq. (10) is integro-differential equation with respect to COD Dui and the singular integrals converge in CPV sense, if a priory smoothness requirements Dui2C1Ca (G) are fulfilled, see Rangelov et al. [40].
3. Fundamental solution, its derivatives and asymptotics Since Eq. (5) is with constant coefficients, there exists a fundamental solution U Z fUkj g, see John [41] that is a solution of the equation ðCij ðvÞ C dij ru2 ÞUjk ðx; x; uÞ Z Kdik dðx K xÞ
(11)
where Cij ðvÞZ Cimjn vxm vxn . Let’s note that in the isotropic case: C11 ðvÞZ ðlC 2mÞv2x1 C mv2x2 ; C12 ðvÞZ C21 ðvÞZ ðlC mÞvx1 vx2 ; C22 ðvÞZ ðlC 2mÞv2x2 C mv2x1 ; where l, m are the Lame´ constants. The Radon transform, see Ludwig [42] and Gelfand and Shilov [43] is of fundamental importance in order to work out the analytic solution of the problem at hand. Let f(x) be a function defined in R2 and s is a real number. The Radon transform R of function f(x) is defined as: ð Ð ^f ðs; nÞ Z R½f ðxÞ Z f ðxÞ dU Z f ðxÞdðs K hn; xiÞ dx; hn;xiZs
(12) where d is Dirac delta function, h,i means scalar product in R2. The Radon transform is an integration of f(x) over all planes defined by hn, xiZs. The inverse Radon transform in
692
P. Dineva et al. / Engineering Analysis with Boundary Elements 29 (2005) 689–698
2D can be written as: ð 1 Kðf^ðs; nÞÞj dn; f ðxÞ Z 2 4p
where: ðCN
vs f^ðs; nÞ ds: Kðf^Þ Z KN s K s
jnjZ1
(13) The next Radon transform properties will be used: R(d(x))Zd(s) and f^ðas; anÞZ a1 f^ðs; nÞ; Rða1 f1 C a2 f2 ÞZ a1 f^1 C a2 f^2 ; Rðvj f ðxÞÞZ nj vs f^ðs; nÞ: Let us apply the Radon transform to both side of Eq. (11). For simplicity xZ(0,0) is considered:
ðCij ðnÞv2s C ru2 dij ÞU^ jk ðs; n; uÞ Z Kdik dðsÞ
(14)
For example, in the orthotropic case the form of Cij(n) is: ! c11 n21 C c66 n22 ðc12 C c66 Þn1 n2 fCij ðnÞg Z (15) ðc12 C c66 Þn1 n2 c22 n22 C c66 n21 Eq. (14) is a linear system of ordinary differential equations. For solving it we transform it into the canonical form. We can consider n, jnjZ1 since only such n is used in the inverse Radon transform (13). The 2!2 matrix CZ{Cij(n)} is positive and symmetric due to the condition (4). It is satisfied: (TrC)2K4 det CO0 and the matrix C has two different positive eigenvalues ai and the corresponding unit eigenvectors gi. Denoting by ! g11 g12 TZ g21 g22 the orthogonal matrix and assuming that ^ n; uÞ U^ ðs; n; uÞ Z T Vðs;
(18)
t
sinðtÞ dt t
(21)
^ V~ mk ðsÞ Z KðV^ mk Þ V~ Z KðVÞ; ikm s Z Kikm am gm K 2½ciðkm sÞcosðkm sÞ k fipe
C siðkm sÞsinðkm sÞ g
(23)
~ the After applying inverse Radon transforms to U^ Z T V; fundamental solution and its derivatives are achieved: ! 1 ! ð g1 u~ 1 g12 u~ 1 g11 g12 1 U ðx;x; uÞ ¼ 2 dn 4p g21 g22 g21 u~ 2 g22 u~ 2 s¼jhn;xKxij jnj¼1
(24) here: u~ j ðs; uÞ Z Aj fipeikm s K 2½ciðkm sÞcosðkm sÞ 1 C siðkm sÞsinðkm sÞ g; Aj Z 2aj ð
0
g11 g12
(25)
1
1 @ A 4p2 2 2 g g 1 2 jnj¼1 0 1 1 1 g1 vs u~ 1 g2 vs u~ 1 A nk sgnhn;xKxidn !@ s¼jhn;xKzij 2 2 g1 vs u~ 2 g2 vs u~ 2
(17)
the following equation is obtained:
ðN
where qffiffiffiffiffi ikm jsj V^ mk Z am gm ; km Z arm u and am Z 2ami km :Using (20) ke it holds:
Uk ðx;x;uÞ ¼
here I2 is the 2!2 unit matrix. Multiplying from the left Eq. (17) with TK1 and since ! a1 0 K1 T CT Z A Z 0 a2
ðAv2s C ru2 I2 ÞV^ Z KdðsÞT K1
t
cosðtÞ dt; siðtÞ Z K t
are integral cosine and integral sine functions, respectively. Then: ! ! V^ 11 V^ 12 g11 g12 U^ Z T V^ Z (22) g21 g22 V^ 21 V^ 22
(16)
Eq. (14) in a matrix form becomes: ðCTv2s C ru2 TÞV^ Z KdðsÞI2 ;
ciðtÞ Z K
ðN
(26) and
2 vs u~ j ZAj Kpkj eikj s K C2kj ½ciðkj sÞsinðkj sÞ s Ksiðkj sÞcosðkj sÞ ;
(27)
Eq. (18) are four ordinary differential equations of the type: ðav2s C ru2 Þn^ Z KdðsÞg
(19)
The solution of the Eq. (19) is well known, see pffiffiffi ^ ageikjsj ; kZ ra u; aZ Vladimirov [44] and it is: nZ i 2ak : Using the calculus with generalized functions, see [43] it is obtained: ^ Z nðsÞ ~ Z ikagfipeiks K 2½ciðksÞcosðksÞ KðnÞ C siðksÞsinðksÞ g
(20)
Near-field asymptotic of Uij and Uij;k are: zdpjq Uij zbij ln r and Upj;q
1 at r/0 r
(28)
where bij and dpjq depend on the elastic constants and on the density, but not on the frequency. The asymptotic of the fundamental solution Uij as r/0 is presented by integrals of linear combinations of u˜j (s, u) over unit circle and its leading term is evaluated from ci(kjs)/O(ln r). The asymptotic of the derivatives
P. Dineva et al. / Engineering Analysis with Boundary Elements 29 (2005) 689–698 Uij;k as r/0 are presented by integrals of linear combinations of vsu˜j(s, u) over unit circle and their leading term is O(1/r). As an example we will show a fundamental solution and its derivatives for the time-harmonic elasto-isotropic case derived with the Radon transform ! ð n1 n2 1 U ðx; x; uÞ Z 2 4p Kn2 n1 jnjZ1 ! n1 u~1 ðs; uÞ Kn2 u~1 ðs; uÞ sZjhrKr ;nij dn ! 0 n2 u~2 ðs; uÞ n1 u~ 2 ðs; uÞ (29)
Note that only in the isotropic case qffiffiffiffia1, a2 do not depend on n and a1 Z lC 2m; a2 Z m; kj Z arj u:The derivative of the is: fundamental solution Uij;k ! ð n n 1 1 2 Uk ðx;x;uÞZ 2 4p Kn2 n1 jnjZ1 ! n1 vs u~ 1 ðs;uÞ Kn2 vs u~ 1 ðs;uÞ ! jsZjhrKr0 ;nij nk n2 vs u~ 2 ðs;uÞ n1 vs u~ 2 ðs;uÞ !sgnhrKr0 ;nidn ð30Þ Correspondingly: u~ j ðs;uÞzK
1 1 1 ln jsj and vs u~ j ðs;uÞzK 2 2 4p aj 4p aj s
(31)
and that gives the asymptotic of the fundamental solution Uij and its derivatives Uij;k : The difference between the fundamental solution derived through Fourier transform and composed by Bessel function UijF ; and the fundamental solution, derived through Radon transform UijR satisfies the Eq. (1) and it is different from zero. On the other side the asymptotic of UijF and UijR coincides. Due to the uniqueness of the solution of (1), (3) it follows that the solution of Eq. (10) with boundary condition (3) obtained by using by both fundamental solutions gives one and the same COD Dui over Scr and one and the same displacement and traction field in the exterior of the crack. The transient and frequency-dependent fundamental solution using Radon transforms is obtained in Wang and Achenbach [32], [33] for 2D and 3D cases, respectively and here we follow their approach. The obtained fundamental solution is based on the Radon transform properties and ellipticity of the differential operator in (11) due to the properties of the elastic tensor, see condition (4). The physical sense of the differential operator ellipticity corresponds to the requirement that the strain energy density EZ1/2Cijkl3ij3kl must remain positive since energy must be minimal in a state of stable equilibrium. The presented here fundamental solution is available for implementation in the non-hypersingular traction BIE (10). Such solution is applied recently in Rangelov [45]
693
for solving by BIEM direct and inverse scattering problems in anisotropic cracked plane.
4. Numerical solution procedure The numerical procedure for solution of the problem defined in Section 2 follows the procedure developed in [40,46,47] for isotropic material by the same authors. The non-singular traction BIE is collocated on one side of the crack boundary by using COD as unknowns. The approximations of the displacement and traction are with parabolic shape functions and they satisfy the following principles: Ho¨lder continuity at least at the collocation points and asymptotic pffiffi behaviour of the displacement near the crack tip as O r : QP-BE are implemented in a quadratic boundary element discretization. The disadvantages of the standard quadratic approximation concerning the smoothness at all irregular points as crack-tips and odd nodes, is overcome by the shifted point method, proposed by the authors in [40]. After the discretization, the obtained integrals are at least CPV integrals. The regular integrals are computed employing the Gaussian quadrature scheme for one-dimensional integrals and Monte Carlo integration scheme for twodimensional integrals. All singular integrals are solved analytically on the small neighbourhood of the field point, by using the approximation of the fundamental solution for small argument. After the discretization of the non-hyper-singular traction BIEs, overcoming of the week and strong singularities in the integrals and satisfaction of the boundary conditions, an algebraic complex system of equations according to the crack opening displacement is obtained and solved. The program codes based on Mathematica and FORTRAN are created following the above described procedure.
5. Validation of the proposed numerical procedure The above presented formulation is general and it allows to solve dynamic fracture problems for anisotropic materials. In order to validate the present approach, three test examples for a line crack in an isotropic, orthotropic and transversely-isotropic plane under incident longitudinal time-harmonic waves are analyzed. The results are compared with those obtained by other authors. Parametric study for ‘in-plane’ wave propagation and wave scattering under different incidence angles by a crack in anisotropic plane is presented. Let us mention that the existing results in the literature in Ohyoshi [7] and in Dhawan [9] are only for incidence angle qZp/2. The crack has a length 2a and it lies on Ox1 on the interval (Ka,Ca). The number of boundary elements on the crack is 5. The first boundary element is left QP-BE, the second, the third and the fourth elements are ordinary parabolic boundary elements and the fifth BE is right
694
P. Dineva et al. / Engineering Analysis with Boundary Elements 29 (2005) 689–698
QP-BE. The boundary elements lengths are correspondingly 0.15a, 0.56a, 0.58a. 0.56a, 0.15a. It is used the shifted points numerical scheme with xZ0.05 on the first and on the second BE, xZK0.05 on the fourth and fifth BE. SIF computation is obtained directly from the traction values ahead the crack-tip, see Saez et al. [25] pffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffi (32) KI Z lim t2 2pr and KII Z lim t1 2pr ; r/0
r/0
where ti is the i-th component of the traction at the point along the crack at the distance r out of the crack-tip. The SIF is normalized by its static value.
Table 1 Comparison of SIFs for L-wave scattering by a crack in an infinite isotropic plane obtained by BIEM based on two different elasto-dynamic fundamental solutions k1 a
Normalized SIF-I obtained by U*F
Incident angle p/2 0.1 1.04337 0.3 1.21845 0.5 1.32026 0.7 1.01096 1.0 0.49946 1.7 0.36855 3.0 0.17388
5.1. Example 1. Diffraction of incident longitudinal timeharmonic wave by a line crack in an infinite elasto-isotropic plane
k1 a
We consider a material with nZ0.25, mZ2.2165! 1010 N/m2, rZ2.4!103 kg/m3. The same test example is solved: analytically by Chen and Sih [48]; numerically with multi-domain displacement BIEM by Chirino and Dominguez [49]; numerically with non-hyper-singular traction BIEM by the authors in Rangelov et al. [40]. In both papers [40,49] it is used the frequency dependent fundamental solution ffiffiffiffiffiffiffi is pffiffiffiffiffiffibased on Fourier transform.pSIF normalized by mk22 pa; where k2Zu/C2 and C2 Z m=r is the shear wave velocity. Here the same example is solved by the usage of the fundamental solution (24)–(27) obtained by Radon transforms in isotropic case and numerical scheme given in Section 4. Table 1 shows SIFs obtained by both: (*) fundamental solution derived through Fourier transform and composed by modified Bessel functions U*F and (**) fundamental solution obtained by Radon transforms U*R. The presented results are for the incident angles p/2, p/4 and 0. Table 1 shows that the values of the normalized SIFs obtained by the nonhypersingular traction BIEM based on two different fundamental solutions coincide.
Incident angle p/4 0.1 0.69556 0.3 0.82063 0.5 0.93436 0.7 0.81340 1.0 0.61289 1.7 0.26908 3.0 0.34026
5.2. Example 2. Diffraction of incident longitudinal timeharmonic wave by a line crack in an infinite orthotropic plane Consider an infinite orthotropic plane Ox1x2 with the axis of material symmetry Ox2. The incident displacement uin i and the incident traction tiin Z sin ij nj components with incidence angle q and at a point xpZ(x1p, x2p) are: in u1 Z g11 ðqÞexp½Kik1 ðqÞðx1p cosðqÞ C x2p sinðqÞÞ (33) uin Z g2 ðqÞexp½Kik ðqÞðx cosðqÞ C x sinðqÞÞ 1 1p 2p 2 1 in sin (34) il Z Cilkm uk;m : pffiffiffiffiffiffiffiffiffiffiffiffiffiffi Here k1 ðqÞZ u r=a1 ðqÞ; a1 ðqÞ is the greater eigenvalue of the matrix (15), g1 ðqÞZ ðg11 ðqÞ; g21 ðqÞÞ is the corresponding unit eigenvector and n1 Z cosðqÞ; n2 Z sinðqÞ:
1.04323 1.21618 1.31335 1.00611 0.49825 0.36875 0.16565
Normalized SIF-I obtained by U*F
k1 a Incident angle 0 0.1 0.3 0.5 0.7 1.0 1.7 3.0
Normalized SIF-I obtained by U*R
Normalized SIF-I obtained by U*R
Normalized SIF-II obtained by U*F
Normalized SIF-II obtained by U*R
0.69547 0.81918 0.92973 0.80970 0.61062 0.26686 0.33596
0.34224 0.37095 0.40272 0.41131 0.37274 0.29225 0.33430
0.34219 0.37040 0.40101 0.40820 0.36899 0.28735 0.32709
Normalized SIF-I obtained by U*F
Normalized SIF-I obtained by U*R
0.34765 0.41057 0.47347 0.42703 0.35364 0.24071 0.05016
0.34760 0.40984 0.47115 0.42513 0.35251 0.23988 0.04935
Ohyoshi [7] solved the problem only using the incident field with angle qZp/2. For this case Eq. (33) is: in K1=2 x2 C tÞ (35) uin 1 Z 0; u2 Z exp Kiuððd22 Þ pffiffiffiffiffiffiffiffiffiffi where d22Zc 22/r and k1 Z u r=c22 Z uðd22 ÞK1=2 : The boundary conditions are js22 jx2Z0 Z Kiru d1=2 22
at jx1 j! a
js12 jx2Z0 Z 0
at jx1 jR a
js22 jx2Z0 Z 0
at jx1 jR a
(36)
The numerical results are obtained by the material characteristics given as a ratio of c11 following Ohyoshi [7]. Let us define cijZsijc, where cZ6.6495 GPa and rZ 2.4!103 kg/m3. Let us consider the following cases with respect to the ratio sij: † † † †
Case Case Case Case
1, 4, 6, 7,
s11Z1, s12Zs21Z1/3, s22Z1, s66Z1/3; s11Z1, s12Zs21Z1/3, s22Z1, s66Z1/30; s11Z1, s12Zs21Z1/3, s22Z1, s66Z1/6; s11Z1, s12Zs21Z1/30, s22Z1, s66Z1/3.
P. Dineva et al. / Engineering Analysis with Boundary Elements 29 (2005) 689–698 (a) 1.5
Ohyoshi results-curve 6 BIEM results-curve 6 Ohyoshi results-curve 7 BIEM results-curve 7
Normalized SIF -I
Normalized SIF-I
1.5
1
0.5
695
curve 1
curve 4
0.5
1 k1a
curve 1
curve 4
curve 6
curve 7
1
0.5
0
0
0
0
0.5
1
1.5 k1a
2
3
2.5
Fig. 2. Comparison between normalized SIF-I for a line crack in an infinite orthotropic plane for L-wave at incidence angle p/2 obtained by Ohyoshi [7] and the authors.
(b) 0.5
1.5
curve 6
2
curve 7
The numbering of the cases is the same as in the Ohyoshi paper [7]. The material in case 1 is isotropic. The orthotropic materials in case 4 and 6 have different coefficients only for c66. The material in case 7 is similar to the material in the case 1, but with different coefficient c12. The dynamic SIF is normalized by its static value ru pffiffiffiffiffiffi dK1=2 ap: The obtained BIEM results for normalized SIFs 22 in Figs. 3–5 are for incidence angle p/2, p/4 and 0, respectively and for different orthotropic elastic constants used in the curves 1, 4, 6 and 7. A comparison is done for the cases 6 and 7 with Ohyoshi results [7], who solved the problem by Fredholm integral equations of second kind. Fig. 2 shows that the present BIEM solution is very accurate. The curve 1 is for an isotropic solid. The effect of orthotropy at incident angle p/2 obtained by the dual integral equations in [7] and by the BIEM (Fig. 3) is that with decreasing the ratio c66/c11, the point of maximum curve 1
curve 4
curve 6
0.3
0.2
0.1
0
0
0.5
1
1.5
2
k1a
Fig. 4. (a,b) Normalized SIFs for a line crack in an infinite orthotropic plane for L-wave at incidence angle p/4.
1.5
curve 1
curve 4
curve 6
curve 7
curve 7
Normalized SIF-I
Normalized SIF-I
1.5
Normalized SIF SIF-II
0.4
1
0.5
1
0.5
0
0 0
0.5
1
1.5 k1a
2
2.5
3
Fig. 3. Normalized SIF-I for a line crack in an infinite orthotropic plane for L-wave at incidence angle p/2.
0
0.5
1 k1a
1.5
2
Fig. 5. Normalized SIF-I for a line crack in an infinite orthotropic plane for L-wave at incidence angle 08.
P. Dineva et al. / Engineering Analysis with Boundary Elements 29 (2005) 689–698
is shifted to lower frequency region, see curve 6 and curve 4. The effect of c12 is small, see curve 1 and curve 7. Figs. 4 and 5 show new results as far as the method used in [7] solves only the case of incident angle p/2. The conclusion made in [7] that the ratio c12/c11 does not influence on the diffraction is no more valid for an incidence angle different from p/2, see Figs. 4 and 5. It is observed that with decreasing of the incidence angle, the influence of the ratio c12/c11 increases. Fig. 4b presents the normalized SIF-II and it is seen that curves 4 and 6 with smaller value of the ratio c66/c11 have smaller values of SIF-II in comparison with curves 1 and 7. All Figs. 3–5 demonstrate the combined effect of the anisotropy and of the incident wave direction on the diffraction field.
BIEM authors reslts Dhawan results
1.6
Normalized SIF-I
696
1.2
0.8
0.4
0 0
0.5
1.5
1
2
2.5
3
k2a Fig. 7. Normalized SIF-I for a line crack in an infinite transversely-isotropic plane for L-wave at incidence angle p/2.
5.3. Example 3. Diffraction of incident longitudinal time-harmonic wave by a line crack in an infinite transversely-isotropic plane
(a)
curve 1
2
curve 2
curve 3
1.8
Normalized SIF-I
1.6 1.4 1.2 1 0.8 0.6 0.4 0.2 0 0
1
3
2
5
4
6
7
8
9
10
q = p N, N = 0, ...,10 10
(b) 0.8
curve 1
curve 2
curve 3
0.7 0.6 Normalized SIF-II
In the orthotropic case there exist only three principal directions that are mutually perpendicular. If isotropy prevails only in one plane, then the body is said to be transversely isotropic [48]. Consider the following three cases, see Fig. 6 for which the axis of material symmetry coincides: (a) with the coordinate axis Ox3 and the plane Ox1x2 is isotropic; (b) with the coordinate axis Ox1 and the plane Ox2x3 is isotropic; (c) with the coordinate axis Ox2 and the plane Ox1x3 is isotropic. Dhawan solved in [9] by the method of singular integral equation the problem for a line crack in an infinite transversely-isotropic plane under incident normal longitudinal time-harmonic wave. It is considered in [9] that the axis of material symmetry coincides with Ox2 axis and the crack lies on the Ox1 axis (case (c)). The material used in [9] is Cadmium crystal. The relation between stress and strains for the considered in [9] plane strain problem is represented by the elastic constants c11, c12, c22, c66 that can be described in terms of Young’s moduli, Poisson’s ratios and shear moduli. Fig. 7 shows that the BIEM results obtained by the authors and Dhawan’s results are very close. The good
0.5 0.4 0.3 0.2 0.1 0
0
1
2
3
4
5
6
7
8
9
10
q = p N, N = 0, ...,10 10
Fig. 6. The geometry of the problem in test example 3; x3 axis of symmetry: (b) x1 axis of symmetry; (c) x2 axis of symmetry.
Fig. 8. (a,b) Normalized SIFs versus incident angle at frequency k2aZ0.5 for a crack situated on axis Ox1 while the axis of material symmetry coincides with axis: (a) Ox3; (b) Ox1; (c) Ox2.
P. Dineva et al. / Engineering Analysis with Boundary Elements 29 (2005) 689–698
comparison confirms that the proposed method works accurate. Fig. 8 shows BIEM solution for normalized SIFs versus incidence angles qZ(p/10) N,NZ0,1,2,.,10 at k2aZ0.5. Curve 1 is for the case (a)—a crack is situated on the axis Ox1 and the axis of material symmetry coincides with the axis Ox3. It is clear that curve 1 corresponds to the elastoisotropic in-plane Ox1x2 deformation case. Curve 2 is for the case (b)—a crack is situated on the axis Ox1 and the axis of material symmetry coincides with the axis Ox1, while curve 3 is for the case (c)—a crack is situated on the axis Ox1 and the axis of material symmetry coincides with the axis Ox2. Fig. 8 demonstrate that the adequate evaluation of the wave-crack interaction field can be obtained after accounting for the direction of the incident wave and the mutual disposition of the crack-line and the axis of the material anisotropy.
6. Conclusions 2D elasto-dynamic analysis of an in-plane finite crack in an anisotropic infinite domain is presented by nonhypersingular traction BIEM in the frequency domain. The dynamic fundamental solution derived by Radon transform in the frequency domain is obtained in a closed form for 2D in-plane deformation state for general anisotropy. The proposed numerical scheme for SIF calculation is validated by the comparison with results achieved by other authors with different methods. Parametric study for the diffraction of waves under different incidence angles and for materials with a different type of anisotropy is presented for the cases in which the material anisotropy axis coincides or does not coincide with the crack boundary. The numerical results show that the wave-crack interaction in anisotropic continua depends on the combined influence of the incidence wave direction, the wave length and frequency and the mutual disposition of the crack-line and the material symmetry axis. The proposed fundamental solution, the presented numerical scheme and the developed program codes can be used as a base for a solution of dynamic anisotropic problems with a more complex geometry, mechanics and type of the dynamic load. Dynamic coupled problems for the piezoelectric materials can be studied with the proposed method as well.
Acknowledgements The authors acknowledge the support of the DFG under the grant number 436BUL 113/118/0-1.
697
References [1] Sih GC, Paris PC, Irwin GR. On cracks in rectilinear anisotropic bodies. Int J Fract 1965;1:189–203. [2] Sih GC, Liebowitz H, editors. Mathematical theory of brittle fracture. New York: Academic Press; 1968. [3] Lekhnitski SG. Theory of elasticity of an anisotropic body. San Francisco: Holden-Day; 1963. [4] Lekhnitski SG. Anisotropic plates. New York: Gordon and Breach; 1968. [5] Savin GN. Stress concentration around holes. New York: Pergamon Press; 1961. [6] Ohyoshi T. Effect of orthotropy on singular stresses produced near a crack tip by incident SH-waves. Z Angew Math Mech 1973;53: 409–11. [7] Ohyoshi T. Effect of orthotropy on singular stresses for a finite crack. ASME J Appl Mech 1973;40:491–7. [8] Dhawan GK. Interaction of elastic waves by a Griffith crack in an infinite transversely-isotropic medium. Int J Fract 1982;19:29–37. [9] Dhawan GK. Interaction of SV-waves by a Griffith crack in an infinite transversely-isotropic medium. Int J Fract 1983;20:103–10. [10] Karim MR, Kundu T. Transient surface response of layered isotropic and anisotropic half-spaces with interface-cracks: SH case. Int J Fract 1988;37:245–62. [11] Karim MR, Kundu T. Dynamic response of an orthotropic half-space with a subsurface crack: in-plane case. ASME J Appl Mech 1991;58: 988–95. [12] Kundu T, Bostrom A. Elastic wave scattering by a circular crack in a transversely isotropic solid. Wave Motion 1992;15:285–300. [13] Sarkar J, Mandal SC, Ghosh M. Diffraction of elastic waves by three coplanar Griffith cracks in an orthotropic medium. Int J Eng Sci 1995; 33(2):163–77. [14] Hua Z, Tian-You F, Lan-Quao T. Composite materials dynamic fracture studies by generalized Shmuely difference algorithm. Eng Fract Mech 1996;54:869–77. [15] Brebbia CA, Venturini WS, editors. Boundary element techniques: applications in stress analysis and heat transfer. Southampton: Computational Mechanics Publications; 1987. [16] Pageau SS, Biggers SB. Finite element evaluation of free-edge singular stress fields in anisotropic materials. Int J Numer Methods Eng 1995;38:2225–39. [17] Pageau SS, Joseph PF, Biggers SB. Finite element analysis of anisotropic materials with singular in-plane stress fields. Int J Solids Struct 1995;32:571–91. [18] Snyder MD, Cruse TA. Boundary integral equation analysis of cracked anisotropic plates. Int J Fract 1975;11:315–28. [19] Chan KS, Cruse TA. Stress intensity factors for anisotropic compacttension specimens with inclined cracks. Eng Fract Mech 1986;23: 863–74. [20] Zhao MH, Shen YP, Lin YJ, Liu GN. The method of analysis of cracks in three-dimensional transversely-isotropic media: boundary integral equation approach. Eng Anal BE 1988;21:169–78. [21] Tan CL, Gao YL. Boundary element analysis of plane anisotropic bodies with stress concentrations and cracks. Comput Struct 1992;20: 17–28. [22] Aliabadi MH, Brebbia CA, editors. Advances in BEM for fracture mechanics. Southampton: Computational Mechanics Publication; 1990. [23] Chandra A, Hu KX, Huang Y. A hybrid BEM formulation for multiple cracks in orthotropic elastic components. Comput Struct 1995;56: 785–97. [24] Pan E, Amadei B. Fracture mechanics analysis of cracked 2D anisotropic media with a new formulation of the boundary element method. Int J Fract 1996;77:161–74. [25] Saez A, Ariza MP, Dominguez J. Three-dimensional fracture analysis in transversely isotropic solids. Eng Anal BE 1997;20:287–98.
698
P. Dineva et al. / Engineering Analysis with Boundary Elements 29 (2005) 689–698
[26] Ang WT, Clements DL, Cooke T. A hyper-singular boundary integral equation for anti-plane crack problems for a class of inhomogeneous anisotropic elastic materials. Eng Anal BE 1999;23:567–72. [27] Zhang Ch. Transient elasto-dynamic anti-plane crack analysis of anisotropic solids. Int J Solids Struct 2000;37:6107–30. [28] Saez A, Dominguez J. Dynamic crack problems in three-dimensional transversely isotropic solids. Eng Anal BE 2001;25:203–10. [29] Saez A, Dominguez J. Far-field dynamic Green’s functions for BEM in transversely isotropic solids. Wave Motion 2000;32:113–23. [30] Saez A, Dominguez J. BEM analysis of wave scattering in transversely isotropic solids. Int J Numer Methods Eng 1999;44(9): 1283–300. [31] Abuquerque EL, Sollero P, Aliabadi MH. The BEM applied to time dependent problems in anisotropic materials. Int J Solids Struct 2002; 39:1405–22. [32] Wang CY, Achenbach JD. Elastodynamic fundamental solutions for anisotropic solids. Geophys Int J 1994;118:384–92. [33] Wang CY, Achenbach JD. Three-dimensional time-harmonic elastodynamic Green’s functions for anisotropic solids. Proc R Soc London A 1995;449:441–58. [34] Wang CY, Achenbach JD. Lamb’s problem for solids of general anisotropy. Wave Motion 1996;24:227–42. [35] Song Ch, Wolf J. Semi-analytical representation of stress singularities as occurring in cracks in anisotropic multimaterials with the scaled boundary finite-element method. Comput Struct 2002;80:183–97. [36] Cruse TA, Swedlow JL. Interactive program for analysis and design problems in advanced composites. Technical report 1971; CarnegieMellon University, Report AFLM-TR-71-268.
[37] Zhang Ch, Gross D. On wave propagation in elastic solids with cracks. Southampton: Computational Mechanics Publication; 1998. [38] Xu LY. Crack surface displacement coupling due to material anisotropy. Eng Fract Mech 1994;49(3):425–34. [39] Su RKL, Sun HY. Numerical solution of two-dimensional anisotropic crack problems. Int J Solids Struct 2003;40:4615–35. [40] Rangelov T, Dineva P, Gross D. A hyper-singular traction BIEM for stress intensity factor computation in a finite cracked body. Eng Anal BE 2003;27:9–21. [41] John F. Plane waves and spherical means applied to partial differential equations. New York: Wiley International Science; 1955. [42] Ludwig D. The Radon transform on Euclidean space. Commun Pure Appl Math 1966;29:49–81. [43] Gelfand M, Shilov G. Generalized functions. New York: Academic Press; 1964. [44] Vladimirov V. Equations of mathematical physics. Moscow: Mir; 1984. [45] Rangelov T. Scattering from cracks in elasto-anisotropic plane. Theor Appl Mech 2003;33:55–72. [46] Dineva P, Gross D, Rangelov T. Dynamic behaviour of a bi-material interface-cracked plate. Eng Fract Mech 2002;69:1193–218. [47] Dineva P, Gross D, Rangelov T. Dynamic behaviour of a cracked solder joint. J Sound Vibration 2002;256:81–102. [48] Chen EP, Sih GC. Scattering wave about stationary and mooving cracks. In: Sih GC, editor. Mechanics of fracture: elasto-dynamic crack problems. Leyden: Noordhoff; 1977. p. 119–212. [49] Chirino Fr, Dominguez J. Dynamic analysis of cracks using BEM. Eng Fract Mech 1989;34:1051–61.
BIEM for 2D steady-state problems in cracked anisotropic materials P. Dinevaa, T. Rangelovb,*, D. Grossc b
a Institute of Mechanics, Bulgarian Academy of Sciences, 1113 Sofia, Bulgaria Institute of Mathematics and Informatics, Bulgarian Academy of Sciences, 1113 Sofia, Bulgaria c Institute of Mechanics, Darmstadt University of Technology, 64289 Darmstadt, Germany
Received 28 October 2004; revised 17 February 2005; accepted 17 February 2005 Available online 21 April 2005
Abstract The two-dimensional ‘in-plane’ time-harmonic elasto-dynamic problem for anisotropic cracked solid is studied. The non-hypersingular traction boundary integral equation method (BIEM) is used in conjunction with closed form frequency dependent fundamental solution, obtained by Radon transform. Accuracy and convergence of the numerical solution for stress intensity factor (SIF) is studied by comparison with existing solutions in isotropic, transversely-isotropic and orthotropic cases. In addition a parametric study for the wave field sensitivity on wave, crack and anisotropic material parameters is presented. q 2005 Elsevier Ltd. All rights reserved. Keywords: Steady-state elastodynamics; Cracked anisotropic material; BIEM; SIF computation AMS Subject Classification: 74E10; 74S15; 74R10; 74H35
1. Introduction In many cases, engineering and construction materials such as modern fibre-reinforced composite, many crystals, piezoelectric materials, ice, wood, diamond, etc. possess both anisotropy and crack-like defects induced due to manufacturing or in-service conditions. Elastodynamic anisotropic crack analysis has attracted the attention of scientists in such fields as ultrasonic non-destructive evaluation for material quality inspection or seismology and exploration geophysics where observation of elastic waves generated by earthquakes or explosions may reveal the presence of structural discontinuities in the earth’s interior. Studies on crack problems in anisotropic materials have been done since many years. The initial groundwork was proposed by Sih et al. [1] and Sih and Liebowitz [2]. Lekhnitski’s [3,4] and Savin’s [5] systematic formulations of anisotropy are regarded as a foundation for the anisotropic elasticity. The methods used for the solution
* Corresponding author. Tel.: C359 2 979 2845; fax: C359 2 971 3649. E-mail address: [email protected] (T. Rangelov).
0955-7997/$ - see front matter q 2005 Elsevier Ltd. All rights reserved. doi:10.1016/j.enganabound.2005.02.006
of anisotropic crack problems are semi-analytical in Ohyoshi [6,7], Dhawan [8,9], Karim and Kundu [10,11], Kundu and Bostrom [12], Sarkar et al. [13]; numerical as finite-difference schemes in Hua et al. [14]; finite-element method in Brebbia and Venturini [15], Pageau and Biggers [16], Pageau et al. [17], BIEM in Snyder and Cruse [18], Chan and Cruse [19], Zhao et al. [20], Tan and Gao [21], Aliabadi and Brebbia [22], Chandra et al. [23], Pan and Amadei [24], Saez et al. [25], Ang et al. [26], Zhang [27], Saez and Dominguez [28–30], Abuquerque et al. [31], Wang and Achenbach [32–34] and hybrid methods in Song and Wolf [35]. In most of the semi-analytical papers various types of integral equation techniques are used. Fourier or Laplace transforms reduce the boundary-value problem to dual integral equations that are expressed in terms of Fredholm integral equations. Using this approach Ohyoshi [6,7] solved the anti-plane and the in-plane diffraction problems of waves in a cracked orthotropic plane. Dhawan [8,9] treated the same problem but for the transversely-isotropic case. Karin and Kundu [10,11] considered the dynamic response of a layered anisotropic half-space with anti-plane interfacecracks and an orthotropic half-space with a subsurface inplane crack. Kundu and Bostrom [12] studied wave scattering by a circular crack in a transversely-isotropic
690
P. Dineva et al. / Engineering Analysis with Boundary Elements 29 (2005) 689–698
solid and the dynamic response of three co-planar cracks in an infinite orthotropic medium in [13]. Although BIEM is a very well suited approach for the type of problems considered here, only a few works can be found for dynamic anisotropic crack problems. The main difficulty when applying BIEM to anisotropic dynamic crack problems is the availability of a fundamental solution that can be efficiently evaluated and implemented in a BEM code. The fundamental solution for 2D elastostatics obtained by Cruse and Swedlow [36] is applied by Abuquerque et al. [31] to transient anisotropic problems. Using multi-domain BIEM formulation, dual reciprocity method and traction singular quarter-point boundary elements (TSQP-BE) in [31] is determinated SIFs for a crack in a quasi-isotropic infinity long strip and for a rectangular orthotropic plate with a central crack and with a slanted edge crack. Wang and Achenbach [32–34] derived fundamental solutions by Radon transforms for general anisotropic solids. These fundamental solutions are related to transient and time-harmonic cases of 3D and 2D out of plane deformation states and they have been used also in Zhang [27] and in Saez and Dominguez [28–30]. The main aim of the present paper is the development of a non-hypersingular traction-based BIEM for solving 2D steady-state elastic in-plane wave scattering problems in plane anisotropic solids. The procedure is based on the non-hypersingular traction BIE proposed by Zhang and Gross [37] for the isotropic case, in conjunction with a closed form 2D fundamental solution, obtained by Radon transform and with crack-tip quarter-point boundary elements (QP-BE) implemented in the used parabolic approximation. The novelty of the present work lies in the generality of the proposed methodology and in the parametric study results revealing how the anisotropy of the material, the existence of cracks, the wave parameters as frequency and incident angle and mutual disposition of the crack line, axis of the material symmetry and wave propagation direction influence on the dynamic stress concentration fields. The advantages of the proposed method as being BIEM include all its important features as reduction of the dimension of the problem, implicit satisfaction of radiation condition at infinite, semi-analytical character of the method, high accuracy at solution of problems with stress gradients, in addition: (a) The method is valid for general type of anisotropy. It allows considering the cases when the coordinate axes do not coincide with the axis of material symmetry. (b) The proposed numerical procedure works for arbitrary incident wave angle, while in the most cases the used for the purpose methods as singular integral method has been applied only for the normal incident wave. (c) ‘In-plane’ crack is used in the proposed mechanical model, while in the most cases in the literature the simpler ‘anti-plane’ crack case is considered.
(d) The method has a potential to be applied for transient dynamic loads, after applying the Laplace or Fourier transformation. The paper is organized as follows. The setting of the boundary-value problem and steady-state non-hypersingular traction BIE formulation is given in Section 2. Fundamental solution in a closed form is presented in Section 3. Paragraph 4 focuses on the numerical solution procedure. Some results concerning the accuracy and the convergence of the proposed numerical scheme together with a parametric study on wave-crack interaction in anisotropic media are discussed in Section 5. The paper is finishing by some concluding remarks in Section 6.
2. Problem statement and traction BIM formulation We consider an infinite, homogeneous, anisotropic and linearly elastic solid containing a finite line crack G with a length 2a, see Fig. 1. The solid is subjected to an incident time-harmonic plane L-wave, and the deformation of the solid is as in a state of a plane strain. The non-zero quantities are the displacements components ui and the stresses components sij, i,jZ1,2. In cases of general anisotropy, the in-plane and the out-of-plane displacement are usually coupled. In contrast, we will assume, see Xu [38], that the material possesses a mid-plane symmetry such as the inplane (x1, x2 coordinate plane) and the out-of-plane deformations are decoupled. In the following we consider the in-plane problem, described by the equation of motion sij;j þ ru2 ui ¼ 0 in R2 nG;
(1)
Hooke’s law sij Z Cijkl uk;l
(2)
and the traction-free boundary condition on the crack-faces tj Z 0 on G:
Fig. 1. Incident L-wave in a cracked anisotropic plane.
(3)
P. Dineva et al. / Engineering Analysis with Boundary Elements 29 (2005) 689–698
Here r is the mass density, u is the frequency, tiZsijnj is the traction on the crack boundary GZGCgGK, ni are the components of the unit normal vector on GC and Cijkl are the elastic constants that satisfy the conditions
691
and the traction is obtained by using the Hooke’s law (2) ð tpsc ðx; uÞ ZKCpqkl nq sijk;l ðx; x; uÞDui ðx; uÞnj dGx ; x;GC GC
Cijkl Z Cjikl Z Cklij ; Cijkl gij gkl O 0 for every non zero real symmetric tensor gij :
(9) (4)
All subscripts vary 1, 2 and the conventional summation rules over double indices is implied. Finally, it is assumed that the displacement ui satisfies Sommerfeld type radiation condition at infinity. Substituting Eq. (2) into Eq. (1) the equation of motion can be written as c11 u1;11 C c66 u1;22 C 2c16 u1;12 C c16 u2;11 C c26 u2;22 Cðc12 C c66 Þu2;12 C ru2 u1 Z 0 c u C c u C ðc C c Þu C c u 16 1;11 26 1;22 12 66 1;21 66 2;11 Cc22 u2;22 C 2c26 u2;12 C ru2 u2 Z 0:
GC
(5)
sc ui ðx; uÞ Z uin i ðx; uÞ C ui ðx; uÞ
ð6Þ
where xZ(x1, x2) and u is the frequency. The functions in uin i ; ti are the displacement and traction components of the incident wave field in the absence of the crack, while sc functions usc i ; ti are the corresponding displacement and traction of the scattered field due to the interaction of the incident wave with the crack. It is assumed the incident wave to be known, while the scattered wave field is unknown and has to be determined. The scattered wave would satisfy the equation of motion (5) and the boundary condition (3), which can be rewritten as tisc Z Ktiin on G
(7)
The scattered displacement usc i ðx; uÞ can be represented as in Zhang and Gross [37] by the boundary integral ð ðx; uÞ Z K skij ðx; x; uÞDuk ðx; uÞni dGx ; x ;GC (8) usc j GC
Kru2 Cijkl nj ðxÞ KCijkl nj ðxÞ
Constants cmn (m,nZ1,2,6) denote the anisotropic material parameters in abbreviated notation, obtained from Cijkl following the rules (11)41, (22)42, (12Z21)46. There are altogether six independent constants. In the orthotropic case with the elasticity axes parallel to the coordinate axes c16Zc26Z0 is fulfilled and the independent material constants are four: c11, c12, c22, c66, see Su and Sun [39]. When the coordinate axis and the elasticity axis do not coincide there are 6 material constants. The interaction of an incident wave with the crack induces scattered waves. The total wave field can be presented as a sum of the incident wave field and the scattered wave field:
and ti ðx; uÞ Z tiin ðx; uÞ C tisc ðx; uÞ
Here sijk ðx; x; uÞ is the stress Green’s function and Dui Z ui jGC K ui jGK is the crack-opening displacement (COD). Following the approach in Zhang and Gross [37], for the isotropic case, the BIE for the anisotropic case is obtained by using the boundary condition (7) and the radiation condition at infinity: ð in ti ðx;uÞZCijkl nj ðxÞ sphk ðx;x;uÞDup;h ðx;uÞdll nl dGx
ð
ð
Udk ðx;x;uÞDud ðx;uÞdll nl dGx
GC
smlk ðx;x;uÞDum;l ðx;uÞnl dGx ; x2GC
GC
(10) Eq. (10) is integro-differential equation with respect to COD Dui and the singular integrals converge in CPV sense, if a priory smoothness requirements Dui2C1Ca (G) are fulfilled, see Rangelov et al. [40].
3. Fundamental solution, its derivatives and asymptotics Since Eq. (5) is with constant coefficients, there exists a fundamental solution U Z fUkj g, see John [41] that is a solution of the equation ðCij ðvÞ C dij ru2 ÞUjk ðx; x; uÞ Z Kdik dðx K xÞ
(11)
where Cij ðvÞZ Cimjn vxm vxn . Let’s note that in the isotropic case: C11 ðvÞZ ðlC 2mÞv2x1 C mv2x2 ; C12 ðvÞZ C21 ðvÞZ ðlC mÞvx1 vx2 ; C22 ðvÞZ ðlC 2mÞv2x2 C mv2x1 ; where l, m are the Lame´ constants. The Radon transform, see Ludwig [42] and Gelfand and Shilov [43] is of fundamental importance in order to work out the analytic solution of the problem at hand. Let f(x) be a function defined in R2 and s is a real number. The Radon transform R of function f(x) is defined as: ð Ð ^f ðs; nÞ Z R½f ðxÞ Z f ðxÞ dU Z f ðxÞdðs K hn; xiÞ dx; hn;xiZs
(12) where d is Dirac delta function, h,i means scalar product in R2. The Radon transform is an integration of f(x) over all planes defined by hn, xiZs. The inverse Radon transform in
692
P. Dineva et al. / Engineering Analysis with Boundary Elements 29 (2005) 689–698
2D can be written as: ð 1 Kðf^ðs; nÞÞj dn; f ðxÞ Z 2 4p
where: ðCN
vs f^ðs; nÞ ds: Kðf^Þ Z KN s K s
jnjZ1
(13) The next Radon transform properties will be used: R(d(x))Zd(s) and f^ðas; anÞZ a1 f^ðs; nÞ; Rða1 f1 C a2 f2 ÞZ a1 f^1 C a2 f^2 ; Rðvj f ðxÞÞZ nj vs f^ðs; nÞ: Let us apply the Radon transform to both side of Eq. (11). For simplicity xZ(0,0) is considered:
ðCij ðnÞv2s C ru2 dij ÞU^ jk ðs; n; uÞ Z Kdik dðsÞ
(14)
For example, in the orthotropic case the form of Cij(n) is: ! c11 n21 C c66 n22 ðc12 C c66 Þn1 n2 fCij ðnÞg Z (15) ðc12 C c66 Þn1 n2 c22 n22 C c66 n21 Eq. (14) is a linear system of ordinary differential equations. For solving it we transform it into the canonical form. We can consider n, jnjZ1 since only such n is used in the inverse Radon transform (13). The 2!2 matrix CZ{Cij(n)} is positive and symmetric due to the condition (4). It is satisfied: (TrC)2K4 det CO0 and the matrix C has two different positive eigenvalues ai and the corresponding unit eigenvectors gi. Denoting by ! g11 g12 TZ g21 g22 the orthogonal matrix and assuming that ^ n; uÞ U^ ðs; n; uÞ Z T Vðs;
(18)
t
sinðtÞ dt t
(21)
^ V~ mk ðsÞ Z KðV^ mk Þ V~ Z KðVÞ; ikm s Z Kikm am gm K 2½ciðkm sÞcosðkm sÞ k fipe
C siðkm sÞsinðkm sÞ g
(23)
~ the After applying inverse Radon transforms to U^ Z T V; fundamental solution and its derivatives are achieved: ! 1 ! ð g1 u~ 1 g12 u~ 1 g11 g12 1 U ðx;x; uÞ ¼ 2 dn 4p g21 g22 g21 u~ 2 g22 u~ 2 s¼jhn;xKxij jnj¼1
(24) here: u~ j ðs; uÞ Z Aj fipeikm s K 2½ciðkm sÞcosðkm sÞ 1 C siðkm sÞsinðkm sÞ g; Aj Z 2aj ð
0
g11 g12
(25)
1
1 @ A 4p2 2 2 g g 1 2 jnj¼1 0 1 1 1 g1 vs u~ 1 g2 vs u~ 1 A nk sgnhn;xKxidn !@ s¼jhn;xKzij 2 2 g1 vs u~ 2 g2 vs u~ 2
(17)
the following equation is obtained:
ðN
where qffiffiffiffiffi ikm jsj V^ mk Z am gm ; km Z arm u and am Z 2ami km :Using (20) ke it holds:
Uk ðx;x;uÞ ¼
here I2 is the 2!2 unit matrix. Multiplying from the left Eq. (17) with TK1 and since ! a1 0 K1 T CT Z A Z 0 a2
ðAv2s C ru2 I2 ÞV^ Z KdðsÞT K1
t
cosðtÞ dt; siðtÞ Z K t
are integral cosine and integral sine functions, respectively. Then: ! ! V^ 11 V^ 12 g11 g12 U^ Z T V^ Z (22) g21 g22 V^ 21 V^ 22
(16)
Eq. (14) in a matrix form becomes: ðCTv2s C ru2 TÞV^ Z KdðsÞI2 ;
ciðtÞ Z K
ðN
(26) and
2 vs u~ j ZAj Kpkj eikj s K C2kj ½ciðkj sÞsinðkj sÞ s Ksiðkj sÞcosðkj sÞ ;
(27)
Eq. (18) are four ordinary differential equations of the type: ðav2s C ru2 Þn^ Z KdðsÞg
(19)
The solution of the Eq. (19) is well known, see pffiffiffi ^ ageikjsj ; kZ ra u; aZ Vladimirov [44] and it is: nZ i 2ak : Using the calculus with generalized functions, see [43] it is obtained: ^ Z nðsÞ ~ Z ikagfipeiks K 2½ciðksÞcosðksÞ KðnÞ C siðksÞsinðksÞ g
(20)
Near-field asymptotic of Uij and Uij;k are: zdpjq Uij zbij ln r and Upj;q
1 at r/0 r
(28)
where bij and dpjq depend on the elastic constants and on the density, but not on the frequency. The asymptotic of the fundamental solution Uij as r/0 is presented by integrals of linear combinations of u˜j (s, u) over unit circle and its leading term is evaluated from ci(kjs)/O(ln r). The asymptotic of the derivatives
P. Dineva et al. / Engineering Analysis with Boundary Elements 29 (2005) 689–698 Uij;k as r/0 are presented by integrals of linear combinations of vsu˜j(s, u) over unit circle and their leading term is O(1/r). As an example we will show a fundamental solution and its derivatives for the time-harmonic elasto-isotropic case derived with the Radon transform ! ð n1 n2 1 U ðx; x; uÞ Z 2 4p Kn2 n1 jnjZ1 ! n1 u~1 ðs; uÞ Kn2 u~1 ðs; uÞ sZjhrKr ;nij dn ! 0 n2 u~2 ðs; uÞ n1 u~ 2 ðs; uÞ (29)
Note that only in the isotropic case qffiffiffiffia1, a2 do not depend on n and a1 Z lC 2m; a2 Z m; kj Z arj u:The derivative of the is: fundamental solution Uij;k ! ð n n 1 1 2 Uk ðx;x;uÞZ 2 4p Kn2 n1 jnjZ1 ! n1 vs u~ 1 ðs;uÞ Kn2 vs u~ 1 ðs;uÞ ! jsZjhrKr0 ;nij nk n2 vs u~ 2 ðs;uÞ n1 vs u~ 2 ðs;uÞ !sgnhrKr0 ;nidn ð30Þ Correspondingly: u~ j ðs;uÞzK
1 1 1 ln jsj and vs u~ j ðs;uÞzK 2 2 4p aj 4p aj s
(31)
and that gives the asymptotic of the fundamental solution Uij and its derivatives Uij;k : The difference between the fundamental solution derived through Fourier transform and composed by Bessel function UijF ; and the fundamental solution, derived through Radon transform UijR satisfies the Eq. (1) and it is different from zero. On the other side the asymptotic of UijF and UijR coincides. Due to the uniqueness of the solution of (1), (3) it follows that the solution of Eq. (10) with boundary condition (3) obtained by using by both fundamental solutions gives one and the same COD Dui over Scr and one and the same displacement and traction field in the exterior of the crack. The transient and frequency-dependent fundamental solution using Radon transforms is obtained in Wang and Achenbach [32], [33] for 2D and 3D cases, respectively and here we follow their approach. The obtained fundamental solution is based on the Radon transform properties and ellipticity of the differential operator in (11) due to the properties of the elastic tensor, see condition (4). The physical sense of the differential operator ellipticity corresponds to the requirement that the strain energy density EZ1/2Cijkl3ij3kl must remain positive since energy must be minimal in a state of stable equilibrium. The presented here fundamental solution is available for implementation in the non-hypersingular traction BIE (10). Such solution is applied recently in Rangelov [45]
693
for solving by BIEM direct and inverse scattering problems in anisotropic cracked plane.
4. Numerical solution procedure The numerical procedure for solution of the problem defined in Section 2 follows the procedure developed in [40,46,47] for isotropic material by the same authors. The non-singular traction BIE is collocated on one side of the crack boundary by using COD as unknowns. The approximations of the displacement and traction are with parabolic shape functions and they satisfy the following principles: Ho¨lder continuity at least at the collocation points and asymptotic pffiffi behaviour of the displacement near the crack tip as O r : QP-BE are implemented in a quadratic boundary element discretization. The disadvantages of the standard quadratic approximation concerning the smoothness at all irregular points as crack-tips and odd nodes, is overcome by the shifted point method, proposed by the authors in [40]. After the discretization, the obtained integrals are at least CPV integrals. The regular integrals are computed employing the Gaussian quadrature scheme for one-dimensional integrals and Monte Carlo integration scheme for twodimensional integrals. All singular integrals are solved analytically on the small neighbourhood of the field point, by using the approximation of the fundamental solution for small argument. After the discretization of the non-hyper-singular traction BIEs, overcoming of the week and strong singularities in the integrals and satisfaction of the boundary conditions, an algebraic complex system of equations according to the crack opening displacement is obtained and solved. The program codes based on Mathematica and FORTRAN are created following the above described procedure.
5. Validation of the proposed numerical procedure The above presented formulation is general and it allows to solve dynamic fracture problems for anisotropic materials. In order to validate the present approach, three test examples for a line crack in an isotropic, orthotropic and transversely-isotropic plane under incident longitudinal time-harmonic waves are analyzed. The results are compared with those obtained by other authors. Parametric study for ‘in-plane’ wave propagation and wave scattering under different incidence angles by a crack in anisotropic plane is presented. Let us mention that the existing results in the literature in Ohyoshi [7] and in Dhawan [9] are only for incidence angle qZp/2. The crack has a length 2a and it lies on Ox1 on the interval (Ka,Ca). The number of boundary elements on the crack is 5. The first boundary element is left QP-BE, the second, the third and the fourth elements are ordinary parabolic boundary elements and the fifth BE is right
694
P. Dineva et al. / Engineering Analysis with Boundary Elements 29 (2005) 689–698
QP-BE. The boundary elements lengths are correspondingly 0.15a, 0.56a, 0.58a. 0.56a, 0.15a. It is used the shifted points numerical scheme with xZ0.05 on the first and on the second BE, xZK0.05 on the fourth and fifth BE. SIF computation is obtained directly from the traction values ahead the crack-tip, see Saez et al. [25] pffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffi (32) KI Z lim t2 2pr and KII Z lim t1 2pr ; r/0
r/0
where ti is the i-th component of the traction at the point along the crack at the distance r out of the crack-tip. The SIF is normalized by its static value.
Table 1 Comparison of SIFs for L-wave scattering by a crack in an infinite isotropic plane obtained by BIEM based on two different elasto-dynamic fundamental solutions k1 a
Normalized SIF-I obtained by U*F
Incident angle p/2 0.1 1.04337 0.3 1.21845 0.5 1.32026 0.7 1.01096 1.0 0.49946 1.7 0.36855 3.0 0.17388
5.1. Example 1. Diffraction of incident longitudinal timeharmonic wave by a line crack in an infinite elasto-isotropic plane
k1 a
We consider a material with nZ0.25, mZ2.2165! 1010 N/m2, rZ2.4!103 kg/m3. The same test example is solved: analytically by Chen and Sih [48]; numerically with multi-domain displacement BIEM by Chirino and Dominguez [49]; numerically with non-hyper-singular traction BIEM by the authors in Rangelov et al. [40]. In both papers [40,49] it is used the frequency dependent fundamental solution ffiffiffiffiffiffiffi is pffiffiffiffiffiffibased on Fourier transform.pSIF normalized by mk22 pa; where k2Zu/C2 and C2 Z m=r is the shear wave velocity. Here the same example is solved by the usage of the fundamental solution (24)–(27) obtained by Radon transforms in isotropic case and numerical scheme given in Section 4. Table 1 shows SIFs obtained by both: (*) fundamental solution derived through Fourier transform and composed by modified Bessel functions U*F and (**) fundamental solution obtained by Radon transforms U*R. The presented results are for the incident angles p/2, p/4 and 0. Table 1 shows that the values of the normalized SIFs obtained by the nonhypersingular traction BIEM based on two different fundamental solutions coincide.
Incident angle p/4 0.1 0.69556 0.3 0.82063 0.5 0.93436 0.7 0.81340 1.0 0.61289 1.7 0.26908 3.0 0.34026
5.2. Example 2. Diffraction of incident longitudinal timeharmonic wave by a line crack in an infinite orthotropic plane Consider an infinite orthotropic plane Ox1x2 with the axis of material symmetry Ox2. The incident displacement uin i and the incident traction tiin Z sin ij nj components with incidence angle q and at a point xpZ(x1p, x2p) are: in u1 Z g11 ðqÞexp½Kik1 ðqÞðx1p cosðqÞ C x2p sinðqÞÞ (33) uin Z g2 ðqÞexp½Kik ðqÞðx cosðqÞ C x sinðqÞÞ 1 1p 2p 2 1 in sin (34) il Z Cilkm uk;m : pffiffiffiffiffiffiffiffiffiffiffiffiffiffi Here k1 ðqÞZ u r=a1 ðqÞ; a1 ðqÞ is the greater eigenvalue of the matrix (15), g1 ðqÞZ ðg11 ðqÞ; g21 ðqÞÞ is the corresponding unit eigenvector and n1 Z cosðqÞ; n2 Z sinðqÞ:
1.04323 1.21618 1.31335 1.00611 0.49825 0.36875 0.16565
Normalized SIF-I obtained by U*F
k1 a Incident angle 0 0.1 0.3 0.5 0.7 1.0 1.7 3.0
Normalized SIF-I obtained by U*R
Normalized SIF-I obtained by U*R
Normalized SIF-II obtained by U*F
Normalized SIF-II obtained by U*R
0.69547 0.81918 0.92973 0.80970 0.61062 0.26686 0.33596
0.34224 0.37095 0.40272 0.41131 0.37274 0.29225 0.33430
0.34219 0.37040 0.40101 0.40820 0.36899 0.28735 0.32709
Normalized SIF-I obtained by U*F
Normalized SIF-I obtained by U*R
0.34765 0.41057 0.47347 0.42703 0.35364 0.24071 0.05016
0.34760 0.40984 0.47115 0.42513 0.35251 0.23988 0.04935
Ohyoshi [7] solved the problem only using the incident field with angle qZp/2. For this case Eq. (33) is: in K1=2 x2 C tÞ (35) uin 1 Z 0; u2 Z exp Kiuððd22 Þ pffiffiffiffiffiffiffiffiffiffi where d22Zc 22/r and k1 Z u r=c22 Z uðd22 ÞK1=2 : The boundary conditions are js22 jx2Z0 Z Kiru d1=2 22
at jx1 j! a
js12 jx2Z0 Z 0
at jx1 jR a
js22 jx2Z0 Z 0
at jx1 jR a
(36)
The numerical results are obtained by the material characteristics given as a ratio of c11 following Ohyoshi [7]. Let us define cijZsijc, where cZ6.6495 GPa and rZ 2.4!103 kg/m3. Let us consider the following cases with respect to the ratio sij: † † † †
Case Case Case Case
1, 4, 6, 7,
s11Z1, s12Zs21Z1/3, s22Z1, s66Z1/3; s11Z1, s12Zs21Z1/3, s22Z1, s66Z1/30; s11Z1, s12Zs21Z1/3, s22Z1, s66Z1/6; s11Z1, s12Zs21Z1/30, s22Z1, s66Z1/3.
P. Dineva et al. / Engineering Analysis with Boundary Elements 29 (2005) 689–698 (a) 1.5
Ohyoshi results-curve 6 BIEM results-curve 6 Ohyoshi results-curve 7 BIEM results-curve 7
Normalized SIF -I
Normalized SIF-I
1.5
1
0.5
695
curve 1
curve 4
0.5
1 k1a
curve 1
curve 4
curve 6
curve 7
1
0.5
0
0
0
0
0.5
1
1.5 k1a
2
3
2.5
Fig. 2. Comparison between normalized SIF-I for a line crack in an infinite orthotropic plane for L-wave at incidence angle p/2 obtained by Ohyoshi [7] and the authors.
(b) 0.5
1.5
curve 6
2
curve 7
The numbering of the cases is the same as in the Ohyoshi paper [7]. The material in case 1 is isotropic. The orthotropic materials in case 4 and 6 have different coefficients only for c66. The material in case 7 is similar to the material in the case 1, but with different coefficient c12. The dynamic SIF is normalized by its static value ru pffiffiffiffiffiffi dK1=2 ap: The obtained BIEM results for normalized SIFs 22 in Figs. 3–5 are for incidence angle p/2, p/4 and 0, respectively and for different orthotropic elastic constants used in the curves 1, 4, 6 and 7. A comparison is done for the cases 6 and 7 with Ohyoshi results [7], who solved the problem by Fredholm integral equations of second kind. Fig. 2 shows that the present BIEM solution is very accurate. The curve 1 is for an isotropic solid. The effect of orthotropy at incident angle p/2 obtained by the dual integral equations in [7] and by the BIEM (Fig. 3) is that with decreasing the ratio c66/c11, the point of maximum curve 1
curve 4
curve 6
0.3
0.2
0.1
0
0
0.5
1
1.5
2
k1a
Fig. 4. (a,b) Normalized SIFs for a line crack in an infinite orthotropic plane for L-wave at incidence angle p/4.
1.5
curve 1
curve 4
curve 6
curve 7
curve 7
Normalized SIF-I
Normalized SIF-I
1.5
Normalized SIF SIF-II
0.4
1
0.5
1
0.5
0
0 0
0.5
1
1.5 k1a
2
2.5
3
Fig. 3. Normalized SIF-I for a line crack in an infinite orthotropic plane for L-wave at incidence angle p/2.
0
0.5
1 k1a
1.5
2
Fig. 5. Normalized SIF-I for a line crack in an infinite orthotropic plane for L-wave at incidence angle 08.
P. Dineva et al. / Engineering Analysis with Boundary Elements 29 (2005) 689–698
is shifted to lower frequency region, see curve 6 and curve 4. The effect of c12 is small, see curve 1 and curve 7. Figs. 4 and 5 show new results as far as the method used in [7] solves only the case of incident angle p/2. The conclusion made in [7] that the ratio c12/c11 does not influence on the diffraction is no more valid for an incidence angle different from p/2, see Figs. 4 and 5. It is observed that with decreasing of the incidence angle, the influence of the ratio c12/c11 increases. Fig. 4b presents the normalized SIF-II and it is seen that curves 4 and 6 with smaller value of the ratio c66/c11 have smaller values of SIF-II in comparison with curves 1 and 7. All Figs. 3–5 demonstrate the combined effect of the anisotropy and of the incident wave direction on the diffraction field.
BIEM authors reslts Dhawan results
1.6
Normalized SIF-I
696
1.2
0.8
0.4
0 0
0.5
1.5
1
2
2.5
3
k2a Fig. 7. Normalized SIF-I for a line crack in an infinite transversely-isotropic plane for L-wave at incidence angle p/2.
5.3. Example 3. Diffraction of incident longitudinal time-harmonic wave by a line crack in an infinite transversely-isotropic plane
(a)
curve 1
2
curve 2
curve 3
1.8
Normalized SIF-I
1.6 1.4 1.2 1 0.8 0.6 0.4 0.2 0 0
1
3
2
5
4
6
7
8
9
10
q = p N, N = 0, ...,10 10
(b) 0.8
curve 1
curve 2
curve 3
0.7 0.6 Normalized SIF-II
In the orthotropic case there exist only three principal directions that are mutually perpendicular. If isotropy prevails only in one plane, then the body is said to be transversely isotropic [48]. Consider the following three cases, see Fig. 6 for which the axis of material symmetry coincides: (a) with the coordinate axis Ox3 and the plane Ox1x2 is isotropic; (b) with the coordinate axis Ox1 and the plane Ox2x3 is isotropic; (c) with the coordinate axis Ox2 and the plane Ox1x3 is isotropic. Dhawan solved in [9] by the method of singular integral equation the problem for a line crack in an infinite transversely-isotropic plane under incident normal longitudinal time-harmonic wave. It is considered in [9] that the axis of material symmetry coincides with Ox2 axis and the crack lies on the Ox1 axis (case (c)). The material used in [9] is Cadmium crystal. The relation between stress and strains for the considered in [9] plane strain problem is represented by the elastic constants c11, c12, c22, c66 that can be described in terms of Young’s moduli, Poisson’s ratios and shear moduli. Fig. 7 shows that the BIEM results obtained by the authors and Dhawan’s results are very close. The good
0.5 0.4 0.3 0.2 0.1 0
0
1
2
3
4
5
6
7
8
9
10
q = p N, N = 0, ...,10 10
Fig. 6. The geometry of the problem in test example 3; x3 axis of symmetry: (b) x1 axis of symmetry; (c) x2 axis of symmetry.
Fig. 8. (a,b) Normalized SIFs versus incident angle at frequency k2aZ0.5 for a crack situated on axis Ox1 while the axis of material symmetry coincides with axis: (a) Ox3; (b) Ox1; (c) Ox2.
P. Dineva et al. / Engineering Analysis with Boundary Elements 29 (2005) 689–698
comparison confirms that the proposed method works accurate. Fig. 8 shows BIEM solution for normalized SIFs versus incidence angles qZ(p/10) N,NZ0,1,2,.,10 at k2aZ0.5. Curve 1 is for the case (a)—a crack is situated on the axis Ox1 and the axis of material symmetry coincides with the axis Ox3. It is clear that curve 1 corresponds to the elastoisotropic in-plane Ox1x2 deformation case. Curve 2 is for the case (b)—a crack is situated on the axis Ox1 and the axis of material symmetry coincides with the axis Ox1, while curve 3 is for the case (c)—a crack is situated on the axis Ox1 and the axis of material symmetry coincides with the axis Ox2. Fig. 8 demonstrate that the adequate evaluation of the wave-crack interaction field can be obtained after accounting for the direction of the incident wave and the mutual disposition of the crack-line and the axis of the material anisotropy.
6. Conclusions 2D elasto-dynamic analysis of an in-plane finite crack in an anisotropic infinite domain is presented by nonhypersingular traction BIEM in the frequency domain. The dynamic fundamental solution derived by Radon transform in the frequency domain is obtained in a closed form for 2D in-plane deformation state for general anisotropy. The proposed numerical scheme for SIF calculation is validated by the comparison with results achieved by other authors with different methods. Parametric study for the diffraction of waves under different incidence angles and for materials with a different type of anisotropy is presented for the cases in which the material anisotropy axis coincides or does not coincide with the crack boundary. The numerical results show that the wave-crack interaction in anisotropic continua depends on the combined influence of the incidence wave direction, the wave length and frequency and the mutual disposition of the crack-line and the material symmetry axis. The proposed fundamental solution, the presented numerical scheme and the developed program codes can be used as a base for a solution of dynamic anisotropic problems with a more complex geometry, mechanics and type of the dynamic load. Dynamic coupled problems for the piezoelectric materials can be studied with the proposed method as well.
Acknowledgements The authors acknowledge the support of the DFG under the grant number 436BUL 113/118/0-1.
697
References [1] Sih GC, Paris PC, Irwin GR. On cracks in rectilinear anisotropic bodies. Int J Fract 1965;1:189–203. [2] Sih GC, Liebowitz H, editors. Mathematical theory of brittle fracture. New York: Academic Press; 1968. [3] Lekhnitski SG. Theory of elasticity of an anisotropic body. San Francisco: Holden-Day; 1963. [4] Lekhnitski SG. Anisotropic plates. New York: Gordon and Breach; 1968. [5] Savin GN. Stress concentration around holes. New York: Pergamon Press; 1961. [6] Ohyoshi T. Effect of orthotropy on singular stresses produced near a crack tip by incident SH-waves. Z Angew Math Mech 1973;53: 409–11. [7] Ohyoshi T. Effect of orthotropy on singular stresses for a finite crack. ASME J Appl Mech 1973;40:491–7. [8] Dhawan GK. Interaction of elastic waves by a Griffith crack in an infinite transversely-isotropic medium. Int J Fract 1982;19:29–37. [9] Dhawan GK. Interaction of SV-waves by a Griffith crack in an infinite transversely-isotropic medium. Int J Fract 1983;20:103–10. [10] Karim MR, Kundu T. Transient surface response of layered isotropic and anisotropic half-spaces with interface-cracks: SH case. Int J Fract 1988;37:245–62. [11] Karim MR, Kundu T. Dynamic response of an orthotropic half-space with a subsurface crack: in-plane case. ASME J Appl Mech 1991;58: 988–95. [12] Kundu T, Bostrom A. Elastic wave scattering by a circular crack in a transversely isotropic solid. Wave Motion 1992;15:285–300. [13] Sarkar J, Mandal SC, Ghosh M. Diffraction of elastic waves by three coplanar Griffith cracks in an orthotropic medium. Int J Eng Sci 1995; 33(2):163–77. [14] Hua Z, Tian-You F, Lan-Quao T. Composite materials dynamic fracture studies by generalized Shmuely difference algorithm. Eng Fract Mech 1996;54:869–77. [15] Brebbia CA, Venturini WS, editors. Boundary element techniques: applications in stress analysis and heat transfer. Southampton: Computational Mechanics Publications; 1987. [16] Pageau SS, Biggers SB. Finite element evaluation of free-edge singular stress fields in anisotropic materials. Int J Numer Methods Eng 1995;38:2225–39. [17] Pageau SS, Joseph PF, Biggers SB. Finite element analysis of anisotropic materials with singular in-plane stress fields. Int J Solids Struct 1995;32:571–91. [18] Snyder MD, Cruse TA. Boundary integral equation analysis of cracked anisotropic plates. Int J Fract 1975;11:315–28. [19] Chan KS, Cruse TA. Stress intensity factors for anisotropic compacttension specimens with inclined cracks. Eng Fract Mech 1986;23: 863–74. [20] Zhao MH, Shen YP, Lin YJ, Liu GN. The method of analysis of cracks in three-dimensional transversely-isotropic media: boundary integral equation approach. Eng Anal BE 1988;21:169–78. [21] Tan CL, Gao YL. Boundary element analysis of plane anisotropic bodies with stress concentrations and cracks. Comput Struct 1992;20: 17–28. [22] Aliabadi MH, Brebbia CA, editors. Advances in BEM for fracture mechanics. Southampton: Computational Mechanics Publication; 1990. [23] Chandra A, Hu KX, Huang Y. A hybrid BEM formulation for multiple cracks in orthotropic elastic components. Comput Struct 1995;56: 785–97. [24] Pan E, Amadei B. Fracture mechanics analysis of cracked 2D anisotropic media with a new formulation of the boundary element method. Int J Fract 1996;77:161–74. [25] Saez A, Ariza MP, Dominguez J. Three-dimensional fracture analysis in transversely isotropic solids. Eng Anal BE 1997;20:287–98.
698
P. Dineva et al. / Engineering Analysis with Boundary Elements 29 (2005) 689–698
[26] Ang WT, Clements DL, Cooke T. A hyper-singular boundary integral equation for anti-plane crack problems for a class of inhomogeneous anisotropic elastic materials. Eng Anal BE 1999;23:567–72. [27] Zhang Ch. Transient elasto-dynamic anti-plane crack analysis of anisotropic solids. Int J Solids Struct 2000;37:6107–30. [28] Saez A, Dominguez J. Dynamic crack problems in three-dimensional transversely isotropic solids. Eng Anal BE 2001;25:203–10. [29] Saez A, Dominguez J. Far-field dynamic Green’s functions for BEM in transversely isotropic solids. Wave Motion 2000;32:113–23. [30] Saez A, Dominguez J. BEM analysis of wave scattering in transversely isotropic solids. Int J Numer Methods Eng 1999;44(9): 1283–300. [31] Abuquerque EL, Sollero P, Aliabadi MH. The BEM applied to time dependent problems in anisotropic materials. Int J Solids Struct 2002; 39:1405–22. [32] Wang CY, Achenbach JD. Elastodynamic fundamental solutions for anisotropic solids. Geophys Int J 1994;118:384–92. [33] Wang CY, Achenbach JD. Three-dimensional time-harmonic elastodynamic Green’s functions for anisotropic solids. Proc R Soc London A 1995;449:441–58. [34] Wang CY, Achenbach JD. Lamb’s problem for solids of general anisotropy. Wave Motion 1996;24:227–42. [35] Song Ch, Wolf J. Semi-analytical representation of stress singularities as occurring in cracks in anisotropic multimaterials with the scaled boundary finite-element method. Comput Struct 2002;80:183–97. [36] Cruse TA, Swedlow JL. Interactive program for analysis and design problems in advanced composites. Technical report 1971; CarnegieMellon University, Report AFLM-TR-71-268.
[37] Zhang Ch, Gross D. On wave propagation in elastic solids with cracks. Southampton: Computational Mechanics Publication; 1998. [38] Xu LY. Crack surface displacement coupling due to material anisotropy. Eng Fract Mech 1994;49(3):425–34. [39] Su RKL, Sun HY. Numerical solution of two-dimensional anisotropic crack problems. Int J Solids Struct 2003;40:4615–35. [40] Rangelov T, Dineva P, Gross D. A hyper-singular traction BIEM for stress intensity factor computation in a finite cracked body. Eng Anal BE 2003;27:9–21. [41] John F. Plane waves and spherical means applied to partial differential equations. New York: Wiley International Science; 1955. [42] Ludwig D. The Radon transform on Euclidean space. Commun Pure Appl Math 1966;29:49–81. [43] Gelfand M, Shilov G. Generalized functions. New York: Academic Press; 1964. [44] Vladimirov V. Equations of mathematical physics. Moscow: Mir; 1984. [45] Rangelov T. Scattering from cracks in elasto-anisotropic plane. Theor Appl Mech 2003;33:55–72. [46] Dineva P, Gross D, Rangelov T. Dynamic behaviour of a bi-material interface-cracked plate. Eng Fract Mech 2002;69:1193–218. [47] Dineva P, Gross D, Rangelov T. Dynamic behaviour of a cracked solder joint. J Sound Vibration 2002;256:81–102. [48] Chen EP, Sih GC. Scattering wave about stationary and mooving cracks. In: Sih GC, editor. Mechanics of fracture: elasto-dynamic crack problems. Leyden: Noordhoff; 1977. p. 119–212. [49] Chirino Fr, Dominguez J. Dynamic analysis of cracks using BEM. Eng Fract Mech 1989;34:1051–61.