3D BEM for general anisotropic elasticity

International Journal of Solids and Structures 44 (2007) 7073–7091 www.elsevier.com/locate/ijsolstr

3D BEM for general anisotropic elasticity C.-Y. Wang a, M. Denda a b

b,*

Mathematics and Modelling Department, Schlumberger-Doll Research, Old Quarry RD, Ridgefield, CT 06877-4108, USA Mechanical and Aerospace Engineering Department, Rutgers University, 98 Brett Road, Piscataway, NJ 08854-8058, USA Received 3 November 2006; received in revised form 21 March 2007 Available online 31 March 2007

Abstract A new method for computing the system matrices G and H of the three-dimensional boundary element method is presented for the elastostatic problems of general anisotropy. Green’s functions, expressed by the line integrals over a semicircle, have a simple structure that allows analytical integration over the boundary elements. The computation for the system matrices G and H is reduced to numerical integration of regular line integrals over the unit semi-circle. The proposed method has been implemented and tested with numerical examples.  2007 Elsevier Ltd. All rights reserved. Keywords: 3D BEM; General anisotropic solids; Line integral representation of Green’s functions

1. Introduction The boundary element method (BEM) is one of the most widely used numerical methods along with the finite element method (FEM) and finite difference method (FDM). To a large extend, however, applications of the BEM in solving elasticity problems is limited to isotropic materials. This is due to the lack of Green’s functions needed for the BEM implementation. Among the six fundamental solutions, namely, static, timeharmonic and time-transient solutions in 2D and 3D, the 2D static Green’s functions are the simplest (Eshelby et al., 1953; Stroh, 1958; Ting, 1996; Wang, 1994). The 3D static Green’s functions (Wang, 1997; Sales and Gray, 1998; Tonon et al., 2001), and the 2D time-transient Green’s functions (Wang and Achenbach, 1992) can be derived, albeit complicated, explicitly by use of Cauchy’s residue theorem. The 2D and 3D time-harmonic and 3D time-transient Green’s functions, which cannot be expressed explicitly in general, are given in terms of integrals over a unit circle or a unit sphere (Wang and Achenbach, 1994). The BEM has been applied to the 2D static problems of general anisotropy successfully (Berger and Tewary, 1997; Clements and Haselgrove, 1983; Denda and Marante, 2004; Snyder and Cruse, 1975; Sollero and Aliabadi, 1993), and to the 3D static problems of transversely isotropy where relatively simple explicit solutions of Green’s functions are available (Pan and Chou, 1976). For the 3D problems of general anisotropy, the BEM has been developed *

Corresponding author. Tel.: +1 732 445 4391; fax: +1 732 445 3124. E-mail address: [email protected] (M. Denda).

0020-7683/$ - see front matter  2007 Elsevier Ltd. All rights reserved. doi:10.1016/j.ijsolstr.2007.03.026

7074

C.-Y. Wang, M. Denda / International Journal of Solids and Structures 44 (2007) 7073–7091

using either the integral expressions of Green’s functions (Barnett, 1972; Gray et al., 2003; Synge, 1957; Vogel and Rizzo, 1973; Wilson and Cruse, 1978) or the explicit expressions (Chen and Lin, 1995; Ting and Lee, 1997; Nakamura and Tanuma, 1997; Pan and Yuan, 2000; Sales and Gray, 1998; Tonon et al., 2001; Lee, 2003). The latter involves tedious calculations of the derivatives of the displacement Green’s functions. Most of the work mentioned above followed the conventional BEM scheme for the numerical computation of the integrations over the boundary elements by calculating Green’s functions at each integration point. Up until now, the 3D BEMs for the anisotropic materials, in comparison to their isotropic counterparts, have been slow in computation and difficult in implementation due to the complexity of Green’s functions. Progress has been made in the 2D dynamic analysis using the integral expressions of Green’s functions derived by the Radon transform (Wang and Achenbach, 1994). These solutions have a simple structure that enables the analytical evaluation of the integrals over the boundary elements. This approach has been successfully applied to solve the 2D problems of the time-transient wave scattering (Wang et al., 1996), time-transient analysis of crack opening and stress concentration (Hirose et al., 2000), and the time-harmonic eigenvalue analysis in elastic (Denda et al., 2003) as well as in piezoelectric solids (Denda et al., 2004). In performing the analytical integrations over the boundary elements, straight(in 2D)/flat(in 3D) elements have been used. Since a large curved element can always be substituted by several small straight segments, the boundaries of most practical problems can be approximated by straight elements quite well. The extension of this approach to the 3D dynamic problems requires some additional manipulation of Green’s functions. As has been shown in Wang and Achenbach (1994), the time-transient and time-harmonic Green’s functions can split to a singular static part and a regular dynamic part. Since the treatment of the regular part is straightforward, it comes down to the problem of how to treat the 3D static problem, which is the goal of this paper. We exploit the unique structure of Green’s functions and subsequent analytical integration over the boundary elements to reduce the computation for the system matrices G and H to the regular line integrations over the unit semi-circle. Furthermore, the integrands of the line integrals have a simple structure that facilitates the coding procedure and computation for very large systems. The proposed method has been implemented and tested with numerical examples. While the proposed approach shares some features with the techniques used by Gray et al. (2003), and with Stokes’ theorem techniques by Li et al. (1998), it is a distinctly different algorithm. The paper is organized as follows. Section 2 presents the boundary integral equations, definition of Green’s functions, and the standard BEM discretization procedure. Section 3 follows the framework of Wang and Achenbach (1995) to derive Green’s functions in terms of a line integral over a unit semi-circle, using an integral representation of the delta function and a subsequent application of the Cauchy residue theorem. The integrand is decoupled into a product of two terms: one is a function of the material properties only and the other is a simple algebraic function that depends on the source and observation points. In Section 4, the unique decoupling properties of Green’s functions are utilized to analytically perform the integrations over the flat triangular boundary elements with the piecewise linear shape functions. This analytical evaluation of the surface integrals is the key results of this paper that enables a saving in computation time. The final expressions for the system matrices G and H are given by the regular line integrals of these analytical terms over the unit semi-circle. In Sections 5 and 6, the proposed method is implemented and demonstrated through numerical examples. The numerical accuracy is tested with a uniformly loaded cube whose analytical solution is known. Then, the boundary value problem of a spherical cavity is solved, first for an isotropic solid, where the results are compared against the analytical solution and the solution by the conventional BEM. Numerical solutions are produced for a number of anisotropic materials, including aluminum polycrystal, aluminum crystal, barium sodium niobate, and sapphire. Section 7 gives concluding remarks. 2. BEM formulation In a fixed Cartesian coordinate system, relationship between the displacement vector u and traction vector t on the boundary oD of a three dimensional linearly elastic body D is provided by the boundary integral equation Z Z cjp ðyÞuj ðyÞ ¼ gjp ðx; yÞtj ðxÞdx  hjp ðx; yÞuj ðxÞdx; ð1Þ oD

oD

C.-Y. Wang, M. Denda / International Journal of Solids and Structures 44 (2007) 7073–7091

7075

where cjp are the free term constants, and x and y are position vectors usually referred to as the observation and source points of Green’s functions gjp(x, y) and hjp(x, y). In this paper, a Roman subscript takes values of 1, 2, and 3. Summation convention is applied over the range of the suffixes. For a homogeneous linearly elastic solid of general anisotropy, the displacement Green’s functions gjp are defined by the partial differential equation Cjp ðox Þgpk ðx; yÞ ¼ djk dðx  yÞ;

ð2Þ

where djk is the Kronecker delta, d(.) is the delta function, and Cjp ðox Þ ¼ cijpq oi oq :

ð3Þ

The elastic stiffness tensor cijpq is fully symmetric and positive definite, i.e, cijpq = cjipq = cijqp = cpqij and cijpqeijepq > 0 for all non-zero real symmetric tensors eij. The traction Green’s functions hjp are determined by hjp ðx; yÞ ¼ ni ðxÞcijlq glp;q ðx; yÞ;

ð4Þ

where n is the unit outward normal to the boundary oD. Here and in the sequel, derivatives with respect to xi are denoted by oi or ,i. The basic idea of the BEM is to approximate the integral Eq. (1) with a finite system of linear equations through a discretization process. The first step of the discretization is to replace the boundary surface oD by a collection of small boundary elements S(m) oD 

M X

S ðmÞ

ð5Þ

m¼1

with M being the total number of elements. The second step is to approximate the boundary variables in each element by the interpolations between nodes. For example, for x 2 S(m) we may write uj ðxÞ 

  ðnÞ ðnÞ ðmnÞ uðnÞ n1 ; n2 uj ;

N X

tj ðxÞ 

n¼1

N X

  ðnÞ ðnÞ ðmnÞ uðnÞ n1 ; n2 tj ;

x 2 S ðmÞ ;

ð6Þ

n¼1 ðmnÞ

ðmnÞ

where u(n) are the interpolation functions (or the shape functions), uj and tj denote the displacement and ðnÞ ðnÞ traction vectors at the nth node in the mth element. In (6), (n1 , n2 ) are the local coordinates which will be specified later. Substituting (6) into (1), we obtain cjp up ðyÞ ¼

M X N h X

ðmnÞ

ðmnÞ

gjp ðyÞtj

ðmnÞ

ðmnÞ

 hjp ðyÞuj

i ;

ð7Þ

m¼1 n¼1

where ðmnÞ gjp ðyÞ ðmnÞ

¼

hjp ðyÞ ¼

Z ZS

ðmÞ

S ðmÞ

      ðnÞ ðnÞ ðnÞ ðnÞ ðnÞ ðnÞ gjp x n1 ; n2  y uðnÞ n1 ; n2 dn1 dn2 ;

ð8Þ

      ðnÞ ðnÞ ðnÞ ðnÞ ðnÞ ðnÞ hjp x n1 ; n2  y uðnÞ n1 ; n2 dn1 dn2 :

ð9Þ

Locating y to every nodes of all elements, (7) provides a system of regular linear algebraic equations to comðmnÞ ðmnÞ ðmnÞ ðmnÞ pute the boundary values uj and/or tj . Obviously, calculations for gjp and hjp are the core part of the BEM. In the next two sections, we shall first derive integral expressions for Green’s functions gjp and hjp, and then calculate analytical solutions for integrals (8) and (9). 3. Green’s functions In this section we follow the framework of Wang and Achenbach (1995) to derive a line integral expression for Green’s functions defined by (2). Starting with the integral representation of d(x)

7076

C.-Y. Wang, M. Denda / International Journal of Solids and Structures 44 (2007) 7073–7091

1 D 8p2

Z

d½d  ðx  yÞdd ¼ dðx  yÞ;

ð10Þ

jdj¼1

where D is the Laplace operator in 3D. It is easily verified that (2) is satisfied by Z 1 Ajp ðdÞ d½d  ðx  yÞdd; gjp ðx; yÞ ¼ 2 8p jdj¼1 DðdÞ

ð11Þ

where Ajp ðdÞ ¼ adj Cjp ðdÞ;

DðdÞ ¼ det Cjp ðdÞ;

ð12Þ

are the adjoint matrix (transpose of the cofactor matrix) and the determinant of matrix Cjp(d), respectively. Let (e1, e2, n) be a set of unit orthogonal base vectors (recall n is the outward normal to the boundary oD). We express the unit sphere |d| = 1 by d ¼ sin / cos he1 þ sin / sin he2 þ cos /n

ð13Þ

and change integral (11) to Z pZ p 1 Ajp ðdÞ d½d  ðx  yÞ sin / d/ dh: gjp ðx  yÞ ¼ 2 4p 0 0 DðdÞ

ð14Þ

In deriving Eq. (14), symmetry with respect to h has been used. A subsequent substitution of g = ctan(/) leads to Z pZ 1 1 Ajp ðsÞ d½s  ðx  yÞdg dh; ð15Þ gjp ðx  yÞ ¼ 2 4p 0 1 DðsÞ where s ¼ cos he1 þ sin he2 þ gn:

ð16Þ

Following Gel’fand et al. (1966)   1 1 ¼ }  ipsignðcÞdðxÞ; x þ ic x

for c ! 0;

ð17Þ

where } indicates the principal part,   signðcÞ 1 I dðxÞ ¼ ; for c ! 0: p x þ ic

ð18Þ

Hence, we have gjp ðx; yÞ ¼

signðcÞ I 4p3

Z

p 0

Z

1

1

Ajp ðsÞ dg dh ; DðsÞ s  ðx  yÞ þ ic

for c ! 0:

ð19Þ

For a fixed h, D(s(g)) is a polynomial of g of order six with real coefficients. Furthermore, the positive definiteness of the elastic constants cijpq guarantees that the system of partial differential equations (2) is strictly elliptical, or equivalently, the polynomial D(s(g)) does not have real roots. Therefore we may write DðsÞ ¼ d 6 g6 þ d 5 g5 þ    þ d 0 ¼ d 6

3 Y

ðg  gk Þðg  gk Þ;

ð20Þ

k¼1

where di(h) are the real valued coefficients, gk(h) are roots of D = 0 with positive imaginary parts (i.e, Iðgk Þ > 0Þ, and  gk are complex conjugates of gk. Based on Cauchy’s residue theorem, integration of (19) with respect to g yields residues at g = gk. The result is Z pX 3 sign½n  ðx  yÞ dh; ð21Þ gjp ðx; yÞ ¼ R gkjp ðhÞ sk  ðx  yÞ 0 k¼1

C.-Y. Wang, M. Denda / International Journal of Solids and Structures 44 (2007) 7073–7091

7077

where sk ¼ cos he1 þ sin he2 þ gk n

ð22Þ

and gkjp ðhÞ ¼

1 Ajp ðsk Þ : 2p2 D;g ðsk Þ

ð23Þ

In the above expression,

" # 6 X o ðm1Þ D;g ðs Þ  DðsÞðg¼gk Þ ¼ md m g og m¼1 k

¼ d 6 ðgk  gk Þ

3 Y

ðg  gq Þðg  gq Þ:

ð24Þ

ðq¼1Þðq6¼kÞ

ðg¼gk Þ

Following the procedure presented above, we can also obtain the traction Green’s functions Z pX 3 sign½n  ðx  yÞ  hjp ðx; yÞ ¼ R dh; hkjp ðhÞ 2 ½sk  ðx  yÞ 0 k¼1

ð25Þ

where  hkjp ðhÞ ¼ ni cijlq skq gklp ðhÞ:

ð26Þ

We notice that the integrand of the line integral (21) is given by the product of two terms. The first one gkjp ðhÞ given by (23) is a function of the material properties only, i.e., independent of location vectors x and y. The second term, given by sign [n Æ (x  y)]/sk Æ (x  y), is a simple algebraic function. It is in fact a plane wave function which is essentially a one-dimensional function. The interesting and simple structure of the integrands in (21) and (25) is exploited in this paper to analytically evaluate surface integrals over the boundary elements. Also notice that, since g is complex-valued (i.e. Ig 6¼ 0), the integrands in (21) and (25) are regular functions when n Æ (x  y) 5 0. Special situations when n Æ (x  y) = 0 are elaborated in some detail in Appendix A. 4. Piecewise linear triangular elements Substituting solutions (21) and (25) into (8) and (9), and exchanging integral orders, we obtain Z p X Z p X 3 3 ðmnÞ ðmnÞ ðmnÞ hk ðhÞwðmnÞ ðy; hÞdh; gjp ðyÞ ¼ R R gkjp ðhÞ/k ðy; hÞdh; hjp ðyÞ ¼ k jp 0

where ðmnÞ

/k

ðy; hÞ ¼ 

Z S ðmÞ

and ðmnÞ

wk

ðy; hÞ ¼

0

k¼1

Z S ðmÞ

ð27Þ

k¼1

  ðnÞ ðnÞ ðnÞ ðnÞ sign½n  ðx  yÞuðnÞ n1 ; n2 dn1 dn2 h   i ðnÞ ðnÞ s k  x n1 ; n2  y

  ðnÞ ðnÞ ðnÞ ðnÞ sign½n  ðx  yÞuðnÞ n1 ; n2 dn1 dn2 : n h   io2 ðnÞ ðnÞ s k  x n1 ; n2  y

ð28Þ

ð29Þ

In this section, we shall derive analytical solutions of integrals (28) and (29) for triangular elements with piecewise linear interpolation functions. 4.1. Geometry in local coordinates Consider a flat triangular element as shown in Fig. 1. Its three nodes (x(1), x(2), and x(3)) are arranged counterclockwise around the outward normal vector n so that

7078

C.-Y. Wang, M. Denda / International Journal of Solids and Structures 44 (2007) 7073–7091

X

n

(3)

e1 e2

(2)

X

x3

X

(1)

x2 x1 Fig. 1. A triangular element x(1)x(2)x(3) with the local orthogonal base vectors (e1, e2, n), where e1, e2 are an arbitrary pair of orthogonal vectors on the element plane and n is the unit normal to the element.

 ð2Þ     ð3Þ     ð1Þ    x  xð1Þ  xð3Þ  xð2Þ x  xð2Þ  xð1Þ  xð3Þ x  xð3Þ  xð2Þ  xð1Þ ¼ ¼ ; n¼ 2M 2M 2M

ð30Þ

where n is the area of the triangular element. This equation can be put into the following compact form  ðnÞ    x  xðn1Þ  xðnþ1Þ  xðnÞ n¼ ; ð31Þ 2M upon introducing the convention xð4Þ  xð1Þ ;

xð0Þ  xð3Þ :

ð32Þ

Define ðnÞ

q1 ¼

xðn1Þ  xðnþ1Þ ; jxn1  xðnþ1Þ j

ðnÞ

ðnÞ

q2 ¼ n  q1 : ðnÞ

ðnÞ

ð33Þ ðnÞ

Fig. 2 illustrates the geometry of q1 and q2 when n = 3. Notice that q1 is along the element edge across from ðnÞ ðnÞ the nth node. We use (q1 ; q2 ; n) as the unit orthogonal base vectors for the local coordinates. The coordinate transformation is given by ðnÞ

ðnÞ ðnÞ

ðnÞ ðnÞ

x ¼ x0 þ n1 q1 þ n2 q2 þ n3 n;

ð34Þ

ðnÞ x0

where the coordinate origin is given by (see Figs. 1 and 2) h h  ðnÞ i ðnÞ  ðnÞ i ðnÞ ðnÞ x0 ¼ xðnþ1Þ þ xðnÞ  xðnþ1Þ  q1 q1 ¼ xðnÞ  xðnÞ  xðnþ1Þ  q2 q2 :

ð35Þ

In such local coordinates, the interpolation function u(n) yields a simple form   nðnÞ ðnÞ ðnÞ uðnÞ n1 ; n2 ¼ 2ðnÞ ; h

ð36Þ

which, in turn, simplifies integration of (28) and (29). Some quantities we will use later on are (see Fig. 2)   ðnÞ hðnÞ ¼ xðnÞ  xðnþ1Þ  q2 ;

  lðnÞ ¼ xðn1Þ  xðnþ1Þ ;

  ðnÞ ðnÞ l1 ¼ xn  xðnþ1Þ  q1 ;

ðnÞ

ðnÞ

l2 ¼ lðnÞ  l1 :

ð37Þ

C.-Y. Wang, M. Denda / International Journal of Solids and Structures 44 (2007) 7073–7091

7079

ξ (3) 4 (3)

X

(3)

h

(3)

q2

(3)

q1

(3)

l1 (1)

(2)

X (3) 0

X ðnÞ

ξ(3) 3

(3)

l2

X

ðnÞ

ðnÞ

ðnÞ

Fig. 2. Local coordinate system ðn1 ; n2 ; n3 Þ with the unit vectors q1 and q2 and n, where n = 3 in this figure. The axis n3 and the unit normal n, perpendicular to the element, are not shown.

In the local coordinates, the source point y is expressed by ðnÞ

ðnÞ ðnÞ

ðnÞ ðnÞ

y ¼ x0 þ y 1 q1 þ y 2 q2 þ y 3 n;

ð38Þ

where ðnÞ

ðnÞ

ðnÞ

ðnÞ

y 1 ¼ ðy  x0 Þ  q1 ;

ðnÞ

ðnÞ

y 2 ¼ ðy  x0 Þ  q2 ; ðnÞ

ðnÞ

y 3 ¼ ðy  x0 Þ  n:

ð39Þ

ðnÞ

Recall that (e1, e2, n) and ðq1 ; q2 ; nÞ are both unit orthogonal base vectors. Henceforth, using Eqs. (22), (34) and (38), we get     ðnÞ ðnÞ ðnÞ ðnÞ sk  ðx  yÞ ¼ n1  y 1 cosðh  an Þ þ n2  y 2 sinðh  an Þ þ ðn3  y 3 Þgk ; ð40Þ ðnÞ

where an, the angles of q1 measured positive counterclockwise relative to e1, are determined by ðnÞ

ðnÞ

ðnÞ

cos an ¼ e1  q1 ¼ e2  q2 ;

ðnÞ

sin an ¼ e2  q1 ¼ e1  q2 :

ð41Þ

4.2. Analytical solutions of / and w Substitution of (36) and (40) into integral (28) yields ðmnÞ /k ðy; hÞ

signðy 3 Þ ¼ hðnÞ

Z 0

hn

Z

ðnÞ

ðnÞ

ðnÞ

ðnÞ ðnÞ

l2 l2 n2 =hðnÞ

l1 þl1 n2 =hðnÞ

ðnÞ

ðnÞ

ðnÞ

n2 dn1 dn2 ðnÞ

ðnÞ

ðnÞ

n1 cos hn þ n2 sin hn  1k

;

ð42Þ

where hn ¼ h  an

ð43Þ

and ðnÞ

ðnÞ

ðnÞ

1k ¼ y 1 cos hn þ y 2 sin hn þ y 3 gk :

ð44Þ

In writing expression (42), n3 = 0 for x 2 S(m) has been considered. By a straightforward integration of (42), we obtain ðmnÞ

/k

ðy; hÞ ¼

signðy 3 ÞhðnÞ ff2 ðhn Þ  f1 ðhn Þg; 2 cos hn

ð45Þ

7080

C.-Y. Wang, M. Denda / International Journal of Solids and Structures 44 (2007) 7073–7091

where

2 c2 ðhn Þ c ðhn Þ  2 ½log fc3 ðhn Þg  log fc2 ðhn Þg; b2 ðhn Þ b2 ðhn Þ

2 c1 ðhn Þ c1 ðhn Þ f1 ðhn Þ ¼  ½log fc3 ðhn Þg  log fc1 ðhn Þg; b1 ðhn Þ b1 ðhn Þ

f2 ðhn Þ ¼

ð46Þ

with ðnÞ

b1 ðhn Þ ¼ hðnÞ sin hn  l2 cos hn ; c1 ðhn Þ ¼

ðnÞ l1

cos hn 

ðnÞ 1k ;

ðnÞ

b2 ðhn Þ ¼ hðnÞ sin hn þ l1 cos hn ;

c2 ðhn Þ ¼

ðnÞ l2

cos hn 

ðnÞ 1k

ðnÞ 1k ;

c3 ðhn Þ ¼ h

ð47Þ ðnÞ

sin hn 

ðnÞ 1k :

ð48Þ

ðmnÞ /k ðy; hÞ

We note that in (42) is complex-valued if y3 5 0. Therefore given by (45) should be regular functions of h. There are a few locations of h where (45) is in the apparent singular form of 0/0, but these singularities can be removed rather easily by rearranging terms in (45). When y3 = 0, that is when y is on the element plane, (45) cannot be used directly. Special treatment for cases when y3 = 0 is given in Appendix B. Similarly, we can evaluate the integral (29) to find Z Z ðnÞ ðnÞ ðnÞ ðnÞ ðnÞ ðnÞ signðy 3 Þ hn l2 l2 n2 =h n2 dn1 dn2 ðmnÞ wk ðy; hÞ ¼ ð49Þ  2 ; ðnÞ ðnÞ ðnÞ ðnÞ ðnÞ ðnÞ hðnÞ 0 l1 þl1 n2 =hðnÞ n1 cos hn þ n2 sin hn  1k which has the following analytical solution

signðy 3 ÞhðnÞ f2 ðhn Þ f1 ðhn Þ ðmnÞ wk ðy; hÞ ¼  : c2 ðhn Þ c1 ðhn Þ cos hn

ð50Þ

Eqs. (45), (50) together with Eq. (27) are the key results of this paper. They show that the construction of the system equations (7) can be achieved by numerically computing the line integrals in (27). The integrands of the ðmnÞ ðmnÞ line integrals in (27) are given by gkjp ðhÞ/k ðy; hÞ and hkjp ðhÞwk ðy; hÞ, among which gkjp and hkjp are functions ðmnÞ ðmnÞ of the material properties only, while /k and wk are simple algebraic functions of the location and geometry. This simple structure facilitates the coding procedure and the computation for very large systems. ðnÞ ðmnÞ ðmnÞ We note again that 1k in (42) and (49) is complex-valued if y3 5 0. Therefore /k ðy; hÞ and wk ðy; hÞ that enter the line integrals (27) as part of the integrands are regular functions. When y3 = 0, formulas (45) and (50) cannot be used without a careful evaluation of the limits as y3 ! 0. Although the direct derivation of the limits for (45) and (50) as y3 ! 0 is possible, the derivation is somewhat complicated. As a simpler alternative, we have used an indirect derivation, shown in Appendix A and B, that first calculates the limits for Green’s functions as y3 ! 0. Having to deal with the case y3 = 0 is an inelegant and complicating aspect of the formulation. This situation occurs not only when y is on the element plane but also when it lies on its extended plane outside the element. Therefore it is necessary to keep track of this condition even when y is outside the element concerned. The verification results for a cube, for which this situation occurs frequently, are shown successfully in Section 6. 5. Implementation In solving the discretized system of equations (7), we follow the standard procedure as presented in Snyder and Cruse (1975) to construct the global system of equations that takes into account of possible discontinuities of traction at edges and corners of the boundary. Let N be the total number of global nodes and M the total number of elements. The system of linear equations are given by 2 11 38 ð1Þ 9 2 11 38 ð1Þ 9 u > t > G H H12 . . . H1N > G12 . . . G1M > > > > > > > > ð2Þ > > 6 21 6 21 22 2M 7> 22 2N 7> ð2Þ > = < < u H H . . . H G G . . . G 6 7 6 7 t = 6 . 7 6 7 ð51Þ .. 7> .. > ¼ 6 .. .. .. .. 7> .. >; .. .. 6 . 4 . . 5> . . . 5> . . > 4 . > > . > > . > > > > > : ðN Þ ; : ðMÞ ; u t HN 1 HN 2 . . . HNN GN 1 GN 2 . . . GNM

C.-Y. Wang, M. Denda / International Journal of Solids and Structures 44 (2007) 7073–7091 ðP 1Þ

ðP 1Þ

ðP 1Þ

ðP 3Þ

ðP 3Þ

7081

ðP 3Þ T

where tðP Þ ¼ ft1 ; t2 ; t3 ; . . . ; t1 ; t2 ; t3 g is a 1 · 9 element traction vector, which has three traction ðQÞ ðQÞ ðQÞ T components at each of the three nodes of an element P, and uðQÞ ¼ fu1 ; u2 ; u3 g is a 1 · 3 nodal displacePQ PR ment vector of node Q. The sub-matrices H and G (P, Q = 1, . . . , N and R = 1, . . . , 3N) have the dimensions 3 · 3 and 3 · 9, respectively. While the nodal displacement components are continuous, the traction components can be, in general, discontinuous at locations where the boundary is not smooth, such as corners or edges. The algorithm to handle discontinuous boundary traction components, as discussed by Snyder and Cruse (1975), is to allow only three unknown traction components at a node. Other traction components at the same node must be either treated as prescribed known quantities or assumed the same as the three unknown traction components. Eventually the system of equations (51) can be reduced to the following form: ½Afxg ¼ fyg;

ð52Þ

where [A] is a 3N · 3N matrix and {x} and {y} are the unknown and known vectors of dimension 3N. 6. Numerical results 6.1. Verification of accuracy Consider a unit cube (|xi| 6 1, i = 1, 2, 3) shown in Fig. 3 subjected to the uniform unit stress field r3 = 1 and introduce the boundary element mesh (14 nodes and 24 triangular elements) shown in Fig. 3. The material selected is aluminum crystal (cubic) with the material constants given by (56). The boundary conditions used are u1 ð0:0; 1:0; 0:0Þ ¼ u1 ð0:0; 1:0; 0:0Þ ¼ 0:0; u2 ð0:0; 0:0; 1:0Þ ¼ u2 ð0:0; 0:0; 1:0Þ ¼ 0:0;

ð53Þ

u3 ð1:0; 0:0; 0:0Þ ¼ u3 ð1:0; 0:0; 0:0Þ ¼ 0:0: The analytical displacement solution, with the boundary condition (53), is given by u1 ¼ s13 r3 x1 ;

u2 ¼ s23 r3 x2 ;

u3 ¼ s33 r3 x3 ;

ð54Þ

where s13, s23 and s33 are compliance constants and r3 = 1. The numerical evaluation of the line integrals in (27) is performed by dividing the range (0, p) into n sub-intervals, in each of which eight point Gauss quadrature is used. Notice that the analytical displacement solution is linear and the geometry involves straight lines only. Therefore, in the BEM solution with the linear interpolation and the flat triangular element, the only source of error is the line integral. With n = 4 for the line integral, the numerical results for the nodal displacement components are in agreement with the analytical values up to six significant digits at all 14 nodes giving the utmost confidence in the BEM developed. In this cube example, the case y3 = 0, discussed in Section 4, occurs frequently. The successful verification of the results confirms that the scheme developed to deal with this case is working. X3 σ3

X2

X1

σ3 Fig. 3. Cube under uniform stress r3. The BEM mesh has 14 nodes and 24 triangular elements.

7082

C.-Y. Wang, M. Denda / International Journal of Solids and Structures 44 (2007) 7073–7091

6.2. Spherical cavity Consider in an infinite elastic body a spherical cavity surface with radius a subjected to uniform pressure p. For an isotropic solid the radial displacement on the cavity surface is given by a uanal p; ð55Þ ¼ r 4G where G is the shear modulus. We should solve for this isotropic case by our proposed method and compare the solution against the analytical solution (55) and the numerical solution obtained by the standard BEM code based on the isotropic fundamental solutions. This isotropic BEM, also using triangular elements and linear interpolations, is identified as the IBEM below; while the BEM developed in this paper is identified as ABEM. Computations were carried out using a coarse mesh with 50 nodes and 96 elements, and with a fine mesh with 146 nodes and 288 elements. Three dimensional views of these meshes are shown in Fig. 4. Results for the radial displacement components are (IBEM = 0.860 vs AEBM = 0.875 with the coarse mesh. IBEM = 0.949 vs AEBM = 0.947 with the fine mesh. Analytical solution = 1.00). We observe good agreements between IBEM and ABEM results in both the coarse and fine meshes. With regard to the comparison to the analytical solution (55), results with the fine mesh model is much better than that with the coarse mesh, indicating a convergence to the analytical solution. We further carried out computations for a number of anisotropic materials with the fine mesh model as shown in Fig. 4(b). In those cases there are no analytical solutions to compare with. The materials considered are: aluminum polycrystal (isotropic), aluminum crystal (cubic), barium sodium niobate (orthorhombic), and sapphire (trigonal). The non-dimensional stiffness constants of these materials are given 3 2 11:1; 6:1; 6:1; 0:0; 0:0; 0:0 6 11:1; 6:1; 0:0; 0:0; 0:0 7 7 6 7 6 6 11:1; 0:0; 0:0; 0:0 7 7; ½cij =c0 al pcrystal ¼ 6 6 2:5; 0:0; 0:0 7 7 6 7 6 4 2:5; 0:0 5 2

½cij =c0 al crystal

2

½cij =c0 bsn

23:9;

6 6 6 6 ¼6 6 6 6 4 2

½cij =c0 sapphire

10:8;

6 6 6 6 ¼6 6 6 6 4

6 6 6 6 ¼6 6 6 6 4

2:5 3 6:13; 6:13; 0:0; 0:0; 0:0 10:8; 6:13; 0:0; 0:0; 0:0 7 7 7 10:8; 0:0; 0:0; 0:0 7 7; 2:85; 0:0; 0:0 7 7 7 2:85; 0:0 5 10:4; 24:7;

5:0; 5:2;

0:0; 0:0; 0:0; 0:0;

13:5; 0:0; 0:0; 6:5; 0:0; 6:6; 49:4;

15:8; 11:4; 2:3; 49:4; 11:4; 2:3; 49:6;

0:0; 14:5;

2:85 3 0:0 0:0 7 7 7 0:0 7 7; 0:0 7 7 7 0:0 5

7:6 3 0:0; 0:0 0:0; 0:0 7 7 7 0:0; 0:0 7 7: 0:0; 0:0 7 7 7 14:5; 2:3 5 16:8

ð56Þ

C.-Y. Wang, M. Denda / International Journal of Solids and Structures 44 (2007) 7073–7091

a

7083

b Z

Y X

Fig. 4. (a) Coarse and (b) fine meshes on the spherical cavity surface. The normals point into the cavity.

Given the displacement u, length x, pressure p, and stiffness c, non-dimensional quantities are introduced as bellow u x p c  p ¼ ; c ¼ ; ð57Þ u ¼ ; x ¼ ;  u0 a p0 c0 where u0 = 103a (m), p0 = 107 (N/m2), and c0 = 1010 (N/m2) are the reference displacement, pressure, and stiffness, respectively. A non-dimensional pressure p ¼ 1 is applied to the spherical cavity surface with a non dimensional radius  a = 1. The displacement components are sampled at 24 element vertices located on the yz, xz, and xy-planes as shown in Fig. 5(a). Shown in Fig. 6 are the results for polycrystal aluminum, which is an isotropic material. Fig. 6(a) shows the variation of the displacement components ur, ut, and uo, which are the radial, tangential, and out-of-plane components, respectively, on the cavity surface measured at the element vertices located on the xy-plane. The horizontal axis of Fig. 6 given by the angle h is defined in Fig. 5(b). For the xy-plane, the out-of-plane direction is along the direction of the z-axis. The displacement plots on the yz-plane are given in Fig. 6(b). As expected from the material isotropy, these two figures are identical. Similar results for the crystal aluminum are shown in Fig. 7. Notice the difference in the magnitudes of the radial displacement between the polycrystal and crystal aluminums. The results in the xz-plane are identical to those in the yz-plane. The barium sodium niobate is an orthorhombic material with c33 drastically different from c11  c22. We expect significant difference between the xy-plane and yz-plane plots, while the xz-plane plot will be almost identical to the yz-plane plot. Fig. 8, showing the results in the xy-plane and yz-plane, verifies our expectations. The xz-plot, not shown, is almost identical to the yz-plot, also as expected.

z

y

φ

y

x

θ

x

Fig. 5. (a) Top view of the fine mesh for the spherical cavity. The displacement components are plotted for the element vertices on the yz, xz, and xy-planes. (b) Spherical coordinate system.

7084

C.-Y. Wang, M. Denda / International Journal of Solids and Structures 44 (2007) 7073–7091 0.1

0.1

0.09

0.09

ur

0.07

0.07

0.06

0.06

0.05

0.05

0.04

0.04

0.03

0.03

0.02

0.02

ut

0.01

ur

0.08

ur u t u o

ur ut u o

0.08

ut

0.01

uo

uo

0

0

-0.01

-0.01 -100

0

-100

100

θ

0

100

φ

Fig. 6. Aluminum polycrystal: Radial ur, tangential ut, and out-of-plane uo displacement components on the xy-plane (a) and on the yzplane (b).

0.1

0.1

0.09

0.09 0.08

ur

0.07

0.07

0.06

0.06

ur ut uo

ur ut uo

0.08

0.05 0.04 0.03

0.05 0.04 0.03

0.02

0.02

ut

0.01

uo

0 -0.01

ur

ut

0.01

uo

0 -100

0

θ

100

-0.01 -100

0

φ

100

Fig. 7. Aluminum crystal: Radial ur, tangential ut, and out-of-plane uo displacement components on the xy-plane (a) and on the yz-plane (b).

Finally we consider sapphire, which is a trigonal material. Fig. 9 shows the plots for the xy-plane and yzplane. Notice, in the yz-plane plot of Fig. 9(b), the presence of the out-of-plane component of the cavity surface displacement. This is expected from the trigonal property of sapphire. 7. Concluding remarks In this paper we have presented a new scheme of three-dimensional boundary element method for the elastic materials of general anisotropy. Green’s functions, expressed by the line integrals over a semi-circle, have a simple structure that allows analytical integration over the boundary elements. The computation for the system matrices G and H is reduced to numerical integration of regular line integrals over the unit semi-circle. This reduction from the surface to line integrations makes the proposed BEM an attractive alternative to

C.-Y. Wang, M. Denda / International Journal of Solids and Structures 44 (2007) 7073–7091

7085

0.05

0.035 0.03

0.04

ur

ur

0.025 0.03

urutuo

urutuo

0.02 0.015 0.01

0.02

0.01

0.005

ut

ut

uo 0

0

uo -0.01

-0.005 -100

0

-100

100

θ

0

100

φ

Fig. 8. Barium sodium niobate: Radial ur, tangential ut, and out-of-plane uo displacement components on the xy-plane (a) and on the yzplane (b). 0.02

0.02

0.018

0.018

ur

ur

0.016

0.014

0.014

0.012

0.012

urutuo

urutuo

0.016

0.01 0.008 0.006

0.01 0.008 0.006

0.004

ut

0.002

0.004

uo

ut

0.002

0

uo

0 -100

0

θ

100

-100

0

φ

100

Fig. 9. Sapphire: Radial ur, tangential ut, and out-of-plane uo displacement components on the xy-plane (a) and on the yz-plane (b).

the existing 3D BEMs that use the surface integrals of the explicit Green’s functions. The integrands of the line ðmnÞ ðmnÞ integrals in (27) are given by gkjp ðhÞ/k ðy; hÞ and  hkjp ðhÞwk ðy; hÞ, among which gkjp and hkjp are functions of the ðmnÞ ðmnÞ material properties only, while /k and wk are simple algebraic functions of the location and geometry. This simple structure facilitates the coding procedure and the computation for very large systems. The analytðmnÞ ðmnÞ ical expressions for /k ðy; hÞ and wk ðy; hÞ, given by Eqs. (45) and (50), along with the integrals in (27) are the key results in this paper. The proposed method has been implemented and tested successfully with numerical examples. Appendix A. Green’s functions when [~ y3 ¼ n  ðx  yÞ ! 0] Let ~y 3 ¼ n  ðx  yÞ; and let

ðA:1Þ

7086

C.-Y. Wang, M. Denda / International Journal of Solids and Structures 44 (2007) 7073–7091

p ¼ ðcos he1 þ sin he2 Þ  ðx  yÞ ¼ e1  ðx  yÞ cos h þ e2  ðx  yÞ sin h:

ðA:2Þ

It follows sign½n  ðx  yÞ ¼ signð~y 3 Þ;

ðA:3Þ

sk  ðx  yÞ ¼ ðcos he1 þ sin he2 þ gk nÞ  ðx  yÞ ¼ p  gj ~y 3 :

ðA:4Þ

and

Based on Eq. (17) and Iðgk Þ > 0, we find sign½n  ðx  yÞ ¼  lim ~y 3 !0 sk  ðx  yÞ sign½n  ðx  yÞ o lim ¼ lim ~y 3 !0 fsk  ðx  yÞg2 ~y 3 !0 op lim

~y 3 !0

signð~y 3 Þ signð~y 3 Þ  ipdðpÞ; ¼} ðp  gj ~y 3 Þ p signð~y 3 Þ signð~y 3 Þ ¼} þ ipd0 ðpÞ: ðp  gj ~y 3 Þ p2

ðA:5Þ

Substituting the two equations above to the fundamental solutions (21) and (25), we obtain ) ) Z p (X Z p (X 3 3 dh k k p lim gjp ðx; yÞ ¼ þsignð~y 3 Þ} R I gjp ðhÞ gjp ðhÞ dðpÞdh; ~y 3 !0 p 0 0 k¼1 k¼1 ) ) Z p (X Z p (X 3 3 dh k k  h ðhÞ d0 ðpÞdh: R p I lim hjp ðx; yÞ ¼ signð~y 3 Þ} hjp ðhÞ jp 2 ~y 3 !0 p 0 0 k¼1 k¼1 We can further simplify these equations to ) Z p (X 3 k I gjp ðhÞ dðpÞdh; lim gjp ðx; yÞ ¼ p ~y 3 !0

0

lim hjp ðx; yÞ ¼ p

~y 3 !0

Z

p

( I

0

k¼1 3 X

)  hkjp ðhÞ

d0 ðpÞdh 

k¼1

ðA:6Þ djp dðe1  ðx  yÞÞdðe2  ðx  yÞÞ; 2

using the following three equalities: ( ) ( ) 3 3 X X 1 k k  hjp ðhÞ ¼ 2 djp ; gjp ðhÞ ¼ 0; R R 4p k¼1 k¼1 and 1 } 2p2

Z 0

p

dh ¼ dðe1  ðx  yÞÞdðe2  ðx  yÞÞ: p2

We shall first prove (A.8). From (23) and (26) we observe, based on Cauchy’s residue theorem ( ) ( ) Z 3 3 X X 1 Ajp ðsk Þ 1 1 Ajp ðsÞ k dg; R gjp ðhÞ ¼ R ¼ 2 D ðsk Þ 2 DðsÞ 2p 4pi 2p ;g C k¼1 k¼1 ( ) ( ) Z 3 3 X X 1 ni cijlq sq Alp ðsÞ k k k  dg; R ni cijlq sq glp ðhÞ ¼ hjp ðhÞ ¼ R 4pi 2p2 DðsÞ C k¼1 k¼1

ðA:7Þ

ðA:8Þ

ðA:9Þ

ðA:10Þ ðA:11Þ

where C is a closed Jordan contour in the complex g-plane that includes all the poles at g = gj and g ¼ gj . Now, let g = |g|eih and let |g| ! 1. Thus dg = i|g|eih dh and Z Z 2p Z 2p 4 i4h 1 1 Ajp ðsÞ 1 1 ajp 1 1 ajp ih 4 jgj e 4 dh dg ¼ lim lim jgje dh ¼ ¼ 0; ðA:12Þ 6 2 2 2 i6h 4pi C 2p DðsÞ 4p jgj!1 0 2p d 6 jgj e 4p jgj!1 0 2p d 6 jgjeih Z Z 2p 4 i4h 1 ni cijlq sq Alp ðsÞ 1 ni cijlq nq jgjeih alp ni cijlq nq alp 4 jgj e 4 ih dg ¼ lim jgje dh ¼ ; ðA:13Þ 4pi C 2p2 DðsÞ 4p jgj!1 0 2p2 4p2 d 6 d 6 jgj6 ei6h

C.-Y. Wang, M. Denda / International Journal of Solids and Structures 44 (2007) 7073–7091

7087

where d6 is the leading coefficient of the polynomial D(g) as shown in (20), and ajp 4 are the leading coefficients of Ajp(s(g)). It can be shown using (12) and (22) that alp 4 ¼ adj Clp ðnÞ;

d 6 ¼ det Cjp ðnÞ:

ðA:14Þ

Therefore, (A.13) reduces to ni cijlq nq alp 4 =d 6 ¼ Cjl ðnÞadj Clp ðnÞ=d 6 ¼ djp det Cmn ðnÞ=d 6 ¼ djp :

ðA:15Þ

With this, we conclude the proof of the two equalities in (A.8). Some of the arguments used in the above between (A.10)–(A.13) can be found in Wang and Achenbach (1995) and Wang et al. (1996). We now proceed to prove (A.9). Let ~x1 ¼ e1  ðx  yÞ and let ~x2 ¼ e2  ðx  yÞ. Then, (A.9) reads Z p Z p Z p 1 dh 1 dh 1 } ¼ } ¼  D logð~x1 cos h þ ~x2 sin hÞdh; ðA:16Þ 2 2p2 p2 2p2 2p2 x1 cos h þ ~x2 sin hÞ 0 0 ð~ 0 where D  o~2x1 þ o~2x2 is the Laplace operator. Now, introducing the substitution z = ctan(h), we have dh = dz/(1 + z2) and Z p Z 1 1 1 logð~x1 z þ ~x2 Þ  2D logð~x1 cos h þ ~x2 sin hÞdh ¼  2 DR dz: ðA:17Þ 2p 2p 1 þ z2 0 1 A subsequent integration using Cauchy’s residue theorem yields Z 1 qffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 logð~x1 z þ ~x2 Þ 1 ~x21 þ ~x22 ¼ dð~x1 Þdð~x2 Þ: D log  2 DR dz ¼  2p 1 þ z2 2p 1

ðA:18Þ

This concludes the proof of (A.9). Derivations similar to (A.17), (A.18) can be found in Wang and Achenbach (1996). Appendix B. Solutions when y is on the element plane Substitution of (36), (A.6) and (A.7) into (8) and (9) yields ) Z p (X 3 ðmnÞ ðmnÞ lim gjp ðyÞ ¼ p I gkjp ðhÞ /k ðy; hÞdh; y 3 !0

0

k¼1 ðnÞ

lim

y 3 !0

ðmnÞ hjp ðyÞ

djp y 2 ¼ p 2 hðnÞ

Z

p

( I

3 X

) ðmnÞ k  hjp ðhÞ wk ðy; hÞdh;

0

k¼1

ðnÞ

  ðnÞ ðnÞ ðnÞ ðnÞ d n1 cos hn þ n2 sin hn  y 1 cos hn  y 2 sin hn n2 dn1 dn2 ;

ðB:19Þ ðB:20Þ

where ðmnÞ

/k

ðmnÞ

wk

ðy; hÞ ¼ ðy; hÞ ¼

1 hðnÞ 1 hðnÞ

Z

hðnÞ

0

ðnÞ

l2 ð1n2 =hðnÞ Þ ðnÞ

ðnÞ

l1 ð1þn2 =hðnÞ Þ

0

Z

Z

hðnÞ

Z

ðnÞ ðnÞ l2 ð1n2 =hðnÞ Þ ðnÞ

ðnÞ

l1 ð1þn2 =hðnÞ Þ

  ðnÞ ðnÞ ðnÞ ðnÞ d0 n1 cos hn þ n2 sin hn  y 1 cos hn  y 2 sin hn n2 dn1 dn2 :

Apply the coordinate transformation ðnÞ

ðnÞ

ðnÞ

ðnÞ

ðnÞ

ðnÞ

t1 ¼ cos hn n1 þ sin hn n2 ; t2 ¼  sin hn n1 þ cos hn n2 ; and rewrite the above integrals as

ðB:21Þ

7088

C.-Y. Wang, M. Denda / International Journal of Solids and Structures 44 (2007) 7073–7091

ðmnÞ /k ðy; hÞ

1 ¼ ðnÞ h

Z

þ1

d



ðnÞ t1



ðnÞ t1y

Z

1

t

ðnÞ

ðnÞ

t2U ðt1 Þ



ðnÞ ðnÞ ðt1 Þ

ðnÞ

ðnÞ

sin hn t1 þ cos hn t2

h n o ðnÞ ðmnÞ H t1  t1L

2L n oi ðnÞ ðmnÞ ðnÞ ðnÞ H t1  t1U dt2 dt1 ;

ðmnÞ

wk

ðy; hÞ ¼

1 hðnÞ

Z

þ1

 Z ðnÞ ðnÞ d0 t1  t1y

1

t

ðB:22Þ

ðnÞ

ðnÞ

t2U ðt1 Þ



ðnÞ ðnÞ ðt1 Þ

ðnÞ

ðnÞ

sin hn t1 þ cos hn t2

h n o ðnÞ ðmnÞ H t1  t1L

2L n oi ðnÞ ðmnÞ ðnÞ ðnÞ H t1  t1U dt2 dt1 ;

ðB:23Þ

where H is the Heaviside step function and ðnÞ

ðnÞ

ðnÞ

t1y ¼ cos hn y 1 þ sin hn y 2 : ðnÞ

ðnÞ

ðB:24Þ ðnÞ

ðnÞ

ðnÞ

ðnÞ

The integration limits t2L ðt1 Þ and t2U ðt1 Þ depend on t1 and the lower and the upper limits of t1 are given by ðmnÞ ðmnÞ t1L and t1U , respectively, as shown in Fig. 10. Let Z tðnÞ ðtðnÞ Þ   2U 1 ðnÞ ðnÞ ðnÞ ðnÞ F ðt1 Þ ¼ sin hn t1 þ cos hn t2 dt2 ðB:25Þ ðnÞ ðnÞ

t2L ðt1 Þ

¼

ðnÞ sin hn t1

n n o 1 o2 n o2 ðnÞ ðnÞ ðnÞ ðnÞ ðnÞ ðnÞ ðnÞ ðnÞ t2U ðt1 Þ  t2L ðt1 Þ þ cos hn t2U ðt1 Þ  t2L ðt1 Þ : 2

Its derivative is ðnÞ

oF ðt1 Þ ðnÞ

ot1

n o n o ðnÞ ðnÞ ðnÞ ðnÞ ðnÞ ðnÞ ðnÞ ðnÞ ðnÞ ¼ sin hn t2U ðt1 Þ  t2L ðt1 Þ þ sin hn t1 t_2U ðt1 Þ  t_2L ðt1 Þ h i ðnÞ ðnÞ ðnÞ ðnÞ ðnÞ ðnÞ ðnÞ ðnÞ þ cos hn t2U ðt1 Þt_2U ðt1 Þ  t2L ðt1 Þt_2L ðt1 Þ ;

ðB:26Þ

ðnÞ

where the dot indicates differentiation with respect to t1 . Eqs. (B.25) and (B.26) can be simplified to    ðnÞ ðnÞ ðnÞ ðnÞ ðnÞ ðnÞ ðnÞ ðB:27Þ F ðt1 Þ ¼ 12 n2U þ n2L t2U ðt1 Þ  t2L ðt1 Þ ; ðnÞ

oF ðt1 Þ ðnÞ

ot1

n o ðnÞ ðnÞ ðmnÞ ðnÞ ðnÞ ðnÞ ðnÞ ðnÞ ðnÞ ðmnÞ ¼ sin hn t2U ðt1 Þ  t2L ðt1 Þ þ t_2U ðt1 Þn2U  t_2L ðt1 Þn2L ;

ðB:28Þ

ξ (3) 2 t (3) 2

X(3)

(3)

t 2U

t 1(3) (3)

t 1y X

θn

(1)

(3)

y

t 2L

(m3)

t 1U

(m3) t 1L

ðnÞ

ðnÞ

(3)

X(2) ξ 1

ðnÞ

ðnÞ

Fig. 10. Local coordinate system (t1 , t2 , n3) obtained by rotating (n1 , n2 , n3) by an angle hn, where n = 3 in this figure. The axis n3, perpendicular to the element, is not shown.

C.-Y. Wang, M. Denda / International Journal of Solids and Structures 44 (2007) 7073–7091

7089

where ðnÞ

ðnÞ

ðnÞ

ðnÞ

ðnÞ

n2U ¼ sin hn t1 þ cos hn t2U ðt1 Þ;

ðnÞ

ðnÞ

ðnÞ

n2L ¼ sin hn t1 þ cos hn t2L ðt1 Þ:

ðB:29Þ

Evaluating integrals (B.22) and (B.23), we obtain ðmnÞ

/k

ðy; hÞ ¼ ¼

1 hðnÞ

Z

þ1

1

  h n o n oi ðnÞ ðnÞ ðnÞ ðnÞ ðmnÞ ðnÞ ðmnÞ ðnÞ d t1  t1y F ðt1 Þ H t1  t1L  H t1  t1U dt1

ðnÞ F ðt1y Þ h n ðnÞ H t1y ðnÞ

ðmnÞ

 t1L

h

o

n oi F ðtðnÞ Þ 1y ðnÞ ðmnÞ  H t1y  t1U ¼ ðnÞ h

ðB:30Þ

and Z

  h n o n oi ðnÞ ðnÞ ðnÞ ðnÞ ðmnÞ ðnÞ ðmnÞ ðnÞ d0 t1  t1y F ðt1 Þ H t1  t1L  H t1  t1U dt1 1 ( ) ðnÞ n oi F ðtðnÞ Þ h n o n oi 1 oF ðt1y Þ h n ðnÞ ðmnÞ o 1y ðnÞ ðmnÞ ðnÞ ðmnÞ ðnÞ ðmnÞ ¼  ðnÞ  t  t d t  t  t  H t   d t H t 1y 1L 1y 1U 1y 1L 1y 1U ðnÞ hðnÞ h ot1 ( ) ðnÞ ðnÞ n oi F ðt1y Þ h n ðnÞ ðmnÞ o 1 oF ðt1y Þ ðnÞ ðmnÞ d t  t  t ¼  ðnÞ   d t ; ðB:31Þ 1y 1L 1y 1U ðnÞ hðnÞ h ot1

ðmnÞ wk ðy; hÞ ¼

1 hðnÞ

þ1

where ðnÞ

oF ðt1y Þ ðnÞ

ot1

 ðnÞ oF ðt1 Þ ¼  ðnÞ  ot 1

t1 ¼t1y

. Substituting (B.30) and (B.31) into (B.19) and (B.20) we obtain lim

y 3 !0

lim

y 3 !0

ðmnÞ gjp ðyÞ

ðmnÞ hjp ðyÞ

p ¼  ðnÞ h

Z

p

( I

0

3 X

) gkjp ðhÞ

ðnÞ

F ðt1y Þdh;

ðB:32Þ

k¼1

)( Z p (X ðnÞ ðnÞ 3 h n o oF ðt1y Þ djp y 2 p ðnÞ ðnÞ ðmnÞ k  ¼ þ I ðhÞ þ F ðt Þ d t  t h 1y 1y 1L jp ðnÞ 2 hðnÞ hðnÞ 0 ot1 k¼1 n oio ðnÞ ðmnÞ d t1y  t1U dh: ðnÞ

ðB:33Þ ðmnÞ

Note that the delta functions in (B.33) pick up nonzero contributions only when t1y is equal to either t1L or ðmnÞ t1U . This is the case when the observation point y is located on one of the edges of the element or their extenðnÞ sions including the vertices as shown in Fig. 11. Since t1y , given by (B.24), is the function of hn = h-an we can rewrite (B.33) as ) Z p (X ðnÞ 3 djp y 2 p ðmnÞ k  ðhÞ lim hjp ðyÞ ¼  þ I ðB:34Þ h jp y 3 !0 2 hðnÞ hðnÞ 0 k¼1 8 2 39 > > < o dfhn  hnU g7= ðnÞ ðnÞ 6dfhn  hnL g      F ðt Þ þ F ðt Þ  4 5 1y 1y ðnÞ  ðnÞ 0   ðnÞ 0  ;dh; > :ot1 t1y ðhnL Þ t1y ðhnU Þ >

7090

C.-Y. Wang, M. Denda / International Journal of Solids and Structures 44 (2007) 7073–7091 (3)

ξ2 X(3)

t (3) 2

t (3) 3

(3)

t 2U

t(3) 1y

θ 3U

X(1)

X(2)

(3)

ξ1

x0(3) (m3)

y

ðmnÞ

ðmnÞ

t 1U

(3)

t 2L

(m3)

t 1L ðnÞ

Fig. 11. When y is on the element edge or its extension, t1y takes the extreme value, t1L or t1U , at an angle hnL or hnU. Then the integral ðmnÞ ðnÞ ðmnÞ ðnÞ ðmnÞ hjp ðyÞ picks up an additional contribution from the delta function, dðt1y  t1L Þ or dðt1y  t1U Þ, in (B.33). The case when n = 3 and ð3Þ ðm3Þ t1y ¼ t1U is shown in this figure. ðnÞ

ðmnÞ

ðnÞ

ðmnÞ

where hnL and hnU are the roots of t1y ðhn Þ  t1L ¼ 0 and t1y ðhn Þ  t1U ¼ 0 and  ðnÞ dt1y  ðnÞ 0 ðnÞ ðnÞ t1y ðhnL Þ ¼ ¼  sin hnL y 1 þ cos hnL y 2 ;  dhn  h ¼h  n nL ðnÞ  dt1y  ðnÞ 0 ðnÞ ðnÞ ¼  sin hnU y 1 þ cos hnU y 2 : t1y ðhnU Þ ¼  dhn 

ðB:35Þ

hn ¼hnU

Further evaluation of (B.34) gives lim

y 3 !0

ðmnÞ hjp ðyÞ

) Z p (X ðnÞ ðnÞ 3 oF ðt1y Þ djp y 2 p k  ¼ þ I ðhÞ dh h jp ðnÞ 2 hðnÞ hðnÞ 0 ot1 k¼1 8 ( 9 ) ( ) ðnÞ ðnÞ = 3 3 X F ft ðh Þg F ft ðh Þg p < X nL nU 1y 1y  hk ðhU Þ   I  ; þ ðnÞ I hkjp ðhL Þ  0 jp 0   ðnÞ  ðnÞ h : k¼1 k¼1 t1y ðhnL Þ t1y ðhnU Þ ;

ðB:36Þ

where, from (43), hL ¼ hnL þ an ; hU ¼ hnU þ an :

ðB:37Þ

References Barnett, D.M., 1972. The precise evaluation of derivatives of the anisotropic elastic Green’s functions. Phys. Stat. Sol. (b) 49, 741–748. Berger, J.R., Tewary, V.K., 1997. Boundary integral equation formulation for interface cracks in anisotropic materials. Int. J. Numer. Meth. Eng. 20, 261–266. Chen, T., Lin, F.Z., 1995. Boundary integral formulations for three-dimensional anisotropic solids. Comput. Mech. 15, 485–496. Clements, D.L., Haselgrove, M.A., 1983. A boundary integral equation for a class of crack problems in anisotropic elasticity. Int. J. Comp. Math. 12, 267–278. Denda, M., Wang, C.-Y., Yong, Y.C., 2003. 2-D time-harmonic BEM for solids of general anisotropy with application to eigenvalue problems. J. Sound Vib. 261, 247–276. Denda, M., Marante, M.E., 2004. Mixed mode BEM analysis of multiple curvilinear cracks in the general anisotropic solids by the crack tip singular element. Int. J. Solids Struct. 41 (5-6), 1473–1489.

C.-Y. Wang, M. Denda / International Journal of Solids and Structures 44 (2007) 7073–7091

7091

Denda, M., Araki, Y., Yong, Y.C., 2004. Time-harmonic BEM for 2-D piezoelectricity applied to eigenvalue problems. Int. J. Solids Struct. 41 (26), 7241–7265. Eshelby, J.D., Read, W.T., Shockley, W., 1953. Anisotropic elasticity with applications to dislocation theory. Acta Metall. 1, 251–259. Gray, L.J., Griffith, A., Johnson, L., Wawrzynek, P.A., 2003. Evaluation of Galerkin singular integrals for anisotropic elasticity: displacement equation. Electron. J. Bound. Elem. 1 (2), 68–94. Gel’fand, I.M., Graev, M.I., Vilenkin, N.Y., 1966. In: Generalized Functions, vol. 5. Academic Press, New York. Hirose, S., Wang, C.-Y., Achenbach, J.D., 2000. Boundary element method for elastic wave scattering by a crack in an anisotropic solid. In: Proceedings of the 20th Int. Congress Theoret. Appl. Mech. ICTAM 2000, Chicago. Lee, V.G., 2003. Explicit expression of derivatives of elastic Green’s functions for general anisotropic materials. Mech. Res. Commun. 30, 241–249. Li, S., Mear, M.E., Xiao, L., 1998. Symmetric weak-form integral equation method for three-dimensional fracture analysis. Comput. Meth. Appl. Mech. Eng. 151, 435–459. Nakamura, G., Tanuma, K., 1997. A formula for the fundamental solution of anisotropic elasticity. Q. J. Mech. Appl. Math. 50, 179–194. Pan, E., Yuan, F.G., 2000. Boundary element analysis of three-dimensional cracks in anisotropic solids. Int. J. Numer. Meth. Eng. 48, 211–237. Pan, Y.C., Chou, T.W., 1976. Point force solution for an infinite transversely isotropic solid. J. Appl. Mech. 43, 608–612. Sales, M.A., Gray, L.J., 1998. Evaluation of the anisotropic Green’s function and its derivatives. Comput. Struct. 69, 247–254. Snyder, M.D., Cruse, T.A., 1975. Boundary integral equation analysis of cracked anisotropic plates. Int. J. Fract. 11, 315–328. Sollero, P., Aliabadi, M.H., 1993. Fracture mechanics analysis of anisotropic plates by the boundary element method. Int. J. Fract. 64 (4), 269–284. Stroh, A.N., 1958. Dislocations and cracks in anisotropic elasticity. Phil. Mag. 7, 625–646. Synge, J.L., 1957. The Hypercircle in Mathematical Physics. Cambridge University Press, London. Ting, T.C.T., Lee, V.G., 1997. The three-dimensional elastostatic Green’s function for general anisotropic linear elastic solids. Q. J. Mech. Appl. Math. 50, 407–426. Ting, T.C.T., 1996. Anisotropic Elasticity: Theory and Applications. Oxford University Press, Oxford. Tonon, F., Pan, E., Amadei, B., 2001. Green’s functions and boundary element method formulation for 3D anisotropic media. Comput. Struct. 79, 469–482. Vogel, S.M., Rizzo, F.J., 1973. An integral equation formulation of threedimensional anisotropic elastostatic boundary value problems. J. Elast. 3, 203–216. Wang, C.-Y., Achenbach, J.D., 1992. A new look at 2-D time-domain elastodynamic Green’s functions for general anisotropic solids. Wave Motion 16, 389–405. Wang, C.-Y., Achenbach, J.D., 1994. Elastodynamic fundamental solutions for anisotropic solids. Geophys. J. Int. 118, 384–392. Wang, C.-Y., 1994. 2-D elastostatic Green’s functions for general anisotropic solids and generalization of Strohs formalism. Int. J. Solids Struct. 31 (19), 2591–2597. Wang, C.-Y., Achenbach, J.D., 1995. Three-dimensional time-harmonic elastodynamic Green’s functions for anisotropic solids. Proc. R. Soc. Lond. A 449, 441–458. Wang, C.-Y., Achenbach, J.D., Hirose, H., 1996. Two-dimensional time domain BEM for scattering of elastic waves in solids of general anisotropy. Int. J. Solids Struct. 33 (26), 3843–3864. Wang, C.-Y., Achenbach, J.D., 1996. Lamb’s problem in solids of general anisotropy. Wave Motion 24 (3), 227–242. Wang, C.-Y., 1997. Elastic fields produced by a point source in solids of general anisotropy. J. Eng. Math. 32 (1), 41–52. Wilson, R.B., Cruse, T.A., 1978. Efficient implementation of anisotropic three dimensional boundary-integral equation stress analysis. Int. J. Numer. Meth. Eng. 12, 1383–1397.