# Self-consistent estimates for the Rate-Dependent elastoplastic behaviour of polycrystalline materials

Renaud Masson, André Zaoui

► To cite this version:

Renaud Masson, André Zaoui. Self-consistent estimates for the Rate-Dependent elastoplastic behaviour of polycrystalline materials. Journal of the Mechanics and Physics of Solids, 1999, 47, pp.1543-1568. 10.1016/S0022-5096(98)00106-9. hal-00111618

HAL Id: hal-00111618

https://hal.science/hal-00111618v1

Submitted on 12 Nov 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Self-consistent estimates for the rate-dependent elastoplastic behaviour of polycrystalline materials

R. Masson, A. Zaoui

Laboratoire de Mécanique des Solides, CNRS, Ecole Polytechnique, 80017 Palaiseau Cedex, France

## Abstract

This paper aims at proving that, contrary to previous contributions to the subject, Hills conception of the nonlinear self-consistent scheme, which has been applied in the past to elastoplasticity and to viscoplasticity, can still be adopted with success for elastoviscoplasticity. After a qualitative presentation of the main arguments for this statement, a new Hill-type approach is proposed for rate-dependent elastoplastic heterogeneous materials. The associated linearization procedure relies on an affine formulation instead of Hills incremental one and on the use of the correspondence principle to solve the concentration problem; this problem is proved to reduce to a linear viscoelastic one with eigenstrain, i.e. to a linear thermoviscoelastic problem. The full set of equations is reported for the case of the self-consistent scheme and illustrative applications are given for polycrystals: they are shown to be, as expected, always softer than Kröner-type predictions and to take better into account the viscoelastic coupling and the associated long range memory effect. In conclusion, the connection and differences between the present approach and other ones already proposed for viscoplastic materials is emphasized and the limits of Hills conception itself are acknowledged and discussed.

*Keywords:* Creep; Elastic-viscoplastic material; Polycrystalline material; Inhomogeneous material

## 1. Introduction

The self-consistent scheme has been intensively used for forty years in view of the prediction and the interpretation of the mechanical properties of heterogeneous materials, especially polycrystals. Roughly speaking (Zaoui, 1997), this scheme is

<footnote>* Corresponding author</footnote>
based on a specific method of estimating the average mechanical interactions between every constituent phase and all the other ones: these interactions are approximated by the ones which are sustained by an ellipsoidal inhomogeneity, representative of the considered constituent, embedded in a homogeneous infinite matrix, made of the homogeneous equivalent medium (H.E.M.) which is searched for, subjected to the macroscopic loading. This medium is determined by a self-consistency condition which expresses the fact that the H.E.M.'s overall response results from the weighted individual contributions of all the phases when they are considered alternatively as ellipsoidal inhomogeneities.

Whereas this scheme is now widely accepted as a well-defined and efficient model for elastic (Hershey, 1954; Kroner, 1958) or linear viscoelastic (Laws and McLaughlin, 1978) disordered materials such as polycrystals, its extensions to nonlinear behaviour are still open to discussion. During the last decades, two main approaches have been proposed in this field. The first one, initially developed by Kroner (1961) for elastoplastic polycrystals as a straightforward application of Eshelby's solution of the inclusion problem (Eshelby, 1957), was extended later to rate-dependent elastoplastic materials by Weng (1981) for infinitesimal strains and by Nemat-Nasser and Obata (1986) and by Harren (1991) for finite strains. The second one was proposed by Hill (1965) as an alternative model to that of Kroner and then applied and extended to elastoplastic polycrystals (Hutchinson, 1970; Berveiller and Zaoui, 1979; Iwakuma and Nemat-Nasser, 1984; Lipinski et al., 1990) and developed for viscoplastic materials (Hutchinson, 1976; Molinari et al., 1987; Lebensohn and Tome, 1993). It is based on the linearization of the constitutive equations and on the use of Eshelby's solution of the elastic inhomogeneity problem (Eshelby, 1957).

Whereas Hill's version of the self-consistent scheme rapidly gained acceptance for rate-independent elastoplasticity, for definite reasons which are recalled in the following, no attempt has been made up to now to extend this version to rate-dependent elastoplasticity. On the contrary, Kroner's model has been either claimed as adapted better (Weng, 1981) to this case or implicitly considered as such (Nemat-Nasser and Obata, 1986; Harren, 1991). This paper aims at proving the superiority of Hill's conception in this case too.

First, the matter of Hill vs Kroner's controversy for elastoplastic polycrystals is briefly recalled in order to prove that the same arguments which lead to acknowledging the better adequacy of Hill's conception still hold for rate-dependent elastoplasticity (Section 2). Second, a new Hill-type approach relying on an 'affine' formulation is proposed for this latter case (Section 3) and some applications are developed as an illustration and compared to Kroner's model-type predictions (Section 4). It is shown that, as in the case of rate-independent plasticity, Hill's conception leads to softer predictions than Kroner's one. In conclusion (Section 5), the limits and deficiencies of Hill's conception itself are briefly discussed and further developments are suggested.

## 2. Kroner's vs Hill's versions of the nonlinear self-consistent model

Hill vs Kroner's controversy for the self-consistent modelling of elastoplastic polycrystals has been already exposed and discussed for a long time (e.g. Berveiller and
Zaoui, 1979). First we briefly report the matter of this debate, in order to extend the discussion to the new framework of rate-dependent elastoplasticity. Kroner's approach was based on the assimilation of the plastic strain to a 'stress-free strain' (or 'transformation strain', 'eigenstrain'...) as defined by Eshelby in the context of elasticity and on the associated elastic treatment of the ellipsoidal inclusion problem: at any given stage of the plastic flow, while the polycrystal is sustaining the macroscopic plastic strain $E^p$ under the overall stress $\Sigma$, every set of grains ($r$) with the same lattice orientation, average shape and average plastic strain $\varepsilon_r^p$ is supposed to undergo the same average back-stress ($\sigma_r - \Sigma$) as the one sustained by an inclusion (with the uniform eigenstrain $\varepsilon_r^p$) embedded in an infinite matrix (with the uniform eigenstrain $E^p$). For isotropic elasticity, isochoric plastic strains and spherical shapes, this treatment leads to the well-known Kroner's concentration (or interaction) equation:

$$ \underline{\underline{\sigma}}_r - \underline{\underline{\Sigma}} = 2\mu(1 - \beta)(E^p - \underline{\underline{\varepsilon}}_r^p) \quad (1) $$

with $\beta$ an elastic coefficient ranging 0.4–0.6.

Such a treatment of the intergranular interactions is clearly open to the criticisms of being a purely elastic (instead of an elastoplastic) one, according to which the plastic strains in the inclusion and in the matrix considered when solving the concentration problem are not disturbed by these mechanical elastic interactions. Consequently, this leads to a strong overestimation of these interactions and of the resulting overall stiffness of the polycrystal. As a matter of fact, the elastic multiplier $2\mu(1 - \beta)$ in (1), which is of the order of $\mu$, is so high with respect to the yield stress level that it leads to almost null plastic strain deviations ($\varepsilon_r^p - E^p$), which practically makes Kroner's model coincide with that of Taylor (Taylor, 1938), based on the assumption of uniform plastic strains.

This is the reason why Hill (1965) proposed an alternative treatment relying on the incremental linearization of the elastoplastic constitutive equations: at a given stage, the ($r$) phase has the instantaneous (multibranched) elastoplastic moduli $L_r^{ep}$ defined for the average strain $\varepsilon_r$, whereas the H.E.M. has the overall elastoplastic moduli $L^{ep}$ so that the concentration equation derives from the solution of the corresponding Eshelby's inhomogeneity problem for the ellipsoid $\Omega_r$ with the moduli $L_r^{ep}$ in an infinite matrix with the moduli $L^{ep}$. In simple cases it can be put in the form

$$ \dot{\underline{\underline{\sigma}}}_r - \underline{\underline{\Sigma}} = -L_r^*: (\underline{\underline{\varepsilon}}_r - \underline{\underline{E}}), \quad (2) $$

where Hill's constraint tensor $L_r^*$ depends on $L^{ep}$ and on the shape and orientation of $\Omega_r$.

It has been significantly checked (Hutchinson, 1970; Berveiller and Zaoui, 1979; Lipinski et al., 1990) that, in agreement with most experimental data, Hill's incremental model yields much softer predictions than those of Kroner and Taylor and can take into better account the specific hardening properties of the considered crystals. Consequently, Hill's treatment was accepted as better than Kroner's a long time ago, as far as rate-independent elastoplasticity was considered. When the time came to extend the self-consistent scheme to elastoplastic finite strains (Iwakuma and Nemat-Nasser, 1984; Lipinski et al., 1990), Hill's version was naturally chosen as the adequate
basis for this development; moreover, when Hutchinson wanted to derive self-consistent estimates for creep of polycrystalline materials (Hutchinson, 1976) in the framework of viscoplasticity (without elasticity), he naturally founded his development on Hill's version too, by replacing strains by strain rates in the basic inhomogeneity problem (before deriving a 'secant' formulation for the special case of power-law creep).

The question of rate-dependence in the context of the self-consistent scheme was first addressed by Laws and McLaughlin (1978) for linear viscoelasticity: they proposed to use the principle of correspondence and the Laplace transform technique in order to convert the viscoelastic problem into a symbolical elastic one and to derive the viscoelastic response by final inversion of the transformed variables. But such an elegant solution of the problem obviously fails for nonlinear viscoelasticity or elastoviscoplasticity. On the other hand, Hill's linearization method does not work any more, at least in its simplest form, due to the simultaneous occurrence in the constitutive equations of time derivatives of different orders of stress and strain, so that no simple 'tangent' moduli or compliances can be derived straightforwardly.

The existence of such an apparent deadlock on the way to extend Hill's method to rate-dependent elastoplasticity may have been the reason why Kroner's approach was then called for help. Weng (1981) even found some physical reasons to do so: his analysis of the 'time-dependent creep problems' concluded that "... Kroner's model then becomes entirely consistent... The crux of the matter is that creep, unlike plastic deformation, is a truly 'stress-free' process in the sense of Eshelby, because the creep strain rate, at any generic state, depends only on the current stress and deformation history and is independent of the stress rate. The latter independence obviously does not hold in plasticity, and this leads to the presence of tangent moduli in Hill's constraint tensor. This subtle point appears to have not been realized...".

Actually, whether the creep strain rate is dependent or not on the stress rate does not matter: the actual 'crux of the matter' is that, 'at any generic state' and at any point $\underline{x}$ of the heterogeneous material, the anelastic strain and its rate depend on the local mechanical state and are not known in advance. In other words, the simple decomposition of the total local strain into its elastic and viscoplastic parts, say

$$ \underline{\varepsilon}(\underline{x}, t) = \underline{\varepsilon}^{\text{el}}(\underline{x}, t) + \underline{\varepsilon}^{\text{vp}}(\underline{x}, t) \qquad (3) $$

with Hooke's equations ruling the elastic behaviour

$$ \underline{\varepsilon}^{\text{el}}(\underline{x}, t) = \underline{s}(\underline{x}):\underline{\sigma}(\underline{x}, t) \qquad (4) $$

through the elastic compliance $\underline{s}(\underline{x})$ cannot be considered as yielding constitutive equations for elastoviscoplasticity as long as the stress-dependence of $\underline{\varepsilon}^{\text{vp}}(\underline{x}, t)$ has not been specified. Recall here that Eshelby's treatment of the inclusion problem relies on the fact that the 'stress-free strain' is given a priori and does not depend on the mechanical loading, so that the foregoing relations are (elastic) constitutive equations and can be used as such in order to solve the concentration problem. If rate-dependent as well as rate-independent plasticity is considered and if the same equations are used as constitutive equations when solving the concentration problem, this means that
the considered constitutive behaviour is a purely elastic one (with eigenstrains): consequently, the resulting concentration equations will express elastic intergranular interactions instead of elasto(visco)plastic ones, as it was already discussed when criticizing Kroner's original model.

More recently, Weng's formulation was extended to finite strains either within an Eulerian (Nemat-Nasser and Obata, 1986) or a Lagrangian (Harren, 1991) approach. The corresponding analyses did not resume explicitly Weng's arguments about the physical nature of the creep strain rate but they used the classical Green function technique in a way which indeed presupposed Weng's arguments to be valid. If, for the sake of simplicity, we get rid of the technical complication of the finite strain formalism which is not essential for our present purpose and if we convert their treatment into the framework of infinitesimal strain analysis, this way may be summed up as follows:

—consider the decomposition of the total strain rate into its elastic and viscoplastic parts and Hooke's equations as constitutive equations, both at the micro and the macroscale, namely

$$ \begin{aligned} \dot{\underline{\underline{\sigma}}}_r &= \mathbf{c}_r : (\underline{\underline{\varepsilon}}_r - \underline{\underline{\varepsilon}}_r^{\text{vp}}) \\ \underline{\underline{\Sigma}} &= \mathbf{C}^{\text{SC}} : (\underline{\underline{E}} - \underline{\underline{E}}^{\text{vp}}) \end{aligned} \qquad (5) $$

with 'SC' standing for 'self-consistent' and with $\mathbf{c}_r$ and $\mathbf{C}^{\text{SC}}$ the local and overall elastic moduli;

—use the Green technique in order to solve the elastic problem of an ellipsoidal inhomogeneous inclusion $\Omega_r(\mathbf{c}_r, \underline{\underline{\varepsilon}}_r^{\text{vp}})$ embedded in the infinite homogeneous matrix $(\mathbf{C}^{\text{SC}}, \underline{\underline{E}}^{\text{vp}})$ through the strain-Green operator of the (elastic) effective medium $\Gamma^{\text{SC}}$, with $\underline{\underline{\varepsilon}}_r^{\text{vp}}$ and $\underline{\underline{E}}^{\text{vp}}$ uniform. This resolution can make use of the well-known uniformity of the (hereafter mentioned as) 'Green tensor'

$$ \mathbf{P}_r^{\text{SC}}(\underline{x}) = \int_{\Omega_r} \Gamma^{\text{SC}}(\underline{x} - \underline{x}') \, d\Omega' \qquad (6) $$

when $\underline{x} \in \Omega_r$;

—combine the resulting concentration relation with the constitutive equations and the average conditions in order to predict the overall behaviour.

The corresponding equations are a special case of the general equations for the homogenization of a heterogeneous linear thermoelastic medium, i.e. a linear elastic medium with eigenstrains; if $\mathbf{A}_r$ and $\mathbf{B}_r$ are the average strain and stress concentration tensors for the phase ($r$), respectively, in the absence of eigenstrains and if $\underline{\underline{\varepsilon}}_r^T$ and $\underline{\underline{E}}^T$ are the local and overall eigenstrains, we have, as classically:

$$ \mathbf{C}^{\text{eff}} = \langle \mathbf{c} : \mathbf{A} \rangle = \langle {}^t \mathbf{A} : \mathbf{c} : \mathbf{A} \rangle $$

$$ \mathbf{S}^{\text{eff}} = (\mathbf{C}^{\text{eff}})^{-1} = \langle \mathbf{s} : \mathbf{B} \rangle = \langle {}^t \mathbf{B} : \mathbf{s} : \mathbf{B} \rangle $$
$$ \mathbf{A}_r = \mathbf{s}_r : \mathbf{B}_r : \langle \mathbf{s} : \mathbf{B} \rangle^{-1} $$

$$ E_{\underline{\underline{\underline{\mathbf{\epsilon}}}}}^{\mathrm{T}} = \langle {}^{\mathrm{T}}\mathbf{B} : \underline{\underline{\underline{\mathbf{\epsilon}}}}^{\mathrm{T}} \rangle = \mathbf{S}^{\mathrm{eff}} : \langle {}^{\mathrm{T}}\mathbf{A} : \mathbf{c} : \underline{\underline{\underline{\mathbf{\epsilon}}}}^{\mathrm{T}} \rangle, \qquad (7) $$

where ${}^{\mathrm{T}}\mathbf{A}$ stands for the transposed tensor of $\mathbf{A}$ and $\langle \cdot \rangle$ denotes a spatial average; in addition, a (model-dependent) concentration equation and a specific expression for $\mathbf{A}_r$ or $\mathbf{B}_r$ are needed.

When applied to the self-consistent scheme and to the special case which is under consideration here, the resulting procedure is nothing but the self-consistent elastic treatment for an elastic inhomogeneous body with (known) eigenstrain rates. It is easy to show that additional simplifications (uniform elasticity, spherical representative inclusions and isotropy) lead to Kroner's concentration eqn (1), written for rates. An illustration of the fact that such a modelling of rate-dependent elastoplastic heterogeneous materials reduces to the transcription (or the extension to finite strains) of Kroner's model can be found in the comparisons which have been performed by Harren (1991) between its own predictions and those derived from the Taylor model for various responses and texture evolutions corresponding to problems of metal forming: it can be checked that, as expected from what precedes, these predictions almost coincide....

At this point, we can conclude that, for rate-dependent elastoplasticity, Kroner's approach is as inadequate as it was for rate-independent elastoplasticity and that Hill's conception would be more appropriate in this case too. It remains to make efforts in order to accommodate Hill's conception to the specific difficulties associated with the rate-dependence and especially with the critical viscoelastic coupling. In addition, while adopting Hill's general point of view, we can feel free to choose another formulation than Hill's incremental one. This is done in what follows by combining a new 'affine' formulation and the use of the correspondence principle (and the associated Laplace transform technique) in this new linearized quasi-'ther-moviscoelastic' context.

### 3. A new Hill-type formulation for rate-dependent elastoplasticity

If we try to extract the essential features of Hill's approach (Hill, 1965) to the self-consistent modelling of nonlinear heterogeneous materials, we can stress two basic ideas: the first one is the reduction of the nonlinear homogenization procedure to a (quasi-)linear one through some adequate linearization of the constitutive equations; the second one refers to the derivation of self-consistent concentration equations from the resolution of Eshelby-type problems where both the ellipsoidal inhomogeneity and the infinite matrix are given uniform linearized properties. Within this framework, several variants can exist and we can find in the literature several Hill-type formulations, for elastoplasticity as well as for viscoplasticity, which mainly differ from each other in the chosen linearization method. So, we have first to decide which formulation to adopt.
### 3.1. Which formulation?

In the case of elastoplasticity, as briefly mentioned above [eqn (2)], the original Hill method consisted in an incremental procedure based on the definition of instantaneous multibranched elastoplastic moduli $I_r^{ep}$ and $L^{ep}$ relating the local and global stress and strain rates, respectively. The choice of an incremental formulation was closely related to the multibranched nature of the elastoplastic constitutive equations and to the irreversibility of the plastic strain. According to this incremental formulation (Hutchinson, 1970), the overall elastoplastic response to a given loading path is derived through a step by step procedure: at each step, the initial mechanical state of every phase is known from the resolution of the linearized self-consistent problem associated to the previous step; the local elastoplastic moduli $I_r^{ep}$ are then updated as is the global one $L^{ep}$ too from the self-consistent equation. The mechanical state of every phase can be updated from the concentration equations and the procedure may be resumed for the next step.

An alternative approximate treatment which can still be considered as related to Hill's approach was proposed later (Berveiller and Zaoui, 1979) for monotonic proportional loading paths: based on the deformation theory of plasticity, it makes use of isotropic secant elastoplastic moduli and does not differ in nature from a nonlinear elastic treatment. At any stage, the overall secant moduli and the associated local and global mechanical variables can be derived independently of their values at any other stage through the (nonlinear) determination of a self-consistent ensemble of local and global quantities. A secant formulation has also been proposed in the case of viscoplasticity: starting from a strict transcription of Hill's incremental formulation, now relating the stress rate and the second-order time derivative of the strain, Hutchinson (1976) showed that for power-law creep this formulation can rigorously be integrated into a 'total' one which makes use of anisotropic secant creep compliances.

Apart from various extensions to finite strains which can be omitted here, a new Hill-type formulation was proposed by Molinari et al. (1987) for viscoplasticity and developed later by Lebensohn and Tome (1993). This formulation is neither an incremental nor a secant one: it makes use of 'tangent' moduli as well as of initial stresses for the linearization of the local and overall constitutive equations. Roughly speaking, it proposes to assimilate at any stage the nonlinear curve correlating the viscoplastic strain rate and the stress to its tangent at this stage. This is done not only for the constituent phases but also for the overall behaviour. Through the restriction to power-law creeping phases with the same exponent, the authors consider as granted that the overall response obeys a similar power-law and that the overall moduli are 'tangent' to the overall response. That is why this formulation can be referred to as a 'tangent formulation'.

Starting from a similar linearization method at the local level, adapted to the case of rate-dependent elastoplastic materials, Rougier et al. (1994) proposed a different homogenization procedure of the thermo(visco)elastic type, referring to linear viscoelasticity with eigenstrains. We adopt the same general framework in the following, but we modify different aspects of this formulation in order to transform it into a
consistent ‘affine’ one (see Masson, 1998 or Zaoui and Masson, 1998) in the specific context of rate-dependent elastoplasticity.

## 3.2. The linearization procedure

We consider a representative volume element $\Omega$ of a multiphase material. At each point $\underline{x}$ and at any time $t$, the total strain rate is the sum of the elastic and the inelastic (either viscous or viscoplastic) strain rates according to (3). The inelastic strain generally depends on internal parameters which model irreversible or viscous processes such as hardening effects. Furthermore, we assume time effects to dominate (no rate-independent plasticity) so that the inelastic strain is derived from a viscoplastic (or nonlinear viscoelastic) potential $w(\underline{x}, \underline{\sigma}, \underline{\alpha})$:

$$ \dot{\underline{\varepsilon}}^{\text{in}}(\underline{x}, t) = \frac{\partial w}{\partial \underline{\sigma}}(\underline{x}, \underline{\sigma}, \underline{\alpha}) = \underline{g}(\underline{x}, \underline{\sigma}, \underline{\alpha}), \qquad (8) $$

where the internal parameters $\underline{\alpha}$, which are supposed here, for the sake of simplicity, to be vectorial variables, obey evolution equations of the form

$$ \dot{\underline{\alpha}}(\underline{x}, t) = \underline{h}(\underline{x}, \underline{\sigma}, \underline{\alpha}). \qquad (9) $$

Both functions $\underline{g}$ and $\underline{h}$ are nonlinear but continuous in $\underline{\sigma}$ and $\underline{\alpha}$.

The mechanical macroscopic loading is supposed to be null before $t = 0$ and known at any time $t \ge 0$. We want to predict the overall response at some time $t = \tau$. Thus we have to derive linearized local constitutive equations referring to the heterogeneous material at time $t = \tau$ but, due to the viscoelastic coupling, these constitutive equations have to be expressed for any $t$; in addition, for the same reason, we have to be able, at $t = \tau$, to consider as known the whole local stress and strain histories at any point $\underline{x}$ up to $t = \tau$.

Before dealing with this specific difficulty (see Section 3.3), we have to identify the mechanical nature of the local linearized constitutive equations in order to predict the form of the global ones. This will be first discussed, for a better understanding, in the simpler case where the flow function $\underline{g}$ does not depend on internal variables $\underline{\alpha}$. Since the variable $\underline{x}$ remains fixed in this section, it is omitted in the following.

### 3.2.1. Simple case

In this case, the local nonlinear constitutive equations read

$$ \dot{\underline{\varepsilon}}^{\text{in}}(t) = \underline{g}(\underline{\sigma}(t)). \qquad (10) $$

According to our affine approach, the linearized constitutive equations read for $t \ge \tau$:

$$ \begin{aligned} \dot{\underline{\varepsilon}}(t) &= \mathbf{s}:\underline{\sigma}(t) + \mathbf{m}(\tau):\underline{\sigma}(t) + \underline{\varepsilon}^0(\tau, t) \quad (t \ge \tau), \\ \mathbf{m}(\tau) &= \frac{\mathrm{d}\underline{g}}{\mathrm{d}\underline{\sigma}}(\underline{\sigma}(\tau)), \\ \dot{\underline{\varepsilon}}^0(\tau, t) &= \underline{g}(\underline{\sigma}(\tau)) - \mathbf{m}(\tau):\underline{\sigma}(\tau), \end{aligned} \qquad (11) $$
where the fourth-order tangent creep compliance $\mathbf{m}(\tau)$ and the strain rate $\dot{\underline{\varepsilon}}^0(\tau, t)$ express the stress-dependence of the inelastic strain in a linearized way.

As explained above, we need, at $t = \tau$, constitutive equations for any $t$: so an expression similar to (11) must be defined for $0 \le t \le \tau$ too. Instead of a linear relation, it is then better to use the potentially known nonlinear expression $\underline{\varepsilon}^{\text{in}}(t) = \underline{g}(\underline{\sigma}(t))$ for $t \le \tau$. Consequently, denoting by $H(t)$ the unit step function at $t = 0$, the complete set of constitutive equations to be considered in the state corresponding to $\tau$ reads, for any $t$:

$$ 
\begin{aligned}
\dot{\underline{\varepsilon}}(t) &= \mathbf{s}: \dot{\underline{\sigma}}(t) + \mathbf{m}(\tau): \underline{\sigma}(t) + \underline{\varepsilon}^0(\tau, t), \\
\mathbf{m}(\tau) &= \frac{d\underline{g}}{d\underline{\sigma}}(\underline{\sigma}(\tau)), \\
\dot{\underline{\varepsilon}}^0(\tau, t) &= \underline{g}(\underline{\sigma}(\tau)) - \mathbf{m}(\tau): \underline{\sigma}(\tau) \\
&\quad + [\underline{g}(\underline{\sigma}(t)) - \underline{g}(\underline{\sigma}(\tau)) - \mathbf{m}(\tau):(\underline{\sigma}(t) - \underline{\sigma}(\tau))][1 - H(t - \tau)].
\end{aligned} 
\quad (12) $$

The crucial point of the proposed method lies in the nature of these constitutive equations and especially in the nature of the auxiliary strain rate $\underline{\varepsilon}^0(\tau, t)$. Except for this contribution, (12) define an ordinary linear viscoelastic (Maxwellian) behaviour; as for the additional term $\underline{\varepsilon}^0(\tau, t)$, the important thing to notice is that it can be considered as entirely known in advance for any $t$ when using (12), since all the involved quantities will have been determined on $[0, \tau]$ (see Section 3.3). It follows that $\underline{\varepsilon}^0(\tau, t)$ can be considered as a true 'stress-free strain rate' (or an eigenstrain rate) in a sense somewhat generalized from that of Eshelby and it will be dealt with as such in the sequel. Let us emphasize the fact that this eigenstrain rate $\underline{\varepsilon}^0(\tau, t)$ has nothing in common with the one which is considered in Kroner's type models, namely $\underline{\varepsilon}^{\text{in}}(t)$. In this treatment, similarly to Hill's elastoplastic one, the first-order stress dependence of the viscoplastic part of the strain rate is well accounted for through the tangent creep compliance $\mathbf{m}(\tau)$: this compliance, which could be considered as multibranched as well, accounts for the viscoelastic or elastic-viscoplastic flow.

In other words, the affine formulation allows the concentration problem to be transformed into a linear thermo(visco)elastic type problem: we put our linearized constitutive equations in the more classical form used for linear viscoelasticity (with eigenstrains) through the use of a Stieljes-type convolution product, denoted by $\otimes$ and defined as

$$ 
\begin{aligned}
[\mathbf{s}_{\tau} \otimes \underline{\sigma}](\tau, t) &= \frac{d}{dt} \left[ \int_{0}^{t} \mathbf{s}_{\tau}(\tau, t-u): \underline{\sigma}(u) \, du \right] \\
&= \int_{0}^{t} \mathbf{s}_{\tau}(\tau, t-u): \dot{\underline{\sigma}}(u) \, du + \mathbf{s}_{\tau}(\tau, t): \underline{\sigma}(0) \\
&= \int_{0}^{t} \frac{d\mathbf{s}_{\tau}}{du}(\tau, u): \underline{\sigma}(t-u) \, du + \mathbf{s}_{\tau}(\tau, 0): \underline{\sigma}(t),
\end{aligned} 
\quad (13) $$
where possible discontinuities of the involved quantities, considered as generalized functions, have to be taken into account; eqns (12) then read:

$$ \dot{\underline{\varepsilon}}(t) = [\mathbf{s}_{\tau} \otimes \underline{\sigma}](\tau, t) + \dot{\underline{\varepsilon}}^{0}(\tau, t), $$

$$ \mathbf{s}_{\tau}(\tau, t) = \mathbf{s} + \mathbf{m}(\tau)t, $$

$$ \mathbf{m}(\tau) = \frac{d\underline{g}}{d\underline{\sigma}}(\underline{\sigma}(\tau)), $$

$$ \dot{\underline{\varepsilon}}^{0}(\tau, t) = \underline{g}(\underline{\sigma}(\tau)) - \mathbf{m}(\tau)\underline{\sigma}(\tau) + [\underline{g}(\underline{\sigma}(t)) - \underline{g}(\underline{\sigma}(\tau)) - \mathbf{m}(\tau)(\underline{\sigma}(t) - \underline{\sigma}(\tau))][1 - H(t - \tau)]. \quad (14) $$

We can then be sure that the corresponding overall constitutive equations will read:

$$ \underline{E}(t) = [\mathbf{S}_{\tau} \otimes \underline{\Sigma}](\tau, t) + \underline{E}^{0}(\tau, t), $$

$$ \underline{\Sigma}(t) = [\mathbf{C}_{\tau} \otimes (\underline{E} - \underline{E}^{0})](\tau, t), \quad (15) $$

with the homogenized quantities $\mathbf{S}_{\tau}(\tau, t)$, $\mathbf{C}_{\tau}(\tau, t)$ and $\underline{E}^{0}(\tau, t)$ derived from the corresponding local ones through equations similar to (7), extended to viscoelasticity by use of the correspondence principle. Of course, the Laplace transform technique is expected to be of some use when dealing with such calculations (see hereafter).

### 3.2.2. General case

We go back to the more general behaviour defined by (8) and (9) where $\alpha_I(0 \le I \le N)$ could be, in the case of a crystalline material, some reference shear stress on the slip system (I). The linearization method is the same as above. For $t \ge \tau$, we write from (8) and (9):

$$ \dot{\underline{\varepsilon}}(t) = \mathbf{s}:\underline{\sigma}(t) + \underline{g}(\tau) + \mathbf{m}(\tau):[\underline{\sigma}(t) - \underline{\sigma}(\tau)] + \mathbf{n}(\tau)\cdot[\underline{\alpha}(t) - \underline{\alpha}(\tau)], $$

$$ \underline{\alpha}(t) = \underline{h}(\tau) + \mathbf{p}(\tau):[\underline{\sigma}(t) - \underline{\sigma}(\tau)] + \underline{q}(\tau)\cdot[\underline{\alpha}(t) - \underline{\alpha}(\tau)] \quad (t \ge \tau), \quad (16) $$

where $\underline{g}(\tau)$ and $\underline{h}(\tau)$ are shortened notations for $\underline{g}(\underline{\sigma}(\tau), \underline{\alpha}(\tau))$ and $\underline{h}(\underline{\sigma}(\tau), \underline{\alpha}(\tau))$, respectively, $\underline{q}(\tau)$ is a second-order $(N \times N)$ matrix whereas $\mathbf{n}(\tau)$ and $\mathbf{p}(\tau)$ are third order $((3 \times 3 \times N)$ and $(N \times 3 \times 3)$ resp.) matrices. Their components read, with $(i,j = 1, 2, 3)$ and $(I, J = 1, \dots, N)$:

$$ m_{ijkl}(\tau) = \frac{\partial g_{ij}}{\partial \sigma_{kl}}(\underline{\sigma}(\tau), \underline{\alpha}(\tau)), $$

$$ n_{ijl}(\tau) = \frac{\partial g_{ij}}{\partial \alpha_I}(\underline{\sigma}(\tau), \underline{\alpha}(\tau)), $$

$$ p_{Ikl}(\tau) = \frac{\partial h_I}{\partial \sigma_{kl}}(\underline{\sigma}(\tau), \underline{\alpha}(\tau)), $$

$$ q_{IJ}(\tau) = \frac{\partial h_I}{\partial \alpha_J}(\underline{\sigma}(\tau), \underline{\alpha}(\tau)), \quad (17) $$
where $\underline{q}(\tau)$ is assumed to be invertible with $\underline{q}^{-1}(\tau)$ its inverse and $\underline{I}_{\alpha}$ the associated unit $(N \times N)$ matrix.

Instead of (16), different expressions hold for $t \ge \tau$. To make them simpler, we first solve (16) as a differential equation for $\underline{\alpha}(t)$ with a continuity condition at $t = \tau$, namely,

$$ \underline{\alpha}(t) - \underline{\alpha}(\tau) = \underline{\alpha}(\tau, t) + \int_{0}^{t} \underline{e}^{(t-u)\underline{q}(\tau)} \cdot \mathbf{p}(\tau): \underline{\sigma}(u) \, du, \quad (18) $$

where the matrix exponential of $\underline{tq}(\tau)$ and $\underline{\alpha}(\tau, t)$ are defined by

$$ \underline{e}^{\underline{tq}}(\tau) = \sum_{0}^{\infty} \frac{1}{m!} (\underline{tq}(\tau))^m $$

and

$$ \underline{\alpha}(\tau, t) = \underline{q}^{-1}(\tau) \cdot [\underline{e}^{(t-\tau)\underline{q}(\tau)} - \underline{I}_{\alpha}] \cdot [\underline{h}(\tau) - \mathbf{p}(\tau): \underline{\sigma}(\tau)] - \int_{0}^{\tau} \underline{e}^{(\tau-u)\underline{q}(\tau)} \cdot \mathbf{p}(\tau): \underline{\sigma}(u) \, du. $$

Thus, we get for $t \ge \tau$:

$$ \begin{aligned} t \ge \tau \qquad \dot{\underline{\varepsilon}}(t) &= \mathbf{s}: \dot{\underline{\sigma}}(t) + \mathbf{m}(\tau): \underline{\sigma}(t) + \mathbf{n}(\tau) \cdot \int_{0}^{t} \underline{e}^{(t-u)\underline{q}(\tau)} \cdot \mathbf{p}(\tau): \underline{\sigma}(u) \, du + \dot{\underline{\varepsilon}}^{0}(\tau, t) \\ \dot{\underline{\varepsilon}}^{0}(\tau, t) &= \underline{g}(\tau) - \mathbf{m}(\tau): \underline{\sigma}(\tau) + \mathbf{n}(\tau) \cdot \underline{\alpha}(\tau, t). \end{aligned} $$

As before, $\dot{\underline{\varepsilon}}^{0}(\tau, t)$ can be given its actual potentially known value for $t \le \tau$ instead of the linearly back-extrapolated one. For the same reasons as above it can be considered as a true eigenstrain rate and dealt with as such when solving the corresponding linear thermoviscoelastic concentration problem.

Finally the linearized constitutive equations read:

$$ \begin{aligned} \dot{\underline{\varepsilon}}(t) &= [\mathbf{s}_{\tau} \otimes \dot{\underline{\sigma}}](\tau, t) + \dot{\underline{\varepsilon}}^{0}(\tau, t), \\ \mathbf{s}_{\tau}(\tau, t) &= \mathbf{s} + \mathbf{m}(\tau)t - \mathbf{n}(\tau) \cdot \underline{q}^{-1}(\tau) \cdot [\underline{I}_{\alpha}t + \underline{q}^{-1}(\tau) \cdot (\underline{I}_{\alpha} - \underline{e}^{\underline{t}\underline{q}(\tau)})] \cdot \mathbf{p}(\tau), \\ \dot{\underline{\varepsilon}}^{0}(\tau, t) &= \underline{g}(\tau) - \mathbf{m}(\tau): \underline{\sigma}(\tau) + \mathbf{n}(\tau) \cdot \underline{\alpha}(\tau, t) + \dot{\underline{\varepsilon}}(\tau, t)[1 - H(t - \tau)], \\ \dot{\underline{e}}(\tau, t) &= \underline{g}(t) - \underline{g}(\tau) - \mathbf{m}(\tau): (\underline{\sigma}(t) - \underline{\sigma}(\tau)) - \mathbf{n}(\tau) \\ & \qquad \cdot \left[ \underline{\alpha}(\tau, t) + \int_{0}^{t} \underline{e}^{(t-u)\underline{q}(\tau)} \cdot \mathbf{p}(\tau): \underline{\sigma}(u) \, du \right], \end{aligned} \quad (19) $$

with $\mathbf{m}(\tau)$, $\mathbf{n}(\tau)$, $\mathbf{p}(\tau)$, $\underline{q}(\tau)$ and $\underline{\alpha}(\tau, t)$ defined above. In the following, we suppose $\mathbf{s}_{\tau}$ to have the full diagonal symmetry (i.e. $s_{\tauijkl} = s_{\tauklji}$). Clearly, this condition implies some restrictions on the form of the flow and hardening functions $\underline{g}$ and $\underline{h}$. This point will be evoked later (Section 4).

According to the above discussion, we can once again presume that the cor-
responding overall constitutive equations will read as in (15). Equivalent integrated expressions can be used instead of (15) or of (19). For instance, at the local level, they have the form

$$ 
\begin{aligned}
\underline{\varepsilon}(t) &= [\mathbf{s}_{\tau} \otimes \underline{\sigma}](\tau, t) + \underline{\varepsilon}^{0}(\tau, t), \\
\underline{\varepsilon}^{0}(\tau, t) &= \underline{\varepsilon}^{0}(\tau, 0) + \int_{0}^{\tau} \underline{\varepsilon}^{0}(\tau, u) \, du, \\
\underline{\alpha}(t) &= \underline{\alpha}(0) + \int_{0}^{t} \underline{h}(u) \, du,
\end{aligned}
\qquad (20) $$

with $\underline{\varepsilon}^{0}(\tau, 0_{+}) = \underline{\varepsilon}(0_{+}) - \mathbf{s}: \underline{\sigma}(0_{+})$ when the initial response at $t = 0_{+}$ is elastic.

Notice that, instead of the proposed resolution of (16) as a differential equation for $\underline{\alpha}(t)$, we could have used the Laplace–Carson transform technique as well (see Navidi et al., 1996) with the same result. Since this technique is applied in the following to the resolution of the homogenization problem, we now indicate some aspects of the way it will be used hereafter.

### 3.2.3. Laplace–Carson transformation

For any time $\tau$, we are left with a linear viscoelastic homogenization problem which, except for the presence of an eigenstrain, is classical. The usual way to solve it is to convert this problem into an elastic one by use of the correspondence principle. A convenient tool for that is the Laplace–Carson transform technique which transforms any temporal function $f(t)$ into a symbolic function $f^{*}(p)$ of the complex variable $p$ according to

$$ f^{*}(p) = p \int_{0}^{+\infty} f(t) e^{-pt} dt. $$

Consequently, our transformed linearized constitutive equations read

$$ \underline{\varepsilon}^{*}(p) = \mathbf{s}_{\tau}^{*}(\tau, p): \underline{\sigma}^{*}(p) + \underline{\varepsilon}^{0*}(\tau, p) $$

with

$$ 
\begin{aligned}
s_{\tauijkl}^{*}(\tau, p) &= s_{ijkl} + \frac{1}{p}m_{ijkl}(\tau) + \frac{1}{p}n_{ijl}(\tau)w_{IJ}(\tau, p)p_{Jkl}(\tau) \\
(i, j, k, l = 1, 2, 3; I, J = 1, \dots, N; w(\tau, p) = [pI_{\alpha} - q(\tau)]^{-1}).
\end{aligned}
$$

After the resolution of the symbolic elastic homogenization problem, an inversion of the transformed quantities is needed. In the sequel, we use the collocation method according to which any complex function $f^{*}(p)$ is approximated as a finite series of the form
$$f^*(p) \approx a + \frac{b}{p} + \sum_{k=1}^{k=M} f_k \frac{1}{1 + p \theta_k},$$

so that $f(t)$ can be estimated as

$$f^*(p) \approx a + bt + \sum_{k=1}^{k=M} f_k (1 - e^{-\frac{t}{\theta_k}}).$$

Additional details are given in Appendix A.

### 3.3. The homogenization scheme

In order to derive the macroscopic response as well as the local stresses and internal variables at $t = \tau$, we need to know the entire local mechanical history from $t = 0$ to $t = \tau$ in every phase. This is already apparent before specifying any homogenization model when we use the correspondence principle in order to solve the concentration problem within the linearized framework which has been described above: referring to the elastic problem with eigenstrains whose general homogenization relations have been described in (7), we have now to write with obvious notations

$$ 
\begin{aligned}
\mathbf{C}_{\tau}(\tau, t) &= \langle [\mathbf{c}_{\tau} \otimes \mathbf{A}_{\tau}](\tau, t) \rangle, \\
\mathbf{S}_{\tau}(\tau, t) &= \langle [\mathbf{s}_{\tau} \otimes \mathbf{B}_{\tau}](\tau, t) \rangle, \\
\mathbf{A}_{r\tau}(\tau, t) &= [\mathbf{s}_{r\tau} \otimes \mathbf{B}_{r\tau} \otimes \langle \mathbf{s}_{\tau} \otimes \mathbf{B}_{\tau} \rangle^{-1\otimes}](\tau, t), \\
\mathbf{E}^0(\tau, t) &= \langle [{}^t\mathbf{B}_{\tau} \otimes \underline{\mathbf{e}}^0](\tau, t) \rangle = [\mathbf{S}_{\tau} \otimes \langle {}^t\mathbf{A}_{\tau} \otimes \mathbf{c}_{\tau} \otimes \underline{\mathbf{e}}^0 \rangle](\tau, t),
\end{aligned}
\quad (21) $$

where $\mathbf{c}_{r\tau}(\tau, t)$ (resp. $\mathbf{C}_{r\tau}(\tau, t)$) is the convolutive reciprocal of $\mathbf{s}_{r\tau}(\tau, t)$ (resp. $\mathbf{S}_{r\tau}(\tau, t)$) through the operation: $\mathbf{c}_{r\tau}(\tau, t) = [\mathbf{s}_{r\tau}(\tau, t)]^{-1\otimes}$.

Notice first that we need to know $\underline{\mathbf{e}}^0(\tau, t)$ for any $t \le \tau$, which means that we have to know the local stress and internal variables for $t \le \tau$ too. Moreover, except for very crude estimates such as Reuss–Voigt ones for which the concentration tensors $\mathbf{A}_{\tau}$ or $\mathbf{B}_{\tau}$ reduce to unity and the stress or strain concentration equations are time-invariant, the choice of specific concentration tensors and equations will enhance this necessity.

The same conclusion can be drawn when considering the Laplace–Carson transform of (21), namely:

$$ 
\begin{aligned}
\mathbf{C}_{\tau}^*(\tau, p) &= \langle \mathbf{c}_{\tau}^*(\tau, p) : \mathbf{A}_{\tau}^*(\tau, p) \rangle, \\
\mathbf{S}_{\tau}^*(\tau, t) &= \langle \mathbf{s}_{\tau}^*(\tau, p) : \mathbf{B}_{\tau}^*(\tau, p) \rangle, \\
\mathbf{A}_{r\tau}^*(\tau, p) &= \mathbf{s}_{r\tau}^*(\tau, p) : \mathbf{B}_{r\tau}^*(\tau, p) : \langle \mathbf{s}_{\tau}^*(\tau, p) : \mathbf{B}_{\tau}^*(\tau, p) \rangle^{-1}, \\
\underline{\mathbf{E}}^0*(\tau, p) &= \langle {}^t\mathbf{B}_{\tau}^*(\tau, p) : \underline{\mathbf{\epsilon}}^0*(\tau, p) \rangle = \mathbf{S}_{\tau}^*(\tau, p) : \langle {}^t\mathbf{A}_{\tau}^*(\tau, p) : \mathbf{c}_{\tau}^*(\tau, p) : \underline{\mathbf{\epsilon}}^0*(\tau, p) \rangle,
\end{aligned}
\quad (22) $$

with $\mathbf{c}_{r\tau}^*(\tau, p) = [\mathbf{s}_{r\tau}^*(\tau, p)]^{-1}$ and $\mathbf{C}_{r\tau}^*(\tau, p) = [\mathbf{S}_{r\tau}^*(\tau, p)]^{-1}$.

So we have to overcome two kinds of difficulties. On the one hand, an affine (by contrast with an incremental) approach leads to concentration tensors and equations,
which are necessary to determine the local and global mechanical responses at time $t = \tau$, depending themselves on these mechanical responses, so that the whole non-linear problem has an implicit nature. On the other hand, due to the viscoelastic coupling, we even need to know not only these responses at time $t = \tau$, but also at any previous time from $t = 0$. Consequently, we have to proceed in the following general way: the whole concerned time interval $[0, \tau]$ must be discretized with intermediate times $\tau_i$ and the associated unknown variables $\underline{\sigma}_r(\tau_i, \tau)$ and $\underline{\alpha}_r(\tau_i, \tau)$ as derived from an affine model defined at $\tau = \tau_i$; an iterative scheme has then to be elaborated so as to allow us the simultaneous determination of these variables which have, for some given macroscopic loading path, to satisfy at any $\tau_i$ the concentration and constitutive equations referring to $\tau_i$, which depend themselves on the set of values of these variables on the interval $[0, \tau_i]$. On the contrary, an incremental formulation, aiming at deriving at $t = \tau$ time increments of the local stresses and internal variables and of the associated time increments of the overall quantities, would have to use already determined histories of the local variables on $[0, \tau]$, as derived from their incremental (step by step) computation on this interval.

Additional details on this question are given in the next subsection 3.4 which is devoted to the self-consistent scheme. But we emphasize that the proposed treatment is quite general and is not restricted to this scheme (see e.g. Beurthey, 1997 for an application to the generalized self-consistent scheme and to polymer blends).

### 3.4. Application to the self-consistent scheme

Up to now, we have applied the affine linearization procedure at each point $\underline{x}$ of every phase ($r$) so that the linearized ‘comparison medium’ (according to Ponte Castaneda, 1991) is a continuously heterogeneous one. Drastic simplifications are needed now when choosing a specific homogenization scheme, so as to deal with uniform per phase tangent compliances and prestrain rates. These quantities will be defined at the phase stress averages $\underline{\sigma}_r$. Obviously this is a crucial (but usual) assumption which we have to keep in mind (see Section 5).

The concentration tensors $\mathbf{A}_{r\tau}(\tau, t)$ and $\mathbf{B}_{r\tau}(\tau, t)$ which are then uniform per phase are derived from the case of inhomogeneous elasticity with eigenstrains which has already been discussed in Section 2, by simply replacing the involved constant quantities by the Laplace–Carson transformed of the variable ones. According to classical results in elasticity, we have now:

$$ 
\begin{aligned}
\mathbf{A}_{r\tau}^*(\tau, p) &= (\mathbf{I} + \mathbf{P}_{r\tau}^*: (\mathbf{c}_{r\tau}^* - \mathbf{C}_{\tau}^*))^{-1}: \langle (\mathbf{I} + \mathbf{P}_{\tau}^*: (\mathbf{c}_{\tau}^* - \mathbf{C}_{\tau}^*))^{-1} \rangle^{-1}, \\
\mathbf{P}_{r\tau}^*(\tau, p) &= \int_{\Omega_r} \Gamma_{\tau}^*(\tau, p, \underline{x}_r, \underline{x}'_r) \, d\Omega'_r, \\
\mathbf{B}_{r\tau}^*(\tau, p) &= (\mathbf{I} + \mathbf{Q}_{r\tau}^*: (\mathbf{s}_{r\tau}^* - \mathbf{S}_{\tau}^*))^{-1}: \langle (\mathbf{I} + \mathbf{Q}_{\tau}^*: (\mathbf{s}_{\tau}^* - \mathbf{S}_{\tau}^*))^{-1} \rangle^{-1}, \\
\mathbf{Q}_{r\tau}^*(\tau, p) &= \mathbf{C}_{\tau}^* - \mathbf{C}_{\tau}^*: \mathbf{P}_{r\tau}^*: \mathbf{C}_{\tau}^*.
\end{aligned} 
\quad (23) $$

As for the concentration equation, it can be written for the local stresses $\underline{\sigma}_r(\tau)$ instead
of the strains. When restricting oneself, as done in the sequel, to identically shaped and oriented ellipsoidal (or spherical) inclusions, this equation reads simply

$$ \underline{\sigma}_r^*(\tau, p) = \mathbf{B}_{r\tau}^* : \underline{\Sigma}^* + \mathbf{B}_{r\tau}^* : \mathbf{Q}_\tau^* : (\underline{E}_\tau^0 * - \underline{\varepsilon}_{r\tau}^0 * ) \quad (24) $$

or, in the real space $t$:

$$ \underline{\sigma}_r(\tau, t) = [\mathbf{B}_{r\tau} \otimes \underline{\Sigma}](\tau, t) + [\mathbf{B}_{r\tau} \otimes \mathbf{Q}_\tau \otimes (\underline{E}_\tau^0 - \underline{\varepsilon}_{r\tau}^0)](\tau, t), \quad (25) $$

whereas internal variables $\underline{\alpha}_r(\tau, \tau)$ are defined by (20). Notice that, as suggested in the preceding subsection, eqns (23)–(25) will have to be used at previous times $\tau_i$ too.

Consideration of (25) makes the difference between the affine and incremental formulations quite clear for the self-consistent scheme: whereas a concentration equation for an incremental treatment would relate $\underline{\sigma}_r(\tau, t)$ and $\underline{\Sigma}(t)$, $\underline{E}_\tau^0(\tau, t)$, $\underline{\varepsilon}_{r\tau}^0(\tau, t)$ by an equation similar to (25), the time derivative of $\underline{\sigma}_r(\tau)$ from this equation would make additional terms appear, originating from the stress-dependence of the involved quantities.

Different iterative schemes could be conceived in order to determine intermediate values of the local stresses and internal variables at times $\tau_i$ prior to $t = \tau$ as argued in subsection 3.3. The one we have selected looks for the set $(\underline{\sigma}_r(\tau_i, \tau_i), \underline{\alpha}_r(\tau_i, \tau_i))$ by finding the solution through increasing values of $\tau_i$ from 0 to $\tau$: an iterative procedure has then to be used at each $\tau_i$. The corresponding numerical implementation may be summed up as follows. Let us consider e.g. a given prescribed macroscopic stress path $\underline{\Sigma}(t)$. In order to derive the macroscopic response $\underline{E}(\tau)$, we first determine, through an elastic analysis, the initial response $\underline{E}(0_+)$ and the associated initial local stresses, strains, eigenstrains and internal variables. Suppose we have succeeded in deriving these variables at increasing times $\tau_1, \tau_2, \dots, \tau_{i-1}$. In order to do the same at time $\tau_i = \tau_{i-1} + \delta\tau$ close to $\tau_{i-1}$, we first estimate a trial set of values for $\underline{\sigma}_r(\tau_i)$ and $\underline{\alpha}_r(\tau_i)$ from the known state at $\tau_{i-1}$, namely

$$ \begin{aligned} \underline{\sigma}_r^{(1)}(\tau_i) &= \underline{\sigma}_r(\tau_{i-1}, \tau_i) \\ &= [\mathbf{B}_{r,i-1} \otimes \underline{\Sigma}](\tau_{i-1}, \tau_i) + [\mathbf{B}_{r,i-1} \otimes \mathbf{Q}_{i-1} \otimes (\underline{E}_{i-1}^0 - \underline{\varepsilon}_{r,i-1}^0)](\tau_{i-1}, \tau_i) \\ \underline{\alpha}_r^{(1)}(\tau_i) &= \underline{\alpha}_r(\tau_{i-1}, \tau_i) = \underline{\alpha}_r(\tau_{i-1}, \tau_{i-1}) + \underline{h}_r(\tau_{i-1}, \tau_{i-1})\delta\tau. \end{aligned} $$

From these estimates, tentative local creep compliances and eigenstrains as well as overall creep compliance $\mathbf{S}_i^{(1)}(\tau_i, t)$, eigenstrain $\underline{E}_i^{(1)}(\tau_i, t)$, Green and concentration tensors $\mathbf{Q}_i^{(1)}(\tau_i, t)$, $\mathbf{B}_{r,i}^{(1)}(\tau_i, t)$ are evaluated as explained in Appendix B. This defines a new affine model which allows one to get updated estimates for $\underline{\sigma}_r(\tau_i, \tau_i)$ and $\underline{\alpha}_r(\tau_i, \tau_i)$, say

$$ \begin{aligned} \underline{\sigma}_r^{(2)}(\tau_i) &= [\mathbf{B}_{r,i}^{(1)} \otimes \underline{\Sigma}](\tau_i, \tau_i) + [\mathbf{B}_{r,i}^{(1)} \otimes \mathbf{Q}_i^{(1)} \otimes (\underline{E}_i^{0(1)} - \underline{\varepsilon}_{r,i}^{0(1)})](\tau_i, \tau_i), \\ \underline{\alpha}_r^{(2)}(\tau_i) &= \underline{\alpha}_r(\tau_{i-1}, \tau_{i-1}) + [(1-\theta)\underline{h}_r(\tau_{i-1}, \tau_{i-1}) + \theta\underline{h}_r^{(1)}(\tau_i, \tau_i)]\delta\tau, \end{aligned} $$

with $\theta \in [0, 1]$ ($\theta \approx 0.5$). Notice that the use of a linear interpolation for stresses and internal variables associated to the collocation method allows us to compute the
convolution products analytically. If these new estimates deviate too much from the previous ones (say by a relative variation larger than 1%), they are updated once more and the process is repeated until adequate convergence is obtained: a very small number of inner iterations (mostly only one) was found enough at each $\tau_i$.

We emphasize once more that, even in this 'step by step' version, this affine formulation does not reduce to an incremental one, which would deal with different concentration equations at each step. It is easily extended to the case of a prescribed macroscopic strain or strain rate path, as considered in some of the following illustrative applications to the prediction of the mechanical response of rate-dependent elastoplastic face-centred cubic polycrystals.

# 4. Application to FCC polycrystals

The volume sample is now constituted of a large number of single crystals. They have the same ellipsoidal (in practice, spherical) shape and are perfectly bounded at the boundary. They only differ by their lattice orientation which is defined by three Euler angles so that volume averages over the whole sample reduce to averages over these orientations.

## 4.1. Single crystal behaviour

In the following, we assume that the inelastic strain is only due to glide on slip systems. Denoting by $\underline{n}_s$ and $\underline{m}_s$, the unit vector normal to the slip plane and along the slip direction of the slip system (s), we define the orientation tensor of this slip system by

$$R_{s(kl)} = \frac{1}{2} [n_{s(k)} m_{s(l)} + n_{s(l)} m_{s(k)}],$$

so that the resolved shear stress on this system reads

$$\tau_s = \underline{R}_s \cdot \underline{\sigma}.$$

The traceless inelastic strain rate is the sum of the contributions of the shear rates $\dot{\gamma}_s^0$ on all the systems, i.e.

$$\dot{\underline{\epsilon}}_{\underline{\underline{s}}}^{\text{in}} = \sum_{s} \dot{\gamma}_{s} \underline{R}_{s},$$

where the shear rate $\dot{\gamma}_s$ is taken to depend on the stress through the resolved shear stress in the form

$$\dot{\gamma}_s = \dot{\gamma}_s^0 \left( \frac{\tau_s}{\tau_s^0} \right)^n.$$

The power-law exponent $n$ is the inverse of the material rate sensitivity, $\dot{\gamma}_s^0$ is a reference slip rate of system (s) and $\tau_s^0$ its reference shear stress. In general, $\tau_s^0$ depends on temperature and inelastic work. However, intragranular hardening effects are
neglected in the sequel so that the reference shear stress is kept constant under isothermal conditions. Consequently, the creep compliance at each step $\tau$ is given by

$$s_{ijkl}(\tau, t) = s_{ijkl} + n \sum_s \frac{\dot{\gamma}_s^0}{\tau_s^0} \left( \frac{\tau_s(\tau)}{\tau_s^0} \right)^{n-1} R_{s(ij)} R_{s(kl)} t$$

which respects the full diagonal symmetry condition supposed in Section 3.

## 4.2. Polycrystal overall responses

In the following, we consider an isotropic aggregate of face-centred cubic crystals. In each crystal, glide occurs on octahedral slip systems. All these twelve slip systems have the same reference shear strain rates $\dot{\gamma}^0$ and shear stresses $\tau^0$. In addition, the elastic behaviour is taken isotropic (Young modulus $E_y$ and Poisson ratio $\nu$). The material constants are: $n = 5$, $\tau^0/E_y = 10^{-3}$ and $\nu = \frac{1}{3}$.

The polycrystal is untextured: the considered inclusions are spherical and averages are obtained from an addition over eighty four equally weighted orientations. The accuracy of each simulation was evaluated by decreasing the number of orientations down to forty. Relative deviations between the two types of simulations never exceeded one percent. Time integrations of Taylor's model and Kroner-Weng's model were carefully conducted through a classical fourth-order Runge-Kutta method. As for our model, the computation of the Laplace-Carson transformed of the concentration tensors involves evaluations of the Green tensor which are reported in Appendix C.

The polycrystal is subjected to the macroscopic stress $\Sigma$. In the steady regime, Hutchinson (1976) has shown that the steady creep strain rate is a homogeneous function of degree $n$ of the overall stress. In addition, since the polycrystal is isotropic, the equivalent effective steady strain rate reads

$$\dot{E}_{\text{eq}} = \dot{\gamma}^0 \left( \frac{\Sigma_{\text{eq}}}{\Sigma^0} \right)^n. \qquad (26)$$

The reference stress $\Sigma^0$ governs the steady regime of the polycrystal and does not depend on the loading path. Moreover, $\Sigma^0/\tau^0$ depends only on the power-law exponent $n$ for a given radial loading.

### 4.2.1. Tensile stress–strain curves

The polycrystal is subjected to a monotonically increasing uniaxial tensile test at a fixed axial strain rate. The simulations of the polycrystal tensile–stress responses according to Taylor's model, Kroner-Weng's model and our model are shown in Fig. 1.

As previously mentioned in the rate-independent context, both self-consistent schemes lead to quite similar predictions in the early stage of the plastic flow. According to the present model, the behaviour is indeed quasi-elastic in this range so that the concentration tensors are very close to elastic Kroner's ones. But the differences
<figure><img src="image_2.png" /><figcaption>Fig. 1. Normalized tensile stress–strain curve for an isotropic rate-dependent elastoplastic FCC polycrystal as predicted from Taylor’s model, Kroner–Weng’s model and this model ($\tau^0/E_y = 10^{-3}$, $n = 5$, $E/\gamma^0 = 1$ with $E$ the overall prescribed axial strain rate).</figcaption></figure>

are increasing whereas the plastic flow occurs. As in the case of rate-independent elastoplasticity, the use of the elastic–viscoplastic tangent modulus instead of the elastic one increases the overall inelastic strain for a given macroscopic stress. Predictions are then softer than the classical ones.

The axial stress of the polycrystal tends towards a steady value for increasing plastic strains. According to (26), this steady effective axial stress reads

$$ \Sigma = \Sigma^0 \left( \frac{E_{\text{eq}}}{\gamma^0} \right)^{1/n} . $$

Like Bishop–Hill’s limit yield stress for rate-independent elastoplasticity, Taylor’s asymptotic value is an upper bound for this limit stress. For this power-law exponent, Taylor’s reference stress $\Sigma^0$ is

$$ \Sigma^0 \approx 2.76\tau^0, $$

which is in very good agreement with the result directly obtained by Hutchinson (1976). As discussed before, it is not surprising that Kroner–Weng’s steady stress
coincides with that of Taylor since the steady viscoplastic strain rate becomes uniform in the whole polycrystal when the stresses reach their steady values.

### 4.2.2. Creep

The polycrystal is now subjected to a creep tensile test. Under a stepwise application of the axial stress, single crystals are elastically stressed so that the instantaneous response at $t = 0$ is elastic. Consequently, the initial stresses are uniform in the polycrystal so that the effective strain at $t = 0$ reads

$$ \dot{E}(t) = \left\langle \dot{\gamma}^0(\tau^0)^{-n} \sum_s [R_s : \Sigma]^n R_s \mathrm{sg}(R_s : \Sigma) \right\rangle. $$

Since the single crystal responses display transverse isotropy around the axial loading direction, their study can be restricted to one spherical standard triangle. The evolutions of the local stresses for four orientations are reported in Fig. 2 with the corresponding orientations indicated in the standard triangle attached to this figure.

Whereas stress heterogeneities predicted by the present theory are lower than Kroner-Weng's model ones, the variations obtained in the transient regime are qualitatively similar. As before, both predictions deviate in the steady creep regime.

<figure><img src="image_3.png" /><figcaption>Fig. 2. Intragranular stress evolution during a tensile creep test ($\Sigma/2\tau^0 = 2$) for some crystals of an isotropic rate-dependent elastoplastic FCC polycrystal ($\sigma$ is the equivalent stress, $\tau^0/E_y = 10^{-3}$, $n = 5$).</figcaption></figure>
In this regime, Kroner–Weng’s model coincides again with that of Taylor which yields an upper bound for the effective creep strain rate. Once again, (26) is satisfied with very good agreement.

At odds with this result, our model predicts a significant creep strain heterogeneity. Therefore, creep strain rates have been evaluated for 21 equally weighted orientations in the standard triangle. In the resulting Fig. 3, we have reported the local normalized equivalent creep strain rates: the corresponding relative fluctuations can reach 40%. This clearly illustrates the grain-to-grain viscoplastic incompatibility during the steady creep regime. Conversely, the grain-to-grain quasi-compatibility associated with Kroner’s approach leads to overestimating the single crystal stresses while the corresponding steady strains are underestimated (see Fig. 4).

In Fig. 5 the overall axial creep strain response predicted by both self-consistent schemes has been reported. As before, the effective response predicted by the present model is definitely much softer. Of course our treatment needs much (≈ 10–30 times) longer computations than that of Weng but the results obtained look much more realistic.

More recently, Weng (1993) proposed a new version of an approximate self-

<figure><img src="image_4.png" /><figcaption>Fig. 3. Normalized intragranular equivalent steady creep strain rates for grains of an isotropic rate-dependent elastoplastic FCC polycrystal as predicted from this model ($\Sigma/2\tau^0 = 2$; $E$ and $\varepsilon$ denote the overall and local equivalent creep strain rates, respectively; $\tau^0/E_y = 10^{-3}$, $n = 5$).</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>Fig. 4. Intragranular equivalent strain evolution during a tensile creep test ($\Sigma/2\pi\tau^0 = 2$) for some crystals of an isotropic rate-dependent elastoplastic FCC polycrystal as predicted by Kroner-Weng's model (dashed lines) and this model ($\epsilon^e$ is the elastic response at $t = 0$; $\tau^0/E_y = 10^{-3}$, $n = 5$).</figcaption></figure>

consistent scheme taking into account the viscoplastic flow of the polycrystal in a simplified way, but the results obtained from this scheme still lead to uniform creep strains in this steady regime. This can be shown to be due to the simplified form of the interaction law.

**5. Conclusion**

This study shows that Hill's conception for extending the self-consistent model for nonlinear constitutive behaviour has not to be restricted to rate-independent elastoplasticity and viscoplasticity but can indeed be applied successfully to rate-dependent elastoplasticity. This conception relies on a linearization procedure applied to the local constitutive equations which leads to much softer mechanical interactions between the constituents than the ones derived from the elastic Kroner-type analysis adopted in previous treatments. The associated predictions for typical mechanical responses of FCC polycrystals have been shown to be, as expected, more realistic (i.e. softer) than that of Kroner-Weng. Similar conclusions have also been drawn from
<figure><img src="image_6.png" /><figcaption>Fig. 5. Evolution of the overall axial strain during a tensile creep test ($\Sigma/2\tau^0 = 2$) for an isotropic rate-dependent elastoplastic FCC polycrystal as predicted by Kroner-Weng's model and this model ($E^e$ is the overall instantaneous elastic response at $t = 0$; $\tau^0/E_y = 10^{-3}$, $n = 5$).</figcaption></figure>

other current applications (hexagonal crystals, Masson, 1998, polymer blends, coupling to rate-independent plasticity, Beurthey, 1997...).

Instead of Hill's original incremental formulation, a new 'affine' approach has been adopted here: it leads to the reduction of the concentration problem to a linear thermoviscoelastic one which can be solved by use of the correspondence principle and classical results of the homogenization theory of linear thermoelastic heterogeneous bodies. The viscoelastic coupling and the associated long range memory effect are responsible for a rather complex (though tractable) numerical implementation: this complexity is related to the dependence of the current response on the whole former load history. The same affine formulation can be applied more easily to rigid-viscoplasticity or nonlinear elasticity (Masson and Zaoui, 1997; Masson, 1998; Zaoui and Masson, 1998). In these cases, it has proved to yield softer responses than Hill's incremental and Hutchinson's secant formulations which have been shown (Gilormini, 1995) to lead to the violation of a nonlinear upper bound for the moduli. These apparent merits of the affine formulation which, contrary to the tangent method of Molinari et al. (1987) and Lebensohn and Tome (1993), does not lead to a Reuss-
type behaviour with an increasing nonlinearity for power-law creep and can deal with more general constitutive behaviour, are likely to be correlated with its relationship to the 'second-order procedure' recently proposed by Ponte Castaneda (1996). The main advantage of the affine formulation with respect to this procedure lies in its ability to be extended to elastoviscoplasticity which cannot be defined by one single potential. Though no variational approach leading to nonlinear bounds still exists in this field, the affine method is expected to yield more realistic predictions than the incremental one in this case too.

Nevertheless, we cannot expect the affine version of the self-consistent scheme to be immune to bounds violation since it is based on an approximate resolution of the basic underlying matrix/inclusion problem where both moduli are assumed to be uniform whereas they should vary from point to point with the local stresses and strains. According to the 'linear composite comparison medium' point of view, this intrinsic limit of Hill's conception itself, whatever formulation one has chosen, originates in the rough approximation of the substitution of the actual continuously heterogeneous linear comparison medium by a simplified body with tangent moduli which are uniform per phase and defined at some reference (usually phase average) stress or strain. Obviously, this suggests the need for further developments which, similarly to the 'modified secant moduli theory' (Suquet, 1995), would refer the phase properties to refined average quantities in connection with adequate variational procedures....

## Acknowledgements

Yves Rougier is gratefully acknowledged for his decisive initial contribution to this subject. We are also grateful to Pirouz Navidi for his valuable help for the numerical treatment.

## Appendix A: Numerical inversion of the Laplace–Carson transformation

Let us consider e.g. the Laplace–Carson transform of the effective compliance tensor $S_{\tau}^*(\tau, p)$. The asymptotic behaviour of $S_{\tau}(\tau, t)$ is first determined from one of the local compliances which obey

$$ \lim_{|p| \to \infty} \mathbf{s}_{\tau}^*(\underline{x}, \tau, p) = \mathbf{s}(\underline{x}), \quad \lim_{p \to 0} p \cdot \mathbf{s}_{\tau}^*(\underline{x}, \tau, p) = \mathbf{s}_{\tau}^{v}(\underline{x}, \tau), $$

with the purely viscous part $\mathbf{s}_{\tau}^{v}(\underline{x}, \tau)$ given by

$$ \mathbf{s}_{\tau}^{v}(\underline{x}, \tau) = \mathbf{m}_{\tau}(\underline{x}, \tau) - \mathbf{n}_{\tau}(\underline{x}, \tau) \cdot [\mathbf{q}_{\tau}(\underline{x}, \tau)]^{-1} \cdot \mathbf{p}_{\tau}(\underline{x}, \tau). $$

Assuming the effective symbolic properties to be continuous functions of the local compliances, we can define the effective elastic and purely viscous overall compliances $S$ and $S_{\tau}^{v}$ by
$$ \lim_{|p| \to \infty} S_{\tau}^{*}(\tau, p) = S \quad \lim_{p \to 0} p \cdot S_{\tau}^{*}(\tau, p) = S_{\tau}^{v}, \qquad (27) $$

which is equivalent (limit value theorem) to

$$ \lim_{t \to 0} S_{\tau}(\tau, t) = S, \quad \lim_{t \to \infty} S_{\tau}(\tau, t) \frac{1}{t} = S_{\tau}^{v}. \qquad (28) $$

According to the collocation method, we are looking for an approximation of $S_{\tau}(\tau, t)$ in the form of a Dirichlet series, namely

$$ S_{\tau}(\tau, t) = S + S_{\tau}^{v} t + \sum_{k=1}^{k=M} b_{\tau k} (1 - e^{-\frac{t}{\theta_{k}}}). $$

The $M$ relaxation times $\theta_k$ are chosen equispaced on a logarithmic scale over six decades, with $M \approx 20$. From the computed values of $S_{\tau}^{*}(\tau, p)$ at $M$ real points $p_l = (1/\theta_l)(1 \le l \le M)$, we are left with a linear system of $N$ equations, namely

$$ 1 \le l \le M \quad S_{\tau}^{*}(\tau, p_l) = S + S_{\tau}^{v} \frac{1}{p_l} + \sum_{k=1}^{k=M} b_{\tau k} \frac{1}{1 + p_l \theta_k}. $$

The resolution of this system then leads to the determination of the set of the unknown fourth-order tensors $(b_{\tau 1}, \dots, b_{\tau M})$ and to the expected Dirichlet series.

A similar treatment is used for the overall moduli, Green and concentration tensors. As pointed out by Turner and Tome (1993) for linear viscoelasticity, the elastic and purely viscous self-consistent effective compliance tensors are given by substituting $s_{\tau}^{*}$ in (22) by $s$ or $s_{\tau}^{v}$, respectively. This asymptotic analysis allows the collocation method to be focused on the transient regime, which greatly improves the quality of the numerical inversions.

## Appendix B: Numerical resolution of the self-consistent equation

Given an intermediate time $\tau_i$, the creep and purely viscous compliances and eigenstrains are first computed in each phase. Note once more that these eigenstrains have to be evaluated from $t = 0$ to $t = \tau_i$ for each phase.

To get the effective tangent behaviour at $\tau_i$, $M$ scalar collocation points $[p_1 < \dots < p_M]$ are chosen and the Laplace–Carson transform of every creep compliance is expressed at each collocation point. The resolution of the nonlinear equations (22) yields the purely viscous effective compliance $S_{\tau_i}^v$. The following iterative process is used to solve this equation: we start with the first estimate $S_{\tau_i}^{v(1)} = \langle s_{\tau_i}^v \rangle$. The stress concentration tensors are then derived in each phase, which improves the estimate of the purely viscous effective compliance $S_{\tau_i}^{v(2)}$ by using (22). This iterative process is repeated until adequate convergence is achieved. Starting from $S_{\tau_i}^{v(1)*}(\tau_i, p_l) = S + (1/p_l)S_{\tau_i}^v$, the same iterative scheme is used to evaluate the Laplace–Carson transform $S_{\tau_i}^{*}(\tau_i, p_l)$ at each collocation point $p_l$. Once these calculations are
performed, the creep compliances and concentration tensors are evaluated by the collocation method. The overall eigenstrain is then deduced from (21).

## Appendix C: Calculation of the concentration tensors

Any self-consistent model involves an evaluation of the Green tensor $P_r^{SC}$ for the effective medium. For the applications considered in Section 4, we restrict ourselves to a simple tension loading, so that the Laplace–Carson transform of the effective compliance tensor $S_{\tau}^*(\tau, p)$ displays transverse isotropy with respect to the tensile axis. In this case, evaluation of $P_r^{SC}$ only needs the use of ordinary integrals. Corresponding formulae are, for instance, detailed in Hutchinson (1970) (with some misprints rectified in Hutchinson, 1976). Then the concentration tensors are evaluated in each phase through (23).

However, due to the viscous incompressibility, these formulae do not strictly apply to the derivation of the asymptotic purely viscous effective compliance needed for the Laplace–Carson inversions. We then proceeded in an approximate manner by using the former calculation of $S_{\tau}^*(\tau, p)$ for a very small value $p_{\varepsilon}$ of $p$. As Hutchinson already did in the viscoplastic context (Hutchinson, 1976), this amounts to introducing a small compressibility for the derivation of $P_r^{SC}$. The effective purely viscous compliance $S_{\tau}^{v}$ is then given by

$$S_{\tau}^{v} \approx p_{\varepsilon} S_{\tau}^{*}(\tau, p_{\varepsilon}).$$

Numerical tests were performed to ensure that decreasing this small compressibility does not significantly change the final result.

## References

* 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. Modelisations du comportement d'alliages de polymeres. Ph.D. thesis, Ecole Poly-technique, France.

* Eshelby, J.D., 1957. The determination of the elastic field of an ellipsoidal inclusion and related problems. Proc. Roy. Soc. London A241, 376–396.

* Gilormini, P., 1995. Insuffisance de l'extension classique du modele autocoherent au comportement non lineaire. C.R. Acad. Sci. Paris 320 II, 115–122.

* Harren, S.V., 1991. The finite deformation of rate dependent polycrystals I and II. J. Mech. Phys. Solids 39, 345–383.

* Hershey, A.V., 1954. The elasticity of an isotropic aggregate of anisotropic cubic crystals. A.S.M.E. J. Appl. Mech. 21, 236–240.

* Hill, R., 1965. Continuum micro-mechanics of elastoplastic polycrystals. J. Mech. Phys. Solids 13, 89–101.

* Hutchinson, J.W., 1970. Elastic–plastic behaviour of polycrystalline metals and composites. Proc. Roy. Soc. London A319, 247–272.

* Hutchinson, J.W., 1976. Bounds and self-consistent estimates for creep of polycrystalline metals. Proc. Roy. Soc. London A348, 101–127.
Iwakuma, T., Nemat-Nasser, S., 1984. Finite elastic-plastic deformation of polycrystalline metals. Proc. Roy. Soc. London A394, 87-119.

Kroner, E., 1958. Berechnung der elastischen Konstaten des Vielkristalls aus den Konstaten des Einkristalls. Z. Phys. 151, 504-518.

Kroner, E., 1961. Zur plastischen Verformung des Vielkristalls. Acta. Metall. Mater. 9, 155-161.

Laws, N., McLaughlin, R., 1978. Self-consistent estimates for the viscoelastic creep compliance of composite materials. Proc. Roy. Soc. London A359, 251-273.

Lebensohn, R., Tome, C.N., 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.

Lipinski, P., Krier, J., Berveiller, M., 1990. Elastoplasticite des metaux en grandes deformations: comportement global et evolution de la structure interne. Rev. Phys. Appl. 25, 361-388.

Masson, R., 1998. Estimations non lineaires du comportement global de materiaux heterogenes en formulation affine. Ph.D. thesis, Ecole Polytechnique, France.

Masson, R., Zaoui, A., 1997. From rate-dependent to rate independent self-consistent modeling of elastoplastic multiphase materials. In: Khan, A.S. (Ed.), Physics and Mechanics of Finite Plastic and Viscoplastic Deformation. Neat Press, Maryland, U.S.A., pp. 209-210.

Molinari, A., Canova, G.R., Ahzi, S., 1987. A self-consistent approach of the large deformation polycrystals viscoplasticity. Acta. Metall. Mater. 35, 2983-2994.

Navidi, P., Rougier, Y., Zaoui, A., 1996. Self-consistent modelling of elastic-viscoplastic multiphase materials. In: Pineau, A., Zaoui, A. (Eds.), Micromechanics of Plasticity and Damage of Multiphase Materials. Kluwer Academic, Dordrecht, pp. 123-130.

Nemat-Nasser, S., Obata, M., 1986. Rate-dependent finite elastoplastic deformation of polycrystals. Proc. Roy. Soc. London A497, 343-375.

Ponte Castaneda, P., 1991. The effective mechanical properties of nonlinear composite materials. J. Mech. Phys. Solids 39, 45-71.

Ponte Castaneda, P., 1996. Exact second order estimates for the effective mechanical properties of nonlinear composite materials. J. Mech. Phys. Solids 44, 827-862.

Rougier, Y., Stolz, C., Zaoui, A., 1994. Self-consistent modelling of elastic-viscoplastic polycrystals. C.R. Acad. Sci. Paris 381 II, 145-151.

Suquet, P., 1995. Overall properties of nonlinear composites: a modified secant moduli theory and its link with Ponte Castaneda's nonlinear variational procedure. C.R. Acad. Sci. Paris 320 II, 563-571.

Taylor, G.I., 1938. Plastic strain in metals. J. Inst. Metals 62, 307-324.

Turner, P.A., Tome, C.N., 1993. Self-consistent modeling of viscoelastic polycrystals: application to irradiation creep and growth. J. Mech. Phys. Solids 41, 1191-1211.

Weng, G.J., 1981. Self-consistent determination of time-dependent behaviour of metals. J. Appl. Mech. 48, 41-46.

Weng, G.J., 1993. A self-consistent relation for the time-dependent creep of polycrystals. Int. J. Plast. 9, 181-198.

Zaoui, A., 1997. Structural morphology and constitutive behaviour of micro-heterogeneous materials. In: Suquet, P. (Ed.), Continuum Micromechanics, CISM Lecture Notes. Springer Verlag, Berlin, pp. 291-347.

Zaoui, A., Masson, R., 1998. Micromechanics-based modeling of plastic polycrystals: an affine formulation. In: NSF-IMM Symposium on Micromechanical Modeling of Industrial Materials. Seattle, July 1998, to appear.