INHOMOGENEITY PROBLEM REVISITED VIA THE MODULUS PERTURBATION APPROACH H. Department
of Mechanical
FAN, L.
M.
Engineering.
KEER
and T.
Northwestern
MURA
University.
Evanston,
IL 60208, U.S.A.
r\btnlct--The present paper reconsiders the inhomo~eneity problem via an elastic modulus perturbation approach. By using the eigenstrain concept, the inhomo~eneity problem is converted to a series of inclusion problems in the present modulus perturbation procedure. A homogeneous state is chosen as reference state or the 0th order solution. Some extensions of the modulus perturbation approach by example are also discussed in the paper.
I. INTROD~~lON The investigation
of the field of an inhomogcncity cmbeddcd in a matrix has long been an
attrnctivc subject in micromcchanics and the mechanics of composites. However. nitions
of the inclusion
and the inhomogcncity
have not been unified.
investigation
the terminology of Mura ( 1982) for characterizing
gcncity will
hc adopted. The
inclusion
the def-
In the present
both inclusion and inhomo-
is dofincd as a s~;b~I~~rn~linhaving cigcnstrains
prcscribcrl. while the inll~~Ill~~~cncity is dcIincd as a subdomain with a prescribed modulus that is tIill+crcnt from the rcmaindcr of’ tflc m;ttcri;ll.
It is noted that the inhomogencity can
bc cquivalcnt to an inclusion with a proper cigcnstrain that dcpcnds on the external loading and gcomctry of the body. 1Jnfortunatcly. is too rlitlicult
irll~~~t!~~)~cflcitycmbcdclcd in an infinite Esholby solution. For cxamplr. Moschovidis
the cquivalcnt cigcnstrain of an inhomogcneity
to he obtained for conligurations
other than that of a single ellipsoidal
matrix,
as solved by Eshclby (1957).
thcrc arc only a few other configurations
the interaction and Mura (1973,
of two inhomogeneities
Besides the
that have analytical solutions.
in an infinite matrix was studied by
and a semi-spherical inhomogeneity near a free surface has
recently been provided by Tsuchida
Ed al. (1990). Generally, solutions
of this type are too
complicated to be used in further applications. As another branch of micromechanics, scribed cigenstrains
the study of the inclusion
(1976) dcvelopcd a cuboidal inclusion solution,
of an inclusion inclusion
is found in metallurgical
research, such as in studies of phase trans-
whcrc the cigcnstrain is known. Also, it is expected that the inclusion solutions
may be used for the study of inhomo~cneities I%?.
and Laird
and Sro and Mura (1979) solved the case
near the free surface of a half-space. A direct engineering application of the
solution
formation,
problem with pre-
has attracted many researchers. Among others, Snnkaran
by the equivalent inclusion method (Mura,
Section 22). In thr prrscnt paper, thr inhomogeneity
perturbation problem to
problem is formulated via the elastic modulus
approach. It is noted that this procedure (Section 2) turns the inhomogeneity a11
confi~tlr;ition.
inclusion problem. The 0th order solution say a homogcnsous
is obtained in a properly selected
nlcdium. which is chosen to be as simple as possible.
The higher order solutions tt z I show the gcncral properties of an inclusion problem. The corrcsponrling cquivalcnt cigcnstrains for tt z I arc obtained from the solutions of previous orders. Through the present perturbation approach, the existing inclusion solutions (Mura, I%?) will have found a new class of applications. which brings additional importance to thcsc solutions. Meanwhile. as it will be pointed out. the perturbation procedure itself. i~ltl~~ol~~hgreatly simplifying
the inhomoEeneity problem. restricts the solution somewhat. It is easy to see that the diffcrencc between C,,&, (matrix modulus) and C$, (inhomogeneity modulus) is required to bc small compared to C,,,, itself. Other remarks on the advantages and disadvantages of the proccdurc will bc made in the last section of the paper. 25R3
2%
H. FAY.
ct ul.
A similar technique had been used by Walpole (1967) to tnzat an inclusion in an anisotropic medium. He used an isotropic solution as the reference state. the 0th order solution. which was obtained by Eshelby (1957). Other related works. such as Gao (1991). will be reviewed in the following sections.
1. PERTURBATION
PROCEDURE
Consider an elastic medium D with elastic modulus Co, in which there is a subdomain Q contained in D with a different modulus C*. We call R the domain of the inhomogeneity and D-R as the domain of the matrix. One may write Hooke’s law in the matrix as, *,, =
and in the inhomogeneity
D-R
(14
Q.
(lb)
as. 01,
The perturbation
in
C,,kt&kI
=
in
C;k,&k,
is carried out as follows :
(2) where the parameterj’may
be chosen by the ratio : ./‘ = (C’ -
C”,/C”.
whcrc C,, iIt'd c'* are the representative components in c$, and c;,,, rcspcctively. All the stress and strain tields are expitndcd as follows : 61) =
cJp/+/u;, + O(/2).
t:,,= E:;+/‘t,‘) +O(J’?,.
(31. 3b)
By identifying terms with the same power of,/: the first two order approximations R
in
are:
O(J),
CT:;= c:;,,t:,:.
(4a)
a,‘)= C:;trE:,+ c,‘&:.
(4b)
1st order solution :
In general, one has : O(f"),nth order solution, n > I : 6;; =
c,$,i:;, + c,;,,t$; ’ .
(4c)
On the other hand, in D-R, (5) By using the eigenstrain concept (Mum. 1982). eqn (4b) may be rewritten as:
’ = cp,kh:t--E:;). or/ where ~1; is the so-called eigenstrain, which is obtained by equating (4b) to (6).
(6)
Inhomogeneity problem 8
1.
=
2585
-(CO)-‘(C’)cO.
(7)
In terms of the eigenstrain and Green’s function of the corresponding order approximation
is expressed as,
u!(x) =
where ( ).i = a(
domain, the first
C;h~L(~')G,,,(~,
x’) dr’
(8a)
)/a.~;.
Green’s function G,,(x. x’) in the above equations is the displacement in the xi direction at point x when the unit point force is acting in the direction .vj at x’. The integrals in eqn (8b) do not exist in the sense of Riemann Appendix)
when a numerical calculation
integrals.
It is seen that the higher order approximations Although
A special treatment
is needed (see
is carried out. can be performed
in the same fashion.
in principle the higher order (n 3 I) can be obtained. most practical applications
only need the first order modification. In the following first order formulation most of the time. Combining cqn (7) with (Ha). the perturbation hc rcwrittcn as :
If,’(x) = -
I
c,:&,(x’)G
section. WC will mention only the
of the displaccmcnt
equation
,,,,(x, x’) dx’
(8a) may
(94
11
,
G,I.,(x.x’) + &G&x. Expression (9a) may be obtained by omitting
x’)
) dx’.
the cigcnstrain concept (Gao.
order to connect the inhomogrncity
and inclusion
problems,
mention and calculate the eigrnstrain
as necessary in the following
1991). In
the present analysis shall sections.
3. EXAMPLES
In order to have a quantitative the Eshclby (1957)
sense of the accuracy of the perturbation
inclusion solution is re-examined.
An infinitely
approximation,
long solid cylindrical
inhomogeneity is embedded in an isotropic infnite matrix. The far field is loaded by pure shear stress, u::. By using the equivalent inclusion concept (Mura, 1982, Chapter I I), the disturbance due to the inhomogeneity
can be expressed as :
01, = C$,L,, =
C,“,,,(Et,-GA
(10)
whcrc o,, and c,, are the disturbed fields. and the E: is the so-called eigenstrain. From the Eshclby solution, the strcsscs in the inhomogeneity arc known to be uniform. The only non-zero componrnt under shear is o,!. The exact solution of the eigenstrain is :
whcrc Ai1 = IL* - 11. The constant S, :, ? for the ellipsoidal (1957). On the other hand. an approximated
inclusion was given by Eshelby
solution can be obtained
from the perturbation
proccdurc given in Section 2 above. The 0th order solution is taken as the uniform shear ” (7,: = ZCCE:~.In the homogeneity. the first order stress and strain are related as :
H. FA-.
2586
ZI
ul.
(17)
L&1 where
‘I-v
co =
$1
-_
I-2w
l-71
I’
I-21,
,
0
I-1’ -I-71.
0
0
I,
0
I--V*
w
__
*
c* = 7p+
1’*
I -7y*
1-p
\'* I_zy*
I -VI’+ 1 -'y*
0
0
0 (13)
0 I J
and
(CO)-1
=;
i For simplicity,
--Y
0
V
I-w
0
0
0
Ii
v = \** is taken in the following
order is obtained by substituting
For the case
I-V
tbrmulatioii.
the above Ihrmtk~o
The cigcnstrain
of’ pure shear. eqn (14) is the first term of the Taylor
By substituting
the eigenstrain
into eqn (t(b). the first
obtained. Having noted that the intrgrnl
(or the lirst
into crln (7) :
expansion ofeqn (I I).
order correction
(Hb) with a uniform
cigsnstrain
on the strain
is
has been carried
out by Eshelby (1957), E,) can be expressed in terms of the so-called Eshelby tensor S as:
In most cases, the 1st order solution present approximation
With
the Eshclby
approximation
provides II good approximation.
can be obtained by comparing (IS)
tensor the discussion
The error
can be extcndcd as follows.
The next order
is read as :
By following
ot‘ the
with the exact Eshclby solution.
the same procedure, the nth order eigcnstrain is obtained as :
Inhomogeneity
2587
problem
(16) where n > I.
It is seen that the sum of all the orders of eigenstrains converges to the exact solution, eqn (I I ). in the range of convergence. i.e.
ET: =
=&
.
(17)
I *leaact)
A numerical range of Ajl can be reached if the shape of the inhomogeneity ratio are given. For a circular cylinder with v = 0.25. the above condition
and Poisson’s
is read as :
which will guarantee that the series (I 7) will convcrgc to the exact solution (I I). Information approximation
can bc dorived from this cxamplc concerning the error of the perturbation at the first and higher orders. Usually.
obtained. while the higher order approximations the perturbation
only the first order can bc easily
bccomc more complicated.
solution does not show the extra complexity
It is clear that
when the inhomogcneity
is
not of elliptical shape. The cfl&t of the shape on the stress and strain fields may be discussed by using this approximation pcrturbntion
schcmc. It is noted that this benefit,
gained through
procedure, is at the cxpcnsc of the restriction that IACJ = IC* -Cl
As another cxamplc without a closed form solution, consider an inhomogeneity a free surface, as shown in Fig.
I. where the inhomogeneity
is taken again as the 0th order approximation.
near
has a radius of I and the
distance between free surface and the center of the inhomogeneity half-plane
the
is small.
is h. The homogeneous
The first order eigenstrains are
easily obtained as the following :
E:‘, =
”
E::
=
E,,,
E:)* = 0,
for biaxial tension loading. The eigenstrains are obtained through
Fig. I. An inhomogeneity
near ;Lsurface.
(18)
Fig. Z(a). Stress distributions along x2 = 0 with ir T 1.S fdashcufline is uj: and sofid one is al,).
(14)
provided that the Poisson’s ratios are the same in the homogeneity and matrix. It is obvious that Green’s function used in the integration [eqn (8)] is for a half-plane only. For the sake of completeness, the Green’s function is listed in the Appendix of the present paper. Stress distributions for the biaxial case are shown in Fig. 2. It is shown that the Esheiby solution for an inhomogeneity embedded in an infinite space can be used as long as the inhomogeneity is beneath the free surface at roughty twice its radius. Another fact revealed by this calculation is that there is a high tangential stress at the free surface when the inhomogeneity is close to the free surface, This resutt indicates that fracture failure may first occur at the free surface for a soft inhomogeneity near the surface. On the other hand, the fracture may be observed at the interface of the soft inhomogeneity and the matrix when the i~homogeneity is deeply embedded in the matrix, 4. BOUNDARY
CONDITIONS
AND
MfXED
BOUNDARY
VALUE
PROBLEMS
Beside the perturbation of the governing equations in the inhomogeneity, boundary conditions are also required. It is easy to develop solutions for the stress prescribed condition as:
Inhomogeneity
t = to
2589
problem
on
s,
(19a)
S,.
(19b)
and for the displacement prescribed conditions as : u =
u” on
The higher order prescribed data are set to zero. Additional care must be taken when prescribing the boundary conditions for the mixed boundary value problem [see Erdogan (1978)]. It will be noted that the higher order boundary conditions also depend on the previous order solutions. The tractions and displacements in this case must be expanded as infinite series : t = t”+ft’+...
and
u = uO+fu’+...
where to and u” are known, while t’ and u’ are to be determined from a 0th order solution. In order to clearly present this procedure, the case of a rigid indenter pressed on a semiinfinite plane with a near surface inhomogeneity is considered (as shown in Fig. 3). By assuming the Hertzian contact area c >>d (inhomogeneity size). the 0th order boundary conditions are given as : to = p”e:,
w
u” = ttze,.
P’ Order and Olh Order
Fig. 3. An inhomogeneity
near a contacted
region.
259u
fl.
They arc uniform
in the .r-direction.
F*\
‘I/
C’f
The stress and strain ticldz nar
the surf~e
arc [WC
Johnson ( 19X5)] :
f7,
=
6,
p”.
=
IT,,
r,,
=
=
CT,,
=
0
(21)
and
(‘2)
The corresponding cipcnstrain is obtaincd from the equation
(23)
From these cigcnstrains.
the corresponding
displaccmcnt dis[urbancc at z = 0 (harrctl
variables) is given. cqn (&I). as :
li I (x) =
11;(x)
=
I c;:,,,,,,:/,;,
(s’)G,,,t(s. x’)tls’
LI
.
c ,t,,,,,~:l!,,c/,.t (s. s’ ) dx’. *I
=-
(24)
i0
This
is ;I displacement that is not allowed at the contact surface tluc to the prcscncc of
the rigid indenter. Thus, The
correction
equation
a negative distribution
to the surface traction
of the above displacement at : = 0 is
at this order is gikcn by the L)llowing
added.
intcgr;tl
:
‘. I” Cc)‘,,s =
I
,
.\‘-
.Y
nE ‘(I
-\a?)
f’li 1 I?.\.
(25)
where c >> L>> J is properly selected. The stress traction at z = 0 disturbed by the inhomogeneity is expected as in Fig. 3.
The stress or displacement fields of the lirst order consists of two parts. The lirst part is computed directly from eqn (8). The second part is obtained by considering p’(.r) acting at z = 0 with ;I homogeneous half plant. Thcexact formulation ofthis prohlcm iscomplicatcd hcc~~usc ofthc interaction between the inhomogcncity and stress distribution within the cont;lct area. Miller and Kccr (19x3) and Bryant CI trl. (IWJ) forniulatctl the prohluni of twv-tlinicnsion;ll contact ol';~ri indenter interacting with ;I near surface inhomogeneity. They provitlcd numerical results for the inlcrxlion bclwccn 3 contact intknlcr al the surfax and ;I nwr surface circular void or rigid inclusion beneath the indcntcr. Marc gcncrally shapccl inhomo~cncitics arc too tlillicult to be solved by an exact formulation. On the contr;iry. lhc pcrturhation approxh can trcal arbitrarily shapctl inhomogcncitics. cvcn for the three-dimcllsion~ll USC. Thu present pcrturhalion
schcmc has ;Ilso dccouplctl
the inhomogencity
and surface traction
lhc tlctcrniination
distribution.
of tlic cquiv;~lcnt cigcnslrain
ot
2591
Inhomogeneity problem 5. OTHER
The modulus
APPLICATIONS
perturbation
VIA THE
approach
MODULUS
PERTURBATION
is not limited
to the inhomogeneity
problem
discussed above. Next. two problems are given, which are extensions of the perturbation approach.
Elastic fiulds in an elastic material of variable modulus When the elastic modulus is a function of spatial co-ordinates. dramatically
more complicated
compared
to the homogeneous
and Erdogan
(1983) considered a crack in an elastic medium
exponentially
as a function of a co-ordinate.
the problem becomes
case. For example.
Delale
with its modulus changing
It is noted that the perturbation
solution can
be obtained for this kind of problem under some restrictions. The 0th order solution is still taken as the solution in the homogeneous the first order equation is written as :
where C’(x)
is a known
function.
Section 2 and no special treatment A systematic formulation moduli perturbation
The remainder
of the formulation
is the same as in
is required for this problem.
of fracture analysis of nonhomogeneous
approach
body. while
was recently provided by Gao (1991).
stress intcnsitiy factor of a crack in the nonhomogeneous
materials via the He emphasized
elastic medium. Combined
his puhlishcd works. hc has cxtcnsively studied the perturbation
the with
stress intensity factors of
a crack in the elastic medium with C(x).
AS iI furthrr cxtcnsion of the modulus perturbation [Fig. 4(a)]
under a certain
rcfcrcncc material
is considered
approach,
a nonlinear c&tic
hcrc. In Fig. 4(a),
curve (I)
which will bs used to construct the 0th order solution.
the real stress-strain contiguration
loading
relation.
The 0th order solution
under consideration,
Curvs (2) gives
is the linear elastic solution
while the first order correction
Furthermore
(27)
theory is applied, then one may have C ’ = C ‘(J&
Jz = Is,,s,,. s,] = an elosto-plastic
a,, -
for a
is
4, = c/rrt::,+ c,:t,h,)4L If JZ deformation
body
shows a
where
fotk6,,.
(28)
analysis with a stress-strain
be treated by the same procedure, although an extra calculation
curve as in Fig. 4(b) can in the 1st order formulation
is required to determine the plastic zone and the region where the eigenstrain is prescribed. It is noted that there is (are) some subdomain yield limit rr,,. The determination Mises, as :
R where the effective stress exceeds the
of R is achieved according
to a yield criterion,
e.g. von
whcrc Q,~,is the yield stress in uniaxial tension. With the above cigsnstrain and $2 known, one obtains the first order stress and strain through the integration equation (8). Finally the plastic zone Qt, of the 1st order where the total effective stress exceeds d,,., is dctcrmincd
by
1(hy+rs,: ) (sp, +f&‘,) 2 u,!,. It is obvious that higher order perturbations zone R,.
will further
(30) correct 52 and the plastic
H. Fx\\u t’r ul.
0
T /’
th order solution
6. CONCLUDING
REMARKS
We have presented ;I modulus perturbation and stress fields disturbed typos of inhomogeneity tages by the simplicity
by inhomogeneities.
scheme for determining Although
problems, the present perturbation and analytical form of the results.
approach. some insight is obtained for many dificult be very troublesome from both the computational required in FEM
FEM
and BEM
numerical
surface traction distribution,
the displacement can handle these
procedure still shows its advan-
By using the modulus perturbation
inhomogeneity
problems which may
complexity and the amount of CPU time
procedures.
problem discussed in Section 4. the interaction
or BEM
For
instance. the mixed boundary
between ;I near surface inhomogeneity and
can be quite costly in CPU time [see a similar caIcuI;ltion on the
elastic-plastic contact analysis by Komvopoulos (lC,SU)]. Under the present perturbation procedure. an inhomogcncity problem is convortcd to ;I series of inclusion problems. The inhomogcncitics other than those of ellipsoidal shape can be considered through the perturbation schcmc. Extensions of the ;lppiication of the modulus perturbation have been made for two examples. The elastic-.plastic analysis through the modulus perturbation scheme is easy to bc ;lpplicd to cnginccring applications. The perturbation procedure may also bc helpful for numerical c;llculations even for the GISC when the modulus dilfcrcncc of the inhomogcncity/matrix system is out of the perturbation range. For instance. in an invcrsc problem [see Gao and Muru (lY89)] the initial guess of the solution is vital for the convergence of the numerical iteration scheme. The perturbation solution may provide a rcasonablc initial guess for the itoration. It should also ht emphasircd that there arc some limitations on the application of the perturbation schcmc. First, the present approach is based on the assumption that AC,,,, is
Inhomogeneity
1593
problem
small compured to C,,,,. A quantitative restriction is shown in the first example in Section 1 where the exact solution acts as a benchmark. The appiicable range of the ~rturbation -. solution also depends on the problems and parameters under consideration. Gao (1991) showed that the perturbation solution of stress intensity factor of a Mode III crack surrounded by inhomogeneity is in good agreement with the exact solution even though Ap/p 2 2. However. for those problems where AC,,&,is small. this assumption may not prove too iimitin~ for engineering appli~tions. Since the proposed modulus perturbation scheme is a regular perturbation procedure. the approximated solution may lose some important information present in the original problem. For example, as interface crack tip singularity cannot be derived from a 0th order solution which is chosen in a homogeneous space. The singularity at the interface crack tip is a special characteristic of the original problem which must appear in the 0th order con~~L[~~tion. This example illustrates that the application of the perturbation approach should be carried out with great care in terms of the original problem and application. Acknon,lc,c/~enlt,~f~-This work was supported by the U.S. Army Research Office under contract 0019. and NSF under contract MSM-8817969.
DAALO?-BP-K-
REFERENCES Brchbia, C. A.. Tellcs. J. C. F. and Wrohcl. L. C. (IY83). B~wtrtltrr~ Ekmc*nrTechniqurs.Springer. New York. Bryant. M. 3.. Miller. G. R. and Kecr. L. M. (19X-l).Lint contxt between ;\ rinid -Fm~ elns11c .z indenter and a damavcd &dy. Q, J. ,bldr. ,.lppl. Murh. 37. 467 ~478.’ Dclalc. F. and Erdogan, F. (I’IR.3). Crock in ;I varied claslic motlulus body. J. Appl. .Clc*ch. 50, MN -614. Erdopan. F. ( 10%). Mitcd h~~inll~~ry-v~~l~~eproblems in mech;tnics. In .~~l.rh~f~i~..~ ?‘otfri~ (Edited by S. NematNasser). Vol. 4, pp. I -84. Pergamon Press, New York. Eshclhy, J. I). f 1057). The dctcrmin;&m of the elastic lield of’xn ellipsoidal inclusion ;mrl related prohlan. Frrrc. K. .SO(~.AZJI. 376 .I’%.
Tachida. E.. Kouris. D. and Jasiuk. 1. (1990). The hemispherical inhomogeneity subjected to 3 conaznrrated force. In ,~fic.ronrrc/tcrrri a& f!jhunruymrirp (Edited by C. J. Weng, M. Taya and H. Abe). Springer. New York. Walpole. L. J. (1967). The ei:tslic tield of an inclusion in an anisotropic medium. from. K. Sot. A3UU. 27&288.
APPENDIX The C&en’s
:
IIALF-PLANE
GREEN
FUNCTION
function for J half plx~e can be wriusn as:
G,,(r,x’)
= G’:,(r,x’)cG:,fx.x’),
whcrc G:,(x. x’) corresponds to the whole-plant Green function and G:,(x, x’) is called the complementary of the h;df-plume Green function. They are given by
and
G;, = Kd -[8(1-~)‘-(3-4vIjIn
Rc
Rz
(Ai) part
294
H.
FAN cr 01
(A3)
H here r, =
.A-,
-.r:.
R, = _x, cx’,.
r=(r,r,)‘:.
R, = .r: -.v’:.
R=(R,R,)“. 0 = arctan
.i = x.‘,.
c=_v,. R.
( l>
. K,, =
fi
I
8np(I -I)
The first order displacement disturbance caused by the inhomogeneity fun&on and the eigenstrain formulatton as follows:
u;(x) = Furthermore.
the stress disturbance
“,, = -
’
--A-y:
(
,: =
1::: L
= -K,r:
Ri
by taking the derivatives
Z,,,(x. ~‘)C:,_E:,,(X’)~~‘-C:I,.E::.
16clR;
Z[_~‘-r’-2c.~+2.~R,(I-2~)J f---RI
Hi
-.-.-.
with respect to x. i.c.
s
+
( i + 3c)( I - Zr) + _‘[R,(~i+2~~)-2~.ri+2Sri(I
r’:,=-K,
r;
(-
(I -2V)
is obtained
(A4)
x’) +G:,(x. x’))dx’.
L)
whcrc
r;:, =
C;,,,&(G,:(x.
is obtalned by using the above Green
R” -‘v)]
Ra 2(r’-.~:+6c.~-21R,(I-Zv)l - --_~____-~-_._..___
164 + ___
RI
3( I - 2v) l[ri-2‘.‘-4~.\‘-2.iR,(I ..-.. ~-. + _-_- __._ - ..-___
R!
R’
R’
-2v)j
>’ I6rlR,ri
+-yr-’
>
>’
I6c.\‘Ri
+7.
>
(A7)
whcrc K. = I In( I -v). f%r ;I gcncral distribution of cigcnstrain, ;1 spczial trcatmcnt on strcsr evaluation from cqn (A5) is nczdcd (Brcbbia er al.. 1983. Chapter 6). In the second example of Section 3, the uniformly distributed first order cigcnstrain can simplify the present numerical calculation procedure as:
cr;, = i”
(Z:,(x.x’)+Z;U(x.
x’))C:,_&,,(x’)dx’.
(As)
Note that the complementary part of the Green function in the above integral is not singular at x’ I x. Thus, the second part of eqn (AS) can be evaluated by using a conventional numerical integration technique. The first part of eqn (AR) is just the Eshclby solution for an inclusion embedded in an infinite matrix which was obtained in a closed form expression.