INTERNATIONAL JOURNAL FOR NUMERICAL METHODS IN ENGINEERING, VOL. 40, 89–109 (1997) TREATMENT OF BODY-FORCE VOLUME INTEGRALS IN BEM BY EXACT TRANSFORMATION FOR 2-D ANISOTROPIC ELASTICITY J. J. ZHANG, C. L. TAN AND F. F. AFAGH Department of Mechanical and Aerospace Engineering, Carleton University, 3135 Mackenzie Bldg, 1125 Colonel By Drive, Ottawa, Canada K1S 5B6 SUMMARY In the Boundary Element Method (BEM) based on the direct formulation, body-force eects manifest themselves as an additional volume integral term in the Boundary Integral Equation (BIE). The numerical solution of the integral equation with this term destroys the notion of the BEM as a truly boundary solution method. This paper discusses the treatment of this volume integral for two-dimensional anisotropic elasticity with bodyforces present. The analytical basis for transforming this integral exactly into boundary ones is presented for geometrically convex regions. This restores the application of the BEM to such problems as a truly boundary solution technique. Numerical examples are presented to demonstrate the veracity of the transformation and implementation. KEY WORDS: BEM; numerical fracture mechanics 1. INTRODUCTION The Boundary Element Method (BEM), also known as the Boundary Integral Equation (BIE) method, is an ecient numerical tool for engineering stress analysis. In elastostatics, the analytical basis of the so-called direct BIE formulation, which the present work is focused on, is the existence of the fundamental solution to the governing dierential equations, and the reciprocal work theorem. From these, and following the usual limiting process, Somigliana’s identity and the BIE, which relate the boundary displacements to the boundary tractions, may be obtained. The fact that integrals present in the BIE and Somigliana’s identity are surface integrals establishes the BEM as a boundary solution technique as distinct from domain solution techniques, such as the nite element method. The implication here is that only the surface of the domain needs to be modelled or discretized, instead of the entire domain as in FEM. However, at the basic level of formulation, this holds true only if body forces and temperature changes are absent. When body forces and thermal eects are considered, the basic form of Somigliana’s identity and the BIE will not only include surface integrals, but volume integrals as well. The direct numerical evaluation of these integrals necessitates the discretization of the entire domain into ‘cells’, so that integration over each of these ‘cells’ may be carried out. Although these integration ‘cells’ are not ‘nite elements’, since they do not contribute to the coecients corresponding to the unknowns in the nal equation matrix, the notion of BEM as a boundary solution technique is destroyed in such cases. Thus the solution to the ‘volume integral problem’ is of fundamental importance to the development of BEM. CCC 0029–5981/97/010089–21 ? 1997 by John Wiley & Sons, Ltd. Received 27 July 1995 Revised 8 February 1996 90 J. J. ZHANG, C. L. TAN AND F. F. AFAGH There are evidently two dierent general schemes to resolve the volume integral problem in the BIE formulation. The rst is to attempt to perform the volume integration without any explicit discretization of the solution domain, such as the Monte Carlo Method, developed by Gipson and Camp,1 or the domain Fanning method, proposed by Camp and Gipson.2 Although a very imaginative idea, the Monte Carlo method is not ecient because a large number of calculation points must be randomly chosen, and the accuracy of the solution is not always assured. The domain Fanning method takes advantage of the properties of the fundamental solution and improves the eciency of evaluating the volume integral through an implicit domain discretization scheme. However, for domains with certain boundary shapes, diculties will arise. The second and more attractive procedure is to further transform the volume integrals to surface or boundary integrals. The exact transformation of volume integrals to boundary integrals, by Rizzo and Shippy,3 the dual reciprocity method, by Nardini and Brebbia,4 the multiple reciprocity method, by Nowak,5 provide examples of this general scheme. These methods have met with tremendous success in isotropic elasticity. However, the developments for anisotropic elasticity have, hitherto, not been that successful. This is primarily due to the signicantly more complicated form of the fundamental solutions in anisotropic cases. A limited amount of work in this regard has been reported recently, using the ‘particular integral approach’. It was rst discussed by Watson,6 and later by Banerjee and Buttereld,7 for isotropic elasticity. This approach has recently been further extended to anisotropic elasticity by Deb and Banerjee,8 and Deb et al.9 Yet another approach, developed by Ao17 very recently, is to reformulate the elasticity equations into a BIE which is in terms of unknown scalar potentials to be solved for. The method, as presented in Reference 17 and Ao and Shippy,18 applies, however, to only orthotropic materials. The particular integral approach mentioned above has been successful in solving problems where the stress distribution is smooth in anisotropic media. However, it lacks generality to treat problems involving stress singularities such as those encountered in cracks. The motivation of the present work is to overcome these diculties, and an attempt to develop a method which transforms exactly the volume integrals to boundary ones in 2-D anisotropic elasticity is presented here. Although the approach discussed is restricted to simply connected convex regions, it nevertheless represents, to the authors’ knowledge, the rst reported successful eort in the use of this scheme for treating volume integrals in the BIE method with body force eects for generally anisotropic domains. 2. THE BIE FOR PLANE ANISOTROPIC ELASTICITY For a plane stress state of a 2-D homogeneous anisotropic elastic region, the generalized Hooke’s law can be written as (see Reference 12) 11 12 16 11 ”11 ”22 = 12 22 26 22 (1) 12 16 26 66 12 where (”11 ; ”22 ; 12 ) and (11 ; 22 ; 12 ) are strains and stresses, respectively, and constants ij are the elastic compliances of the materials. In the case of plane strain conditions, the coecients ij in equation (1) are replaced by ij , where ij = ij − i3 j3 33 (2) TREATMENT OF BODY-FORCE VOLUME INTEGRALS IN BEM 91 According to Lekhnitskii, the governing dierential equation for a two-dimensional problem of anisotropic elasticity can be obtained in terms of Airy’s stress function F as follows: 22 @4 F @4 F @4 F @4 F @4 F − 2 + (2 + ) − 2 + =0 26 12 66 16 11 @41 @21 @22 @42 @31 @2 @1 @32 (3) The corresponding characteristic equation of (3) is 22 − 226 + (212 + 66 )2 − 216 3 + 11 4 = 0 (4) It can be shown that for anisotropic materials, equation (4) has only complex roots, and they are distinct (see Reference 12). In addition, because the coecients of equation (4) are all real constants, the roots are in complex conjugate pairs. The four roots are denoted by 1 ; 1 ; 2 ; 2 In terms of generalized complex variables, i.e. z1 = (1 − y1 ) + 1 (2 − y2 ) z2 = (1 − y1 ) + 2 (2 − y2 ) (5) where (1 ; 2 ) is the eld point, and (y1 ; y2 ) is the source point, the fundamental solution for an anisotropic 2-D elastic body can be obtained as (see Reference 11) Uji (; y) = 2 Re[ri1 Aj1 log z1 + ri2 Aj2 log z2 ] (6) where Re [·] is an operator which takes the real part of a complex function. The constants rij are r11 = 11 12 + 12 − 16 1 r12 = 11 22 + 12 − 16 2 r21 = 12 1 + 22 − 26 1 r22 = 12 2 + 22 − 26 2 (7) The explicit equations for constants Aij can be found in Tan et al.10 The corresponding traction fundamental solutions of equation (6) can be obtained as Tj1 = 2n1 Re[12 Aj1 =z1 + 22 Aj2 =z2 ] − 2n2 Re [1 Aj1 =z1 + 2 Aj2 =z2 ] Tj2 = − 2n1 Re [1 Aj1 =z1 + 2 Aj2 =z2 ] (8) + 2n2 Re [Aj1 =z1 + Aj2 =z2 ] where nj are unit outward normal components at a boundary point (x1 ; x2 ). The direct formulation of BEM for anisotropic materials may be developed by following the same steps as in the isotropic cases. The resulting boundary integral equation provides a direct 92 J. J. ZHANG, C. L. TAN AND F. F. AFAGH integral relation between boundary tractions, ti , and boundary displacements, ui , and is given as (see Reference 11) I Cji ui (y) + ui (x)Tji (x; y) dSx @ I = ti (x)Uji (x; y) dSx Z@ + Xi ()Uji (; y) d ; x; y ∈ @ ∈ (9) where Xi are body forces, and @ and represent the boundary and the domain, respectively. In elastostatics, the body forces can usually be expressed as the gradient of a scalar potential , i.e. Xi = ; i ; ii = k0 (10) where k0 is a constant. The distinct advantage of the BEM as a numerical tool over domain solution methods is that only boundary discretization of the problem domain is needed. However, with the volume integral (VI), i.e. Z (VI)j = Xi ()Uji (; y) d ; y ∈ @ ; ∈ (11) due to the body forces Xi , appearing in the boundary integral equation, the notion of the method being a boundary solution technique is now violated. Because closed-form solutions to the above VI cannot be obtained generally, numerical evaluation is usually required. But the direct numerical evaluation of equation (11) needs discretization of the problem domain. Here, an attempt to transform the VI into surface or boundary integrals is made so as to restore equation (9) as a boundary-only integral equation. 3. EXACT TRANSFORMATION FOR CONVEX REGIONS This paper discusses the method for transforming the VI to boundary integrals for convex regions. The term convex region applies to a region whose tangent at every boundary point lies outside the region. For example, regions shown in Figure 1 are convex regions, while those shown in Figure 2 are not. With reference to the co-ordinate system shown in Figure 1, one can always nd a highest point A (with maximum 2 ), and a lowest point C (with minimum 2 ) along the boundary of a convex region. In some cases, these points may lie along straight lines of constant 2 , as shown in Example II, Figure 1. The points A and C divide the entire boundary of a convex region into ˙ ˙ two parts, namely the arcs ABC and CDA . For a case as in Example II, Figure 1, one can always specify any point along line A as the point A, and any point along line C as the point C. In dealing with the VI in equation (11), one has to deal with the logarithmic function of a generalized complex variable, i.e. z = (1 − y1 ) + (2 − y2 ) (12) TREATMENT OF BODY-FORCE VOLUME INTEGRALS IN BEM 93 Figure 1. Typical convex regions where = a + ib represents the characteristic roots, and b never vanishes for anisotropic materials since equation (4) has only complex roots. A two-step mapping of equation (12) follows. First, the (1 ; 2 ) plane is mapped to an intermediate (1 0 ; 2 0 ) plane through 01 = 1 − y1 02 = 2 − y2 (13) Then, the intermediate (01 ; 02 ) plane is mapped to the nal (1 ; 2 ) plane through 1 = 01 + a02 2 = b02 (14) In the (1 ; 2 ) plane, the variable z can be written as standard complex variable, z = 1 + i2 (15) Through equations (13) and (14), the original domains in Figure 1 are mapped into the new domain , as shown in Figure 3. 94 J. J. ZHANG, C. L. TAN AND F. F. AFAGH Figure 2. Typical non-convex regions From equation (13), one may note that the origin of the (01 ; 02 ) system is always on the boundary and is coincident with the source point. Furthermore, it is evident that for source points ˙ along arc ABC , the positive 01 axis cuts across the domain, while for source points along arc ˙ CDA , the negative 01 axis cuts across the domain. The mapping depicted by equation (14) has the distinctive feature of keeping the points along the 01 axis unchanged. It can be easily seen that for points 02 = 0, one always has 1 = 01 . In addition, it can be noted that the image of point A remains the highest point if b ¿ 0, and it becomes the lowest point if b ¡ 0. Conversely, the image of point C in the (1 ; 2 ) plane remains the lowest point if b ¿ 0, and it becomes the highest point if b ¡ 0. Hence, the images of points A and C remain as the dividing points for the boundary @ of the mapped region in the (1 ; 2 ) plane. Also, the image of the arc ˙ ˙ ABC remains on the left side of the region, and the image of the arc CDA remains on the right side of the region. Thus, one can conclude that if from a source point, the negative 1 0 axis cuts across the region , then the negative 1 axis will cut across the corresponding region ; and if from a source point, the positive 01 axis cuts across the region , then the positive 1 axis will cut across the corresponding region . Similarly, if the 01 axis does not cut across the region , then the 1 ˙ axis does not cut across the corresponding region . Thus, for any point along the arc ABC , the TREATMENT OF BODY-FORCE VOLUME INTEGRALS IN BEM 95 Figure 3. The images of original regions in the 1 2 plane ˙ positive 1 axis cuts across the region ; and for any point along arc CDA , the negative 1 axis cuts across the region . By denition, the complex function log(z) is log (z) = log | z | + i arg (z) (16) In computing the fundamental solution, the principal value of arg (z) is employed, i.e. − ¡ arg (z)6 (17) It is clear from equation (17) that along the negative 1 axis, arg (z) is not continuous. Thus, ˙ for any source point along ABC , the function log (z) is analytic in the problem domain; but for ˙ any source point along CDA , the function log (z) is not analytic in the problem domain because the negative 1 axis cuts across the region. However, this diculty can be overcome by simply ˙ redening the principal value of arg(z) for the source points along CDA of a convex region as 0 ¡ arg (z)62 (18) 96 J. J. ZHANG, C. L. TAN AND F. F. AFAGH In equation (18), arg(z) is now discontinuous along the positive axis, but continuous along the negative 1 axis, and log (z) becomes analytical in the problem region for source points along ˙ CDA . With the understanding that log (z) is made to be analytic in the problem domain for every source point along the boundary by redening the limits of its argument, a method for transforming the volume integral in equation (11) to boundary integrals may now be formulated. First, the following notation is introduced 1 1 [km ] = (19) 1 2 so that with the application of the summation convention, equation (5) can be simplied as zm = km (k − yk ); m; k = 1; 2 (20) Next, two functions G1 (; y) and G2 (; y) are introduced such that they satisfy the following dierential equations respectively, Gm; jj (; y) = = 1 zm 1 km (k − yk ) ; m; k = 1; 2 (21) where dierentiation is taken to be with respect to the eld point . To nd a solution to equation (21), it is important to note that the second derivative of a real function log () with respect to variable is 1=. This suggests that the following are particular solutions of equation (21), G1 = C1 z1 log (z1 ) G2 = C2 z2 log (z2 ) (22) where Cm (m = 1; 2) are constants to be determined. Applying the Laplacian operator to equation (22), and noting that log (z) is now analytic in the problem domain, gives G1; jj = C1 (k1 kj log (z1 ) + k1 kj ); j 1 = C1 j1 j1 · z1 1 G2; jj = C2 j 2 j 2 · z2 (23) where ij is the Kronecker delta function. By substituting equation (23) into equation (21), the constants Cm can be obtained as C1 = 1 ; j1 j1 C2 = 1 j 2 j 2 (24) Thus, the functions Gm ; m = 1; 2, can be nally determined as G1 = 1 1 z1 log (z1 ) = 2 z log (z1 ) 2 1 j1 j1 11 + 21 1 1 G2 = z2 log (z2 ) = 2 z log (z2 ) 2 2 j2 j2 12 + 22 (25) TREATMENT OF BODY-FORCE VOLUME INTEGRALS IN BEM 97 From functions Gm and by using a similar structure for the fundamental solution, two important functions can be constructed, i.e. Hj (; y) = 2 Re[ri1 i1 Aj1 G1 (; y) + ri2 i2 Aj 2 G2 (; y)]; j = 1; 2 (26) Applying the Laplacian operator to functions Hj , and noting equation (21), one can obtain Hj; kk = 2 Re[ri1 i1 Aj1 G1; kk (; y) + ri2 i2 Aj 2 G2; kk (; y)] 1 1 = 2 Re ri1 i1 Aj1 + ri2 i2 Aj 2 z1 z2 (27) The relation between the functions Hj ; j = 1; 2, and the fundamental solution Uji (x; y); i; j = 1; 2, may now be established. Applying the divergence operator to the fundamental solution gives Uji; i (; y) = 2 Re[ri1 Aj1 log (z1 ) + ri2 Aj 2 log (z2 )]; i 1 1 = 2 Re ri1 Aj1 k1 ki + ri2 Aj 2 k2 ki z1 z2 1 1 = 2 Re r i1 Aj1 i1 + ri2 Aj 2 i2 z1 z2 (28) Comparing equation (27) to equation (28), it is evident that Hj; ii (; y) = Uji; i (; y) (29) Two vector functions E1 = (E11 ; E21 )T , and E2 = (E12 ; E22 )T are introduced next. They are related to functions Gm (m = 1; 2) by Ei1; i = G1 (; y) = 1 z1 log (z1 ) j1 j1 Ei2; i = G2 (; y) = 1 z2 log (z2 ) j 2 j 2 (30) That is, the function Gm (m = 1; 2) is the divergence of vector Em . To nd a particular solution to the above dierential equations, it is important to note that the derivative of a real function 2 log () with respect to is 2 log () + . This suggests that the vectors Em may have the following structure, Ei1 = Di1 z 21 log (z1 ) + Fi1 z 21 Ei2 = Di2 z22 log (z2 ) + Fi2 z22 (31) where constants Dim and Fim ; i; m = 1; 2, are to be determined. By applying the divergence operator to the rst of equation (31), i.e. Ei1; i = Di1 (z12 log (z1 )); i + Fi1 (z12 ); i (32) (z12 ); i = 2z1 k1 ki = 2z1 i1 (33) and noting that (z12 log (z1 )); i = 2z1 k1 ki log (z1 ) + z1 k1 ki = 2i1 z1 log (z1 ) + i1 z1 (34) 98 J. J. ZHANG, C. L. TAN AND F. F. AFAGH one can reduce equation (32) to Ei1; i = 2Di1 i1 z1 log (z1 ) + Di1 i1 z1 + 2Fi1 z1 i1 By comparing equation (35) to equation (30), one obtains 1 2Di1 i1 = j1 j1 (35) (36) Di1 i1 + 2Fi1 i1 = 0 (37) From equations (36) and (37), the constants Di1 and Fi1 can be chosen as 1 Di1 = 4i1 j1 j1 1 1 Fi1 = − Di1 = − 2 8i1 j1 j1 (38) (39) Similarly, Di2 and Fi2 are determined to be 1 4i2 j 2 j 2 1 =− 8i2 j 2 j 2 Di2 = (40) Fi2 (41) Thus the vectors E1 (; y) and E2 (; y) are found as 1 Ei1 (; y) = z 2 (log z1 − 12 ) 4i1 j1 j1 1 1 z 2 (log z2 − 12 ) Ei2 (; y) = 4i2 j 2 j 2 2 (42) (43) By employing the vector functions E1 (; y) and E2 (; y), two other vector functions W1 = (W11 ; W12 ) and W2 = (W21 ; W22 ) can be constructed by using a structure similar to that of the fundamental solution, i.e. Wji = 2 Re[rk1 Aj1 k1 Ei1 + rk2 Aj 2 k2 Ei 2 ]; i; j; k = 1; 2 (44) Applying the divergence operator to the vectors Wj yields Wji; i = 2 Re[rk1 Aj1 k1 Ei1; i + rk2 Aj 2 k2 Ei2;i ] = 2 Re[rk1 Aj1 k1 G1 + rk2 Aj 2 k2 G2 ] (45) Comparing equation (26) to equation (45), one can immediately recognize that Wji; i (; y) = Hj (; y) (46) With the help of functions Hj ; j = 1; 2 and vectors Wi ; i = 1; 2, it is now possible to transform the volume integral equation (11) to boundary ones. Substituting equation (10) into equation (11) and applying Green’s theorem gives Z (VI)j = I ; i Uji (; y) d = @ Z Uji (x; y)ni dSx − Uji; i (; y) d ; ∈ ; y; x ∈ @ (47) TREATMENT OF BODY-FORCE VOLUME INTEGRALS IN BEM 99 where n = (n1 ; n2 ) is the outward normal of the boundary at point (x1 ; x2 ). From equation (29), the divergence of the fundamental solution Uji can be expressed in terms of the functions Hj . Thus, by applying Green’s theorem consecutively, one has I Z (VI)j = Uji (x; y)ni dSx − Hj; ii (; y) d I @ I Z = Uji (x; y)ni dSx − Hj; i (x; y)ni dSx + ; i Hj; i (; y) d I @ I @ = Uji (x; y)ni dSx − Hj; i (x; y)ni dSx @ @ I Z + ; i Hj (x; y)ni dSx − ; ii Hj (; y) d (48) @ From equation (10) and (46), the last volume integral in the above equation can be further reduced to Z Z ; ii Hj (; y) d = k0 Wji; i (; y) d I = k0 Wji (x; y) ni dSx (49) @ Thus, the volume integral in equation (11) nally can be expressed as boundary integrals, i.e. I I (VI)j = Uji (x; y)ni dSx − Hj; i (x; y)ni dSx @ @ I I + ; i Hj (x; y)ni dSx − k0 Wji (x; y)ni dSx (50) @ @ Equation (50) successfully transforms the volume integral associated with body forces in the BIE method into boundary ones. Therefore, the basic integral equation of the BEM for the 2-D anisotropic elastostatic problem can be written as I Cji ui (y) + ui (x)Tji (x; y) dSx I @ = ti (x)Uji (x; y) dSx @ I I + Uji (x; y)ni dSx − Hj; i (x; y)ni dSx @ @ I I + ; i Hj (x; y)ni dSx − (51) k0 Wji (x; y)ni dSx @ @ Equation (51) is now a boundary-only integral equation, which, at least for convex regions, overcomes the volume integral problem of the BEM for 2-D anisotropic elastostatics. The numerical solution of equation (51) may be obtained using, for example, the conventional quadratic isoparametric element formulation, the details for which may be found in e.g. Reference 6. The evaluation of the transformed surface integrals containing functions H and W due to the body force term presents no diculties as they are not singular. The method developed above can also be applied to the case of some non-convex regions. For instance, in the case of Example III, Figure 2, the principal value of the argument of a complex variable can be redened to be within the limits 3 (52) − ¡ arg (z)6 2 2 100 J. J. ZHANG, C. L. TAN AND F. F. AFAGH for all the source points along the inner semi-circular boundary. In a similar manner, log (z) can be made analytic in the problem domain for all source points along the whole boundary. It should perhaps also be remarked that, in general, a non-convex region can be judiciously divided into several convex sub-regions, or sub-regions like the case of Example III, Figure 2, on each of which the method described herein may be applied. This concept of sub-regioning is, of course, quite common in BEM analysis. 4. NUMERICAL EXAMPLES Three test problems are treated to demonstrate the validity of the exact transformation procedure developed in this paper. In all these examples, the properties of materials are described in terms of the engineering elastic constants, which are related to the elastic compliances ij as follows (see, Reference 12), 1 E1 12 21 =− = − E1 E2 1 = E2 12;1 1;12 = = E1 G12 1 = G12 12; 2 2; 12 = = E2 G12 11 = 12 22 16 66 26 (53) where Ek is the Young’s modulus in the k direction, G12 is the shear modulus in the (1 ; 2 ) plane, and jk is the Poisson’s ratio which is dened as the compressive strain in the k direction due to unit extensional strain in the j direction. Constants ij; k and i; jk are referred to as the coecients of mutual inuence of the rst and second kind, respectively. The quadratic isoparametric element formulation was employed in the present study. In the rst test problem, the volume integral associated with the body force was also evaluated for each boundary node by means of direct two-dimensional numerical integration and the values are compared with those obtained from evaluating the transformed boundary integrals. The second test problem involves the determination of the stress distribution in a rotating, plane, circular, orthotropic disc for which the exact analytical solution exists. Finally, a rotating, plane, circular, anisotropic disc with a central crack is analysed. The 8-point Gauss quadrature was used throughout for evaluating the boundary integrals in equation (51) Problem 1 A square anisotropic region of 2 m × 2 m is considered with the BEM mesh discretization as shown in Figure 4. The material constants are chosen to be as follows: E1 (GPa) E2 (GPa) G12 (GPa) 12 12;1 12;2 345 516 173 0·131 0·0 0·0 TREATMENT OF BODY-FORCE VOLUME INTEGRALS IN BEM 101 Figure 4. Problem 1—A square domain and the BEM mesh Figure 5. Body force potential for Problem 1 The body forces are taken to be X1 = C [(e−2 − e2 ) cos(1 ) + (e−1 − e1 ) sin(2 )] X2 = C [−(e−2 − e2 ) sin(1 ) − (e−1 − e1 ) cos (2 )] (54) with the corresponding potential = C [(e−2 − e2 ) sin(1 ) − (e−1 − e1 ) sin(2 )] (55) which is shown graphically in Figure 5. C is a constant, and is chosen to be 106 in the numerical calculations, assigning both X1 and X2 to have units of N (Newton). Although this may not be a physically encountered situation, the choice of this potential is to demonstrate the mathematical soundness of the approach developed. In the analysis, the boundary of the square is divided into 8 elements with a total of 16 nodes (Figure 4). Table I shows the results of the body force volume integral (VI) as obtained by the Direct Integration Method (DIM) of evaluating equation (11), and by the ‘Exact Transformation Method’ (ETM) of evaluating equation (50). In the former method, the solution domain was 102 J. J. ZHANG, C. L. TAN AND F. F. AFAGH Table I. Values of the volume integral for the square region—Problem 1 (VI)1 (mm) Node 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 (VI)2 (mm) ETM DIM ETM DIM 0·624095E − 02 0·762926E − 02 − 0·149011E − 07 − 0·762925E − 02 − 0·624094E − 02 − 0·567059E − 02 − 0·383093E − 02 − 0·567054E − 02 − 0·624099E − 02 − 0·762930E − 02 0·163917E − 07 0·762930E − 02 0·624091E − 02 0·567060E − 02 0·383092E − 02 0·567057E − 02 0·624251E − 02 0·733150E − 02 0·174622E − 09 − 0·733150E − 02 − 0·624251E − 02 − 0·566292E − 02 − 0·383115E − 02 − 0·566292E − 02 − 0·624251E − 02 − 0·733150E − 02 0·000000E + 00 0·733150E − 02 0·624251E − 02 0·566292E − 02 0·383115E − 02 0·566292E − 02 0·533362E − 02 0·514822E − 02 0·402040E − 02 0·514822E − 02 0·533363E − 02 0·629162E − 02 0·167638E − 07 − 0·629167E − 02 − 0·533360E − 02 − 0·514820E − 02 − 0·402040E − 02 − 0·514823E − 02 − 0·533361E − 02 − 0·629162E − 02 0·223517E − 07 0·629164E − 02 0·533267E − 02 0·515314E − 02 0·402029E − 02 0·515314E − 02 0·533267E − 02 0·645052E − 02 − 0·116415E − 09 − 0·645052E − 02 − 0·533268E − 02 − 0·515314E − 02 − 0·402029E − 02 − 0·515314E − 02 − 0·533267E − 02 − 0·645052E − 02 0·174622E − 09 0·645052E − 02 Figure 6. A circular rotating disc and the co-ordinate system—Problem 2 divided into 64 ‘cells’ and a 10 × 10 point Gauss quadrature was used per cell. As can be seen in the table, even with the small number of boundary elements, the results for the VI from ETM and the DIM are in very good agreement indeed. Problem 2 The stress distribution in a thin circular orthotropic disc which is rotating about its centroidal axis x3 is analysed next (Figure 6). The results of this BEM analysis can be compared to those obtained from the analytical closed-form solutions by Lekhnitskii.13 For such a solid disc with radius R and rotating at a constant angular velocity !, the stress components depend only TREATMENT OF BODY-FORCE VOLUME INTEGRALS IN BEM 103 Figure 7. BEM mesh—Problem 2 on the distance r from the centre of the rotation and, in a polar co-ordinate system, they are given as13 r = 12 !2 R2 (1 − A)(1 − (r=R)2 ) = 12 !2 R2 [(1 − A)(1 − (r=R)2 ) + 2A(r=R)2 ] (56) r = 0 A = 11 + 212 + 22 311 + 212 + 66 + 322 where the constant A is an invariant with the transformation of the co-ordinate system. For the purpose of analysis, the material used in this example is a glass–epoxy composite with the following elastic constants: E1 (GPa) E2 (GPa) G12 (GPa) 12 12;1 12;2 48·26 17·24 6·89 0·29 0·0 0·0 The body forces due to the rotation are X1 = x1 !2 (57) X2 = x2 ! (58) 2 Therefore, the body force potential is = 12 !2 (x12 + x22 ) (59) Since constant A is an invariant with the transformation of the co-ordinate system, for simplicity, the material principal axes are also taken as x1 and x2 here. Due to the symmetry with respect to 104 J. J. ZHANG, C. L. TAN AND F. F. AFAGH Figure 8. Variation of radial and tangential stresses in the rotating disc with radial distance—Problem 2 Table II. Percentage errors for radial and tangential stresses—Problem r=R 0·0000 0·0625 0·1250 0·1875 0·2500 0·3750 0·5000 % error of r % error of r=R % error of r % error of 0·2058 1·7651 1·1089 1·5382 0·4369 1·9147 − 0·0156 0·2229 0·4814 0·2103 0·7569 − 0·4029 0·6843 − 0·5445 0·6250 0·7500 0·8125 0·8750 0·9375 1·0000 2·6281 0·9407 2·9434 1·6480 3·5840 0·6703 0·4282 − 0·9390 1·9364 − 2·7796 0·2496 0·9766 x1 and x2 axes, only the rst quarter of the disc need be considered with the appropriate boundary condition. The corresponding boundary discretization is shown in Figure 7, where 20 boundary elements with a total of 40 nodes are employed. The BEM results of radial and tangential stresses, normalized by !2 R2 , together with the exact solution from Lekhnitskii13 are shown in Figure 8, where it can be seen that the agreement between the two sets of results is excellent. The corresponding percentage errors are listed in Table II. Problem 3 Next, a thin circular anisotropic rotating disc with a central crack as shown in Figure 9 is investigated. It is assumed that the crack can have any arbitrary orientation with respect to the material principal axis. The origin of the co-ordinate system is taken to be at the centre of the disc with the x1 axis aligned along the plane of the crack. The stress intensity factors KI and KII are to be determined. As the closed-form solution for this problem is not readily available, Bueckner’s superposition method14 was used to obtain an independent set of results for comparison with those obtained with the present approach. Using this method, the stress intensity factor for the crack in the TREATMENT OF BODY-FORCE VOLUME INTEGRALS IN BEM 105 Figure 9. Circular rotating disc with a central crack—Problem 3 Figure 10. Circular non-rotating disc with normal tractions along crack faces rotating disc is equal to that for a crack in an identical, but non-rotating disc for which the crack faces are subjected to normal tractions equal to the circumferential stress acting on the prospective crack plane in the uncracked rotating disc (Figure 10). This circumferential stress distribution is readily available in the literature.13 In this manner, the original crack problem with body forces is converted to a simple crack problem without any body forces. This problem can now be solved using any standard BEM code, the solution of which may thus be used for verication of the results from the exact volume to boundary integral transformation approach. Due to the material anisotropy, the problem is not symmetric about either the x1 or x2 axis in general. Hence the boundary of the whole domain needs to be discretized. First, the domain is divided into two sub-regions along the plane of the crack; a typical mesh used is shown in 106 J. J. ZHANG, C. L. TAN AND F. F. AFAGH Figure 11. Typical BEM mesh—Problem 3 Figure 11. Traction-singular quarter-point boundary elements, see e.g. Reference 15, were employed adjacent to the crack tips and stress intensity factors were computed directly from the computed nodal traction coecients there, as is well established in the literature for BEM fracture mechanics analysis.15; 16 The material used in this example is the T300 / N5208 composite which has the following elastic constants:8 E1 (GPa) E2 (GPa) G12 (GPa) 12 12;1 12; 2 131 13 6·4 0·038 0·0 0·0 Several material property orientations, as measured by the angle (see Figure 9), namely, = 0◦ , 15◦ , 30◦ , 45◦ , 60◦ , 75◦ and 90◦ , as well as normalized crack length, a=R = 0·125, 0·25, 0·375, 0·5, 0·625, 0·75, 0·875 were considered. Figure 12 shows the results for the non-dimensionalized stress intensity factor KI =Ka at the crack tips as functions of crack length, where √ Ka = 0 a (60) and 0 is the tangential stress at r = 0 in the absence of the crack. The percentage deviations between the values obtained by Bueckner’s method and the present method are generally less than 1 per cent, the greatest discrepancy being only 1·2 per cent. The corresponding results for KII =Ka as a function of crack length a=R are shown in Figure 13 for the dierent orientations of the material TREATMENT OF BODY-FORCE VOLUME INTEGRALS IN BEM 107 √ Figure 12. Variation of normalized stress intensity factor KI =Ka with relative crack length a=R; Ka = 0 a axes. The results of KII for = 0 and 90◦ of the material axes are zero for any crack length, as is to be expected. It is evident again that there is excellent agreement between the results obtained using the BEM formulation presented in this paper and Bueckner’s superposition approach. The relative dierences for KII are however larger than the ones for KI . This is attributed to the fact that the KII are numerically much smaller in magnitude compared to KI ; hence their values are more easily inuenced by the calculation round-o errors in numerical analysis. 108 J. J. ZHANG, C. L. TAN AND F. F. AFAGH √ Figure 13. Variation of normalized stress intensity factor KII =Ka with relative crack length a=R; Ka = 0 a 5. CONCLUSION In the direct BEM for elastostatics, the presence of body force eects results in a volume integral term in the boundary integral equations. In isotropic elasticity, this volume integral can be transformed exactly into boundary ones. This enables the BEM to maintain itself as a truly boundary solution computational technique. Hitherto, however, the same cannot be said of the BEM for general anisotropic elasticity involving body forces. In this study, the exact transformation of the volume integral to boundary integrals, for two dimensional anisotropic regions, has been presented. The validity of the analysis has been demon- TREATMENT OF BODY-FORCE VOLUME INTEGRALS IN BEM 109 strated by three examples. Although this exact transformation has been successfully implemented in the BEM for convex geometric regions, the study nevertheless suggests the feasibility of the approach to treat body forces and, in principle, thermoelastic eects. Research to remove this limitation on the geometry of the anisotropic region analysed and to extend the procedures to the thermoelastic problem is currently in progress. REFERENCES 1. G. S. Gipson and C. V. Camp, ‘Eective use of Monte Carlo quadrature for body force integrals occurring in integral form of elastostatics’, Proc. 7th Int. Conf. on Boundary Elements, 1985, pp. 17–26. 2. C. V. Camp and G. S. Gipson, Boundary Element Analysis of Nonhomogeneous Biharmonic Phenomena, Springer, Berlin, 1992. 3. F. J. Rizzo and D. J. Shippy, ‘An advanced boundary integral equation method for three dimensional thermoelasticity’, Int. j. numer. methods eng., 11, 1753–1768 (1977). 4. D. Nardini and C. A. Brebbia, ‘A new approach to free vibration using boundary elements’, in C. A. Brebbia (ed.), Boundary Element Methods in Engineering, springer, Berlin, 1982, pp. 312–326. 5. A. J. Nowak, ‘Temperature elds in domains with heat sources using boundary-only formulation’, Proc. 10th BEM Conf., Springer, Berlin 1988, pp. 233–247. 6. J. O. Watson, ‘Advanced implementation of the boundary element method for two and three-dimensional elastostatics’, in P. K. Banerjee and R. Buttereld (eds.), Developments in BEM, Applied Science Publisher, 1979. 7. P. K. Banerjee and R. Buttereld, Boundary Element Methods in Engineering Science, McGraw-Hill, London, 1981. 8. A. Deb and P. K. Banerjee, ‘BEM for general anisotropic 2D elasticity using particular integrals’, Comm. appl. numer. methods, 6, 111–119 (1990). 9. A. Deb, D. P. Henry and R. B. Wilson, ‘Alternate BEM formulations for 2- and 3-D anisotropic thermoelasticity’, Int. J. Solids Struct., 27, 1721–1738 (1991). 10. C. L. Tan, Y. L. Gao and F. F. Afagh, ‘Anisotropic stress analysis of inclusion problems using the boundary integral equation method’, J. Strain Anal., 27, 67–76 (1992). 11. T. A. Cruse, Boundary Element Analysis in Computational Fracture Mechanics, Kluwer Academic, Dordrecht, 1988. 12. S. G. Lekhnitskii, Theory of Elasticity of an Anisotropic Body, Mir, Moscow, 1981. 13. S. G. Lekhnitskii, Anisotropic Plates, Gordon and Breach, London, 1968. 14. H. F. Bueckner and N. Y. Schenectady, ‘The propagation of cracks and the energy of elastic deformation’, Trans. ASME, 80E, 1225–1230 (1958). 15. J. Martinez and J. Dominguez, ‘On the use of quarter-point boundary elements for stress intensity factor computations’, Int. j. numer. methods eng., 20, 1941–1950 (1984). 16. C. L. Tan and Y. L. Gao, ‘Boundary element analysis of plane anisotropic bodies with stress concentration and cracks’, Composite Struct., 20, 17–28 (1992). 17. Q. Ao, ‘Development of boundary methods of solution of anisotropic problems’, Ph.D. Dissertation, University of Kentucky, 1994. 18. Q. Ao and D. J. Shippy, ‘BEM treatment of body forces in plane orthotropic elastostatics’, J. Eng. Mech., 121, 1130–1135 (1995).