Micromechanics of inelastic composites with misaligned inclusions: Numerical treatment of orientation

Comput. Methods Appl. Mech. Engrg. 195 (2006) 1387–1406 www.elsevier.com/locate/cma

Micromechanics of inelastic composites with misaligned inclusions: Numerical treatment of orientation Issam Doghri a

a,*

, Laurent Tinel

b

Universite´ Catholique de Louvain (UCL), CESAME, Baˆtiment Euler, 4 Avenue G. Lemaıˆtre, B-1348 Louvain-la-Neuve, Belgium b e-Xstream engineering SA, 4 Avenue A. Einstein, B-1348 Louvain-la-Neuve, Belgium Received 30 November 2004; received in revised form 20 March 2005; accepted 9 May 2005

Paper dedicated to Professor Thomas J.R. Hughes on the occasion of his 60th birthday

Abstract This paper is concerned with the micromechanics of multi-phase fiber- or inclusion-reinforced inelastic composites with a special emphasis on the numerical treatment of misaligned orientation. Mean-field homogenization of these composites involves volume and orientation averaging. The latter usually needs average orientation measures known as the second- and fourth-rank orientation tensors a and A, respectively. Usually the only orientation data available is a and a first issue is to reconstruct A using closure approximations. We propose a new weighted regularization method which deals with orientations which are strictly neither 1D, 2D nor 3D and ensures a smooth transition between the 3 cases. A second issue is that there are situations where the inclusionsÕ orientation distribution function (ODF) is needed for homogenization but is unavailable. The ODF is recovered using a and A and numerically computed at a number of discrete orientations. A third issue is the actual computation of orientation averaging integrals. An efficient algorithm approximates the integrals with sums over a finite number of orientations. Orientation and volume averaging concepts are applied to the mean-field homogenization of two classes of multi-phase composites: linear thermo-elastic and rateindependent inelastic. In the latter case, an incremental formulation is proposed which enables the simulation of unloading, cyclic and otherwise non-proportional loadings. Implicit time-discretization and consistent tangent operators are employed. The procedures and the corresponding algorithms were implemented in the [DIGIMAT version 1.4, 2004. Linear and nonlinear multi-scale material modeling software, e-Xstream engineering SA (http://www.e-Xstream.com), Louvain-la-Neuve, Belgium.] software and several numerical simulations show their accuracy and efficiency.
2005 Elsevier B.V. All rights reserved. Keywords: Composites; Micromechanics; Homogenization; Orientation distribution function; Orientation tensors; Inelasticity

*

Corresponding author. Tel.: +32 10 478042/472350; fax: +32 10 472180. E-mail addresses: [email protected] (I. Doghri), [email protected] (L. Tinel).

0045-7825/$ - see front matter
2005 Elsevier B.V. All rights reserved. doi:10.1016/j.cma.2005.05.041

1388

I. Doghri, L. Tinel / Comput. Methods Appl. Mech. Engrg. 195 (2006) 1387–1406

1. Introduction Fiber- or inclusion-reinforced composite materials enjoy a wide use in various industries. The modern way to model them is via multi-scale or micro–macro methods which predict the influence of the microstructure on the macroscopic properties. Scale-transition methods for these composites can be classified into four categories. (1) Direct finite element (FE) simulations of a unit cell (assuming periodic microstructures) or a representative volume element (RVE). A pioneering contribution was made by Adams [3] who realized transverse load simulations of an elasto-plastic matrix reinforced with long and aligned elastic fibers. A large number of other papers followed, e.g., [31,58,8,51,32,49,11,30]. (2) Transformation field analysis or subcell method (where a unit cell or a RVE is subdivided into a number of 2D pixels or 3D voxels); see [23,2,6,48,27]. The latter authors subdivided each cell into 4 subcells and enforced traction and strain compatibility in terms of average stresses and strains in each subcell. (3) Homogenization based on asymptotic expansion of the displacement field (it assumes a periodic microstructure and ends up with a unit cell problem to be solved by FE); examples are [35,31]. (4) Mean-field homogenization based on assumed relations between average values of microstrain and stress fields in each phase; examples are [55,36,47,42,26,45,19,18,43]. In this paper, we study the latter scale transition method. There are several phenomena it is unable to describe (clustering, percolation, strain localization and size effects—except when combined with a Hall– Petch-type model) and for homogenization models based on the Eshelby [24] results, the inclusionsÕ shape can only be ellipsoidal. However, when the conditions of their application are met, and when one is interested only in the effective properties and the per-phase averages of stress and strain, mean-field homogenization models provide by far the most cost-effective solution (in terms of computer time and also ease-of-use). The models have been applied with success to a large variety of composites based on polymer, metal, ceramic or concrete matrices. Most of the literature on micro–macro modeling of composites (including the above-mentioned references) deals with fixed orientation fibers or orientation-irrelevant spherical inclusions, while this paper deals with multi-phase composites. These are such that a matrix material is reinforced with multiple phases of misaligned inclusions. Each inclusion is supposed to be a spheroid (that is an ellipsoid of revolution) with unit axis vector p and aspect ratio Ar > 0 (boldface symbols denote tensors or matrices, the rank of which is indicated by the context). Fibers and platelets can be approximated with prolate (Ar > 1) and oblate (Ar < 1) spheroids, respectively. Mean-field homogenization of such multi-phase composites has been studied in linear thermo-elasticity (e.g., [10,22,44]) and elasto-plasticity (e.g., [21,20]). However, some important issues related to the numerical treatment of misaligned inclusionsÕ orientation have not been addressed and it is the objective of the present paper to study them. A first issue is that the inclusionsÕ orientation distribution function (ODF) w(p) is usually unavailable (e.g., injection molding of short fiber reinforced polymer composites). The only information one has is a second-rank orientation tensor a which is the ODF-weighted average of p
p (the symbol
designates a tensor or dyadic product). However, in the computation of effective stiffness operators, one needs a fourth-rank orientation tensor A which is the ODF-weighted average of p
p
p
p. An important problem is to approximate A from a. This is the subject of closure approximations. We propose a new weighted formula for cases (most often encountered in practice) where the orientation is strictly neither 2D nor 3D. Numerical simulations show that unless the new weighted scheme is used, unrealistic predictions can be obtained in some instances.

I. Doghri, L. Tinel / Comput. Methods Appl. Mech. Engrg. 195 (2006) 1387–1406

1389

A second issue is that in many cases (e.g., inelastic composites) the ODF although unavailable is needed in order to homogenize the composite. The ODF must be recovered from orientation tensors a and A. We present numerical algorithms for the evaluation of the reconstructed ODF at discrete orientations. Analytical fitting is proposed either for computational efficiency or for cases where the recovered ODF has physically unacceptable results (e.g., negative) for some orientations. A third issue is that in order to compute the compositeÕs homogenized response, one needs to compute orientation averages of appropriate microfields. Based on the ODF values at discrete orientations, we present efficient algorithms for the approximation of orientation integrals by finite sums. The algorithms must be accurate yet computationally cheap, because the cost increases rapidly (especially in 3D) with the number of discrete orientations. The mathematical and numerical tools for orientation treatment are applied to the homogenization of multi-phase composites within a general two-step procedure in two cases: linear thermo-elasticity and rate-independent inelasticity. The latter is based on the incremental formulation and algorithms proposed by Doghri and Ouaar [19], Doghri and Friebel [18] and Doghri and Tinel [20] which enable the simulation of unloadings, cyclic and otherwise non-proportional loadings, unlike a secant or total deformation formulation (e.g., [55]). The implementation of the incremental formulation is based on implicit time-discretization and algorithmic or consistent tangent operators [29,52,53]. The paper has the following outline. In Section 2, basic averaging results for multi-phase composites are presented, and a two-step homogenization procedure is described and then applied to linear thermo-elasticity and rate-independent inelasticity. Section 3 shows how A is reconstructed from a (exact linear and quadratic expressions, closure approximations and a new weighting formula). In Section 4 both the theory of ODF reconstruction from a and A and its implementation are presented (including an analytical fitting of the ODF when necessary). Section 5 deals with the numerical computation of orientation averaging (based on an algorithm for iso-size facets in 3D and the computation of appropriate weights). Several numerical simulations are given in Section 6. They illustrate the topics covered in this work and the accuracy and efficiency of the proposed numerical methods. The paper ends with conclusions and a few directions for future work in Section 7.

2. Homogenization of multi-phase composites 2.1. General results The main objective of mean-field homogenization is to relate the volume averages of stress hriX and strain hiX over a RVE occupying a domain X (see for instance [40]). Consider a RVE containing a matrix (domain X0, volume fraction v0) and a large number of randomly located spheroids, each characterized by its unit axis vector p, aspect ratio Ar and material properties. With respect to (w.r.t.) a fixed Cartesian frame (O, x, y, z), p is defined by two spherical angles: h 2 [0, p] between p and (O, ez) and / 2 [0, 2p] between the projection of p on the (O, x, y) plane and (O, ex). The inclusions are classified into N phases (i) of volume fractions vi, v0 þ

N X

vi ¼ 1.

ð1Þ

i¼1

The matrix and each inclusion of each phase are supposed to be homogeneous. Each inclusionsÕ phase is characterized by a given material, a constant aspect ratio Ari and a differentiable ODF wi(p) defined such that the probability of finding a fiber belonging to phase (i) and whose orientation is between p and (p + dp) is w(p) dp. The ODF satisfies the following conditions:

1390

I. Doghri, L. Tinel / Comput. Methods Appl. Mech. Engrg. 195 (2006) 1387–1406

wi ðpÞ ¼ wi ðpÞ;

I

wi ðpÞ dp ¼ 1;

ð2Þ

where the first equality expresses the fact that inclusions oriented at p and (p) are indistinguishable, and the second is a normalization condition. The RVE is decomposed into a set of infinitesimal pseudo-grains. Each pseudo-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 all inclusions (in concentration (1  v0)) of phase i of orientation between p and (p + dp). Let l(x) represent an arbitrary microfield (e.g., microstrain or stress fields). It can be shown (e.g., [44]) that the volume average of l(x) over the RVE can be written as follows: hliX ¼

N X i¼1

vi hhlixi;p iwi  hhlixi;p ii;wi ; ð1  v0 Þ

ð3Þ

where hlixi;p represents the volume average over a pseudo-grain xi,p and hiwi a wi-weighted average. The orientation average of a function l(p) is its ODF-weighted average defined as follows: I ð4Þ hlðpÞiwi  lðpÞwi ðpÞ dp. In the compact format (3), hii;wi is the orientation average over all pseudo-grains, that is over all phases i and all orientations. Result (3) shows that homogenization of the RVE can be performed in two steps: (1) homogenization of each pseudo-grain individually, and (2) homogenization of the pseudo-grains. The twostep homogenization procedure is illustrated in Fig. 1. A direct (one-step) homogenization of the composite (top of Fig. 1) is not equivalent in general to a two-step homogenization, but we prefer the latter, as discussed in Section 2.2. Note that result (3) and therefore the two-step approach are very general and do not depend on the microfield l(x) or the material model for any phase. 2.2. Linear thermo-elastic composites For simplicity, a pseudo-grainÕs domain xi,p is designated by x, with subdomains x0 for the matrix (of stiffness c0) and x1 for the inclusionsÕ phase (c1). In isothermal linear elasticity, any mean-field homogenization model is identified by its strain concentration tensors B
and A
defined as follows:

DECOMPOSITION

FIRST STEP

SECOND STEP

Fig. 1. Two-step homogenization. First step: the RVE is decomposed into a set of pseudo-grains and each one is homogenized. Second step: homogenization over all pseudo-grains.

I. Doghri, L. Tinel / Comput. Methods Appl. Mech. Engrg. 195 (2006) 1387–1406

hix1 ¼ B
: hix0 ¼ A
: hix ; 1

A
¼ B
: ½ð1  v0 ÞB
þ v0 I s  ;

1391

ð5Þ

where a colon designates a tensor product contracted over two indices and Is is the fourth-rank symmetric identity tensor. The knowledge of B
and A
allows to compute the effective stiffness and coefficients of thermal expansion (CTE) in thermo-elasticity; see [10,7,44]. In general, B
depends on c0, c1, v0 and (except for the simplest homogenization models such as Voigt-uniform strain or Reuss-uniform stress) on the shape and orientation of inclusions (I) through EshelbyÕs tensor EðI; c0 Þ [24]. A model which works particularly well for moderate values of v1 is the Mori and Tanaka [39] model (M–T) for which B
is equal to an isolated inclusionÕs strain-concentration tensor H
, n h io1 . ð6Þ B
ðI; c0 ; c1 ; v0 Þ ¼ H
ðI; c0 ; c1 Þ ¼ I s þ EðI; c0 Þ : ðc0 Þ1 : c1  I s More models are considered in [44]. For 3D or 2D randomly oriented fibers in an isotropic matrix, Christensen and Waals [14] computed the effective stiffness based on a two-cylinder model. Tandon and Weng [54] extended the M–T model to two-phase composites with misaligned inclusions. Qiu and Weng [46] found that this direct M–T homogenization does not give the same predictions as a two-step M–T/self-consistent method. In this work, we adopt a two-step homogenization procedure: M–T for each two-phase pseudo-grain and Voigt for the ‘‘aggregate’’; see [10,44]. This M–T/Voigt procedure ensures physically acceptable results and a symmetric effective stiffness tensor. Indeed, Benveniste et al. [7] showed that M–T model gives a diagonally-symmetric effective stiffness in 2 cases. First, when the composite contains exactly 2 materials (one for the matrix and another for the inclusions) and the inclusions have the same stiffness matrix in a fixed Cartesian frame (which for misaligned inclusions restricts their material to being isotropic). Second, for a composite containing 3 or more materials, when all inclusions have the same shape, aspect ratio and orientation. For isotropic materials, each pseudo-grain xi,p is transversely isotropic and a closed-form expression of its effective stiffness cxi;p exists which depends on the invariants and the direction of anisotropy (which for each pseudo-grain is a unit vector p along the spheroidsÕ revolution axis). It can then be shown (e.g., [10]) that for each inclusionsÕ phase ðiÞ, the orientation averages hcxi;p iwi only depend on the above-mentioned invariants and the orientation tensors of the phase (Section 3.1). So in this case, one does not need a precise knowledge of neither each pseudo-grainÕs orientation p nor each phaseÕs ODF wi(p). Similar results hold for the CTE. 2.3. Rate-independent inelastic composites For simplicity, x designates a pseudo-grainÕs domain, and subscripts 0 and 1 refer to the matrix and inclusionsÕ phases, resp. If a phase (r) is inelastic but rate-independent, one can write in the isothermal case _ _ tÞ; rðx; tÞ ¼ cr ðx; tÞ : ðx;

8x 2 xr ;

ð7Þ

where cr(x, t) designates a tangent operator. So in principle, we could extend the linear elastic results by writing them in a rate form; we thus obtain the so-called incremental formulation. However, there are at least two major differences with the linear case. First, tangent moduli cr(x, t) are not uniform per phase. This makes impossible the generalization of the solutions of Eshelby and the isolated inclusion problem, which are the foundation of useful mean-field homogenization models. This problem is circumvented by defining so-called reference tangent operators bc r ðtÞ which are uniform per phase _ tÞ; _ rðx; tÞ ¼ bc r ðtÞ : ðx;

8x 2 xr .

ð8Þ

In this work, bc r ðtÞ for each phase is computed by sending per-phase averages of strain and strain increments to the constitutive box of the real material. A second difficulty in the nonlinear regime is that strain concentration tensors depend on reference tangent operators which in turn depend on per-phase averages of

1392

I. Doghri, L. Tinel / Comput. Methods Appl. Mech. Engrg. 195 (2006) 1387–1406

strain. So, in order to compute the latter, we have to solve a nonlinear problem. Rate relations such as (5a) are discretized over a time interval [tn, tn+1] as follows: hDix1 ¼ B
ðtnþa Þ : hDix0 ;

tnþa ¼ ð1  aÞtn þ atnþ1 ; a 20; 1;

ð9Þ

where D designates an increment. For the M–T model, the strain concentration tensor is given by the extension of Eq. (6) as follows: B
ðtnþa Þ ¼ H
ðI; bc 0 ðtnþa Þ; bc 1 ðtnþa ÞÞ.

ð10Þ

The data for the pseudo-grain homogenization are: h(tn)ix, hDix and per-phase history variables at tn. The problem is to compute the homogenized stress increment hDrix and tangent moduli cx(tn+a). Algorithms are presented in [19,18,20]. Notice that B
does not remain constant w.r.t. time, but changes with time steps, and within each time step with each iteration of the homogenization algorithm. For two-phase composites made of an elasto-plastic matrix reinforced with elastic spheres or long and aligned fibers, an algorithm was proposed by [34] based on backward Euler time integration (a = 1) and a ep numerical computation of EshelbyÕs tensor using the ‘‘continuum’’ tangent operator bc 0 (relating infinitesimal stress and strain increments). In the present work, we use a mid-point rule (a = 1/2) and algorithmic alg moduli bc r ¼ bc r obtained by consistent linearization of the time-discretized constitutive equations, as exalg plained in [29,52,53,17]. Those operators bc r are anisotropic. In this work, we use them in all computations alg iso except that of EshelbyÕs tensor which is approximated with an isotropic part of bc 0 , i.e., EðI; bc 0 Þ. Otherwise, the predicted macroresponse is usually too stiff, as found out by [25,26,37,9,45]. There are some partial but still no definitive explanations for the problem or its workaround. According to Molinari et al. [38], the incremental formulation leads to the cumulation over successive time steps of the error made on the ep alg iso assumption of uniform reference tangent operators. Doghri and Ouaar [19] prove that bc 0 P bc 0 P bc 0 (where P means: ‘‘stiffer than’’) and remark that the overall response computed with those 3 operators follows the same order. Doghri and Friebel [18] notice that even the often-used secant or total-deformation theory computes EshelbyÕs tensor with an isotropic secant operator. Chaboche and Kanoute´ [12] and Chaboche et al. [13] give plastic flow-based arguments to explain why an isotropic projection being more compliant delivers better predictions than an anisotropic tangent operator. For a multi-phase RVE containing possibly misaligned inclusions, the final macroresponse is obtained via homogenization of all pseudo-grains xi,p. The algorithm is detailed in [20]. Unlike the linear case, even if every material is isotropic, the pseudo-grainsÕ tangent operators are not transversely isotropic. Therefore, the knowledge of the orientation tensors of each inclusionsÕ phase (i) is not sufficient for the determination of the RVEÕs overall properties. One needs to reconstruct the ODF (see Sections 3 and 4) and discretize the orientation integrals into sums over a finite number of orientations (see Section 5).

3. Reconstruction of the fourth-rank orientation tensor 3.1. Motivation Usually, the ODF is unavailable and the only given orientation data is the so-called second-rank orientation tensor a. A first issue is to find an estimate for a fourth-rank orientation tensor A. Those tensors are defined as the following ODF-weighted averages [4,5]: a  hp
piw ;

A  hp
p
p
piw .

ð11Þ

From the definition of a and kpk = 1, the following properties are obvious: aij ¼ aji ;

a11 ; a22 ; a33 P 0;

aii ¼ 1.

ð12Þ

I. Doghri, L. Tinel / Comput. Methods Appl. Mech. Engrg. 195 (2006) 1387–1406

1393

EinsteinÕs summation convention over repeated indices is used throughout unless otherwise indicated. As for A, it is symmetric w.r.t. any permutation of indices and satisfies Aijij P 0 ðno sumsÞ;

Aijll ¼ aij .

ð13Þ

In order to reconstruct A from a, exact formulae exist for aligned or random (2D or 3D) orientations (Section 3.2). For other cases, only approximate so-called closure formulae exist (Section 3.3). We present in Section 3.4 a new weighted formula for cases which are neither strictly 2D nor 3D. In general, orientation averaging is not unique. Indeed, the reconstruction of A from a is not unique but based on various approximate models which might lead to quite different results. If the ODF needs to be recovered (e.g., inelastic composites) than there are approximations also. Finally, when discrete orientation algorithms are used, additional inaccuracies are introduced. Some appreciation of the overall precision of orientation and volume averaging can be gained by examining the examples of Section 6. 3.2. Linear and quadratic closures For random orientation, the so-called linear closure Al(a) is exact. It is given in index notation in [4,5]. In this paper, we re-write the expression of Al(a) in a tensorial format which is both very compact and readily translated into matrix operations Al ðaÞ ¼ al 1
1 þ 2ðal  bl ÞI s þ bl ð1
a þ a
1Þ þ 2bl ½Iða þ 1Þ  IðaÞ;

ð14Þ

where 1 and Is designate the symmetric second- and fourth-rank identity tensors, resp. and I(a)ijkl = 1/2(aikajl + ailajk). Note that Is = I(1). The coefficients in Eq. (14) are given by al ¼ 

1 35

and

bl ¼

1 in 3D; 7

al ¼ 

1 24

and

bl ¼

1 in 2D. 6

ð15Þ

Now, for aligned (fixed) orientation, the following quadratic expression is exact: Aq ðaÞ ¼ a
a.

ð16Þ

3.3. Closure approximations If the fiber orientation is neither random nor fixed, then estimating A from a is no longer exact, and one must use a closure approximation such as the hybrid [5], the orthotropic [15] and the natural [56]. In this work, we implemented the hybrid closure which is defined as the following weighted average between the linear and quadratic expressions: Ah ðaÞ ¼ ð1  f ÞAl ðaÞ þ f Aq ðaÞ;

f ¼ 1  N h det a;

ð17Þ

where the number Nh is chosen as follows: N h ¼ 27 for 3D orientation;

N h ¼ 4 for 2D orientation.

The hybrid closure reduces to the exact linear expression if the orientation is random in 3D or 2D

 1 1 1 1 1 f ¼ 0 if a ¼ DIAG ; ; or a ¼ DIAG ; ; 0 . 3 3 3 2 2 It also reduces to the exact quadratic expression for fixed orientation f ¼ 1 if a ¼ DIAGð1; 0; 0Þ.

ð18Þ

1394

I. Doghri, L. Tinel / Comput. Methods Appl. Mech. Engrg. 195 (2006) 1387–1406

3.4. A new weighting formula One important problem with closure approximations is that they have expressions for either 2D or 3D orientations, while in practice there are numerous cases where one cannot decide for sure whether the orientation is 1D, 2D or 3D. For instance, consider the following cases: a ¼ DIAGð0:8; 0:1; 0:1Þ;

a ¼ DIAGð0:98; 0:01; 0:01Þ.

ð19Þ

The first case is typical of injection molding where many fibers are aligned in the flow direction, while the second situation represents a small perturbation about the fully aligned 1D case. For simulation purposes, one can very well assume that both examples (19) are basically 1D or 2D in the (x, y) or (z, x) planes or 3D. However, it is shown in Table 1 (Section 6.1) that depending on the properties of each phaseÕs material and the contrast between them, computing different approximations of A assuming either 2D or 3D can lead to completely different predictions of the effective CTE and those predictions are either physically wrong or too far away from the aligned 1D case results. This is why we propose a new scheme which computes the relative weights of 1D, 2D and 3D orientations and ensures a smooth transition between these 3 cases. The new approach is based on the following two steps. First we extract 1D, 2D, and 3D projections of a, and we compute the corresponding three closure approximations A1D, A2D and A3D. Second, those three expressions contribute to the final approximation of A via a weighting formula. We now explain those steps in detail. Since a is symmetric, we find its eigenvalues (arranged in decreasing order) and eigenvectors (forming an orthonormal basis). In that basis, a is now represented by a diagonal matrix, adiag ¼ DIAGða1 ; a2 ; a3 Þ;

a1 P a2 P a3 P 0;

a1 þ a2 þ a3 ¼ 1.

ð20Þ

Table 1 Predictions of effective CTE (106/K) and stiffness (first three diagonal terms) (GPa) for a highly contrasted composite

Ref.: a = (1, 0, 0) a = (0.8, 0.1, 0.1) hybrid 3D a = (0.8, 0.1, 0.1) new weighting a = (0.98, 0.01, 0.01) hybrid 2D a = (0.98, 0.01, 0.01) quadratic 2D a = (0.98, 0.01, 0.01) quadratic 3D a = (0.98, 0.01, 0.01) new weighting Ref.: a ¼ ð13 ; 13 ; 13Þ a = (0.343, 0.323, 0.333) hybrid 3D a = (0.343, 0.323, 0.333) new weighting Ref.: a ¼ ð12 ; 12 ; 0Þ a = (0.49, 0.51, 0) hybrid 2D a = (0.49, 0.51, 0) new weighting

a11

a22

a33

C 11

C 22

9.453 8.378 7.287 8.195 13.31% 8.737 7.57% 8.102 14.29% 9.339 1.21%

91.127 96.293 92.084 94.410 3.60% 90.645 0.53% 91.421 0.32% 91.263 0.12%

91.127 96.293 92.084 90.018 1.22% 92.036 1.00% 91.421 0.32% 91.263 0.15%

7.454 5.143 7.154 7.299 2.08% 7.327 1.70% 7.206 3.33% 7.448 0.08%

1.157 1.212 1.169 1.129 2.42% 1.158 0.09% 1.158 0.09% 1.157 0.00%

1.157 1.212 1.169 1.157 0.00% 1.158 0.09% 1.158 0.09% 1.157 0.00%

39.610 37.887 4.35% 37.636 4.98%

39.610 41.398 4.51% 41.612 5.05%

39.610 39.593 0.04% 39.644 0.09%

2.425 2.478 2.20% 2.488 2.60%

2.425 2.370 2.27% 2.366 2.43%

2.425 2.424 0.04% 2.424 0.04%

19.099 19.582 2.53% 19.792 3.63%

19.099 18.631 2.45% 18.427 3.52%

130.20 130.19 0.01% 130.19 0.01%

3.526 3.463 1.79% 3.495 0.88%

3.526 3.589 1.79% 3.618 2.61%

1.157 1.157 0.00% 1.157 0.00%

C 33

Study of three reference cases (aligned, 3D and 2D random) and perturbations around them with different closure approximations and our new weighting scheme. Predictions and relative differences between perturbed and reference results are shown.

I. Doghri, L. Tinel / Comput. Methods Appl. Mech. Engrg. 195 (2006) 1387–1406

1395

The absence of shear terms makes it easy to extract 1D, 2D and 3D projections. The 1D part is simply a1D ¼ DIAGð1; 0; 0Þ.

ð21Þ

The corresponding A1D is the exact quadratic closure (Section 3.2). The 2D projection of a is defined as follows:
 a1 a2 a2D ¼ DIAG ; ;0 . ð22Þ a1 þ a2 a1 þ a2 The corresponding closure approximation A2D is then computed unambiguously using the 2D expressions of Section 3.3. Similarly, the 3D closure expression A3D is found from Section 3.3 by using a3D = adiag. The final approximation to A is given by the proposed weighting formula A ¼ a1 A1D þ a2 A2D þ a3 A3D ;

a1 ¼

a1  a2 a2  a3 a3 ; a2 ¼ ; a3 ¼ . a1 a1 a1

ð23Þ

It is seen that in 1D (a2 = a3 = 0), the formula is correct (a1 = 1, a2 = a3 = 0). In 2D (a3 = 0), a3 = 0, and the closure a2 is to a1, the larger a2 (and if a2 = a1 = 1/2—random 2D—then a1 = 0 and a2 = 1). If a3 5 0, then in general the 1D, 2D and 3D expressions all contribute, but with weights depending on how close we are to each case (and if a1 = a2 = a3 = 1/3—random 3D—then a1 = a2 = 0 and a3 = 1). Finally, A is rotated back from the basis of the eigenvalues of a to the original basis of a.

4. Reconstruction of the ODF Usually, only orientation tensors a and A (typically reconstructed from a as explained in Section 3) are available and the ODF is not. However, the latter is needed in some cases, such as for the mean-field homogenization of inelastic composites. In Section 4.1, the Onat and Leckie [41] formulae for reconstructing the ODF from a and A are presented, while Section 4.2 deals with numerical implementation. 4.1. Theory Onat and Leckie [41] define tensor basis functions of p as follows (see also Advani and Tucker [4]). A second-rank tensor function f(p) is the deviatoric part of (p
p), i.e., f ðpÞ ¼ devðp
pÞ.

ð24Þ

A fourth-rank tensor function F(p) is defined as follows (we give it in a very concise tensor format which can be immediately translated into matrix implementation): FðpÞ ¼ Aq ðp
pÞ  Al ðp
pÞ; q

ð25Þ

l

where A (p
p) and A (p
p) are given by expressions (14) and (16) when one replaces argument a with (p
p). Now define orientation averages b and B as follows: hf ðpÞiw ¼ devðaÞ  bðaÞ; l

hFðpÞiw ¼ AðaÞ  A ðaÞ  BðaÞ;

ð26Þ ð27Þ

where A(a) is the fourth-order orientation tensor—Eq. (11b)—and Al(a) its linear closure approximation— Eq. (14). Onat and Leckie [41] show that the ODF can be reconstructed as follows: wðpÞ  w1 þ w2 bðaÞ : f ðpÞ þ w3 BðaÞ :: FðpÞ;

ð28Þ

1396

I. Doghri, L. Tinel / Comput. Methods Appl. Mech. Engrg. 195 (2006) 1387–1406

where a double colon indicates a tensor product contracted over 4 indices. In 3D, the coefficients are given by w1 ¼

1 ; 4p

w2 ¼

15 ; 8p

w3 ¼

315 . 32p

ð29Þ

While for 2D orientation, w1 ¼

1 ; 2p

2 w2 ¼ ; p

8 w3 ¼ . p

ð30Þ

For random orientation in 3D (a = DIAG(1/3, 1/3, 1/3)) or 2D (a = DIAG(1/2, 1/2,0)), we have dev (a) = 0 and A = Al, thus b(a) = 0 and B(a) = 0. Therefore, the exact ODF expressions (w = w1) are recovered in random 3D or 2D. 4.2. Numerical implementation In practice, we proceed as follows. Orientations tensors a and A (usually reconstructed from a—Section 3) are known data. Thus, tensors b(a) and B(a)—Eqs. (26) and (27)—are known. Tensors f(p) and F(p)— Eqs. (24) and (25)—and therefore the ODF-formula (28)—are computed at discrete values of p. In 2D, these correspond to discrete angle values /i 2 [0, 2p]. For computational efficiency (or for cases when some reconstructed discrete ODF values are negative), an analytical fit to the discrete values w(/i) can be sought. For instance, we often use the following Gaussian: "
2 # 1 1 /  /M wG ð/Þ ¼ pffiffiffiffiffiffi exp  ; ð31Þ 2 S S 2p where /M is the mean value and S the standard deviation. In 3D, the algorithm of Section 5.1 computes a finite number of iso-size facets, each with a unit normal p defined by spherical angles (hi, /i). Tensors f(p) and F(p) and the ODF w(p) are computed for each discrete orientation (hi, /i) 2 [0, p] · [0, 2p]. As in 2D, whenever possible, an analytical fit is sought, e.g., with a Gaussian wðh; /Þ  wG ðhÞwG ð/Þ.

ð32Þ

5. Discrete orientation averaging In this section, we show how orientation averaging integrals such as (4) which are needed in homogenization are discretized into sums over a finite number of discrete orientations. 5.1. Algorithm for iso-size facets in 3D We implemented an algorithm which was proposed by Weber et al. [57] for critical plane criteria in high cycle fatigue. Accordingly, a unit vector p defined by two spherical angles h (meridien) and / is viewed as the outer unit normal to a unit sphere. Equal-sized increments of h are considered, and the corresponding values of /-increments are computed such that the surface of the unit sphere is subdivided into facets of almost equal areas (see Figs. 3 and 4 in [57]). In practice, it suffices to consider half the unit sphere, (h, /) 2 [0, p] · [0, p]. The input is Dh, the constant increment of h. A spherical ring situated between angles (h  Dh/2) and (h + Dh/2) has the following area:
 Z hþDh2 Dh Sh ¼ p sinðcÞ dc ¼ 2p sinðhÞ sin . ð33Þ Dh 2 h 2

I. Doghri, L. Tinel / Comput. Methods Appl. Mech. Engrg. 195 (2006) 1387–1406

1397

On the ring h = p/2 ± Dh/2 (h = p/2 corresponds to the (O, x, y) plane) the choice D/ = Dh is made, therefore the number of equal-sized facets on that ring is simply N p=2 ¼

p . Dh

ð34Þ

The equal-size Seq of the facets is then computed by applying formula (33) to h = p/2, which gives
 S p=2 2p Dh S eq ¼ ¼ sin . 2 N p=2 N p=2

ð35Þ

For any ring h ± Dh/2, the number of iso-size facets, and therefore of /-increments is then given by the simple formula Nh ¼

Sh ¼ N p=2 sin h. S eq

ð36Þ

At each pole h = 0 and h = p, there is a single facet of area 
 Z Dh 2 Dh S0 ¼ p sinðcÞ dc ¼ p 1  cos . 2 0

ð37Þ

5.2. Algorithm for orientation averaging In 3D, orientation integral (4) is computed as follows: I

3D

¼

Z

p

h¼0

Z

2p

lðh; /Þwðh; /Þ sin h dh d/  n /¼0

Nf X

lðhif ; /if Þwðhif ; /if ÞSðhif ; /if Þ;

ð38Þ

if ¼1

where n = 2 for the discrete ODF and n = 1 for the fitted Gaussian. Nf designates the total number of facets on half the unit sphere computed by the algorithm of Section 5.1, and for each facet if defined by angles hif and /if , wðhif ; /if Þ is the ODF computed in Sections 4.1 and 4.2 and Sðhif ; /if Þ the area computed in Section 5.1. In 2D, orientation integral (4) is computed as follows: I2D ¼

Z

2p

lð/Þwð/Þ d/  nD/

/¼0

N/ X

lð/i/ Þwð/i/ Þ;

ð39Þ

i/ ¼0

where n = 2 for the discrete ODF and n = 1 for the fitted Gaussian, D/ is a constant /-angle increment, with a total number of increments N/ = p/D/, and for each discrete angle /i/ ¼ i/ D/, wð/i/ Þ is the corresponding ODFÕs value computed in Sections 4.1 and 4.2.

6. Numerical simulations The procedures presented in this paper were implemented in the DIGIMAT [16] software. DIGIMAT stands for Digital Materials and enables the numerical simulation of linear and nonlinear composites based on mean-field homogenization. The following examples illustrate the accuracy and usefulness of the methods described in this article.

1398

I. Doghri, L. Tinel / Comput. Methods Appl. Mech. Engrg. 195 (2006) 1387–1406

6.1. Effective CTE of a highly contrasted polymer matrix composite We consider a polymer matrix with YoungÕs modulus E0 = 0.35 GPa, PoissonÕs ratio m0 = 0.4 and CTE a0 = 104/K. It is reinforced with v1 = 30% volume fraction of glass fibers of aspect ratio Ar = 20 and thermo-elastic properties E1 = 73 GPa, m1 = 0.25 and a1 = 5 · 106/K. This is a severe test because the phasesÕ contrast is very important (E1/E0 = 208.6 and a0/a1 = 20). The effective stiffness and CTE predictions are reported in Table 1 w.r.t. a fixed Cartesian frame for several cases. We start with aligned fibers, a ¼ DIAGð1; 0; 0Þ. This is the first reference case. The M–T prediction of the stiffness corresponds to one of the Hashin and Shtrikman [28] bounds, and we checked that using those bounds, the predicted CTE values lie within the Schapery [50] bounds. Next, we consider misaligned fibers described by a ¼ DIAGð0:8; 0:1; 0:1Þ. It is seen in Table 1 that the 3D hybrid closure (Section 3.3) predicts a negative CTE in direction 1, which is unexpected and cannot be explained physically. The corresponding stiffness value seems also too far off. However, our proposed weighting procedure (Section 3.4) predicts positive CTE and stiffness values close to the aligned case. We point out that closure approximations influence both the effective stiffness and CTE (e.g., [44]). In order to better quantify the behaviour of different closure approximations, we consider a smaller perturbation around the aligned case a ¼ DIAGð0:98; 0:01; 0:01Þ. Table 1 shows that different closures (2D hybrid, 2D and 3D quadratic) predict CTE values whose relative differences with the reference aligned case are far too important to be reliable, while our proposed weighting formula gives much smaller and more intuitively acceptable differences with aligned fibers results. In order to see how the weighting formula performs in pure 3D cases, we consider 3D random orientation and a slight perturbation around it, i.e.,
 1 1 1 a ¼ DIAG ; ; ; a ¼ DIAGð0:343; 0:323; 0:333Þ. 3 3 3 Table 1 shows that both the purely 3D hybrid closure and our regularization formula give acceptable predictions. Finally, consider purely 2D cases, namely a 2D random orientation and a slight perturbation,
 1 1 a ¼ DIAG ; ; 0 ; a ¼ DIAGð0:49; 0:51; 0Þ. 2 2 From Table 1, it appears that the predictions of both the purely 2D hybrid closure and our weighting method are acceptable. The main advantage of the latter is to ensure a smooth transition between 1D, 2D and 3D cases and therefore deal with orientations which are not strictly in any of those cases. 6.2. Accuracy of orientation algorithms A first example checks the accuracy of ODF reconstruction in 2D. Consider a Gaussian, Eq. (31) with /M = p/2 and S = 0.25. First, we compute the corresponding a directly: Eq. (4) with lij = pipj, p1 = cos / and p2 = sin /. It is found that the only non-zero values are a11 = 0.9615 and a22 = 0.0384. Second, A is reconstructed from a (Section 3). Third, an ODF is reconstructed from a and A (Section 4). The initial Gaussian ODF (multiplied by 1/2) and the reconstructed ODF are plotted together in Fig. 2, and a good agreement can be seen.

I. Doghri, L. Tinel / Comput. Methods Appl. Mech. Engrg. 195 (2006) 1387–1406

1399

A second example shows the necessity of analytical fitting of a reconstructed ODF in some cases. Consider a 2D second-rank orientation tensor a11 = 0.88, a22 = 0.12 and a12 = 0.037. First, A is reconstructed from a (Section 3). Second, the ODF is reconstructed (Section 4). Fig. 3 shows that the recovered ODF 0.8

Given ODF Reconstructed ODF 0.7

0.6

ODF

0.5

0.4

0.3

0.2

0.1

0

0

0.5

1

1.5

φ [rad]

2

2.5

3

3.5

Fig. 2. Accuracy of ODF reconstruction. Starting with a given ODF, aij and Aijkl are computed and then an ODF is recovered. 0.6

Reconstructed ODF Gauss fitting of reconstructed ODF 0.5

0.4

ODF

0.3

0.2

0.1

0

–0.1 –2

–1.5

–1

–0.5

0

0.5

1

1.5

2

φ [rad]

Fig. 3. Analytical fitting of a reconstructed ODF. Starting with aij, Aijkl is computed, then an ODF is recovered. It has unacceptable negative values and is fitted with a Gaussian.

1400

I. Doghri, L. Tinel / Comput. Methods Appl. Mech. Engrg. 195 (2006) 1387–1406

have unacceptable negative values at some angles. Therefore, the relevant discrete data is fitted with a Gaussian, Eq. (31) whose parameters are /M = 0.05 rad and S = 0.3978. The good fitting is illustrated in Fig. 3 where the reconstructed ODF is superposed to the fitted Gaussian (multiplied by 1/2). A third example illustrates the accuracy of the discretization of the unit orientation sphere into a number of facets of equal areas. Consider a given fixed value Dh of the meridien angle increment. For each discrete angle value hi, the algorithm of Section 5.1 computes angle values /j over half the unit orientation sphere such that the latter is subdivided into facets of approximately equal areas. A measure of accuracy is to compare the sum of all facet areas to the exact area of half the unit sphere (i.e., 2p). This is illustrated in Table 2 which shows that acceptable accuracy is achieved for Dh = 5, i.e., 846 facets. A classical algorithm subdividing the half sphere with equal values Dh = D/ = 5 would create 1296 discrete orientations instead of 846. A fourth example illustrates the combination of all orientation algorithms in 3D and the final resultant accuracy. Consider a given second-rank orientation tensor a in 3D. First, the fourth-rank orientation tensor A is reconstructed via the hybrid closure approximation (Section 3). Second, the ODF is reconstructed from a and A using the formulae of Section 4. Actually, the ODF is computed at discrete orientations given by the iso-size facetsÕ algorithm of Section 5.1. Third, a is recomputed from the ODF using the discretization formulae of Section 5.2 with lij = pipj, p1 = sin h cos /, p2 = sin h sin / and p3 = cos h. The original and reconstructed values of aij are given in Table 3 for different values of Dh. Note that acceptable accuracy

Table 2 A measure of accuracy of the iso-size facets algorithm Number of h-increments h-increment Number of facets Sum of facet areas Accuracy (%)

12 15 100 5.9728 4.94

24 7.5 382 6.1195 2.60

36 5 846 6.1619 1.93

45 4 1312 6.1710 1.78

90 2 5196 6.2204 1.00

For each fixed value of Dh, a number of facets on half the unit sphere is computed and the sum of facet areas is compared to the exact value (2p).

Table 3 Accuracy of orientation treatment in 3D a11

a22

a33

a12

a23 0.005

a13

Given aij

0.338333

0.333333

0.328333

0.005

12 h-increments

0.324131 4.20%

0.319274 4.22%

0.312544 4.81%

0.004827 3.47%

0.002698 46.03%

0.001230 75.39%

0.005

24 h-increments

0.330951 2.18%

0.326050 2.18%

0.318073 3.12%

0.004897 2.06%

0.004326 13.48%

0.003943 21.14%

36 h-increments

0.332331 1.77%

0.327419 1.77%

0.321442 2.10%

0.004911 1.78%

0.004676 6.48%

0.004505 9.90%

45 h-increments

0.332840 1.62%

0.327920 1.62%

0.321737 2.01%

0.004920 1.60%

0.004814 3.72%

0.004741 5.18%

90 h-increments

0.335852 0.73%

0.330875 0.74%

0.323458 1.48%

0.004966 0.68%

0.004912 1.76%

0.004884 2.32%

Starting with a given aij, Aijkl is approximated, then an ODF is recovered and aij is recomputed using the iso-size facets algorithm, the ODF values for each facet and discrete orientation averaging. Initial and reconstructed values of aij, and their relative differences are shown.

I. Doghri, L. Tinel / Comput. Methods Appl. Mech. Engrg. 195 (2006) 1387–1406

1401

is achieved with 45 increments of Dh = 4, i.e., 1312 facets (Table 2), while a classical subdivision of half the sphere would have created 2025 discrete orientations. 6.3. Elasto-plastic 6061 Al matrix reinforced with 3D randomly oriented short SiC fibers An aluminum alloy (6061 Al) matrix is reinforced with vf = 20% of elastic SiC fibers of aspect ratio l/d = 10 and properties Ef = 455 GPa and mf = 0.15. The matrixÕ behaviour is elastic-J2 plastic with Em = 72.5 GPa, mm = 0.34, initial yield stress rY = 250 MPa and power-law isotropic hardening (R(p) = kpm) with k = 173 MPa and m = 0.46. A 3D random fiber orientation is assumed in the numerical simulations. Macro stress–strain curves under uniaxial macro tension predicted with our two-step M–T/ Voigt homogenization procedure are plotted in Fig. 4 for several values of the number of facets used in orientation discretization. The figure shows that the predicted response stiffens with the number of facets and compares the results with those of Dunn and Ledbetter [21]. 6.4. Elasto-plastic Al–5.5Mg matrix reinforced with 2D randomly oriented short d-Al2O3 fibers The matrix is an aluminum alloy (Al–5.5Mg) reinforced with vf = 10% of d-Al2O3 elastic fibers with l/ d = 20 and elastic properties Ef = 300 GPa and mf = 0.2. The matrix is elastic-J2 plastic with Em = 70.2 GPa, mm = 0.33, rY = 100 MPa, and power law isotropic hardening with k = 479 MPa and m = 0.36. The numerical simulations assume 2D random orientation. The predicted macro stress–strain response under uniaxial macro stress states is plotted in Fig. 5. Under monotonic tension, our M–T/Voigt predictions can be compared in Fig. 5 to experimental results reported by Kang et al. [33]. There is a good agreement especially when taking into account that several modeling assumptions (e.g., 2D randomness and constant l/d) are not well verified experimentally. The simulation of a first cycle under imposed macro strain

400

350

Macro Stress [MPa]

300

250

200

150

100 DIGIMAT, Number of facets : 100 DIGIMAT, Number of facets : 382 DIGIMAT, Number of facets : 1312 DIGIMAT, Number of facets : 5196 Dunn and Ledbetter (1997)

50

0 0

0.005

0.01

0.015

0.02

0.025

0.03

Macro Strain

Fig. 4. Macrotension tests of an elasto-plastic matrix reinforced with 3D randomly oriented fibers. Influence of the number of discrete facets on the Mori–Tanaka/Voigt predictions, and comparison with Dunn and Ledbetter [21] results.

1402

I. Doghri, L. Tinel / Comput. Methods Appl. Mech. Engrg. 195 (2006) 1387–1406 Al-5.5Mg matrix composite with 2D randomly oriented δ-Al2 O3 fibers (l/d = 20) under tension/compression cycle.

Experiment, Kang et al. (2002), Fig. 9 DIGIMAT, 2D random, Mori-Tanaka/Voigt

300

Macro stress [MPa]

200 100 0 –100 –200 –300 –0.01

–0.005

0

0.005

0.01

Macro strain Fig. 5. Macrotension/compression tests of an elasto-plastic matrix reinforced with 2D randomly oriented fibers. Predictions with Mori–Tanaka/Voigt homogenization. Comparison with experimental results of Kang et al. [33] in the monotonic part and simulation of a first cycle under imposed macrostrain.

(±1 · 102) is also plotted in Fig. 5 in order to illustrate the capability of our incremental formulation to handle unloading and cyclic loadings unlike a secant or total deformation theory such the one used in [21]. 6.5. Coupled two-scale elasto-plastic analysis of a turbine blade The procedures presented in this paper enable the simulation of products made of composite materials using a two-scale micro-to-macro approach. The example of the turbine blade of a jet engine (Figs. 6 and 7) is studied here. The blade is made of Titanium reinforced with silicon carbide short fibers. A J2 elasto-plastic model is used for Ti with the following properties: E0 = 110 GPa, m0 = 0.25, rY = 300 MPa and powerlaw isotropic hardening (k = 332 MPa and m = 0.18). The fibersÕ data are E1 = 410 GPa, m1 = 0.17, Ar = 10 and v1 = 8%. The blade is subject to a centrifugal load under imposed rotational velocity W having a triangular history. W increases linearly with time up to a maximum value WM = 900 rad/s (8595 RPM), then decreases

Fig. 6. View of a jet engine. Each turbine blade is supposed to be made of a Ti matrix reinforced with short SiC fibers.

I. Doghri, L. Tinel / Comput. Methods Appl. Mech. Engrg. 195 (2006) 1387–1406

1403

Fig. 7. Two-scale simulation of a turbine blade. FE analysis at the macro scale and homogenization at the microlevel for each macro integration point. Axial stress (MPa) versus axial strain at the most stressed point for the homogeneous and reinforced matrix.

linearly to a zero value. At the macroscopic level, a FE analysis is carried out with the ABAQUS [1] program. The macro mesh uses 24205 tetrahedron elements (of type C3D4) and 6344 nodes. A coupled two-scale procedure was developed by linking the homogenization software (DIGIMAT) to ABAQUS through its user-defined function UMAT. For each time step and each iteration of the FE equilibrium conditions, and at each integration (Gauss) point of the mesh, ABAQUS calls DIGIMAT with strains at that point as input. DIGIMAT views the point as the center of a RVE and the strains as macro strains. It carries out homogenization computations and returns macro values of stress and tangent stiffness to ABAQUS. Some results are illustrated in Fig. 7 where stress–strain curves in the axial (centrifugal) direction are plotted at the most stressed point for 4 cases: the unreinforced matrix and 3 fiber orientations. It is found that when fibers are aligned in the axial direction, the response remains elastic, but in all 3 other cases, plastic strains develop. Also, as expected the compositeÕs stiffness for 3D randomly oriented fibers lies between the axial and transverse fixed orientation results and again it is seen that our incremental formulation is able to simulate non-monotonic loadings.

7. Conclusions and prospects This paper focused on the numerical treatment of misaligned inclusionsÕ orientation in the mean-field homogenization of multi-phase elastic and inelastic composites. Some important conclusions are listed hereafter. Mean-field homogenization of multi-phase composites is carried out via a two-step procedure based on general orientation and volume averaging results. The procedure was worked out for linear thermo-elastic and rate-independent inelastic composites. In the latter case, our proposed incremental formulation is developed based on implicit time-discretization algorithms and consistent tangent operators. Unloading,

1404

I. Doghri, L. Tinel / Comput. Methods Appl. Mech. Engrg. 195 (2006) 1387–1406

cyclic and otherwise non-proportional loadings can be simulated, to the contrary to secant or total deformation theories. Usually, orientation averaging involves second- and fourth-rank orientation tensors a and A, resp., but only the former is available. The latter is reconstructed (approximately, except in special cases) from the closure formulae in 1D, 2D or 3D. We proposed a new weighting scheme which ensures a smooth transition between those cases. Numerical simulations show that without this formula, physically unrealistic predictions can be obtained for some composites. In several cases such as inelastic composites, the ODF is needed but is unavailable. It is therefore reconstructed from a and A. We implemented an efficient algorithm which computes the ODF at a number of discrete orientations. Analytical fitting of the discrete ODF values is proposed either for computational efficiency or for cases where physically unacceptable (e.g., negative) discrete ODF values are obtained. Again for inelastic behaviour and other instances, the computation of orientation averages for suitable microfields is needed. In 3D, an efficient algorithm approximates the integrals with sums over a finite number of iso-size facets on the surface of a unit sphere. The methods and algorithms proposed in this paper were implemented in the DIGIMAT [16] software and their accuracy and efficiency were illustrated through several numerical simulations. Some issues can be explored further in future work: (1) closure approximations other than the hybrid; (2) ODF reconstruction for cases which are strictly neither 2D nor 3D; (3) analytical fitting of the ODF with functions other than the Gaussian; (4) algorithms for orientation averaging involving smaller numbers of discrete orientations than the iso-size facets algorithm studied in this paper.

Acknowledgments The authors would like to thank the following institutions and initiatives for their support: (1) e-Xstream engineering SA (http://www.e-Xstream.com); (2) IUAP P5/08 project ‘‘From microstructure towards plastic behaviour of single- and multi-phase materials’’; (3) Re´gion Wallonne through FIRST projects.

References [1] ABAQUS, A general purpose finite element program, Abaqus, Inc., Pawtucket, RI, USA, 2004. [2] J. Aboudi, M.-J. Pindera, S.M. Arnold, Higher-order theory for periodic multiphase materials with inelastic phases, Int. J. Plast. 19 (6) (2003) 805–847. [3] D.F. Adams, Inelastic analysis of a unidirectional composite subjected to transverse normal loading, J. Compos. Mater. 4 (1970) 310–328. [4] S.G. Advani, C.L. Tucker, The use of tensors to describe and predict fiber orientation in short fiber composites, J. Rheol. 31 (8) (1987) 751–784. [5] S.G. Advani, C.L. Tucker, Closure approximations for three-dimensional structure tensors, J. Rheol. 34 (3) (1990) 367–386. [6] B.A. Bednarcyk, S.M. Arnold, J. Aboudi, M.-J. Pindera, Local field effects in titanium matrix composites subject to fiber-matrix debonding, Int. J. Plast. 20 (8–9) (2004) 1707–1737. [7] Y. Benveniste, G.J. Dvorak, T. Chen, On diagonal and elastic symmetry of the approximate effective stiffness tensor of heterogeneous media, J. Mech. Phys. Solids 39 (7) (1991) 927–946. [8] H.J. Bo¨hm, W. Han, Comparisons between three-dimensional and two-dimensional multi-particle unit cell models for particle reinforced metal matrix composites, Modell. Simul. Mater. Sci. Engrg. 9 (2001) 47–65. [9] M. Bornert, T. Bretheau, P. Gilormini (Eds.), Homoge´ne´isation en me´canique des mate´riaux, 2. Comportements non line´aires et proble`mes ouverts, HERMES Science, Paris, 2001. [10] C.W. Camacho, C.L. Tucker III, S. Yalvac, R.L. McGee, Stiffness and thermal expansion predictions for hybrid short fiber composites, Polym. Compos. 11 (4) (1990) 229–239. [11] N. Carrere, R. Valle, T. Bretheau, J.-L. Chaboche, Multiscale analysis of the transverse properties of Ti-based matrix composites reinforced by SiC fibres: from the grain scale to the macroscopic scale, Int. J. Plast. 20 (4–5) (2004) 783–810.

I. Doghri, L. Tinel / Comput. Methods Appl. Mech. Engrg. 195 (2006) 1387–1406

1405

[12] J.L. Chaboche, P. Kanoute´, Sur les approximations ‘‘isotrope’’ et ‘‘anisotrope’’ de lÕope´rateur tangent pour les me´thodes tangentes incre´mentale et affine, C.R. Me´canique (2003) 857–864. [13] J.L. Chaboche, P. Kanoute´, A. Roos, On the capabilities of mean-field approaches for the description of plasticity in metal matrix composites, Int. J. Plast. 21 (6) (2005) 1409–1434. [14] R.M. Christensen, F.M. Waals, Effective stiffness of randomly oriented fibre composites, J. Compos. Mater. 6 (1972) 518–532. [15] J.S. Cintra, C.L. Tucker, Orthotropic closure approximations for flow-induced fiber orientation, J. Rheol. 39 (6) (1995) 1095– 1122. [16] DIGIMAT version 1.4, Linear and nonlinear multi-scale material modeling software, e-Xstream engineering SA (http://www. e-Xstream.com), Louvain-la-Neuve, Belgium, 2004. [17] I. Doghri, Mechanics of Deformable Solids—Linear, Nonlinear, Analytical and Computational Aspects, Springer-Verlag, Berlin, 2000. [18] I. Doghri, C. Friebel, Effective elasto-plastic properties of inclusion-reinforced composites. Study of shape, orientation and cyclic response, Mech. Mater. 37 (2005) 45–68. [19] I. Doghri, A. Ouaar, Homogenization of two-phase elasto-plastic composite materials and structures: study of tangent operators, cyclic plasticity and numerical algorithms, Int. J. Solids Struct. 40 (7) (2003) 1681–1712. [20] I. Doghri, L. Tinel, Micromechanical modeling and computation of elasto-plastic materials reinforced with distributed-orientation fibers, Int. J. Plasticity 21 (2005) 1919–1940. [21] M.L. Dunn, H. Ledbetter, Elastic–plastic behavior of textured short-fiber composites, Acta Mater. 45 (8) (1997) 3327–3340. [22] D. Duschlbauer, H.E. Pettermann, H.J. Bo¨hm, Mori–Tanaka based evaluation of inclusion stresses in composites with nonaligned reinforcements, Scripta Mater. 48 (2003) 223–228. [23] G.J. Dvorak, Transformation-field analysis of inelastic composite materials, Proc. R. Soc. London, Ser. A 431 (1992) 89–110. [24] J.D. Eshelby, The determination of the elastic field of an ellipsoidal inclusion, and related problems, Proc. R. Soc. London, Ser. A 241 (1957) 376–396. [25] P. Gilormini, Insuffisance de lÕextension classique du mode`le autocohe´rent au comportement non line´aire, CRAS 320 (Se´rie IIb) (1995) 115–122. [26] C. Gonzalez, J. Llorca, A self-consistent approach to the elasto-plastic behaviour of two-phase materials including damage, J. Mech. Phys. Solids 48 (2000) 675–692. [27] R.M. Haj-Ali, A.H. Muliana, A multi-scale constitutive formulation for the nonlinear viscoelastic analysis of laminated composite materials and structures, Int. J. Solids Struct. 41 (2004) 3461–3490. [28] Z. Hashin, S. Shtrikman, A variational approach to the theory of the elastic behaviour of multiphase materials, J. Mech. Phys. Solids 11 (1963) 127–140. [29] T.J.R. Hughes, K.S. Pister, Consistent linearization in mechanics of solids and structures, Comput. Struct. 8 (1978) 391–397. [30] T. Iwamoto, Multiscale computational simulation of deformation behavior of TRIP steel with growth of martensitic particles in unit cell by asymptotic homogenization method, Int. J. Plast. 20 (4–5) (2004) 841–869. [31] S. Jansson, Homogenized nonlinear constitutive properties and local stress concentrations for composites with periodic internal structure, Int. J. Solids Struct. 17–29 (1992) 2181–2200. [32] B. Ji, T. Wang, Plastic constitutive behavior of short-fiber/particle reinforced composites, Int. J. Plast. 19 (5) (2003) 565–581. [33] G.-Z. Kang, C. Yang, J. Zhang, Tensile properties of randomly oriented short d-Al2O3 fiber reinforced aluminum alloy composites. I. Miscrostructure characteristics fracture mechanisms and strength prediction, Composites: Part A 33 (2002) 647– 656. [34] D.C. Lagoudas, A.C. Gavazzi, H. Nigam, Elastoplastic behavior of metal matrix composites based on incremental plasticity and the Mori–Tanaka averaging scheme, Comput. Mech. 8 (1991) 193–203. [35] F. Le´ne´, Damage constitutive relations for composite materials, Engrg. Frac. Mech. 25 (5–6) (1986) 713–728. [36] G. Li, P. Ponte Castan˜eda, Variational estimates for the elastoplastic response of particle-reinforced metal–matrix composites, in: M. Ostoja-Starzewki, I. Jasiuk (Eds.), Micromechanics of random media, Appl. Mech. Rev. 47 (1) (1994) S77–S94 (Part 2). [37] R. Masson, M. Bornert, P. Suquet, A. Zaoui, An affine formulation for the prediction of the effective properties of nonlinear composites and polycrystals, J. Mech. Phys. Solids 48 (2000) 1203–1227. [38] A. Molinari, S. Ahzi, A. Koudane, On the self-consistent modeling of elastic–plastic behavior of polycrystals, Mech. Mater. 26 (1997) 43–62. [39] T. Mori, K. Tanaka, Average stress in matrix and average elastic energy of materials with misfitting inclusions, Acta Mater. 21 (1973) 571–574. [40] S. Nemat-Nasser, M. Hori, Micromechanics: Overall Properties of Heterogeneous Materials, second ed., Elsevier Science Publishers, Amsterdam, 1999. [41] E.T. Onat, F.A. Leckie, Representation of mechanical behavior in the presence of changing internal structure, ASME J. Appl. Mech. 55 (1988) 1–10. [42] H.E. Pettermann, A.F. Planskensteiner, H.J. Bo¨hm, F.G. Rammerstorfer, A thermo-elasto-plastic constitutive law for inhomogeneous materials based on an incremental Mori–Tanaka approach, Comput. Struct. 71 (1999) 197–214.

1406

I. Doghri, L. Tinel / Comput. Methods Appl. Mech. Engrg. 195 (2006) 1387–1406

[43] O. Pierard, I. Doghri, An enhanced affine formulation and the corresponding numerical algorithms for the mean-field homogenization of elasto-viscoplastic composites, Int. J. Plast. 22 (2006) 131–157. [44] O. Pierard, C. Friebel, I. Doghri, Mean-field homogenization of multi-phase thermo-elastic composites: a general framework and its validation, Compos. Sci. Technol. 64 (2004) 1587–1603. [45] P. Ponte Castan˜eda, P. Suquet, Nonlinear composites and microstructure evolution, Mechanics for a New Millenium (Proc. ICTAM 2000-Chicago), Kluwer Academic Publishers, Dordrecht, The Netherlands, 2001, pp. 253–273. [46] Y.P. Qiu, G.J. Weng, Elastic constants of a polycrystal with transversely isotropic grains, and the influence of precipitates, Mech. Mater. 12 (1991) 1–15. [47] Y.P. Qiu, G.J. Weng, An energy approach to the plasticity of a two-phase composite containing aligned inclusions, J. Appl. Mech. 62 (4) (1995) 1039–1046. [48] A. Roos, J.-L. Chaboche, L. Ge´le´bart, J. Cre´pin, Multiscale modelling of titanium aluminides, Int. J. Plast. 20 (4–5) (2004) 811– 830. [49] D. Saraev, S. Schmauder, Finite element modeling of Al/SiCp metal matrix composites with particles aligned in stripes—a 2D–3D comparison, Int. J. Plast. 19 (6) (2003) 733–747. [50] R.A. Schapery, Thermal expansion coefficients of composite materials based on energy principles, J. Compos. Mater. 2–3 (1968) 380. [51] J. Segurado, J. Llorca, C. Gonzalez, On the accuracy of mean-field approaches to simulate the plastic deformation of composites, Scripta Mater. 46 (2002) 525–529. [52] J.C. Simo, R.L. Taylor, Consistent tangent operators for rate-independent elastoplasticity, Comput. Methods Appl. Mech. Engrg. 48 (1985) 101–118. [53] J.C. Simo, T.J.R. Hughes, Computational Inelasticity second ed., Springer-Verlag, New York, 2000. [54] G.P. Tandon, G.J. Weng, Average stress in the matrix and effective moduli of randomly oriented composites, Compos. Sci. Technol. 27 (1986) 111–132. [55] G.P. Tandon, G.J. Weng, A theory of particle-reinforced plasticity, J. Appl. Mech.—Trans. ASME 55 (1988) 126–135. [56] V. Verleye, F. Dupret, Prediction of fibre orientation in complex injection moulded parts, in: C.E. Altan, D.A. Siginer, W.E. Van Arsdale, A.N. Alexandrou (Eds.), Development in Non-Newtonian Flows, AMD-vol. 175, ASME, New York, 1993, pp. 139–163. [57] B. Weber, B. Kenmeugne, J.C. Clement, J.L. Robert, Improvements of multiaxial fatigue criteria computation for a strong reduction of calculation duration, Comput. Mater. Sci. 15 (1999) 381–399. [58] E. Weissenbek, H.J. Bo¨hm, F.G. Rammerstorfer, Micromechanical investigations of arrangement effects in particle reinforced metal matrix composites, Comput. Mater. Sci. 3 (1994) 263–278.