The solution of an inhomogeneity in a finite plane region and its application to composite materials

Composites Science and Technology 60 (2000) 75±82

The solution of an inhomogeneity in a ®nite plane region and its application to composite materials $

J. Wang a, J.H. Andreasen b, B.L. Karihaloo c,* a

Department of Mechanics and Engineering Science, Peking University, Beijing 100871, People's Republic of China b Institute of Mechanical Engineering, Aalborg University, Pon 101, DK-9220 Aalborg East, Denmark c Division of Civil Engineering, University of Wales, Cardi€, Queen's Buildings, PO Box 686, Newport Road, Cardi€ CF2 3TB, UK Received 3 August 1998; received in revised form 24 May 1999; accepted 30 July 1999

Abstract A method is described for determining the elastic ®eld in a ®nite isotropic region containing an isotropic circular inhomogeneity whose elastic constants di€er from those of the rest of the region. The method is based on a combination of Muskhelishvili complex potentials and boundary collocation. The method is particularly useful for solving problems of practical interest in the ®eld of composite materials. Two examples are given to illustrate this. # 2000 Elsevier Science Ltd. All rights reserved. Keywords: Circular inhomogeneity; Finite region; Boundary collocation; Composite materials; Transverse and shear moduli

1. Introduction In 1957, Eshelby proposed an elegant method for determining the elastic ®eld in an in®nite body containing a ®nite inhomogeneity whose elastic constants di€er from those of the rest of the body [1]. He showed how the ®eld inside and outside the inhomogeneity could be constructed by using the notion of an equivalent inclusion whose elastic constants are the same as those of the rest of the body but which is subjected to unknown eigenstrains. Eshelby also suggested a method for constructing the elastic ®eld if the ®nite inhomogeneity is contained in a body of ®nite dimensions by using the notion of image ®elds that cancel tractions on the traction-free surfaces of the body [2]. The construction of these image ®elds is a formidable task. In two-dimensions, the Eshelby formalism can of course be expressed equally elegantly in terms of Muskhelishvili [3] complex potentials (see, e.g. Ref. 4). We shall exploit this duality between the Eshelby formalism and complex potentials on the one hand and on the other, have recourse to the $ This work was completed while JW and BLK were at Aalborg University, Denmark. * Corresponding author. Tel.: +44-1222-874-597; fax: +44-1222874-597. E-mail address: karihaloob@cardi€.ac.uk (B.L. Karihaloo).

boundary collocation method [5] which is known to solve two-dimensional elastostatic problems for ®nite domains with a high degree of accuracy. A judicious choice of the complex potentials allows satisfaction of the traction and displacement continuity conditions on the interface between the inhomogeneity and surrounding body, simultaneously with the prescribed traction and displacement conditions on the outer boundary of the body. This technique was used ®rst by Sendeckyj [6] who applied it to an inhomogeneity in an in®nite plane and who anticipated its application to a ®nite plane without actually developing the corresponding potentials. In this manner, we shall develop a method for the solution of the problems of an inhomogeneity in a ®nite-plane region. We shall demonstrate the development of the solution method on a circular inhomogeneity, but the method can be easily generalised to any inhomogeneity which can be mapped on to a circle, as has been shown by Sendeckyj [6]. We shall then illustrate the use of the method on two examples commonly encountered in the design of composite materials. The ®rst example concerns the e€ect of the relative values of moduli of the inhomogeneity and the surrounding material on the stress ®eld in a rectangular region under unidirectional loading. The results of this example are useful for the prediction of the post-impact compressive strength of a ®bre-reinforced quasi-isotropic

0266-3538/00/$ - see front matter # 2000 Elsevier Science Ltd. All rights reserved. PII: S0266-3538(99)00103-7

76

J. Wang et al. / Composites Science and Technology 60 (2000) 75±82

composite laminate that has su€ered a low-velocity impact resulting in a circular damage area at its centre. It is of great interest to know how this damaged area a€ects the in-plane compressive strength of the composite laminate as it is barely visible on the surface, thus posing a hidden menace to the in-plane compressive strength of the laminate (e.g. Refs. 7±10). The second example concerns the prediction of transverse and shear moduli and of the stress ®eld in composites reinforced by continuous circular ®bres distributed in a square or hexagonal pattern. The theoretical prediction of the moduli using the general method of solution are shown to be in very good agreement with experimental results. 2. Method of solution Consider a ®nite plane region containing a circular inhomogeneity at its centre, as shown in Fig. 1. The shear moduli and Kolosov constants of the region outside and inside the inhomogeneity are 1 and 1 , and 2 and 2 , respectively. In the sequel, index 1 will refer to the region outside the inhomogeneity and 2 to that inside, as indicated in Fig. 1. For the plane problem shown in Fig. 1, the stresses xx , yy and xy , can be expressed in terms of the complex potentials …z† and …z† through [3] xx ‡ yy ˆ 2f0 …z† ‡ 0 …z†g; yy ÿ xx ‡ 2ixy ˆ 2fz00 …z† ‡

0

…z†g;

…1†

where a prime denotes di€erentiation with respect to z ˆ x ‡ iy and an overbar denotes the complex conjugate. The displacements u and v in the x- and y-directions, respectively, can be written as

2…u ‡ iv† ˆ …z† ÿ z0 …z† ÿ …z†;

…2†

where  is the shear modulus and  is the Kolosov constant related to the Poisson ratio  through  ˆ 3 ÿ 4 for plane strain and  ˆ …3 ÿ †=…1 ÿ † for plane stress. It is shown in Ref. 11 that the solution of twodimensional problems in the xy-plane for isotropic media as shown in Fig. 1 can be readily applied to transversely isotropic (in the xy-plane) media through the introduction of a parameter similar to the Kolosov constant for isotropic media. The complex potentials …z† and …z†, like  and , will be di€erent inside and outside the inhomogeneity. The complex potentials for the region outside the inhomogeneity can be expressed in Laurent series for any general loading on the ®nite outer boundary  z n  X Rn an ‡An ; 1 …z† ˆ z R nˆ1 …3†  z n  X Rn bn ‡Bn ‡ B0 ; 1 …z† ˆ z R nˆ1 whereas for the region occupied by the inhomogeneity, the potentials can be expressed in Taylor series X  z n cn ; 2 …z† ˆ R nˆ1 …4† X  z n dn : 2 …z† ˆ R nˆ0 The choice of the series representations (3) and (4) ensures that the potentials are holomorphic in their respective regions. The constants B0 and d0 permit rigid body displacements of the inhomogeneity and the surrounding matrix. Additional constant terms in the potentials would be redundant. In the above series representations, the four groups of complex coecients an , bn , cn , and dn are eliminated analytically by satisfying the continuity conditions of displacement and traction across the boundary of the inhomogeneity. This leaves only the coecients An and Bn to be determined from the boundary conditions on the outer boundary that can be of arbitrary geometry and be subjected to arbitrary loading. The radial and tangential displacements ur and u corresponding to (2) are 2…ur ‡ iu † ˆ eÿi ……z† ÿ z0 …z† ÿ …z††;

…5†

and the net force increment as z traverses an arc is Py ‡ iPx ˆ ÿ…z† ÿ z0 …z† ÿ …z†: Fig. 1. A ®nite region containing a circular inhomogeneity of radius R.

…6†

The integral representation of the tractions (6) can also be derived from the stresses (1). Constant terms

J. Wang et al. / Composites Science and Technology 60 (2000) 75±82

which represent rigid body displacements and have no bearing on stresses have been ignored. Continuity of (6) ensures continuity of the tractions. Substituting (3) and (4) in turn into (5) and equating the resulting expressions for z ˆ Rei gives the condition for continuity of displacements across the boundary of the inhomogeneity. Similarly, substituting (3) and (4) in turn into (6) and equating the resulting expressions for z ˆ Rei ensures continuity of tractions across the inhomogeneity. Collecting like terms in the powers of k…k ˆ 0; 1; 2; . . .† in eik in the two continuity conditions gives the following sets of equations: 0 ˆ 2A2 ‡ B0 ÿ ÿ…2c2 ‡ d0 †; 0 ˆ A2 ‡ b 2 ÿ c2 ;

…7†

0 ˆ 1 A2 ÿ b 2 ÿ ÿ2 c2 ; 0 ˆ …A1 ‡ A 1 † ‡ b1 ÿ …c1 ‡ c1 †; 0 ˆ …A1 ÿ 1 A 1 † ‡ b1 ÿ ÿ…c1 ÿ 2 c1 †:

…8†

and for n51:

…9† where ÿ ˆ 2 =1 , and an overbar again denotes complex conjugate. Solving (7)±(9) for an , bn , cn , and dn in terms of An and Bn , and substituting the resulting expressions into (3) and (4) gives the following potentials for the region outside the inhomogeneity

R  A1 ‡ A 1 …1 ‡ l4 † ‡ A 1 …l1 ÿ l2 † z  n M ‡2 X R A n …l1 ‡ n…n ÿ 2†l2 † ‡ z nˆ1 " #   M  z n X R …n‡2†  ‡ Bn ‡nBn l2 ; R z nˆ0 ˆÿ



and for the inhomogeneity itself

2 …z†

…11†

ˆ ÿ2A2 …1 ‡ l2 † ‡ …2A2 ‡ B0 †ÿ M‡2 M  z nÿ2 X  z n X nAn …l1 ÿ l2 † ‡ Bn …1 ‡ l2 † : ÿ R R nˆ2 nˆ1

The constants l1;...;4 are given by ‡ ; 1ÿ ÿ ; l2 ˆ 1‡ 1‡ ; l2 ˆ 1ÿ 1‡ ; l4 ˆ 2 ÿ  ÿ 1

l1 ˆ

…12†

ÿ…1 ‡ 1† ÿ …2 ‡ 1† ; ÿ…1 ‡ 1† ‡ …2 ‡ 1† ÿ…k1 ÿ 1† ÿ …2 ÿ 1† : ˆ ÿ…1 ‡ 1† ‡ …2 ‡ 1†

ˆ

0 ˆ 1 An‡2 ‡ na n ÿ b n‡2 ÿ ÿ2 cn‡2 ;

1 …z†

2 …z† ˆ

where  and are the Dundurs' parameters [12]

0 ˆ …n ‡ 2†An‡2 ‡ Bn ‡ a n ÿ ……n ‡ 2†cn‡2 ‡ dn †; 0 ˆ …n ‡ 2†An‡2 ‡ Bn ÿ 1 a n ÿ ÿ……n ‡ 2†cn‡2 ‡ dn †; 0 ˆ An‡2 ÿ na n ‡ b n‡2 ÿ cn‡2 ;

z 1 …z† ˆ ÿ2A 2 l2 ÿ A 1 l2 R " #  nÿ2 M ‡2  z n X R  nAn l2 ‡An ‡ z R nˆ1   n M X R B n l2 ; ‡ z nˆ1

z 1 …A1 ‡ A 1 †l4 ÿ …A1 ÿ A 1 †l3 2 R M ‡2  z n X An …1 ‡ l1 † ; ‡ R nˆ2

77

…10†

…13†

The use of a truncated series in the potentials (10) and (11) is consistent with the equations determining the coecients (9) and does not violate any conditions on the boundary of the inhomogeneity. Substituting (10) and (11) into (1) and (2) gives the complete stress and displacement ®elds outside and inside the inhomogeneity, respectively. The coecients An and Bn are to be determined from the boundary conditions on the outer boundary. The choice of the number of terms M is dictated by the desired accuracy. In the limiting case, when the inhomogeneity degenerates into a hole (2 ˆ 0) the constants simplify to l1 ˆ l2 ˆ ÿ1, and l4 ˆ 0. At the other extreme, when the inhomogeneity is perfectly rigid …2 ! 1† the constants simplify to l1 ˆ 1 ; l2 ˆ 1=1 , and l4 ˆ ÿ…1 ÿ 1†=2. In both extreme cases, 2 and 2 for the inhomogeneity are meaningless. The displacements and force resultants outside the inhomogeneity are obtained by substituting (10) into (2) and (6), respectively z  21 …u ‡ iv† ˆ …A1 ÿ 1 A 1 †l2 ‡ …A1 ‡ A 1 †…1 ‡ l4 † R R ‡A1 …l1 ÿ l2 †† ÿ 21 A 2 l2 ÿ B 0 z

78

J. Wang et al. / Composites Science and Technology 60 (2000) 75±82

       n  z n zz R ‡ An 1 ‡ n…n ÿ 2†l2 2 ÿ 1 ÿ l1 R R z nˆ1 "  nÿ2  n #) R z z ÿ ‡nA n 1 l2 z z R (   n‡2 M X zz R nBn l2 2 ÿ 1 ‡ R z nˆ1  n  n   R z ÿ ; ‡B n 1 l2 z R …14† and M ‡2 X

z  Py ‡ iPx ˆ l2 …A1 ‡ A 1 † ‡ …A1 ‡ A 1 †…1 ‡ l4 † R R ‡A 1 …l1 ÿ l2 † ‡ 2A2 l2 ÿ B0 ( "z   # M‡2 X R nÿ2 z  z n ‡ ÿnAn l2 ‡ z z R nˆ1   n  n    zz R z  ‡An n…n ÿ 2†l2 2 ÿ 1 ÿ l1 ÿ R z R   n    M  X R z n ÿ Bn l2 ‡ z R nˆ1   n‡2 ) zz R : ÿnB n l2 2 ÿ 1 R z …15†

In the illustrative examples to follow, we shall be particularly interested in a general symmetric or skewsymmetric loading. For the former the following symmetry conditions prevail xy ˆ u ˆ 0; along x ˆ 0;

3. Illustrative examples 3.1. E€ect of moduli mismatch on stress distribution As the ®rst example, consider a rectangular region 2b2c containing a circular inhomogeneity of diameter 2a at its centre. A Cartesian coordinate system is set up for this region such that the origin is at the centre of the circular inhomogeneity. The x- and y-axis are parallel with the sides of lengths 2b and 2c, respectively. Suppose that the region is subjected to a uniform normal stress yy ˆ  0 on the two sides of length 2b. From a strength point of view, the most relevant stresses under uniaxial loading are the normal stress yy along the section ÿ b4x4b, y ˆ 0 (or by symmetry along 04x4b, y ˆ 0) and the interfacial stresses along the interface between the matrix and the inhomogeneity. In this example, we only show the stresses xx and yy along the section 04x4b, y ˆ 0. The interfacial stresses will be examined in the second example. Figs. 2 and 3 show the distribution of the major stress yy along the section 04x4b, y ˆ 0 for two values of the ratio ÿ ˆ 2 =1 . In Fig. 2, this ratio is 0.1, so that the inhomogeneity is softer than the surrounding material, whereas in Fig. 3, ÿ is 10, so that the inhomogeneity is much sti€er than the matrix. The Poisson ratios are assumed to be equal, 1 ˆ 2 ˆ 0:3. The ratio c=b is 1.5. Also shown in Fig. 2 is the distribution of yy = 0 , when the inhomogeneity degenerates into a hole …ÿ ˆ 0†. The stress ®eld for an in®nitely long plate with a central hole (i.e. c=a ! 1) can be found in the book by Savin [13]. For large values of c=a, the present method gives results which are very close to those given by Savin. For example, when c=a ˆ 10:0 and b=a ˆ 2:0, the values of

…16†

xy ˆ v ˆ 0; along y ˆ 0; whereas for the latter, we have the following skew-symmetry conditions xx ˆ v ˆ 0; along x ˆ 0;

…17†

yy ˆ u ˆ 0; along y ˆ 0: For both the symmetric (16) and skew-symmetric (17) loading the even order terms in z must vanish, so that An ˆ Bn ˆ 0 for n even. For the symmetric case, the remaining coecients are real for odd n, whereas for the skew-symmetric case they are imaginary. In the latter case, the imaginary coecient A1 gives a rigid body rotation. In the limit as 2 ! 0, the expressions (14) and (15) reduce to those in Ref. 5.

Fig. 2. Variation of yy = 0 along the section 4x4b; y ˆ 0 for a rectangular region (c/b=1.5) containing a circular inhomogeneity. The solid lines show the results for ÿ ˆ 0:10, and the dashed lines show the results for 2 =1 ˆ 0. The (horizontal) dotted line shows the (uniform) stress in the inhomogeneity when both b/a and c/a tend to in®nity.

J. Wang et al. / Composites Science and Technology 60 (2000) 75±82

79

Fig. 3. Variation of yy = 0 along the section 04x4b; y ˆ 0 for a rectangular region containing a circular inhomogeneity with ÿ ˆ 10.

yy = 0 at the points x ˆ a, y ˆ 0 and x ˆ b, y ˆ 0 are 4.34 and 0.75 (Fig. 2), respectively, whereas the values given by Savin for an in®nitely long plate are 4.32 and 0.73, respectively [13]. When b=a and c=a tend to in®nity, we know from the Eshelby formalism [1] that the stresses in the inhomogeneity should be uniform, if the remote applied stress is uniform. The uniform stress yy inside the inhomogeneity normalized by  0 for the two ratios of the shear moduli is shown in the ®gures by the dotted lines. The stress ®eld in the inhomogeneity remains almost uniform despite the boundary e€ect. For example, for the most severe case when the inhomogeneity occupies half the width of the region at y ˆ 0 (i.e. a=b ˆ 0:5) and ÿ ˆ 10:0, yy = 0 increases from 1.25 at the centre of the inhomogeneity to 1.32 at its edge x ˆ a, y ˆ 0. The corresponding value (uniform) for an in®nite matrix is 1.41. Fig. 2 shows that even a very soft inhomogeneity can signi®cantly reduce the stress concentration in the matrix. For instance when ÿ ˆ 0:10 and a=b ˆ 0:3, the stress concentration factor in the matrix at the point x ˆ a, y ˆ 0 is 2.76, i.e. a reduction of 20% in comparison with the stress concentration factor 3.46 for a hole of the same radius. Figs. 4 and 5 show the variations of the normalized minor stress xx = 0 along the section 04x4b, y ˆ 0 for the two values of ÿ. The sign of xx = 0 changes when ÿ increases from 0.1 to 10. As seen from Figs. 4 and 5, xx takes a substantial non-zero maximum value. It appears that a small inhomogeneity or cutout gives rise to a large xx along the section 04x4b, y ˆ 0 in the material than does a large inhomogeneity or cutout. 3.2. Transverse Young and shear moduli and stresses in ®bre-reinforced composites 3.2.1. Transverse Young and shear moduli As the second illustrative example, we shall use the above method to predict the transverse and shear moduli and the stress ®eld in a unidirectional ®bre-reinforced composite by assuming that the circular continuous ®bres are doubly periodically distributed in the matrix. Two distributions of ®bre are considered Ð the square and the hexagonal con®gurations. Figs. 6(a) and 7(a) show the

Fig. 4. Variation of xx = 0 along the section y ˆ 0; 04x4b. The solid lines show the results for ÿ ˆ 0:10. The dashed line above the datum shows the uniform stress in the inhomogeneity, when both b/a and c/a tend to in®nity.

Fig. 5. Variation of xx = 0 along the section y ˆ 0; 04x4b for ÿ ˆ 10. The dashed line below the datum shows the uniform stress in the inhomogeneity, when both b/a and c/a tend to in®nity.

normal cross-sections of the composite for these con®gurations of ®bre. The ®bres are directed normal to the xy-plane (i.e. in the z-direction). For the prediction of the transverse Young modulus, it is assumed that the composite is subjected to a unidirectional stress in the y-direction with an average 0 , whereas for the determination of its shear value yy modulus, it is assumed to be subjected to a shear stress in 0 . Due to periodithe xy-plane with an average value xy city, only the two unit cells shown in Figs. 6(b) and 7(b) need to be studied. Moreover, we need only consider a quarter of each unit cell denoted ACDE in the ®gures. Under transverse loading, the edge DE of the unit cell in Fig. 6(b) will remain straight, but that in Fig. 7(b) will be distorted. Under shear loading all the edges of the unit cells will be distorted. The boundary conditions for the

80

J. Wang et al. / Composites Science and Technology 60 (2000) 75±82

two unit cells under transverse and shear loading are listed in Tables 1 and 2, respectively. These boundary conditions are prescribed in terms of the overall equilibrium conditions of the unit cells and the symmetry of the deformation and loading states. The overall Young and shear moduli of the composite are calculated in such a way that they are the parameters that relate the corresponding average stress and strain. Both the average stress and strain are calculated over the outer boundary of the unit cell. The transverse Young modulus Ec of the composite is calculated using the following relation, assuming a plain-strain deformation state,

Ec ˆ

0 1 ‡ …2zy yy †=…Ez "yy †

Fig. 7. A hexagonal con®guration of ®bres in a composite.

Table 1 Boundary conditions for the unit cell under uniform transverse loading Fibre con®guration

On side CD

On side DE

Square

u=constant [Px]jD C =0  xy=0 u=constant A [Px]jD C ÿ[Px]jE =0  xy=0

u=constant [Py]jED= 0yyb  xy=0 u(x)ÿu(0)=u(b)ÿu(bÿx) v(x)ÿv(0)=v(b)ÿv(bÿx) [Py]jED= 0yyb  yy(x)= yy(bÿx)  xy(x)= xy(bÿx)

Hexagonal

;

…18†

where "yy is given by "yy ˆ


: c

…19†

 is the average vertical displacement of the Here,
edge DE. In Eq. (18), the longitudinal Young modulus Ez and the Poisson ratio zy are estimated using the law of mixtures. The shear modulus is calculated using the following relation, Gc ˆ

Fig. 6. A square con®guration of ®bres in a composite.

0 yy ="yy

0 xy ;
u=c ‡
v=b

…20†

where
u and
v are the displacements of point D in the x- and y-directions, respectively. b and c are equal for the square con®guration. For a glass/epoxy unidirectional ®bre-reinforced composite studied in Ref. 14, the variation of the transverse Young modulus Ec with the volume fraction of ®bre Vf …Vf ˆ a2 =…4b2 †, square con®guration; p Vf ˆ a2 =…4bc†, hexagonal con®guration …c=b ˆ 3=2†† is shown in Fig. 8 for the two con®gurations. The Young moduli and the Poisson ratios of the ®bre and the matrix are 73.1 GPa and 3.45 GPa, and 0.22 and 0.35, respectively. It is found that for the same volume fraction of ®bre, the hexagonal con®guration of ®bres gives a lower value of Ec than does the square con®guration. Also shown in Fig. 8 are the experimental results given in Ref. 15. It is seen that the experimental data fall by and large around an area bounded by the two theoretical curves. The upper bound of this area is the curve predicted using the square con®guration, and the lower bound that predicted using the hexagonal con®guration. Fig. 9 shows the variation of the transverse shear modulus of a carbon/epoxy composite with the volume Table 2 Boundary conditions for the unit cell under uniform shear loadinga Fibre con®guration On side CD

On side DE

Square

u=constant [Px]jED= 0xyb  yy=0 u(x)ÿu(0)=u(b)ÿu(bÿx) v(x)ÿv(0)=v(b)ÿv(bÿx) [Px]jED= 0xyb  yy(x)= yy(bÿx)  xy(x)= xy(bÿx)

Hexagonal

a

v=constant 0 [Py]jD C =sxyc  xx=0 v=constant 1 E 0 [Py]jD C +2[Py]jD=  xyc  xx=0

Rigid body rotation is prevented by setting A1 equal to zero in Eq. (14).

J. Wang et al. / Composites Science and Technology 60 (2000) 75±82

81

0 Fig. 10. Variation of rr =yy . The solid lines are for the square con®guration and the dashed lines for the hexagonal con®guration.

Fig. 8. Variation of the transverse Young modulus with the volume fraction of ®bre for a glass/epoxy composite.

0 Fig. 11. Variation of r =yy . The solid lines are for the square con®guration and the dashed lines for the hexagonal con®guration.

Fig. 9. Variation of the shear modulus with the volume fraction of

fraction of ®bre. The transverse shear modulus of this composite was also predicted in Ref. 14 using the boundary element method. The shear moduli and the Poisson ratios of the ®bre and the matrix are 6.60 GPa and 1.25 GPa, and 0.348 and 0.366, respectively. Also shown in this ®gure are the experimental results in Ref. 16. It is seen that the experimental data also fall by and large around the two theoretical curves.

0 Fig. 12. Variation of  =yy at the interface on the matrix side. The solid lines are for the square con®guration and the dashed lines for the hexagonal con®guration.

3.2.2. Stress ®eld Using the method developed above, the stress and displacement ®elds are accurately determined, when the boundaries CD and DE are discretized into, say, 12 equally-spaced collocation points. In the following, we only show the stresses at the interface between the ®bre and matrix when the composite is subjected to a unidirec0 . tional stress in the y-direction with an average value yy Figs. 10 and 11 show the variations of the normalized 0 0 and shear stress r =yy at the ®breradial stress rr =yy matrix interface with the volume fraction of ®bre. Each ®gure shows the results for the square and hexagonal con®gurations. This information on stresses is extremely

useful to the designers in investigating potential debond sites between the ®bre and matrix. For a low volume fraction of ®bre, the interfacial radial stresses for the two ®bre con®gurations are very close in magnitude, indicating an insigni®cant e€ect of the interaction among the ®bres. For the square con®guration, the radial stress attains its maximum at  ˆ 90 for all volume fractions of ®bre. The radial stress varies more smoothly when the ®bres are distributed in a hexagonal pattern than when they are distributed in a 0 appears square pattern. The maximum value of rr =yy  at  ˆ 90 when the ®bres are dilutely distributed. The volume fraction of ®bre does not signi®cantly a€ect the

82

J. Wang et al. / Composites Science and Technology 60 (2000) 75±82

magnitude and trend of the interfacial shear stress under the considered loading condition. This stress attains its maximum when  is about 45 , with the exception of a large fraction of ®bre in the square con®guration. Fig. 0 12 shows the normalized circumferential stress  =yy i  in the matrix at the interface z ˆ ae … 2 ‰0;2 Š† for the two con®gurations. It exhibits a similar trend to that of 0 . rr =yy 4. Discussion and conclusions A method for the solution of two-dimensional elastostatic problems for ®nite domains containing a circular inhomogeneity was developed using complex potentials and boundary collocation. The complex potentials inside and outside the inhomogeneity were chosen in a series form whose unknown complex coecients were determined by satisfying the traction and displacement continuity conditions at the interface and the traction and/or displacement conditions on the outer boundary of the ®nite domain. Although the method was developed for a circular inhomogeneity, there is no reason why it could not be easily extended to any inhomogeneity which can be mapped on to a circle. The loading on, and geometry of, the outer boundary can be arbitrary. The accuracy of the solution can be controlled by an appropriate choice of the number of terms in the series expansions, i.e. the number of boundary collocation points. The application of the method of solution was illustrated on two examples borrowed from the ®eld of composite materials. Highly accurate stress ®eld and overall moduli were obtained with the least amount of computational e€ort. Acknowledgements The support of the Danish Research Academy to JW through the award of a visiting post-doctoral fellowship (DANVIS) is gratefully acknowledged.

References [1] Eshelby JD. The determination of the elastic ®eld of an ellipsoidal inclusion, and related problems. Proc Roy Soc London 1957;A241:376±96. [2] Eshelby JD. Elastic inclusions and inhomogeneities. In: Sneddon IN, Hill R, editors. Progress in solid mechanics, vol. 2. Amsterdam: North-Holland, 1961. p. 89±140. [3] Muskhelishvili NI. Some basic problems of the mathematical theory of elasticity. Leyden: Noordho€ International Publishing, 1954 [English Edition, translated by JRM Radok]. [4] Karihaloo BL, Andreasen JH. Mechanics of transformation toughening and related topics. Amsterdam: Elsevier, 1996. [5] Isida M, Nemat-Nasser S. A uni®ed analysis of various problems relating to circular holes with edge cracks. Eng Fract Mech 1987;27:571±91. [6] Sendeckyj GP. Elastic inclusion problems in plane elastostatics. Int J Solids Structures 1970;6:1535±43. [7] Starnes Jr. JH, Rhodes MD, Williams JG. E€ect of impact damage and holes on the compressive strength of a graphite/ epoxy laminate. In: Pipes RB, editor. Nondestructive evaluation and ¯aw criticality for composite materials (ASTM STP 696). ASTM 1979. p. 145±71. [8] Xiong Y, Poon C, Straznicky PV, Vietingho€ H. A prediction method for the compressive strength of impact damaged composite laminates. Compos Struct 1995;30:357±67. [9] Soutis C, Curtis PT. Prediction of the post-impact compressive strength of CFRP laminated composites. Compos Sci Technol 1996;56:677±84. [10] Wang J. Prediction of post-impact compressive strength of composite laminates using an inhomogeneity model. J Compos Mater, in press. [11] Wang J, Karihaloo BL. Mode II and mode III stress singularities and intensities at a crack tip terminating on a transversely isotropic-orthotropic bimaterial interface. Proc Roy Soc London 1994;A444:447±60. [12] Dundurs J. Edge-bonded dissimilar orthogonal elastic wedges. J Appl Mech 1969;36:650±2. [13] Savin GN. Stress concentration around holes. London: Pergamon Press, 1961 [English Edition, translated by W Johnson]. [14] Shan HZ, Chou T-W. Transverse elastic moduli of unidirectional ®bre composites with ®bre/matrix interfacial debonding. Compos Sci Technol 1995;53:383±91. [15] Tsai SW, Hahn HT. Introduction to composite materials. Westport: Technomic Publishing Co, Inc., 1980. [16] Zhu X, Qi Z, Wu R, Leung WP, Choy CL. Elastic moduli of unidirectional carbon ®bre composite. In: Loo TT, Sun CT, editors. Proc Int Symp Composite Materials and Structures. Westport: Technomic Publishing Co, Inc., 1986. p. 135±40.