Micromechanics application to materials-related topics

Tetsu ICHITSUBO

Institute for Materials Research
Tohoku University

Contents

1 Eshelby’s ellipsoidal inclusion theory
 1.1 Eigenstrain
 1.2 Elastic-equilibrium equation for elastic homogeneity
  1.2.1 General solution
  1.2.2 Solution for ellipsoidal inclusion
 1.3 The Eshelby tensor
 1.4 Equivalent inclusion method
2 Phase-field microelasticity theory
 2.1 Macroscopic shape change and total strain
 2.2 Elastic-equilibrium equation for elastic inhomogeneity
 2.3 Application of Discrete Fourier Transformation technique
  2.3.1 Definition of DFT
  2.3.2 On singular point
  2.3.3 Constraint conditions
3 Mechanical energies
 3.1 Definitions and theorem
  3.1.1 Strain energy
  3.1.2 Elastic interaction energy
  3.1.3 Colonnetti’s theorem
  3.1.4 Total mechanical energy and mechanical interaction
 3.2 For ellipsoidal inclusion
  3.2.1 Strain energy
  3.2.2 Mechanical interaction energy
 3.3 External stress applied to elastic inhomogeneity
4 Effective elastic constants
 4.1 Definition of effective elastic constants
 4.2 Strain-concentration factor
  4.2.1 Equivalent inclusion method
  4.2.2 Utilization of Mori-Tanaka’s mean-field theory
  4.2.3 Multi-component composite
 4.3 Application to special composites
  4.3.1 Porous material
  4.3.2 Liquid-inclusion composite
  4.3.3 Rigid-body-inclusion composite
 4.4 Strain-energy balance
 4.5 Effective mean-field (EMF) theory

1 Eshelby’s ellipsoidal inclusion theory

Eshelby established the basis of the micromechanics theory for isotropic body[1]. Afterward it has been extended for anisotropic cases by other researchers, and has become more sophisticated theory; for example, one can see more systematized theory in Mura’s book[2]. In this section, we introduce Eshelby’s ellipsoidal inclusion theory.

1.1 Eigenstrain

Suppose that only a small part Ω in an originally uniform substance is expanded because of thermal dilatation or a certain transformation. Then, the substance becomes in an internal-stress state. If the surroundings of Ω does not exist, the displacement due to such an intrinsic dilatation of Ω will be larger; the strain expressed by this displacement is not elastic but plastic, which is called eigenstrain or stress-free transformation strain, which is denoted as ϵkl*. However, the surroundings prohibits such an intrinsic dilatation; Ω is shrinked elastically because of the constraint by the surroundings. This elastic strain is denoted as ϵkl. When the total displacement measured referring to the initial state is denoted as u = (u1,u2,u3), the total strain can be written as γkl = ∂uk∕∂xl, which is given by sum of the eigenstrain and elastic strain;

γkl = ϵ*kl + ϵkl,                                                        (1)
where ϵkl*0 in Ω, and ϵ kl* = 0 in the exterior matrix D-Ω. This situation is illustrated in Fig. 1.


PIC

Figure 1: Schematic figure indicating the concept of eigenstrain.


Next suppose that a number of inclusions with an eigenstrain ϵkl* exist in a substance; the spatial distribution of ϵkl* is denoted as ϵ kl*(x). When x is located inside Ω (i.e., x Ω), ϵkl*(x) = ϵ kl*, and when x is outside Ω (i.e., x Ω), ϵ kl*(x) = 0. In the case where the elastic constants Cijkl of matrix are equal to those of inclusions, Hooke’s law can be written as

              (              )       (∂u  (x )        )
σij(x) = Cijkl γkl(x ) - ϵ*kl(x) = Cijkl---k----  ϵ*kl(x) ,                 (2)
                                        ∂xl
which holds at each position x.

1.2 Elastic-equilibrium equation for elastic homogeneity

1.2.1 General solution

When the total displacement u or total strain γ are known, all the stress/strain fields can be analyzed. To do this, we have to solve the following elastic-equilibrium equation:

∂σij(x)-=  ∂σi1(x)+  ∂σi2(x)-+ ∂-σi3(x-)=  0.
  ∂xj       ∂x1        ∂x2       ∂x3  (3)

Hereafter, the summation convention is applied when the same subscripts appear. Substituting Eq. (70) into Eq. (3), we obtain

    ∂2uk-(x)         ∂ϵ*k′l′(x-)
Cijkl ∂xjxl  = Cijk′l′  ∂xj   .  (4)

To solve this equation, we transform uk(x) and ϵkl*(x) into the Fourier forms: ũk(g) and ˜ϵ kl*(g). Then, in the Fourier spcace, Eq. (4) is written as

                           *
Cijklgjgl˜uk(g) = - igjCijk′l′˜ϵk′l′(g ).  (5)

The solution of Eq. (5) is given by

˜uk(g) = - iK -ik1gjCijk′l′˜ϵ*k′l′(g) = - iK -ik1 gjCijmn ˜ϵ*mn(g),  (6)

where Kik Cijklgjgl and the subscripts klare changed to mn. When g is denoted as gg (g is the unit vector, and g is the magnitude of the vector),

                  2     ----    2                  ----
Kik = Cijklgjgl = g Cijklgjgl = g Gik.  (Gik ≡ Cijklgjgl)  (7)

Hence, Eq. (6) is written as

             --  -2  -1      *           - 1 - 1-       *
˜uk(g) = - (iggj)g  G ik Cijmn˜ϵmn (g ) = - ig  GikgjCijmn ˜ϵmn(g),  (8)

and hence the expression of ∂uk∂xl in the Fourier space is

{    }
 ∂uk-   = - (igg-)ig -1G- 1g-C    ˜ϵ* (g ) = G - 1 g-g-C    ˜ϵ*  (g),
  ∂xl g         l      ik  j ijmn mn        ik  lj  ijmn mn  (9)

where {A(x)}g means the Fourier transform of A(x). Therefore, the symmetric total strain γkl(x) is expressed as

         1 (∂uk    ∂ul )
γkl(x) = -- ----+  ----                                                (10 )
         2∫ ∂xl    ∂xk
       =     1C    (G -1 g-+ G- 1g-)g-˜ϵ*  (g)exp (ig ⋅ x )-dg-.          (11 )
             2  ijmn   ik l    il  k  jmn               (2π)3
Moreover, since
 *       ∫   *    ′            ′  ′
˜ϵmn(g) =    ϵmn(x )exp (- ig ⋅ x )dx ,  (12)

from Eq. (8)

               ∫  ∫
           -∂--       -2       - 1*    ′               ′    ′-dg---
uk (x) = - ∂xj  g  x′ g CijmnG ik ϵmn (x )exp(ig ⋅ (x - x ))dx (2π)3   (13 )
            ∫  ∫
       = - i      g-1CijmnG - 1 g-ϵ*  (x′)exp(ig ⋅ (x - x′))dx′-dg-,     (14 )
             g  x′          ik  jmn                         (2π)3
and from Eq. (11)
         ∫  ∫
               1-        -1--    -1-- -- *    ′               ′
γkl(x) =  g  x′2 Cijmn(G ik gl + G ilgk)gjϵmn (x )exp(ig ⋅ (x - x ))    (15 )

                                                       dx ′-dg--.      (16 )
                                                          (2π )3
Equations (11), (14) and (16) are the general solutions of the elastic-equilibrium equation (3) for an elastically homogeneous substance.

1.2.2 Solution for ellipsoidal inclusion

Hereafter, we consider the case wehre there is one inclusion Ω that has the eigenstrain ϵ*(x). Since ϵ*(x) = 0 for x′∋ Ω, Eq. (14) is

            ∂  ∫     ∫                                        dg
uk (x) = - ----   dx′   g-2CijmnG -ik1ϵ*mn (x′)exp(ig ⋅ (x - x′))---3,    (17 )
           ∂xj  Ω     g                                      (2π)
                                                                       (18 )
                 ∂  ∫     ∫
     uk (x ) = - ----   dx′   g-2CijmnG -i1k ϵ*mn (x′)exp(- ig ⋅ (x - x′))   (19 )
                ∂xj  Ω     g
                                                              dg
                                                             ----3,    (20 )
                                                             (2π)
because g-2C ijmnGik-1ϵ mn*(x) is an even function of g, and the integral range g-space is invariable when dg →-dg. By using the relation of dg = dg1dg2dg3 = dS(g)dg = g2dS(g)dg, Eq. (20) is expressed as
                    ∫      ∫ ∞ ∫
uk (x) = - --1---∂--   dx ′         CijmnG -1ϵ* (x ′)dS (g)dg,          (21 )
           (2π)3∂xj   Ω     0   |S2|       ik  mn
                                                                       (22 )
                    ∫      ∫ ∞ ∫
u  (x) = - --1---∂--   dx ′         C    G -1ϵ* (x ′)dS (g)dg,          (23 )
  k        (2π)3∂xj   Ω     0   |S2| ijmn  ik  mn

                     ∫     ∫    ∫                                      (24 )
           --1----∂--     ′  0             -1 *    ′    --
       = - (2π)3 ∂x    dx         2 CijmnG ik ϵmn (x )dS (g)dg.         (25 )
                   j  Ω     -∞   |S |
From the upper and lower equations in Eq. (25), we obtain
                      ∫     ∫    ∫
             1     ∂       ′  ∞             -1 *    ′       --       ′
uk (x) = - -----3-----  dx           CijmnG ik ϵmn (x )exp(igg ⋅ (x - x ))(26 )
           2(2π)  ∂xj  Ω     -∞   |S2|                             --
                                                               dS(g)dg,(27 )
where the integral |S2|dS(g) is performed on the surface of the unit sphere S2, with g the vector variable moving on the surface |S2|. By using the characteristics with regard to the delta function
∫  ∞

 - ∞ exp(igx)dg =  2πδ(x),  (28)

we obtain

                   ∫             {∫                          }
u  (x) = - -1---∂--    C     G- 1    ϵ* (x ′)δ(g-⋅ (x - x′))dx ′ dS (g-). (29 )
  k        8π2 ∂xj  |S2|  ijmn  ik   Ω  mn
Since the delta function has a following characteristics
∫  a2
     f(x)δ(x - a)dx =  f(a)   (a1 < a < a2),
 a1  (30)

we have to seek xin Ω so as to satisfy g (x - x) = 0.


PIC

Figure 2: Relation between the ellipsoidal inclusion and unit sphere. The vectors x, xand g are converted to y, yand h through the transformations Eq. (33). When x and g (or y and h) are fixed, and x(or y) is the variable vector, the hatched plane satisfies g (x - x) = hh (y - y) = 0. Namely, δ(hh (y - y)) has a value only when yexists on the hatched plane.


We first deal with the integral Ωϵmn*(x)δ(g (x - x))dxin Eq. (29), where we ragard that x and g are fixed, and consider the case where x = (x1,x2,x3) is located in the ellipsoidal inclusion Ω with radii a1, a2, and a3:

           x21-  x22   x23-
x ∈ Ω :    a21 +  a22 + a23 ≤ 1.  (31)

Transformations of variables are made as follows:

                                    ′  ′   ′
 y =  (x1, x2-, x3),        y ′ = (x-1, x2, x3),                       (32 )
       a1 a2  a3                  a1  a2  a3
        --   --   --         --   h1- h2- h3-
h  = (a1g1,a2g2,a3g3),       h = ( h ,h , h ),                         (33 )
where |h| = 1 and h = |h| =     -- 2      --2      -- 2
{(a1g1) +  (a3 g2) + (a3g3) }12. When the vector variable xmoves inside the ellipsoidal inclusion Ω, ymoves inside the unit sphere S2 (see Fig. 2):
x′ ∈ Ω    ⇐ ⇒     y ′ ∈ S2 :   y2+  y2+ y2 ≤ 1.
                                1    2   3  (34)

Since dx = a1a2a3dy and g (x - x) = hh (y - y), Eq. (29) is written as

                      ∫
           a1a2a3--∂--           -1
uk (x) = -  8π2   ∂xj  |S2|CijmnG ik                                    (35 )
                        { ∫             --             }
                             ϵ*  (x′)δ (h h ⋅ (y - y′))dy ′ dS (g-).       (36 )
                           S2 mn
Here, yis expressed using the cylindrical coordinate (z,r,ϕ); the origin O of the coordinate is the cross point of the vector h and the plane that is normal to h and is passing on the point y, as shown in Fig. 2. Then, the following relations hold:
dy ′ = rdrd ϕdz,

    --
z = h ⋅ (y′ - y),

   --       ′
δ(hh ⋅ (y - y )) = δ(- hz ) = δ(hz).

In Eq. (36), S2ϵmn*(x)δ(hh (y - y))dyis performed under the fixed y and h. Therefore,

           a1a2a3  ∂  ∫
uk (x) = - ----2------    CijmnG -i1k                                    (37 )
            8π    ∂xj  |S2|{∫                           }
                                *    ′      r-              --
                              2 ϵmn(x )δ(hz) hdrdϕd (hz) dS(g)         (38 )
                  ∫          S
                       CijmnG  -1h-1ϵ*                                  (39 )
                   |S2|       ik     mn
                                      {∫  2π   ∫ R    }
                                  -∂--      dϕ     rdr  dS(g),         (40 )
                                  ∂xj    0      0
where the latter formula is obtained under the condition that the eigenstrain is uniform in Ω (i.e., ϵmn*(x) = ϵ mn*). By using the following geometric relation, R2 = 1 - (y h)2 = 1 - (x g)2∕h2 (see Fig. 2), we obtain
   { ∫ 2π   ∫ R    }            2
-∂--     dϕ     rdr  = 2π -∂--R--=  - 2π(x ⋅ g)gjh-2.                 (41 )
∂xj   0      0            ∂xj  2
Substituting Eq. (41) into Eq. (40), we obtain
         a1a2a3 ∫      --          --          --
uk (x) = -------   (x ⋅g)CijmnG -ik1gjϵ*mnh -3dS (g).                     (42 )
           4π    S2
After all, we can express the symmetric total strain as
         (                                            )
           a1a2a3 ∫          - 1--    -1-- -- - 3   --   *
γkl(x) =   -------   Cijmn (Gik gl + G il gk)gjh dS (g)  ϵmn.           (43 )
             8π    S2
It is noted that the total strain in Ω is uniform, being independent of x.

When x is located in the exterior matrix (i.e., x Ω), it is quite complicated to calculate the outside fields because the hatched area displayed in Fig. 2, which satisfies h (y - y) = 0, does not exist for any ywhen y h > 1. Taking into account this matter, Cheng and Mura has derived the outside field of ellipsoidal inclusion[3].

1.3 The Eshelby tensor

As found from Eq. (43), the total strain γik in the inclusion Ω can be written in the following form:

            *
γik = Sikmnϵmn.                                                        (44 )
The tensors Sikmn are called the Eshelby tensor. The Eshelby tensors were first derived for elastically isotropic body[1], which can be calculated analytically. Afterward, the Eshelby tensors for anisotropic cases has been extended by Kinoshita and Mura[4], which requires numerical calculations. The formulae for both cases are summarized in literature[2].

General ellipsoidal shape As seen in Eq. (43), the most general form of the fourth-rank Eshelby tensors for the arbitrary shape of ellipsoidal inclusion are given by

                ∫
         a1a2a3-            -1--     -1-- -- -3    --
Sikmn =    8π     2 Cjlmn (G ijgk + G kjgi)glh  dS (g),                 (45 )
                 |S|
           ----
Gij = Cipjqgpgq,

    {   --        --       --  }1∕2
h =  (a1g1)2 + (a3g2)2 + (a3g3)2   ,

where the surface integral is performed over the unit sphere |S2| that is made by unit vector g, Cjlmn denotes the elastic constants of matrix, gp is a component of g, and (a1, a2, a3) denote the radii of the ellipsoidal inclusion.

Cylindrical shape For cylindrical (needle) shape, i.e., a3∕a1(= a3∕a2) →∞, S tensors reduce to

            ∫
S     =  -1-    C    (G -1g- + G -1 g-)g-dL (g),                        (46 )
  ikmn    4π  SL  jlmn   ij  k     kj i  l
where the line integral is performed over the unit circle SL normal to the x 3 (i.e., a3) direction.

Plate-disc shape For the plate-disc inclusion with zero aspect ratio i.e., a3∕a1(= a3∕a2) 0, S tesors further reduce to

         1           --      -- --
Sikmn =  -Cjlmn (G-ij1gk + G -k1jgi)gl,                                   (47 )
         2
where the vector g is taken to be normal to the disc basal plane.

In general, the principle axes of ellipsoidal inclusion are not always consistent with the crystallographic axes (in which the elastic constants, strain and stress etc are defined). In such a case, it is convenient to transform the elastic constants of crystal-coordinate system into those of the inclusion-coordinate system:

  (I)                (C)
C ijkl = aipajqakralsCpqrs,                                              (48 )
where aij denotes a component of the rotation tensor for the crystalinclusion coordinate transformation.

Moreover, there are often cases that we want to use the Eshelby tensors of 6 × 6 matrix form, because the elastic-constant tensors can be reduced to 6 × 6 matrix form, and the tensor calculation becomes more simple. As found from Eqs. (45)-(46), there is a following symmetry Sijkl = Sjikl = Sijlk and SijklSklij. Therefore, the Eshelby-tensor matrix S can be expressed as

     (  S     S      S                            )
     |   1111    1122   1133                        |
     ||  S2211  S2222  S2233           O            ||
     ||  S3311  S3322  S3333                        ||
S  = ||                     2S2323                 || ,                  (49 )
     |(          O                  2S1313         |)

                                           2S1212
when the matrix has an elastic symmetry higher than orthorhombic system. By using this S matrix, γ = Sϵ* becomes equivalent to the tensor calculation γik = Sikmnϵmn*.

1.4 Equivalent inclusion method

When the elastic constants of inclusion are different from those of matrix, we have to reconsider the elastic-equilibirum equation (this will be discussed in the later section). Here, a convenient method that does not require the reconsideration of the elastic-equilibirum equation is presented. If the elastic constants of inclusion can be somehow regarded as the same as those of matrix by modifying the eigenstrain, we can deal with the same elastic-equilibrium equation as well as in the previous sections. There is a well-known method, so-called equivalent inclusion method, to solve the elastic inhomogeneity problem.

Type-I equivalent inclusion When an ellipsoidal inclusion has no eigenstrain but elastic constants Cijkl different from those of matrix Cijkl, stress and strains due to an applied external stress becomes not uniform because of the elastic inhomogeneity; the stress and strain disturbances, which are denoted as σij and γkl, are tried to be reproduced using an equivalent inclusion.

When there is no inclusion, a uniform strain ϵkl0 is caused by external stress:

  ext        0
σ ij  = Cijklϵkl.                                                        (50 )
When an elastically inhomogeneous inclusion exsits, the external stress is disturbed due to the inclusion. The stress inside the inclusion is expressed as
σexijt + σ′ij = C′ijkl(ϵ0kl + γkl) = Cijkl(ϵ0kl + γkl - ϵ*kl),                  (51 )
            *
γkl = Sklmnϵmn.

We have to seek the equivalent eigenstrain ϵkl* so as to satisfy Eq. (51). Note that if the elastic constants are the same, the stress and strain associated with the external stress are uniform, even when the inclusion has an eigenstrain. Furthermore, the intrinsic (original) eigenstrain ϵkl* of the inclusion does not appear in Eq. (51). This is because the initial internal-stress state due to ϵkl* in the absence of external stress is regarded as a standard. These matters are related with Colonnetti’s theorem, which is described in Sec. 3.

Type-II equivalent inclusion We consider the case where an ellipsoidal inclusion has both eigenstrain ϵkl* and different elastic constants Cijkl. The internal stress due to the inhomogeneous inclusion is tried to be reproduced using an equivalent inclusion.

The stress inside the inclusion is expressed as

                   *
σinijc = C ′ijkl(γkl - ϵ′kl) = Cijkl(γkl - ϵ*k*l),                             (52 )
γ   = S    ϵ** .
 kl    klmn mn

We have to seek the equivalent eigenstrain ϵkl** so as to satisfy Eq. (52).

2 Phase-field microelasticity theory

In the previous section, we derived Eqs. (11), (16) and (14) as solutions of linear elastic-equilibrium equation (3) for elastic homogeneity, and obtained Eq. (43) for the special case where there is one ellipsoidal inclusion in the substance. However, the following two points has not been discussed yet.

  1. In the elastic-equilibrium equation, we did not consider the case where the elastic constants of inclusions are different from those of matrix. As stated in previous section, there is a well-known method, so-called equivalent inclusion method, to solve the elastic inhomogeneity problem. If we do not use the equivalent inclusion method, how do we write the elastic-equilibrium equation, and how should we deal with this problem?
  2. The integrand in Eq. (11), Cijmn(Gik-1g l + Gil-1g k)gj˜ϵ mn*(g) exp(ig x), shows a singularity at g = 0, because Gik = Cipkqgpgq has no inverse matrix. In the derivation of Eq. (43) for ellipsoidal inclusion, we have avoid calculating the integrand at g = 0, by utilizing the characteristics of delta function δ(x). If we treat Eq. (11) more generally, how should we deal with this singularity?

In this section, to overcome the above problems, the phase-field microelasticity theory of Khachaturyan[5] is introduced. The theory adopts Eshelby’s concept, but the ellipsoidal inclusion method is no longer used and the discrete Fourier transformation (DFT) technique is utilized in the theory.

2.1 Macroscopic shape change and total strain

The total strain is redefined as

γkl(x ) = γkl + δγkl(x),  (53)

where γkl and δγkl(x) represent the average and deviation of the total strain, respectively. From the definition,

∫

 V δγkl(x)dx =  0.  (54)

Thus, the macroscopic shape change of the substance is represented by γkl. Since the Eshelby theory treats the disturbances (of stress, strain, displacement fields) caused by the small region with eigenstrain, the shape change of the whole substance is not taken into consideration. Namely, the total strain γkl(x) uesd throughout the previous sections is a quantity measured referring to the sunstance without any macroscopic-shape change. According to Eq. (53), the above can be written as γkl = 0 and γkl(x) = δγkl(x). Also in the case where the macroscopic shape of substance changes, we have to use

δγ  (x ) = ∂uk-,
  kl      ∂xl  (55)

where the displacement u is measured referring to the macroscopical shape with a deformation represented by γkl.

External strain Without macroscopic change due to the eigenstrain, when an external stress σklext, which is given by

       1 ∫
σekxlt=  --   Cklij(x )ϵexijt(x)dx,
       V  V  (56)

is applied to a substance, the macroscopic average γkl of the total strain is equal to the average strain ϵklext cuased by the external stress;

        ∫
γ- =  1-   ϵext(x )dx = ϵext.
 kl   V   V kl          kl  (57)

Internal strain Consider a structural phase transformation, for example, the substance with crystal lattice of cubic symmetry transforms into that of tetragonal symmetry; its eigenstrain ϵkl* is measured referring to the cubic lattice. Then, the macroscopic average of the total strain is

         ∫                ∫
γ-  = -1    γkl(x )dx = -1    (ϵkl(x) + ϵ*(x))dx                         (58 )
  kl  V   V            V   V           kl
         1 ∫                     1 ∫
      = --    C -kl1ij(x )σkl(x)dx +  --   ϵ*kl(x)dx.                         (59 )
        V   V                    V
On the other hand, in the elastic equilibrium (i.e., ∂σik(x)∂xk = 0 holds) and in the absence of external force (i.e., Xi = σiknk = 0 holds, where nk is the normal vector of the surface), the sum of the internal stresses σkl(x) caused by the eigenstrain ϵkl*(x) throughout the whole substance is zero:
∫              ∫                ∫
                                         ∂xj-
    σij(x)dx =     σik(x)δjkdx  =    σik(x)∂x  dx                        (60 )
  V     ∫       V          ∫     V          k
     =      σ (x )n  x dS -     ∂σik(x)x dx =  0.                       (61 )
         |V | ik     k j       V  ∂xk    j
Therefore, when the elastic constants are spatially homogeneous (i.e., Cklij(x) = Cklij), Eq. (59) is writen as
--     1 ∫   *
γkl = --    ϵkl(x )dx.                                                   (62 )
      V   V

External and internal strain The strains in both cases are summed up:

         ∫                ∫               ∫
--    -1     ext         1-    *         1-     -1
γkl = V     ϵkl (x)dx + V    ϵkl(x)dx +  V    Cklij(x)σkl(x)dx.          (63 )
          V

2.2 Elastic-equilibrium equation for elastic inhomogeneity

We first define a new parameter Φ(x) to prescribe which material exists at a position x; Φ(x) = 0 or Φ(x) = 1 mean the matrix or inclusion, respectively. When the difference of the elastic constants between inclusion and matrix is denoted as ΔCijkl, the elastic constants at x are written as

Cijkl(x) = Cijkl + ΔCijkl Φ(x ).  (64)

Similarly, the eigenstrain ϵkl* representing the inclusion is defined as

ϵ*(x) = ϵ* Φ(x ).
 kl      kl  (65)

According to Eq. (70), the internal stress field is given by

 σij(x ) = (Cijkl + ΔCijklΦ (x ))(γkl + δγkl(x) - ϵ*klΦ (x ))                 (66 )

         --                                                            (67 )
  =  Cijklγkl + Cijklδγkl(x) - Cijklϵ*klΦ (x)                               (68 )
                                                                       (69 )
             --
+ ΔCijklΦ (x)γkl + ΔCijklΦ (x )δγkl(x ) - ΔCijklϵ*klΦ2(x )                (70 )
Then, the elastic-equilibrium equation is
      ∂σij(x)-       ∂δγkl(x)-       * ∂Φ-(x)-
  0 =   ∂xj   = Cijkl  ∂xj    - Cijklϵkl ∂xj                            (71 )
                                                       2
+ ΔC    γ- ∂-Φ(x-)+ ΔC     ∂Φ-(x)δγkl(x)-  ΔC    ϵ* ∂Φ-(x-)            (72 )
      ijkl kl ∂xj        ijkl     ∂xj           ijkl kl ∂xj
          ∂2u (x )          ∂Φ (x )            ∂Φ (x)
   = Cijkl---k-----  Cijklϵ*kl------+  ΔCijklγkl-------                  (73 )
          ∂xj ∂xl            ∂xj               ∂xj
                        ∂ (     ∂uk (x))           ∂Φ2 (x)
              + ΔCijkl ----Φ (x)------- - ΔCijklϵ*kl-------.            (74 )
                       ∂xj        ∂xl                ∂xj
Thus, we obtain the elastic-equilibrium equation to be solved for the elastically inhomogeneous substance:
      2                   (            )
C    ∂-uk-(x-)+  ΔC    -∂-- Φ(x )∂uk(x-)                                (75 )
  ijkl ∂xj∂xl       ijkl∂xj         ∂xl
    (                       )                       2
 =   Cijmnϵ*mn - ΔCijmn  γmn  ∂Φ-(x)-+ ΔCijmn ϵ*mn ∂Φ--(x).              (76 )
                              ∂xj                  ∂xj
The approximate solution of this equation has been given by Chen et al.[6].

Zeroth-order approximation When the effect of ΔCijkl is neglected (i.e., ΔCijkl = 0), Eq. (76) reduces to

      2
Cijkl∂-uk-(x-)=  Cijmnϵ*mn ∂Φ-(x),                                       (77 )
      ∂xj∂xl              ∂xj
which is totally the same as Eq. (3). Therefore, from Eq. (8), the solution is given by
        ∫                   dg
uk(x) =    ˜uk(g) exp(ig ⋅ x)-----
                           (2π)3   (78)

˜uk(g) = - ig-1G -ik1gjCijmnϵ*mn ˜Φ(g ).  (79)

Therefore, from Eq. (55), we obtain

          ∫
δ γ (x) =    δ˜γ  (g)exp (ig ⋅ x)-dg--,                                  (80 )
   kl          kl              (2 π)3
where
          1-        -1--    - 1- -- *  ˜
δ ˜γkl(g) ≡ 2 Cijmn(G ik gl + Gilgk )gjϵmn Φ (g ).                          (81 )

First-order approximation We write Eq. (76) as

     ∂2uk (x )   (                   --  )∂Φ (x)
Cijkl--------=   Cijmnϵ*mn -  ΔCijmn γmn  -------                       (82 )
      ∂xj∂xl                              ∂xj
                   *  ∂Φ2(x-)          -∂--(     ∂uk(x-))
         + ΔCijmn ϵmn  ∂x    -  ΔCijmn ∂x   Φ(x ) ∂x     .             (83 )
                          j               j          l
In the Fourier space, Eq. (83) is written as
                  (                    -- )
- Cijklgjgl˜uk(g ) =  Cijmnϵ*mn - ΔCijmn  γmn (igj)˜Φ(g )                  (84 )
                                               {            }
         + ΔCijmn ϵ*mn (igj){Φ2 }g - ΔCijmn (igj) Φ(x )∂uk(x-)  .        (85 )
                                                       ∂xl   g
Since Eq. (85) is nonlinear, it cannot be soloved analytically. Thus, we obtain an approximate solution with the zeroth-order solution. The solution of zeroth-order equation uk(0) is substituted into the right hand side of Eq. (85), and a more accurate solution uk(1) is newly obtained; since C ijklgjgl = g2C ijklgjgl = g2G ik,
                      (
  (1)          -1--  -1 (       *           --  )˜
u˜k (g) = - ig   gjG ik   Cijmn ϵmn - ΔCijmn γmn  Φ (g)                  (86 )
                                            {        (0)   } )
               + ΔC     ϵ*  {Φ2 } -  ΔC      Φ (x)∂u-m-(x)    ,        (87 )
                     ijmn mn      g      ijmn         ∂xn    g
and we obtain the (non-symmetric) total strain as
                   ((                       )
δ ˜γ(1k)l (g) = gjglG -ik1 Cijmn ϵ*mn - ΔCijmn γmn  ˜Φ (g)                     (88 )

                                         {             } )
             + ΔCijmn ϵ*mn{Φ2 }g - ΔCijmn   Φ(x)δγm(0n)(x)    ,           (89 )
                                                        g
where δγmn(0)(x) = ∂u m(0)(x)∂x n.

Higher-order approximation Similarly, we obtain the higher-order solution through the iterative calculations:

                        (
  (N+1 )            --     (                   --  )
u˜k    (g ) = - ig- 1gjG -ik1 Cijmnϵ*mn -  ΔCijmn γmn  ˜Φ (g )               (90 )
                                                               )
                           *     2            {     ∂u (Nm )(x )}
                  +ΔCijmn  ϵmn {Φ  }g - ΔCijmn  Φ (x)--------- g  ,     (91 )
                                                      ∂xn
and
                     (
   (N+1 )      ---- - 1 (      *            --  )˜
δ ˜γkl   (g) = gjglG ik   Cijmnϵmn -  ΔCijmn γmn  Φ(g )                  (92 )
                                            {             } )
                +ΔC      ϵ*  {Φ2 } -  ΔC      Φ (x)δγ(N)(x)    .        (93 )
                     ijmn mn      g      ijmn         mn     g

2.3 Application of Discrete Fourier Transformation technique

2.3.1 Definition of DFT

We try to apply the discrete Fourier transformation (DFT) technique for the calculation of Eq. (80) with Eqs. (81) [0th-order solution], (89) [1st-order solution], and (93) [Higher-order solution]. For the sake of simplicity, as an example, we here treat the solution of zeroth-order approximation, i.e., Eq. (81).


Table 1: Physical quantities corresponding to the indices, xi and gi, used in the discrete Fourier transformation. Note that 2πgi∕NΔx in DFT corresponds to gi = 2π∕λi of conventional Fourier transformation.








Index xi 0 1 2 ⋅⋅⋅ N - 1 (N)
L [m] xiΔx 0 Δx x⋅⋅⋅(N - 1)Δx(NΔx)








Index gi 0 1 2 ⋅⋅⋅ N - 1 (N)
λ [m] NΔgxi- NΔxN-Δ2x⋅⋅⋅ N(NΔ-x1) x)
2π
λ [1/m]2πgi
NΔx 0 -2π-
N Δx -4π-
N Δx⋅⋅⋅ 2(N--1)π
 N Δx (-2π-
Δx)

















In order to perform DFT, the substance is divided into Nd elements, where d denotes the dimension, and the division length is denoted as Δx [m]. Figure 3 indicates the relation between real and reciprocal spaces. With regard to physical quantities A(x) and (g), the definition of DFT is given by

        N∑d
A˜(g) =    A (x)exp (- i2πx ⋅ g ),                                      (94 )
         x              N
            N∑d
A (x) = -1-    A˜(g)exp (i2πx ⋅ g ),                                    (95 )
        N  d g            N
        Nd                    d
δ (g ) = ∑  exp(i2π-x ⋅ g) = N    (g = 0 )                              (96 )
        x       N           0    (g ⁄=  0)
where 2π∕N appears in the exponential term because xi and gi are integer (xi,gi = 0, 1, 2,⋅⋅⋅,N - 1) in the practical DFT calculation. Then, xiΔx and 2πgi∕NΔx correspond to gi(= 2π∕λi) and xi in the conventional Fourier transformation e.g., Eq. (11). Table 1 gives the physical quantities such as wavelength and wavenumber corresponding to the indices, xi and gi, used in the discrete Fourier transformation. Thus, the exponential term can be understood as
                ( 2πgi )          2πxigi   2π
g ⋅ x   =⇒       ------ (xiΔx  ) = -------= ---g ⋅ x.
                 N Δx               N      N  (97)


PIC

Figure 3: Relation between the real and reciprocal spaces. Intensity of O point of the reciprocal spcae is given by summation throughout the real system.


2.3.2 On singular point

We consider the g = 0 point in Eq. (81). The DFT forms of eigenstrain ϵ˜ kl*(g) is given by

 *        * ˜      ∑   *             2π-
˜ϵkl(g) = ϵklΦ (g) =  x ϵklΦ (x) exp(- iN g ⋅ x ).                         (98 )
Hence, at g = 0,
˜ϵ* (0) = ϵ*Φ˜(0) = ∑  ϵ*Φ (x).                                         (99 )
 kl       kl        x  kl
The intensity of origin in the reciprocal space, ˜ϵ kl*(0), equals the sum of ϵ kl*(x) in the real space (see Fig. 3).

In the similar way with this DFT characteristics, the following relation

          ∑
δ ˜γkl(0) =    δγkl(x) = 0                                             (100 )
           x
should be satisfied from Eq. (54). By utilizing Eq. (100), we can avoid calculating the g = 0 point in Eq. (81). Equation (100) is the crucial key for all the calculations.

2.3.3 Constraint conditions

Next let us consider the constraint conditions with regard to the macroscopic-shape change. For this purpose, we rewrite Eq. (63) in the discrete form as

--     1  ∑            1  ∑            1 ∑
γkl = ---d   ϵekxlt (x ) +--d    ϵ*kl(x ) + --d   C -k1lij(x)σkl(x ).          (101 )
      N    x          N    x          N   x
Furthermore, applying DFT formulae and using the characteristics of delta function Eq. (96), the total and elastic strains are expressed as
                      [                 ]
γ-  + δγ  (x) = -1- ∑  γ- δ(g) + δ˜γ  (g ) exp(i2π-g ⋅ x )             (102 )
  kl    kl      N d  g   kl         kl         N
and
         -1- ∑  (--                  *˜    )     2π-
ϵkl(x) = N d     γklδ(g ) + δ˜γkl(g) - ϵklΦ (g) exp(iN g ⋅ x),          (103 )
              g
respectively, where note that ϵkl*(x) = ϵ kl*Φ(x). In the above equations, γ klδ(0) means the constraint condition, and the key is to calculate somehow the average total strain γkl.

(1) Without macroscopic shape change If the substance is fixed in the rigid box, γkl = 0 should be used for the calculations.

(2) Without external stress If no external field is applied and the substance is elastically homogeneous, the average total strain is given by the spatial average of eigenstrain:

  d--    ∑   *         * ˜
N  γkl =    ϵklΦ(x ) = ϵklΦ(0 ).                                      (104 )
          x
The last equal sign is obtained by comparing with Eq. (99). Using Eq. (100) and this equation, we find
--                   *˜
γklδ(0) + δ˜γkl(0) - ϵklΦ (0) = 0,                                     (105 )
for g = 0 in the summation of Eq. (103). After all, the summation in Eq. (103) can be performed excluding the singular (g = 0) point.

If the substance is elastically inhomogeneous, since

         ∑            ∑
N dγkl =    ϵ*klΦ(x ) +   C -k1lij(x)σkl(x ),                             (106 )
          x            x
Eq. (105) does not hold due to the extra term xCklij-1(x)σ kl(x).

(3) With uniform strain Origin of the uniform strain (e.g., mechanical, electrical, magnetical sources) is not limited. Then, γkl is given as a constraint condition.

For example, when the external stress σijext is applied to the elastically homogeneous system, the external uniform strain

 ext       ext    -1  ext
ϵkl (x) = ϵkl = C klijσij  (107)

is caused. Then, when calculating Eq. (103), we use as a constraint condition

--          d--    ∑   ext  ∑   *
γklδ(0) = N  γkl =    ϵkl +    ϵklΦ(x ).                              (108 )
                    x        x

(4) With external stress Since the eigenstrain effect ϵkl*Φ(x) can be added later, we distinguish the effects of the external stress and eigenstress (eigenstrain), and here discuss the only effect of external stress.

This condition is the most general, but the most complicated in the case of the elastic inhomogeneity, which is totally different from the uniform strain condition. Because the uniform strain is usually not caused, when an external stress σklext is applied to the elastically inhomogeneous substance. The relation that necessarily holds is

-ext    1 ∑           ext       1  ∑          -ext    ext
σkl =  --d   Cklij(x)ϵij (x) = --d    Cklij(x )(ϵkl + δ ϵkl (x)),
       N   x                  N    x  (109)

where

          -                       ∑
ϵekxlt(x ) = ϵexklt+  δϵexktl (x)  and        δϵekxtl (x) = 0.
                                   x

Since the eigenstrain effects can be added later, this is now ignored, and then the macroscopic average γkl of the total strain equals the average strain ϵklext cuased by the external stress;

--    -ext   -1- ∑    -1     ext
γkl = ϵkl =  N d    Cklij(x)σ ij (x).
                 x  (110)

Usually, we do not know the detail of σijext(x), and therefore this constraint condition is somewhat complicated. If you know the effective (average) elastic constants Cklij in advance, we obtain

--    -ext   ---1--ext
γkl = ϵkl =  C klijσ ij .  (111)

Therefore, solving this constraint problem is converted to obtaining the effective elastic constants.

The effective elastic constants Cklij obtained by the other methods are generally not self-consistent. Thus, the following iterative calculations are needed.

  1. First assume Cklij arbitrarily (with the other method), and calculate the average total strain using γmn = Cmnij-1σ ijext.
  2. Obtain the equilibrium solution with higher-order approximation:
                                  (            {             } )
δ˜γ(kNl+1 )(g) = - gjglG -ik1ΔCijmn   γmnΦ˜(g) +  Φ (x)δγ(mNn)(x)    ,
                                                          g
    where δγmn(0)(x) = 0 and G ik = Cipkqgpgq. Note that the eigenstrain is not related to the effective elastic constants in the regime of linear elasticity theory.
  3. Check the sum of local stresses using Eq. (70), i.e.,
        ∑               ∑
-1-    σij(x) = -1-   (Cijkl + ΔCijkl Φ(x))(γ +  δγkl(x))
N d  x          N d x                       kl
    If this sum equals the external stress σijext, the assumed C klij is valid, because the sum of internal stresses (eigen-stresses) is zero; see Eq. (61).
  4. After obtaining the self-consistent Cklij, γmn is recalculated and then the effect of the eigenstrain is considered with higher-order approximation.

Remarks This method for evaluating the effective elastic constants is implicitly based on Colonnetti’s theorem, which is described next section. According to this theorem, the external stress does not interact with the elastic strain due to eigenstrain. When the substance includes the elastic inhomogeneous inclusions with an eigenstrain, the substance is in an internal-stress state, regardless of the existence of an external stress. We regard this (equilibrium) state as a standard state, and consider the case where an external stress is further applied to the equilibrium state. Colonnetti’s theorem says that the total strain energy Estr stored in the substance is given by the simple sum of the strain energy due to the eigenstrain Estrint and that due to the external stress Estrext, i.e., E str = Estrint + E strext. If the inclusion has no eigenstrain but has only different elastic constants, the strain energy is caused only by the external stress, which equals Estrext stated above. Namely, we may ignore the eigenstrain effect when considering the effective elastic costants. Once we obtain the accurate or self-consistent effective elastic constants, δγmn is to be calculated with Eq. (93) using the average total strain γmn = Cmnij-1σ ijext.

3 Mechanical energies

In this section, we discuss various kinds of mechanical energies. There are two types of energies to be discussed here. One is the elastic strain energy Estr, and the other is the potential energy Epot of external stress (load). The sum of these energies is called the total mechanical energy:

Emec  = Estr + Epot = (Einsttr + Eexsttr + Ei-ster ) + (Eipnott + Eepxott).          (112 )
Estr consists of three parts: the internal strain energy Estrint caused by the eigenstrain, the energy Estrext by the external stress, and the elastic interaction energy E stri-e between internal stress due to eigenstrain and external stress. According to Colonnetti’s theorem described later, Estri-e = 0 always holds. E pot consists of two parts: Epotint represents the potential energy caused by the interaction between external stress and eigenstrain, and Epotext represents that caused by interaction between the external stress and external strain. Here, we present how to express those energies.

3.1 Definitions and theorem

3.1.1 Strain energy

Let us consider the case where the substance has only one inclusion with eigenstrain ϵij*. The strain energy stored in the whole substance is

       1 ∫                   1 ∫        (              )     1 ∫        (∂ui (x )        )
Estr = --   σij(x)ϵij(x)dx =  --   σij(x ) γij(x) - ϵ*ij(x ) dx = --   σij(x) --------  ϵ*ij(x) dx(113 )
       2  V                 ∫2  V                   ∫        2  V          ∂∫xj
                          1-                      1-   ∂-σij(x-)          1-          *
                        = 2     σij(x )ui(x )njdS - 2      ∂xj   ui(x)dx - 2    σij(x)ϵij(x)dx(114 )
                             |V|                     V                      V
where σijnj = 0 on the surface |V | and ∂σij∂xj = 0 everywhere. Therefore,
           ∫
         1       *
Estr = - 2-   σijϵijdx.                                               (115 )
            V

3.1.2 Elastic interaction energy

Let us consider the case where two inclusions exist in the substance, and we define the elastic interaction energy. The total strain energy due to the eigen-stresses (1) and (2) is given by

       1 ∫  (                )(               )
Estr = --    σ(i1j) (x ) + σ (2ij)(x) ϵ(i1)j (x) + ϵ(i2)j (x) dx                  (116 )
       2∫ V                     ∫
       1-    (1)    (1)         1-    (2)    (2)
    =  2    σij (x)ϵij (x)dx + 2    σij (x)ϵij (x )dx                   (117 )
          V  ∫  (                V            )
          + 1-   σ (1)(x)ϵ(2)(x) + σ(2)(x )ϵ(1)(x ) dx.                   (118 )
            2  V   ij     ij        ij    ij
The third term can be regarded as the elastic interaction energy. Here, the following relation
σ (1ij)(x)ϵ(i2)j (x) = Cijkl(x)ϵ(1k)l (x)ϵ(2ij)(x) = Cklij(x )ϵ(i2j)(x )ϵ(1k)l (x)       (119 )
                                               (2)    (1)
                                            = σij (x )ϵkl (x)         (120 )
holds, so that
 ∫  (                             )
1-   σ (1)(x)ϵ(2)(x) + σ(2)(x )ϵ(1)(x ) dxσ (2)(x)ϵ(1)(x)dx                (121 )
2  V   ij     ij        ij    ij        ij     kl
                    ∫         (   (1)             )
                  =    σ(2)(x ) ∂u-i-(x) - ϵ*(1)(x) dx                (122 )
                     V  ij       ∂xj       ij
       ∫                        ∫     (2)
     =     σ (2)(x)u(1)(x )n  dS -    ∂-σij-(x)u (1)(x)dx                (123 )
         |V|  ij     i      j      V    ∂xj    i
        ∫                        ∫
      -    σ (2ij)(x)ϵ*i(1j)(x)dx =  -    σ(i2j) (x )ϵ*i(j1)(x)dx.                (124 )
          V                       V
After all, we can defeine elastic interaction energy Estr(I) as
         ∫                        ∫
E (sI)tr = -    σ (2ij)(x)ϵ*ij(1)(x)dx =  -    σ(i1j) (x)ϵ*i(j2)(x)dx.                (125 )
           V                       V
To calculate the interaction energy Estr(I), the stress fields σ ij(1)(x) or σ ij(2)(x) outside inclusions have to be calculated. Such an exterior field has been given by Mura and Cheng[3].

3.1.3 Colonnetti’s theorem

When the eigenstress/eigenstrain and the external stress/strain coexist, the elastic strain energy is expressed as

       1 ∫  (               )(               )
Estr = --    σij(x ) + σexijt(x) ϵij(x) + ϵeixjt(x) dx                     (126 )
       2∫ V                   ∫
       1-                   1-    ext    ext
    =  2    σij(x)ϵij(x )dx + 2    σij (x )ϵij (x)dx                     (127 )
          V  ∫  (              V            )
          + 1-   σij(x)ϵext(x) + σext(x)ϵij(x) dx.                     (128 )
            2  V        ij       ij
The third term can be regarded as the interaction energy Estri-e between the external and internal stresses. Here,
       ext                   ext              ext
σij(x)ϵij (x ) = Cijkl(x)ϵkl(x )ϵij (x) = Cklij(x)ϵij (x )ϵkl(x)          (129 )
                           =  σekxlt (x )ϵkl(x) = σeixjt(x )ϵij(x)            (130 )
holds, so that
          ∫
  i-e   1-  (       ext       ext         )         ext
E str =  2     σij(x)ϵij (x ) + σ ij (x)ϵij(x ) dxσij(x)ϵij (x)dx          (131 )
           V                           ∫           ext
                                     =    σ  (x)∂u-i-(x)dx           (132 )
                                        V  ij     ∂xj
               ∫                       ∫  ∂ σ (x)
             =     σij(x )ueixt(x)njdS  -    ---ij----ueixt(x)dx           (133 )
                |V|                     V   ∂xj
eventually this equals 0, because σij(x)nj = 0 and ∂σij(x)∂xj = 0. Thus, external stress and internal stress (eigenstress or eigenstrain) do not interact with each other. This is called Colonnetti’s theorem. Note that there is no limitation on the elastic constants, i.e., the elastic inhomogeneity is taken into consideration in this argument.

Furthermore, using the other hand, σijext(x)ϵ ij(x) = σijext(x)[γ ij(x) -ϵij(x)], in the expansion,

                               ∫  (                           )
                      Ei- e=  1-   σ  (x)ϵext(x) + σext(x)ϵ (x) dx    (134 )
                        str    2  V  ij    ij        ij     ij
                                    ∫         (∂u  (x)        )
                                  =    σexijt(x) ---i----  ϵ*ij(x) dx    (135 )
                                      V          ∂xj
   ∫    ext               ∫  ∂ σeixjt (x )          ∫   ext    *
=      σij (x )ui(x )njdS -   ---------ui(x)dx -    σij (x)ϵij(x)dx    (136 )
    |V|                    V   ∫∂xj             ∫V
                                                    ext    *
                             =     Xiui (x )dS -    σij (x)ϵij(x)dx    (137 )
                                |V|              V
Since Estri-e = 0 from Eq. (133), we obtain
   ∫                  ∫
                          ext    *
-      Xiui(x)dS =  -    σij (x)ϵij(x)dx,                             (138 )
    |V |                V
which is so-called mechanical interaction energy Epotint, discussed next.

3.1.4 Total mechanical energy and mechanical interaction

Let us consider the case where the eigenstrain and external stress coexist. Considering the potential energy due to external stress (potential energy of weight in the system consisting of spring and weight), the total mechanical energy Emec, which corresponds to the Gibbs free energy, can be defined as

          ∫  (               )(               )
Emec  = 1-    σij(x ) + σext(x) ϵij(x ) + ϵext(x) dx                   (139 )
        2  V            ij               ij
                       ∫      (              )
                     -     Xi  ui(x ) + uexit(x) dS,                   (140 )
                        |V|
where the first term means the total strain energy of the eigenstrain and external strain, and the latter term represents the potential energy of external load. Taking notice of Colonnetti’s theorem, this is summarized as
          ∫                  {                  ∫               }
E     = 1-   σ  (x)ϵ (x)dx +   1σext(x)ϵext(x) -     X  uext(x)dS      (141 )
  mec   2  V  ij    ij          2  ij     ij       |V|  i i
                                                   ∫
                                                -      Xiui(x)dS.    (142 )
                                                    |V |
In Eq. (142), the first term is the internal strain energy Estrint caused by the eigenstrain, and the second term, {⋅⋅⋅}, is the total mechanical energy caused by the external stress, where the former term in {⋅⋅⋅} means the external strain energy Estrext, and the latter means the potential energy, Epotext, of the external load, which is expressed as
          ∫                    ∫         ∂uext(x)
Eexptot = -     Xiueixt(x)dS =  -    σeixjt (x )--i-----dx.                 (143 )
           |V |                  V          ∂xj
The third term can be defined as the mechanical interaction energy between external stress and eigenstrain, which is denoted as Epotint in the analogy of E potext. From Eq. (138), Epotint is given by
          ∫                  ∫
Einptot = -     Xiui(x)dS  = -    σeixjt(x)ϵ*ij(x)dx.                      (144 )
           |V |                V

Remarks In the above arguments, there is no limitation or constraint on the elastic constants Cijkl, i.e., the case of elastic inhomogeneity can also be dealt with, and furthermore the shape of inclusions are not specified.

3.2 For ellipsoidal inclusion

When the inclusion shape is assumed to be ellipsoidal, the calculations including volume integrals are fairly reduced, because the internal stress and strain inside the inclusion become uniform.

3.2.1 Strain energy

Elastic homogeneity From Eq. (115), the elastic strain energy is given by

         ∫   1   *
Estr = -     -σijϵijdx.
           V 2
When the inclusion is ellipsoidal, the internal stress inside inclusion, which is expressed using the Eshelby tensor Sklmn as
                  *     *
σij = Cijkl(Sklmn ϵmn - ϵkl),                                          (145 )
is uniform, so that
Estr = - Ω-σijϵ*ij,                                                    (146 )
         2
where Ω is the volume of inclusion.

Elastic inhomogeneity If an elastically inhomogeneous (Cijkl) inclusion has an eigenstrain ϵkl*, using the equivalent eigenstrain ϵkl* (see section 1.4), the internal stress inside the inhomogeneous inclusion is given by

σij = Cijkl(Sklmn ϵ*mn - ϵ*kl),                                          (147 )
and the strain energy is to be calculated as
Estr = - Ω-σijϵ′* = - Ω-Cijkl(Sklmn ϵ*  - ϵ*)ϵ′*.                      (148 )
         2     ij     2            mn    kl  ij

3.2.2 Mechanical interaction energy

From Eq. (144), the mechanical interaction energy Epotint is

          ∫
Eint =  -    σext(x )ϵ* (x)dx.
  pot      V  ij    ij

Elastic homogeneity In this case, the above equation is straightforwardly applicable, i.e.,

Eint =  - Ω σextϵ*.                                                  (149 )
  pot       ij  ij

Elastic inhomogeneity In the case of elastically inhomogeneous (Cijkl) inclusion, if it has no eigenstrain, Epotint is not caused. When the elastically inhomogneneous inclusion has an eigenstrain ϵkl*, E potint should be considered. The applied stress to the inclusion, σijext(x), is different from the originally uniform stress σ ijext(= C ijklϵkl0), i.e., σijext(x) = σ ijext + σ ij(x) σijext. In order to know the stress inside inclusion, Type-I equivalent inclusion method, Eq. (51), is used:

σext(x) = σ ′ext= C ′  (ϵ0 + γ)
  ij         ij     ijkl  kl
Note that ϵkl* does not appear in the equation. This is because the external stress field σijext(x) is not disturbed with the existence of eigenstrain region, in the regime of the linear elasticity theory. This is also a conclusion of Colonnetti’s theorem. After all, in the case where an inclusion possesses both Cijkl and ϵkl*, using σ ijext, the mechanical interaction energy is calculated as
          ∫
  int         ext   ′*              ′ext ′*
E pot = -    σij (x )ϵij(x )dx = - Ωσ ij ϵij                           (150 )
           V       0         **    **  ′*
       = - ΩCijkl(ϵkl + Sklmnϵmn - ϵkl)ϵij.                          (151 )

3.3 External stress applied to elastic inhomogeneity

In the previous sections, the general definitions of the mechanical energies are presented. Here, we focus especially on the case where the external stress is applied to an elastically inhomogeneous substance. We are not concerned about the existence of eigenstrain regions; for the sake of simplicity, the substance has no eigenstrain regions.

In the case where an external stress is applied to an elastically inhomogeneity, the total mechanical energy (the residual parts of Eq. (142) after removing the eigenstrain terms) is expressed as

          ∫                    ∫
        1     ext    ext                ext
Emec  = 2-   σij (x)ϵij (x)dx -     Xiu i (x )dS                      (152 )
           V                     |V|
When σijext(x) and u iext(x) are expressed as σ ijext + σ ij(x) and uiext + u i(x) (σij(x) and ui(x) are the disturbances due to elastically inhomogeneous inclusions), the total mechanical energy is
          ∫                               ∫
        1-     ext    ′  ∂-(uexit-+-u′i)              ext    ′
Emec  = 2    (σij  + σ ij)     ∂xj    dx -      Xi(ui  + u i)dS        (153 )
         ( V∫        ext     ∫            )|V |(  ∫        ′
   = 0 +  1-   σext∂ui--dx -     X  uextdS  +  1-   σext∂u-idx        (154 )
          2  V  ij  ∂xj       |V|  i i        2  V  ij ∂xj
             ∫        ext       ∫        ′      ∫           )
          + 1-  σ ′ij∂ui--dx + 1-   σ ′ij∂u-idx -     Xiu ′idS  .       (155 )
            2 V      ∂xj      2  V     ∂xj       |V|
The first term “0” means the energy in the case where the elastic inhomogeneity only exists without external stress, and the second term represents the total mechanical energy without elastic inhomogeneity. In definition, the third term can be regarded as the mechanical interaction energy (which is denoted as Emecinho) between the external stress and the elastic inhomogeneity:
         1∫      ∂u ′i      1 ∫     ∂uext      1∫      ∂u ′i
Einmheoc =  --  σexijt----dx +  --   σ′ij---i-dx +  --  σ ′ij----dx         (156 )
         2  V    ∂xj       2  V     ∂xj       2∫ V    ∂xj
                                                      ′ext
                                             -     Xiu i dS.         (157 )
                                                |V |
                                                                     (158 )
Since σij∂xj = 0 and σijnj = 0, we can write Emecinho as
         1∫      ∂u ′i      ∫
Einmheoc =  --  σexijt----dx -      Xiu′idS ≡ ΔE  +  ΔV,                  (159 )
         2  V    ∂xj        |V |
where Emecinho was distinguished into the elastic-strain-energy change ΔE and the potential-energy change ΔV . Since ΔE can be written as
         ∫         ′        ∫                  ∫     ext
ΔE   ≡ 1-   σext∂u-idx =  1-   σextu′ n dS - 1-   ∂-σij u′dx         (160 )
       2  V  ij ∂xj       2  |V | ij   i j     2  V  ∂xj   i
                                               1 ∫
                                            =  --   Xiu ′idS,         (161 )
                                               2  |V |
we find the following relation
         ∫
ΔV   ≡ -     Xiu ′idS  = - 2ΔE.                                       (162 )
           |V|
After all, the mechanical interaction energy Emecinho is expressed as
Einmheoc = ΔE   + ΔV  =  - ΔE  = ΔV--.                                  (163 )
                               2

After all, the total mechanical energy is expressed as

        (1 ∫      ∂uext         )  (  ∫                  )
Emec  =  --   σeixjt---i-dx + ΔE   +  -     XiuexitdS +  ΔV             (164 )
         2  V      ∂xj                  |V|
         (1 ∫      ∂uexit         )  (∫                   )
       =  --   σeixjt-----dx + ΔE   -      XiuexitdS +  2ΔE             (165 )
          2  V      ∂xj               |V|
        1 ∫      ∂uext     ∫
Emec  = --   σeixjt--i--dx -     XiuexitdS  - ΔE.                       (166 )
        2  V      ∂xj        |V|

Application of equivalent inclusion The above argument holds regardless of the inclusion shape. We here assume an ellopsoidal-shape inclusion. Referring to Eq. (144), we apply the Type-I equivalent inclusion method to this case:

         ∫                ∫
ΔV   = -     Xiu ′idS  = -    σextϵ*dV  = - Ωσextϵ*,                  (167 )
           |V|              V  ij  ij           ij  ij
where ϵij* is the equivalent eigenstrain that is defined in Eq. (51). Then,
         ΔV--   Ω- ext *
ΔE   = -  2   =  2σij ϵij,                                            (168 )
and
  inho             Ω- ext*
E mec = - ΔE  =  - 2 σij ϵij.                                        (169 )

By substituting Eq. (168) into the first term of the right-hand side of Eq. (165), we can discuss macroscopic effective elastic constants. Next section presents more systematic treatment for deriving the effective elastic constants, and this strain-energy balance will be discussed in Section 4.4.

4 Effective elastic constants

The general formulae for the effective elastic constants are to be derived based on the Eshelby’s ellipsoidal inclusion theory[1] and Mori-Tanaka’s meand-field approximation (MTMF) theory[7]. On this issue, the mathmatical treatment has been first presented by Benveniste[8]. His approach is very simple and widely applied to other relatated topics. In this section, following the procedure of Benveniste, we will present the derivation of the effective elastic constants.

4.1 Definition of effective elastic constants

Consider a composite material consisting of a matrix with the elastic constants C0 and a type of inclusion with C1, whose volume fractions are denoted by 1 - p and p, respectively. The average stress σ and average strain ϵ of the composite are approximated as σ = (1 - p)σ0 + pσ1 and ϵ = (1 - p)ϵ0 + pϵ1, where σ0 = C0ϵ0 and σ1 = C1ϵ1. The average stress σ must be equal to an external stress, because the sum of internal stresses have to be zero: V σij(x)dV = 0. Then, the average elastic constants C of the composite material can be defined as σ = Cϵ, i.e.,

                       --
(1 - p)C0ϵ0 + pC1 ϵ1 = C [(1 - p)ϵ0 + p ϵ1].  (170)

When expressing ϵ1 = Aϵ0, the average elastic constants C of the composite can be written as

--                                     -1
C =  [(1 - p)C0 +  pC1A ][(1 - p )I + pA ] ,  (171)

where I is the unit matrix. The matrix A is called the strain-concentration factor. In order to obtain the effective elastic constants C, we must find A.

4.2 Strain-concentration factor

The strain-concentration factor A can be derived on the basis of Eshelby’s equivalent-inclusion[1] and Mori-Tanaka’s mean-field theories[7].

4.2.1 Equivalent inclusion method

We suppose first that the external stress σext is applied to an elastically uniform (elastic homogeneous) substance:

σext = C0ϵ0,  (172)

where ϵ0 denotes the uniform strain in the absence of any inclusions. Next, we consider a substance containing, in a matrix with the elastic constants C0, an infinitesimally small inclusion that has different elastic constants C1 but has no intrinsic eigenstrain such as a lattice misfit, and consider the case where an external stress σext is applied to such an elastically inhomogeneous substance. Then, the internal stress σ1 inside the inclusion is given by

                           ext    ∞
σ1 = C1 ϵ1 = C1(ϵ0 + γ) = σ   + σ  ,  (173)

where σ and γ represent the disturbances of internal stress and strain. which are caused by the elastic inhomogeneity under an applied stress. Such a situation is depicted in Fig. 4. When the inclusion is ellipsoidal, according to Eshelby’s concept, Eq. (173) can be expressed using the eigenstrain ϵ* of “the equivalent inclusion” as

σ  = σext + σ ∞ = C (ϵ +  γ - ϵ*),
 1                 0  0  (174)

where the extra strain γ is given by γ = Sϵ* with S the Eshelby tensor. Hence, σ can be expressed as

σ∞  = C  (γ - S-1γ ).
        0  (175)

If σ is expressed without using the Eshelby tensor S,

σ∞  = (C1 - C0 )ϵ0 + C1γ =  (C1  - C0 )C -01σext + C1 γ.  (176)

When ϵ0 = 0 or σext = 0, the stress disturbance is not caused i.e., σ = 0 and therefore γ = 0. When ϵ00 or σext0, σ0 and therefore γ = 0 because of SI in Eq. (175). Thus, σ and γ are functions of σext or ϵ 0.

4.2.2 Utilization of Mori-Tanaka’s mean-field theory

We deal with a composite material with a large number of inclusions. Under an external stress σext(σext), the average stresses and strains of matrix and inclusion are denoted as σ0,ϵ0 and σ1,ϵ1, respectively. (Then, σext = (1 - p)C 0ϵ0 + pC1ϵ1 holds.) The differences between σ1 and σ0 and between ϵ1 and ϵ0 are denoted by σand γ, i.e., σ1 = σ0 + σand ϵ1 = ϵ0 + γ.

On the other hand, Mori-Tanaka’s mean-field (MTMF) theory[7] allows us to express the average internal stress of inclusions as

σ1 = σ0 + σ∞.  (177)

Equation (177) can be understood as follows; see Fig. 4. We suppose that new one inclusion is further added into the matrix of the composite with numerous inclusions. As stated in the previous paragraph, in the case of single inclusion, the stress disturbance due to the inclusion embedded in a matrix with the external uniform strain ϵ0 (caused by σext) is given by σ. Since the new inclusion is to be embedded into the arleadly stressed (or strained) matrix with σ0 (or ϵ0), the internal stress of new inclusion is approximately given by sum of σ and σ 0, i.e., σ1 σ0 + σ. Moreover, as far as σ0(x) is not so much deviated from σ0, we may regard that the internal stress σ1 of the newly added inclusion is virtually equal to the average internal stress σ1 of inclusions in the pre-added original composite, i.e., σ1 σ1. Therefore we obtain σ1 σ0 + σ. It should be noted here that we have to use σ in the situation that ϵ 0 equals ϵ0. After all, we find that σ= σ(ϵ 0) and γ= γ(ϵ0).

Then, by using Eq. (175) based on the equivalent inclusion concept, we can write Eq. (177) as

   -       -     ′      -         ′   - 1 ′
C1 ϵ1 = C1 (ϵ0 + γ) = C0 ϵ0 + C0(γ - S  γ ).  (178)

Here, expressing A as

ϵ1 = ϵ0 + γ ′ = A ϵ0,  (179)

and substituting Eq. (179) into eq. (178), we obtain

         -1              - 1
A =  [SC  0 (C1 - C0 ) + I] .  (180)


PIC

Figure 4: Schematic figure indicating the concept of Mori-Tanaka’s mean field theory with the type-I equivalent inclusion.


4.2.3 Multi-component composite

Essentially, the similar treatment is applicable to this case. The species of the component phases are denoted i. Then, the effective elastic constants can be expressed as

--   (∑         )(∑       )-1
C =      piCiAi      piAi    ,
       i           i  (181)

where ipi = 1. The average strain of i phase is denoted as ϵi, and we write ϵi in the similar way to Eq. (179) as

-   -           -
ϵi = ϵ0 + γi = Aiϵ0,  (182)

where Ai represents the strain-concentration factor of i phase, which is given (in the regime of Mori-Tanaka’s mean-field theory) by

Ai = [SC -0 1(Ci -  C0) + I]-1.  (183)

Note that Eq. (183) satisfies A0 = I.

4.3 Application to special composites

4.3.1 Porous material

We here consider the case of porous material; pore can be regarded as inclusion with C1 = 0. Here, the slightly different way of derivation of A is presented. First, let us consider the case where the single pore exists in an infinite matrix. Since σ1 = C1ϵ1 = 0 holds in Eq. (174), the following relation

 ext      ∞              - 1
σ   =  - σ  = - C0 (γ - S  γ )  (184)

must hold. Namely, according to the magnitude of an external stress σext, γ is determined to be

γ = - (S - I)-1SC - 1σext = - (S - I)-1Sϵ0,
                  0  (185)

so as to guarantee that σ1 = 0. As mentioned abobe, one can find that γ is a function of σext or ϵ 0.

Next a porous composite is considered. Since under an external stress σext,

σext′ = (1 - p)C  ϵ + pC  ϵ
                00      1 1   (186)

holds, substituting C1 = 0 into above equation, we obtain σext = (1 - p)C 0ϵ0. In order to refer to the case of an infinitesimal single inclusion, we suppose the case where the average strain ϵ0 of the matrix is equal to ϵ0;

-      -1σext′          -1 ext
ϵ0 = C 0 1---p-= ϵ0 = C 0 σ   ,  (187)

i.e., σext = (1 - p)σext in the case of ϵ 0 = ϵ0. Referring to ϵ0 + γ = Aϵ0 [Eq. (179)], we obtain

      ext′                    ext′           ext′
C -1-σ----- (S - I)-1SC  -1σ-----= AC  -1-σ----,
  0 1 - p                0 1 - p       0 1 - p  (188)

or from Eqs. (179), (185) and ϵ0 = ϵ0, we directly obtain

ϵ - (S -  I)- 1Sϵ =  A ϵ .
 0              0      0  (189)

Eliminating σext or ϵ 0 in both hand sides, after all we obtain

               - 1           -1
A =  I - (S -  I)  S = (I - S)  ,  (190)

which is consistent with Eq. (180) in the case of C1 = 0.

Other derivation of A From the equations,

 ext′         --     --
σ   =  (1 - p )σ0 + p σ1

and

σ- = σ- + σ∞,
 1    0

we obtain

σ- = - pσ∞  + σext′   and    σ- =  (1 - p )σ ∞ + σext′.
 0                            1  (191)

For porous materials, since σ1 = 0,

    --           ∞     ext′                   -1      ext′
0 = σ1 = (1 - p)σ   + σ    = (1 - p)C0 (γ -  S  γ) + σ   .  (192)

Hence,

 ext′                     -1
σ    = - (1 - p)C0(γ - S   γ),  (193)

and from this equation we obtain

                     σext′
γ = - (S - I)-1SC -0 1-----.
                     1 - p  (194)

Therefore, from ϵ0 + γ = Aϵ0 [Eq. (179)],

     σext′                  σext′          σext′
C -01------- (S - I)-1SC  -01------= AC  -01------,
    1 - p                  1 - p         1 - p  (195)

which is the same as Eq. (188).

4.3.2 Liquid-inclusion composite

When the shear modulus of liquid is assumed to be zero, using the Lame constant (λ and μ = 0), the elastic constants of liquid is expressed as

      (                     )
      |  λ  λ   λ           |
      ||  λ  λ   λ     O     ||
      ||  λ  λ   λ           ||
C1  = ||            0        || .                                      (196 )
      |                     |
      (     O         0     )
                          0
As far as liquid inclusion has the above elastic constants, when a strain of any arbitrary forms is applied, a hydrostatic stress/pressure is only yielded, because
(  σ1 )       (  ϵ1)    (  λ(ϵ1 + ϵ2 + ϵ3) )
|  σ  |       |  ϵ |    |  λ(ϵ + ϵ  + ϵ ) |
||   2 ||       ||   2||    ||     1   2    3  ||
||  σ3 ||  = C1 ||  ϵ3||  = ||  λ(ϵ1 + ϵ2 + ϵ3) || .                       (197 )
||  σ4 ||       ||  ϵ4||    ||        0        ||
|(  σ5 |)       |(  ϵ5|)    |(        0        |)
   σ6            ϵ6              0
Thus, even if the strain of inclusion is controlled by the eigenstrain of equivalent inclusion, the internal stress of inclusion can retain a hydrostatic pressure, which means that the equivalent inclusion method is applicable to liquid-phase composites as well as to porous composites.

4.3.3 Rigid-body-inclusion composite

Rigid body can be regarded as inclusion with C1-1 0; even when an external stress is applied to the rigid body, the external strain is not caused, i.e.,

ϵ1 = ϵ0 + γ = 0,  (198)

and therefore γ = -ϵ0. Since C1 →∞ and ϵ1 = 0, σ1 = C1ϵ1 formally has an indeterminate form, and σ1 represents the virtual stress applied to the rigid-body inclusions. This virtual stress is, of course, not a real stress of rigid body, but is the quantity that guarantees the average stress of matrix satisfies (1 -p)σ0 = σext-pσ 1 in the case of composite.

First, we consider the case where a single rigid body exists in an infinite matrix. Therefore, Eq. (174) can be written as

σ  = σext + σ ∞ = C (ϵ -  ϵ + S -1ϵ ) = C S -1ϵ
 1                 0  0    0       0     0     0   (199)

must hold, and therefore

σ∞  = C0 (S-1 - I)ϵ0.  (200)

Next a composite with numerous rigid-body inclusions is considered. In this case, it is convenient to consider the stress-concentration factor, not the strain-concentration factor A, because A 0 from Eq. (198). The average elastic compliances C-1 of the composite material can be generally defined as ϵ = C-1σ, i.e.,

(1 - p)C- 1σ0 + pC -1σ1 = C--1[(1 - p )σ0 + pσ1].
        0          1  (201)

When expressing σ1 = Bσ0, the average elastic compliances C-1 of the composite can be written as

--- 1
C    = [(1 - p)C -01+  pC -11B ][(1 - p)I + pB ]-1,                     (202 )
where the matrix B is called the stress-concentration factor. From Eq. (177) and σ1 = Bσ0, we obtain
 ∞           --
σ   = (B - I)σ0.  (203)

Comparing with Eq. (200) obtained for a rigid body, we find

          -1      - 1
B =  C0(S   -  I)C0  + I,  (204)

where B = S-1 when C 0S-1 = S-1C 0 holds. Since we can regard that C1-1 = 0 for rigid bodies, from Eq. (202), the average elastic constants for composite with rigid-body inclusions are expressed as

C-=  (I + -p---B )C  .
          1 - p    0  (205)

Thus, when p 1, the effective elastic constants are diverged. However, since p∕(1 - p) ~ O(1) when p < 0.99, this divergence takes place steeply in the vicinity of p = 1.

4.4 Strain-energy balance

In the previous arguments, we considered the effective elastic constants on the basis of Hooke’s law, (1 - p)σ0 + pσ1 = C[(1 - p)ϵ0 + pϵ1], with the average stress and strain. Here we consider this matter in terms of the external strain energy caused by an external stress σijext.

When the substance is uniform, the external strain energy is given by

Eexsttr(0)=  1C -1 σextσext.                                             (206 )
          2  ijkl ij   kl
On the other hand, in the case where there is one inhomogeneous inclusion, from Eq. (168), the increase of the external strain energy due to one inclusion is given by ΔE = (Ω2)σijextϵ ij*(), where ϵ kl*() is the equivalent eigenstrain (Type-I) of one isolated inclusion, and depends on both the elastic constants Cijkl(M) of matrix and the external stress σklext. Namely, the external strain energy in the case of one inhomogeneous inclusion is given by
Eexsttr(1)=  Eesxttr(0)+  ΔE  = 1-C -1σextσext+ Ω-σextϵ*ij(∞ ).                (207 )
                         2  ijklij  kl    2  ij
Similarly, when there exist N inclusions in the substance, the external strain energy may be expressed as
  ext(N )    ext(0)      (N ) ~ 1- -1  ext ext   1-∑N  Ω- ext *(i)
E str    = E str   + ΔE     =  2C ijklσij σ kl +  V     2 σkl ϵkl ,       (208 )
                                                i
where ϵkl*(i) means the equivalent eigenstrain of ith inclusion, which is dependent on the sorroundings of the ith inclusion, and is usually different from ϵkl*(j) (ji) and ϵ ij*(). Under a certain condition, ΔE(N) is approximated to be
           1( N Ω         )   p
ΔE  (N) ≈ --  ---σexktl ϵ*k(∞l ) = -σekxltϵ*k(l∞ ),                             (209 )
          V    2              2
where p = NΩ∕V means the volume fraction of inclusions. The condition for satisfying Eq. (209) is that p (or N) is small, and each inclusion is located apart from each other, since this equation does not consider the elastic interaction among inclusions. In this condition, the total strain energy is given by
  ext(N )    ext(0)      (N )   1- -1  ext ext   p- ext *(∞ )
E str    = E str   + ΔE     ≈  2C ijklσij σ kl +  2σ kl ϵkl  .             (210 )
When the macroscopic elastic constants are denoted as Cijkl, the total strain energy is expressed as
  ext   1---1  ext ext
E str =  2C ijklσij σ kl .                                              (211 )
Thus, the following equation
1---1  ext ext   1- - 1 ext ext   p- ext *(∞)
2C ijklσij σ kl  ≈ 2 Cijklσij σkl + 2 σij ϵij  ,                         (212 )
holds from the strain-energy balance. Equation (212) will hold fairly exactly when the fraction of inclusion p is small. By eliminating the external stress σijext in both hand sides of equation, we can obtain the macroscopic elastic constants Cijkl.

In the matrix form, we can rewrite Eq. (212) as

1-ext(T)---1  ext   1- ext(T) - 1 ext   p- ext(T)*(∞ )
2σ     C   σ   =  2σ     C 0 σ    + 2σ      ϵ   ,                   (213 )
where σext(T) means the transposed matrix of σext. Eliminating (12)σext(T),
--- 1 ext     -1  ext     *(∞ )
C   σ    = C 0 σ   +  pϵ   .                                         (214 )
By substituting the following equations,
              --    --             -        -
σext = (1 - p)σ0 + pσ1 = (1 - p)C0 ϵ0 + pC0 (ϵ0 + γ - S-1γ ),

ϵ*(∞ ) = S-1γ,

into the right-hand side of Eq. (214), and taking notice of ϵ1 = ϵ0 + γ and (1 - p)ϵ0 + pϵ1 = ϵ, we obtain

  - 1 ext     *(∞ )     -1          -        -         - 1        -1
C 0 σ    + pϵ    =  C 0 [(1 - p)C0 ϵ0 + pC0 (ϵ0 + γ - S γ )] + pS  γ  (215 )
                                                                     (216 )
              -      -         -1       - 1           -     -
     = (1 - p)ϵ0 + p(ϵ0 + γ - S  γ) + pS   γ = (1 - p)ϵ0 + p(ϵ0 + γ) (217 )
                                                                     (218 )
                                              = (1 - p)ϵ +  pϵ = ϵ.  (219 )
                                                        0     1
After all, we find that Eq. (213) based on the strain-energy balance is completely equivalent to Hooke’s law, σext = Cϵ, i.e.,
       --    --    --       -     -
(1 - p)σ0 + pσ1 =  C [(1 - p)ϵ0 + p ϵ1],
which is the basic definition for deriving the effective elastic constants.

4.5 Effective mean-field (EMF) theory

In general, the MTMF theory does not yield accurate elastic constants for large p, because Eq. (212) does not consider the elastic interaction between inclusions. One of the essential reasons of this problem is that the important conclusion of Mori-Tanaka’s theory, σ1 = σ0 + σ, does not hold when the fraction of inclusions is fairly large, where σ0 and σ1 are the average stress fields of the matrix and inclusions, and σ denotes the stress disturbance caused by a single inclusion formed in an infinite matrix. This is because the elastic-interaction effect between inclusions cannot be neglected any longer. This means that we should avoid deriving directly the elastic constants for a composite with a high inclusion fraction from those of the non-inclusion material with the conventional MTMF theory. With this in mind, we here consider the evaluation method of effective elastic constants for large p.[9]

Equation (212) is rewritten as

1---1              1---1              p
-C ijkl(p )σexijtσekxlt≈  -C ijkl(0)σeixjt σexktl + -σexijtϵ*i(0j),                    (220 )
2                  2                  2
where the superscript (0) means p = 0, Cijkl-1(0) and ϵ ij*(0) represent C ijkl-1 (the elastic compliance of matrix) and ϵij*() in the previous section. Namely, ϵ ij*(0) is the equivalent eigenstrain that is prepared referring to the substance without any inclusions. Equation (220) is fairly accurate when p 1. We consider the case where the fraction of inclusion is further increased, for example, to 2p. According to Eq. (220), the effective elastic constants for 2p is given by
1--                 1--                2p
-C -ij1kl(2p )σeijxtσekxlt ≈ -C -ij1kl(0)σeixjtσekxtl + ---σeijxtϵ*i(j0).                  (221 )
2                   2                   2
However, there is no longer guarantee that Eq. (221) is accurate, because we cannot judge whether 2p 1 or not (even if p 1). If we believe that Cijkl-1(p) obtained for p 1 from Eq. (220) is accurate, we had better use the following equation, instead of Eq. (221):
1---1      ext  ext   1--- 1    ext ext   1---p---ext *(p)
2C ijkl(2p )σij σkl  ≈ 2C ijkl(p)σij σkl + 2 1 - pσij ϵij  ,             (222 )
where ϵij*(p) is the equivalent eigenstrain that is prepared referring to the homogeneous substance with the elastic constants Cijkl-1(p), in which influences of the inclusions are renormalized. The elastic constants Cijkl-1(2p) obtained from Eq. (222) are more reliable than those obtained from Eq. (221), because the condition of p 1 is still used in the equation. Thus, in the case of (N + 1)p, we can write successively
 --                  --
1C -1 (3p )σextσext ≈ 1C - 1(2p)σextσext+ 1---p---σextϵ*(2p)           (223 )
2  ijkl     ij  kl    2  ijkl     ij  kl   2 1 - 2p  ij  ij
                                                                    (224 )
1--                 1--                 1   p
-C -ij1kl(4p )σeijxtσekxlt ≈ -C -ij1kl(3p)σeixjtσekxlt+ --------σexijtϵ*ij(3p)           (225 )
2                   2                   2 1 - 3p
                                                                    (226 )

                                    ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅           (227 )
and after all we obtain the following fomula, using small p,
1---1                      1---1                1   p
-C ijkl((N  + 1)p)σeixjtσekxlt≈  -C ijkl(N p )σexijtσekxlt+  ---------σeijxtϵ*i(jNp), (228 )
2                          2                    21 - N p
                                                                    (229 )
                            --                      N∑
                        ≈  1C -i1jkl(0)σeixjt σexktl + 1σexijt    ---p---ϵ*(ijnp)
                           2                  2    n=0 1 - np
                                                                    (230 )
Since ϵij*(np) is a function of the average elastic constants C ijkl-1(np), it seems difficult to obtain the explicit formula for the effective elastic constants. However, we can calculate the effective elastic constants for high fraction of inclusions successively. Since the influence of inclusions is renormalized into the previous averaged elastic constants in each successive calculation for next higher fraction, we will refer this method as effective mean-field (EMF) method.

As seen in the previous section, the strain-energy balance equation (209) is totally equivalent to Eq. (171) with Eq. (180) derived on the basis of Hooke’s law σext = Cϵ. Therefore, Eq. (171) is still used in the EMF treatment, but is applied only in the low fraction range. Namely, C0 in Eq. (171) is sequentially superseded by C (effective elastic constants) for a lower fraction obtained by a previous calculation step, and then the effective elastic constants of a renewed composite with a slightly higher fraction are to be calculated. When the increment of fraction with regard to the original non-porous material is given by Δp, the division number Q of the entire fraction range (0 p 1) is Q = 1Δp p should be chosen in order for Q to be integer). Then, the following relation,

        Q∑-1             Δp
Q Δp =     (1 - nΔp )---------=  1,
        n=0          1 - nΔp  (231)

holds, where 1 - nΔp means the fraction of the previous (pre-renewed) porous sample and Δp∕(1 - nΔp) represents the virtual fraction measured referring to the previous porous sample that is regarded as a matrix. Therefore, the fraction p of a renewed porous sample is expressed as

                  N
p = (N +  1)Δp =  ∑  (1 - n Δp )--Δp---,    (0 ≤ N ≤  Q - 1 ),
                  n=0          1 - nΔp  (232)

Based on the above concept, Eq. (171) is modified as a following formula:

          (                                     )
--(N+1)     1---(N-+--1)Δp---(N )  ---Δp----
C       =      1 - N Δp    C    + 1 - N Δp  C1A                      (233 )
                              (                               )- 1
                                1---(N-+--1)Δp-    ---Δp----
                                   1 - N Δp    I + 1 - N Δp A     ,  (234 )
where C(N+1) and C(N) are the effective elastic constants of the renewed and previous porous samples, respectively. It should be noted that the strain concentration factor A and the Eshelby tensor S have to be calculated using C(N). Furthermore, C(N) is not always isotropic even if the initial matrix is elastically isotropic. For an anisotropic medium one has to calculate the Eshelby tensor using Eqs. (45)-(47).

References

[1]   Eshelby, J.D., Proc. Roy. Soc. London, A241 (1957), 376.

[2]   Mura, T., Micromechanics of Defects in Solids, Second, Revised, Edition (Martinus Nijhoff Pul., The Hegue, 1987).

[3]   Mura, T. and Cheng, PC., J. Appl. Mech., 44 (1977), 591.

[4]   Kinoshita N. and Mura T., Phys. Status Solidi, A5 (1971), 759.

[5]   Khachaturyan AG, Theory of Structural Transformation in Solids (Wiley, New York, 1983).

[6]   Hu SY. and Chen LQ., Acta Mater. 2001, 49 1879.

[7]   Mori, T. and Tanaka, K., Acta Metal., 21 (1973), 573.

[8]   Benveniste, Y., Mech. Mater., 6 (1987), 147.

[9]   Tane, M. and Ichitsubo, T., Appl. Phys. Lett. 85 (2004), 197.