Homogenization of elastic–viscoplastic heterogeneous materials: Self-consistent and Mori-Tanaka schemes

International Journal of Plasticity 25 (2009) 1024–1048

Contents lists available at ScienceDirect

International Journal of Plasticity journal homepage: www.elsevier.com/locate/ijplas

Homogenization of elastic–viscoplastic heterogeneous materials: Self-consistent and Mori-Tanaka schemes S. Mercier *, A. Molinari Laboratoire de Physique et Mécanique des Matériaux, UMR CNRS 7554, Université Paul Verlaine Metz, Ile du Saulcy, 57045 Metz, France

a r t i c l e

i n f o

Article history: Received 6 February 2008 Received in final revised form 22 August 2008 Available online 18 September 2008

Keywords: Interaction law Mori-Tanaka model Self-consistent model Elastic–viscoplastic materials Composite materials

a b s t r a c t This paper deals with the prediction of the macroscopic behavior of a multiphase elastic–viscoplastic material. The proposed homogenization schemes are based on an interaction law postulated by Molinari et al. [Molinari, A., Ahzi, S., Kouddane, R. 1997. On the self-consistent modelling of elastic–plastic behavior of polycrystals. Mech. Mater., 26, 43–62]. Self-consistent schemes are developed to describe the behavior of disordered aggregates. The Mori-Tanaka approach is used to capture the behavior of composite materials, where one phase can be clearly identified as the matrix. The proposed schemes are developed within a general framework where compressible elasticity and anisotropy of the materials are taken into account. Inclusions can have various shapes and orientations. Illustrations of the homogenization procedure are given for a two-phase composite materials. Comparisons between results of the literature and predictions based on the interaction law are performed and have demonstrated the efficiency of the proposed homogenization schemes. Ó 2008 Elsevier Ltd. All rights reserved.

1. Introduction Homogenization is an efficient technique to predict the macroscopic behavior of composite materials or aggregates. Two stages are needed to capture the macroscopic response of heterogeneous materials. The first step is named localization and consists in solving the problem of an inclusion embedded into an infinite reference homogeneous medium (Eshelby problem). The second step is averaging and consists in connecting local fields to global ones. When disordered materials are

* Corresponding author. Tel.: +33 03 87 31 54 09; fax: +33 03 87 31 53 66. E-mail address: [email protected] (S. Mercier). 0749-6419/$ - see front matter Ó 2008 Elsevier Ltd. All rights reserved. doi:10.1016/j.ijplas.2008.08.006

S. Mercier, A. Molinari / International Journal of Plasticity 25 (2009) 1024–1048

1025

considered, the self-consistent scheme is appropriate. In that case, the reference medium used for solving the non-linear Eshelby problem has the behavior of the aggregate (the so-called homogeneous equivalent medium). When dealing with a composite material, where a matrix phase can be clearly identified, the matrix material is usually adopted for the reference medium and the well-known Mori-Tanaka scheme is obtained. The Eshelby problem has an exact solution when dealing with linear elasticity, see Eshelby (1957). But it has no analytical solution when materials have non-linear behavior. To overcome the difficulties, the non-linear behavior of the individual components can be linearized. In that case, an approximate solution based on the original Eshelby solution can be found. From this general principle, }ner numerous homogenization schemes have been proposed for different classes of materials, Kro (1961); Hill (1965) and Berveiller and Zaoui (1979) for elastoplasticity, Hutchinson (1976); Molinari et al. (1987) and Lebensohn and Tomé (1993) for viscoplasticity. Note that in all contributions cited above, only average information within phases are obtained. Indeed, an important aspect when dealing with materials having non-linear behavior is the fact that the strain (stress and strain rate) distribution within an inclusion is non-uniform. Some efforts have been devoted to account for second order moment of strain. Fluctuation of stress and strain within homogeneous phases have been analysed via second order moment approaches by Ortiz and Molinari (1988), Kreher (1990) and Kreher and Molinari (1993) for linear elastic heterogeneous materials and by Ponte-Castañeda and Suquet (1998) and Ponte-Castañeda (2002) for non-linear rigid viscoplastic materials (no elastic deformation). In this paper, the fluctuations of strain within phases are not accounted for. Nevertheless, to the authors knowledge, those theories dealing with second moment of strain and stress, have not been yet extended to elastic–viscoplastic materials. Elastic–viscoplastic materials are considered in the present contribution. The main problem with this type of behavior has been illustrated by Suquet (1985). He has demonstrated that the macroscopic behavior of an aggregate of incompressible linear viscoelastic phases whose behaviors are represented by a Maxwell law does not have in general a Maxwell behavior, although, as demonstrated by Li and Weng (1994), the overall response could be of Maxwell type if the constituent phases have identical relaxation time. Therefore, except for some cases, the mathematical description of the macroscopic behavior is not known a priori. This leads to difficulties in the homogenization procedure. Only few exact solutions exist, even with linear behavior. Hashin (1969) considered the Eshelby problem where the inclusion and the matrix have incompressible linear viscoelastic behaviors described by a Maxwell law. He was able to determine analytically the strain inside the inclusion as a function of the remote loading. The closed form solution was derived by using the Laplace transform and the correspondence principle. Using also Laplace transform, Rougier et al. (1994) have obtained the explicit solution for the overall behavior of a two-phase incompressible linear viscoelastic material. The aggregate was assumed to be perfectly disordered and a self-consistent scheme was used for averaging. In the same spirit, an extension to non-linear elastic–viscoplastic behavior has been proposed by Masson and Zaoui (1999). The non-linear behavior is first linearized via the affine formulation at each step of the deformation process. The Laplace transform is applied to the linearized problem and a symbolic problem is solved in the Laplace space. The inverse transform to the real time space is performed numerically with a collocation technique, see Schapery (1962). Note that the inverse transform has a high computational time cost. In addition, the approach by Laplace transform is only exact for non-ageing linear viscoelastic behavior. In the non-linear range where material parameters are time evolving, this induces approximations which are difficult to quantify. Recently, Pierard and Doghri (2006) adopted the affine method together with the Laplace transform to predict the overall behavior of two-phase elastic–viscoplastic composites. Predictions obtained with the Mori-Tanaka scheme were compared to finite element calculations performed on a 2D axisymmetric unit cell. The unit cell contained only one inclusion. The authors observed that the results are more accurate when the affine stiffness modulus of the reference medium, used for the evaluation of the Eshelby tensor, are chosen isotropic. When the anisotropy of the stiffness tensor is considered, the overall response is too stiff. Pierard et al. (2007) extended the previous work by considering a 3D unit cell. The representative volume element contained 30 spheres randomly and isotropically dispersed. Calculations were performed for two different volume fractions of inclusions, 0.15 and 0.3. The macroscopic response was well reproduced for all considered loading paths (uniaxial and biaxial tension, shear, cyclic

1026

S. Mercier, A. Molinari / International Journal of Plasticity 25 (2009) 1024–1048

loading). Nevertheless, the authors mentioned that the affine approximation underestimated the stress within inclusions. To deal with elastic–viscoplastic behaviors, other approaches have been explored. Asaro and Needleman (1985) have adopted the Lin model to predict texture development in polycrystals. Zouhal et al. (1996) and Tòth et al. (2000) used a Taylor–Lin approach to model the mechanical behavior of polycrystalline copper in cyclic tension or compression. In that case, no interaction law is needed since the macroscopic strain rate is applied to all grains. Recently, Mercier et al. (2007) have analysed the behavior of nanocrystalline materials with a Taylor–Lin approach. The inverse Hall Petch behavior was captured with this simple model. Despites its simplicity, the Taylor–Lin approach is able to reproduce the overall response of aggregates when the contrast between phases is not too pronounced. However, based on the assumption of uniform strain rate, this model is not able to give any information concerning the distribution of strain rate and strain within the different material phases. To overcome this weakness and improve the predictions of the overall response of the aggregate, several authors have proposed more refined homogenization schemes (for instance Mori-Tanaka or self-consistent models) for non-linear elastic–viscoplastic materials, Weng (1993); Kouddane et al. (1993); Turner et al. (1994); Molinari et al. (1997); Paquin et al. (1999); Molinari (2002); Sabar et al. (2002); Berbenni et al. (2007). To be more precise, Weng (1993) developed a self-consistent approach for polycrystalline materials devoted to creep. In his approach, the behavior of the homogeneous medium representing the polycrystal is replaced by a linear Maxwell solid with time dependent secant viscoplastic modulus. More recently, Li and Weng (1997) and Li and Weng (1998) have modeled the behavior of an elastic–viscoplastic two-phase composite under creep and constant strain-rate loading. For these two contributions, the overall response of the composite is derived using a Mori-Tanaka homogenization scheme. The non-linear behavior of the matrix is approximated by a linear viscoelastic comparison material, the shear viscosity of the comparison material is set equal to the shear viscosity of the non-linear viscoplastic matrix for the considered deformation stage. The problem being linearized, the correspondence principle is used. Note that the same type of approach has been adopted recently for the creep prediction of nanocrystalline materials, see Barai and Weng (2007). Paquin et al. (1999) and Sabar et al. (2002) have proposed different self-consistent schemes based on projection operator and translated fields. Both schemes are able to reproduce with good accuracy, the behavior of a two-phase viscoelastic composite. Applications to polycrystalline materials and polymeric composites are presented and enlight the relevance of their approaches when non-linear behaviors are considered. Kouddane et al. (1993) and Molinari et al. (1997) dealt with the response of an aggregate where all phases have a non-linear Maxwell behavior. They proposed an interaction law which includes both the instantaneous elastic response during a strain rate jump and the long range viscous behavior. To carry out the homogenization stage, they also postulated that the global behavior of the aggregate can be of the same type but with time evolving material parameters. This type of approach was used to predict the cyclic response of polycrystals, see for instance Abdul-Latif (2004) where the non-linear behavior of FCC polycrystals is analysed. In this paper, the model was tested for different complex cyclic loading paths under strain and stress controlled conditions. The main outcome of this work was the prediction of the ratcheting behavior of stainless steel 316L. To validate the interaction law proposed by Molinari et al. (1997), the inclusion problem for elastic–viscoplastic materials has been addressed by Mercier et al. (2005). In this paper, incompressible and compressible elasticity were considered. The theory was developed for anisotropic material responses, but applications were restricted to isotropic materials. The average strain rate in the inclusion was related to the far field loading prescribed on the matrix through the interaction law. Comparisons between predictions derived from the interaction law and finite element results were performed. A good agreement with numerical results was observed for various inclusion shapes, material parameters and loading conditions, as far as average strain rates and stresses are considered. The paper is organized as follows. Using the interaction law of Molinari et al. (1997), homogenization schemes (self-consistent or Mori-Tanaka) are proposed in Section 2, for multiphase materials having elastic–viscoplastic behaviors. The general framework is presented. Compressible elasticity, material anisotropy and ellipsoidal inclusions of different shapes and orientations are taken into account. The rigid viscoplastic case is treated separately, so that notations and basis of the theory

S. Mercier, A. Molinari / International Journal of Plasticity 25 (2009) 1024–1048

1027

are first presented in a simple case. In Section 3, a comparison is made between the two self-consistent schemes proposed in Section 2, for a two-phase incompressible linear viscoelastic material with perfect disorder. Self-consistent calculations are also performed for disordered non-linear elastic–viscoplastic materials. For composite materials with non-linear responses, the predictions of the MoriTanaka scheme are shown to be in good correlation with the finite element results of Pierard and Doghri (2006) and Pierard et al. (2007). Finally, the behavior of a 6.6 polyamide reinforced by spherical polypropylene inclusions is simulated using the proposed Mori-Tanaka approach and compared to experimental data. It is worth noting that the calculations are straightforward and need low computational time. 2. Homogenization The main purpose of the paper is to propose homogenization procedures for multiphase materials having elastic–viscoplastic behavior. Nevertheless, the foundation of the homogenization procedure is first presented for rigid viscoplastic materials so that it will become clear later how this simple case has been extended to elastic–viscoplastic materials. 2.1. Homogenization of rigid viscoplastic materials Consider an ellipsoidal inclusion embedded into an infinite homogeneous matrix. The strain rate prescribed to the matrix at infinity is D. The inclusion (respectively, the matrix) has a rigid viscoplastic incompressible behavior (elasticity neglected) defined by the strain-rate potential g (respectively, g o ). The average strain rate d within the inclusion can be estimated by the following interaction law, see Molinari et al. (1987):

  1 1 tg d  D ¼ Atg : ðs  SÞ o  Po

ð1Þ

where S ¼ ogodo ðDÞ is the deviatoric stress at infinity and s ¼ og ðdÞ is the deviatoric stress within the od inclusion associated to the average strain rate d. The fourth order tensor Ptg o is calculated with use of Green functions related to the tangent modulus of the matrix Atg o calculated at D, see for example Molinari (2002). The relationship (1) provides an approximate solution of the non-linear Eshelby inclusion problem for rigid viscoplastic materials. This analytical solution was obtained by replacing the matrix material response:

D1 ! S1 ¼

og o ðD1 Þ od

ð2Þ

by a first order Taylor expansion (tangent approximation) at D:

D1 ! Atg o ðDÞ : D1 þ So where the tangent modulus

ð3Þ Atg o

and So are given by:

2

Atg o ¼

o go og ðDÞ and So ¼ o ðDÞ  Atg o : D odod od

ð4Þ

The relationship (1) has a form similar to other interaction laws developed for different classes of behavior, as the incremental formulation of Hill (1965) for elastoplasticity and the incremental theory of Hutchinson (1976) for viscoplasticity. Nevertheless, the linearization upon which Eq. (1) is based appears to be different from the incremental approach developed by Hutchinson (1976) and the model results are different, see for instance, Molinari et al. (2004). When the remote loading D, inclusion shape and strain-rate potentials are given, the interaction law (1) provides a non-linear equation for the average strain rate d within the inclusion. Validation using finite element calculations has been performed in Molinari et al. (2004). It has been shown that the average strain rate d calculated in the inclusion using Eq. (1) is in good agreement with finite element results.

1028

S. Mercier, A. Molinari / International Journal of Plasticity 25 (2009) 1024–1048

The solution of the Eshelby problem obtained from Eq. (1) is the basic stone for homogenization schemes. A second important point is the definition of the salient reference medium. For a composite material with inclusions embedded in a connected matrix, the Mori-Tanaka approach is adopted. Each phase, represented by an ellipsoidal inclusion, is considered as surrounded by the reference medium chosen as the matrix material and subjected at infinity to the remote loading Dm . Dm and Sm are, respectively, the average strain rate and the associated deviatoric Cauchy stress tensors in the matrix material. We denote by (p) (p ¼ 1 . . . N) the label of phases included in the matrix and p ¼ 0 is the label of the matrix. The average strain rate dðpÞ in a given phase is obtained by solving the non-linear equation derived from Eq. (1):

   1 1 og   p tg dðpÞ  Dm ¼ Atg : dðpÞ  Sm o  Po od

ð5Þ

where g p is the potential of phase ðpÞ. The macroscopic deviatoric stress tensor Sm is derived from the o ðDm Þ. The tangent modulus Atg potential g o by: Sm ¼ og o of Eq. (4) is evaluated at Dm . Note that the relaod tionship (5) is also valid when the matrix phase is considered, providing dðoÞ ¼ Dm . By volume averaging of Eq. (5) over all phases and by satisfying the consistency relation D ¼ hdi, the macroscopic strain rate in the matrix Dm can be linked to the loading D applied to the composite material:

* Dm ¼ D 

+  1 1 og   p tg ðpÞ  Atg  P : d S m o o od

ð6Þ ðpÞ

For a given loading D, Eqs. (5) and (6) provide a set of N þ 1 relationships for the unknowns d , ðp ¼ 1 . . . NÞ and Dm . Let us consider now an heterogeneous and disordered material. The appropriate homogenization approaches are now the self-consistent schemes. For these models, the reference medium is the homogeneous equivalent medium (HEM) having the macroscopic response of the aggregate. Each inclusion representing a phase is embedded in the HEM. The response of the homogeneous equivalent medium is described by the macroscopic potential g hom . As mentioned previously, the response of the HEM is linearized by the tangent approximation specified in Eq. (3). The solution of the corresponding Eshelby problem is obtained from Eq. (1) where Atg o is the tangent modulus of the HEM defined from g hom . Note that g hom cannot be always defined in an analytical way. This leads to difficulties for the definition of the tangent modulus of the HEM. Therefore, different estimates of Atg o can be considered. For instance, one may adopt the definition given by the incremental self-consistent approach of Hutchinson (1976), also used in the affine model of Masson and Zaoui (1999):

Atg o ¼



h  i1 h  i1 1 tg tg tg tg tg tg atg : K þ P : a  A K þ P : a  A : p p p 0 0 0 0

ð7Þ

where atg p is the tangent modulus of phase ðpÞ given by:

atg p ¼

o2 g p  ðpÞ  d odod

ð8Þ

K is the fourth order tensor operating on symmetric traceless second order tensors: K ijkl ¼ 12 ðdik djl þ dil djk Þ  13 dij dkl . The working of the self-consistent model is the following. Consider a material made up of N different phases, subjected to a remote loading D. Each phase, labeled ðpÞ, is represented by an ellipsoidal inclusion embedded in the infinitely extended HEM subjected to the remote strain rate D (a priori difðpÞ ferent from D). For each phase (p), the average strain rate d is provided by the interaction law derived from Eq. (1):

dðpÞ  D ¼



  1 1 og   p tg Atg : dðpÞ  S o  Po od

ð9Þ

S. Mercier, A. Molinari / International Journal of Plasticity 25 (2009) 1024–1048

1029

The macroscopic deviatoric stress tensor is defined as the average of the deviatoric stress tensor in all phases: S ¼ hsi. Using the linearization of the response of the HEM at D, see Eq. (3), one obtains:

S  S ¼ hsi  S ¼ Atg o ðDÞ : ðD  D Þ

ð10Þ

The introduction of Eq. (10) into Eq. (9) leads to:

   1 1 og p tg ðpÞ tg dðpÞ  D ¼ Atg  P : Þ  hsi þ A : ðD  D Þ ðd  o o o od

ð11Þ

From the consistency equation hdi ¼ D and Eq. (9), the strain rate tensor D is linked to the strain rates d in the local phase, see Molinari and Mercier (2004):

D ¼

  1   1 1 tg P tg : P : d o o

ð12Þ

The N unknowns dðpÞ are calculated as solutions of Eq. (11) where Eq. (12) has been used. In the case of ellipsoidal inclusions with the same orientation and shape ratio, the tensor P tg o is the same for all phases. Then, from Eq. (12), D ¼ D and the interaction law (11) simplifies into:

dðpÞ  D ¼



 1 1   Atg  Ptg : sðpÞ  hsi o o

ð13Þ

Homogenization methods based on a tangent linearization of the matrix response (first order Taylor expansion with tangent modulus Atg o ) are named tangent approaches, see Molinari et al. (1987), Molinari (2002). The affine model of Masson and Zaoui (1999) is based on the same linearization by a first order Taylor expansion of the non-linear material response. Concerning the Mori-Tanaka approaches, both the tangent model or the affine model provide identical results. In addition, it has been shown in Molinari and Mercier (2004), that the self-consistent scheme based on the tangent interaction law with Atg o defined by Eq. (7) provides the same results as the affine method. However, the tangent approach developed here is broader as various estimates for the tangent modulus of the HEM can be proposed, see Molinari (1999). 2.2. Position of the problem for elastic–viscoelastic materials Consider an heterogeneous medium V made of different material phases. The global material is an aggregate of phases (distributed in a disordered or in an ordered manner). All phases have a non-linear elastic–viscoplastic behavior. The total strain rate tensor d at any material point of the medium V, is split into an elastic part de and an anelastic (viscoplastic) part dvp :

d ¼ de þ dvp The elastic part is linked to the Cauchy stress tensor

ð14Þ

r by the incremental elastic law:

de ¼ ðae Þ1 : r_

ð15Þ

e

where a is the fourth order tensor of elastic moduli and ð_Þ refers to the material time derivative of ð Þ. The present contribution is restricted to a small deformation theory. For finite deformation, an objective time derivative, for instance Jaumann or Green Naghdi, can be used in Eq. (15) instead of the material rate of the Cauchy stress tensor r_ . The viscoplastic contribution dvp , assumed to be volume preserving, is related to the deviatoric Cauchy stress tensor s by:

dvp ¼

of og ðsÞ or s ¼ ðdvp Þ os od

where f and g are conjugated stress or strain-rate potentials.

ð16Þ

1030

S. Mercier, A. Molinari / International Journal of Plasticity 25 (2009) 1024–1048

As for individual phases, the total strain rate D is split into elastic De and anelastic Dvp contributions:

D ¼ De þ Dvp

ð17Þ

The most important step in the homogenization process is to determine the local strain rate d in terms of the macroscopic strain rate D. 2.3. Interaction law Let us consider the Eshelby problem where an inclusion is embedded into a reference medium. The inclusion and reference materials have an elastic–viscoplastic behavior described by Eqs. (15) and (16). The remote loading prescribed on the reference medium at infinity is D . The strain rate in the inclusion d can be related to the remote loading D through an interaction law. In the proposed work, we adopt the one defined by Molinari (2002):

   1 1  1 1 tg d  D ¼ Atg : ðs  S Þ þ Aeo  Peo : ðr_  R_  Þ o  Po

ð18Þ

e where Atg o is the tangent viscoplastic stiffness tensor and Ao the elastic stiffness tensor of the reference e medium. The fourth order tensor Ptg o (respectively, P o ) is calculated with use of Green functions related e to Atg o (respectively, Ao ). R represents the Cauchy stress tensor at the remote boundary of the reference medium. S denotes the corresponding deviatoric Cauchy stress tensor. Eq. (18) provides the exact solution of the Eshelby problem for linear incompressible viscoelastic materials when inclusions are spherical and the materials (inclusion and matrix) are obeying to a Maxwell behavior. Note that this interaction law (18) has been validated by Mercier et al. (2005) for various inclusion shapes, non-linear material responses and loading conditions.

2.4. Effective properties 2.4.1. Disordered materials 2.4.1.1. Non-linear material responses. For an heterogeneous medium with disordered phases, a selfconsistent estimate is able to better represent the effective properties of the overall medium. As mentioned previously, the reference medium is chosen to have the properties of the effective material. In a first approach, the effective elastic stiffness Aeo is defined by Aeo ¼ Ae where Ae is given by the classical elastic self-consistent scheme for a purely elastic material:

Ae ¼

 h  i1 h  i1 1 ae : I þ Pe : ae  Ae : I þ Pe : ae  Ae

ð19Þ

The tangent modulus tensor Atg o is formulated using the incremental self-consistent scheme of tg with: Hutchinson (1976). Then, Atg o ¼ A

Atg ¼



h  i1 h  i1 1 : K þ Ptg : atg  Atg atg : K þ Ptg : atg  Atg

ð20Þ

where atg is derived from Eq. (16) by the relationship:

atg ¼

o2 g vp ðd Þ odod

ð21Þ

I is the fourth order unit tensor defined by Iijkl ¼ 12 ðdik djl þ dil djk Þ. The macroscopic stresses are defined by volume averaging over all phases:

S ¼ hsi R ¼ hri

ð22Þ

In this section, we restrict our attention to the special case where all inclusions have the same ellipsoidal shape with the same orientation of principal axes. In that case, D ¼ D, R ¼ R and the interaction law Eq. (18) becomes:

S. Mercier, A. Molinari / International Journal of Plasticity 25 (2009) 1024–1048

 1  1 d  D ¼ Atg  ðPtg Þ1 : ðs  SÞ þ Ae  ðPe Þ1 : ðr_  R_ Þ

1031

ð23Þ

The general case, where inclusion shapes or orientations of principal axes are different, is discussed in Appendix A. The numerical procedure used for this self-consistent scheme, named SC1, is presented in Appendix B. It will be shown in Section 3 that this approach is too stiff for linear viscoelastic materials when the behavior of the different phases is strongly different. Nevertheless, for polycrystals with moderate heterogeneities, Molinari et al. (1997) and Abdul-Latif (2004) have shown that this type of self-consistent approach provides accurate results. Another self-consistent scheme will be defined next to account for strong material heterogeneities. 2.4.1.2. Non-linear elastic–viscoplastic materials with same strain rate sensitivity. A second way of homogenization is proposed for non-linear elastic–viscoplastic materials (named SC2) for which the viscoplastic behavior of each phase is described by a strain potential of the form: eq

g ¼ hðp Þ

ðd Þmþ1 mþ1

ð24Þ

pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi R eq eq where d ¼ 2=3dvp : dvp denotes the equivalent plastic strain rate. p ¼ d dt is the cumulated p plastic strain. The function hð Þ represents the hardening behavior of the material and m is the strain rate sensitivity. For the present SC2 model, all phases have the same strain rate sensitivity m, which is supposed to be lower than unity (m 6 1). According to Eq. (24), the tensor of viscoplastic tangent moduli of each individual phase is:

atg ¼

  2 2 dvp dvp eq hðp Þðd Þm1 K þ ðm  1Þ eq  eq 3 3 d d

ð25Þ

The overall response of the multiphase material is assumed to be defined by the following constitutive law:

R_ ¼ Aeo ðtÞ : ðD  Aso ðtÞ1 : RÞ

ð26Þ

where the elastic modulus Aeo and the viscous modulus Aso are in general time dependent. The term Aso ðtÞ1 : R is not equal to the macroscopic viscoplastic strain rate. For a Maxwell material, Aeo would be time independent and Aso would be a function of R. Thus, the material response (26), with timedependent moduli, is not in general of Maxwell type. In addition, we assign to the modulus of the reference medium the following definition:

Aeo ðtÞ ¼ aðtÞAe

and Aso ðtÞ ¼ As

ð27Þ

where a is a real parameter and Ae is the elastic tensor given by Eq. (19). As is defined as the secant modulus of the homogenized viscous material with purely viscous phases characterized by the strain-rate potential g p ðp ¼ 1 . . . NÞ. For materials having the same strain rate sensitivity, Hutchinson (1976) has shown that secant and tangent macroscopic moduli As and Atg of an heterogeneous material are linked by:

As ¼

1 tg A m

ð28Þ

As a consequence, the following constitutive law is proposed for the overall material:



R_ ¼ aAe : D  mðAtg Þ1 : R



ð29Þ

Note that the constitutive law (29) is compatible with the limiting cases of pure elastic composite material (a ¼ 1 and m ¼ 0) and pure non-linear viscous composite material (a ¼ þ1). For the particular set of parameters: a ¼ 1, m ¼ 1, Ae ¼ 2le K and Atg ¼ 2ge K, Eq. (29) simplifies and leads to:

R_ þ

le S ¼ 2le D ge

ð30Þ

1032

S. Mercier, A. Molinari / International Journal of Plasticity 25 (2009) 1024–1048

For this given set of parameters, one observes that a linear Maxwell behavior is retrieved. In this section, we consider the general case where the ellipsoidal inclusions can have different shapes. The micro–macro transition is postulated to be governed by the interaction law (18) with Aeo ¼ aAe and Aso ¼ m1 Atg . We obtain:

 1 1 1 e d  D ¼ Atg  ðPtg Þ1 : ð s  S Þ þ A  ðPe Þ1 : ðr_  R_ Þ

a

ð31Þ

_ are reference mechanical fields used to solve the inclusion problem. The quantities D , S and R  Since the HEM has the mechanical behavior of the overall material (for a self-consistent approach), _ are linked by Eq. (29): D , S and R 



R_  ¼ aAe : D  mðAtg Þ1 : R



ð32Þ

The homogenization scheme SC2 is defined by the interaction law (31), the overall response (32) and the consistency equation D ¼ hdi together with the relation R ¼ hri. Before entering the details of the resolution of the proposed scheme, it is useful to verify that this approach restitutes the appropriate material response in the case of loading at high stress rates and for low loading rates. For high stress rates, the term with ðs  S Þ in Eq. (31) becomes negligible with respect to the other contributions involving the rates of stress and strain. Thus, Eq. (31) simplifies into:

d  D ¼

1

a

  1 1 Ae  P e : ðr_  R_ Þ

ð33Þ

The constitutive laws of the overall material (32) and of the different phases reduce, for the same reason, to:

R_  ¼ aAe : D r_ ¼ ae : d

ð34Þ

The solution corresponds to a pure elastic homogenization problem with a ¼ 1. For low stress rates, the interaction law (31) reduces to the rigid viscoplastic law (1). The solution corresponds to a pure non-linear viscous homogenization problem where the strain rates are purely anelastic. Thus, our homogenization procedure is consistent with the material responses expected at high and low stress rates. The algorithm to solve the SC2 scheme is presented in the following. We assume that the representative volume element contains N phases. At time t, the stress tensors r, R and R and the strain rate tensor D are supposed to be known. From the local behavior of the different phases Eqs. (15) and (16), the viscoplastic strain rate dvp , the tensors of elastic moduli ae and of viscoplastic tangent moduli atg are evaluated. Using relationships (19) and (20), the fourth order tensors Ae and Atg are obtained. The remaining unknowns are d, D and a. We have to determine N þ 2 quantities, some of them being tensors. The introduction of Eqs. (15) and (32) into the relationship (31) provides a new expression for the strain rate d in individual phases:

h  i1 1 d ¼ ae  a Ae  Pe " e

vp

: a :d

e

tg 1

þ aA : ðD  mðA Þ

  1 : : R  Þ  a Ae  P e

(

)#   1 1 tg tg D þ A  P : ðs  S Þ ð35Þ

The latter relation (35) provides N relationships. The consistency hdi ¼ D and the relation R ¼ hri _ and R _ permit the closure of equations. d, D and a being determined, the rate of the stress tensors r_ , R  are calculated from Eqs. (15), (29) and (32). Therefore, the stress tensors can be updated at time t þ Dt by using an Euler forward scheme with Dt being the time increment. The procedure is pursued up to the end of the deformation process. When all inclusions have the same ellipsoidal shape with the same principal axes, the algorithm for the restitution of this SC2 scheme can be simplified and is presented in Appendix C.

S. Mercier, A. Molinari / International Journal of Plasticity 25 (2009) 1024–1048

1033

2.4.1.3. Range of applications of the two self-consistent approaches. The first self-consistent scheme (SC1) provides a general frame for non-linear elastic–viscoplastic materials. It has already been adopted for the description of polycrystalline materials. The second self-consistent scheme is more specific since it is devoted to heterogeneous materials with phases having the same strain rate sensitivity. 2.4.2. Composite materials A Mori-Tanaka estimate is used to derive the effective properties of composite materials where inclusions are embedded in a matrix material. In that case, the reference medium is the matrix material. The interaction law becomes:

d  Dm ¼

   1 1  1 1 tg e Atg  P : ðs  S Þ þ A  Pem : ðr_  R_ m Þ m m m m

ð36Þ

where the superscript m represents quantities related to the matrix phase. Note that the matrix phase behavior is given by Eqs. (15) and (16). From hdi ¼ D, the strain rate in the matrix Dm is linked to D:

Dm ¼ D 

* + * +  1 1  1 1 tg e e _ _ Atg  P : ðs  S Þ  A  P : ð r  R Þ m m m m m m

ð37Þ

As for self-consistent schemes, the macroscopic stress tensor is obtained from volume averaging over all phases: R ¼ hri. The procedure for the Mori-Tanaka homogenization is presented in Appendix D. 3. Results In Section 2, different homogenization schemes for multiphase materials have been proposed. The associated framework is general in a sense that compressible elasticity and non-linear viscoplasticity can be handled by the present models. In this section, for illustration purpose, we restrict our attention to two-phase elastic–viscoplastic (or viscoelastic) materials. When considering disordered two-phase material, self-consistent approaches will be adopted to determine the overall behavior of the aggregate. For a composite material with inclusions embedded in a matrix material, the Mori-Tanaka scheme will be used. Only spherical inclusions will be considered. 3.1. Linear viscoelastic material, self-consistent modelling The material is a disordered aggregate made up of two viscoelastic phases, with isotropic and incompressible behavior, obeying to a Maxwell law. Therefore only deviatoric stress and strain rate tensors are considered. From Eq. (15), the rate of the deviatoric stress tensor is linked to the elastic strain rate tensor by:

s_ i ¼ 2li dei

ð38Þ

where li (i ¼ 1; 2) denotes the shear modulus of the phase i. The viscous behavior is Newtonian, given by the relationship:

si ¼ 2gi dvp i

ð39Þ

where gi (i ¼ 1; 2) represents the viscosity of the phase i. The aggregate is assumed to be isotropic. The effective tensors for the elastic and viscous moduli are depending upon the single parameters leff and geff :

Ae ¼ 2leff K;

Atg ¼ 2geff K

ð40Þ

For spherical inclusions, expressions for the tensors Pe and P tg are:

Pe ¼

1 K; 5leff

Ptg ¼

1 K 5geff

ð41Þ

1034

S. Mercier, A. Molinari / International Journal of Plasticity 25 (2009) 1024–1048

Both phases are assumed perfectly disordered in the aggregate, so a self-consistent approach is proposed to capture the global response of the material. Based on Eqs. (40) and (41), the effective shear modulus leff is solution of a quadratic equation, deduced from Eq. (19):

3





2



leff þ 2l2 f  3l2 ð1  f Þ þ 2l1 ð1  f Þ  3l1 f leff  2l1 l2 ¼ 0

ð42Þ

where f stands for the volume fraction of phase (1) in the global material. The viscosity geff is solution of Eq. (42) by replacing l1 by g1 and l2 by g2 . For the SC2 type of homogenization, combining Eqs. (29), (31), (38), (40), (41) and (68), the average total strain rate tensor di in each phase is given by:



 leff leff 2li þ 3aleff di ¼ 5aleff D  a eff S  a eff ðsi  SÞ þ 2li dvp i

g

ð43Þ

g

The parameter a is evaluated at each step of the loading by minimizing the difference between the volume average of the local strain rate hdi and the applied strain rate D. Under monotonic loading in tension or in compression, for a material with no initial internal stresses, the consistency hdi ¼ D can be fulfilled exactly. Indeed, one can easily prove that the conditions hd11 i ¼ D11 or hd22 i ¼ D22 or hd33 i ¼ D33 lead to the same definition of a. From hd11 i ¼ D11 together with Eq. (43), a second order algebraic equation is obtained for a: 2

 eff 2

a ðl Þ 6D11 

6Dvp 11



eff

þ al

þ 2l1 ð1  f Þðs2  S11 ÞÞ þ 6

( 

  1 5D11  2Dvp 11 2l2 f þ 2l1 ð1  f Þ  eff ð2l2 f ðs1  S11 Þ:

g



l

1vp 1 d11 f

2vp 2 d11 ð1

þl

 f Þ  ðl1 þ l2 ÞD11

) 

 vp  1 2vp þ 4l1 l2 d11 f þ d11 ð1  f Þ  D11 ¼ 0

ð44Þ

where s1 (respectively, s2 ) represents the component of the deviatoric Cauchy stress tensor s11 in phase 1 (respectively, in phase 2). Note that a similar equation can be derived when the loading is monotonic in shear. The update of all mechanical quantities is detailed in Appendix C. When dealing with the SC1 model, the interaction law (23) can be written in a form that leads to a set of two equations for the stress rate tensors in the two phases:

  vp 2l1 2l 1 2l 1 1 þ ð1  f Þ s_  ð1  f Þ s_ ¼ 2l1 ðD  d1 Þ  ð1  f Þðs1  s2 Þ 3leff 1 3leff 2 3geff   vp 2l2 2l2 2l2 f s_ þ 1 þ f s_ ¼ 2l2 ðD  d2 Þ  eff f ðs2  s1 Þ 3leff 1 3leff 2 3g

ð45aÞ ð45bÞ

To obtain Eq. (45), the relation: S ¼ f s1 þ ð1  f Þs2 and its time derivative: S_ ¼ f s_ 1 þ ð1  f Þs_ 2 have been used. The two ways of homogenization SC1 and SC2 have been tested for different sets of material parameters (li , gi ) under cyclic uniaxial loading. The applied strain rate tensor D at the remote boundary corresponds to uniaxial traction or compression:

2

1 0 6 D ¼ Do 4 0 0:5 0

0

3 0 7 0 5 0:5

ð46Þ

The predictions based on the two schemes SC1 and SC2 are compared to the exact solution obtained by Rougier et al. (1994) and presented in Paquin (1998) and Paquin et al. (1999). Fig. 1 presents the cyclic response of an aggregate made of two linear viscoelastic materials. The volume fraction of the two materials is f ¼ 0:5. The material parameters are: l1 ¼ 50 MPa, g1 ¼ 10 MPa s, l2 ¼ 250 MPa, g2 ¼ 1000 MPa s. The applied strain rate is Do ¼ 104 s1 . It is observed that the SC2 model with the adjusted parameter a provides accurate results. Predictions are almost superimposed with the exact solution obtained via the self-consistent approach of Rougier et al. (1994). The model SC1 proposed by Molinari (2002) overestimated the stress during the viscoelastic transition. Results obtained by Paquin (1998) and Paquin et al. (1999) are also close to the exact solution. Fig. 1 presents

S. Mercier, A. Molinari / International Journal of Plasticity 25 (2009) 1024–1048

1035

Deviatoric stress s11 (MPa)

0.04

0.02

Self consistent SC2, α parameter Self consistent SC1, Molinari (2002) Self consistent, Paquin et al. (1999) Maxwell behavior Exact solution, Rougier et al (1994)

0

-0.02

-0.04 0

0.001 0.002 Longitudinal strain ε11

0.003

Fig. 1. Predictions of the stress–strain response of a two-phase material obtained with various self-consistent schemes. The material is subjected to cyclic tension/compression. Results are compared to the exact solution of the literature. Both materials are linear viscoelastic. Material parameters are: l1 ¼ 50 MPa, g1 ¼ 10 MPa s, l2 ¼ 250 MPa, g2 ¼ 1000 MPa s. The applied strain rate is Do ¼ 104 s1 . The volume fraction of each phase is f ¼ 0:5.

also the mechanical response of the aggregate when the macroscopic material is assumed to have a Maxwell behavior, with two constant moduli leff and geff . In that case, the macroscopic response is given by:

 S_ ¼ 2leff D 

S 2geff

 ð47Þ

For this particular behavior, in uniaxial tension or compression, the component of the macroscopic deviatoric Cauchy stress tensor S11 can be deduced in a closed form by integration of Eq. (47). Adopting a Maxwell description for the response of the aggregate overestimates the stress level and the viscoelastic transition is not well captured. These results clearly show that the two homogeneous schemes SC1 and SC2 do not finally lead to a Maxwell behavior for the aggregate. Fig. 2 shows the evolution of the parameter a during straining. At the beginning of the loading, the response in all phases is purely elastic; then a ¼ 1. As the deformation goes on, a decreases continuously and saturates. When the loading is reversed, a peak is observed to occur at the value a ¼ 1. From Fig. 1, one observes that the slope of the stress–strain curve when the loading is reversed, is identical for the four models. It can be easily shown that the slope is twice the initial slope observed when the material is stress-free (S11 ¼ 0, 11 ¼ 0). Therefore, even if a ¼ 1 when the loading is reversed, the material does not respond in a purely elastic way but in a viscoelastic one. Fig. 3(a) presents other results for the following parameters: l1 ¼ 50 MPa, g1 ¼ 1000 MPa s, l2 ¼ 1000 MPa, g2 ¼ 1000 MPa s. The applied strain rate is still Do ¼ 104 s1 . Both phases have the same volume fraction f ¼ 0:5. In the present test, the shear moduli are strongly different, the viscosity being identical for the two phases. It is observed that accurate results are obtained with all models. In Fig 3(b), the material parameters are: l1 ¼ 50 MPa, g1 ¼ 10 MPa s, l2 ¼ 50 MPa, g2 ¼ 1000 MPa s. The elastic shear moduli are identicals for both phases. Phase 2 has a viscosity hundred times larger than phase 1. In that case, the SC2 model gives good predictions while the SC1 model overestimates the stress level. An interesting conclusion (results not presented here) is that for a moderate heterogeneity on viscous and elastic properties, the SC1 and SC2 schemes provide both good results. This feature will be confirmed for non-linear behaviors. In a single phase polycrystal, the viscous and elastic moduli of crystals have frequently a small contrast and the SC1 scheme will provide accurate results.

1036

S. Mercier, A. Molinari / International Journal of Plasticity 25 (2009) 1024–1048

1.1 1 0.9

α parameter

0.8 0.7 0.6 0.5 0.4 0.3 Self consistent SC2, α parameter

0.2 0.1 0

0

0.001

0.002

0.003

Longitudinal strain ε11 Fig. 2. Evolution of the a parameter during the loading of Fig. 1.

3.2. Non-linear viscoplastic behavior, self-consistent modelling In this section, the viscoplastic behavior of the two phases is described by Eq. (24) with a strain rate sensitivity below unity. Both phases present a strain hardening behavior described by:

hðp Þ ¼ jðp þ  Þn

ð48Þ

where  is a reference strain and n is the strain hardening exponent. The two materials are elastically incompressible. The general framework of the SC1 and SC2 schemes for non-linear behavior is similar to the one developed for the linear case, see Section 3.1. The only difference lies in the definition of Atg . No simple analytical expression is obtained for Atg from the relationship (20). This equation is solved numerically and the associated tensor P tg is also calculated by numerical means. Fig. 4 presents the cyclic response of a two-phase aggregate. The volume fraction is f ¼ 0:5. The applied loading is Do ¼ 104 s1 . The material parameters are: l1 ¼ 81 GPa, j1 ¼ 1200 IS, n1 ¼ 0:05, l2 ¼ 40:5 GPa, j2 ¼ 600 IS, n2 ¼ 0:1, 1 ¼ 2 ¼ 105 . The strain rate sensitivity of both phases is m ¼ 0:3. No isotropic estimate of the viscoplastic tangent modulus tensor has been used, as in Pierard et al. (2007). The prediction of the two self-consistent schemes SC1 and SC2 are in good correlation. In the example presented in Fig. 5, a low strain rate sensitivity is adopted, m ¼ 0:03, consistently with the value of real strain rate sensitivity in standard metallic materials at room temperature. The material parameters are similar to those of Fig. 4 except for the shear modulus l2 ¼ l1 ¼ 81 GPa (corresponding to the shear modulus of a steel). The elastic response and the long term viscoplastic behavior are identically predicted by the two SC1 and SC2 schemes. A difference exists just after yielding, see Fig. 5. This situation is clearly depicted in Fig. 6. Since the two phases have the same shear modulus and since j1 ij2 , the phase 2 plastifies first while the phase 1 remains elastic during a certain time. During this stage (plastification of phase 2 only), the macroscopic response seems to be also linear. In fact, a transition occurs until the phase 1 presents also some plastic activity. Afterwards, the macroscopic response of the aggregate presents a low hardening rate, owing to the low value of n1 and n2 , see Fig. 6. 3.3. Non-linear composite, Mori-Tanaka scheme The predictions of the proposed elastic–viscoplastic Mori-Tanaka scheme are compared to finite element results obtained on 2D or 3D unit-cell models, see Pierard and Doghri (2006) and Pierard et al. (2007). First, a Perzyna-type elastic–viscoplastic constitutive model is used to describe the

1037

S. Mercier, A. Molinari / International Journal of Plasticity 25 (2009) 1024–1048

0.2

Deviatoric stress s11 (MPa)

0.15 0.1 0.05 Self consistent SC2, α parameter Self consistent SC1, Molinari (2002) Self consistent, Paquin (1998) Maxwell behavior Exact solution, Rougier et al (1994)

0 -0.05 -0.1 -0.15 -0.2 0

0.001

0.002 0.003 0.004 0.005 Longitudinal strain ε11

0.006

0.007

Deviatoric stress s11 (MPa)

0.04

0.02

Self consistent SC2, α parameter Self consistent SC1, Molinari (2002) Self consistent, Paquin (1998) Maxwell behavior Exact solution, Rougier et al (1994)

0

-0.02

-0.04 0

0.002

0.004

0.006

0.008

0.01

Longitudinal strain ε11 Fig. 3. Cyclic response of a two-phase material. Predictions are based on various self-consistent models. Materials have linear viscoelastic behavior. (a) The two phases have different shear moduli: l1 ¼ 50 MPa, l2 ¼ 1000 MPa and the same viscosity g1 ¼ g2 ¼ 1000 MPa s. The applied strain rate is Do ¼ 104 s1 . The volume fraction of each phases is f ¼ 0:5. (b) The two phases have the same shear moduli: l1 ¼ l2 ¼ 50 MPa and different viscosities g1 ¼ 10 MPa s and g2 ¼ 1000 MPa s.

behavior of both phases. The incremental elastic law is governed by Eq. (15) and the material is assumed to be elastically compressible. The viscoplastic response is described as in Perzyna (1986):

dvp ¼ U

oF or

ð49Þ

where U is the viscoplastic function. F is the quasi-static yield function, defined by:

F ¼ req  ry  Rðp Þ

ð50Þ

pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi where req ¼ 3=2s : s is the equivalent von Mises stress, ry the initial yield stress and Rðp Þ the hardeq ening function. Since a J 2 flow theory is adopted, it appears that U ¼ d . We assume that U is given by:

1038

S. Mercier, A. Molinari / International Journal of Plasticity 25 (2009) 1024–1048

40

Self consistent SC2, α parameter Self consistent SC1, Molinari (2002)

Deviatoric stress s11 (MPa)

30 20 10 0 -10 -20

m=0.3

-30 0

0.01 Longitudinal strain ε11

0.02

Fig. 4. Comparison of the predictions obtained with the proposed self-consistent schemes SC1 and SC2. A two-phase material is considered. The aggregate is subjected to cyclic tension/compression. The behavior of each phase is non-linear and presents a strain hardening. The material parameters are: f ¼ 0:5, Do ¼ 104 s1 . l1 ¼ 81 GPa, j1 ¼ 1200 IS, n1 ¼ 0:05, l2 ¼ 40:5 GPa, j2 ¼ 600 IS, n2 ¼ 0:1, 1 ¼ 2 ¼ 105 . The common strain rate sensitivity is m ¼ 0:3.

500

Self consistent SC2, α parameter Self consistent SC1, Molinari (2002)

Deviatoric stress s11 (MPa)

400 300 200 100 0 -100 -200 -300

m=0.03

-400 0

0.01 Longitudinal strain ε11

0.02

Fig. 5. Comparison of the predictions obtained with the self-consistent schemes SC1 and SC2. A two-phase material is considered. The aggregate is subjected to cyclic tension/compression. The behavior of each phase is non-linear and presents a strain hardening, identical to Fig. 4 except for the strain rate sensitivity which is here m ¼ 0:03.



U¼j

M F p ry þ Rð Þ

if

F i0;

U ¼ 0 otherwise

1=M represents the strain rate sensitivity of the material. (49) to (51), the flow stress of the material is:

ð51Þ

j is a scalar parameter. Combining Eqs.

1039

S. Mercier, A. Molinari / International Journal of Plasticity 25 (2009) 1024–1048

Self consistent SC2, α parameter, S11 Viscoplastic strain rate, phase 1 Viscoplastic strain rate, phase 2

Deviatoric stress s11 (MPa)

400

2

300 1

200 100

0

0 -100

-1

-200 -300

m=0.03

-400 0

100

200

300

400

500

600

700

Viscoplastic strain rate d11 (10-4 s-1)

500

-2 800

time (s) Fig. 6. Time evolution of the predicted stress for the aggregate defined in Fig. 5. The variation in the macroscopic stress is linked to the viscoplastic activities within the two phases. The SC2 self-consistent scheme is adopted.

p

req ¼ ðry þ Rð ÞÞ 1 þ

 eq M1 ! d

j

ð52Þ

In the following calculations, a power law model is adopted to characterize the hardening behavior of the material:

Rðp Þ ¼ kðp Þn

ð53Þ

n is the hardening exponent and k a reference stress, representing the stress increase when a plastic strain of one is cumulated. Note that Pierard et al. (2007) adopted the Perzyna-type of behavior since this model is implemented in the finite element software ABAQUS standard. In this part, we have used the same material description for validation purpose. Note that only spherical inclusions will be considered. Nevertheless, the present formulation of the Mori-Tanaka scheme based on the interaction law (36) is applicable to any elastic–viscoplastic constitutive response and to any ellipsoidal inclusion shapes. Fig. 7 presents the cyclic response of a composite material with 10% spherical inclusions. The behavior of the matrix and of the inclusions is described by the set of relationships (15) and (49) to (53). Material parameters of the matrix and of the inclusions are listed in Table 1. Note that the strain rate sensitivity of both phases is one. The composite material sustains cyclic uniaxial tension/compression up to 1% strain. The uniaxial strain rate D11 is equal to 103 s1 . Since the material is compressible, the component of the strain rate tensor D22 and D33 are adjusted at each step of the deformation so as to ensure that the stress components R22 and R33 are equal to zero. A very good agreement is observed between finite element results of Pierard and Doghri (2006) and the calculations obtained by the present Mori-Tanaka method based on the interaction law of Molinari (2002). The unloading and also the elastic–viscoplastic transitions are well captured by the proposed scheme. The approach based on the affine method and developed in Pierard and Doghri (2006) provides also accurate results. The main advantages of our homogenization scheme is that all calculations are performed in the real space-time without recoursing to any Laplace transform. In the second example, a composite material with a larger volume fraction of inclusions is considered: f ¼ 0:3. Both phases have an elastic–viscoplastic behavior described by the Perzyna model. The parameters are given in Table 2, see Pierard et al. (2007). The strain rate sensitivity 1=M is 0.667 for both materials. The Poisson’s ratio of the two phases are different: m ¼ 0:33 for the matrix and 0:286

1040

S. Mercier, A. Molinari / International Journal of Plasticity 25 (2009) 1024–1048

1000 Present Model, Molinari (2002) Abaqus, Pierard and Doghri (2006)

Stress (MPa)

500

0

-500 f=0.1 Uniaxial Tension/Compression at D11=10-3 s-1 -1000 -0.01

-0.005

0

0.005

0.01

Strain Fig. 7. Cyclic uniaxial tension/compression test. Results obtained using the proposed elastic–viscoplastic Mori-Tanaka scheme are close to F.E. results of Pierard and Doghri (2006). The material parameters are listed in Table 1. The applied strain rate is D11 ¼ 103 s1 . The volume fraction of inclusions is f ¼ 0:1.

Table 1 Material parameters for the matrix and the inclusions used for calculations of Fig. 7, see Pierard and Doghri (2006) Material

Young modulus (GPa)

Poisson’s ratio

ry (MPa)

k (GPa)

n

j (s1 )

M

Matrix Inclusion

100 500

0.3 0.3

100 500

5 5

1 1

0:3  103 0:3  103

1 1

Table 2 Material parameters for the matrix and the inclusions used for calculations of Fig. 9, see Pierard et al. (2007) Material

Young modulus (GPa)

Poisson’s ratio

ry (MPa)

k (GPa)

n

j (s1 )

M

Matrix Inclusion

70 400

0.33 0.286

70 400

4 8

0.4 0.4

0:3  103 0:2  103

1.5 1.5

for the inclusions. The effect of the applied strain rate on the uniaxial tensile response of the composite is investigated in Fig. 8. For the two strain rates, 103 s1 and 104 s1 , a perfect agreement is observed between finite element results and the proposed homogenization scheme. The scheme based on the affine method can also capture with good accuracy the mechanical response at low strain rate. A second advantage of our homogenization scheme is that the calculations are carried out with a full anisotropic tangent tensor. No isotropic estimate of the matrix tangent modulus is involved to get a better accuracy of the results as in Pierard et al. (2007). The same composite material as in Fig. 8 is now facing cyclic loading in tension and compression, see Fig. 9. The maximum strain is 0:05. For the two strain rates 106 s1 and 103 s1 , both schemes (affine method or Molinari’s interaction law) capture accurately the mechanical response. It is important for a micro–macro model to be able to reproduce the overall behavior of the aggregate. It is, also of interest, to propose good estimates of the stress and strain levels in all phases. In Fig. 10, the average stress in the inclusion and in the matrix predicted by the present scheme are compared to finite element results and to results based on the affine method and developed by Pierard et al. (2007). During uniaxial tension, at 103 s1 , the mean stress in the matrix material is clearly well

1041

S. Mercier, A. Molinari / International Journal of Plasticity 25 (2009) 1024–1048

3500

Present model, Molinari (2002) Affine method, Pierard et al (2007) Abaqus, Pierard et al (2007)

Equivalent stress (MPa)

3000 D11 =10 -3 s -1

2500 2000 1500

D11 =10

-4

s

-1

1000 500

Uniaxial tension, f=0.3

0 0

0.01

0.02

0.03

0.04

0.05

Longitudinal strain Fig. 8. Tensile response of a two-phase composite, effect of the strain rate. Results obtained using the proposed elastic– viscoplastic Mori-Tanaka scheme are close to F.E. results of Pierard et al. (2007). The material parameters are listed in Table 2. Two strain rates have been adopted D11 ¼ 104 s1 and D11 ¼ 103 s1 . The volume fraction of inclusions is f ¼ 0:3.

predicted for both schemes. For the von Mises stress in the inclusion, a 6% underestimation is provided by our approach. A 21% difference is observed with the method of Pierard et al. (2007). The difference in the stress estimate is linked to an underestimation of the cumulated plastic strain. Table 3 presents the plastic strain in the matrix and in the inclusions for two strain rates 103 s1 and 106 s1 , at the end of the cyclic loading of Fig. 9. For the larger loading rate 103 s1 , the present approach predicts a mean cumulated plastic strain of 8.02% in the matrix (to be compared to 8.23% for finite element calculations) and of 4.75% in the inclusions (5.08% for finite element calculations). As already mentioned in Fig. 10, the response of the matrix is accurate for both schemes. For the inclusion, our approach provides better results. For the fastest loading, a 7% underestimation is obtained for the plastic strain in the inclusion while Pierard et al. (2007) obtained a 23% difference. Tests depicted in Fig. 7 to Fig. 10 have validated the proposed elastic–viscoplastic Mori-Tanaka scheme. Without Laplace transform and without isotropic estimate of the matrix viscoplastic tangent stiffness, the results obtained with the present approach are accurate. Finally, our homogenization scheme is used to predict the behavior of a polymeric composite. The matrix is made of polyamide 6.6 (PA). Spherical inclusions of polypropylene are isotropically dispersed in the matrix material. The volume fraction of inclusions is f ¼ 0:25. Both phases are assumed incompressible. Beurthey (1997) has identified the behavior of both phases and of the composite, at different strain rates ranging from 103 s1 to 1 s1 . Based on experiments, the flow stress of the polyamide 6.6 is found to be rate dependent with strain hardening and can be described by a powerlaw function:

Table 3 Average cumulated plastic strain in the matrix and in the inclusions at the end of the loading of Fig. 9 Strain rate (s1 )

Method

Matrix (%)

Inclusions (%)

106 106 106

Finite element Present model Affine

18.9 18.74 19.5

9.8 8.47 6.53

103 103 103

Finite element Present model Affine

8.23 8.02 8.18

5.08 4.75 3.94

1042

S. Mercier, A. Molinari / International Journal of Plasticity 25 (2009) 1024–1048

Abaqus, Pierard et al (2007) Present model, Molinari (2002) Mori-Tanaka, Pierard et al (2007)

3000

Longitudinal stress (MPa)

2000

1000

0

-1000

f=0.3

-2000

Uniaxial tension/compression at D11 = 10 -0.04

-0.02

0.02

s

-1

0.04

Abaqus, Pierard et al (2007) Present model, Molinari (2002) Mori-Tanaka, Pierard et al (2007)

4000 Longitudinal stress (MPa)

0 Longitudinal strain

-6

2000

0

-2000

f=0.3

-4000

Uniaxial tension/compression at D11 = 10 -0.04

-0.02

0

0.02

-3

s -1

0.04

Longitudinal strain Fig. 9. Effect of the applied strain rate on the cyclic response of an elastic–viscoplastic two-phase material. The material parameters are listed in Table 2. The volume fraction is f ¼ 0:3. For both applied strain rates 106 s1 and 103 s1 , the Mori Tanaka schemes based on the affine method Pierard et al. (2007) or on the interaction law proposed by Molinari (2002) give accurate predictions. eq m p n rPA y ¼ kðd Þ ðo þ  Þ

ð54Þ

with k ¼ 103  106 IS, m ¼ 0:022, n ¼ 0:089 and o ¼ 1010 . The shear modulus of polyamide 6.6 is 800 MPa. The viscoplastic response of polypropylene presents a saturation of the hardening with strain. The following phenomenological relationship is adopted for the flow stress: eq m p rPP y ¼ ½ro þ ðrs  ro Þð1  a expðb ÞÞðd Þ 6

6

ð55Þ

with ro ¼ 15  10 IS, rs ¼ 48  10 IS, a ¼ 0:5, b ¼ 60 and m ¼ 0:04. The shear modulus of the polypropylene material is 340 MPa. Note that the plastic behavior of both phases (PA and PP) is well

1043

S. Mercier, A. Molinari / International Journal of Plasticity 25 (2009) 1024–1048

5000 Present model, Molinari (2002) Affine method, Pierard et al (2007) Abaqus, Pierard et al (2007) Equivalent stress (MPa)

4000

Inclusion

3000

2000

Matrix

1000

f=0.3 Uniaxial tension at D11 =10 -3 s -1

0 0

0.01

0.02

0.03

0.04

0.05

Longitudinal strain Fig. 10. Evolution of the average equivalent stress in the matrix and in the inclusion during tensile loading. Material parameters and configurations are similar to Fig. 9. It is shown that the proposed Mori-Tanaka scheme is able to capture the heterogeneity of stress and strain between phases.

80 70

Uniaxial stress (MPa)

60 50 40 30 -1

1s , 10-2-1 s-1 , -1 10 s , -3 -1 10 s ,

20 10 0

0

1 s-1, Present model 10-1 s-1, Present model 10-2 s-1, Present model -3 10 s-1, Present model Experiments, Beurthey (1997) Experiments, Beurthey (1997) Experiments, Beurthey (1997) Experiments, Beurthey (1997)

0.05

0.1

0.15

Uniaxial strain Fig. 11. Tensile response of a polymeric composite at different strain rates, comparisons with experiments of Beurthey (1997). The matrix is made of Polyamide 6.6. Spherical inclusions are made of polypropylene. The volume fraction of inclusions is 0.25.

described by relationships (54) and (55) for strain rates in the range ½103 s1  1 s1  (results not presented here). Since inclusions of PP are embedded into a PA matrix, the present Mori-Tanaka approach based on the interaction law (36) is adopted to model the macroscopic stress–strain curve. Fig. 11 shows that the overall response of the polymeric composite is well captured up to 0.15 uniaxial strain for the three smallest strain rates (103 s1 to 101 s1 ). A certain discrepancy is met for the fastest strain rate 1 s1 . From Fig. 11, it is observed that the proposed scheme is also efficient for composite material having different strain rate sensitivities and different strain hardening behaviors.

1044

S. Mercier, A. Molinari / International Journal of Plasticity 25 (2009) 1024–1048

4. Conclusion Homogenization for elastic–viscoplastic materials is a difficult problem that has received in the last decade lots of attention. Some of the contributors adopt the affine method together with the use of Laplace transform to replace the original problem by a symbolic thermo-elastic problem in the Laplace space. We have adopted a different point of view. Based on the interaction law, defined by Molinari (2002) and validated on the Eshelby problem by Mercier et al. (2005), we have proposed two self-consistent schemes for perfectly disordered materials and a Mori-Tanaka model for composite materials. The main advantages of the proposed approaches by comparison with schemes based on the affine method are twofold. First, the method is not computer time consuming since no Laplace transform is needed. Besides, the inverse Laplace transform induces approximations that are difficult to quantify, when non-linear behaviors are considered. Second, Pierard et al. (2007) have observed that the quality of the results is depending on the choice of the isotropic restriction for the tangent modulus of the reference medium. In our approach, the regular anisotropic tangent modulus is used. The first self-consistent scheme SC1 is valid for any non-linear behavior. The second one SC2 is dedicated to aggregates with phases having the same strain rate sensitivity. It is observed that predictions based on the second self-consistent model SC2 are accurate as compared to exact solutions obtained by Rougier et al. (1994) for a two-phase linear viscoelastic aggregate. Both schemes SC1 and SC2 provide similar results when non-linear behavior are considered, at least for the cases presented in this paper. For composite materials, the proposed Mori-Tanaka scheme compares very closely to numerical results obtained via ABAQUS, see Pierard et al. (2007). The overall response of the composite material is predicted accurately. In addition, the scheme is able to capture the strain and stress histories in the phases. This is crucial when dealing for example with residual stresses after metal forming operation. In the present paper, the tensile response of a polymeric composite has been modeled, both phases having low strain rate sensitivities. This shows that our Mori-Tanaka scheme based on the interaction law of Molinari (2002) is valid for various non-linear behaviors. Acknowledgements This work was supported in part by the EU Network of Excellence project Knowledge-based Multicomponent Materials for Durable and Safe Performance (KMM-NoE) under the Contract No. NMP3CT-2004-502243. Appendix A. General framework for the first self-consistent scheme (SC1) The SC1 scheme has been described in the main core of the paper when all inclusions have same shape and same orientation of principal axes. Let us consider an aggregate of N different phases. In the homogenization scheme, each phase is represented by an ellipsoidal inclusion whose shape ratio and orientation is varying from one phase to another. The overall material is subjected to a remote strain rate D. Adopting the self-consistent scheme SC1, one seeks the overall stress tensor R, the stress _ and the macroscopic viscoplastic strain rate tensor Dvp . rate R The interaction law is given by Eq. (18):

   1 1  1 1 d  D ¼ Atg  Ptg : ðs  S Þ þ Ae  Pe : ðr_  R_  Þ

ð56Þ

where Atg , Ae , Ptg , Pe denote fourth order tensors referring to the homogeneous equivalent medium. D , R and R_  are uniform mechanical fields considered to be applied at the remote boundaries associated to each inclusion problem. The viscous behavior of the HEM (having the properties of the overall material) is linearized in the vicinity of Dvp , Eq. (3): vp

vp

D0 ! S0 ¼ Atg : D0 þ So

ð57Þ

S. Mercier, A. Molinari / International Journal of Plasticity 25 (2009) 1024–1048

1045

where Atg is the tangent modulus of the HEM at Dvp and So ¼ S  Atg : Dvp . From Eq. (57), the difference in deviatoric stress tensors S  S is linked to the difference in viscoplastic strain rate tensors Dvp  Dvp  by:

S  S ¼ Atg : ðDvp  Dvp  Þ

ð58Þ

The incremental elastic law for the aggregate gives the following relationship for the stress rate _ : _ and R tensors R 



R_ ¼ Ae : ðD  Dvp Þ R_  ¼ Ae : D  Dvp 



ð59Þ

Note that the rate of the Cauchy stress r_ is given by the local elastic law (15). The last two equations of the SC1 scheme are the consistency equation and the definition of the macroscopic stress:

hdi ¼ D hri ¼ R

ð60Þ

Let’s assume that at time t, the stress tensors r, R and R are given. The strain rate tensors d are defined from the local behavior of the phases (16). The fourth order tensors Ae , Atg are obtained via Eqs. (19) and (20). Since the shape of inclusions are known a priori, the tensors Ptg , P e can be evaluated. The relationships (56) and (58) to (60) provide a complete set of equations for the unknowns: d, D , Dvp , _ _ Dvp  , R and R . When all ellipsoidal inclusions have the same shape ratio and same orientation of the principal axes, Ptg and P e are identical for all inclusions. By volume averaging of Eq. (56), one obtains: vp

   1 1  1 1 D  D ¼ Atg  P tg : ðS  S Þ þ Ae  Pe : ðR_  R_  Þ

ð61Þ

In addition, one has, from relationship (59):

R_  R_  ¼ Ae : ðD  D  Dvp þ Dvp  Þ

ð62Þ

An obvious solution for the set of Eqs. (58), (61) and (62) is: vp D ¼ D S ¼ S R_  ¼ R_ Dvp  ¼ D

ð63Þ

This solution has been adopted for the SC1 scheme developed in the main core of the paper in the case where the ellipsoidal inclusions have same shape ratio and orientation of principal axes.

Appendix B. Numerical resolution for the first self-consistent scheme (SC1) The procedure to determine the evolution of the mechanical fields r, d and dvp in all phases during loading are presented here for the self-consistent scheme (SC1). The material is an aggregate of disordered phases. A representative volume element of the aggregate has been defined and contains N phases or grains. Each phase is represented by an ellipsoidal inclusion. We assume here that the shape ratio and the orientation of the principal axes of inclusions are identical for all phases. Therefore, R ¼ R and D ¼ D. The general case where shape ratios or orientations of the principal axes of inclusions are different could be also treated by using equations of Appendix A. The material is subjected to a remote strain rate D whose time evolution is known. Assume that at time t, the Cauchy stress tensor r has been calculated in all phases. By volume averaging over the representative volume element, the macroscopic stress is obtained: R ¼ hri. Based on the local behavior of each phase defined by Eq. (16), the anelastic strain rate tensors dvp ðtÞ are found and atg ðtÞ can be calculated. As a consequence, Atg ðtÞ is obtained as the solution of the non-linear Eq. (20). The knowledge of the elastic stiffness of each phase defines the elastic modulus of the aggregate, Eq. (19). Knowing Atg ðtÞ and Ae ðtÞ, the related tensors Ptg ðtÞ and P e ðtÞ are calculated. To this point, the total strain rate tensor dðtÞ and the rate of the Cauchy stress tensor r_ ðtÞ in each phase are unknown. From the incremental elastic law (15): r_ ðtÞ ¼ ae : ðdðtÞ  dvp ðtÞÞ, the definition of the rate of the macroscopic stress: R_ ¼ hr_ i and the interaction law (23), the strain rate d is defined by:

1046

S. Mercier, A. Molinari / International Journal of Plasticity 25 (2009) 1024–1048

  D E  1 1  1 1  d  D ¼ Atg  P tg : ðs  SÞ þ Ae  Pe : ae : ðd  dvp Þ  ae : ðd  dvp Þ

ð64Þ

This latter relationship provides N equations for the N unknowns dðtÞ. The set of equations is linear. The local strain rate tensors dðtÞ in each phase are found. Subsequently, the rate of the Cauchy stress tensors r_ ðtÞ are obtained. This enables the update of the stress state in each phase for time t þ Dt, where Dt is the time increment:

rðt þ DtÞ ¼ rðtÞ þ r_ ðtÞDt

ð65Þ

Thus, the evolution of the mechanical fields can be determined during the whole loading. This algorithm is straightforward. Compared to technique which needs Laplace transform, the computing time is much reduced. Appendix C. Numerical resolution for the second self-consistent scheme (SC2) The procedure adopted to update, from time t to time t þ Dt, the mechanical fields in each phase is explained in the following for the second self-consistent scheme SC2. We assumed here that all inclusions have the same shape ratio and same principal axes. The general case could be treated as well by using the equations of Section 2.4.1. We first consider a volume of heterogeneous medium subjected to uniform strain rate D. At time t, the Cauchy stress tensors in all phases rðtÞ and the macroscopic stress tensor RðtÞ are supposed to be known. From the local behavior of the different phases, dvp ðtÞ can be calculated in terms of rðtÞ, and the tensors for local elastic and viscoplastic tangent moduli can be obtained. Using relationships (19) and (20), the fourth order tensors Ae (t) and Atg ðtÞ are evaluated. By volume averaging over all phases and using hdi ¼ D and hri ¼ R, one has from Eq. (31):

   1 1  1 1   1 D  D ¼ Atg  P tg : ðS  S Þ þ Ae  Pe : R_  R_ 

a

ð66Þ

From Eqs. 29, 32 and 66 and noticing that:



  1 1  1 1 Atg  Ptg : ðS  S Þ ¼ Atg  Ptg : ðR  R Þ;

the strain rate difference D  D is linked to the stress difference R  R by:

"

#   1 1 e e e I A  P : A : ðD  D Þ ¼

" #   1 1  1 1  1 Atg  Ptg  m Ae  P e : Ae : Atg

: ð R  R Þ

ð67Þ

Initially, the material is undeformed without any internal stresses. Thus, at t ¼ 0, r ¼ R ¼ R ¼ 0. Let us consider that at time t, the stress tensors R and R are related by R ¼ R. From Eq. (67), one _ ¼R _ . Therefore, at time t þ Dt, one also obtains that D ¼ D. With Eqs. (29) and (32), it follows: R  has: R ¼ R. Finally, at any time, when inclusions have same ellipsoidal shape ratio with same orientation of the principal axes, one gets:

D ¼ D;

R ¼ R;

R_  ¼ R_

ð68Þ

In that case, the relationship (31) is rewritten as:

h  i1 d ¼ ae  a Ae  Pe1 " ( )#    1     1 1 e vp e tg e e1 tg tg : Dþ A  P : a : d þ aA : D  m A :R a A P : ðs  SÞ

ð69Þ

1047

S. Mercier, A. Molinari / International Journal of Plasticity 25 (2009) 1024–1048

Note that for linear behavior, applying the volume averaging operator on Eq. (69) leads to Eq. (44). At this stage, only the strain rate tensors dðtÞ and aðtÞ have to be determined. The parameter aðtÞ is calculated so as to satisfy the consistency equation hdðtÞi ¼ DðtÞ. For complex loading, the consistency condition can not be always rigorously fulfilled. Thus, aðtÞ is evaluated in order to minimize the relpffiffiffiffiffiffiffiffiffi kDhdik ative error kDk where the norm kxk ¼ x : x can be used for instance. Note that in uniaxial tension, the consistency condition is fulfilled exactly. From the incremental Hooke’s law, r_ ðtÞ can be calculated and the stress tensors can be updated by using an Euler forward scheme: rðt þ DtÞ ¼ rðtÞ þ r_ ðtÞDt. Using Eq. (29), the rate of macroscopic _ ðtÞ can be evaluated and the stress tensor Rðt þ DtÞ updated. The quantities rðt þ DtÞ stress tensor R and Rðt þ DtÞ are now defined at time t þ Dt. The iterative process can be pursued for the next time increment. From volume averaging performed on Eq. (31) and using relationship (68), one shows that the difference between the volume average of the local stresses hri and the macroscopic stress R satisfies a first order differential equation:

   1 1  1 1 1 Atg  P tg : ðhri  RÞ þ Ae  Pe : ðhr_ i  R_ Þ ¼ 0

ð70Þ

a

The consistency equation hdi ¼ D has been used to derive Eq. (70). Since initially, proved from Eq. (70) that by construction, one has hri ¼ R.

r ¼ R ¼ 0, it is

Appendix D. Numerical resolution for the Mori-Tanaka scheme The update procedure for the Mori-Tanaka scheme is quite similar to the algorithm presented in Appendix C. At time t, the Cauchy stress tensors in the matrix Rm ðtÞ and in other phases rðtÞ are supposed to be known. From the local behavior of phases, the viscoplastic strain rate tensors dvp ðtÞ and vp Dvp and Avp m ðtÞ are obtained, and the viscoplastic tangent stiffness tensors a m can be calculated. The elastic stiffness tensors ae and Aem are given. From the interaction law (36) and from the incremental elastic law for all phases, the total strain rate in each phase is obtained at time t by:

   1  1 d ¼ ae  Aem  Pem " ( )#    1   1 1   e vp e vp e e tg tg : a : d þ Am : Dm  Dm  Am  Pm : ðs  Sm Þ : D m þ Am  P m

ð71Þ

From the consistency equation hdðtÞi ¼ DðtÞ, the strain rate in the matrix Dm ðtÞ is found and henceforth, the strain rate in all other phases dðtÞ. Using the incremental Hooke’s law and an Euler forward scheme, the quantities rðt þ DtÞ, Rm ðt þ DtÞ are defined at time t þ Dt. The process can be continued for the next time increment. References Abdul-Latif, A., 2004. Pertinence of the grains aggregate type on the self-consistent model responses. Int. J. Solids Struct. 41, 305–322. Asaro, R., Needleman, A., 1985. Overview. Texture development and strain-hardening in rate dependent polycrystals. Acta Metall. 33, 923–953. Barai, P., Weng, G., 2007. Mechanics of creep resistance in nanocrystalline solids. Acta Mech. 24, 327–348. Berbenni, S., Favier, V., Berveiller, M., 2007. Impact of the grain size distribution on the yield stress of heterogeneous materials. Int. J. Plast. 23, 114–142. Berveiller, M., Zaoui, A., 1979. An extension of the self-consistent scheme to plastically flowing polycrystals. J. Mech. Phys. Solids 26, 325–344. Beurthey, S., 1997. Modélisations du comportement mécanique d’alliages de polymères. Ph.D. Thesis, Ecole Polytechnique, Palaiseau, France. Eshelby, J.D., 1957. The determination of the elastic field of an ellipsoidal inclusion and related problems. Proc. Roy. Soc. Lond. A 241, 376–396. Hashin, Z., 1969. The inelastic inclusion problem. Int. J. Eng. Sci. 7, 11–36. Hill, R., 1965. Continuum micro-mechanisms of elastoplastic polycrystals. J. Mech. Phys. Solids 13, 89–101.

1048

S. Mercier, A. Molinari / International Journal of Plasticity 25 (2009) 1024–1048

Hutchinson, J.W., 1976. Bounds and self-consistent estimates for creep of polycrystalline materials. Proc. Roy. Soc. Lond. A 348, 101–127. Kouddane, R., Molinari, A., Canova, G.R., 1993. Self-consistent modelling of heterogeneous viscoelastic and elastoviscoplastic materials. In: Teodosiu, C., Raphanel, J., Sidoroff, F. (Eds.), Large Plastic Deformation. Balkema, Rotterdam, pp. 129–141. Kreher, W., 1990. Residual stresses and stored elastic energy of composites and polycrystals. J. Mech. Phys. Solids 38, 115–128. Kreher, W., Molinari, A., 1993. Residual stresses in polycrystals as influenced by grain shape and texture. J. Mech. Phys. Solids 41, 1955–1977. }ner, E., 1961. Zur plastischen verformung des vielkristalls. Acta Metall. 9, 155–161. Kro Lebensohn, R., Tomé, C., 1993. A self-consistent anisotropic approach for the simulation of plastic deformation and texture development of polycrystals: application to zirconium alloys. Acta Metall. Mater. 41, 2611–2624. Li, J., Weng, G., 1994. Strain-rate sensitivity relaxation behavior and complex moduli of a class of isotropic viscoelastic composites. ASME J. Eng. Mater. Technol., 494–504. Li, J., Weng, G., 1997. A secant-viscosity approach to the time-dependent creep of an elastic–viscoplastic composite. J. Mech. Phys. Solids 45, 1069–1083. Li, J., Weng, G., 1998. A unified approach from elasticity to viscoelasticity to viscoplasticity of particle-reinforced solids. Int. J. Plast. 14, 193–208. Masson, R., Zaoui, A., 1999. Self-consistent estimates for the rate-dependent elastoplastic behavior of polycrystalline materials. J. Mech. Phys. Solids 47, 1543–1568. Mercier, S., Jacques, N., Molinari, A., 2005. Validation of an interaction law for the Eshelby inclusion problem in elastoviscoplasticity. Int. J. Solids Struct. 42, 1923–1941. Mercier, S., Molinari, A., Estrin, Y., 2007. Grain size dependence of strength of nanocrystalline materials as exemplified by copper: an elastic–viscoplastic modelling approach. J. Mater. Sci. 42, 1455–1465. Molinari, A., 1999. Extensions of the self-consistent tangent model. Model. Simul. Mater. Sci. Eng. 7, 683–697. Molinari, A., 2002. Averaging models for heterogeneous viscoplastic and elastic–viscoplastic materials. J. Eng. Mater. Technol. 124, 62–70. Molinari, A., Ahzi, S., Kouddane, R., 1997. On the self-consistent modelling of elastic–plastic behavior of polycrystals. Mech. Mater. 26, 43–62. Molinari, A., Canova, G.R., Ahzi, S., 1987. A self consistent approach of the large deformation polycrystal viscoplasticity. Acta Metall. 35, 2983–2994. Molinari, A., El-Houdaigui, F., Tóth, L.S., 2004. Validation of the tangent formulation for the solution of the non-linear Eshelby inclusion problem. Int. J. Plast. 20, 291–307. Molinari, A., Mercier, S., 2004. Homogeneization of viscoplastic materials. In: Ahzi, S., Cherkaoui, M., Khalel, M., Zbib, H., Zikry, M., LaMatina, B. (Eds.), Multiscale Modeling and Characterization of Elastic–Inelastic Behavior of Engineering Materials. Kluwer Academics Publishers. Ortiz, M., Molinari, A., 1988. Microstructural thermal stresses in ceramic materials. J. Mech. Phys. Solids 36, 385–400. Paquin, A., 1998. Modélisation micromécanique du comportement élastoviscoplastique des matériaux hétérogènes. Ph.D. Thesis, University of Metz, Metz, France. Paquin, A., Sabar, H., Berveiller, M., 1999. Integral formulation and self-consistent modelling of elasto-viscoplastic behavior of heterogeneous materials. Arch. Appl. Mech. 69, 14–35. Perzyna, P., 1986. Internal state variable description of dynamic fracture of ductile solids. Int. J. Solids Struct. 22, 797–818. Pierard, O., Doghri, I., 2006. An enhanced affine formulation and the corresponding numerical algorithms for the mean-field homogenization of elasto-viscoplastic composites. Int. J. Plast. 22, 131–157. Pierard, O., Llorca, J., Segurado, J., Doghri, I., 2007. Micromechanics of particle-reinforced elasto-viscoplastic composites: Finite element simulations versus affine homogenization. Int. J. Plast. 23, 1041–1060. Ponte-Castañeda, P., 2002. Second-order homogenization estimates for nonlinear composites incorporating field fluctuations: Itheory. J. Mech. Phys. Solids 50, 737–757. Ponte-Castañeda, P., Suquet, P., 1998. Non linear composites. Adv. Appl. Mech. 34, 171–302. Rougier, Y., Stoltz, C., Zaoui, A., 1994. Self-consistent modelling of elastic–viscoplastic polycrystals. C. R. Acad. Sci. Paris, Série II 318, 145–151. Sabar, H., Berveiller, M., Favier, V., Berbenni, S., 2002. A new class of micro–macro models for elastic–viscoplastic heterogeneous materials. Int. J. Solids Struct. 39, 3257–3276. Schapery, R., 1962. Approximate methods of transform inversion in viscoelastic stress analysis. In: Proceedings of 4th US National Congress on Applied Mechanics, vol. 2, pp. 1075–1085. Suquet, P., 1985. Elements of homogeneization for inelastic solid mechanics. In: Sanchez-Palencia, E., Zaoui, A. (Eds.), Homogeneization Techniques for Composite Media. Springer-Verlag, Berlin, pp. 193–278. Tòth, L., Molinari, A., Zouhal, N., 2000. Cyclic plasticity phenomena as predicted by polycrystal plasticity. Mech. Mater. 32, 99– 113. Turner, P.A., Tomé, C.N., Woo, C.H., 1994. Self-consistent modelling visco-elastic polycrystals: an approximate scheme. Philos. Mag. A 70, 689–711. Weng, G., 1993. A self-consistent relation for the time-dependent creep of polycrystals. Int. J. Plast. 9, 181–198. Zouhal, N., Molinari, A., Toth, L., 1996. Elastic–plastic effects as predicted by the Taylor–Lin model of polycrystal elastoviscoplasticity. Int. J. Plast. 12, 343–360.