Numerical equivalent inclusion method: a new computational method for analyzing stress fields in and around inclusions of various shapes

Materials Science and Engineering A285 (2000) 229 – 238 www.elsevier.com/locate/msea

Numerical equivalent inclusion method: a new computational method for analyzing stress fields in and around inclusions of various shapes Yuji Nakasone a,*, Hirotada Nishiyama b, Tetsuharu Nojiri b a

Department of Mechanical Engineering, Faculty of Engineering, Science Uni6ersity of Tokyo, 1 -3 Kagurazaka, Shinjuku-ku, Tokyo 162 -8601, Japan b Graduate School of Engineering, Science Uni6ersity of Tokyo, 1 -3 Kagurazaka, Shinjuku-ku, Tokyo 162 -8601, Japan

Abstract The present study has attempted to develop a new computational method for the elastic stress analysis of inclusions based on the equivalent inclusion method. The proposed method can avoid the complexity of mathematics required for the analysis of non-uniform eigenstrain distributions within inclusions having various shapes. The paper is focused on the formulation for two-dimensional case. The fundamental integral equation is shown first to have a kernel with the 1/r-singularity. The two-dimensional equations are then discretized by using the triangle polar coordinates. It is shown that the adoption of this coordinate system can eliminate the singularity. Eigenstrain distributions within inclusions having various shapes were calculated by the present method in order to obtain stress distributions within them as well as those in the vicinity of the matrix–inclusion interfaces. The shapes of the inclusions described here are ellipse, circle, triangle and rectangle. The results obtained by the present method were compared and showed good agreements with those obtained by the theories and/or by the FEM analyses except for the sharp corner points of the triangular inclusions where the outside normal vectors can not be determined uniquely and thus the stress becomes singular. Published by Elsevier Science S.A. Keywords: Micromechanics; Equivalent inclusion method; Green’s function; Computational mechanics; Finite elements; Discretization; Triangular polar coordinates

1. Introduction Advanced materials such as superalloys and metal– matrix composites tend to have complex microstructures constrained by design and process requirements from recent advanced technologies. These materials are no longer homogeneous even from the viewpoint of the conventional continuum mechanics; therefore, an effective micromechanics approach is required to analyze their mechanical behavior. The equivalent inclusion method was first proposed by J.D. Eshelby, and is considered to be one of the most powerful tools for analyzing the mechanics of microscopic inhomogeneities often found in advanced materials; e.g. strengthening dispersoids, inclusions, precipi* Corresponding author. Tel.: +81-3-32604272, ext. 3356; fax: +81-3-32604291. E-mail address: [email protected] (Y. Nakasone)

tates, micro-cracks, etc. [1–3]. The method, however, requires advanced mathematics, and thus is difficult to treat or even is theoretically unsolvable in a closed form for the cases of non-elliptic and non-ellipsoidal inclusions. Little work has yet been done for non-elliptic and non-ellipsoidal inclusions. Recently, Mura et al. [4] claimed that the elastic strain and stress distributions within a pentagonal star-shaped inclusion are uniform. Being inspired by their work, Rodin [5] has derived a closed-form solution for polygonal and polyhedral inclusions with uniform eigenstrains and found that vertex singularities do exist in these inclusions except for incompressible polygonal inclusions whose angles are either p/2 or 3p/2. Nozaki and Taya [6] have investigated a regular n polygon-shaped inclusion with uniform eigenstrains. They have obtained the logarithmic singularity of the total strain fields in polygonal inclusions. The strain fields approach to the uniform distribution as the number of vertices n increases.

0921-5093/00/$ - see front matter Published by Elsevier Science S.A. All rights reserved. PII: S 0 9 2 1 - 5 0 9 3 ( 0 0 ) 0 0 6 3 7 - 7

230

Y. Nakasone et al. / Materials Science and Engineering A285 (2000) 229–238

This paper presents a new computational method for calculating the eigenstrain and elastic stress distributions within inclusions of arbitrary shapes. Inclusions to analyze need not be isotropic. The proposed method can be obtained by directly discretizing the fundamental integral equations of the type-I equivalent inclusion method. This paper is focussed on the formulation for two-dimensional case, but the three-dimensional formulation can be easily derived in a similar way. The fundamental integral equation is shown first to have a kernel with the 1/r-singularity. The two-dimensional equations are then discretized by using the triangle polar coordinates [7,8]. It is shown that the adoption of this coordinate system can eliminate the singularity. Eigenstrain distributions within elliptic, triangular and rectangular inclusions are calculated by the present method in order to obtain stress distributions inside them and those in the vicinity of the matrix–inclusion interfaces. The results obtained by the present method are compared with those obtained by the theories or by the finite element method (FEM).

2. Fundamental equations; type-I equivalent inclusion method The equivalent inclusion method can be classified into two types according to whether or not the applied external forces exist [2,3]. The type-I equivalent inclusion method can be defined as the case where the external forces are applied to an infinite body having inclusions in it, while the type-II for the case where no external forces are applied, thus the problem of internal stress field only. This chapter will describe the computerization of the type-I equivalent inclusion problem. Figure 1 depicts a schematic explanation of the type-I equivalent inclusion problem. By using the concept of the eigenstrain, the method simulates the stress or strain disturbance caused by an inhomogeneity with the elastic moduli C%ijkl occupying the domain V in an infinite matrix

Fig. 1. Schematic view of the type-I equivalent inclusion problem.

D–V with the elastic moduli Cijkl. The inhomogeneity is substituted for the domain V with the same elastic moduli as the matrix but with eigenstrains o *kl which brings about the equivalent stress or strain disturbance [2,3]. The eigenstrain o *kl is chosen so as to satisfy the following equation: A C%*ijkl (u A k,l + uk,l )= Cijkl (u k,l + uk,l − ok,l ) in V

(1)

where u A k,l is the displacement gradient induced by the applied stress s A ij at infinity, uk,l the displacement gradient disturbance caused by the presence of the inclusion V, and ( · )k,l (( · )k /(xl. Equation 1 is the necessary and sufficient condition for the equivalency of the stress and strain fields in the type-I equivalent inclusion problems. The strain gradient disturbance uk,l can be expressed in terms of the static Green’s function, Gik (x−x%). The function Gik (x−x%) can be regarded as the xi component of the displacement at point x, when a unit body force in the xk -direction is applied at point x% in the infinitely-extended material [2,3]. By using its partial derivative with respect to xl, Gijl (x−x%), uij can be written by the following equation [2,3]. uij = −

( (xj

&

V

Cklmnomn (x%)Gik,l (x−x%) dx%.

(2)

The Green’s function Gik,l (x−x%) in Eq. (2) can be obtained by the following formula: Gik,l (x−x%) =



n

−1 (3−4n)dikx¯l − (dklx¯l + dilx¯k ) 2x¯ix¯jx¯k + 8p(1− n)m x¯ 2 x¯ 4

(3) where m is shear modulus, n Poisson’s ratio, dik the Kronecker delta, and x¯l xl − x%l, x¯ x¯lx¯l. The summation convention for the repeated indices is employed in the above equation. In this chapter, the fundamental equations Eq. (1) will be discretized with the singularity taken into consideration. The description will be focused upon the formulation for two-dimensional problems, since the three-dimensional formulation is analogous to the twodimensional one, and thus is easily deductible from the procedure mentioned below. The three-dimensional formulation will be presented elsewhere. Firstly, let the inclusion V be divided into Ne triangular elements Ve (e= 1, 2, … , Ne ) by using Ne nodal points. The eigenstrain o *(x) within the eth element Ve ij can then be approximated by the values of the eigenstrain at each node of the element o *ij (I) and the interpolation function L (I) e (I=1, 2, 3; see Fig. 2a) for the standard linear triangular element [9]; i.e. (I) (I) o*(x)=L for xVe ij e o* ij

(4)

Y. Nakasone et al. / Materials Science and Engineering A285 (2000) 229–238

Fig. 2. Discretization of an inclusion embedded in a two-dimensional infinite elastic body. (a) Non-singular points (x
Q Ve ); (b) singular points (x
 Ve ).

where, for the index I, summation is taken from 1 to 3 for the two-dimensional triangular finite element. The index I in Eq. (4) denotes the local nodal number of the element Ve. The interpolation function L (I) e for node I is expressed as (I) (I) (I) L (I) e =a e +b e x1 +g e x2.

(5)

The coefficients a , b , and g in Eq. (5) are given in terms of the coordinates of the other two nodes J and K; i.e. (I) e

(I) e

(K) (K) (J) x (J) 1 x2 −x1 x2 , 2Ae

a (I) e = b (I) e = g (I) e =

(I) e

(K) x (J) 2 −x2 , 2Ae

x x x

x à à x Ã. x à à (I) 2 (J) 2 (K) 2

(6b)

For non-singular points (x Q Ve ), the integrand in Eq. (2) does not have a singularity; thus, the order of the integration and the differentiation can be changed to integrate Eq. (2) numerically by using the standard triangular elements. After changing the order of the integration and the differentiation in Eq. (2) and substituting Eq. (4) into Eq. (2), the following discretization formula can be obtained by the Gaussian quadrature procedure [10]. ui, j = −

( (xj Ni

=− %

&

Ve

+ 2Dijmno*mn (1) + (1− 2n) + (1− 2n)(2Gjmo*im (3) − Gijo*mm (3))

where Ae is the area of the triangular elements Ve and given by the following formula (I) 1 (J) 1 (K) 1

(ui −1 = [(1− 2n)(2Ajmo*im (1) − Aijo*mm (1)) (I) 4p(1− n) (x j × (2Bjmo*im (2) − Bijo*mm (2))+ 2Eijmno*mn (2)

(J) x (K) 1 −x 1 , 2Ae

Ã1 Ã Ae =det Ã1 Ã1 Ã

and J the Jacobean determinant for the transformation of the global coordinate system into the local coordinate one. Ni is chosen as 16 in the present calculations. For singular points (xVe ), i.e. in the case where the reference point x falls on point x% (see Fig. 2b), the integrand of Eq. (2) has the 1/r-singularity. In this particular case, the order of the integration and the differentiation is not interchangeable; therefore, the above discretization can not be applied. The discretization, however, can be made by the use of a finite element having the triangle polar coordinates [7,8]. This finite element is employed in the direct boundary element analyses in which body forces, and hence interior points, need be investigated [7,8]. The singularity can be eliminated in Cauchy’s sense by the adoption of the finite element having the triangle polar coordinate system for the discretization of the singular element, Vs, i.e. Ve x. Analytical differentiation of ui (x) with respect to xj, uij (x) can be obtained explicitly by known values of the nodal point coordinates x (iI ) and unknown values of eigenstrains o *ij (I) on the node x (I) having the local node number I (i, j =1, 2; I= 1, 2, 3). For two-dimensional cases, the resultant expression of uij (x) can be written as follows: ui, j =

(6a)

231

Cklmno*mn (x%)Gik,l (x − x%) dx

3

(K) % wlCklmnL (K) Gik,lj (x, hI ) J e (hI )o* mn

(7)

I=1 K=1

where hI is the local coordinates of the Ith integration point, wI the weights, Ni the number of Gauss points

+2Zijmno*mm (3)]

(8)

where x (I) (i.e. I= 1; j =1, 2) represents the cartesian j coordinates of the singular point [7,8]. The coefficients of the nodal values of eigenstrains Aij, Bij, Gij, Dijkl, Eijkl and Zijkl are expressed only by the coordinates of the three nodal points of the triangular element. The detailed expressions of the coefficients are listed in Appendix to the present paper. Substitution of Eqs. (7) and (8) into Eq. (1) gives the following simultaneous equations for unknown variables o *kl(N) in the element Ve : (N) (Cijkl − C%ijkl ) % B (N) − Cijklo*kl (M) klmno* mn N  Vre

1 A = (C%ijkl − Cijkl )C − klmns mn,

(9)

where the superscripts N and M represent the global number of each node of the element and that of the x (M) of interest, respectively. N and M range from 1 to Nn where Nn is the total number of the nodal points. Equation 8 form some parts of the diagonal components of the coefficient matrix for the unknown values of o *kl(M) in the above Eq. (9). The eigenstrain distribu-

232

Y. Nakasone et al. / Materials Science and Engineering A285 (2000) 229–238

(l+ 2m)dkm − (l+ m)nknm m(l+ 2m)

tion o *(x) within the inclusion V can be obtained by kl calculating Eq. (9) for all the Nn nodal points in the Ne triangular elements and then solving a set of 3Nn linear equations with respect to the unknown variables o *kl(M) on the nodal points x (M) within V (M=1, 2, … , Nn ). j Namely, the resultant simultaneous equations can be written in the matrix notation as follows:

Nkm (n)D − 1(n)=

[K]{o*}= [C − 1]{s A},

The computer code for the present numerical equivalent inclusion method has been developed based on the formulation described in the previous chapter. Figure 3 shows the flow chart of the present method. The code calculates first the nodal values of the eigenstrains o *ij (M) in the inclusion V (M= 1, 2, … , Nn ) by solving the algebraic equation Eq. (11). Then, the code calculates the internal stress distribution s in ij (x) in the inclusion and the interfacial stress distribution s out ij (x) in the vicinity of the periphery of the inclusion. Since Eq. (2) is the solution for the inclusion in the infinitely-extended two-dimensional elastic body, the present method requires only the inclusion region to be meshed. (M) has a useful The coefficient matrix M  V B (M) klmn for o * kl property similar to that of the Eshelby tensor Sklmn, i.e. the value of each component of the coefficient matrix is the same when the matrix material, the shape of the inclusion and the discretization of the inclusion region are the same. The present code adopts the parametric least square method, a modified version of the parametric spline interpolation method, in order to obtain the outside normal ni on the circumference of an inclusion having an arbitrary shape. In the case of a triangular inclusion illustrated in Fig. 8, which will be shown later, for example, the outside normal ni at the corner point A can be approximately obtained as a unit vector having components of (− 1, 0). In a strict sense, the outside normal ni can not be defined at the point A; nevertheless, the present code approximately calculates the outside unit normal vector by rounding off the sharp corner via the parametric least square method. Interpolating the periphery of an inclusion, the code automatically calculates the interfacial stress s out ij (x) by the use of Eqs. (12) and (13). The present computer code can calculate the stress field in and around inclusions having different shapes. The inclusions are not necessarily isotropic, so that the code can also analyze anisotropic inclusions. The elliptic inclusions and holes having different aspect ratios were analyzed first to evaluate the present code. Then polygonal inclusions such as triangular and rectangular inclusions were analyzed by the present computer code. Table 1 shows the mechanical properties of the matrix material and the inclusion materials investigated in the present study. Steel was chosen as the matrix material. SiC and Ti–6Al–4V were chosen as the inclusion materials in order to examine the effects of the matrix material and inclusion material

(10)

where [K] is the coefficient matrix for the unknown nodal values of eigenstrains {o*}, [C − 1] the elastic compliance and {s A} the applied stress at infinity. Consequently, the resultant eigenstrain distributions within the inclusion V can be obtained by the following formulae. {o*}=[K] − 1[C − 1]{s A},

(11)

Substitution of Eq. (2) into Eq. (1) brings about the formula for the interior stress distribution s in ij in the inclusion of interest. The following equation Eq. (12) is used to calculate the stress distribution in the vicinity of the interface of the inclusion and the matrix: −1 in (n) + o*mn ] s out ij =s ij + Cijmn [ −Cklpqo* pqnlnnNkm (n)D (12)

where nl is the unit outside normal on the interface and Nkm (n)D − 1(n) given by Eq. (13) for the isotropic matrix material with Lame’s constants l, m and the Kronecker delta dkm [2,3].

Fig. 3. Flow chart of the numerical inclusion method developed by the present study.

(13)

3. Computer code for the numerical equivalent inclusion method

Y. Nakasone et al. / Materials Science and Engineering A285 (2000) 229–238 Table 1 Elastic properties of the matrix and inclusion material used in the present analysis Matrix material

Young’s modulus, E (GPa) Poisson’s ratio, n

210 0.30

Inclusion material SiC

Ti–6Al–4V

410

110

0.16

0.34

233

between the maximum and the minimum values obtained by the present calculations. The error bar becomes longer for smaller and larger values of the aspect ratio a/b. The cause of longer error bar may be attributed to the flatness of the triangular element. Table 2 lists the relative errors, namely, Žcalculated average value −Žtheoretical value Žtheoretical value × 100% obtained for different aspect ratios and different mesh discretizations of the SiC elliptic inclusion. The magnitudes of the relative errors for the average values of internal stress disturbance Ds in 22 and eigenstrains o * 11 and o *22 become smaller for finer mesh discretization for each value of a/b. The concentration factor is defined by the ratio of the maximum interfacial stress to the applied stress; A out i.e. a [s out 22 ]max/s 22. The maximum stress [s 22 ]max appears at points B and B% for the SiC inclusions, whereas it appears at points A and A% for the Ti– 6Al–4V inclusions. Closed symbols in Figs. 5a and b

Fig. 4. The present elliptic inclusion model (b/a = 0.5) and its finiteelement discretization.

combination on the stress distribution in and around the inclusion. The former inclusion material is harder, whereas the latter is softer than the matrix material.

4. Results and discussion

4.1. Elliptic inclusions and holes Figure 4 shows the present model of an elliptic inclusion in a two-dimensional infinite medium and its finite-element discretization. The total number of nodal points is 120 and that of linear triangular elements is 192. The value of the aspect ratio a/b was changed from 0.25 to 4. The remote uniaxial stress sA 22 = 200 MPa is applied to the infinite body having the elliptic inclusion in it. Figures 5a and b show the internal stress s in 22 versus a/b diagram and the stress concentration factor a versus a/b diagram, respectively. Although the theory claims that the value of s in 22 should be constant, the value of s in numerically obtained by the present 22 method slightly varies within the elliptic inclusions. The solid circle symbols in Fig. 5a indicate the average value of s in 22 within the inclusions, showing a good agreement with the theoretical curve. The error bars attached to the solid circles indicate the range

in Fig. 5. Internal stress s in 22 and stress concentration factor a [s 22]max/ sA vs. aspect ratio a/b diagrams for elliptic inclusions of SiC and 22 Ti – 6Al – 4V. (a) Internal stress s in 22 vs. aspect ratio a/b; (b) concentraA tion factor a [s out 22 ]max/s 22 versus aspect ratio a/b.

Y. Nakasone et al. / Materials Science and Engineering A285 (2000) 229–238

234

Table 2 Relative errors obtained for different mesh discretizations and aspect ratios of SiC elliptic inclusions Aspect ration, a/b

0.33

0.50

0.75

1.0

1.33

2.0

3.0

Total number of nodes, Nn

72 120 180 252 72 120 180 252 72 120 180 252 72 120 180 252 72 120 180 252 72 120 180 252 72 120 180 252

Total number of elements, Ne

108 192 300 432 108 192 300 432 108 192 300 432 108 192 300 432 108 192 300 432 108 192 300 432 108 192 300 432

indicate the results obtained by the present method, showing that they agree well with the theoretical curves indicated by the solid lines. Both Figs. 5a and b illustrate that the behavior of the internal stress s in 22 and that of the stress concentration factor a are different between the harder (SiC) and the softer (Ti – 6Al–4V) inclusions. Figure 6 shows the stress concentration factor a versus aspect ratio a/b diagram. The solid line depicts the theoretical line expressed by the well-known formula, a= 1+2(a/b). The calculated values of a vary linearly with the aspect ration a/b and show a good agreement with the theoretical line. Due to the flatness of the elements, however, the discrepancy between the calculated and the theoretical values becomes large for the small and the large values of a/b except for the point of a/b= 3.0. Figure 7 compares the calculated and the theoretical interfacial tangential stress distributions s out u around a circular hole. Open circular and triangular symbols indicate the values of s out calculated by the u present method and those analytically obtained by the conventional equivalent inclusion method, respectively. The solid line indicates the theoretical s out distribution u A expressed as s out /s =1 + 2 cos(2u). The calculated u 22 values of s out agree well with the two theoretical results. u

Relative error (%) Ds in 22

o *11

o *22

−3.34 −2.88 −1.76 −1.07 −4.02 −2.88 −2.09 −1.10 −0.976 −0.698 −0.518 −0.265 −7.93 −5.77 −4.13 −3.91 −11.2 −8.28 −5.90 −5.70 −19.2 −14.5 −9.99 −9.26 −1.78 −1.35 −1.01 −0.860

−1.94 −1.69 −1.11 −0.821 −1.84 −1.27 −0.851 −0.646 −1.73 −1.24 −0.863 −0.666 −1.78 −1.29 −0.929 −0.622 −1.74 −1.29 −0.922 −0.633 −1.842 −1.395 −0.988 −0.697 −2.08 −1.71 −1.15 −0.797

−3.82 −3.54 −3.10 −2.87 1.33 0.965 0.663 0.452 1.38 0.965 0.663 0.452 −1.45 −1.05 −0.756 −0.576 1.54 1.10 0.782 0.614 1.65 1.22 0.854 0.668 2.02 1.56 1.09 0.823

4.2. Triangular inclusions Figure 8 shows the triangular inclusion model investigated in the present and its finite-element discretization. The total number of nodal points is 135 and that of linear triangular elements is 220. The aspect ratio a/b of the triangle is 2.0. The remote uniaxial tensile stress

A Fig. 6. Stress concentration factor a [s out 22 ]max/s 22 as a function of aspect ratio a/b for elliptic holes.

Y. Nakasone et al. / Materials Science and Engineering A285 (2000) 229–238

Fig. 7. Comparison of the calculated and the theoretical interfacial tangential stress distributions s out around a circular hole. u

235

calculated by the parametric least square method. In the case of the hard SiC inclusion, the internal A stress s in 22 is less than the applied tensile stress s 22 in the vicinity of the point A and is increased monotonically along the x1-axis. The value of s in 22, however, is larger than the value of s A 22 in the most parts of the SiC inclusion. The tendency of the stress distribution along the x1-axis is completely reversed for the soft Ti–6Al– 4V inclusion: i.e., s in 22 is tensile in the vicinity of the point A and is decreased monotonically along the x1A axis. The value of s in 22 is smaller than s 22 in the most parts of the Ti alloy inclusion. The analysis was also made for the triangular inclusion having a smaller aspect ratio a/b= 1.0. The tendency of the results obtained was the same as mentioned above but showed a better agreement with the FEM results.

4.3. Rectangular inclusions Figure 10 shows the rectangular inclusion model investigated in the present study and its finite-element discretization. The total number of nodal points is 231 and that of linear triangular elements is 400. The aspect

Fig. 8. The present triangular inclusion model (b/a= 0.5) and its finite-element discretization.

sA 22 = 200 MPa is applied to the infinite body having the triangular inclusion in it. Figure 9a and b compare the stress distributions along the x1-axis within the inclusion and at the inclusion–matrix interface obtained by the present method and by one of the commercial FEM code called ANSYS for the SiC and the Ti – Al – 4V inclusions, respectively. Open symbols indicate the results obtained by the present method NEIM and closed symbols by the FEM. The results obtained by the two methods agree well with each other for the both inclusion materials expect for the vicinity of the sharp corner point A and for the point B on the base of the triangular inclusion. The stress discontinuity occurring at the points A and B is simulated by the present method. The order of the singularity ls at the point A is expressed in real number; i.e. ls =0.09 for the SiC inclusion and ls = 0.02 for the Ti – 5Al – 4V inclusion [11]. Thus, in reality, the stress should be infinite at the point A; nevertheless, the results obtained by the present method remain finite values, because the interfacial stress at the sharp corner can be approximately

Fig. 9. Stress distributions s22 along the x1-axis within and at the periphery of triangular inclusions of SiC and Ti – 6Al –4V. (a) SiC inclusion; (b) Ti – 6Al – 4V inclusion.

236

Y. Nakasone et al. / Materials Science and Engineering A285 (2000) 229–238

Fig. 10. The present rectangular inclusion model (b/a = 0.5) and its finite-element discretization.

Al–4V inclusions, respectively. Open symbols indicate the results obtained by the present method NEIM and closed symbols by the FEM. The results obtained by the two methods agree well with each other for the both inclusion materials. The internal stress s in 22 is larger than the applied stress s A 22 everywhere in the hard SiC inclusion. The internal stress minimum at the center of the inclusion, whereas it becomes maximum in the vicinity of the inclusion–matrix interfaces. The stress just out of the inclusion–matrix interfaces s out 22 is lower than the applied stress s A 22. The tendency of the stress distribution along the x1-axis is completely reversed for the soft Ti–6Al–4V rectangular inclusion as observed in the case of the triangular inclusion, i.e. s in 22 is smaller than s A 22 everywhere in the Ti alloy inclusion. The internal stress becomes maximum at the center of the inclusion, whereas it becomes minimum in the vicinity of the inclusion–matrix interfaces. The stress just out side of the interfaces s out 22 is higher than the applied stress s A . 22

5. Conclusions

Fig. 11. Stress distributions s22 along the x1-axis within and at the periphery of rectangular inclusion of SiC and Ti–6Al–4V. (a) SiC inclusion; (b) Ti – 6Al –4V inclusion.

ratio a/b of the rectangle is 2.0. The remote uniaxial tensile stress s A 22 =200 MPa is applied to the infinite body having the rectangular inclusion in it. Figures 11a and b compare the stress distribution along the x1-axis within the inclusion and at the inclusion–matrix interface obtained by the present method and the FEM code for the SiC and the Ti–

The present paper investigates the computerization of the type-I equivalent inclusion problems in order to avoid the complexity of the mathematics required for the analysis of non-uniform eigenstrain distributions within such inhomogeneities as inclusions, precipitates, etc. This paper is focused upon the formulation for two-dimensional cases. The fundamental integral equations were shown first to have kernels with the l/r-singularity. The integral equations were then discretized by using the finite element having the triangle polar coordinates. The discretization was shown to be very similar to the one, which is often encountered in the direct boundary element analysis where interior points need be investigated. The present paper then showed some of the examples of the analyses done by the numerical equivalent inclusion method; e.g. stress distributions in and around inclusions of various shapes such as ellipse, circle, triangle, rectangle, etc. Eigenstrain distributions within these inclusions are calculated to obtain stress distributions within them as well as those in the vicinity of the interfaces between the matrices and the inclusions. The results obtained by the present method were compared and showed good agreements with those obtained by the theories and by the FEM analyses expect for the corner points of polygonal inclusions where the outside normal vectors can not be determined uniquely and thus the stress becomes singular.

Y. Nakasone et al. / Materials Science and Engineering A285 (2000) 229–238

Appendix A. Appendix

=

The coefficients of the nodal values of eigenstrains Aij, Bij, Gij, Dijkl, Eijkl and Zijkl described in Chapter 2 are expressed by the following equations: Aij

 !



1 (I1 (I2 − dij I1 +x (21) +x (32) i i (1) 2 (x (1) (x j j

 "

(A1)





1 1 I7 (4a1a3 − a 22)2 2

=

(I1 1 (I2 Bij − dij (I1 −I2) +x (21) − (1) i (1) (x j (x j 2 + x (32) i Gij





(I2 (I3 − (1) (x (1) (x j j



(I2 (I3 1 − dij I2 +x (21) +x (32) i i (1) (x j (x (1) 2 j

Dklpq

(A3)

1 ( (X 1kpqI4 +X 2kpqI5 +X 3kpqI6 +X 4kpqI7) 2 (x (1) l

(A4)

+X 3kpqI6(I6 −I7) + X 4kpqI7(I7 −I8)}

(A5)

1 ( Zklpq (X 1kpqI5 +X 2kpqI6 +X 3kpqI7 +X 4kpqI8) 2 (x (1) l

(A6)



=tan − 1

1 1 I2 (4a1a3 −a 22)2 2

=

1

1 dr2 2 +a r a 2 2 +a3r 2 0 1



 

a2 +2a3 a −tan − 1 2 , 4A 4A

 

&

(A7)



& 

1 1 I3 (4a1a3 −a 22)2 2



!

r 22 dr2 2 0 a1 +a2r2 +a3r 2

(A8)

=



"



1 1 I6 (4a1a3 −a 22)2 2

&

r 22 dr2 2 2 0 (a1 +a2r2 +a3r 2) 1

1 kpq

2 kpq

3 kpq

(A16) 4 kpq

and X appearing in Eqs. (A4) to X ,X ,X (A6) are defined as follows: (21) (21) X 1kpq = x (21) k xp xq

(A17)

(21) (21) (32) (21) (21) (32) + x (21) + x (21) X 2kpq = x (32) k xp xq k xp xq k xp xq (A18) (32) (32) (21) (32) (32) (21) + x (32) + x (32) X 3kpq = x (21) k xp xq k xp xq k xp xq (A19)

(A20)

Finally, a1, a2 and a3 in the solutions of the definite integrals, Eqs. (A7) to (A14), are expressed in terms of x (21) and x (32) . i i (A21) (A22)

2

)2 a3 % (x (32) i

(A23)

i=1

The area of the triangular element can be written by these parameters. (A10)

r2 dr2 2 2 0 (a1 +a2r2 +a3r 2)



(A15)

i=1

1

1 a2 + 2a3 a − 2I , 8A a1 +a2 +a3 2A 1

(1) = x (2) x (21) i i −xi

2

1 dr2 2 2 0 (a1 +a2r2 +a3r 2)

&

1 {2A − 2a2a3I7 − (a 22 + 2a1a3)I6 − 2a1a2I5 −a 21I4}. a 23 (A14)

a2 2 % (x (21) x (32) ) i i

1



r 42 dr2 2 2 0 (a1 + a2r2 + a3r 2) 1

i=1

(A9)

1 a2 + 2a3 a a − 2 + 3I , 8A a1 +a2 +a3 a1 A 1

1 1 I5 (4a1a3 −a 22)2 2

=

&

"

2

A a a 2 −2a1a3 a +a2 +a3 = 2− 2log 1 + 2 I1 , a3 a3 a1 2a3A 1 1 I4 (4a1a3 −a 22)2 2

&

)2 a1 % (x (21) i

1



!

(32) (32) X 4kpq = x (32) k xp xq

r2 dr2 2 a +a r 0 1 2 2 +a3r 2 1

a + a2 +a3 A a log 1 − 2 I1 , a3 a1 2A



(2) x (32) = x (3) i i −xi

where I1, I2, …, I8 are definite integrals solved analytically and expressed only by the coordinates of the three nodal points of the triangular element.

&

r 32 dr2 2 2 0 (a1 + a2r2 + a3r 2) 1

1 a 22 − 2a1a3 + a1a2 a2(6a1a3 −a 22) I1 , − a1 + a2 + a3 4a3A 8a3A (A13)

1 1 I8 (4a1a3 − a 22)2 2

=

(A12)

and x (32) in Eqs. (A1) to (A3) are defined as the x (21) i i difference of the coordinates of the nodes; i.e.,

1 ( Eklpq {X 1kpqI4(I4 −I5) +X 2kpqI5(I5 −I6) 2 (x (1) l

1 1 I1 (4a1a3 −a 22)2 2



&

a1 + a2 + a3 A log a 23 a1 +

(A2)



1 a2 + 2a1 a1 − + I1 , a1 + a2 + a3 A 8A

237

(A11)

1 1 A= (4a1a3 − a 22)2 4

(A24)

References [1] J.D. Eshelby, in: I.N. Sneddon, R. Hill (Eds.), Progress in Solid Mechanics, North-Holland, Elsevier, 1961, p. 89. [2] T. Mura, T. Mori, Micromechanics, Baifukan, Tokyo, 1976.

238

Y. Nakasone et al. / Materials Science and Engineering A285 (2000) 229–238

[3] T. Mura, Micromechanics of Defects in Solids, second ed., Martinus Nijhoff, Dordrecht, 1987. [4] T. Mura, H. Shoja, H.M. Lin, Y. Safadi, A. Makkawy, Bull. Tech. Univ. Istanbul 47 (1994) 267. [5] G.J. Rodin, J. Mech. Phys. Solids 44 (12) (1996) 1977. [6] H. Nozaki, M. Taya, Trans. ASME: J. Appl. Mech. 64 (1997) 495. [7] H.-B. Li, G.-M. Han, H.A. Mang, Int. J. Num. Methods Eng. 21 (1984) 2071.

[8] R. Yuuki, T. Matsumoto, M. Sato, Trans. JSME 53A (496) (1987) 2428. [9] O.C. Zienkiewicz, K. Morgan, Finite Elements and Approximation, Wiley, New York, 1983. [10] A. Stroud, H.D. Secrest, Gaussian Quadrature Formulas, Englewood Cliffs, NJ, Prentice-Hall, 1966. [11] R. Yuuki (Ed.), Mechanics of Interfaces, Baifukan, Tokyo, 1993.