# Structural Morphology and Constitutive Behaviour of Microheterogeneous Materials

André Zaoui

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

## To cite this version:

André Zaoui. Structural Morphology and Constitutive Behaviour of Microheterogeneous Materials. GAMM, 1999, Metz, France. pp.1-57, 10.1007/978-3-7091-2662-2_6. hal-00117859

HAL Id: hal-00117859

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

Submitted on 17 Oct 2022

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_2.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
STRUCTURAL MORPHOLOGY
AND CONSTITUTIVE BEHAVIOUR
OF MICROHETEROGENEOUS MATERIALS

A. Zaoui

Ecole Polytechnique, Palaiseau, France

ABSTRACT

One of the main specific aspects of continuum micromechanics is related to the fact that one
has generally to deal with ill-defined bodies: only partial information on the statistical
distribution of the constituent phases of the considered random inhomogeneous materials is
available.

The first chapter briefly reports the main classical ways to use such an information in
the context of linear elasticity for the derivation of bounds and estimates; attention is then
focused on isotropic particulate composites and especially on Hashin's composite spheres
assemblage. New estimates are proposed for the overall shear modulus which include
Hashin's bounds as well as the three-phase model estimate : improved Hashin-Shtrikman-
type bounds are conjectured from this analysis.

The second chapter presents a proof of this guess as well as a generalization of the
method used for arbitrary "morphological representative patterns"; the classical Hashin-
Shtrikman variational procedure is adapted to such a morphological analysis. Applications to
derivation of bounds and estimates for the overall elastic moduli are developed both for
isotropic and anisotropic distributions of the patterns.

In the third chapter, we lay the stress on behavioral aspects of this approach, with
special emphasis on the viscoelastic coupling. Whereas it is easy to study the influence of
morphological characteristics on the overall behaviour of linearly viscoelastic materials, this
is more difficult to do for nonlinear behaviours: Hill's treatment of the rate-independent
classical self-consistent model is highlighted and extended to rate-dependent behaviours but
this remains an approximation and new ways have still to be explored.
# 1. INTRODUCTION

## 1.1 Continuum micromechanics and structural morphology

In what follows, we assume the usual conditions which define the framework of continuum micromechanics to be satisfied: a random inhomogeneous material is considered which obeys macrohomogeneity requirements, which implies that the pertinent scale lengths of the body (referring to the size of the elementary "microscopic" inhomogeneities, of the "mesoscopic" volume elements and of the "macroscopic" structure and mechanical loading respectively) differ by one order of magnitude at least from each other. This separation of the scales allows "representative volume elements" (R.V.E.) to be defined : their details differ from one to the other but they lead to the same "homogeneous equivalent medium" (H.E.M.).

Nevertheless, the detailed microstructure of each R.V.E. is so complicated that it cannot be described completely: only some properties of the statistical distribution of the constituent phases are known. These properties may be quantitative, such as the knowledge of the volume fractions of the phases, of the two-point correlation functions of the elastic moduli etc. They may also be qualitative, so as e.g. to specify the isotropy or the symmetries of the anisotropy of the spatial phase distribution, the connectedness of one "matrix" phase, the average shape of inclusions belonging to another phase... We refer to all such informations as the "morphological" description of the inhomogeneous material.

The question is now: how to deal with such informations in order to derive either optimal bounds or pertinent estimates for the overall mechanical parameters? The main classical ways to answer this question are first reported briefly in what follows (§2) in the context of linear elasticity. The well-known solution of Eshelby's inclusion problem [1] has been used to derive several classical estimates whose morphological meaning is discussed. A general "systematic theory" for random elastic media [2] has been developed too which allows to derive a set of bounds for the overall elastic moduli associated with correlation functions of various orders.

These methods are shown to be unsatisfactory when materials exhibiting one continuous phase, such as the matrix of reinforced composites, are considered. For the rest of the chapter, we focus on this problem with special attention paid to Hashin's [3] composite spheres (or cylinders) assemblage (C.S.A. or C.C.A.). We build up a continuous set of estimates for the overall shear modulus of the isotropic C.S.A. which lie between the older bounds derived by Hashin by using a Voigt-Reuss type approach adapted to the C.S.A. whereas the three-phase model [4] is recovered as the self-consistent treatment of the problem. We can guess that two of the obtained estimates correspond to new Hashin-Shtrikman-type bounds: such improved bounds with respect to Hashin's ones should be obtained by taking better advantage of the isotropic distribution of the composite spheres.

This guess is demonstrated and developed systematically in the next chapter (§3) which combines the definition of "morphologically representative patterns" (M.R.P.) with the use of Hashin-Shtrikman variational procedure [5]. It is shown that, for an isotropic distribution of the patterns, explicit bounds can be derived from the solution of elementary composite spheres problems: this solution may be approximated analytically as closely as wanted. When the centers of the patterns obey an "ellipsoidal distribution", the numerical resolution of similar composite inclusions problems [6] lead to refined bounds whose comparison with older ones is meaningful. Finally new estimates can be defined too, especially in the context of the self-consistent procedure: M.R.P.-based generalized self-consistent schemes [7] using finite element computations for the solution of the associated inclusion elementary problems are suggested to yield an interesting alternative to the "unit cell" or periodic homogenization techniques.
## 1.2 Behavioral aspects

In the last chapter (§4), we lay the stress on behavioral aspects of the proposed approach. This can be done easily in the case of linear (non ageing) viscoelasticity through the Laplace transform technique. An illustration of the influence of morphology on the "long range memory effect" is given by comparing the relaxation spectra of the shear modulus of an isotropic two-phase material whose constituents obey a Maxwellian behaviour as predicted either by the classical or by the generalized self-consistent scheme: whereas in both cases the relaxation spectrum is a continuous one, which proves that the overall behaviour is no more Maxwellian, this spectrum is split up into two parts when the three-phase model is used: this property is characteristic of the continuous morphology of the matrix phase [8].

The connection between morphology and overall behaviour is much harder to be modelled in case of nonlinearity. For the rest of the chapter, attention is mainly focused on the nonlinear classical self-consistent scheme which already confronts us with most of the difficulties to be overcome. The advantages of Hill's treatment of elastoplasticity [9] with respect to former ones are first recalled. The difficulty to extend it to rate-dependent plasticity arises from the fact that, since stress and strain time derivatives of different orders are simultaneously present in the constitutive equations, Green techniques cannot be used directly. This difficulty can be solved [10] by combining the same kind of linearization procedure along the loading path and the use of the Laplace transform technique: at each incremental step, we are left with a linear viscoelastic behaviour with eigenstrains which can be converted into a (symbollically) classical elastic problem with eigenstrains.

Nevertheless an important basic problem of nonlinear micromechanics has not been solved by Hill's self-consistent treatment: it is associated to the fact that, at any stage, a nonlinear phase which is not homogeneously deformed cannot be defined by a unique set of instantaneous moduli (this is the case for the matrix for the self-consistent scheme). This difficulty is still increased with the different versions of the generalized self-consistent scheme which consider heterogeneously deformed phases. Some ways to appreciate this difficulty in connection with the recent theories of nonlinear bounds are suggested in conclusion.

# 2. CONTINUUM MICROMECHANICS AND MORPHOLOGICAL ANALYSIS

## 2.1 The need of a morphological analysis

It is well known that, except for quite special cases, e.g. when the assumption of a periodic microstructure can be made, the material content of a Representative Volume Element (R.V.E.) of a microinhomogeneous material cannot be described in a deterministic way. Even if a statistical point of view is adopted, the statistical description of the spatial distribution of the constituent phases cannot be performed completely. Consequently, we cannot aim more than at bounding or estimating the overall mechanical characteristics of the considered material. This is performed the more efficiently, the more completely the available information on the phase distribution is used or the more pertinently this information is obtained.

As a matter of fact, it is not easy to decide which kind of information, in addition to the minimal one consisting in the determination of the volume fractions of the constituents and their own mechanical characteristics, will lead to the closest bounds or to the most adequate estimates. It will depend on the specific morphology of the considered material: in some cases, the two-point correlation functions of some material parameters will be more efficient; in other cases, it will be the symmetry of the overall anisotropy or the fact that one phase is continuous and the other ones are not...

In any case, some morphological analysis has to be performed before any attempt to derive bounds or estimates. In what follows, we shall briefly recall first the main current
ways to build up such bounds and estimates in order to illustrate the way they could take account of such a morphological analysis. This will be done for linear elasticity in order to concentrate on these morphological aspects: behavioral ones will be considered later. Let us start with the most classical estimates of the overall elastic moduli which have been derived from the classical solution of Eshelby's inclusion problem.

## 2.2 Inclusion-based estimates of the overall elastic moduli

The principle of any method for the derivation of estimates of the overall elastic moduli $C^{est}$ or compliances $S^{est}$ of macrohomogeneous heterogeneous media is quite simple: since the local values $c_r$ and $s_r$ are uniform per phase, we only have to estimate the average stress $\overline{\sigma}_r$ or strain $\overline{\varepsilon}_r$ in each phase ($r$) for a given macroscopic stress $\Sigma$ or strain $E$, i.e. to estimate the average stress $\overline{B}_r^{est}$ or strain $\overline{A}_r^{est}$ concentration tensors, from which we derive:

$$ 
\begin{gathered}
E = \langle \overline{\varepsilon} \rangle = \sum_r f_r \overline{\varepsilon}_r = \sum_r f_r \overline{A_r^{est}} : E \\
C^{est} = \sum_r f_r c_r : \overline{A_r^{est}} = \langle c : A^{est} \rangle
\end{gathered}
\qquad (2.1) $$

or:

$$ 
\begin{gathered}
\Sigma = \langle \overline{\sigma} \rangle = \sum_r f_r \overline{\sigma}_r = \sum_r f_r \overline{B_r^{est}} : \Sigma \\
S^{est} = \sum_r f_r s_r : \overline{B_r^{est}} = \langle s : B^{est} \rangle
\end{gathered}
\qquad (2.2) $$

where $f_r$ is the volume fraction of phase ($r$) and $\langle \cdot \rangle$ stands for the volume average in the considered R.V.E..

Inclusion-based estimates may be defined by reference to Eshelby's solution [1] of the problem of an ellipsoidal inhomogeneity $H$ whose moduli $c$ differ from those, $C$, of the surrounding infinite matrix submitted to homogeneous strain $E^0$ at infinity. The strain tensor $\varepsilon_H$ in the inclusion is known to be uniform and given by

$$ \varepsilon_H = [I + S^E : C^{-1} : (c - C)]^{-1} : E^0 = (I + P : \delta c)^{-1} : E^0 \qquad (2.3) $$

where $S^E$ is the Eshelby tensor, $P = S^E : C^{-1}$, $\delta c = c - C$ and $I$ is the unit fourth-order tensor. We may now estimate $\overline{A_r^{est}}$ by considering that the mechanical state of each phase ($r$) in the R.V.E. submitted to the global strain $E$ is, in average, the same as the one of an ellipsoidal homogeneity $H_r$ with the same moduli $c_r$ embedded in a fictitious homogeneous medium with moduli $C^0$, submitted to some homogeneous strain $E^0$ at infinity. The shape and orientation of $H_r$ can be specified from what we know of the ($r$) phase domains geometry whereas $E^0$ is determined by the average condition

$$ \langle \varepsilon \rangle = \sum_r f_r \overline{\varepsilon}_r = E \qquad (2.4) $$

As for the moduli $C^0$, they can be chosen at will in order to express at best the specific morphology of the considered material; each choice of $C^0$ leads to a specific estimate $C_0^{est}$ through the corresponding $P^0$ value.
The associated set of equations reads:

$$
\begin{align}
\bar{\epsilon}_r &= (\mathbf{I} + \mathbf{P}_r^0 : \delta\mathbf{c}_r^0)^{-1} : \mathbf{E}^0 \\
\mathbf{E} &= \langle \bar{\epsilon} \rangle = \langle (\mathbf{I} + \mathbf{P}^0 : \delta\mathbf{c}^0)^{-1} \rangle : \mathbf{E}^0 \nonumber \\
\mathbf{A}_{r0}^{\text{est}} &= (\mathbf{I} + \mathbf{P}_r^0 : \delta\mathbf{c}_r^0)^{-1} : \langle (\mathbf{I} + \mathbf{P}^0 : \delta\mathbf{c}^0)^{-1} \rangle^{-1} \nonumber \\
\mathbf{C}_0^{\text{est}} &= \langle \mathbf{c} : (\mathbf{I} + \mathbf{P}^0 : \delta\mathbf{c}^0)^{-1} \rangle : \langle (\mathbf{I} + \mathbf{P}^0 : \delta\mathbf{c}^0)^{-1} \rangle^{-1} \nonumber \\
&= \sum_r (f_r \mathbf{c}_r : (\mathbf{I} + \mathbf{P}_r^0 : \delta\mathbf{c}_r^0)^{-1}) : (\sum_s f_s (\mathbf{I} + \mathbf{P}_s^0 : \delta\mathbf{c}_s^0)^{-1})^{-1} \nonumber
\end{align}
$$

Of course, similar results could be obtained for the compliances, namely, with obvious
notation:

$$
\begin{align*}
\bar{\sigma}_r &= (\mathbf{I} + \mathbf{Q}_r^0 : \delta\mathbf{s}_r^0)^{-1} : \Sigma^0 \\
\mathbf{Q}_r^0 &= \mathbf{C}^0 - \mathbf{C}^0 : \mathbf{P}_r^0 : \mathbf{C}^0 \\
\Sigma &= \langle \bar{\sigma} \rangle = \langle (\mathbf{I} + \mathbf{Q}^0 : \delta\mathbf{s}^0)^{-1} \rangle : \Sigma^0 \\
\mathbf{B}_{r0}^{\text{est}} &= (\mathbf{I} + \mathbf{Q}_r^0 : \delta\mathbf{s}_r^0)^{-1} : \langle (\mathbf{I} + \mathbf{Q}^0 : \delta\mathbf{s}^0)^{-1} \rangle^{-1} \\
\mathbf{S}_0^{\text{est}} &= \langle \mathbf{s} : (\mathbf{I} + \mathbf{Q}^0 : \delta\mathbf{s}^0)^{-1} \rangle : \langle (\mathbf{I} + \mathbf{Q}^0 : \delta\mathbf{s}^0)^{-1} \rangle^{-1} \\
&= \sum_r (f_r \mathbf{s}_r : (\mathbf{I} + \mathbf{Q}_r^0 : \delta\mathbf{s}_r^0)^{-1}) : (\sum_s f_s (\mathbf{I} + \mathbf{Q}_s^0 : \delta\mathbf{s}_s^0)^{-1})^{-1}
\end{align*}
\tag{2.6}
$$

Several classical estimates enter this framework for different C⁰ or S⁰ values. Voigt's (C<sup>V</sup> = <c>) and Reuss's (S<sup>R</sup> = 
) ones correspond to vanishing values of S⁰ and C⁰ respectively, which conforms with their extremal character and their underlying assumptions (uniform strains, as transmitted to the inclusion by a rigid matrix, and uniform stresses, as transmitted to the inclusion by an infinitely soft matrix, respectively).

Two other usual estimates, namely Mori-Tanaka's [11] and the self-consistent [12, 13] ones, may be obtained too. The first one corresponds to the choice $C^0 = C^{\text{Mat}}$, where $C^{\text{Mat}}$ are the elastic moduli of the matrix phase of a composite material. This is an approximate way for expressing the continuity of the matrix of such reinforced composites: in this case, the reinforcements as well as the matrix itself are supposed to be surrounded by the continuous matrix, which corresponds obviously to an overestimation of the mechanical role of this phase since interactions between particles or fibres as well as their strengthening effect on the matrix response are neglected. From Walpole's interpretation [14] of Hashin-Shtrikman [15] bounds in terms of matrix/inclusion situations, it is easy to show that for spherical particles and isotropy, this estimate coincides with Hashin-Shtrikman's lower bounds when, as usual, the matrix is softer than the reinforcements. Note that this derivation only insures the reciprocity of the overall moduli and compliances ($S^{\text{MT}} = (C^{\text{MT}})^{-1}$) but not their symmetry when particles with different aspect ratios and orientation are considered.

As for polycrystals, their phases (i.e. grains with identical lattice orientation and shape) are rather disordered so that they play a similar morphological role: each grain with a given lattice orientation and shape is surrounded by many other grains and, as a whole, the set of similarly shaped and oriented grains is surrounded by almost all the other phases. Consequently, a more adequate choice for the moduli Cº would be to take them as the researched moduli Cest themselves. This is the basic assumption of the "classical self-
consistent scheme" (C.S.C.S.) whose moduli $C^{CSCS}$ and compliances $S^{CSCS}$ satisfy, from (2.5) and (2.6), the equations

$$ 
\begin{align} 
C^{CSCS} &= \langle c : (I + P^{CSCS} : \delta c^{CSCS})^{-1} \rangle : \langle (I + P^{CSCS} : \delta c^{CSCS})^{-1} \rangle^{-1} \nonumber \\
\text{or} \quad & \langle \delta c^{CSCS} : (I + P^{CSCS} : \delta c^{CSCS})^{-1} \rangle \ge 0 \tag{2.7} \\
S^{CSCS} &= \langle s : (I + Q^{CSCS} : \delta s^{CSCS})^{-1} \rangle : \langle (I + Q^{CSCS} : \delta s^{CSCS})^{-1} \rangle^{-1} \nonumber \\
\text{or} \quad & \langle \delta s^{CSCS} : (I + Q^{CSCS} : \delta s^{CSCS})^{-1} \rangle \ge 0 \nonumber
\end{align} 
$$

with $P^{CSCS}$, $\delta c^{CSCS}$ and $Q^{CSCS}$, $\delta s^{CSCS}$ formed from $C^{CSCS}$ and $S^{CSCS}$.

Many other estimates, of course, could be defined with other choices for $C^0$ or $S^0$, for instance $C^V$, $S^R$, $C^{MT}$ etc., but without any clear morphological meaning. An illustration is given on Fig. 1 for the estimation $\mu^{est}$ of the overall shear modulus of a two-phase isotropic incompressible material at some given volume fraction according to the foregoing method. The ruling formula reads here:

$$ \mu^{est} = \frac{2\mu_1\mu_2 + 3\mu^0 <\mu>}{2\mu_1\mu_2 <\mu^{-1}> + 3\mu^0} \quad (2.8) $$

where $\mu^0$ can vary from zero to infinity. Note that Hashin-Shtrikman's bounds are obtained for $\mu^0 = \mu_1$ and $\mu^0 = \mu_2$ and that the successive estimates $\mu^{(0-)} (= \mu_1)$, $\mu^{(2-)} (= \mu^{HS-} = \mu^{MT})$, ..., $\mu^{(2n-)}$ ... as well as $\mu^{(0+)} (= \mu_2)$, $\mu^{(2+)} (= \mu^{HS+})$, ..., $\mu^{(2n+)}$ converge towards the self-consistent estimate, as expressed by Kröner's recursive formula for this estimate [16], which could be put in connection with his theory of "graded disorder".

<figure><img src="image_3.png" /><figcaption>Fig. 1 An infinite set of estimates for the overall shear modulus of an isotropic incompressible two-phase elastic material, for a given volume fraction (the inclined straight line is the first bisector of the positive quadrant)</figcaption></figure>
## 2.3 A systematic bound theory

There is a systematic way to obtain bounds for the overall elastic moduli of random media which get sharper and sharper by integrating more and more information on the spatial distribution of the local moduli. Without entering into the details of their derivation, let us emphasize some aspects of this "systematic theory" [2] in view of our further developments:

- the fundamental equation which is considered is the so-called "Lippman-Schwinger" equation which derives from the application of the Green technique to the problem of a given inhomogeneous body with moduli $c(\bar{x})$: if $G^0(\bar{x}, \bar{x}')$ is the Green tensor for an arbitrary reference homogeneous medium with the same geometry, the moduli $C^0$ and fixed boundary, the strain field in the heterogeneous body satisfies the integral equation

$$ \varepsilon(\bar{x}) + \Gamma^0(\bar{x}, \bar{x}') * \delta c^0(\bar{x}'): \varepsilon(\bar{x}') = \varepsilon^0 \quad (2.9) $$

where $\varepsilon^0$ is the uniform strain field which would exist if the body was homogeneous and the sign $*$ indicates a convolution-type product. Here $\Gamma^0(\bar{x}, \bar{x}')$ is the "modified Green strain operator", defined from $G^0(\bar{x}, \bar{x}')$ by:

$$ \Gamma^0_{ijkl}(\bar{x}, \bar{x}') = \frac{\partial^2}{\partial x_j \partial x_l'} G^0_{ik}(\bar{x}, \bar{x}')_{(ij)(kl)} \quad (2.10) $$

where the index $(ij)$ stands for symmetrization with respect to $i$ and $j$. Note that $\Gamma^0(\bar{x}, \bar{x}')$ is singular at $\bar{x} = \bar{x}'$ and behaves as a Dirac function at this point.

- if the distribution $c(\bar{x})$ was known, (2.9) could be solved, at least formally, in order to derive the effective moduli $C^{\text{eff}}$:

$$ \begin{gathered} C^{\text{eff}} = \langle c : (\tilde{I} + \Gamma^0 * \delta c^0)^{-1} \rangle : \langle (\tilde{I} + \Gamma^0 * \delta c^0)^{-1} \rangle^{-1} \\ \text{or} \quad \langle \delta c^{\text{eff}} : (\tilde{I} + \Gamma^0 * \delta c^0)^{-1} \rangle = 0 \end{gathered} \quad (2.11) $$

$$ \text{with} \quad \delta c^{\text{eff}} = c - C^{\text{eff}} $$

where the exponent $(-1)$ refers to an operator inversion and $\tilde{I}$ is the Dirac unit tensor. Similar expressions are valid for the effective compliances $S^{\text{eff}}$ by changing $c(\bar{x})$ into $s(\bar{x})$ and $\Gamma^0(\bar{x}, \bar{x}')$ into its stress counterpart $\Delta^0(\bar{x}, \bar{x}')$ defined by:

$$ \Delta^0 = C^0 - C^0 * \Gamma^0 : C^0 \quad (2.12) $$

- the same Green technique may be used to build up sets of kinematically (or statically) admissible strain $\varepsilon'(\bar{x})$ (or stress $\sigma^*(\bar{x})$) fields from any polarization (symmetric) stress field $p(\bar{x})$ through the equations:

$$ \begin{gathered} \varepsilon'(\bar{x}) + \Gamma^0(\bar{x}, \bar{x}') * p(\bar{x}') = \varepsilon^0 \\ \langle \varepsilon' \rangle = E \end{gathered} \quad (2.13) $$

$$ \sigma^* = C^0 : \varepsilon' + p $$

Such trial fields may be used in the potential and complementary energy theorems in order to bound the effective moduli and compliances, according to the available information on the distribution of the local quantities.
when such an information is given under the form of the first *n*-point correlation functions of the local moduli or compliances, it can be used to yield *n*-th order bounds which are the closer, *n* the larger.

- a special case which is worth-mentioning is concerned with materials obeying a *n*-rank disorder condition [17]: if $E^0$ and $F^0(\bar{x}, \bar{x}')$ are the (singular) local part and the (regularized) long-range part of the Green operator respectively, so that this operator may be put in the form

$$ \Gamma^0(\bar{x}, \bar{x}') = E^0 \delta(\bar{x}, \bar{x}') + F^0(\bar{x}, \bar{x}') \quad (2.14) $$

this condition reads:

$$ \langle c^\dagger : (F^0 * c^\dagger)^p \rangle = 0, \quad \forall p = 1 \text{ to } (n-1) \quad (2.15) $$

An interesting property of such a class of materials is that well-known bounds or estimates are recovered as special cases: Voigt-Reuss's bounds correspond to $n = 1$, Hashin-Shtrikman's ones to $n = 2$, whereas the self-consistent estimate is obtained for $n \to \infty$, the so-called "condition of perfect disorder". In this case, all the contributions of the Green operator reduce to their local part, so that (2.11) reduces to:

$$ \begin{gathered} C^{CSCS} = \langle c : (I + E^{CSCS} : \delta c^{CSCS})^{-1} \rangle : \langle (I + E^{CSCS} : \delta c^{CSCS})^{-1} \rangle^{-1} \quad (2.16) \\ \text{or} \quad \langle \delta c^{CSCS} : (I + E^{CSCS} : \delta c^{CSCS})^{-1} \rangle = 0 \end{gathered} $$

which coincides with (2.7a) for spheres and isotropy; (2.7b) could be obtained in the same way by use of the local part of $\Delta^0(\bar{x}, \bar{x}')$ instead of $\Gamma^0(\bar{x}, \bar{x}')$.

- one could think this systematic theory to yield the final solution of the problem of taking morphology into account in homogenizing elastic random media. What makes the situation less optimistic is that, in addition to the mathematical complexity of the practical computation of bounds from given correlation functions of the elastic moduli, these functions can hardly be extracted from an experimental investigation beyond $n = 3$ whereas many primary morphological characteristics, such as the geometrical connectedness of one phase or its inclusive nature could only be expressed through correlation functions of quite high an order...

## 2.4 Treatments of connectedness

This question has been for a long time a bone of contention between people interested in composites and those concerned with polycrystals. For the former, any homogenization model should obviously express, for instance, the fact that a porous medium can have some mechanical strength up to a porosity of 100% and they found the self-consistent scheme to be physically meaningless because of its prediction of vanishing elastic moduli for a porosity larger than .5 or .6; on the contrary, metallurgists have considered for a long time Hashin-Shtrikman's bounds as useless for their usual (polycrystalline) materials whereas they found the self-consistent prediction of a rigid overall behaviour of a two-phase (duplex) metal as soon as the volume fraction of the quasi-rigid second phase is larger than .4 or .5 as quite sensible and meaningful...

If one has in mind any care for the role of morphology in micromechanical modelling, one may consider such a dispute as out-of-date: it only reflects the fact that composites and polycrystals refer to basically different morphological types and cannot be satisfactorily described by the same models, as already discussed. In this section, we focus attention on the treatment of connectedness according to various models proposed in literature before proposing a new method to do it in an improved manner.
From what we have seen before, we can conclude that the question is still open. Among the different inclusion-based estimates reported in section 2.2, the Mori-Tanaka's one offers a rather crude and indirect way to express the geometrical continuity of the matrix of a composite. The fact that it coincides with Hashin-Shtrikman's lower bound for an isotropic particle-reinforced composite illustrates its overestimation of this morphological property and we could hope that this property could be accounted for in a less extremal and indirect way. At the opposite, the self-consistent estimate is clearly inadequate to this purpose and is rather adapted to polycrystal-type disordered and dispersed morphologies.

As for the systematic theory of bounds, which includes Hashin-Shtrikman's ones as a special case, it does not seem adequate anymore since it is centered around the idea of morphological description through point correlation functions. As above-mentioned, such a description would need correlation functions of a quasi-infinite order to have a chance to express the connectedness of one particular phase, which is untractable in practice.

Two other ways have been proposed in literature in view of a better expression of the matrix continuity of a composite. The first one refers to the so-called "differential self-consistent scheme" (D.S.C.S.) [18, 19] and the second one to Hashin's "composite sphere assemblage" (C.S.A.) [3]. We now briefly discuss the first method before paying sharper attention to the second one.

## 2.4.1 The differential self-consistent scheme

Instead of considering a particle (or fiber)-reinforced composite in its actual state, with the reinforcement volume fraction $f$, we take it as the final result of a differential process which starts from the matrix phase alone. The second phase is added progressively, so as to have its current volume fraction $c$ raised from 0 to $f$; at each stage, the medium is re-homogenized according to the self-consistent procedure; the latter may be defined within the dilute approximation, since the considered second phase volume fraction is infinitesimal.

This differential process may be derived as follows: at a current stage, we start from a homogeneous mixture with the overall moduli $C(c)$ containing the matrix and the second phases with the volume fractions $(1-c)$ and $c$ respectively. We add some infinitesimal volume of the second phase with moduli $c_2$, so that its final volume fraction be $c+dc$: it is easy to see that, in this new mixture, the volume fraction of the newly added second phase is $dc/(1-c)$ whereas the one of the previous homogeneous mixture is now $(1-c-dc)/(1-c)$. The overall moduli of this modified mixture, say $C(c+dc)$, are calculated from (2.7). If (2.7a) is specified for a two-phase material, with volume fractions $f_1$ and $f_2$, the overall moduli $\bar{C}$ can be derived from the equation

$$f_1(c_1 - \bar{C}):[I + \bar{P}_1:(c_1 - \bar{C})]^{-1} + f_2(c_2 - \bar{C}):[I + \bar{P}_2:(c_2 - \bar{C})]^{-1} = 0 \quad (2.17)$$

where $\bar{P}_1$ and $\bar{P}_2$ are defined for the searched medium with moduli $\bar{C}$. If now (2.17) is specified for the considered situation, i.e. for $f_1 = (1-c-dc)/(1-c)$, $c_1 = C(c)$, $f_2 = dc/(1-c)$, $\bar{C} = C(c+dc) = C(c) + dC$ and if only terms of first order in $c$ are kept, we find the following differential equation:

$$dC = \frac{dc}{1-c}(c_2 - C):(I + P_2:(c_2 - C))^{-1} \quad (2.18)$$

where $C = C(c)$ and $P_2$ may be defined from $C$ instead of $C(c+dc)$. Integration of (2.18) has to use the initial value $C(0) = c_1$. This is not easy to be performed in closed form in the general case, but it can be achieved better in special cases.

For instance, for an isotropic composite with spherical particles and isotropic constituents, we have the classical relations:
$$C = 2\mu K + 3kJ$$

$$c_2 = 2\mu_2 K + 3k_2 J$$

$$I = J + K \qquad I_{ijkl} = \frac{1}{2}(\delta_{ik}\delta_{jl} + \delta_{il}\delta_{jk}) \qquad J_{ijkl} = \frac{1}{3}\delta_{ij}\delta_{kl} \qquad (2.19)$$

$$P_2 = \frac{\beta}{2\mu} K + \frac{\alpha}{3k} J \qquad \beta = \frac{6(k + 2\mu)}{5(3k + 4\mu)} \qquad \alpha = \frac{3k}{3k + 4\mu}$$

Consequently, the differential equation (2.18) is split into two scalar equations, for the shear and the bulk moduli $\mu$ and $k$, namely:

$$d\mu = \frac{dc}{1-c} \frac{\mu_2 - \mu}{1 + \frac{\beta(\mu_2 - \mu)}{\mu}} \qquad dk = \frac{dc}{1-c} \frac{k_2 - k}{1 + \frac{\alpha(k_2 - k)}{k}} \qquad (2.20)$$

These are highly nonlinear equations and, due to the dependence of $\alpha$ and $\beta$ on $k$ and $\mu$, they are coupled too. In order to illustrate simply the possibilities of this scheme, let us consider the case of incompressibility for both phases, and so for the composite too. In the first equation (2.20), we can set $\beta = 2/5$, so that the equation for the shear modulus reads now:

$$\frac{(3\mu + 2\mu_2) d\mu}{5\mu (\mu_2 - \mu)} = \frac{dc}{1-c} \qquad (2.21)$$

After integration, we get

$$\frac{\mu^{\text{DSCS}}}{\mu_1} \left[ \frac{\mu_2 - \mu^{\text{DSCS}}}{\mu_2 - \mu_1} \right]^{-5/2} = (1-f)^{-5/2} \qquad (2.22)$$

This result can be compared, for $\mu_1 < \mu_2$, with the (classical) self-consistent estimate $\mu^{\text{CSCS}}$ as well as with Hashin-Shtrikman's lower bound (which coincides here with Mori-Tanaka's estimate $\mu^{\text{MT}}$). They are given by the following relations:

$$\frac{\mu^{\text{CSCS}}}{\mu_1} = \frac{(5f-2)}{6} \frac{\mu_2}{\mu_1} + \frac{3-5f}{6} + \frac{1}{6} \sqrt{[(2-5f)\frac{\mu_2}{\mu_1} + 5f-3]^2 + \frac{24\mu_2}{\mu_1}} \qquad (2.23)$$

$$\frac{\mu^{\text{MT}}}{\mu_1} = \frac{3(1-f)\mu_1 + (2+3f)\mu_2}{(3+2f)\mu_1 + 2(1-f)\mu_2}$$

It can be checked that, as expected, the three considered estimates are ordered according to: $\mu^{\text{MT}} < \mu^{\text{DSCS}} < \mu^{\text{CSCS}}$. In order to pay special attention to the treatment of connectedness, it is interesting to make the inclusions harder and harder (i.e. $\mu_2 \to \infty$). It is easy to get in this case (Fig. 2):

$$\frac{\mu^{\text{DSCS}}}{\mu_1} \to (1-f)^{-5/2} \qquad \frac{\mu^{\text{MT}}}{\mu_1} \to \frac{2+3f}{2(1-f)} \qquad (2.24)$$

$$\frac{\mu^{\text{CSCS}}}{\mu_1} \to \frac{2}{2-5f}, \quad f \le \frac{2}{5} \qquad \frac{\mu^{\text{CSCS}}}{\mu_1} \to \infty, \quad f \ge \frac{2}{5}$$
<figure><img src="image_4.png" /><figcaption>Fig. 2 Comparison of the CSCS, DSCS and MT estimates of the inverse of the effective shear modulus for an isotropic incompressible composite with rigid particles.</figcaption></figure>

Whereas the three estimates behave similarly at low volume fractions $f$ of the rigid inclusions, namely as $\mu/\mu_1 \approx 1 + 5f/2$, they differ more and more when $f$ is increasing: the two first ones remain bounded for any $f$ and they tend to infinity only when $f \to 1$ whereas the third one is unbounded as long as $f \ge .4$. These conclusions clearly express the fact that the two first estimates take some geometrical continuity of the matrix into account whereas the classical self-consistent scheme deals with some more disordered morphology so as to allow a rigid skeleton to be formed beyond some percolation-like threshold.

As far as composite materials with a well-connected matrix are concerned, we now have two models at disposal. Can we conclude that the differential self-consistent scheme, which is more complex than Mori-Tanaka's one, is more relevant to express the matrix connectedness? On the one hand, we can suspect that it does not present the exaggerated character we have already underlined when discussing Mori-Tanaka's model. On the other hand, it is not beyond any reproach: it can be developed in closed form for special cases only so that its applications are more limited; moreover, though yielding the expected qualitative properties for particle (or fiber)-reinforced composites, its morphological meaning is not so clear.

As a matter of fact, the step by step re-homogenization procedure after every increase of the particle (or fiber) content is rather artificial, despite its practical efficiency. It has been argued that this procedure is a quite natural one when a wide range of particle sizes is to be taken into account: the first steps would allow to deal with the smallest particles and the
following ones would be concerned with larger and larger ones. Nevertheless, this cannot be
understood literally since the elementary matrix/inclusion problem (with an infinite matrix)
has no absolute length scale and cannot distinguish by itself smaller and larger inclusions...
That is why we cannot consider that the problem of matrix connectedness representation has
received a final solution through the differential self-consistent scheme and other ways to do
so have better to be explored.

2.4.2 The composite spheres assemblage

A much more direct and simple way to represent the matrix connectedness of a
composite material has been proposed by Hashin [3] as soon as 1962. Dealing first with an
isotropic particulate material, the "composite spheres assemblage" (C.S.A.), which has been
followed shortly after by the "composite cylinders assemblage" (C.C.A.) for fiber-reinforced
transversely isotropic materials, consists in an unbounded set of contiguous similar
composite spheres of all sizes, including vanishing ones, so as to be able to fill up the whole
space (see Fig. 3). Each composite sphere, with some radius b, has a spherical core, with
the radius a, made of the reinforcement material and a concentric spherical shell made of the
matrix material. The ratio (a/b) equals f^{1/3} where f is the particle volume fraction (for the
C.C.A., this ratio is f^{1/2}). Since the composite spheres touch each other and all the particles
are surrounded by some matrix material, the matrix is unambiguously connected, just as
wished.

The point is now to derive the overall moduli from those of the constituents and the
particle volume fraction f. Let us consider an isotropic C.S.A. and isotropic phases with the
shear and bulk moduli μ₁, μ₂ and k₁, k₂, where the indices 1 and 2 refer to the matrix and
the particles respectively. Obviously, there is an infinity of ways to realize an isotropic
C.Ş.A. with the parameters μ₁, μ₂, k₁, k₂ and f. Thus the overall moduli μ<sup>C.S.A.</sup> and k<sup>C.S.A.</sup>
must be either bounded or estimated. The first approach to the problem, as proposed by
Hashin, aimed at bounding the overall moduli. More than fifteen years later, an estimate was
derived by Christensen and Lo [20] through the so-called "three-phase model".

<figure><img src="image_5.png" /><figcaption>Fig. 3 The composite spheres assemblage and Hashin's assumption for the derivation of bounds for its overall moduli.</figcaption></figure>
2.4.2.1. Hashin's bounds for the C.S.A.

In order to build up admissible strain and stress fields to be used in the energy
theorems, Hashin considered homogeneous conditions at infinity of the form

$$
\bar{u}^g = \mathbf{E}^\infty \cdot \bar{x} \quad \text{or} \quad \bar{T}^g = \Sigma^\infty \cdot \bar{n} \tag{2.25}
$$

and he imagined to transfer these conditions to the boundary of every composite sphere of
the assemblage and to solve the corresponding elastic problems. It is noticeable that, when
doing so, the required continuity conditions at the interface between any two contiguous
spheres are automatically satisfied so that, if the two elementary problems can be solved, we
are left with a kinematically admissible strain field and a statically admissible stress field
respectively, which are defined throughout the whole assemblage.

From a practical point of view, isotropy allows to split each elementary problem into
two ones, corresponding to purely radial or shear loading conditions, so that the overall bulk
and shear moduli can be bounded separately. The solution of each problem derives from
basic Love's results [21]. The conclusion is twofold: the overall bulk modulus k<sup>C.S.A.</sup> is
unambiguously determined, since the upper and lower bounds coincide with the value

$$
\frac{k^{\text{CSA}}}{k_1} = 1 + \frac{f}{\frac{1}{k_2/k_1 - 1} + \frac{3k_1(1-f)}{3k_1 + 4\mu_1}} \quad (2.26)
$$

whereas the overall shear modulus $\mu^{C.S.A.}$ is not. For reasons which will appear later, the
corresponding Hashin's bounds $\mu^{H+}$ and $\mu^{H-}$ may be conveniently written in the following
form:

$$
\frac{\mu^{H-}}{\mu_1} = 1 + f\varphi(0) \qquad \frac{\mu^{H+}}{\mu_1} = 1 + f\varphi(\infty)
$$

$$
\varphi(x) = \left( A(1-f) + \frac{1}{m_2/m_1 - 1} - \frac{f(1-f^{2/3})^2}{Bf^{7/3} + C + \frac{25(1-n_1)^2}{6(x-1)}} \right)^{-1}
$$

with

$$
A = \frac{2(4 - 5v_1)}{15(1 - v_1)}
$$

$$
B = \frac{10(1 - v_1)}{21} \frac{(7 - 10v_2)(7 + 5v_1)\mu_1 - (7 - 10v_1)(7 + 5v_2)\mu_2}{4(7 - 10v_2)\mu_1 + (7 + 5v_2)\mu_2} \quad (2.27)
$$

$$
C = \frac{10}{21} (7 - 10v_1)(1 - v_1)
$$

In these equations, $v_1$ and $v_2$ are the Poisson ratios of the considered phases.

One can see on Fig. 4 an illustration of these results for a particular case. Whereas it
can be checked that these bounds are closer to each other than Voigt-Reuss's ones, it may
also be noticed that they happen to be more distant than Hashin-Shtrikman's ones. This
remark may seem somewhat disappointing since Hashin-Shtrikman's bounds are known to
deal with unspecified (isotropic) morphologies, whereas the C.S.A. has especially been
defined in order to specify a matrix/inclusion-type morphology... We only note this apparent
inconsistency at the moment and shall go back to this point later.

<figure><img src="image_6.png" /><figcaption>Fig 4 Hashin's bounds for the overall shear modulus for the C.S.A., compared with Voigt-Reuss's and Hashin-Shtrikman's bounds ($\mu_2/\mu_1 = 6$, incompressibility)</figcaption></figure>

2.4.2.2 The three-phase model

Aiming at "solving" the problem of the isotropic C.S.A. for the determination of the
overall shear modulus, Christensen and Lo actually proposed an estimate of this (not
uniquely defined) modulus. This estimate was derived from the condition that, under given
homogeneous (stress or strain) loading conditions at infinity, the strain energy density in
some homogeneous isotropic elastic infinite body remains unchanged if part of this body is
replaced by a composite sphere of the C.S.A.: the elastic moduli of such a homogeneous
body would yield an estimate of $\mu^{\text{C.S.A.}}$ and $k^{\text{C.S.A.}}$. In other words, this so-called "three-
phase model" (see Fig. 5) deals with a composite sphere embedded in a "third phase", the
homogeneous equivalent medium (H.E.M.) whose moduli $\mu^{\text{3PM}}$ and $k^{\text{3PM}}$ are derived from
the foregoing energy condition.

Here again, the elementary inclusion problem is solved from classical Love results as a
function of the unknown parameters $\mu^{3PM}$ and $k^{3PM}$ which are determined from the energy
condition, conveniently transformed into a quite tractable equation by use of Eshelby's
integral [22]. It happens that the resulting value of $k^{3PM}$ coincides with Hashin's finding of
eqn (2.26) whereas $\mu^{3PM}/\mu_1$ is the positive solution of a second order equation. It can be
checked that this equation may be put in the following form, which refers to the notations of
formula (2.27):

$$
x = 1 + f \varphi(x) \tag{2.28}
$$

It is easy to verify that this estimate always lies between Hashin's lower and upper
bounds. Nevertheless, it is not quite clear whether the underlying model still respects the
wished matrix connectedness and so reflects the corresponding morphology. As a matter of fact, the treated configuration (Fig. 5) according to which a representative composite sphere has been embedded in a fictitious homogeneous medium does not express any longer the geometrical continuity of the matrix phase, at least as directly as initially done in Hashin's assemblage. So, the morphological meaning of the three-phase model would improve on clarification, which could be done by comparing this estimate of the C.S.A. elastic moduli to some others to be derived. Such a derivation can be performed easily by applying to the C.S.A. the method proposed in a simpler form in section 2.2: instead of dealing with simple ellipsoidal inclusions, we would just have now to consider composite spheres.

<figure><img src="image_7.png" /><figcaption>Fig. 5 The three-phase model</figcaption></figure>

### 2.4.2.3 A set of estimates for the C.S.A.

Let us consider the "phases" which constitute the C.S.A. in a new sense: instead of still referring to the two mechanical ones, namely the matrix and the particles, this term is applied now to composite spheres with the same diameter. Let ($\lambda$) be such a "phase" defined by the diameter $d_\lambda$, with the volume fraction $f_\lambda$: according to the proposed approach adapted to the present situation, its average strain $\bar{\epsilon}_\lambda$ and stress $\bar{\sigma}_\lambda$ state when the homogeneous strain $E$ (or stress $\Sigma$) is applied to the C.S.A. at infinity is estimated to be the same as the one of a composite sphere belonging to ($\lambda$) embedded in an infinite elastic matrix with moduli ($\mu^0, k^0$), under some prescribed uniform strain $E^0$ (or stress $\Sigma^0$) at infinity. The searched estimates $C_0^{est}$ are derived from

$$ \begin{gathered} \langle \bar{\epsilon} \rangle = \sum_{\lambda} f_{\lambda} \bar{\epsilon}_{\lambda} = E \\ (\text{or } \langle \bar{\sigma} \rangle = \sum_{\lambda} f_{\lambda} \bar{\sigma}_{\lambda} = \Sigma) \end{gathered} \qquad (2.29) \\ \langle \bar{\sigma} \rangle = C_0^{est} : \langle \bar{\epsilon} \rangle \end{gathered} $$

Obviously, $\bar{\epsilon}_\lambda$ and $\bar{\sigma}_\lambda$ do not depend on $d_\lambda$, so that $\bar{\epsilon}_\lambda = E$ (or $\bar{\sigma}_\lambda = \Sigma$). One calculation only is needed in order to derive $\bar{\epsilon}$ (or $\bar{\sigma}$) from $E^0$ (or $\Sigma^0$) for any composite sphere with the volume fraction $f$. This can be made easily from the basic Love's solutions.
If $\bar{e}, \bar{e}_1, \bar{e}_2, \mathcal{E}^0$ and $\bar{\theta}, \bar{\theta}_1, \bar{\theta}_2, \Theta^0$ denote the average strain deviators and traces in the composite inclusion, the spherical shell, the spherical core and at infinity respectively, one can write the solution in the form

$$ 
\begin{aligned}
\bar{e}_1 &= a_1 \mathcal{E}^0, \quad \bar{e}_2 = a_2 \mathcal{E}^0, \quad \bar{e} = [(1-f)a_1 + fa_2] \mathcal{E}^0 \\
\bar{\theta}_1 &= b_1 \Theta^0, \quad \bar{\theta}_2 = b_2 \Theta^0, \quad \bar{\theta} = [(1-f)b_1 + fb_2] \Theta^0,
\end{aligned}
\qquad (2.30) $$

with scalar values $a_1, a_2, b_1$ and $b_2$ depending on $\mu^0$ and $k^0$. We get formally

$$ 
\begin{aligned}
\frac{k_0^{est}}{k_1} &= 1 + \frac{f(k_2/k_1 - 1)b_2/b_1}{1 - f + f b_2/b_1} \\
\frac{\mu_0^{est}}{\mu_1} &= 1 + \frac{f(\mu_2/\mu_1 - 1)a_2/a_1}{1 - f + f a_2/a_1}.
\end{aligned}
\qquad (2.31) $$

It is found easily that $b_2/b_1$ does not depend on $\mu^0$ and $k^0$, so that $k_0^{est}$ has a unique value which coincides with Hashin's solution (2.26). For $a_2/a_1$, we find

$$ \frac{a_2}{a_1} = \frac{(1-f) \varphi(\mu^0/\mu_1)}{\mu_2/\mu_1 - 1 - f \varphi(\mu^0/\mu_1)} \qquad (2.32) $$

where $\varphi(x)$ has been defined in (2.27). From (2.31), (2.32) and (2.27) we get finally:

$$ \frac{\mu_0^{est}}{\mu_1} = 1 + f \varphi(\mu^0/\mu_1) \qquad (2.33) $$

<figure><img src="image_8.png" /><figcaption>Fig. 6 A set of estimates for the shear modulus of an isotropic C.S.A. for a given volume fraction (the inclined straight line is the first bisector of the positive quadrant)</figcaption></figure>
A graphic representation, similar to the one of Fig. 1, can be derived for a given f value: Fig. 6 shows such a plot of $\mu_0^{est}/\mu_1$ as a function of $\mu^0/\mu_1$ which yields a continuous set of estimates for the overall shear modulus of the isotropic C.S.A., lying between Hashin's lower and upper bounds. Several comments can be made from this plot:

- Hashin's bounds $\mu^H$ and $\mu^H_-$ for the shear modulus of the C.S.A. are two "bounds" for the present set of estimates: they correspond to $\mu^0 = 0$ and $\mu^0 \to \infty$ respectively. Referring to the single inclusion-based estimates of section 2.2 and to the analogy with classical Voigt-Reuss's bounds, one could say that Hashin's bounds are the transcription of Voigt-Reuss's ones to the specific morphology of the C.S.A.: as a matter of fact, Hashin's transfer of the homogeneous stress or strain boundary conditions of the R.V.E. to the boundary of every composite sphere is similar to Voigt-Reuss's procedure of prescribing uniform stress or strain conditions throughout the R.V.E. (i.e. of transferring the macroscopic boundary conditions to any point of the R.V.E.). In some sense to be made more precise, we could say that Hashin's bounds are the composite sphere-based extension of Voigt-Reuss's ones or that Voigt-Reuss's bounds are the "point" version of Hashin's ones.

- the three-phase model estimate is nothing but the "self-consistent" element of the present set, corresponding to $\mu^0 = \mu^{est}$, as already apparent from (2.28): this is a justification of the denomination of "generalized self-consistent scheme" (G.S.C.S.) which is given frequently to the three-phase model. Now again, we could say that the three-phase model is the composite sphere-based extension of the classical self-consistent scheme or that the classical self-consistent scheme is the "point" version of the generalized one. Similarly to the single inclusion-based estimates of section 2.2, we could define an iterative procedure and recursive formulae converging towards the three-phase model estimate: this would allow this model to appear as the proper treatment of a hypothetical "perfectly disordered" composite spheres assemblage. Incidentally, it can be checked that Christensen-Lo's "energy condition" coincides with a direct strain or stress condition ($E^0 = E$ or $\Sigma^0 = \Sigma$), a coincidence which can be proved to be valid in more general situations [23].

- going further along this speculative line of thinking, we can focus attention on two special estimates of the considered set, corresponding to $\mu^0 = \mu_1$ and $\mu^0 = \mu_2$ respectively. Note that, when $(k_1 - k_2)(\mu_1 - \mu_2) \ge 0$, the first one is nothing but the classical Hashin-Shtrikman's lower bound (or Mori-Tanaka's estimate) for the considered two-phase isotropic material. Their derivation by embedding the composite sphere in the softer and the stiffer phases looks like an extension of Walpole's interpretation of the classical Hashin-Shtrikman's bounds adapted to the C.S.A., so that one could guess them to coincide with new Hashin-Shtrikman-type bounds for the C.S.A.. If so, they would represent improved bounds, with respect to Hashin's ones, for the shear modulus of the C.S.A., without any new information except for a better use of the property of an isotropic spatial distribution of the C.S.A.. Anticipating the following we could say that such new bounds are the composite sphere-based extension of the classical Hashin-Shtrikman's ones or that Hashin-Shtrikman's bounds are the "point" version of these new ones.

Before demonstrating this guess in the next chapter and generalizing the corresponding "composite spheres approach" to arbitrary geometrical patterns, let us briefly mention that the foregoing analysis of the C.S.A. can be generalized easily both to n-layered isotropic spheres [24] and to transversely isotropic cylinders [25] assemblages, in view of applications to coated particle or fiber-reinforced composites, as well as to approximate treatments of nonlinear constitutive behaviours (see chapter 4).
### 3. MORPHOLOGICALLY REPRESENTATIVE PATTERN-BASED BOUNDING AND ESTIMATING IN ELASTICITY

#### 3.1 Introduction

We have already stressed the fact that most of the classical micromechanical approaches use a description by points of the phase spatial distribution: this is obvious for Kröner's systematic theory which uses intensively point-correlation functions of various orders for the elastic moduli; it is also the case for single inclusion-based estimates (one single inclusion in an infinite matrix has no internal microstructure so that it behaves like a point) as well as for the classical Voigt and Reuss or Hashin and Shtrikman bounds. In principle, point-correlation functions are a powerful tool for the description of complicated microstructure but they reveal themselves rather limited in practice: their experimental determination cannot be achieved beyond an order larger than 2 or 3 while the description of as elementary morphological properties as the connectedness of a matrix phase looks practically out of reach.

Hashin's composite spheres (or cylinders) assemblage and its by-products open a quite novel method of morphological description by using finite inhomogeneous elements (instead of points) as basic "bricks": first, the complete knowledge of the material content of such bricks includes a much richer information that usual; in addition, something must be said on the spatial distribution of the bricks, which still enriches the morphological description of the material (one must keep in mind that both informations, namely the constitution of a composite sphere and the spatial distribution of the spheres, were necessary to make the C.S.A. express the matrix continuity).

Looking back upon Hashin's approach, we could find his definition of the "bricks" rather sufficient for his purpose but his treatment of their spatial distribution somewhat improvable (no explicit used was made of the implicitly assumed isotropy of this distribution). In addition, one family of composite spheres was considered only and the need to fill up the whole space with this spheres made necessary to use an artificial quasi-fractal process, so that no local fluctuation of the particle volume fractions was allowed.

In what follows, we shall aim at saving the basic idea of dealing with finite composite elements, without initial restriction to spheres or cylinders, and at improving the geometrical description of their spatial distribution. By doing so, we want to combine a deterministic description of small, but finite, well-chosen "composite patterns" and a statistical representation of their distribution. Such a morphological representation will then be used to adapt the classical Hashin-Shtrikman's variational procedure in order to yield bounds for the overall moduli of a heterogeneous material [5]. Moreover, this "pattern approach" will be applied to the definition of new estimates for these moduli, including a generalized "morphologically pattern-based self-consistent scheme" [7] which uses numerical computations for the resolution of the basic underlying composite inclusions problems.

#### 3.2 A variational pattern-based approach to heterogeneous elasticity

##### 3.2.1 Simplified approach

Let us first consider a multiphase E.V.R., whose extension can be considered as infinite, submitted to homogeneous strain conditions *E* at infinity; the volume *V* of this domain *Ω* is supposed to have been entirely decomposed into several "morphological phases" (λ) with the volume fraction $c_λ = V_λ/V$, consisting of *N*λ identical composite domains $D_{λ,i}$ centered at $\bar{X}_{λ,i}$, translated from the same domain $D_λ$ centered at the origin.
This domain $D_\lambda$ will be considered as the "morphologically representative pattern" (say M.R.P.) of the $(\lambda)$ phase (Fig. 7).

<figure><img src="image_9.png" /><figcaption>Fig. 7 A morphologically representative pattern</figcaption></figure>

Note that any inhomogeneous material may be described in this manner, at least with as many patterns, reducing to points, as mechanical phases. The classical "point approach" will then be a special case of the present one and classical bounds have to be recovered as by-products. In addition, such a pattern decomposition is not unique for a given material: some descriptions may contain more informations than others and lead to more accurate bounds or estimates. In the same way, for one given pattern description, some variational approaches may be more powerful than others if they take more available information into account. Let us start with the simplest one, namely the Voigt-Reuss-type approach, before adopting the more efficient Hashin-Shtrikman variational procedure.

### 3.2.1.1 M.R.P.-based Voigt-Reuss bounds

Let us transfer the homogeneous strain boundary conditions at the boundary of each pattern, namely

$$ \bar{u}'_\lambda(\bar{x}) = E \cdot \bar{x}, \quad \forall \bar{x} \in \partial D_{\lambda_i}, \quad \forall D_{\lambda_i} \subset \Omega \qquad (3.1) $$

If one can compute the corresponding strain field $\varepsilon'_\lambda(\bar{x})$ inside each pattern $D_\lambda$, one has at disposal, by simply collecting these fields throughout the whole body $\Omega$, a clearly kinematically admissible strain field $\varepsilon'(\bar{x})$ which can be used as a Voigt-type trial field in the potential energy theorem. This leads to a "M.R.P.-based Voigt bound" $C_{MRP}^V$ [6] such that

$$ C_{MRP}^V \ge C^{\text{eff}}, \quad C_{MRP}^V: \langle \varepsilon' \rangle = \langle c: \varepsilon' \rangle = \sum_{\lambda} \frac{c_{\lambda}}{D_{\lambda}} \int_{D_{\lambda}} c_{\lambda}: \varepsilon'_{\lambda} d\omega, \quad \forall E. \qquad (3.2) $$

Similarly, a "M.R.P.-based Reuss bound" is obtained by applying homogeneous stress conditions $\Sigma^*$ at the boundary of each pattern and by computing the associated stress field $\sigma^*_\lambda(\bar{x})$ inside the patterns such that
$$ \sigma_{\lambda}^{*}(\bar{x}).\bar{n}_{\lambda}(\bar{x}) = \Sigma^{*}.\bar{n}_{\lambda}(\bar{x}), \quad \forall \bar{x} \in \partial D_{\lambda_i}, \quad \forall D_{\lambda_i} \subset \Omega \quad (3.3) $$

with $\bar{n}_\lambda(\bar{x})$ denoting the outward unit normal vector to $\partial D_{\lambda_i}$. The optimal choice $\Sigma^*$ for a given $E$ leads to the bound $S_{MRP}^R$ such that

$$ (S_{MRP}^R)^{-1} \le C^{\text{eff}}, \quad S_{MRP}^R: \langle \sigma^* \rangle = \langle s: \sigma^* \rangle = \sum_{\lambda} \frac{c_{\lambda}}{D_{\lambda}} \int_{D_{\lambda}} s_{\lambda} : \sigma_{\lambda}^* d\omega, \quad \forall \Sigma \quad (3.4) $$

It can be checked that $C_{MRP}^V$ and $S_{MRP}^R$ have the right required properties of symmetry and positiveness. In the general case, the computation of the 21 components of these tensors requires the (numerical) resolution of as many problems as different patterns under six linearly independent homogeneous applied strain and stress tensors. The simplest case has been solved in closed form by Hashin for the isotropic C.S.A. (see § 2.4.2.1.) and by Hashin and Rosen [26] for the transversely isotropic C.C.A.

### 3.2.1.2 M.R.P.-based Hashin-Shtrikman bounds

The Hashin-Shtrikman variational procedure can be applied to this pattern approach too. It makes use of trial fields derived from the application of an arbitrary polarization stress field $p(\bar{y})$ to a homogeneous reference medium with the elastic moduli $C^0$ but with the same geometry and boundary conditions as the considered inhomogeneous body. This polarization field of symmetric second order tensors is defined throughout the body as

$$ \mathbf{p}(\bar{y}) = \sum_{\lambda} \mathbf{p}_{\lambda}(\bar{y})\phi_{\lambda}(\bar{y}) \quad (3.5) $$

where $\phi_\lambda(\bar{y})$ is the characteristic function of $V_\lambda$. As for $\mathbf{p}_\lambda(\bar{y})$, we choose it as translated by $\bar{X}_{\lambda_i}$ from some $\mathbf{P}_\lambda(\bar{x})$ field to be defined on $D_\lambda$, so that we get

$$ \mathbf{p}_{\lambda}(\bar{x} + \bar{X}_{\lambda_i}) = \mathbf{P}_{\lambda}(\bar{x}), \quad \bar{x} \in D_{\lambda_i}, \quad \forall \lambda, \forall i. \quad (3.6) $$

<figure><img src="image_10.png" /><figcaption>Fig. 8 Homologous points of a M.R.P.</figcaption></figure>
In other words, we prescribe identical values of the polarization stress field at "homologous points", i.e. at points of the same morphological phase with the same relative position inside the phase domains they belong to (Fig.8).

We start from the classical expression of the Hashin-Shtrikman functional $HS^0(p)$ which reads

$$2HS^0(p) = E: C^0: E + E: <p> + <\epsilon>: p> - <p>: H^0: p> \qquad (3.7)$$
$$\text{with } H^0 = (\delta c^0)^{-1} = (c - C^0)^{-1}$$

where $\epsilon'$ is the (kinematically admissible) strain field deriving from $p$. From (3.6) we get:

$$2HS^0(p) = E: C^0: E + \sum_{\lambda} \frac{N_{\lambda}}{V} \int_{D_{\lambda}} (E + \epsilon_{\lambda}^{M} - P_{\lambda}: H_{\lambda}^{0}): P_{\lambda} \, d\omega \qquad (3.8)$$

where we have used for the strain field $\epsilon'$ the following definition of "pattern-averaged" fields $f(\bar{x})$:

$$f_{\lambda}^{M}(\bar{x}) = \frac{1}{N_{\lambda}} \sum_{i=1}^{N_{\lambda}} f(\bar{x} + \bar{X}_{\lambda_i}) \qquad (3.9)$$

From (2.13), we can derive $\epsilon'$ from $p$ and then $\epsilon_{\lambda}^{M}(\bar{x})$:

$$\begin{align*}
\epsilon_{\lambda}^{M}(\bar{x}) &= \frac{1}{N_{\lambda}} \sum_{i=1}^{N_{\lambda}} \epsilon'(\bar{x} + \bar{X}_{\lambda_i}) \\
&= \epsilon^0 - \frac{1}{N_{\lambda}} \sum_{i=1}^{N_{\lambda}} \int \Gamma^0(\bar{x} + \bar{X}_{\lambda_i}, \bar{x}'): p(\bar{x}') d\omega' \\
&= \epsilon^0 - \int_{V} \Gamma_{\lambda}^{0M}(\bar{x}, \bar{x}'): p(\bar{x}') d\omega' \\
&= \epsilon^0 - \sum_{\mu} \int_{V_{\mu}} \Gamma_{\lambda}^{0M}(\bar{x}, \bar{x}'): p_{\mu}(\bar{x}') d\omega' \\
&= \epsilon^0 - \sum_{\mu} \sum_{j=1}^{N_{\mu}} \int_{D_{\mu}} \Gamma_{\lambda}^{0M}(\bar{x}, \bar{x}'+\bar{X}_{\mu_j}): P_{\mu}(\bar{x}') d\omega' \tag{3.10} \\
&= \epsilon^0 - \sum_{\mu} N_{\mu} \int_{D_{\mu}} \Gamma_{\lambda\mu}^{0MM}(\bar{x}, \bar{x}'): P_{\mu}(\bar{x}') d\omega'
\end{align*}$$

In this equation, (3.9) has been applied twice to the Green operator (once for each variable $\bar{x}$ and $\bar{x}'$). The auxiliary uniform strain $\epsilon^0$ has to be determined through the average property:

$$\langle \epsilon' \rangle = \sum_{\lambda} \frac{N_{\lambda}}{V} \int_{D_{\lambda}} \epsilon_{\lambda}^{M}(\bar{x}) d\omega = \sum_{\lambda} \frac{c_{\lambda}}{D_{\lambda}} \int_{D_{\lambda}} \epsilon_{\lambda}^{M}(\bar{x}) d\omega = E \qquad (3.11)$$

so that we have:
$$ \sum_{\lambda} \frac{N_{\lambda}}{V} \int_{D_{\lambda}} [\varepsilon^{0} - \sum_{\mu} N_{\mu} \int_{D_{\mu}} \Gamma_{\lambda\mu}^{0MM}(\bar{x}, \bar{x}') : \mathbf{P}_{\mu}(\bar{x}') d\omega'] d\omega = \mathbf{E} \qquad (3.12) $$

$$ \text{or } \varepsilon^0 = \mathbf{E} + \sum_{\mu} \sum_{\lambda} \frac{N_{\lambda} N_{\mu}}{V} \int_{D_{\lambda}} \int_{D_{\mu}} \Gamma_{\lambda\mu}^{0MM}(\bar{x}, \bar{x}') : \mathbf{P}_{\mu}(\bar{x}') d\omega' d\omega $$

Note that the difference between $\varepsilon^0$ and $\mathbf{E}$ stems from the infinite extension of the body and from the corresponding definition of the Green operator in that case (an equivalent alternative treatment would have been to deal with polarization fields deviations, i.e. fields with null average).

After integration of (3.10) into (3.8), we get:

$$ \begin{aligned} 2HS^0(\mathbf{P}_\lambda) &= \mathbf{E} : \mathbf{C}^0 : \mathbf{E} + \sum_{\lambda} \frac{N_\lambda}{V} \int_{D_\lambda} \mathbf{P}_\lambda : (\mathbf{E} + \varepsilon_\lambda^{\prime M} - \mathbf{H}_\lambda^0 : \mathbf{P}_\lambda) d\omega \\ &= \mathbf{E} : \mathbf{C}^0 : \mathbf{E} + \sum_{\lambda} \frac{N_\lambda}{V} \int_{D_\lambda} \mathbf{P}_\lambda : (\mathbf{E} + \varepsilon^0 - \mathbf{H}_\lambda^0 : \mathbf{P}_\lambda) d\omega \dots \end{aligned} \qquad (3.13) $$

$$ \dots - \sum_{\mu} \sum_{\lambda} \frac{N_{\lambda} N_{\mu}}{V} \int_{D_{\lambda}} \int_{D_{\mu}} \mathbf{P}_{\lambda}(\bar{x}) : \int_{D_{\mu}} \Gamma_{\lambda\mu}^{0MM}(\bar{x}, \bar{x}') : \mathbf{P}_{\mu}(\bar{x}') d\omega' d\omega. $$

We know that if $\mathbf{C}^0$ is such that $\mathbf{H}_\lambda^0$ is positive or negative everywhere, then $HS^0(\mathbf{p})$ yields bounds for the overall moduli. Thus, the polarization fields $\mathbf{P}_\lambda$ have to be optimized in order to get optimal bounds. The stationarity of the functional $HS^0$ is obtained by:

$$ \frac{\partial HS^0}{\partial \mathbf{p}} : \delta \mathbf{p} = 0, \quad \forall \delta \mathbf{P}_\lambda. \qquad (3.14) $$

From (3.13) and taking account of (3.12) for the dependence of $\varepsilon^0$ on $\mathbf{P}_\lambda$, we get, after some elementary algebra:

$$ \frac{\partial HS^0}{\partial \mathbf{p}} : \delta \mathbf{p} = \sum_{\lambda} \frac{N_{\lambda}}{V} \int_{D_{\lambda}} \delta \mathbf{P}_{\lambda}(\bar{x}) : [\varepsilon_{\lambda}^{\prime M}(\bar{x}) - \mathbf{H}_{\lambda}^{0}(\bar{x}) : \mathbf{P}_{\lambda}(\bar{x})] d\omega = 0. \qquad (3.15) $$

This condition can be satisfied, due to the fact that the patterns are disjoint and the variables $\mathbf{P}_\lambda$ are independent on each other, by

$$ \varepsilon_{\lambda}^{\prime M}(\bar{x}) = \mathbf{H}_{\lambda}^{0}(\bar{x}) : \mathbf{P}_{\lambda}(\bar{x}) \qquad (3.16) $$

$$ \text{or } \mathbf{P}_{\lambda}(\bar{x}) = \delta \mathbf{c}_{\lambda}^{0}(\bar{x}) : \varepsilon_{\lambda}^{\prime M}(\bar{x}). $$

In other words, $\varepsilon_{\lambda}^{\prime M}(\bar{x})$ and $\mathbf{P}_{\lambda}(\bar{x})$ are the solutions of the following integral equations:

$$ \varepsilon_{\lambda}^{\prime M}(\bar{x}) = \varepsilon^{0} - \sum_{\mu} N_{\mu} \int_{D_{\mu}} \Gamma_{\lambda\mu}^{0MM}(\bar{x}, \bar{x}') : \delta \mathbf{c}_{\mu}^{0}(\bar{x}') : \varepsilon_{\mu}^{\prime M}(\bar{x}') d\omega' \qquad (3.17) $$

with the condition (3.11), or
$$
\begin{equation}
\begin{aligned}
\mathbf{H}_{\lambda}^{0}(\bar{x}): \mathbf{P}_{\lambda}(\bar{x}) &= \mathbf{E} + \sum_{\mu} \sum_{\lambda} \frac{N_{\lambda} N_{\mu}}{V} \int_{D_{\lambda}} \int_{D_{\mu}} \Gamma_{\lambda\mu}^{0MM}(\bar{x}, \bar{x}'): \mathbf{P}_{\mu}(\bar{x}') d\omega' d\omega \dots \\
&\quad - \sum_{\mu} N_{\mu} \int_{D_{\mu}} \Gamma_{\lambda\mu}^{0MM}(\bar{x}, \bar{x}'): \mathbf{P}_{\mu}(\bar{x}') d\omega'.
\end{aligned}
\tag{3.18}
\end{equation}
$$

Referring to (2.9), we can guess from (3.17) that $\epsilon_{\lambda}^{MM}(\bar{x})$ could be the solution of an
appropriate problem of interacting inhomogeneities $D_{\mu}$ in an infinite matrix, but such an
interpretation would need additional informations on the spatial distribution of the patterns to
be made more apparent.

3.2.1.3 The case of isotropy

We consider here an isotropic distribution of patterns in the following sense: let
$\Psi_{\lambda\mu}(\bar{u})d\omega_u$ be the number of couples of domains $D_{\lambda_i}$ and $D_{\mu_j}$ such that their centers $\bar{X}_{\lambda_i}$ and $\bar{X}_{\mu_j}$ lie at a distance $\bar{u}$ within $d\omega_u$. The assumption of an isotropic distribution of patterns then implies that $\Psi_{\lambda\mu}(\bar{u})$ depends on $u = \|\bar{u}\|$. In addition, this implies too that the external shape of the patterns must be spherical: consider two patterns $\lambda$ and $\mu$ such that they contain at least two subdomains $D_{\lambda_i}$ and $D_{\mu_j}$ in contact with each other (such a situation happens necessarily since the domains have to map the whole space). As two subdomains cannot lie at the same place, $\Psi_{\lambda\mu}(u)$ must vanish over a finite neighbourhood of 0 (except for 0 itself for $\lambda = \mu$). Let $u_{min} = \|\bar{u}_{min}\|$ be the minimal (non zero) distance between the centers of (thus contiguous) subdomains $\lambda$ and $\mu$: due to the isotropic distribution of the centers, this concerns couples of domains with an equiprobable orientation of $\bar{u}_{min}$ and it is easy to prove that such a situation cannot occur unless the subdomains be spherical.

Thus, taking account also of the translation invariance of the Green operator in an
infinite medium, we can write

$$
\begin{align}
& N_{\mu} \int_{D_{\mu}} \Gamma_{\lambda\mu}^{0MM}(\bar{x}, \bar{x}'): \mathbf{P}_{\mu}(\bar{x}') d\omega' = N_{\mu} \int_{D_{\mu}} \Gamma_{\lambda\mu}^{0MM}(\bar{x} - \bar{x}'): \mathbf{P}_{\mu}(\bar{x}') d\omega' \nonumber \\
& \quad = \frac{1}{N_{\lambda}} \sum_{i=1}^{N_{\lambda}} \sum_{j=1}^{N_{\mu}} \int \Gamma^0(\bar{X}_{\lambda_i} - \bar{X}_{\mu_j} + \bar{x} - \bar{x}'): \mathbf{P}_{\mu}(\bar{x}') d\omega' \tag{3.19} \\
& \quad = \frac{1}{N_{\lambda}} \int d\omega_u \int_{D_{\mu}} \Psi_{\lambda\mu}(u) \Gamma^0(\bar{u} + \bar{x} - \bar{x}'): \mathbf{P}_{\mu}(\bar{x}') d\omega'. \nonumber
\end{align}
$$

In this integral, $\bar{x}$ and $\bar{x}'$ belong to the spheres $D_\lambda$, with the radius $r_\lambda$, and $D_\mu$, with the radius $r_\mu$, respectively so that $\bar{x} - \bar{x}'$ belongs to a sphere with a radius ($r_\lambda + r_\mu$). With use of Fubini's theorem, the integration may be performed first on $u$ over spherical concentric shells with infinitesimal width $du$ where $\Psi_{\lambda\mu}(u)$ is constant. These terms vanish for $u \ge r_\lambda + r_\mu$ due to the well-known property of the Green operator which is responsible for the uniformity of the solution in the ellipsoidal inclusion of Eshelby's problem; they do so too
for $u \le r_\lambda + r_\mu$ when $\lambda \neq \mu$ since $\Psi_{\lambda\mu}(u)$ vanishes in this case (the distance between the centers of two spheres belonging to two different patterns cannot be smaller than the sum of their radii). A non zero contribution comes only from $\lambda = \mu$ and $u = 0$ (with $\Psi_{\lambda\lambda}(0) = N_\lambda$) so that we are led to the result: $N_\mu \Gamma_{\lambda\mu}^{0MM}(\bar{x}, \bar{x}') = \delta_{\lambda\mu} \Gamma^0(\bar{x}, \bar{x}')$, $\forall \bar{x} \in D_\lambda, \bar{x}' \in D_\mu$.

We get finally

$$ 
\begin{aligned}
N_\mu \int_{D_\mu} \Gamma_{\lambda\mu}^{0MM}(\bar{x}, \bar{x}'): \mathbf{P}_\mu(\bar{x}') d\omega' &= \delta_{\lambda\mu} \int_{D_\mu} \Gamma^0(\bar{x} - \bar{x}'): \mathbf{P}_\mu(\bar{x}') d\omega' \\
&= \int_{D_\lambda} \Gamma^0(\bar{x} - \bar{x}'): \mathbf{P}_\lambda(\bar{x}') d\omega', \quad \forall \bar{x} \in D_\lambda
\end{aligned} 
\quad (3.20) 
$$

so that (3.10) reads now:

$$ \epsilon_{\lambda}^{M}(\bar{x}) = \epsilon^{0} - \int_{D_{\lambda}} \Gamma^{0}(\bar{x} - \bar{x}'): \mathbf{P}_{\lambda}(\bar{x}') d\omega', \quad \forall \bar{x} \in D_{\lambda} \quad (3.21) $$

and the integral equation (3.17) becomes:

$$ \epsilon_{\lambda}^{M}(\bar{x}) = \epsilon^{0} - \int_{D_{\lambda}} \Gamma^{0}(\bar{x}, \bar{x}'): \delta c_{\lambda}^{0}(\bar{x}'): \epsilon_{\lambda}^{M}(\bar{x}') d\omega' \quad (3.22) $$

with the condition (3.11) to be used for the determination of $\epsilon^0$. This proves that, due to the assumption of an isotropic distributions of the patterns, the pattern-averaged strain trial fields $\epsilon_{\lambda}^{M}(\bar{x})$ which can lead to optimal bounds for such an assumption can be derived as the solution of elementary problems of spherical composite inhomogeneities in an infinite homogeneous matrix (Fig. 9).

<figure><img src="image_11.png" /><figcaption>Fig. 9 Interpretation of M.R.P.-based Hashin-Shtrikman bounds in terms of composite inclusions problems</figcaption></figure>

The elastic moduli $C^0$ of the matrix have to be optimized too in order to make the condition that $H_\lambda^0$ be positive or negative everywhere "just" satisfied: roughly speaking, this
means that these moduli must be those of the softest or the stiffest mechanical phases belonging to the considered composite. This allows finally to derive optimal expressions $HS^{opt}$ of the Hashin-Shtrikman functional and then Hashin-Shtrikman-type bounds $C_{MRP}^{HS}$

through: $HS^{opt} = \frac{1}{2} E : C_{MRP}^{HS} : E$.

Note that in the foregoing analysis nothing has been specified about the mechanical anisotropy of the constituents and the material content of the patterns: the analysis is still valid when the spherical patterns exhibit no internal spherical symmetry or when elasticity is locally arbitrarily anisotropic so that, despite the isotropic distribution of the centers of the patterns, the overall behaviour may be anisotropic (see examples of such situations on Fig. 10). Note also that, in such cases, the Green operator need not be known in closed form: as a matter of fact, we only have to solve composite inclusions / matrix elastic problems, which can be performed by other techniques (including numerical, e.g. F.E.M., ones) than Green's ones.

<figure><img src="image_12.png" /><figcaption>Fig. 10 Examples of isotropic distribution of patterns</figcaption></figure>

### 3.2.1.4 Going back to the isotropic C.S.A.

A straightforward application of what precedes concerns the isotropic C.S.A. [27]. We can consider the "morphological phases ($\lambda$) as consisting of the composite spheres $S_\lambda$ with given external and internal radii $b_\lambda$ and $a_\lambda$ respectively, satisfying $f = (a_\lambda/b_\lambda)^3$. The optimal value of the Hashin-Shtrikman functional is, from (3.16) and (3.8):

$$HS^0(\varepsilon_{\lambda}^{M}) = \frac{1}{2} E : C^0 : E + \frac{1}{2V} \sum_{\lambda} N_{\lambda} E : \int_{S_{\lambda}} \delta c^0(\bar{x}) : \varepsilon_{\lambda}^{M}(\bar{x}) d\omega \quad (3.23)$$

with $\varepsilon_{\lambda}^{M}(\bar{x})$ given by (3.22) and (3.11) and $\delta c^0(\bar{x})$ positive or negative everywhere. If $\bar{\varepsilon}_{\lambda}^{(2)}$ and $\bar{\varepsilon}_{\lambda}^{(1)}$ are respectively the average strain tensors in the core and the shell of $S_\lambda$ when
embedded in an infinite matrix with moduli Cº, submitted to eº at infinity, we can write
(3.23) as:

$$
HS^0(\epsilon_\lambda^{\prime M}) = \frac{E: C^0: E}{2} + \frac{E}{2} \sum_\lambda c_\lambda [f(c_2 - C^0): \bar{\epsilon}_\lambda^{(2)} + (1-f)(c_1 - C^0): \bar{\epsilon}_\lambda^{(1)}]. \quad (3.24)
$$

But we know that $\bar{\epsilon}_{\lambda}^{(2)}$ and $\bar{\epsilon}_{\lambda}^{(1)}$ do not depend on the size of $S_{\lambda}$, so that we have:

$$
HS^0(\bar{\epsilon}_1, \bar{\epsilon}_2) = \frac{1}{2} E : C^0 : E + \frac{1}{2} E : [f(c_2 - C^0) : \bar{\epsilon}_2 + (1-f)(c_1 - C^0) : \bar{\epsilon}_1] \\
\text{with } E = f\bar{\epsilon}_2 + (1-f)\bar{\epsilon}_1, \quad \bar{\epsilon}_\lambda^{(2)} = \bar{\epsilon}_2, \quad \bar{\epsilon}_\lambda^{(1)} = \bar{\epsilon}_1, \quad \forall \lambda \quad (3.25)
$$

and then
$$ HS^0(\bar{\epsilon}_1, \bar{\epsilon}_2) = \frac{1}{2} E : [f c_2 : \bar{\epsilon}_2 + (1-f) c_1 : \bar{\epsilon}_1] . $$

<figure><img src="image_13.png" /><figcaption>Fig 11 Comparison between the bounds derived, for the overall shear modulus of an isotropic C.S.A., according to Hashin's treatment (continuous lines) or to the present one (dotted lines): phase 2 included in phase 1, μ₂/μ₁ = 6, (a) incompressibility and (b) ν₁ = .45 and ν₂ = 0. (from [27])</figcaption></figure>
With the notations of (2.30) and with $\mathcal{E}$ and $\Theta$ the deviator and trace of $\mathcal{E}$, we get easily:

$$HS^0 = \frac{f\mu_2 a_2 + (1-f)\mu_1 a_1}{fa_2 + (1-f)a_1} \mathcal{E} : \mathcal{E} + \frac{fk_2 b_2 + (1-f)k_1 b_1 \Theta^2}{fb_2 + (1-f)b_1} \frac{2}{2} \quad (3.26)$$

This leads directly to prove that the two associated bounds for $k^{\text{C.S.A.}}$ coincide, with the same value as Hashin's one, whereas the bounds for $\mu^{\text{C.S.A.}}$, say $\mu^{\text{HS}^+}$ and $\mu^{\text{HS}^-}$, are given, with the notations of § 2.4.2.3, by:

$$\frac{\mu^{\text{HS}^+}}{\mu_1} = 1 + f \varphi(\mu_2/\mu_1), \quad \frac{\mu^{\text{HS}^-}}{\mu_1} = 1 + f \varphi(1). \quad (3.27)$$

This result is a demonstration of the guess proposed in § 2.4.2.3. It allows to bound the overall shear modulus of the C.S.A. more efficiently than initially proposed by Hashin, as illustrated on Fig. 11. We can add that, provided that the assumption of an isotropic distribution of the pattern centers remains valid, several families of composite spheres, representing either local fluctuations of particle volume fractions or different kinds of particles, could be considered as well; especially, one family of spheres can be constituted with the matrix material only. Similarly, we can deal with several patterns of multilayered spherical inclusions. In addition, the transcription of this approach to composite (possibly coated) cylinders with transverse isotropy is quite straightforward.

We can also guess that the transcription of Kröner's "graded disorder" theory to the present pattern approach could allow us to prove that the three-phase model could be the adequate treatment of a "perfectly disordered" C.S.A. in a way similar to the one which allowed Kröner to prove, within the point approach, that the classical self-consistent scheme is the adequate treatment of a "perfectly disordered" medium. Nevertheless, the question of a proper definition - or even of the statistical feasibility - of such a "perfectly disordered" C.S.A. remains an open one...

3.2.1.5 Case of isotropic pattern distribution and cubic overall symmetry:

<figure><img src="image_14.png" /><figcaption>Fig. 12 Scheme of a composite with isotropically distributed cubic aligned inclusions</figcaption></figure>
The above treatment of the case of an isotropic distribution of the pattern centers can also be applied to situations of macroscopic anisotropy, depending on the material content of the patterns. Such a situation can be achieved for instance by a composite with cubic aligned particles distributed isotropically (Fig. 12). Both phases are considered as elastically isotropic so that the overall behaviour has a cubic symmetry: let $k^{eff}$, $\mu_a^{eff}$ and $\mu_d^{eff}$ the overall bulk and shear moduli respectively, with the definition

$$2\mu_a^{eff} = \frac{\langle \sigma_{xx} - \sigma_{yy} \rangle}{\langle \epsilon_{xx} - \epsilon_{yy} \rangle}, \quad 2\mu_d^{eff} = \frac{\langle \sigma_{xy} \rangle}{\langle \epsilon_{xy} \rangle} \qquad (3.28)$$

where x, y and z are the directions of the cubic symmetry.

Of course the M.R.P.-based bounds for these moduli cannot be derived in closed form anymore but numerical computations aiming at the solution of the basic matrix/inclusion underlying problems can be performed as well. We report here the results of such computations [6] for both shear moduli of (3.28) performed by a 3D F.E.M. technique using the mesh represented on Fig. 13 (351 20-noded quadratic elements and 1528 nodes). The M.R.P.-based Voigt-Reuss and Hashin-Shtrikman bounds have been computed by using moduli of the reference medium either far higher and lower than those of the constitutive phases for the former or equal to those of each of these phases for the latter. Meaningful comparisons may be made especially with the general bounds derived by Milton and Kohn [28] for cubic symmetry in order to appreciate the advantage of specifying a matrix/inclusion-type morphology.

<figure><img src="image_15.png" /><figcaption>Fig. 13 Typical 3D mesh used for the solution of the auxiliary inclusion problem corresponding to Fig 12 (from [6])</figcaption></figure>
In the present case ($k_2 > k_1$ and $\mu_2 > \mu_1$), Milton and Kohn's bounds read:

$$f_1\left(\frac{2}{\mu_2 - \mu_a^{\text{eff}}} + \frac{3}{\mu_2 - \mu_d^{\text{eff}}}\right) \le \frac{5}{\mu_2 - \mu_1} - \frac{6f_2(k_2 + 2\mu_2)}{\mu_2(3k_2 + 4\mu_2)} \quad (3.29)$$

$$f_2\left(\frac{2}{\mu_a^{\text{eff}} - \mu_1} + \frac{3}{\mu_d^{\text{eff}} - \mu_1}\right) \le \frac{5}{\mu_2 - \mu_1} - \frac{6f_1(k_1 + 2\mu_1)}{\mu_2(3k_1 + 4\mu_1)}$$

It can be checked on Fig. 14 that the specification of a matrix/inclusion morphology with isotropically distributed aligned cubic inclusions leads to a drastic tightening of the bounds for the shear moduli. Note that, in accordance with mechanical intuition, the M.R.P.-based bounds for this particle-reinforced composite lie close to the lower Milton-Kohn's one; nevertheless, the lower bounds do not coincide as they do in case of isotropy and "well-ordered" materials (i.e. with $k_2 > k_1$ and $\mu_2 > \mu_1$).

<figure><img src="image_16.png" /><figcaption>Fig. 14 Comparison between Milton-Kohn's and M.R.P.-based Voigt-Reuss and Hashin-Shtrikman's bounds for the shear moduli of a composite with isotropically distributed aligned cubic inclusions $k_1 = 3, \mu_1 = 1, k_2 = 20, \mu_2 = 10$ (arbitrary units) (2) in (1) (from [6])</figcaption></figure>

In what follows, we shall now explore some routes of generalization of the proposed approach, especially in view of dealing with anisotropic distributions as well as of building up new tools for an improved estimation of the effective behaviour.

### 3.2.2 Generalization

The foregoing analysis has been extended according to two different directions: the first one is concerned with the possibility to deal with more general definitions of patterns, which would not need anymore to be identical within a given family nor to fill up the whole space; the second one aims at taking better into account the overall anisotropy, especially
anisotropic distributions of patterns. These two kinds of extensions are briefly reported in what follows within the framework of linear elasticity.

### 3.2.2.1 An extended definition of patterns

Instead of restricting ourselves to morphological phases both composed of identical patterns and forced to map the whole space, we can consider more general situations without sensible increase of complexity [5]. First we allow any morphological phase ($\lambda$) to be constituted of patterns which may be slightly different, both geometrically and mechanically. Let $S_\lambda$ be the geometrical intersection of such patterns $D_{\lambda_k}$ after translation of their center $\bar{X}_{\lambda_k}$ to the origin: the foregoing definition of the polarization stress fields will be unchanged inside one part of the patterns, namely the translated of $S_\lambda$, but modified outside. This new definition will read:

$$ 
\begin{aligned}
\mathbf{p}_{\lambda}(\bar{X}_{\lambda_k} + \bar{x}) &= \mathbf{P}_{\lambda}(\bar{x}) \quad \text{if} \quad \bar{X}_{\lambda_k} + \bar{x} \in (\lambda), \bar{x} \in S_{\lambda} \\
\mathbf{p}_{\lambda}(\bar{X}_{\lambda_k} + \bar{x}) &= \mathbf{p}_{\lambda_0} \quad \text{if} \quad \bar{X}_{\lambda_k} + \bar{x} \in (\lambda), \bar{x} \notin S_{\lambda}
\end{aligned}
\quad (3.30) $$

where $\mathbf{p}_{\lambda_0}$ is a constant symmetrical second order tensor. In addition, the phases ($\lambda$) may not fill up the whole space: we can consider this part of $\Omega$ where no pattern information is available as an additional phase, say ($\nu$), with $S_\nu = \emptyset$.

Referring to the expression of the Hashin-Shtrikman functional, it appears that we need an additional definition of "pattern-average" fields $f(\bar{x})$, namely:

$$ f_{\lambda}^{m} = \frac{1}{V_{\lambda_0} v_{\lambda_0}} \int f(\bar{x}) d\omega \quad (3.31) $$

where $V_{\lambda_0}$ is that part of ($\lambda$) whose points obey the conditions: $\bar{X}_{\lambda_k} + \bar{x} \in (\lambda)$, $\bar{x} \notin S_{\lambda}$. This allows us to write a new version of (3.8), namely:

$$ 
\begin{aligned}
2HS^0(\mathbf{p}) = E:C^0:E + \sum_{\lambda} \frac{N_{\lambda}}{V} \int_{S_{\lambda}} (E + \varepsilon_{\lambda}^{M} - \mathbf{P}_{\lambda}:\mathbf{H}_{\lambda}^{0M}):\mathbf{P}_{\lambda} \, d\omega \dots \\
\qquad + \sum_{\lambda} \frac{V_{\lambda_0}}{V} (E + \varepsilon_{\lambda}^{m} - \mathbf{p}_{\lambda_0}:\mathbf{H}_{\lambda}^{0m}):\mathbf{p}_{\lambda_0}
\end{aligned}
\quad (3.32) $$

Instead of a deterministic $\mathbf{H}_{\lambda}^{0}(\bar{x})$ tensor with known identical values at homologous points $\bar{x}$, we have now to deal with pattern-averages $\mathbf{H}_{\lambda}^{0M}(\bar{x})$ and $\mathbf{H}_{\lambda}^{0m}$, which have to be derivable from some statistical knowledge about the material content of the patterns $D_{\lambda_k}$.

Eqns (3.10) and (3.11) now read:

$$ 
\begin{aligned}
\varepsilon_{\lambda}^{M}(\bar{x}) &= \varepsilon^0 - \sum_{\mu} N_{\mu} \int_{S_{\mu}} \Gamma_{\lambda\mu}^{0MM}(\bar{x}, \bar{x}') : \mathbf{P}_{\mu}(\bar{x}') d\omega' - \sum_{\mu} V_{\mu_0} \Gamma_{\lambda\mu}^{0Mm}(\bar{x}) : \mathbf{p}_{\mu_0} \\
\varepsilon_{\lambda}^{m}(\bar{x}) &= \varepsilon^0 - \sum_{\mu} N_{\mu} \int_{S_{\mu}} \Gamma_{\lambda\mu}^{0mM}(\bar{x}') : \mathbf{P}_{\mu}(\bar{x}') d\omega' - \sum_{\mu} V_{\mu_0} \Gamma_{\lambda\mu}^{0mm}(\bar{x}) : \mathbf{p}_{\mu_0}
\end{aligned}
\quad (3.33) $$
and

$$<\epsilon^{\prime}> = \sum_{\lambda} c_{\lambda} \left( \frac{1-\alpha_{\lambda}}{S_{\lambda}} \int_{S_{\lambda}} \epsilon_{\lambda}^{\prime M}(\bar{x}) d\omega + \alpha_{\lambda} \epsilon_{\lambda}^{\prime\prime m} \right) = E \quad (3.34)$$

with $\alpha_{\lambda} = V_{\lambda_0}/V_{\lambda}$. The stationarity condition on the Hashin-Shtrikman functional has the new form:

$$\begin{aligned} \frac{\partial H S^0}{\partial \boldsymbol{p}}: \delta \boldsymbol{p} = & \sum_{\lambda} \frac{N_{\lambda}}{V} \int_{D_{\lambda}} \delta \boldsymbol{P}_{\lambda}(\bar{x}) : [\epsilon_{\lambda}^{\prime M}(\bar{x}) - \boldsymbol{H}_{\lambda}^{0M}(\bar{x}) : \boldsymbol{P}_{\lambda}(\bar{x})] d\omega \dots \\ & + \sum_{\lambda} \frac{V_{\lambda_0}}{V} \delta \boldsymbol{p}_{\lambda_0} : (\epsilon_{\lambda}^{\prime\prime m} - \boldsymbol{H}_{\lambda}^{0m} : \boldsymbol{p}_{\lambda_0}) = 0 \end{aligned} \quad (3.35)$$

which implies

$$\epsilon_{\lambda}^{\prime M}(\bar{x}) = \mathbf{H}_{\lambda}^{0M}(\bar{x}) : \mathbf{P}_{\lambda}(\bar{x}) \quad \text{or} \quad \mathbf{P}_{\lambda}(\bar{x}) = \mathbf{H}_{\lambda}^{0M-1}(\bar{x}) : \epsilon_{\lambda}^{\prime M}(\bar{x}) \quad (3.36)$$

$$\epsilon_{\lambda}^{\prime\prime m} = \mathbf{H}_{\lambda}^{0m} : \mathbf{p}_{\lambda_0} \quad \text{or} \quad \mathbf{p}_{\lambda_0} = \mathbf{H}_{\lambda}^{0m-1} : \epsilon_{\lambda}^{\prime\prime m}$$

Integration of (3.36) into (3.33) shows these equations to still rule (composite) inclusions / matrix problems' through pattern-averaged Green operators which have to be derived from the available informations on the pattern distribution. Finally, optimal bounding is obtained by choosing $C^0$ so as to make $\mathbf{H}_{\lambda}^{0M}(\bar{x})$ and $\mathbf{H}_{\lambda}^{0m}$ (just) positive or negative. In case of an isotropic distribution, the intersections $S_{\lambda}$ must be spherical and the pattern-averaged Green operators can be determined in closed form: in addition to the result already derived for $\Gamma_{\lambda\mu}^{0MM}(\bar{x}, \bar{x}')$, namely $N_{\mu} \Gamma_{\lambda\mu}^{0MM}(\bar{x}, \bar{x}') = \delta_{\lambda\mu} \Gamma^0(\bar{x}, \bar{x}')$, $\forall \bar{x} \in S_{\lambda}, \bar{x}' \in S_{\mu}$, we find with similar arguments:

<figure><img src="image_17.png" /><figcaption>Fig. 15 Interpretation of generalized M.R.P.-based Hashin-Shtrikman bounds in terms of composite and homogeneous inclusions problems</figcaption></figure>
$$
\begin{equation}
\begin{aligned}
\Gamma_{\lambda\mu}^{0Mm}(\bar{x}) &= \Gamma_{\lambda\mu}^{0mM}(\bar{x}') = 0, \quad \forall \bar{x} \in S_{\lambda}, \bar{x}' \in S_{\mu} \\
\Gamma_{\lambda\mu}^{0mm} &= \frac{1}{V_{\mu_0}} \delta_{\lambda\mu} E^0
\end{aligned}
\tag{3.37}
\end{equation}
$$

where $E^0$ is the local part of the Green operator $\Gamma^0(\bar{x}, \bar{x}')$.

Now again, the problem reduces to several problems of composite spheres with moduli $C^0 + H_\lambda^{0M^{-1}}(\bar{x})$ and homogeneous ones with moduli $C^0 + H_\lambda^{0m^{-1}}$ in an infinite homogeneous matrix, so that Fig. 9 is now changed into Fig. 15. Note that it is still more apparent on this case that when the patterns reduce to points and the morphological phases to the mechanical ones, i.e. when $S_\lambda = S_r = \emptyset, \forall \lambda = r$, the classical Hashin-Shtrikman bounds (and their interpretation in terms of homogeneous spherical inclusions by Walpole) are recovered soon with $p_{\lambda_0} = p_{\lambda_0} = p_r$ uniform per phase. Note also that the case of composites with one phase constituting a continuous matrix and with inclusive patterns which are identical within each morphological phase has been treated more completely by Bornert [29].

### 3.2.2.2 Ellipsoidal distributions of patterns

Another way of generalization of what has been presented as a "simplified approach" in 3.2.1 is the following: still considering morphological phases consisting in identical patterns, we can deal with more general distributions of the centers than the isotropic one in order to explore more completely the sources of anisotropic behaviour of heterogeneous media. A simple way to do so refers to Willis's original idea [30] of "ellipsoidal phase distribution", here adapted to the pattern approach. To do so, we still make use of the distribution function $\Psi_{\lambda\mu}(\bar{u})$ of section 3.2.1.3 such that $\Psi_{\lambda\mu}(\bar{u})d\omega_u$ is the number of couples of domains $D_{\lambda_i}$ and $D_{\mu_j}$ such that their centers $\bar{X}_{\lambda_i}$ and $\bar{X}_{\mu_j}$ lie at the distance $\bar{u}$ within $d\omega_u$. The assumption of an ellipsoidal distribution of the pattern centers implies that $\Psi_{\lambda\mu}(\bar{u})$ can be written as

$$ \Psi_{\lambda\mu}(\bar{u}) = \psi_{\lambda\mu}(\|\mathbf{B}\cdot\bar{u}\|) \qquad (3.38) $$

where $\psi_{\lambda\mu}(\bar{x})$ is any positive real function of $\bar{x}$ and $\mathbf{B}$ any positive definite symmetric second order (appropriately normalized) matrix.

Similarly to what has been shown for an isotropic distribution, it can be proved now [6] that (3.38) implies that the outer surface of any pattern obeys the equation $\|\mathbf{B}\cdot\bar{x}\| = r$ where $r$ determines the pattern extension: this means that all the patterns must be ellipsoids with the same principal directions (those of the matrix $\mathbf{I}\mathbf{B}\cdot\mathbf{B}$, with $\mathbf{I}\mathbf{B}$ transposed from $\mathbf{B}$) and aspect ratios: we can imagine the ellipsoidal distribution as resulting from an isotropic one which has been stretched along the principal directions of $\mathbf{I}\mathbf{B}\cdot\mathbf{B}$ (see Fig. 16 as an illustration of such an operation). Nevertheless the internal geometry of the patterns is still arbitrary and may exhibit quite different anisotropic properties. An interesting possibility of this approach lies in the possibility to investigate the competition between the two corresponding sources of the overall anisotropy, associated with the pattern distribution and their material content, respectively. Examples of such a competition are given in the following.
<figure><img src="image_18.png" /><figcaption>Fig. 16 Ellipsoidal distribution of patterns: (a) external shape of the patterns; (b) possible internal content (from [6])</figcaption></figure>

The integration of such an information into Hashin-Shtrikman's variational procedure is quite straightforward and makes basically use of the same fundamental property of the Green operator as before. The pattern-averaged trial strain field (3.10) can be computed by integration on $\bar{u}$ over ellipsoidal shells with infinitesimal width where this property can be used, namely

$$ 
\begin{align} 
& N_{\mu} \int_{D_{\mu}} \Gamma^{0MM}_{\lambda\mu}(\bar{x} - \bar{x}') : \mathbf{P}_{\mu}(\bar{x}') d\omega' \dots \nonumber \\ 
& = \frac{1}{N_{\lambda}} \sum_{i=l}^{N_{\lambda}} \sum_{j=l}^{N_{\mu}} \int_{D_{\mu}} \Gamma^{0}(\bar{X}_{\lambda_{i}} - \bar{X}_{\mu_{j}} + \bar{x} - \bar{x}') : \mathbf{P}_{\mu}(\bar{x}') d\omega' \tag{3.39} \\ 
& \quad - \frac{\dots}{N_{\lambda}} \int_{D_{\mu}} d\omega_{r} \int_{D_{\mu}} \psi_{\lambda\mu}(r) \Gamma^{0}(\bar{u} + \bar{x} - \bar{x}') : \mathbf{P}_{\mu}(\bar{x}') d\omega' \nonumber 
\end{align} 
$$

By permuting the order of integration due to Fubini's theorem, the contributions vanish as long as $r \ge r_{\lambda} + r_{\mu}$, with $r_{\lambda}$ and $r_{\mu}$ denoting the extension of the ellipsoids $D_{\lambda}$ and $D_{\mu}$ (whose equations read $\|\mathbf{B} \cdot \bar{x}\| \le r_{\lambda}$ and $\|\mathbf{B} \cdot \bar{x}\| \le r_{\mu}$ respectively): this result stems from the evoked property of the Green integral for concentric similar ellipsoids. The contributions
vanish too when $r \le r_\lambda + r_\mu$ and $\lambda \ne \mu$ since $D_\lambda$ and $D_\mu$ cannot overlap; the only non zero
contribution arises from $r=0$ and $\lambda = \mu$ with $\psi_{\lambda\mu}(0) = N_\lambda \delta_{\lambda\mu}$. Consequently we get:

$$
\begin{align}
& N_{\mu} \int_{D_{\mu}} \Gamma_{\lambda\mu}^{0MM}(\bar{x}, \bar{x}') : \mathbf{P}_{\mu}(\bar{x}') d\omega' = \delta_{\lambda\mu} \int_{D_{\mu}} \Gamma^{0}(\bar{x} - \bar{x}') : \mathbf{P}_{\mu}(\bar{x}') d\omega' \tag{3.40} \\
& = \int_{D_{\lambda}} \Gamma^{0}(\bar{x} - \bar{x}') : \mathbf{P}_{\lambda}(\bar{x}') d\omega', \quad \forall \bar{x} \in D_{\lambda} \nonumber
\end{align}
$$

which is strictly the same result as (3.20), derived for an isotropic distribution, except that
$D_\lambda$ and $D_\mu$ are now similar ellipsoids instead of spheres.

The end of the demonstration and the final results and interpretation are unchanged,
except for the replacement of spheres by ellipsoids. No closed form solution can generally be
expected, even for simple internal microstructures of the patterns. But, here again, direct
F.E.M. numerical calculations can be performed in order to solve the elementary underlying
composite ellipsoidal inclusions problems which give access to the derivation of bounds for
the overall moduli. An illustration of the method is reported now for the case of transversely
isotropic composites [6].

3.2.2.3 Case of transversely isotropic materials:

We consider identically shaped and oriented (along the symmetry axis z) hard isotropic
ellipsoidal inclusions distributed in a continuous weak isotropic matrix according to an
ellipsoidal symmetry. This case allows us to compare our M.R.P.-based Hashin-Shtrikman
bounds for the transverse $\mu_{xy}^{eff}$ and longitudinal $\mu_{xz}^{eff}$ overall moduli to explicit analytical
Hashin-Shtrikman-type bounds which have been derived recently by Ponte Castañeda and
Willis [31] for particulate composites. They assumed too an ellipsoidal distribution of the
centers of the inclusions (not of the patterns !) and used piecewise constant polarization
fields in the matrix and the inclusions; as for taking the matrix connectedness into account,
they did not make any explicit assumption and just used average relations in order to make
the matrix apparently disappear from the equations. We guessed that in fact such a procedure
implied some assumption on the matrix spatial distribution, which would be unlikely to
express the matrix connectedness and we expected our M.R.P.-based Hashin-Shtrikman
upper bounds to be lower than theirs.

The Hashin-Shtrikman-type lower bound derived by Ponte Castañeda and Willis is

$$
C_{PW}^{HS} = C_I + f_2 \left[ (C_2 - C_I)^{-1} + \mathbf{P}_{inc}^{(1)} - f_2 \mathbf{P}_{dis}^{(1)} \right]^{-1} \quad (3.41)
$$

with use of the following uniform tensors:

$$
\begin{align}
\mathbf{P}_{\text{inc}}^{(i)} &= \int_{\Omega_{\text{inc}}} \Gamma^{(i)}(\bar{x} - \bar{x}') d\omega', && \forall \bar{x} \in \Omega_{\text{inc}} \nonumber \\
\mathbf{P}_{\text{dis}}^{(i)} &= \int_{\Omega_{\text{dis}}} \Gamma^{(i)}(\bar{x} - \bar{x}') d\omega', && \forall \bar{x} \in \Omega_{\text{dis}} \tag{3.42}
\end{align}
$$

where (i) refers to the used comparison medium (moduli $C_I$ for the matrix and $C_2$ for the inclusions) and $\Omega_{inc}$ and $\Omega_{dis}$ are the characteristic ellipsoids for the inclusions and for the distribution respectively. The upper bound could be derived in the same way:

$$
C_{PW}^{HS+} = C_2 - [f_1(C_2 - C_I)^{-1} - f_2 \mathbf{P}_{inc}^{(2)} + (f_2)^2 \mathbf{P}_{dis}^{(2)}]^{-1} \quad (3.43)
$$
Comparison was made for spheroidal inclusions and distribution with identical aspect ratios $w = b/a$ where $b$ denotes the length along $z$: in this case, $\Omega_{inc}$ and $\Omega_{dis}$ coincide and so do $P_{inc}^{(i)}$ and $P_{dis}^{(i)}$. Typical results are reported in Fig. 17 and Fig. 18 for $\mu_{xy}^{eff}$ and $\mu_{xz}^{eff}$ as a function of $F(w)$ (with $F(w) = w$ for oblate inclusions and $F(w) = 2 - 1/w$ for prolate ones in view of a more symmetrical presentation). In both cases the lower bounds coincide whereas the upper bounds lie far apart from each other, with definite advantage gained by the M.R.P.-based Hashin-Shtrikman bound. This result illustrates the pertinence of generalized Hashin's assemblages to the representation of the matrix connectedness.

<figure><img src="image_19.png" /><figcaption>Fig. 17 Comparison between Ponte Castañeda and Willis's and M.R.P.-based Hashin-Shtrikman's bounds for the transverse shear modulus of a composite with aligned ellipsoidal inclusions $k_1 = 3, \mu_1 = 1, k_2 = 20, \mu_2 = 10$ (arbitrary units), $f_2 = 3$, in (1) (from [6])</figcaption></figure>

We note that for $w = 1$ (isotropic distribution of spheres) Ponte Castañeda-Willis's bounds coincide with the classical Hashin-Shtrikman's ones, which indicates that, at least in this case, the matrix connectedness has not been fully exploited. As for the noticeable (numerical) coincidence of the lower bounds, it can easily be justified analytically by solving the involved problem of an ellipsoid $\Omega_{inc}$ with aspect ratio $w_{inc}$ made of material (2) embedded in the ellipsoid $\Omega_{dis}$ with aspect ratio $w_{dis}$ made of material (1) which is embedded itself in an infinite matrix with moduli $C_1$: this is nothing but the classical Eshelby problem whose well-known solution leads to (3.41) after some simple calculations [6].

We end this chapter by reporting briefly Bornert's recent developments [7, 29] concerning the continuation of the foregoing approach towards the definition of M.R.P.-based self-consistent estimates.
<figure><img src="image_20.png" /><figcaption>Fig. 18 Comparison between Ponte Castañeda and Willis's and M.R.P.-based Hashin-Shtrikman's bounds for the longitudinal shear modulus of the composite of Fig. 17 ([6])</figcaption></figure>

### 3.3 A M.R.P.-based self-consistent scheme

The principle of such a generalization of the so-called "generalized self-consistent scheme" is quite straightforward: we already noticed that the three-phase model was an extension of the classical self-consistent scheme obtained by replacing the homogeneous ellipsoids representative of the mechanical phases by the composite spheres (or cylinders) of the C.S.A. (or of the C.C.A.). We also justified the possibility of dealing simultaneously with different kinds of composite spheres (or cylinders) within an extension of the three-phase model. The next step consists in dealing with several families of morphological representative patterns (instead of composite spheres or cylinders) so as to derive an estimate of the overall moduli from the numerical resolution of as many composite inclusion / matrix problems as representative patterns: the matrix moduli have to be such that they relate the average stress and strain tensors over the different representative patterns embedded in this matrix (Fig. 19).

Another way of understanding such a procedure is the following [7]: when looking for M.R.P.-based Hashin-Shtrikman bounds, we had first to make an optimal choice of the polarization field $p$, say $p_{opt}^0$, for given values of the moduli $C^0$ of the reference medium, which led to the relation (see 3.7, 3.13 and 3.16):

$$ \begin{gathered} \Phi_{opt}^0(p_{opt}^0, E) = \frac{1}{2}E: \langle p_{opt}^0 \rangle \qquad (3.44) \\ \text{with } \Phi^0(p, E) = HS^0(p, E) - W^0(E) = HS^0(p, E) - \frac{1}{2}E:C^0:E \end{gathered} $$
We had then to chose the moduli $C^0$ of the reference medium so as to be sure, at minimal cost, of the sign of $\Phi_{opt}^0(p_{opt}, E)$ in view of bounding the effective moduli. The self-consistent point of view aims at minimizing $|\Phi^0(p, E)|$ with $C^0 = C^{eff}$ so as to have an estimate of $W^{eff}(E)$ through $HS_{opt}^{eff}(p_{opt}^{eff}, E)$, with $\langle p_{opt}^{eff} \rangle = \mathbf{0}$. Such a relation is equivalent to $\langle \delta c^{eff} : \epsilon_{opt}' \rangle = \mathbf{0}$ or $\langle \sigma_{opt}' \rangle = C^{eff} : \langle \epsilon_{opt}' \rangle$.

<figure><img src="image_21.png" /><figcaption>Fig. 19 Principle of the M.R.P.-based self-consistent scheme</figcaption></figure>

From a practical point of view, several difficulties have still to be overcome. We refer to Bornert's work [7, 29] for more details. Suffice it to say here that consistency of the foregoing definition with an energetical one forces the external shape of the patterns to be ellipsoidal. Thus the numerical procedure developed for the derivation of M.R.P.-based Voigt-Reuss or Hashin-Shtrikman bounds in case of en ellipsoidal distribution of patterns can still be used within an iterative procedure which can start from any initial reference medium $C^{(0)}$ to compute a first estimate $C^{(2)}$, then to compute $C^{(4)}$ from $C^{(2)}$ and so on until some convergence criterion $\|C^{(2n)} - C^{(2n-2)}\| \le \epsilon$ be satisfied. Different pattern geometries have been considered up to now [7] (including non similar concentric ellipsoids) in order to investigate different sources of anisotropy in particulate composites.

It can be appreciated that such an approach yields a useful alternative to the "unit cell" or periodic homogenization techniques. At the price of a numerical treatment somehow heavier due the iterative procedure, the assumption of a periodical distribution of patterns which may happen to be embarrassing is given up; in addition several families of patterns may be considered simultaneously, which may offer definite advantages in many cases, especially when local strain or stress investigations are developed, such as in view of analyzing plastic localization or damage initiation and growth processes. Nevertheless, this approach would need to be extended to nonlinear behaviour, which will be seen not to be quite straightforward in what follows...
# 4. BEHAVIORAL ASPECTS

## 4.1 Introduction

Up to now linear elasticity only has been considered. A more complete understanding of the influence of structural morphology on the overall mechanical properties of inhomogeneous materials obviously needs to explore nonlinear constitutive behaviours too. The basic difficulty one has to overcome is twofold:

- on the one hand, most of the micromechanical methods reported hereabove in case of linear elasticity make use intensively of the Green techniques which are basically attached to linearity. The question is: to what extent some linearization procedures could allow us to preserve anything of these techniques when dealing with nonlinearity?

- on the other hand, nonlinearity implies inhomogeneity of the mechanical properties, even in a homogeneous material, as soon as it is stressed or strained in an inhomogeneous manner. This means that we have to deal not only, as we did up to now, with interphase inhomogeneity but also with "intraphase" inhomogeneity of the material parameters. In other words, the morphological description and integration of the phase spatial distribution must be somehow complemented by those of the material parameters within each phase, the latter being basically changing along any given load path (think of the inextricable extension of the correlation functions of the elastic moduli to the case of nonlinearity !...).

That is why the question of nonlinear continuum micromechanics is a matter of intensive current research which will not be tried to be reported here. In order to give an outline of the difficulties and possible issues but still concentrate on the influence of morphology on the overall behaviour, restriction will be made to the extension of two simple but basic models in this field, namely the C.S.C.S. and the G.S.C.S.: the first one is well-suited to a "polycrystal-type morphology" and the second one to a "composite-type morphology". Comparison of their predicted responses to the same loading path can allow to draw significant conclusions on the influence of the phase connectedness.

In view of going from the simpler to the harder, we start with the easiest way to depart from elasticity but still preserving the convenience of linearity, namely with the case of linear nonageing viscoelasticity.

## 4.2 The case of linear nonageing viscoelasticity

In this case, obviously the mechanical characteristics will not be space-dependent inside the same phase, due to the linear behaviour, but the rate-dependence is responsible for many new phenomena and properties which are characteristic of the coupling between elasticity and viscosity. One of the most important manifestations of this coupling in heterogeneous media is the so-called "long-range memory effect" which results from the complex delayed mechanical interactions between the constituent phases. This effect finds expression, for instance, in the fact that the overall behaviour of a mixture of Maxwellian constituents does not obey a Maxwell constitutive behaviour anymore.

Let the local constitutive equations read

$$ \dot{\epsilon} = a: \sigma + b: \dot{\sigma} \qquad (4.1) $$

with $a$ and $b$ fourth order constant tensors. Then the overall behaviour will obey a different law which can be written in the form [32]

$$ \dot{E} = A^{\text{eff}}: \Sigma + B^{\text{eff}}: \dot{\Sigma} + \int_{0}^{t} J(t-s): \dot{\Sigma}(s)ds \qquad (4.2) $$
where the tensorial kernel $J$ is directly related to this "long range memory effect". Its mathematical expression is, of course, model-dependent and then morphology-dependent. So, it would be interesting to focus attention on the comparison between its properties as predicted, other things being equal, by the C.S.C.S. and the G.S.C.S. respectively in order to appreciate the influence of the microstructural morphology on this memory effect.

A convenient way to do so is to derive the relaxation spectrum of a two-phase isotropic material according to both these models. Whereas this spectrum reduces to a single line for each constituent phase, it is expected to be more complex for the two-phase material in a way which could differ according to the morphology. For the sake of simplicity, we assume incompressibility and isotropy for each phase (i) so that (4.1) reduces to

$$ \dot{e} = a_i s + b_i \dot{s} \quad i = 1, 2 \qquad (4.3) $$

where $e$ and $s$ are the strain and stress deviators and $a_i$ and $b_i$ are material constants which define the single relaxation time $T_i = b_i/a_i$ of each phase. The question is: which are the relaxation times of the two-phase material according to the C.S.C.S. and the G.S.C.S.?

The answer [8] can be derived easily through the use of the Laplace transform technique, since it is well-known, due to the "correspondence principle" [33], that this transformation allows us to convert linear (nonageing) viscoelasticity into (symbolical) elasticity: suffice it, in order to get the equations of the viscoelastic problem, to replace the elastic moduli and variables of the elastic self-consistent equations by the Laplace-transform of their viscoelastic counterpart [34]. Transformation of (4.3) yields:

$$ s^L(p) = 2p\mu_i^L(p)e^L(p), \quad \mu_i^L(p) = \frac{1}{2b_i(p + 1/T_i)} \qquad (4.4) $$

where $p$ is the complex variable, $f^L(p)$ the Laplace transform of $f(t)$ and $\mu_i(t)$ the shear relaxation function of phase (i).

In the elastic case, both considered models yield the overall shear modulus $\mu^{eff}$ as the positive root of a quadratic equation which may be written in the following common form:

$$ L(\beta, f)X^2 + 2M(\beta, f)X + N(\beta, f) = 0 \qquad (4.5) $$

where $X = \mu^{eff}/\mu_1$, $\beta = \mu_2/\mu_1$ and $f = f_2$. For the C.S.C.S., we have, from (2.8):

$$ L = 1, \quad M = \frac{2 - 5f}{6}\beta + \frac{5f - 3}{6}, \quad N = -\frac{2\beta}{3} \qquad (4.6) $$

and for the G.S.C.S., with phase (2) included in phase (1), from (2.28):

$$ \begin{aligned} L(\beta, f) &= 4[3(\beta - 1)x^3 - \eta_3](\eta_1 x^7 - 2\eta_2) - 126\eta_2(\beta - 1)x^3(1 - x^2)^2 \\ M(\beta, f) &= 3\eta_1(\beta - 1)x^{10} + 4\eta_1\eta_3 x^7 - \frac{129}{4}\eta_2(\beta - 1)x^3 + \dots \\ & \quad + \frac{3}{4}\eta_2\eta_3 + 126\eta_2(\beta - 1)x^3(1 - x^2)^2 \\ N(\beta, f) &= -[\frac{9}{2}(\beta - 1)x^3 + \eta_3](\eta_1 x^7 + \frac{19}{2}\eta_2) - 126\eta_2(\beta - 1)x^3(1 - x^2)^2, \\ x &= f^{1/3}, \quad \eta_1 = 19(\beta - 1), \quad \eta_2 = \frac{19\beta + 16}{2}, \quad \eta_3 = \frac{3}{2}(2\beta + 3). \end{aligned} \qquad (4.7) $$
The solution of the viscoelastic problem is obtained by replacing in (4.5), (4.6) and (4.7) the quantities X by $X^L(p) = \mu^{eff L}(p)/\mu_1^L(p)$ and $\beta$ by $\beta^L(p)$ defined by:

$$ \beta^L(p) = \frac{\mu_2^L(p)}{\mu_1^L(p)} = k \frac{p + I/T_1}{p + I/T_2} \quad \text{with} \quad k = b_1/b_2. \qquad (4.8) $$

The solution $\mu^{eff L}(p)$ can be written as the sum of two terms, $f_1^L(p)$ and $f_2^L(p)$, namely:

$$ \begin{gathered} f_1^L(p) = -\frac{M(\beta^L(p), c)}{2b_1(p + 1/T_1)L(\beta^L(p), c)} \\ f_2^L(p) = \frac{[M^2(\beta^L(p), c) - L(\beta^L(p), c)N(\beta^L(p), c)]^{1/2}}{2b_1(p + 1/T_1)L(\beta^L(p), c)} \end{gathered} \qquad (4.9) $$

We have now to go back into the real space:

- the Laplace inversion of $f_1^L(p)$ is easy and yields directly the sum of two (for the C.S.C.S.) or three (for the G.S.C.S.) exponential functions : in terms of spectral representation, this corresponds to two (for the C.S.C.S.) or three (for the G.S.C.S.) discrete lines, at times $T_1$ and $T_2$ for the former and $T_1$, $\theta_1$ and $\theta_2$ for the latter, with $\theta_1$ and $\theta_2$ lying between $T_1$ and $T_2$; note that their intensity can be negative for some values of the volume fraction $f$.

- the inversion of the second term needs recourse to the inverse Laplace integral formula

$$ f_2(t) = \frac{1}{2i\pi} \int_{\Delta} f_2^L(p) e^{pt} dp \qquad (4.10) $$

where $\Delta$ is a vertical axis which has to leave on its left all the critical points of $f_2^L(p)$. This is performed after some rewriting of (4.9b). For the C.S.C.S, we write $f_2^{CL}(p)$ as

$$ \begin{align} f_2^{CL}(p) &= k_C \frac{P_C(p)^{1/2}}{(p + 1/T_1)(p + 1/T_2)} \nonumber \\ P_C(p) &= (p + 1/\tau_1)(p + 1/\tau_2) \tag{4.11} \\ k_C &= \frac{\sqrt{(2-5f)^2 k^2 + 2(6-5f)(5f+1)k + (5f-3)^2}}{12b_1} \nonumber \\ k &= b_1/b_2 \nonumber \end{align} $$

where times $\tau_1$ and $\tau_2$ lie between $T_1$ and $T_2$; for the G.S.C.S., we write $f_2^{GL}(p)$ as

$$ \begin{align} f_2^{GL}(p) &= k_G \frac{P_G(p)^{1/2}}{(p + 1/T_1)(p + 1/\theta_1)(p + 1/\theta_2)}, \quad k_G \ge 0 \tag{4.12} \\ P_G(p) &= (p + 1/\tau'_1)(p + 1/\tau'_2)(p + 1/\tau'_3)(p + 1/\tau'_4) \nonumber \end{align} $$

where $\tau'_1$ and $\tau'_2$ lie between $T_1$ and $\theta_1$ and $\tau'_3$ and $\tau'_4$ between $\theta_1$ and $\theta_2$. Application of (4.10), of the theorem of residues and Jordan's lemma and definition of adequate cuts on
the real negative axis (see Fig. 20) lead to the following expressions (with $T_1<T_2$ to help discussion):

$$f_2^C(t) = \pm \frac{k_C}{\pi} \left\{ -\int_{-1/\tau_1}^{-1/\tau_2} \frac{\sqrt{-P_C(x)}}{(x+1/T_1)(x+1/T_2)} e^{\lambda x} dx - \dots \right. \\ \left. \qquad \dots - \pi \frac{\sqrt{P_C(-1/T_1)}}{1/T_2 - 1/T_1} e^{-\lambda/T_1} + \pi \frac{\sqrt{P_C(-1/T_2)}}{1/T_1 - 1/T_2} e^{-\lambda/T_2} \right\}$$

$$f_2^G(t) = \pm \frac{k_G}{\pi} \left\{ \int_{-1/\tau'_1}^{-1/\tau'_2} \frac{\sqrt{-P_G(x)}}{(x+1/T_1)(x+1/\theta_1)(x+1/\theta_2)} e^{\lambda x} dx - \dots \right. \\ \left. \qquad \dots - \int_{-1/\tau'_3}^{-1/\tau'_4} \frac{\sqrt{-P_G(x)}}{(x+1/T_1)(x+1/\theta_1)(x+1/\theta_2)} e^{\lambda x} dx - \dots \right. \\ \left. \qquad \dots - \frac{\pi}{(1/\theta_2 - 1/\theta_1)} \left[ \frac{\sqrt{P_G(-1/\theta_1)}}{(1/T_1 - 1/\theta_1)} e^{-\lambda/\theta_1} - \frac{\sqrt{P_G(-1/\theta_2)}}{(1/T_1 - 1/\theta_2)} e^{-\lambda/\theta_2} \right] \right\} \quad (4.13)$$

<figure><img src="image_22.png" /><figcaption>Fig. 20 Integration paths used for the inversion of the Laplace transform (from [8]):
(a) C.S.C.S., (b) G.S.C.S. (with $T_1 < T_2$)</figcaption></figure>

The involved integrals can be written, after putting $\tau = -1/x$, in the form of a spectral representation, namely
$$ \mu(t) = \int_{0}^{\infty} g(\tau) e^{-t/\tau} d\tau \qquad (4.14) $$

which corresponds to the continuous spectrum $g(\tau)$, whereas the exponential functions outside the integrals are associated to discrete lines.

Addition of $f_I^C(t)$ for the C.S.C.S. or of $f_I^G(t)$ for the G.S.C.S. lead to the following final results:

- the shear relaxation spectrum as predicted by the C.S.C.S. consists of a continuous part and additional discrete lines. The continuous spectrum extends from $\tau_1$ to $\tau_2$ (with $T_1 < \tau_1 < \tau_2 < T_2$) and its intensity is given by

$$ g_C(\tau) = \frac{k_C T_1 T_2 \sqrt{(\tau - \tau_1)(\tau_2 - \tau)}}{\pi \sqrt{\tau_1 \tau_2} (\tau(\tau - T_1)(T_2 - \tau))}, \quad \tau \in [\tau_1, \tau_2]. \qquad (4.15) $$

The discrete lines lie at $\tau = T_1$ for $f \le .6$ (with the intensity $(3-5f)/6b_1$) and at $\tau = T_2$ for $f \ge .4$ (with the intensity $(5f-2)/6b_2$).

- the shear relaxation spectrum predicted by the G.S.C.S. consists of two separate continuous parts and additional discrete lines. The continuous spectrum lies between $\tau'_1$ and $\tau'_2$ and $\tau'_3$ and $\tau'_4$ respectively, with the intensity

$$ g_G(\tau) = \frac{\epsilon(\tau) k_G T_1 \theta_1 \theta_2}{\pi \sqrt{\tau'_1 \tau'_2 \tau'_3 \tau'_4}} \frac{\sqrt{(\tau - \tau'_1)(\tau - \tau'_2)(\tau - \tau'_3)(\tau'_4 - \tau)}}{\tau(\tau - T_1)(\theta_1 - \tau)(\theta_2 - \tau)} \\ \tau \in [\tau'_1, \tau'_2] \cup [\tau'_3, \tau'_4] \qquad (4.16) $$

with $\epsilon(\tau) = 1$ if $\tau \in [\tau'_1, \tau'_2]$ and $\epsilon(\tau) = -1$ if $\tau \in [\tau'_3, \tau'_4]$.

The discrete lines lie at $\tau = T_1$ in a fixed range of volume fraction $f$ and at $\tau = \theta_1$ or $\theta_2$. or both according to the value of $f$, $T_1$, $T_2$, $k$ and $b_1$.

We can see on Fig. 21 and 22 illustrative examples of these results as predicted, for the same values of the material parameters $a_i$, $b_i$ and $f$, by the C.S.C.S. and the G.S.C.S. respectively. Though such results seem hard to be corroborated by direct experimental investigation, they yield useful information on several micromechanical aspects:

- both models predict a non-Maxwellian overall behaviour from Maxwellian constituents. Whether the resultant spectrum is continuous or not, it definitely does not reduce to a single line: the "long-range memory effect" may be strong enough to result, even for a simplistic two-phase material, in a continuous (bounded) spectrum, which means a continuous infinity of relaxation times which are to be related to the mechanical interactions between the phases and their viscoelatic nature. This result must be kept in mind in view of the future treatment of nonlinear viscoelasticity.

- obviously morphology does affect the spectrum shape: one the one hand the "symmetrical" morphology associated with the C.S.C.S. leads to a unique continuous spectrum which reflects the underlying intricate phase distribution; on the other hand, the unsymmetrical morphology inherent to the G.S.C.S. is responsible for a clear spectrum splitting into two parts and for a spectrum shape which reflects the prominent mechanical role played by the connected phase. This qualitative difference between the spectra might be a sensible indicator of morphological changes which can occur during phase transitions, such as the glass transition of polymers.
The foregoing analysis could be extended to more general constitutive equations than Maxwell's ones without too much efforts [35] or to more complex morphologies: especially the existence of an interphase between the matrix and the particles, which can be taken into account with a four-phase model, can be shown to be responsible for a specific spectrum part which could be used as a "morphological signature" [36].

<figure><img src="image_23.png" /><figcaption>Fig. 21 Shear relaxation spectrum as predicted by the C.S.C.S. ($a_1 = a_2 = 1$; $T_1 = 1$; $T_2 = 10$; $f = .5$) (from [8]).</figcaption></figure>

<figure><img src="image_24.png" /><figcaption>Fig. 22 Shear relaxation spectrum as predicted by the G.S.C.S. ($a_1 = a_2 = 1$; $T_1 = 1$; $T_2 = 10$; $f = .5$) (from [8]).</figcaption></figure>
## 4.3 Nonlinear classical self-consistent modelling

### 4.3.1 Introduction

Let us first stress the fact that the problem of nonlinearity in the context of classical self-consistent modelling is still an open question. Though it was first addressed by Kröner [37] thirty-five years ago for the case of elastic-plastic polycrystals, it can be considered as unsolved yet despite the intensive research which has been, and is still being, devoted to this topic.

Two main approaches have been developed:

- the first one has been initiated by Kröner; it is based on the concept of "stress-free strain" or "eigenstrain" applied to the plastic or viscoplastic part of the total strain and makes extensive use of Eshelby's solution of the inclusion problem [1]. Strictly speaking, this problem refers to an ellipsoidal region of an infinite unloaded elastic medium, undergoing a uniform eigenstrain, i.e. a prestrain which could be maintained without stress in an isolated volume element. The basic idea consists in deriving the self-consistent equation from the solution of the straightforward extension of Eshelby's problem corresponding to the case when the matrix itself undergoes the uniform macroscopic plastic or viscoplastic strain and is subjected to the macroscopic strain or stress at infinity: the inclusion is representative of one particular phase (usually a set of identically oriented grains of a polycrystal) while the matrix is constituted of the unknown H.E.M. which is searched for (usually the homogenized polycrystal). This approach can be extended easily to inhomogeneous elasticity by referring to Eshelby's solution of the problem of the "inhomogeneous inclusion".

- the second one refers to a new formulation proposed by Hill [9] for rate-independent plasticity, starting from the criticism of some aspects of the foregoing approach: it relies upon a linearization procedure along the prescribed macroscopic loading path and on the flow theory of plasticity relating the strain and stress rates through the instantaneous elastoplastic moduli. Now the self-consistent equation can be derived from the solution of Eshelby's inhomogeneity problem, which is used within a rate (or incremental) formulation. For proportional loadings, the deformation theory of plasticity can be used as well by dealing with secant (instead of tangent) elastoplastic moduli. This approach has been extended to finite strain formulations as well as to rate-dependent plasticity.

Though Hill's approach considerably improved Kröner's one, it will be shown nevertheless to be not fully self-consistent anymore, so that the problem has still to be considered as an open one. Before drawing this conclusion, let us go into more details about the current state-of-the-art in this field.

### 4.3.2 Linear classical self-consistent modelling with eigenstrains and Kröner's model

In order to appreciate better the limits of Kröner's initial approach of elastoplastic self-consistent modelling, it is useful to start from apparently a different problem, namely the linear elastic classical self-consistent modelling when fixed eigenstrains $\epsilon^T$ are present in the constituent phases. This is a classical question whose solution is well-known as soon as it has been solved by any model when no eigenstrains exist. Let $\overline{A_r^{eff}}$ and $\overline{B_r^{eff}}$ be the average stress concentration tensors for the elastic phase ($r$) in this case, so that, according to (2.1) and (2.2), the overall moduli and compliances are given by

$$ 
\begin{align} 
C^{eff} &= \sum_r f_r c_r : \overline{A_r^{eff}} = <c: A^{eff}> \tag{4.17} \\ 
S^{eff} &= \sum_r f_r s_r : \overline{B_r^{eff}} = <s: B^{eff}> \text{ with } \overline{B_r^{eff}} = c_r : \overline{A_r^{eff}} : <c: A^{eff}>^{-1}. \nonumber 
\end{align} 
$$
When eigenstrains are present, the local and global constitutive equations read

$$ \sigma = c:(\varepsilon - \varepsilon^T) \\ \Sigma = C^{\text{eff}}:(E - E^T) \qquad (4.18) $$

with $E^T = \langle t \overline{B}^{\text{eff}} : \varepsilon^T \rangle$

or, when $\varepsilon^T$ is piecewise uniform per phase, say $\varepsilon_r^T$:

$$ E^T = \sum_r f_r \overline{B_r^{\text{eff}}} : \varepsilon_r^T \qquad (4.19) $$

If we use the self-consistent scheme, $\overline{A_r^{SCSC}}$ and $\overline{B_r^{SCSC}}$ are given by (2.5) and (2.6) with $C^0 = C^{SCSC}$ and the strain concentration relation reads

$$ \varepsilon_r = (I + \overline{P_r^{CSCS}} : \delta \overline{c_r^{CSCS}})^{-1} : \{\varepsilon_0 + \overline{P_r^{CSCS}} : (\overline{c_r} : \varepsilon_r^T - C^{CSCS} : E^T)\} \qquad (4.20) $$

with $\varepsilon_0$ defined from: $\langle \varepsilon \rangle = E$.

Consider now the case of uniform isotropic elasticity, isochoric eigenstrains denoted $\varepsilon_r^p$ and spherical inclusions. We can see easily that $\delta \overline{c_r^{CSCS}} = 0$ and $\varepsilon_0 = E$, whereas $\overline{P_r^{CSCS}}$ is given by (2.19) where the deviatoric part only has to be used, so that (4.20) becomes

$$ \varepsilon_r = E + \beta(\varepsilon_r^p - E^p) \qquad (4.21) $$

with $E^p = \langle \varepsilon^p \rangle$. This may be written in the equivalent form

$$ \sigma = \Sigma + 2\mu(1 - \beta)(E^p - \varepsilon^p) \qquad (4.22) $$

This relation is nothing but Kröner's interaction equation [37] which was the first one to have been proposed to rule the self-consistent approach to polycrystal plasticity.

This means that Kröner's approach to elastoplastic self-consistent modelling actually reduces to an elastic one with eigenstrains and does not really takes into account the (stress-dependent) plastic flow of the polycrystal. Instead of expressing the actual elastoplastic interactions between the phases, this model is considering far too strong elastic ones: this stems clearly from the underlying inclusion problem which deals with a uniformly plastified matrix whose plastic strain $E^p$ is not disturbed by the inclusion. Referring to the set of estimates reported in section 2, we could say that instead of an elastic (i.e. quasi-rigid) matrix, which leads to a Voigt-type overall behaviour, we should have considered a plastically flowing matrix, corresponding to the actual polycrystal in the plastic regime. So, it is not surprising that the stress-strain curves and textures predicted by this model are almost identical to those derived from the Taylor model which assumes uniform strains. This can be understood from the fact that, in (4.22), the stress deviations $(\sigma - \Sigma)$ are of the order of magnitude of the yield stress whereas the term $2\mu(1-\beta)$ is almost equal to $\mu$, i.e. is $10^2$ to $10^3$ larger, so that the strain fluctuations from grain to grain cannot exceed $10^{-3}$ to $10^{-2}$, which coincides practically with Taylor's assumption.

The main point open to criticism in the foregoing approach lies in the fact that the relations (4.18) have been considered (4.18) as constitutive equations for a plastically flowing inhomogeneous body. As a matter of fact, these equations rule the elastic part of the
constitutive behaviour only and do not express at all the stress-dependence of the (local as well as global) plastic flow in the plastic regime. From this point of view, it could be dangerous to consider the plastic strain as an eigenstrain. On the other hand, stress-independent eigenstrains can be dealt with according to the foregoing approach even if they are time-dependent, such as, for instance, thermal strains under variable temperature (we will make use of this comment later).

Note that for such actual eigenstrains which do not depend on the purely mechanical loading, Kröner's approach can be applied to linear viscoelasticity as well, with use of the correspondence principle and the Laplace-transform technique. Suffice it for that to replace in eqns (4.17) to (4.22) the real elastic variables by their Laplace-transformed viscoelastic counterpart. Use of such a possibility will be made later when dealing with nonlinear viscoelasticity.

### 4.3.3 Elastoplasticity and viscoplasticity

An alternative treatment of elastoplastic self-consistent modelling was proposed by Hill, making use of the flow theory of plasticity. Instead of (4.17), the constitutive equations were written in the form

$$ 
\begin{aligned}
\dot{\sigma} &= l_{ep} \cdot \dot{\epsilon} \\
\dot{\Sigma} &= L_{ep}^{CSCS} \cdot \dot{E}
\end{aligned}
\qquad (4.23)
$$

Such a linearization procedure, using the local elastoplastic instantaneous moduli $l_{ep}$ and the overall ones $L_{ep}^{CSCS}$ and still keeping in mind that these are multi-branched quantities, allows us to convert the elastoplastic problem into a (pseudo)elastic one: at any step of a given loading path, the elastoplastic phases of the inhomogeneous material may be considered as elastic ones for the infinitesimal coming step with moduli $l_{ep}$. Instead of the inclusion problem with stepwise uniform plastic strains, the corresponding elementary problem is Eshelby's problem of an ellipsoidal inhomogeneity with uniform moduli $l_{ep}$ embedded in an infinite matrix with uniform moduli $L_{ep}^{CSCS}$ submitted to the prescribed homogeneous stress $\dot{\Sigma}$ or strain $\dot{E}$ rate at infinity (Fig. 23).

<figure><img src="image_25.png" /><figcaption>Fig. 23 Kröner's (a) and Hill's (b) reference Eshelby's problems</figcaption></figure>
The resulting self-consistent equations are nothing but the elastic ones (2.7) where c and C<sup>CSCS</sup> are to be replaced by l<sub>ep</sub> and L<sub>ep</sub><sup>CSCS</sup> respectively. As for the interaction (or concentration) rule, it reads from (2.5) or (2.6):

$$
\bar{\epsilon}_r = (I + P_r^{CSCS} : \delta l_r^{CSCS})^{-1} : \langle (I + P^{CSCS} : \delta l^{CSCS})^{-1} \rangle^{-1} : \dot{E} \quad (4.24)
$$

where $\mathbf{P}_r^{CSCS}$ has been defined from the Green operator associated with $\mathbf{L}_{ep}^{CSCS}$ and with
$\delta l_r^{CSCS} = l_{ep,r} - \mathbf{L}_{ep}^{CSCS}$. In case of ellipsoids with the same shape and orientations, this
expression reduces to

$$
\bar{\epsilon}_r = (I + \mathbf{P}_r^{\mathrm{CSCS}} : \delta \mathbf{l}_r^{\mathrm{CSCS}})^{-1} : \dot{\mathbf{E}} \quad (4.25)
$$

which has been put by Hill in the equivalent form:

$$
\dot{\sigma} = \dot{\Sigma} + L^* : (\dot{E} - \dot{\epsilon}) \quad (4.26)
$$

where the "constraint tensor" L*, depending on L<sup>CSCS</sup><sub>ep</sub> and on the shape and orientation of
the ellipsoids, is given by

$$
\begin{align}
\boldsymbol{L}^* &= \boldsymbol{L}_{ep}^{\mathrm{CSCS}} : (\boldsymbol{S}_{Esh}^{\mathrm{CSCS}})^{-1} - \boldsymbol{I} \tag{4.27} \\
\boldsymbol{S}_{Esh}^{\mathrm{CSCS}} &= \boldsymbol{P}^{\mathrm{CSCS}} : \boldsymbol{L}_{ep}^{\mathrm{CSCS}} \nonumber
\end{align}
$$

with help of Eshelby's tensor $S_{Esh}^{CSCS}$.

Note that, for proportional monotonic loading, the deformation theory could be used as
well [38]: for isotropic plasticity, it yields an interaction law which only differs from
Kröner's one (4.22) by the fact that the term 2μ(1-β) is multiplied by a scalar "elastoplastic
secant accommodation factor" α<sub>ep</sub> which can be shown to be equal to 1 only in the elastic
regime and to rapidly decrease towards 10<sup>-1</sup> to 10<sup>-2</sup> as soon as the overall plastic flow
occurs. This allows the plastic strain deviations to be larger and this model to yield
predictions which can strongly differ from Taylor's ones, frequently in far better agreement
with experimental results.

A quite similar treatment can be used for the case of viscoplasticity [39] by simply
replacing in (4.23) l<sub>ep</sub> and L<sub>ep</sub><sup>CSCS</sup> by their viscoplastic counterpart h<sub>vp</sub> and H<sub>vp</sub><sup>CSCS</sup> and the
first strain time derivatives ε̇ and Ė by the second ones ε̇̇ and Ė̇. Both approaches have
been extended to finite strain formulations and have successfully predicted the formation of
crystallographic textures in a number of metal forming situations.

4.3.4 Rate-dependent plasticity and nonlinear viscoelasticity

Any extension of Hill's formulation to rate-dependent elastoplasticity has to deal with
the difficult coupling between elasticity and viscosity: such a coupling makes stress and
strain time-derivatives of different orders appear simultaneously in the constitutive equations.
So, Hill's linearization procedure cannot be used anymore. Some authors [40] argued that
Hill's criticism of Kröner's approach should not apply anymore to rate-dependent plasticity
since the viscoplastic strain rate depends at any time on the current stress and not of the
stress rate so that, dislike the plastic strain, the viscoplastic one ε<sup>vp</sup> could be treated as an
eigenstrain: as we have seen before, this question refers to the mechanical dependence or
independence of the considered strain whatever the mechanical variable may be concerned.
Other authors [41, 42] did not explicitly discussed this point but they used eqns (4.18) as basic constitutive equations (with $\varepsilon^T = \varepsilon^{vp}$) in a Green analysis at finite strain: it is easy to see [43] that, in the special case of linear viscoelasticity, this would correspond to the presumption that a mixture of Maxwellian constituents obeys a Maxwellian overall behaviour too, which we know to be wrong. By doing so, one is led directly to Kröner's model (eventually extended to finite strain). For reasons similar to the previous ones, one would expect from such models predictions hardly different from Taylor's ones: Fig 24 gives an illustration of this quasi-coincidence...

<figure><img src="image_26.png" /><figcaption>Fig. 24 Comparison of Harren's (S.C.) and Taylor's (T.M.) model predictions (from [41])</figcaption></figure>

On the other hand, we cannot deal with the elastic-viscous coupling as we did in the linear case since nonlinearity prevents to use the convenient Laplace transform technique. Nevertheless, by use of Hill's linearization procedure, this technique can still be utilized. Let us first consider [10] a simple local nonlinear viscoelastic behaviour obeying the equation

$$ \dot{\varepsilon} = s: \dot{\sigma} + g(\sigma) \qquad (4.28) $$

where s are the elastic compliances. Starting from time $t = 0$, we are supposed to have determined already the local and overall responses to some prescribed macroscopic loading path up to $t = t_n$; now we look for these responses during the subsequent infinitesimal time interval $[t_n, t_n + dt]$. So we can approximate (4.28) by:

$$ \begin{gathered} \dot{\varepsilon}(t) = s: \dot{\sigma}(t) + m_n: \sigma(t) + \dot{\varepsilon}_n^0(t, t_n), \quad m_n = \frac{\partial g}{\partial \sigma}[\sigma(t_n)] \\ \dot{\varepsilon}_n^0(t, t_n) = g[\sigma(t_n)] - m_n: \sigma(t_n) + \dots \\ \dots + [g[\sigma(t)] - g[\sigma(t_n)] - m_n: [\sigma(t) - \sigma(t_n)]] [1 - H(t - t_n)], \end{gathered} \qquad (4.29) $$
where $H(t - \tau)$ is the unit step function at time $\tau$. These equations show unambiguously that the considered linearized behaviour is a Maxwellian one with the eigenstrain rate $\dot{\epsilon}_n^0(t, t_n)$. Such a strain is actually a (variable) eigenstrain because its variation is completely known a priori and does not depend on the external loading which is applied beyond $t_n$; its time derivative is constant beyond $t_n$ and variable but known on $[0, t_n]$. Thus we may use the method suggested hereabove for this case by using the Laplace or Laplace-Carson transform technique. The latter is defined by:

$$f^*(p) = p \int_0^\infty f(t)e^{-pt}dt \qquad (4.30)$$

It can apply to $\epsilon_n^0(t, t_n)$ as well so as to convert the problem into a (symbolical) elastic one with eigenstrains. As already proved, the whole set of equations is the following:

$$\epsilon^* = s^* : \sigma^* + \epsilon_n^{0*} \quad \text{or} \quad \sigma^* = c^* : (\epsilon^* - \epsilon_n^{0*}),$$

$$s^*(p) = s + \frac{m_n}{p}, \quad c^*(p) = s^{*-1}(p),$$

$$E^* = S^{\text{eff}*} : \Sigma^* + E_n^{0*} \quad \text{or} \quad \Sigma^* = C^{\text{eff}*} : (E^* - E_n^{0*}), \qquad (4.31)$$

$$C^{\text{eff}*}(p) = \langle c^*(p) : A^*(p) \rangle \quad \text{or} \quad S^{\text{eff}*}(p) = \langle s^*(p) : B^*(p) \rangle, \\ B^*(p) = c^*(p) : A^*(p) : \langle c^*(p) : A^*(p) \rangle^{-1}, \\ E_n^{0*}(p) = \langle t^* B^*(p) : \epsilon_n^{0*}(p) \rangle.$$

Due to the nonlinearity, we need an additional interaction relation between the local and global variables, dependent on the chosen model, in order to determine at each step the mechanical state of each phase. For the C.S.C.S., this can be made easily from what precedes by applying the Laplace-Carson transformation to (2.5) and (4.20):

$$A_r^{\text{CSCS}^*}(p) = [\mathbf{I} + \mathbf{P}_r^{\text{CSCS}^*}(p) : \delta\mathbf{c}_r^{\text{CSCS}^*}(p)]^{-1} : \dots \\ \dots : \langle [\mathbf{I} + \mathbf{P}^{\text{CSCS}^*}(p) : \delta\mathbf{c}^{\text{CSCS}^*}(p)]^{-1} \rangle^{-1}, \\ \epsilon_r^*(p) = [\mathbf{I} + \mathbf{P}_r^{\text{CSCS}^*}(p) : \delta\mathbf{c}_r^{\text{CSCS}^*}(p)]^{-1} : \{\epsilon_0^*(p) + \dots + \mathbf{P}_r^{\text{CSCS}^*}(p) : [\mathbf{c}_r^*(p) : \epsilon_{nr}^{0*}(p) - \mathbf{C}^{\text{CSCS}^*}(p) : \mathbf{E}_n^{0*}(p)]\}, \\ \langle \epsilon^*(p) \rangle = \mathbf{E}^*(p). \qquad (4.32)$$

$$\epsilon_r^*(p) = [\mathbf{I} + \mathbf{P}_r^{\text{CSCS}^*}(p) : \delta\mathbf{c}_r^{\text{CSCS}^*}(p)]^{-1} : \{\epsilon_0^*(p) + \dots + \mathbf{P}_r^{\text{CSCS}^*}(p) : [\mathbf{c}_r^*(p) : \epsilon_{nr}^{0*}(p) - \mathbf{C}^{\text{CSCS}^*}(p) : \mathbf{E}_n^{0*}(p)]\}, \\ \langle \epsilon^*(p) \rangle = \mathbf{E}^*(p).$$

When the ellipsoids have the same aspect ratios and orientation, $\mathbf{P}^{\text{CSCS}}$ is independent of the phase and the "normalization term" $\langle [\mathbf{I} + \mathbf{P}^{\text{CSCS}^*}(p) : \delta\mathbf{c}^{\text{CSCS}^*}(p)]^{-1} \rangle^{-1}$ reduces to unity.

An illustrative example of the application of this method is reported on Fig. 25 and 26 for a shear relaxation test performed on the two-phase material which has been already studied in the linear case, but now with the local constitutive equations

$$\dot{e} = a\sigma_{\text{equ}}^{m-1}s + bs\dot{s} \qquad (4.33)$$
with $\sigma_{equ}$ the von Mises equivalent stress; in this particular case, the problem reduces to a scalar one for which we derive the global stress response $S(t)$ to a stepwise strain loading $E(t) = E_0 H(t)$. Note that use has been made at each step of the closed form solution of the linear case reported hereabove but that collocation and FFT techniques are also necessary. The results clearly show both the nonlinearity of the response (Fig. 25) and the variation of the overall eigenstrain rate $\dot{E}_n^0(t, t_n)$ at each step (Fig 26).

<figure><img src="image_27.png" /><figcaption>Fig. 25 Normalized shear stress relaxation for a two-phase material according to the C.S.C.S. ($m = 2, f = .5, b_1 = b_2 = 1, a_1 = 1, a_2 = 5$) (from [44])</figcaption></figure>

<figure><img src="image_28.png" /><figcaption>Fig. 26 Time variation of the overall eigenstrain rate at different steps $t_n$ for the case of Fig. 25 (from [44])</figcaption></figure>
This method can be extended [45] to more general rate-dependent elastoplastic constitutive equations such as the ones which are used frequently in crystalline plasticity, namely

$$
\begin{align}
\dot{\epsilon} &= s: \dot{\sigma} + \sum_s \dot{\gamma}_s R^s, & R^s &= \frac{1}{2} (\bar{n}_s \otimes \bar{m}_s + \bar{m}_s \otimes \bar{n}_s) \nonumber \\
\dot{\gamma}_s &= \dot{\gamma}_s^0 (\frac{\tau_s}{\tau_s^0})^{\nu_s}, & \tau_s &= R^s: \sigma, & \dot{\tau}_s^0 &= \sum_t H^{st} \dot{\gamma}_t \tag{4.34}
\end{align}
$$

where $\bar{n}_s$ and $\bar{m}_s$ are the unit normal and the unit vector of the slip system (s), $\dot{\gamma}_s$ is the shear rate on this system and $H^{st}$ is the so-called "hardening matrix" which may depend on stress and strain. Such constitutive equations can be put in the form:

$$
\begin{align}
\dot{\epsilon} &= s: \dot{\sigma} + g(\sigma, \bar{\alpha}) \tag{4.35} \\
\bar{\alpha} &= \bar{h}(\sigma, \bar{\alpha}), \nonumber
\end{align}
$$

where the internal parameters $\bar{\alpha}$ have been given the form of vectors, with components $(\alpha_l)$
with $I = 1$ to $N$, for the sake of simplicity: they could be the $N$ reference shear stresses $\tau_s^0$ of
(4.34). The same linearization procedure as before has now to be applied to both $\dot{\epsilon}$ and $\bar{\alpha}$,
which yields:

$$
\begin{equation}
\begin{aligned}
\dot{\epsilon}(t) &= s: \dot{\sigma}(t) + g(t_n) + m_{(n)}: [\sigma(t) - \sigma(t_n)] + n_{(n)}: [\bar{\alpha}(t) - \bar{\alpha}(t_n)] + \dots \\
&\quad + [g(t) - g(t_n) - m_{(n)}: [\sigma(t) - \sigma(t_n)] - n_{(n)}: [\bar{\alpha}(t) - \bar{\alpha}(t_n)]] [1 - H(t - t_n)]
\end{aligned}
\tag{4.36}
\end{equation}
$$

$$
\text{with: } \mathbf{m}_{(n)} = \frac{\partial g}{\partial \sigma} [\sigma(t_n), \bar{\alpha}(t_n)], \quad \mathbf{n}_{(n)} = \frac{\partial g}{\partial \bar{\alpha}} [\sigma(t_n), \bar{\alpha}(t_n)]
$$

and

$$
\begin{equation}
\begin{aligned}
\bar{\alpha}(t) &= \bar{h}(t_n) + \mathbf{p}_{(n)}: [\sigma(t) - \sigma(t_n)] + \mathbf{q}_{(n)}: [\bar{\alpha}(t) - \bar{\alpha}(t_n)] + \dots \\
&\quad + \{\bar{h}(t) - \bar{h}(t_n) - \mathbf{p}_{(n)}: [\sigma(t) - \sigma(t_n)] - \mathbf{q}_{(n)}: [\bar{\alpha}(t) - \bar{\alpha}(t_n)]\} [1 - H(t - t_n)]
\end{aligned}
\tag{4.37}
\end{equation}
$$

$$
\text{with: } \mathbf{p}_{(n)} = \frac{\partial \bar{h}}{\partial \sigma} [\sigma(t_n), \bar{\alpha}(t_n)], \quad \mathbf{q}_{(n)} = \frac{\partial \bar{h}}{\partial \bar{\alpha}} [\sigma(t_n), \bar{\alpha}(t_n)].
$$

We then apply the Laplace-Carson transform to these equations and eliminate $\bar{\alpha}^*(p)$ through
the relation:

$$
\begin{align}
\bar{\alpha}^*(p) &= [I - \frac{l}{p} \mathbf{q}_{(n)}]^{-1} \cdot [\frac{l}{p} \mathbf{p}_{(n)}: \sigma^*(p) + \tilde{\alpha}_{(n)}^*(p)] \nonumber \\
\text{with:} \quad \tilde{\alpha}_{(n)}^*(p) &= \frac{l}{p} [\bar{h}(t_n) - \mathbf{p}_{(n)}: \sigma(t_n) - \mathbf{q}_{(n)}: \bar{\alpha}(t_n)] + \dots \tag{4.38}
\end{align}
$$

$$
\dots + \int_0^{t_n} \{[\bar{h}(t) - \bar{h}(t_n) - \mathbf{p}_{(n)}: [\sigma(t) - \sigma(t_n)] - \mathbf{q}_{(n)}: [\bar{\alpha}(t) - \bar{\alpha}(t_n)]] \exp(-pt) dt.
$$

where *I* is the second order unit tensor in a *NxN* space. This leads to the final relation:
$$\varepsilon^*(p) = s^*(p) : \sigma^*(p) + \varepsilon_{(n)}^0(p)$$

$$\text{with: } s^*(p) = s + \frac{1}{p} \mathbf{m}_{(n)} + \frac{1}{p} \mathbf{n}_{(n)} \cdot (p\mathbf{I} - \mathbf{q}_{(n)})^{-1} \cdot \mathbf{p}_{(n)}$$

$$\begin{align}
\varepsilon_{(n)}^0(p) &= \mathbf{n}_{(n)} \cdot (p\mathbf{I} - \mathbf{q}_{(n)})^{-1} \cdot \tilde{\alpha}_{(n)}^*(p) + \tilde{\varepsilon}_{(n)}^*(p) \tag{4.39} \\
\tilde{\varepsilon}_{(n)}^*(p) &= \frac{1}{p} [\mathbf{g}(t_n) - \mathbf{m}_{(n)} : \sigma(t_n) - \mathbf{n}_{(n)} \cdot \tilde{\alpha}(t_n)] + \dots \nonumber
\end{align}$$

$$\dots + \int_0^{t_n} \{ \mathbf{g}(t) - \mathbf{g}(t_n) - \mathbf{m}_{(n)} : [\sigma(t) - \sigma(t_n)] - \mathbf{n}_{(n)} \cdot [\tilde{\alpha}(t) - \tilde{\alpha}(t_n)] \} \exp(-pt) dt.$$

So, with another definition of $s^*(p)$ and $\varepsilon_{(n)}^0(p)$, the nature of the local constitutive equations is the same as before: the auxiliary strain function $\varepsilon_{(n)}^0(t)$ can be considered as a true eigenstrain and the homogenization treatment is unchanged.

To illustrate the proposed method, we go back to the same simple example except that the number of phases is arbitrary and that the local constitutive equations depend on one scalar internal parameter [45]:

$$\begin{align}
\dot{\mathbf{e}} &= a\dot{s} + b\left(\frac{\sigma_{\text{equ}}}{\alpha}\right)^{m-1}s \tag{4.40} \\
\dot{\alpha} &= h\left(\frac{\sigma_{\text{equ}}}{\alpha}\right)^{\nu} \nonumber
\end{align}$$

where $a, b, h, m$ and $\nu$ are material constants which may differ from phase to phase. For a proportional loading, the problem reduces to a scalar one: we look for some macroscopic stress relaxation function $\mu^{\text{SCSC}}(t, t_n)$ and the macroscopic eigenstrain $E_{(n)}^0(t, t_n)$ to be defined at each step $t_n$. The overall stress relaxation function $\mu^{\text{CSCS}}$ at time $t_n$ can be derived by an iterative scheme with the recursive formula:

$$\begin{gather}
\mu_{(i+1)}^{\text{CSCS}*} = \left( \sum_r \frac{f_r \mu_r^*}{3\mu_{(i)}^{\text{CSCS}*} + 2\mu_r^*} \right) \cdot \left( \sum_r \frac{f_r}{3\mu_{(i)}^{\text{CSCS}*} + 2\mu_r^*} \right)^{-1} \tag{4.41} \\
\text{with: } 2\mu_r^*(p) = \frac{p}{a_r(p+1/T_r)}, \quad T_r = \frac{a_r}{b_r} \nonumber
\end{gather}$$

whereas the interaction law reads

$$\begin{gather}
s_r^* = \frac{5\mu_r^*}{3\mu^{\text{CSCS}*}} S^* + \frac{6\mu_r^*\mu^{\text{SC}*}}{3\mu^{\text{CSCS}*} + 2\mu_r^*} (E_{(n)}^0 - e_{(n)r}^0)^* \tag{4.42} \\
\text{with: } E_{(n)}^0 = \sum_r c_r B_r^* e_{(n)r}^0, \quad B_r^* = \frac{5\mu_r^*}{3\mu^{\text{CSCS}*} + 2\mu_r^*} \nonumber
\end{gather}$$
Typical results are given on Figure 27 for a shear test with the prescribed constant shear strain rate $\dot{E}_0$. The normalized stress response $S/\dot{E}_0$ depends clearly on the strain rate amplitude, as expected from the nonlinear local behaviour.

Similar computations applied to rate-dependent elastoplastic polycrystals are now in progress. Extensions are developed too when instantaneous plasticity is simultaneously present. We conclude by remarking that the above proposed formulation of the rate-dependent elastoplastic self-consistent scheme reduces, when elasticity is neglected, to the one already proposed (at finite strain) by Molinari et al [46] for viscoplastic polycrystals; whereas this treatment does not coincide with Hutchinson's [39] one, it is in better agreement with Ponte Castañeda's recent one [47].

<figure><img src="image_29.png" /><figcaption>Fig. 27: normalized stress evolution at constant strain rate "edot" (time unit arbitrary) (2 phases; $m_r = v_r = 2$; $a_r = \alpha_r(0) = 1$; $f = .5$; $b_l = h_l = 1$; $b_2 = h_2 = 5$) (from [45]).</figcaption></figure>

### 4.3.5 Going back to Hill's assumption

Hill's treatment of the nonlinear classical self-consistent scheme, including its former extension to viscoplasticity and its latest developments for rate-dependent elastoplasticity, can be considered as the most advanced approach to polycrystal plasticity. Nevertheless, one must keep in mind that it is still an approximate one which cannot be thought as fully self-consistent. Replacing the original Kröner's assumption of piecewise uniform plastic straining for the solution of the reference inclusion problem by the one of piecewise uniform instantaneous moduli (or secant ones, or the correspondant Laplace-transformed viscoelastic moduli for rate-dependent plasticity) allows Hill's treatment to take far better nonlinearity into account, but this is not beyond reproach: when we have to solve any given boundary-value problem in elastoplasticity, we know that the instantaneous (or secant) moduli vary from point to point according to the local stress and strain state...

There is some consistency in considering the ellipsoidal inclusion tangent (or secant) moduli as uniform when the matrix ones have been considered so, since the mechanical state
of the inclusion is homogeneous in that case. But the reverse is not true: even if the inclusion
is uniformly strained and stressed, the matrix cannot be so, neither can its moduli. And
conversely, if the matrix has not uniform moduli, the inclusion, even if ellipsoidal, is not
uniformly deformed and so its moduli are no more uniform.

A fully self-consistent treatment would have needed to embed the inclusion in the Homogeneous Equivalent Medium whose constitutive behaviour, as resulting from self-consistent modelling, should make its instantaneous (or secant) elastoplastic moduli vary from point to point in any non trivial B.V.P. problem (and so differently according to the inclusion under consideration and to its current mechanical state too...). As we can appreciate, such a treatment would need at each step a tremendous number of overlapping iterative procedures which, despite the exponential increase in time of the numerical potential and computing power, still seems disheartening.

Nevertheless, we have to make progress in this direction and imagine better and better approximations. One practical objective could be to force an improved nonlinear self-consistent approximation to abide by bounds in the simple problems for which it has been proved to violate them [48]; such improved "variational estimates" could be derived from bounding approaches by building them up from the self-consistent elastic ones [49, 47]. This suggests that nonlinear improved estimates would take advantage of improved nonlinear bounding theories (and vice versa ?).

## 4.4 Nonlinear generalized self-consistent modelling

When dealing with nonlinear generalized self-consistent model, the foregoing problem becomes quite... generalized: even if Hill's approximation of uniform moduli in the matrix is adopted, it can obviously not be so in the composite inclusion. So, the problem must be faced. Several approximations have been proposed: many of them relieve on the use of plasticity deformation theory and secant elastoplastic moduli so as to preserve an isotropic formulation and to avoid the additional complexity attached to anisotropy. Whether the deformation or the flow theory is used, it can be combined with the approximation of stepwise uniform elastoplastic moduli in phase subdomains as small as possible.

Let us quote for instance the approximate treatment of the elastoplastic three-phase model [23, 50, 51] by decomposing the core and the shell into as many as wished concentric (or coaxial) spherical (or cylindrical) shells inside which the secant moduli are assumed uniform: use is then made of the solution [24, 25] of the elastic so-called "n-phase spherical (or cylindrical) model" with different possible definitions of the secant moduli [50, 51, 52]. In such treatments, the stress or strain definition of the three-phase model [23] has rather to be used instead of the energetical one. Finite element numerical computations may be developed as well [50]. The main difficulty indeed is the same as in the case of the classical self-consistent scheme and is concerned with the way to take the inhomogeneity of the matrix into account since it is made of the (a priori unknown) H.E.M.. Iterative procedures dealing with successive known approximations of the H.E.M. could be a possible approach to this additional difficulty.

# 5. CONCLUSION

To sum up, we can appreciate that the integration of structural morphology into micromechanical approaches is still a widely open question. Despite some progress has been achieved in this field for the last thirty years, many questions still remain unsolved. The future expected developments are concerned both with improved treatments of nonlinear behaviour and with efficient practical means to quantify experimentally and to represent conveniently the main morphological characteristics of inhomogeneous material as well as their evolution under straining.
What is at stake with these efforts is not so much to be able to predict better the macroscopic behaviour of given materials as to do so at the mesoscale: several studies [e.g. 50] have shown that morphology plays a major role at this mesoscale in inducing local stress or strains concentrations which rule damage initiation and evolution, so that one can guess that any damage micromechanical theory will have to take to the best morphological parameters into account. In addition, invention of new materials as well improvement of the existing ones in view of optimized properties for a given use obviously have to get morphological aspects under control.

## Acknowledgements

Michel Bornert and Claude Stolz are friendly and gratefully acknowledged for their essential contribution to most of the above reported results. Many thanks are especially due too, between many others, to Éveline Hervé and Yves Rougier for their kind and efficient cooperation.

## REFERENCES

1. Eshelby, J.D.: The determination of the elastic field of an ellipsoidal inclusion and related problems, Proc. Roy. Soc. London, A 241 (1957), 376-396.

2. Kröner, E.: Bounds for effective elastic moduli of disordered materials, J. Mech. Phys. Solids, 25 (1977), 137-155.

3. Hashin, Z.: The elastic moduli of heterogeneous materials, J. Appl. Mech, 29 (1962), 143-150.

4. Christensen, R.M.: Mechanics of composite materials, Wiley-Interscience, New-York 1979.

5. Stolz, C. and Zaoui, A.: Analyse morphologique et approches variationnelles du comportement d'un milieu élastique hétérogène, C. R. Acad. Sci. Paris, II 312 (1991), 143-150.

6. Bornert, M., Stolz C. and Zaoui A.: Morphologically representative pattern-based bounding in elasticity, J. Mech. Phys. Solids, 44 (1996), 307-331.

7. Bornert, M.: A generalized pattern-based self-consistent scheme, Comp. Mat. Science, 5 (1996), 17-31.

8. Rougier, Y., Stolz, C. and Zaoui, A.: Représentation spectrale en viscoélasticité linéaire des matériaux hétérogènes, C. R. Acad. Sci. Paris, II 316 (1993), 1517-1522.

9. Hill, R.: Continuum micro-mechanics of elastoplastic polycrystals, J. Mech. Phys. Solids, 13 (1965), 89-101.

10. Rougier, Y., Stolz C. and Zaoui, A.: Self-consistent modelling of elastic-viscoplastic polycrystals, C. R. Acad. Sci. Paris, II 318 (1994), 145-151.

11. Mori, T. and Tanaka, K.: Average stress in the matrix and average elastic energy of materials with misfitting inclusions, Acta Metall., 21 (1973), 571-574.

12. Hershey, A.V.: The elasticity of an isotropic aggregate of anisotropic cubic crystals, J. Appl. Mech., 21 (1954), 236-240.

13. Kröner, E.: Berechnung der elastischen Konstanten des Vielkristalls aus den Konstanten des Einskristalls, Z. Phys., 151(1958), 504-518.

14. Walpole, L.: On bounds for the overall elastic moduli of inhomogeneous systems - I, J. Mech. Phys. Solids, 14 (1966), 151-162.

15. Hashin, Z. and Shtrikman, S.: A variational approach to the theory of the elastic behaviour of multiphase materials, J. Mech. Phys. Solids, 11 (1963), 127-140.

16. Kröner, E.: Self-consistent scheme and graded disorder in polycrystal elasticity, J. Phys. F., 8 (1978), 2261-2267.
17. Kröner, E.: Graded and perfect disorder in random media elasticity, J. Engng. Mech. Div., 106-EM5 (1980), 889-914.

18. Zimmermann, R.W.: Elastic moduli of a solid containing spherical inclusions, Mech. Mater., 12 (1991), 17-24.

19. Christensen, R.M.: A critical evaluation for a class of micro-mechanical models, J. Mech. Phys. Solids, 38 (1990), 379-404.

20. Christensen, R.M. and Lo, K.H.: Solutions for effective shear properties in three phase sphere and cylinder models, J. Mech. Phys. Solids, 27 (1979), 315-330.

21. Love, A.E.H.: A treatise on the mathematical theory of elasticity, Dover Publications, New York 1944.

22. Eshelby, J.D.: The continuum theory of lattice defects, in: Progress in Solid State Physics 3 (Eds. F. Seitz and D. Turnbull), Academic Press, New York 1956, 79-144.

23. Hervé, É. and Zaoui, A.: Modelling the effective behaviour of nonlinear matrix / inclusion composites, Eur. J. Mech. A / Solids, 9 (1990), 505-515.

24. Hervé, É. and Zaoui, A.: *n*-layered inclusion-based micromechanical modelling, Int. J. Engng Sci., 31 (1993), 1-10.

25. Hervé, É. and Zaoui, A.: Elastic behaviour of multiply coated fibre-reinforced composites, Int. J. Engng Sci., 33 (1995), 1419-1433.

26. Hashin, Z. and Rosen, B.W.: The elastic moduli of fiber-reinforced materials, J. Appl. Mech., 31 (1964), 223-232.

27. Hervé, É., Stolz, C. and Zaoui, A.: À propos de l'assemblage de sphères composites de Hashin, C. R. Acad. Sci. Paris, II 313 (1991), 857-862.

28. Milton, G. and Kohn, R.: Variational bounds on the effective moduli of anisotropic composites, J. Mech. Phys. Solids, 36 (1988), 597-629.

29. Bornert, M.: Ph. D. Dissertation, E.N.P.C., Paris, 1996.

30. Willis, J.R.: Bounds and self-consistent estimates for the overall properties of anisotropic composites, J. Mech. Phys. Solids, 17 (1977), 185-202.

31. Ponte Casteñeda, P. and Willis, J.R.: The effect of spatial distribution on the effective behaviour of composite materials and cracked media, J. Mech. Phys. Solids, 43 (1995), 1919-1951.

32. Suquet, P.: Elements of homogenization for inelastic solid mechanics, in: Homogenization Tecniques for composite media (Eds. E. Sanchez-Palencia and A. Zaoui), Lecture Notes in Physics 272, Springer-Verlag, Berlin 1985, 193-278.

33. Mandel, J.: Mécanique des milieux continus, Gauthier-Villars, Paris 1966.

34. Laws, N. and McLaughlin, R.: Self-consistent estimates for the viscoelastic creep compliances of composite materials, Proc. R. Soc. Lond., A 359 (1978), 251-273.

35. Rougier, Y.: Ph. D. Dissertation, Ecole Polytechnique, Palaiseau, 1994.

36. Beurthey, S.: private communication, 1995.

37. Kröner, E.: Zur plastischen Verformung des Vielkristalls, Acta Metall, 9 (1961), 155-161.

38. Berveiller, M. and Zaoui, A.: An extension of the self-consistent scheme to plastically flowing polycrystals, J. Mech. Phys. Solids, 26 (1979), 325-344.

39. Hutchinson, J.W.: Bounds and self-consistent estimates for creep of polycrystalline materials, Proc. Roy. Soc. London, A 348 (1976), 101-127.

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

41. Harren, S.V.: The finite deformation of rate-dependent polycrystals, J. Mech. Phys. Solids, 39 (1991), 345-383.

42. Nemat-Nasser, S. and Obata, M.: Rate-dependent finite elastoplastic deformation of polycrystals, Proc. Roy. Soc. London, A 407 (1986), 343-375.
43. Zaoui, A. and Raphanel, J.L.: On the nature of the intergranular accommodation in the modelling of elastoviscoplastic behavior of polycrystalline aggregates, in: Large plastic deformations, Fundamentals and applications to metal forming (Eds. C. Teodosiu, J.L. Raphanel and F. Sidoroff), Balkema, Rotterdam 1993, 185-192.

44. Zaoui, A., Rougier, Y. and Stolz, C.: Micromechanical modelling based on morphological analysis; Application to viscoelasticity, in: Microstructure-Property Interactions in Composite Materials (Ed. R. Pyrz), Kluwer Academic Publishers, Dordrecht 1995, 419-430.

45. Navidi, P., Rougier, Y. and Zaoui, A.: Self-consistent modelling of elastic-viscoplastic multiphase materials, in: Micromechanics of plasticity and damage of multiphase materials (Eds. A. Pineau and A. Zaoui), Kluwer Academic Pub., Dordrecht 1996, 123-130.

46. Molinari, A., Canova, G.R. and Ahzi, S.: A self-consistent approach of the large deformation polycrystal viscoplasticity, Acta Metall. Mater., 35 (1987), 2983-2994.

47. Ponte-Castañeda, P.: Exact second-order estimates for the effective mechanical properties of nonlinear composite materials, J. Mech. Phys. Solids, 44 (1996), 827-862.

48. Gilormini, P.: A critical evaluation for various nonlinear extensions of the self-consistent model, in: Micromechanics of plasticity and damage of multiphase materials (Eds. A. Pineau and A. Zaoui), Kluwer Academic Pub., Dordrecht 1996, 67-74.

49. Suquet, P.: Overall properties of nonlinear composites; remarks on secant and incremental formulations, in: Micromechanics of plasticity and damage of multiphase materials (Eds. A. Pineau and A. Zaoui), Kluwer Academic Pub., Dordrecht 1996, 149-156.

50. Bornert, M., Hervé, É., Stolz, C. and Zaoui, A.: Self-consistent approaches and strain heterogeneities in two-phase elastoplastic materials, Appl. Mech. Rev., 47 (1994), S.66-76.

51. Bornert, M.: Morphological effects at the local scale in two-phase materials; some experimental evidence and tentative modelling, in: Micromechanics of plasticity and damage of multiphase materials (Eds. A. Pineau and A. Zaoui), Kluwer Academic Pub., Dordrecht 1996, 27-34.

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