Mean-field homogenization of multi-phase thermo-elastic composites: a general framework and its validation

COMPOSITES SCIENCE AND TECHNOLOGY Composites Science and Technology 64 (2004) 1587–1603 www.elsevier.com/locate/compscitech

Mean-field homogenization of multi-phase thermo-elastic composites: a general framework and its validation O. Pierard, C. Friebel, I. Doghri

*

Universit
e Catholique de Louvain (UCL), CESAME, B^atiment Euler, 4 Avenue G. Lema^ıtre, B-1348 Louvain-la-Neuve, Belgium Received 4 August 2003; accepted 19 November 2003 Available online 28 January 2004

Abstract This paper deals with mean-field Eshelby-based homogenization techniques for multi-phase composites and focuses on three subjects which in our opinion deserved more attention than they did in the existing literature. Firstly, for two-phase composites, that is when in a given representative volume element all the inclusions have the same material properties, aspect ratio and orientation, an interpolative double inclusion model gives perhaps the best predictions to date for a wide range of volume fractions and stiffness contrasts. Secondly, for multi-phase composites (including two-phase composites with non-aligned inclusions as a special case), direct homogenization schemes might lead to a non-symmetric overall stiffness tensor, while a two-step homogenization procedure gives physically acceptable results. Thirdly, a general procedure allows to formulate the thermo-elastic version of any homogenization model defined by its isothermal strain concentration tensors. For all three subjects, the theory is presented in detail and validated against experimental data or finite element results for numerous composite systems. Ó 2004 Elsevier Ltd. All rights reserved. Keywords: A. Particle-reinforced composites; B. Modelling; B. Thermomechanical properties; C. Computational simulation; Homogenization

1. Introduction We consider composites where the reinforcing phase consists of inclusions in a generic sense, which means fibers (short or long), platelets, spheres or spheroids (i.e., ellipsoids possessing an axis of revolution) with an arbitrary aspect ratio. For those composites, Eshelby [11] based mean-field homogenization models provide a very cost-effective way of predicting the influence of the micro-structure (thermo-mechanical properties of each phase; inclusionsÕ volume fraction, aspect ratio and orientation) on the overall properties. For linear thermoelastic composites, although there exists a large body of literature on mean-field homogenization schemes (see for instance the review article of Tucker and Liang [24]) we devote this paper to three aspects which in our opinion deserve more attention than they did. * Corresponding author. Tel.: +32-10-478042/472350; fax: +32-10472180. E-mail addresses: [email protected] (O. Pierard), [email protected] (C. Friebel), [email protected] (I. Doghri).

0266-3538/$ - see front matter Ó 2004 Elsevier Ltd. All rights reserved. doi:10.1016/j.compscitech.2003.11.009

Firstly, for two-phase composites, that is when the inclusions in a given representative volume element (RVE) all have the same material properties, aspect ratio and orientation, then the interpolative double inclusion model (D-I) proposed by Lielens [17] provides perhaps the best mean-field predictions to date, for a wide range of volume fractions and contrast values. The model is based on the family of D-Is proposed by Nemat-Nasser and Hori [20]. Alternatively, it can also be seen as a non-trivial interpolation between the Hashin– Shtrikman bounds [14]. Secondly, for multi-phase composites, that is when the inclusions in a RVE have different material properties, aspect ratios or orientations (including a twophase composite with randomly oriented or otherwise non-aligned inclusions as a particular case), a commonly used approach is to extend the Mori–Tanaka [19] model (M–T) using BenvenisteÕs observation [2]. However, as shown by Benveniste et al. [3], this direct method might lead to non-symmetric overall stiffness operator, which is physically unacceptable. A solution was advocated by Camacho et al. [7], Lielens [17] and

1588

O. Pierard et al. / Composites Science and Technology 64 (2004) 1587–1603

Friebel [12] and consists in a two-step homogenization procedure. Accordingly, the RVE is decomposed into grains, each one being viewed as a two-phase composite. The first step consists in a homogenization of each grain using M–T or D-I, while in the second step homogenization of the aggregate is performed using Voigt or Reuss schemes. Thirdly, for thermo-elastic problems, instead of extending M–T or other models from the isothermal to the thermal-stress case using ad hoc methods, Lielens [17] has proposed a general framework which allows to extend any homogenization model. Generic formulae are derived and can be used for any mean-field scheme defined by its isothermal strain concentration tensors. In this paper, we give a detailed presentation of the three subjects mentioned above and put a special emphasis on an extensive validation of the predicted results against experimental data and direct finite element (FE) simulations. The paper has the following outline. In Section 2, we recall homogenization models for two-phase composites with a special attention to the interpolative D-I model. Section 3 presents a general two-step procedure for multi-phase composites. In Section 4, generic formulae extending homogenization models from the isothermal to the thermo-elastic case are derived. Section 5 contains a large number of numerical mean-field predictions and their comparison with experimental data or direct FE simulations, for several metal–matrix composites (MMCs) and polymer– matrix composites (PMCs). Conclusions are drawn in Section 6, and the paper also contains Appendix A for M–T model applied to multi-phase composites.

hei ¼ e; hri ¼ r;

1 with hf i
V

Z

f ðx; xÞdV ;

ð2Þ

x

where x is the position vector of a micro particle in a local frame (attached to the RVE). It appears that the problem we would like to solve is to link the strain and stress averages over the RVE. For a two-phase composite, the matrix phase has a volume V0 and a volume fraction (or concentration) v0 ¼ V0 =V where V is the volume of the RVE. We define the same quantities for the reinforcing phase (V1 and v1 ¼ V1 =V ¼ 1  v0 ). It is easy to check that the averages over x (the entire RVE), x0 (the matrix phase) and x1 (the inclusionsÕ phase) are related by hf i ¼ v1 hf ix1 þ ð1  v1 Þhf ix0 :

ð3Þ

2. Homogenization of two-phase isothermal composites In this section, we consider two-phase composites in isothermal linear elasticity. All inclusions have the same shape, aspect ratio, orientation and material properties. A RVE is subjected to linear boundary displacements corresponding to a macro strain e. The strain averages per phase are related through a (still unknown) strain concentration tensor Be : heix1 ¼ B e : heix0 ;

heix1 ¼ Ae : e; 1

Ae ¼ B e : ðv1 Be þ v0 IÞ ;

ð4Þ

where I designates the fourth-order symmetric identity tensor. For any homogenization model defined by Be , the macro stiffness is given by

1.1. Notation

C ¼ ½v1 C 1 : B e þ ð1  v1 ÞC 0  : ½v1 Be þ ð1  v1 ÞI :

Consider a macro material point with position vector x in a fixed Cartesian frame. In a continuum mechanics analysis of the body at macro scale, either the macro strain eðxÞ is given and we have to compute the macro stress rðxÞ, or vice-versa. In isothermal linear elasticity, the macro stress and strain are related via the macro stiffness tensor C,

We recall hereafter some well-known homogenization schemes.

r ¼ C : e;

B e ¼ I;

i:e: rij ¼ C ijkl
lk ;

1

ð1Þ

where throughout the paper, boldface symbols designate tensors – the order of which is clear from the context – and EinsteinÕs convention for repeated indices is used unless otherwise indicated. At the micro-level, a macro-point is viewed as the center of a RVE with domain x and boundary ox. At this level, we can see the heterogeneity of the material. It can be shown (see Nemat-Nasser and Hori [20]) that if linear displacements are applied on ox, then the average strain in the RVE is equal to the macro strain. We have a similar result for the stresses if uniform traction is applied to the boundary,

ð5Þ

2.1. Voigt model Assuming uniform strain in the RVE, the following results are immediately found: C ¼ v1 C 1 þ ð1  v1 ÞC 0 :

ð6Þ

2.2. Reuss model This scheme considers uniform stress inside the RVE. Consequently, the following expressions are found:
 1 1 B e ¼ C 1 C ¼ v1 C 1 : ð7Þ 1 : C 0; 1 þ ð1  v1 ÞC 0 It is known that Voigt and Reuss models give upper and lower bounds for C, respectively. However, those bounds are too far apart and the two models do not take into account neither the shape nor the orientation of the inclusions, therefore they present little interest for two-

O. Pierard et al. / Composites Science and Technology 64 (2004) 1587–1603

phase composites (but they can be useful in the second step of the homogenization for multi-phase composites: see Sections 3.5 and 5.2).

The strain concentration tensor for this model is found to be equal to that of the single inclusion problem: B e ¼ H e ðI; C 0 ; C 1 Þ:

ð8Þ

Consequently, Benveniste [2] proposed the following interpretation: ‘‘each inclusion ðIÞ behaves like an isolated inclusion in the matrix seeing heix0 as a far-field strain’’. 2.4. Double inclusion model The D-I proposed by Nemat-Nasser and Hori [20] supposes that each spheroidal inclusion ðIÞ is wrapped with a hollow inclusion ðI0 Þ of stiffness C 0 . The outer material has a stiffness C R . Inclusions ðIÞ and ðI0 Þ have the same aspect ratio, symmetry axis and center, and their volume ratio equals that of inclusions and matrix in the actual composite ðv1 =v0 Þ. By choosing C R ¼ C 0 , one can show that we retrieve the M–T model: B e ¼ H e ðI; C 0 ; C 1 Þ
Bel :

ð9Þ

Another choice, C R ¼ C 1 leads to another strain concentration tensor: 1


B eu :

ð10Þ

This can be called the inverse M–T model, as it corresponds to M–T for a composite where the material properties of the inclusions and the matrix are permuted. It can be shown that B el and B eu correspond to lower and upper stiffness estimates, respectively, which correspond to the Hashin–Shtrikman bounds (Hashin and Shtrikman [14]). Real behavior is close to Bel for small v1 and close to Beu for large values of v1 . Consequently, Lielens [17] proposed the following interpolative homogenization model: B e ¼ ½ð1  nðv1 ÞÞðBel Þ1 þ nðv1 ÞðBeu Þ1 1 ;

ð11Þ

where nðv1 Þ is a smooth interpolation function which satisfies: dn ðv1 Þ > 0; dv1 lim nðv1 Þ ¼ 1:

nðv1 Þ > 0;

lim nðv1 Þ ¼ 0;

v1 !0

ð12Þ

v1 !1

Lielens proposed the following simple quadratic expression for nðv1 Þ: nðv1 Þ ¼ 12v1 ð1 þ v1 Þ:

bounds. In the remainder of this paper, ‘‘D-I’’ means LielensÕ interpolative model thus defined. 2.5. Self-consistent model

2.3. Mori–Tanaka model

B e ¼ ½H e ðI; C 1 ; C 0 Þ

1589

ð13Þ

LielensÕ model can be seen as a properly chosen interpolation between the M–T and inverse M–T estimates, or equivalently, between the Hashin–Shtrikman

The self-consistent model (S-C) can be defined directly by assuming that each inclusion is embedded in a fictitious homogeneous matrix possessing the compositeÕs unknown stiffness C. It can also be shown that the S-C model belongs to the family of D-I models and is such that the outer material has stiffness C, i.e. C R ¼ C. Implementation of S-C requires an additional iterative loop in order to determine C. Generally, the model gives good predictions for polycrystals but is less satisfying in the case of twophase composites.

3. Homogenization of multi-phase isothermal composites Consider now a RVE x containing a matrix and inclusions of different shapes, orientations or material properties. For multi-phase composites (including twophase composites with non-aligned inclusions as a particular case), a commonly used method in the literature consists in a direct (single step) homogenization based on BenvenisteÕs interpretation of M–T model (e.g., [25]). However, Benveniste et al. [3] proved that the macro stiffness tensor has the required symmetries only if all the inclusions are similarly shaped and have the same orientation. In other cases, direct M–T or D-I cannot be used in general as they might lead to physically unacceptable results. This is why we present in this section a two-step homogenization procedure which was advocated by Camacho et al. [7], Lielens [17] and Friebel [12]. The matrix phase and each inclusion are supposed to be homogeneous, linear elastic and isotropic. The inclusions are classified into N families (i); each family being characterized by given aspect ratio Ari , stiffness C i and orientation distribution function (ODF) Wi ðpÞ. In other words, for each family ðiÞ, the inclusions have identical shape, aspect ratio and stiffness but not necessarily the same orientation. Table 1 summarizes the notation. The volume fractions of the matrix and the families in the RVE obey the following obvious relation: v0 þ

N X

vi ¼ 1:

ð14Þ

i¼1

At this stage, we can consider two cases: (a) each familyÕs ODF Wi ðpÞ represents a differentiable density function wi ðpÞ, or (b) each familyÕs ODF represents a finite number of orientations.

1590

O. Pierard et al. / Composites Science and Technology 64 (2004) 1587–1603

Table 1 Notations used for multi-phase isothermal composites

RVE Matrix Family i

Notation

Volume

Volume fraction

Stiffness

Aspect ratio

x x0 xi

V ðxÞ V ðx0 Þ V ðxi Þ

 v0 ¼ V ðx0 Þ=V ðxÞ vi ¼ V ðxi Þ=V ðxÞ

C C0 Ci

  Ari

3.1. Each family’s ODF is a differentiable density function

The average of a quantity on the RVE can be written as

In this case, the RVE is decomposed into a set of infinitesimal grains. Each grain xi;p , of volume dV ðxi;p Þ, is a two-phase composite containing the matrix material in concentration v0 (the same as in the RVE) and those inclusions (in concentration ð1  v0 Þ) of family i whose orientation is between p and ðp þ dpÞ. The condition of volume conservation for such inclusions gives successively:

hix ¼ ¼ ¼

N 1 X Vð[x Þ hið[xi;k Þ V ðxÞ i¼1 k i;k k Ki N X X V ðxi;k Þ hixi;k V ðxÞ i¼1 k¼1 N X i¼1

Ki X vi hi dWi ðpk Þ: ð1  v0 Þ k¼1 xi;k

ð20Þ

ð1  v0 ÞdV ðxi;p Þ ¼ vi V ðxÞdWi ðpÞ ¼ vi V ðxÞwi ðpÞdp:

ð15Þ

This gives the following identity for the concentration of each grain in the RVE: dV ðxi;p Þ vi ¼ w ðpÞdp: V ðxÞ ð1  v0 Þ i

ð16Þ

The average of a quantity on the RVE can be decomposed into a sum of averages over the grains, N 1 X hix ¼ Vð[x Þ hið[xi;p Þ p V ðxÞ i¼1 p i;p I N X dV ðxi;p Þ hixi;p ¼ V ðxÞ i¼1 I N X vi ¼ ð17Þ hixi;p wi ðpÞdp: ð1  v0 Þ i¼1

3.2. Each family’s ODF represents a finite number of orientations In this case, in each family i, the inclusions span a finite number Ki of orientations. The RVE is decomposed into a discrete set of finite grains. Each grain xi;k of volume V ðxi;k Þ is a two-phase composite with a matrix phase (concentration v0 as in the RVE) and those inclusions (in concentration ð1  v0 Þ) of family i whose direction is pk . Writing the conservation of volume for those inclusions leads to: ð1  v0 ÞV ðxi;k Þ ¼ vi V ðxÞdWi ðpk Þ:

ð18Þ

The concentration of each grain in the RVE is thus given by V ðxi;k Þ vi ¼ dWi ðpk Þ: V ðxÞ ð1  v0 Þ

ð19Þ

3.3. General case We can imagine a RVE where some of the ODFs are differentiable – case (a) – and others are not – case (b). The average over a RVE will then involve a combination of expressions (17) and (20). In any case ((a), (b) or this general case (c)), a RVE-averaged quantity can be cast under the following format: I N X vi hix ¼ hixi;p dWi ðpÞ
hhixi;p ii;Wi ; ð21Þ ð1  v0 Þ i¼1 where the integral symbol can also represent a discrete sum. In the latter equation, hixi;p represents the average over the grain xi;p and hii;Wi the average over all grains, that is over all families i and all orientations. Therefore, homogenization of the RVE is performed in two steps: Step 1: Homogenization of each grain individually. Step 2: Homogenization of the grains. The two-step homogenization procedure is illustrated in Fig. 1. 3.4. First step: homogenization of each grain individually In this first step, each grain is a two-phase composite which possesses a matrix material and a number of inclusions having the same shape, aspect ratio, stiffness and orientation. For this step, any homogenization scheme discussed in Section 2 can be used. Typically, one can use M–T or D-I models. In each grain xi;p , the matrix has a stiffness C 0 and volume fraction v0 , and the inclusions the same stiffness C i , aspect ratio Ari and orientation p (or an orientation between p and p þ dp). The total concentration of inclusions in the grain is ð1  v0 Þ. The homogenization of each grain will give a strain concentration tensor B ei;p and a ‘‘macro’’ (for the grain) stiffness C i;p , defined as in Section 2:

O. Pierard et al. / Composites Science and Technology 64 (2004) 1587–1603

1591

tion of M–T model would be: the average strain in the matrix phase of each grain is the same for all grains, i.e.:

DECOMPOSITION

heix0 ¼ heix0 i;p

8i; p;

ð25Þ

where x0i;p is the matrix phase of the grain xi;p . It is shown in Appendix A that using this assumption in the expression of the average strain and average stress gives the following expression for the macro stiffness of the whole aggregate: FIRST STEP

SECOND STEP

Fig. 1. The RVE is decomposed into a set of grains. Each of them is then individually homogenized (step 1). Finally, a second homogenization is performed over all the grains (step 2).

heix1 ¼ B ei;p : heix0 ; i;p

i;p

hrixi;p ¼ C i;p : heixi;p ;

ð22Þ

where upper indices 0 and 1 refer to the matrix phase and the inclusionsÕ phase in each grain xi;p , respectively. Note that each material being isotropic and the inclusionsÕ orientation being fixed in a given grain, each grain xi;p possesses transverse isotropy about its local spheroid revolution axis p. 3.5. Second step: homogenization over all grains In this second and final step, the RVE is viewed as an aggregate, that is a set of different grains xi;p (Fig. 1). Each grain is characterized by its stiffness C i;p and orientation p, both of which were computed in step 1. The aggregate is thus a multi-phase composite which can be homogenized by one of several homogenization schemes, some of which are recalled hereafter. However, in the end of the section, we draw attention to the fact that some homogenization models for multi-phase composites might lead to physically unacceptable macro predictions. 3.5.1. Voigt The basic assumption of this model is that the strain field is uniform inside the aggregate. We easily obtain that the macro stiffness can be written as C ¼ hC i;p ii;wi :

ð23Þ

3.5.2. Reuss This model assumes that the stress field is uniform in the aggregate. Consequently, the macro stiffness is given by  1 1 C ¼ hðC i;p Þ ii;wi : ð24Þ 3.5.3. Mori–Tanaka For the second step of the procedure (homogenization of the aggregate – Fig. 1), BenvenisteÕs generaliza-

 1 C ¼ ðv0 C 0 þ ð1  v0 ÞhC i : B ei;p ii;wi : ðv0 I þ ð1  v0 ÞhB ei;p ii;wi :

3.5.4. Comments We mentioned in the beginning of this section that we cannot use a direct M–T for all multi-phase composites. On the other hand, it can be shown (e.g., [12,17]) that a direct M–T homogenization on a multi-phase composite is equivalent to a two-step procedure where M–T is used in each step. Consequently, M–T model should not be used in the second step of the procedure advocated in this paper, because it might lead to physically unacceptable macro predictions, as illustrated in Section 5.2. One can use M–T or D-I for the first step since each grain is a two-phase composite with inclusions of identical aspect ratio, stiffness and orientation. But for the second step, when the grains have different aspect ratios or orientations, one should use other homogenization models which are suitable for aggregates or polycrystals (e.g., Voigt or Reuss).

4. Homogenization of two-phase thermo-elastic composites Consider a homogeneous material of stiffness C. For a given total strain eðxÞ and a change in temperature hðxÞ, the stress is given by rðxÞ ¼ C : eðxÞ þ bhðxÞ;

b ¼ C : a; ;

ð26Þ

where a is the thermal expansion tensor. In the isotropic case, aij ¼ adij , with dij being KroneckerÕs symbol. Consider a two-phase composite, where the matrix has uniform properties C 0 and a0 , and the inclusions have the same aspect ratio, orientation and properties C 1 and a1 . Relations (26) hold for each phase; our aim is to write similar macro relations for the whole composite. In order to find the macro properties C and a, one can re-derive a homogenization model taking into account thermo-elastic behavior instead of isothermal elasticity. A better alternative is proposed in Lielens [17]. Its major interest is that given any homogenization model which is defined in the isothermal case by its strain concentration tensor B e (or Ae ), general expressions of the macro thermo-elastic properties can be found. Assume that the composite is subjected to linear displacements on the boundary (corresponding to a macro total strain e ¼ hei)

1592

O. Pierard et al. / Composites Science and Technology 64 (2004) 1587–1603

and uniform temperature change h. The proof is based on the following three-step approach. 4.1. First step In this step, the composite is subjected to linear boundary displacements corresponding to the final total strain es1 ¼ e and to zero change in temperature hs1 ¼ 0. This step corresponds to a classical isothermal transformation whose solution is given in Section 2. The strain averages are found as follows: hes1 ix1 ¼ Ae : e;

v1 hes1 ix1 þ v0 hes1 ix0 ¼ e;

ð27Þ

where Ae is given by Eq. (4). The stress averages are given by hrs1 ix1 ¼ C 1 : hes1 ix1 ;

hrs1 ix0 ¼ C 0 : hes1 ix0 :

ð28Þ

heix1 ¼ Ae : ðe þ Des3 Þ þ Des2 ; v1 heix1 þ v0 heix0 ¼ e þ Des2 þ Des3 ;

ð33Þ

where Des2 is given by Eq. (30). At the end of the three steps, the RVE is subjected to a uniform temperature change h and linear boundary displacements corresponding to a macro total strain e þ Des2 þ Des3 . However, the latter value should be equal to e, therefore we have Des3 ¼ Des2 :

ð34Þ

Consequently, the per-phase strain averages can be computed from Eqs. (30) and (33), the per-phase stress averages from Eqs. (28), (29) and (32), and the stress average over the RVE by superposition,

    hri ¼ Drs2 þ v1 hrs1 ix1 þ hDrs3 ix1 þ v0 hrs1 ix0 þ hDrs3 ix0 :

ð35Þ

4.2. Second step 4.5. Final expressions

In this step, the RVE witnesses a uniform temperature increment equal to the final temperature change Dhs2 ¼ h, and incremental boundary linear displacements corresponding to a macro strain increment Des2 are imposed so that uniform strain and stress increments are obtained:

Carrying out the computations and rearranging terms, the following expressions for the total strain averages are found:

Drs2 ¼ C 0 : Des2 þ b0 h ¼ C 1 : Des2 þ b1 h:

ae
ðAe  IÞ : ðC 1  C 0 Þ

ð29Þ

This allows to compute the uniform total strain increment: s2

De ¼ ðC 1  C 0 Þ

1

: ðb1  b0 Þh:

ð30Þ

1

: ðb1  b0 Þ;

ð36Þ

v1 heix1 þ v0 heix0 ¼ e; where Ae is defined in Eq. (4). Similarly to (26), it is found that the macro thermo-elastic response is written under the following format: hri ¼ C : e þ bh;

4.3. Third step In this step, the composite is subjected to linear boundary displacements corresponding to a macro total strain increment Des3 and to zero temperature increment Dhs3 ¼ 0. This is a classical isothermal transformation whose solution is given in Section 2 as follows: hDes3 ix1 ¼ Ae : Des3 ;

v1 hDes3 ix1 þ v0 hDes3 ix0 ¼ Des3 ; ð31Þ

where Ae is given by Eq. (4). The per-phase stress averages are given by s3

heix1 ¼ Ae : e þ ae h;

s3

hDr ix1 ¼ C 1 : hDe ix1 ;

s3

b ¼ v0 b0 þ v1 b1 þ v1 ðC 1  C 0 Þ : ae ;

ð37Þ

where the macro stiffness C is given by the isothermal expression (5), ae by Eq. (36) and the macro thermal expansion a is defined as follows: a ¼ C

1

: b:

ð38Þ

In conclusion, generic expressions for thermo-elastic properties were obtained for any homogenization model defined in the isothermal case by its strain concentration tensor B e (or equivalently Ae ). 4.6. Special case

s3

hDr ix0 ¼ C 0 : hDe ix0 : ð32Þ

4.4. Superposition Using the superposition theorem, it is seen that at the end of the three steps, the per-phase strain averages are given by

When the two phases have identical stiffness operators (C 0 ¼ C 1
C), Eqs. (36) and (37) become invalid. Lielens [17] suggests the following solution. In the previous three-step method, step 1 remains unchanged and steps 2 and 3 are replaced by one single step in which the RVE is subjected to a uniform temperature change hs2 ¼ h and zero boundary displacements corresponding to zero macro total strain Des2 ¼ 0. The average stress is then

O. Pierard et al. / Composites Science and Technology 64 (2004) 1587–1603

hDrs2 i ¼

1 X

5.1. Two-phase isothermal composites

vr hC : Des2 þ br hixr

r¼0 1 X

s2

¼ C : De þ

! vr br h ¼ ðv0 b0 þ v1 b1 Þh:

r¼0

ð39Þ Superposition of steps 1 and 2 shows then that: hri ¼ C : e þ b h;

1593

b ¼ v0 b0 þ v1 b1 ;

ð40Þ

where the macro stiffness C is given by the isothermal expression (5).

5. Numerical simulations and their validation In this section, we compare our predictions of the macroscopic behavior with experimental data and FE results obtained either on a RVE or on a unit cell. Note that in the latter case, predictions may depend on the unit cell considered in the simulation. Several authors have already discussed this subject, e.g. B€ ohm and Han [5] for two-phase isothermal composites and Weissenbek et al. [26] for two-phase thermoelasto(plastic) composites. Tucker and Liang [24] also made an evaluation of several homogenization schemes and different FE calculations based on different unit cells.

5.1.1. First test A polymer matrix made of epoxy (E0 ¼ 3:16 GPa, m0 ¼ 0:35) is reinforced with spherical silica inclusions (E1 ¼ 73:1 GPa, m1 ¼ 0:18). The contrast (i.e., stiffness ratio E1 =E0 ) is 23:13. Fig. 2 shows predictions of macro elastic modulus obtained with various homogenization models and with unit cell FE computations (described in Appendix A in Doghri and Ouaar [10]). It is seen that the D-I model gives a good agreement with FE, and that the difference between M–T and D-I predictions for high values of v1 is rather pronounced. For v1 ¼ 0:389, Table 2 gives results for macro values of elastic modulus (E), PoissonÕs ratio (m) and shear modulus (G), for different models: M–T ¼ HSL (Hashin– Shtrikman lower bound), HSU (Hashin–Shtrikman upper bound) and D-I (interpolative double inclusion). On the other hand, B€ ohm [4] conducted FE simulations on a RVE (not a unit cell) using the commercial package Palmyra [21] and considering four cases: 50 or 80 randomly positioned monodispersely sized spheres, 80 randomly positioned bi-dispersed or polydispersely sized spheres. The results are reported in Table 2 together with the predictions of homogenization models or bounds. The table shows that the predictions of the interpolative D-I model are very close to the Palmyra FE results.

PMC with spherical inclusions. Macro elastic stiffness with different models. E0=3.16GPa, ν0=0.35, E1=73.1GPa, ν1=0.18 25

Macro elastic modulus/E0

Voigt Hashin-Shtrikman upper bound

20

interpolative double inclusion Mori-Tanaka = HSL Reuss 15

Finite elements (ABAQUS)

10

5

0 0

0.2

0.4

0.6

0.8

1

Volume fraction of inclusions Fig. 2. PMC: macro elastic modulus predictions with different models: Voigt (uniform strain), Hashin–Shtrikman upper bound (HSU), interpolative double inclusion, Mori–Tanaka ¼ Hashin–Shtrikman lower bound (HSL), Reuss (uniform stress) and finite elements (unit cell).

1594

O. Pierard et al. / Composites Science and Technology 64 (2004) 1587–1603

of that stress. This in turn explains the increase of the axial macro elastic modulus for h > 60 : fibers maintain the deformation of the matrix in their transverse direction and the matrix is then transversally reinforced by the fibers. Of course, Voigt and Reuss homogenization models are obvious in an off-axis tensile test as they donÕt take into account the orientation of the inclusions. In Fig. 3, this is illustrated by a constant value for all the off-axis angles. The difference between M–T and D-I predictions remains small due to the low value of the volume fraction of the inclusions. We can thus conclude that homogenization schemes catch well the critical angle of an off-axis tension test and the comparison of the predictions against unit cell FE results are excellent.

Table 2 PMC: YoungÕs modulus E, PoissonÕs ratio m and shear modulus G obtained with homogenization schemes or bounds (this work) and Palmyra FE results [4] M–T ¼ HSL E (Gpa) m G (Gpa)

HSU

6.8028 0.3146 2.5875

DI

20.89 0.2247 8.53

FE

7.7979 0.3071 2.9828

7.89–8.06 0.287–0.292 3.03–3.11

5.1.2. Second test We now consider an aluminum alloy matrix (E0 ¼ 70 GPa, m0 ¼ 0:33) reinforced by short fibers of alumina (dAl2 O3 , E1 ¼ 300 GPa, m1 ¼ 0:2, Ar ¼ 20). We simulate an off-axis tensile test in which the direction of the aligned fibers presents an angle h with the direction of traction. Surprisingly, Fig. 3 shows that the axial macro elastic modulus does not decrease monotonically when increasing h but presents a minimum around h ¼ 60 and not 90 . When comparing our numerical predictions with unit cell (3D cylinders) FE results of Kang and Gao [13], we find a good agreement between the two methods and also a minimum located around the same angle, which confirms the special shape of the curve. The presence of this minimum is explained physically by Kang and Gao [13]. When h 6 60 , increasing h decreases the axial (and tensile) stress in the fibers, weakening the reinforcementsÕ effect. But when h > 60 , the axial stress in the fibers becomes compressive and increasing h induces a rise in the absolute value

5.1.3. Third test Consider the same composite as in the second test (aluminum alloy matrix reinforced by short fibers of alumina). In order to evaluate the influence of the aspect ratio Ar of reinforcements on the overall behavior, we consider Ar varying from 5 to 50. Our predictions with homogenization models are compared with unit cell (same as in second test) FE results performed by Kang and Gao [13]. As expected, Fig. 4 shows that the transverse YoungÕs modulus is not influenced by Ar and the predicted values agree with FE results. On the contrary, the longitudinal YoungÕs modulus is much more

Aluminum alloy matrix with δ-Al2O3 alumina short fibers. E0=70GPa, ν0=0.33, E1=300GPa, ν1=0.2, v1=10%, Ar=20 Voigt

Axial macro elastic modulus / E0

1.35

1.3 Homogenization schemes Finite Elements (Kang and Gao)

1.25

1.2 Double inclusion 1.15 Mori-Tanaka 1.1

Reuss 0

10

20

30

40

50

60

70

80

90

Off-axis angle θ°

Fig. 3. Aluminum alloy matrix composite: off-axis tensile test. Prediction of the axial macro elastic modulus for various orientations of the aligned fibers. Comparison between predictions of homogenization models and FE results performed on a 3D unit cell.

O. Pierard et al. / Composites Science and Technology 64 (2004) 1587–1603

1595

Aluminum alloy matrix with δ-Al2O3 alumina short fibers. Longitudinal and transverse tension tests. E0=70GPa, ν0=0.33, E1=300GPa, ν1=0.2, v1=10% 1.34

Macro elastic modulus / E0

1.32

Longitudinal

1.3 1.28 1.26 1.24

Double-Inclusion

1.22

Mori-Tanaka Finite Elements (Kang and Gao)

1.2 1.18 1.16 Transverse 1.14 1.12 5

10

15

20

25

30

35

40

45

50

Fibers’ aspect ratio Fig. 4. Aluminum alloy matrix composite: influence of the aspect ratio of the reinforcing fibers on the macro longitudinal and transverse elastic moduli. Comparison between the double inclusion model, Mori–Tanaka model and FE results on a 3D unit cell.

dependent on Ar . For an increasing value of Ar , the longitudinal macro elastic modulus increases up to a plateau. FE results show the same behavior but the difference with mean-field predictions is more pronounced. One reason might be that in the FE calculation fibers were modeled by cylinders whereas homogenization schemes consider spheroids. 5.1.4. Fourth test We now consider a composite made of a polypropylene matrix (E0 ¼ 1 GPa, m0 ¼ 0:35) reinforced with talc (E1 ¼ 100 GPa, m1 ¼ 0:25) platelets (i.e., the axis of revolution of the spheroid is the smallest one). Numerical simulations are compared with results from van Es [25] performed with the FE code Palmyra [21] on a RVE. Results are illustrated in Fig. 5 for aspect ratios ranging from 0.01 to 1 (spherical inclusions) and for several volume fractions (5%, 10% and 15%). Because of the low volume fraction of inclusions, the difference between M–T and D-I models remains small with slightly better results obtained with the latter model. The elastic modulus in the direction of the revolution axis (Ek ) is rather independent of the Ar of inclusions, which logically contrasts with the fibers in previous tests. On the opposite, the elastic modulus in the direction of the radius of the spheroid (E? ) is strongly dependent on the aspect ratio. This modulus increases with decreasing values of the aspect ratio. The correspondence between homogenization predictions and PalmyraÕs results is rather excellent, either for Ek or E? .

5.2. Multi-phase isothermal composites 5.2.1. First test A three-phase composite made of Ti3 Al matrix (E0 ¼ 95:6 GPa, m0 ¼ 0:2884, v0 ¼ 0:55) is reinforced with two kinds of constituents. The first ones are SiC fibers (E1 ¼ 431:0 GPa, m1 ¼ 0:2529, v1 ¼ 0:20, Ar1 ¼ 100) and the second ones are carbon circular disks (E2 ¼ 34:4 GPa, m2 ¼ 0:2028, v2 ¼ 0:25, Ar2 ¼ 0:01). The revolution axes of the platelets and the fibers are oriented along the traction direction. We thus have a three-phase composite with two reinforcements of different shapes. As mentioned in Section 3, a direct M–T homogenization scheme should be avoided in this case because it leads to a non symmetric homogenized stiffness tensor C. We can circumvent this problem by using the two-step homogenization procedure with a Voigt or Reuss method for the second step. Li [15] suggested another method in order to obtain a symmetric C. His so-called regularization method is based on a combination of M–T and self-consistent (S-C) methods and can be explained as follows: under uniform strain boundary condition, assume that the average strain in reinforcement i equals the strain of a single inclusion with elastic stiffness tensor C i embedded in an infinite matrix with the yet to be determined effective stiffness tensor C (effective medium assumption – S-C), and subjected to a uniform strain at the boundary, which equals the yet unknown average strain e0 in the matrix (effective field assumption – M–T).

1596

O. Pierard et al. / Composites Science and Technology 64 (2004) 1587–1603

Polypropylene matrix with Talc platelets. E0=1GPa, ν0=0.35, E1=100GPa, ν1=0.25 8 PALMYRA results (van Es 2001) v1 = 5 %

Macro elastic modulus / E0

7

v1 = 10 % 6

v1 = 15 % double inclusion

5

4

E⊥

3

2

E//

1 0.01

0.1

1

Aspect ratio Fig. 5. Influence of the shape of the reinforcements (platelets) on the macro elastic modulus in the direction of the revolution axis (Ek ) and along the radius (E? ) of the spheroids. Comparison between the double inclusion model and FE results (Palmyra).

In Table 3, we compare several homogenization schemes. Three are based on our two-step procedure with M–T for the first step. Another one uses the direct S-C method (performed by Li) and the last one is LiÕs regularization method explained above. As expected, only M–T/Voigt, M–T/Reuss and LiÕs regularization show the required symmetries (e.g., C 1133 ¼ C 3311 ). However, it is difficult to compare quantitatively the corresponding three sets of results since we have no FE simulations or experimental data. 5.2.2. Second test Consider again a three-phase composite reinforced with two kinds of inclusions, but this time, they differ only by their shape (Ar1 ¼ 50, Ar2 ¼ 1) and not by the material (E1 ¼ E2 ¼ 50 GPa and m1 ¼ m2 ¼ 0:3). Results with the direct M–T approach [22] are plotted together with those of the two-step homogenization scheme in Fig. 6. We can observe the equivalence between the direct M–T homogenization and the two-step

procedure using M–T at the two steps (a proof of the equivalence is given by Friebel [12]). An advantage of the two-step procedure is that in the second step we can easily switch to other homogenization schemes which lead to acceptable results. Once again it is seen that even when considering the same material for the two kinds of inclusions, the M–T/M–T prediction is out of the range of M–T/Voigt and M–T/Reuss. 5.2.3. Third test We now consider another kind of three-phase composite where reinforcements have different YoungÕs moduli (E1 ¼ 50 GPa and E2 ¼ 100 GPa), different shapes (fibers or spheres) and are aligned with the direction of traction. Predictions of the two-step homogenization method are compared to Taya and Chou [22] direct M–T homogenization results. As plotted in Fig. 7, we are not surprised that the homogenized longitudinal elastic modulus predicted by M–T/M–T (or equivalently Taya

Table 3 Components (in GPa) of the macro stiffness obtained with several homogenization methods for a three-phase composite reinforced by two kinds of inclusions (different materials and shapes) with revolution axis along direction 3 Approach

C 1111

C 1122

C 1133

C 3311

C 3333

C 1212

C 3131

M–T/M–T M–T/Voigt M–T/Reuss LiÕs self-consistent [15] LiÕs regularization [15]

121.06 136.19 110.16 124.0 128.1

41.43 49.96 37.94 44.25 47.24

34.27 45.77 30.49 36.22 35.10

21.84 45.77 30.49 22.34 35.10

123.96 162.47 95.03 124.26 125.47

39.81 43.14 36.11 32.69 33.32

32.99 41.98 30.95 39.86 40.45

O. Pierard et al. / Composites Science and Technology 64 (2004) 1587–1603

1597

3-phase composite with 2 kinds of diferently shaped inclusions. E0=1GPa, ν0=0.35, E1=E2=50GPa, ν1=ν2=0.3, Ar =50, Ar =1 1

Longitudinal macro elastic modulus / E0

16

2

Mori-Tanaka, Voigt

14

Mori-Tanaka, Mori-Tanaka 12

Mori-Tanaka, Reuss Taya and Chou results

10

8

6

4

2

0 0

0.05

0.1

0.15

0.2

0.25

0.3

0.35

0.4

0.45

0.5

Total volume fraction of inclusions v1 + v2 ( v1 = v2 ) Fig. 6. Prediction of the longitudinal macro elastic modulus with several homogenization methods on a three-phase composite with different shapes for the inclusions (made of the same material and with the same orientation).

3-phase composite with 2 kinds of aligned inclusions. E0=1GPa, ν0=0.35, E1=50GPa, E2=100GPa, ν1=ν2=0.3 Longitudinal macro elastic modulus / E0

35

30

Mori-Tanaka, Voigt Ar =10, Ar =50 1

Mori-Tanaka, Mori-Tanaka 25

2

Mori-Tanaka, Reuss Taya and Chou results

20 Ar =50, Ar =10

15

1

10

2

Ar = 1, Ar =10 1

2

5

0 0

0.1

0.2

0.3

0.4

0.5

0.6

Total volume fraction of inclusions v1+v2 (v2=3v1) Fig. 7. Prediction of the longitudinal macro elastic modulus with several homogenization methods for a three-phase composite with different shapes and materials for the inclusions (with the same orientation).

and ChouÕs direct M–T) is not included between the M– T/Voigt and M–T/Reuss curves. The finding is the same as previously: M–T/M–T leads to unacceptable predictions for differently shaped inclusions in a composite.

5.2.4. Fourth test Consider now a three-phase composite with two kinds of fibers which are identically shaped (Ar1 ¼ Ar2 ¼ 10), oriented along the direction of traction but made of

1598

O. Pierard et al. / Composites Science and Technology 64 (2004) 1587–1603

different materials (E1 ¼ 10 GPa and E2 ¼ 50 GPa). On this composite, we compute predictions with our twostep homogenization scheme with three different methods: M–T/Voigt, M–T/Reuss and M–T/M–T. Note that with this kind of composite, conditions formulated by Benveniste et al. [3] in order to obtain the required symmetries on the macro stiffness tensor are respected. In Fig. 8, we plotted the results of these homogenization schemes for different volume fractions of the inclusions. It is seen that once the required symmetry conditions are satisfied, the M–T/M–T results become included between M–T/Voigt and M–T/Reuss predictions, which is more logical. 5.2.5. Fifth test We consider again a three-phase composite with two kinds of inclusions which are identically shaped (long fibers, Ar1 ¼ Ar2 ¼ 1) and oriented along the direction of traction. They only differ by the material (E1 ¼ 750, E2 ¼ 1000, m1 ¼ m2 ¼ 0:3 and v1 ¼ v2 ¼ 0:1). Material properties of the matrix are: E0 ¼ 250 and m0 ¼ 0:3. We predict the transverse macro elastic modulus with various homogenization schemes and compare these results with a FE simulation on a RVE [8]. FE simulations are performed with ABAQUS [1] assuming generalized plane strain. The mesh was generated with the software TRIANGLE [23] which creates a Delaunay triangulation on a cross-section of the RVE as illustrated in Fig. 9. In Table 4, we compare the FE results of the normalized transverse macro elastic moduli with predictions

Fig. 9. Nodes of the mesh generated by the software TRIANGLE [23] on a cross-section of the RVE of a three-phase composite.

made by various homogenization schemes (the M–T/M– T scheme can be used because the two reinforcing phases only differ by the material). One can observe the good agreement between all models. The results obtained by homogenization schemes are very similar to each other, due to this particular composite (low volume fraction of aligned long fibers).

3-phase composite with 2 kinds of identically shaped aligned fibers. E0=1GPa, ν0=0.35, E1=10GPa, E2=50GPa, ν1=ν2=0.3, Ar = Ar =10 1

Longitudinal macro elastic modulus / E0

8

7

2

Mori-Tanaka, Voigt Mori-Tanaka, Mori-Tanaka

6

Mori-Tanaka, Reuss

5

4

3

2

1 0

0.05

0.1

0.15

0.2

0.25

0.3

0.35

0.4

0.45

0.5

Total volume fraction of inclusions v1+v2 (v1=v2) Fig. 8. Prediction of the longitudinal macro elastic modulus with several homogenization methods on a three-phase composite with the same shape and orientation for the inclusions (but different materials).

O. Pierard et al. / Composites Science and Technology 64 (2004) 1587–1603

1599

Table 4 Normalized macro transverse elastic modulus predicted by finite elements (FE) and various two-step homogenization schemes Method

FE

M–T/V

M–T/M–T

M–T/R

D-I/V

D-I/R

ET =E0

1.2674

1.2342

1.2340

1.2335

1.2407

1.2399

M–T: Mori–Tanaka, D-I: double inclusion, V: Voigt, R: Reuss.

5.3. Thermo-elastic composites

106 K1 ) reinforced with aligned inclusions (E1 ¼ 172 GPa, m1 ¼ 0:2, a1 ¼ 106 K1 ). We simulate a traction test in the direction of the axis of revolution. We compare our M–T predictions with those of van Es [25]. His model is based on a simplified two-phase composite which is just a superposition of two layers in which one is the matrix and the other is the reinforcement. The interaction between the inclusions is thus not directly taken into account. However, in this formulation appears the homogenized stiffness tensor of linear elasticity which might be calculated previously through a classical M–T approach. Results of the predictions of the coefficient of thermal expansion (CTE) in the longitudinal direction are shown in Fig. 11 for various volume fractions of the inclusions and various aspect ratios. It is seen that the predictions of the two methods are almost identical and that for a particular aspect ratio of the inclusions, the macro longitudinal CTE becomes independent of the volume fraction of the reinforcements.

5.3.1. First test Consider a thermo-elastic composite made of a polymer matrix (E0 ¼ 3:0 GPa, m0 ¼ 0:35, a0 ¼ 70 

5.3.2. Second test We now consider a short glass fiber reinforced composite (E1 ¼ 72:5 GPa, m1 ¼ 0:2, a1 ¼ 4:9  106 K1 ,

5.2.6. Sixth test The composite studied here is a two-phase composite with randomly oriented inclusions. The matrix is made of epoxy (E0 ¼ 4:5 GPa, m0 ¼ 0:38) reinforced with mica platelets (E1 ¼ 172 GPa, m1 ¼ 0:2 and Ar1 ¼ 0:04). In order to predict the macroscopic elastic modulus of composites with randomly oriented inclusions, van Es [25] suggested the following relation: E ¼ 0:49Ek þ 0:51E? where Ek (respectively E? ) is the YoungÕs modulus in the direction of the revolution axis of the spheroids (respectively perpendicular to it) obtained with corrected Halpin–Tsai equations. In Fig. 10, van EsÕ predictions are compared to our two-step homogenization results and experimental data. We observe that M– T/Voigt gives a good agreement with the results of van Es (who used a specific method for this kind of composite) and experimental data, at least for small volume fractions of inclusions.

Two-phase composite with randomly oriented platelets. E0=4.5GPa, ν0=0.38, E1=172GPa, ν1=0.2, Ar = 0.04 15

Macro elastic modulus [GPa]

14 Mori-Tanaka, Voigt

13

Mori-Tanaka, Mori-Tanaka 12 Mori-Tanaka, Reuss 11

van Es results

10

experimental

9 8 7 6 5 4 0

5

10

15

20

25

Volume fraction of inclusions [%] Fig. 10. Prediction of the macro elastic modulus of a two-phase composite with randomly oriented platelets. The two-step homogenization method is compared to van EsÕ model and experimental data.

1600

O. Pierard et al. / Composites Science and Technology 64 (2004) 1587–1603

Two phase thermoelastic composite −6

E0=3.0GPa, ν0=0.35, α0=70.10

−6

1

K , E1=172GPa, ν1=0.20, α1=1.10

−1

K

van Es, v =1% 1 van Es, v1=2% van Es, v1=5% van Es, v1=10% M T, v1=1% M T, v1=2% M T, v1=5% M T, v1=10%

−6

−1

CTE along direction of traction [10 K ]

120

100

80

60

40

20

0 −4 10

10

−3

10

−2

10

−1

0

10

1

10

2

10

3

10

4

10

Aspect Ratio Fig. 11. Influence of the aspect ratio on the macro thermal expansion coefficient (comparison of M–T with van Es [25]) for different volume fractions of a thermo-elastic composite.

Ar1 ¼ 35:58 and v1 ¼ 8%). The matrix characteristics are: E0 ¼ 1:57 GPa, m0 ¼ 0:335, a0 ¼ 108:3  106 K1 . A tension test is performed in a direction aligned with the fibers. We compare our predictions of the macro longitudinal and transverse YoungÕs moduli and CTE with experimental data found in Lusti et al. [18]. In that paper, experimental data are compared against several simplified methods which we recall hereafter. McCulloughÕs model predicts the macro characteristics directly. Two micro approaches were proposed by Takao and Taya. The first one is a simple second order average of the longitudinal and transverse CTE (noted hereafter Takao–Taya/aggregate). The second one, proposed by Tandon and Weng tried to include elastic

constraints. This approach is noted Tandon–Weng/T– T/laminate hereafter. Comparison of all these approaches with our M–T and D-I predictions of the macro properties is reported in Table 5. Note that for the predictions made with homogenization schemes, we also performed predictions by multiplying the Ar of the fibers by 1.25 in order to try to model the cylinders by spheroids, as proposed by Li and Ponte Casta neda [16]. Table 5 shows that our predictions, either with M–T or D-I, are very close to experimental results. These predictions also agree with Palmyra [21] FE results while our procedure is much less time-consuming (a small fraction of a second of CPU time on an ordinary PC versus several hours). Note that for the

Table 5 Short glass fiber reinforced composite: predictions of the macroscopic thermomecanical properties by using various homogenization schemes and comparison with several predictive formulae [18] EL (GPa) Experiment [18] M–T M–TðAr  1:25Þ D-I D-IðAr  1:25Þ FE [18] Takao–Taya/aggregate [18] Tandon–Weng/T–T/laminate [18] McCullough [18]

5.99 6.097 6.401 6.136 6.432

<5.89

ET (GPa) 1.929 1.932 1.938 1.941

aL (106 K1 )

aT (106 K1 )

27.7  1.7 30.1 28.9 30.0 28.8 29.3  0.1 32.0 29.4 31.4

121  1 120.6 121.0 120.1 120.5 119  0.1 120 119 121

O. Pierard et al. / Composites Science and Technology 64 (2004) 1587–1603

1601

Table 6 Transversely isotropic graphite fibers-reinforced composites: prediction of the longitudinal and transverse macroscopic CTE with direct predictive formulae, finite-element method [6], homogenization scheme (M–T) and experimental results (Exp.) Composite

v1

CTE (106 K1 )

Exp.

FE (Hex)

FE (Sq)

M–T

SH

RH.

P75/934 P75/934 P75/930 P75/930 C6000/PMR15 C6000/PMR15

0.48 0.48 0.65 0.65 0.63 0.63

aL aT aL aT aL aT

)0.584 19.18 )0.598 17.62 )0.118 12.46

)0.512 19.07 )0.627 14.04 )0.104 12.56

)0.512 19.06 )0.625 14.00 )0.099 12.36

)0.512 18.92 )0.627 13.90 )0.104 12.54

)0.537 19.70 )0.644 14.80 )0.125 14.30

)0.512 18.90 )0.627 13.90 )0.104 12.40

Shapery (SH), Rosen and Hashin (RH), FE on a hexagonal (Hex) or square (Sq) fiber arrangement [6], homogenization scheme (M–T) and experimental results (Exp.).

mean-field CTE prediction, multiplying the Ar by 1.25 leads to slightly better results. 5.3.3. Third test We now consider graphite fiber-reinforced composites. The main difference with the previous cases is that the inclusions material is not isotropic anymore but transversely isotropic (subscript L denotes the longitudinal direction – the one of anisotropy – and T a transverse one). We propose hereafter a slight modification of the method in order to take into account this particularity. Since EshelbyÕs tensor only depends on the matrix properties and the shape of the inclusions, it will not be modified when used in M–T. However, when storing the tensors (accordingly with the traditional storage of strain and stress tensors), we have to write the stiffness tensor and the thermal expansion vector under the following form, using the direction 3 as the direction of the aligned fibers:

0

1=ET B mT =ET B B mLT =EL C1 ¼ B B 0 B @ 0 0

mT =ET 1=ET mLT =EL 0 0 0

mLT =EL mLT =EL 1=EL 0 0 0

0 0 0 ð1 þ mT Þ=ðET Þ 0 0

T

E1;T ¼ 3:35 GPa, l1;LT ¼ 1:30 GPa, m1;LT ¼ 0:20, m1;T ¼ 0:40, Ar1 ¼ 1) and P75 (E1;L ¼ 79:8 GPa, E1;T ¼ 1:38 GPa, l1;LT ¼ 1:00 GPa, m1;LT ¼ 0:20, m1;T ¼ 0:40, Ar1 ¼ 1) where mLT represents the ratio of the strain in the transverse direction over the strain in the longitudinal direction in a uniaxial tension test along the fibers. We can now make predictions by using the M– T homogenization scheme, compare them to several experimental results, some well-known analytical predictions for the CTE and FE results assuming general plane strain on square and hexagonal arrays of fibers [6]. This is illustrated in Table 6. Note that the SchaperyÕs predictive formula, which exists for the longitudinal and transverse directions, is only valid for isotropic constituents and so it does not account for the different values of the CTE of the inclusions in different directions. We can note the good coherence of our predictions, especially for the transverse macro CTE which is much more dependent on the predictive technique used.

0 0 0 0 1=lLT 0

11 0 0 C C 0 C C ; 0 C C 0 A 1=lLT

a1 ¼ ð aT aT a L 0 0 0Þ :

6. Conclusions

For this test, two materials are considered for the matrix: 930/934 epoxy (E0 ¼ 0:63 GPa, m0 ¼ 0:37 and a0 ¼ 24:4  106 K1 ) and PMR15 polyimide (E0 ¼ 0:50 GPa, m0 ¼ 0:35 and a0 ¼ 20:0  106 K1 ), and two kinds of graphite fibers: C6000 (E1;L ¼ 33:8 GPa,

In this paper, we focused on three subjects and presented their theory in detail. Firstly, for two-phase composites, a powerful yet simple predictive tool is LielensÕ interpolative (D-I). Secondly, for multi-phase composites, a general two-step homogenization procedure gives physically acceptable results. Thirdly, a gen-

1602

O. Pierard et al. / Composites Science and Technology 64 (2004) 1587–1603

eral method allows to formulate the thermo-elastic version of any homogenization model defined by its isothermal strain concentration tensors. The predictions have been extensively validated against experimental data or FE results for numerous composite systems. From this study, several observations can be made. In the case of fiber-reinforced two-phase composites, we illustrated the influence of the variation of the angle between the direction of traction and the fibersÕ direction. It was found that sometimes the lowest elastic modulus is obtained for an angle around 60 and not 90 . In the case where reinforcements are platelets (aspect ratio Ar < 1), we also obtained an excellent agreement with FE results, for various values of Ar . For multi-phase composites, we tested several homogenization schemes based on a two-step procedure, where the first one was M–T. Direct M–T homogenization (which is equivalent to a two-step M–T/M–T) was also tested. As explained theoretically by Benveniste et al. [3], when the inclusions have different shapes, a direct M–T does not satisfy the required symmetries and might lead to a non-symmetric overall stiffness tensor. We have also found that a direct M–T might predict a macro stiffness which is outside the M–T/Voigt and M– T/Reuss estimates. However, M–T scheme can still be used for the second step when the inclusions are aligned and similarly shaped. For thermo-elastic composites, we compared our predictions of the CTE with experimental data, either when each phase is isotropic or when the fibers present a transversely isotropic behavior. In both cases, the longitudinal and transverse values of the CTE were in excellent agreement with target results. There are many directions for future work. For instance, for two-phase composites, mean-field homogenization models have been recently extended to the elasto-plastic regime in a general setting (unloading, cyclic loadings, any elasto-plastic model for either phase) by Doghri and Ouaar [10] and Doghri and Friebel [9]. Likewise, work on other generalizations (e.g., rate-dependent behavior, multi-phase nonlinear composites, thermally-coupled nonlinear problems) needs to be performed and validated.

Acknowledgements We would like to thank Mr. A. Ouaar for helpful discussions, Prof. H.J. B€ ohm for providing the Palmyra FE results of the first test in Section 5.1 and Mr. C. de Caritat de Peruzzis for the FE results of the fifth test in Section 5.2. We gratefully acknowledge partial financial support of the IUAP P5/08 project ‘‘From microstructure towards plastic behaviour of single- and multi-

phase materials’’ of the Federal Office for Scientific, Technical and Cultural Affairs, Belgian State Prime MinisterÕs Office.

Appendix A. Second-step homogenization of multi-phase composites with M–T BenvenisteÕs interpretation of Mori–Tanaka for multi-phase composites is that each inclusion can be seen as isolated in a matrix which undergoes a strain field equal to the average strain in the matrix. By denoting x0i;p and x1i;p the matrix and inclusions of the grain xi;p (which belongs to a family i with an orientation p), we can write this assumption as heix0 ¼ heix0 i;p

8i; p:

ðA:1Þ

Using this relation in the expression of the average strain, we get heix ¼ hheixi;p ii;Wi ¼ hv0 heix0 þ ð1  v0 Þheix1 ii;Wi i;p i;p  E e ¼ hðv0 I þ ð1  v0 ÞBi;p : heix0 i;p

¼ hv0 I þ ð1 

i;Wi

v0 ÞB ei;p ii;Wi

: heix0  ¼ ðv0 I þ ð1  v0 ÞhBei;p ii;Wi : heix0 :

ðA:2Þ

Doing the same for the average stress, we get hrix ¼ hhrixi;p ii;Wi ¼ hv0 hrix0i;p þ ð1  v0 Þhrix1 ii;Wi i;p  E e ¼ hðv0 C 0 þ ð1  v0 ÞC i : Bi;p : heix0 i;p

¼ hv0 C 0 þ ð1  v0 ÞC i :

i;Wi

Bei;p ii;Wi

: heix0  ¼ ðv0 C 0 þ ð1  v0 ÞhC i : Bei;p ii;Wi : heix0 :

ðA:3Þ

By putting expression (A.2) in (A.3), we get the expression of the macro stiffness tensor:  C ¼ ðv0 C 0 þ ð1  v0 ÞhC i : Bei;p ii;wi 1 : ðv0 I þ ð1  v0 ÞhB ei;p ii;wi :

References [1] ABAQUS. A general-purpose finite element program, Hibbitt. Karlsson & Sorensen, Inc., Pawtucket, RI, USA, 2001. [2] Benveniste Y. A new approach to the application of Mori– TanakaÕs theory in composite materials. Mech Mater 1987;6: 147–57. [3] Benveniste Y, Dvorak GJ, Chen T. On diagonal and elastic symmetry of the approximate effective stiffness tensor of heterogeneous media. J Mech Phys Solids 1991;39(7):927–46.

O. Pierard et al. / Composites Science and Technology 64 (2004) 1587–1603 [4] B€ ohm HJ. Private communication, May 2003. [5] B€ ohm HJ, Han W. Comparisons between three-dimensional and two-dimensional multi-particle unit cell models for particle reinforced metal matrix composites. Modeling Simul Mater Sci Eng 2001;2:47–65. [6] Bowles DE, Tompkins SS. Prediction of coefficients of thermal expansion for unidirectional composites. J Compos Mater 1989;23:370–88. [7] Camacho CW, Tucker III CL, Yalvac S, McGee RL. Stiffness and thermal expansion predictions for hybrid short fiber composites. Polym Compos 1990;11(4):229–39. [8] de Caritat de Peruzzis C. Mod
elisation microm
ecanique de mat
eriaux composites
elastoplastiques biphas
es. Graduation report, Universit
e Catholique de Louvain, Belgium, 2003. [9] Doghri I, Friebel C. Effective elasto-plastic properties of inclusion-reinforced composites. Study of shape, orientation and cyclic response. Mechanics of materials, in press, 2004. [10] Doghri I, Ouaar A. Homogenization of two-phase elasto-plastic composite materials and structures: Study of tangent operators, cyclic plasticity and numerical algorithms. Int J Solids Struct 2003;40(7):1681–712. [11] Eshelby JD. The determination of the elastic field of an ellipsoidal inclusion, and related problems. Proc Roy Soc London, Ser A 1957;241:376–96. [12] Friebel C. Mod
elisation et Simulation micro-macro de mat
eriaux composites. Etude de lÕinfluence des renforts (mat
eriau, forme et orientation). Graduation report, Universit
e Catholique de Louvain, Belgium, 2002. [13] Guo-Zheng Kang, Gao Qing. Tensile properties of randomly oriented short d  Al2 O3 fiber reinforced aluminum alloy composites: II. Finite element analysis for stress transfer, elastic modulus and stress–strain curve. Composite Part A 2002;33:657– 67. [14] Hashin Z, Shtrikman S. A variational approach to the theory of the elastic behaviour of multiphase materials. J Mech Phys Solids 1963;11:127–40.

1603

[15] Li JY. On micro-mechanics approximation for the effective thermoelastic moduli of multi-phase composite materials. Mech Mater 1999;31:149–59. [16] Li G, Ponte Casta~ neda P. Variational estimates for the elastoplastic response of particle-reinforced metal–matrix composites. In: Ostoja-Starzewki M, Jasiuk I, editors. Micromechanics of random media, Appl. Mech. Rev., vol. 47, No. 1, Part 2, S77-S94. [17] Lielens G. Micro-Macro Modeling of Structured Materials. PhD thesis, Universit
e Catholique de Louvain, Belgium, 1999. [18] Lusti HR, Hine PJ, Gusev AA. Direct numerical predictions for the elastic and thermoelastic properties of short fiber composites. Compos Sci Technol 2002;62:1927–34. [19] Mori T, Tanaka K. Average stress in matrix and average elastic energy of materials with misfitting inclusions. Acta Mater 1973;21:571–4. [20] Nemat-Nasser S, Hori M. Micromechanics: overall properties of heterogeneous materials. second ed. Amsterdam: Elsevier; 1999. [21] Palmyra A. A finite-element software for heterogeneous microstructures. Zurich, Switzerland: MatSim GmbH; 2001. [22] Taya M, Chou T. On two kinds of ellipsoidal inhomogeneities in an infinite elastic body: an application to a hybrid composite. Int J Solids Struct 1981;17:553–63. [23] TRIANGLE. A two-dimensional quality mesh generator and Delaunay triangulator. J.R. Shewchuk, School of Computer Science, Carnegie Mellon University, Pittsburgh, Pennsylvania, USA, 1996. [24] Tucker III CL, Liang E. Stiffness predictions for unidirectional short-fiber composites: review and evaluation. Compos Sci Technol 1999;59:655–71. [25] van Es M. Polymer–clay Nanocomposites. The importance of particle dimensions. PhD thesis, Technische Universiteit Delft, The Netherlands, 2001. [26] Weissenbek E, B€ ohm HJ, Rammerstorfer FG. Micromechanical investigations of arrangement effects in particle reinforced metal–matrix composites. Computat Mater Sci 1994;3: 263–78.