Intelligent modeling of materials

Joseph Zarka, Pirouz Navidi

► To cite this version:

Joseph Zarka, Pirouz Navidi. Intelligent modeling of materials. Mechanics of Materials, 1998, 28 (1-4), pp.61-82. 10.1016/S0167-6636(97)00063-X. hal-00111602

HAL Id: hal-00111602

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

Submitted on 21 Mar 2023

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

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

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

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
Intelligent modeling of materials

Joseph Zarka *,<fnref n="1" />, Pirouz Navidi

Laboratoire de Mécanique des Solides, Ecole Polytechnique, F-91128 Palaiseau Cedex, France

**Abstract**

Over several decades, many contributions have been made to understand the physics and to develop models of aggregates of metals and composites. Three books 'Modelling small deformations of polycrystals' (eds. J. Gittus and J. Zarka) and 'Modelling large deformations of solids' (eds. J. Gittus, J. Zarka and S. Nemat-Nasser) and 'New approach of inelastic analysis of structures' (eds. J. Zarka et al.) gave already a general view of the state of the art. Nemat-Nasser, particularly in his impressive book 'Micromechanics: overall properties of heterogeneous materials', improved and is still greatly improving our knowledge in this fundamental area. In this paper, which is dedicated to his 60th anniversary, we want to give again and for the last time, our own view of the problem based on (may be) not well known works and also some unpublished theses of some students which were done between 1964 and now. In the first part, we shall underline what the physical quantitative description of inelastic behavior of single crystals means for us by answering the following questions: (1) are the various glides on each crystal slip plane, as introduced in the multiple plastic potential theory, potential internal variables? (no); (2) can we build a real physical model for the crystal in correlation with the observations of the evolution of the crystalline defects, such as the dislocations, the vacancies, the precipitates ... and the thermally activated processes? (yes); (3) can we define the physical hardening parameters (yes) and is the latent hardening higher than the self-hardening? (no); (4) is it possible to reach a description which fits the experimental results? (may be!). Then, in the second part, assuming a given behavior for the crystals, for the global inelastic behavior of polycrystals, we shall try to answer the questions: (5) how many different orientations of crystals are necessary to reach the description of an initially isotropic polycrystal? (three well selected ones!); (6) are the simplified models, such as the self-consistent model, sufficient to represent the interactions between crystals? (yes); (7) are the actual numerical simulations based on rate-dependant or rate-independent crystal plasticity, useful for the industrial applications? (no, unhappily!). Finally, in the third part, we shall give our new general framework of 'intelligent' modelling of aggregates where we take into account not only the local behaviors of the crystals but also their size, their shape and their relative distribution. In this framework, it is needed: (i) To build a data base i.e. to obtain some experimental, real or simulated, results where the experts identify all variables or descriptors which may relevant to the given problem. This is, at first, done with some primitive descriptors x which are usually in a limited number. Then, the data are transformed into intelligent descriptors XX in a larger number, using the existing knowledge and theories which are still always insufficient. The descriptors may be numbers, boolean, alphanumeric, name of files which gives access to data bases, or treatments of curves, signals and images. The results or conclusions may be classes (good, not good, ...) or numbers (Young modulus, cost, weight, life time, ...). Usually, the data base may contain roughly 30 to 150 examples with 10 to 1000 descriptors and 1 to 20 conclusions. (ii) To generate the rules with any

<footnote>* Corresponding author. E-mail: zarka@athena.polytechnique.fr.</footnote>

<footnote n="1">¹ Also at: Centre of Excellence for Advanced Materials, Department of Applied Mechanics and Engineering Sciences, University of California, San Diego, 9500 Gilman Drive, La Jolla, CA 92093-0416, USA. E-mail: zarka@ceam.ucsd.edu.</footnote>
automatic learning tool. The intelligent descriptors help these learning algorithms. Each conclusion is explained as function or set of rules of some among the input intelligent descriptors with a known reliability or accuracy. If this reliability is too low, it implies that either there is not enough data or there are bad, missing descriptors or the problem was not well described. (iii) To optimize at two levels (inverse problems): (1) Considering the intelligent descriptors as independent; it is possible to get the optimal solution satisfying the special required properties and allowing the discovery of new mechanisms; (2) considering the intelligent descriptors linked to primitive descriptors, it is possible to obtain the optimal solution which is technologically realizable. So, not only a practical optimal solution is obtained, but also the experts may learn the missing parts, may build models or theories based only on the retained intelligent descriptors and guided by the structures of the rules or relationships. We shall illustrate our framework by treating step by step the problem of the global elastic behavior of concrete with a cement matrix and soft or hard inclusions. This problem has been the subject of several papers, providing simplified models, simple bounds of the elastic moduli and sophisticated theories ... were produced. Our aim is to give the results to the engineers in such a practical way that they could: (1) estimate the elastic properties within a few percent of error for any concentration, any shape, or distribution of the inclusions;(2) select the preparation of the aggregate to reach any pre-assigned elastic properties even at the lowest cost or weight).

# 1. Physical modeling of metallic single crystals

## 1.1. Theoretical model (Fig. 1)

It is well known that experiments have shown that:

(1) Irreversible plastic strains, in metallic single crystals, are induced by the motions of some structural defects; for example, during the slip mode, which is the most important one, the crystal structure characterized by a lattice, remains unchanged and the strains are delayed (viscous) and often very large.

(2) Reversible elastic strains are linked to the mean lattice; the strains are instantaneous and generally very small.

(3) This fundamental fact has allowed to differentiate the elastic strains and plastic strains at any time by introducing two families of triads, $\vec{F}_K$ and $\vec{e}_\alpha$ which are defined with the following properties:

* $\vec{F}_K$ are three macroscopic vectors which characterize the total transformation.

* $\vec{e}_\alpha$ are three microscopic vectors linked to the mean lattice which characterize the elastic deformation.

<figure><img src="image_2.png" /><figcaption>Fig. 1. Different states of the crystal.</figcaption></figure>
• At time $t=0$, chosen as an initial unloaded and natural state, relative to a fixed Cartesian frame of reference $\vec{g}_i$ and at the point $M(X^K \delta_i^K)$, the components of these vectors are $\vec{F}_K(0) = \delta_K^i \vec{g}_i$ and $\vec{e}_\alpha(0) = \delta_\alpha^i \vec{g}_i$, where $\delta_i^K$, $\delta_K^i$ and $\delta_\alpha^i$ stand for the Kronecker symbol.

• At the actual time $t$ and at the point $m(x^i)$, the image of $M$ in the total deformation, the components of the above-mentioned vectors are $\vec{F}_K = (\partial x^i / \partial X^K) \vec{g}_i$ and $\vec{e}_\alpha = x_\alpha^i \vec{g}_i$.

• $F = (\partial x^i / \partial X^K)$ is the classical deformation gradient; $e = (2 e_{ij} = \delta_{ij} - \delta_{KL} (\partial X^K / \partial x_i) (\partial X^L / \partial x_j))$ is the total strain tensor which, similar to the total rotation, is generally large.

• $E = (x_\alpha^i)$ gives the elastic deformation; $e' = (2 e'_{ij} = \delta_{ij} - \delta_{\alpha\beta} x_i^\alpha x_j^\beta)$ is the small elastic strain tensor.

• The elastic strains may be instantaneously recovered by a total unloading (the temperature being here assumed constant) and with the help of a convenient rotation to reach an intermediate unloaded state in which the vectors of the microscopic triad $\vec{e}_\alpha$ are equipotent to the initial $\vec{e}_\alpha(0)$ or $\vec{g}_i$ but in which there is a new distribution of defects relative to the initial state; let us define by $\{y\}$ a family of internal parameters, which describes this actual distribution of the defects.

• By the classical elastic theory, it may be possible to calculate the global stress tensor $\Pi$ ($\Pi^{\alpha\beta}$) (referred to the frame $\vec{e}_\alpha$) or the Cauchy stress tensor,

$$ \sigma \left( \sigma^{ij} = \frac{1}{\det(x_\alpha^i) x_\alpha^j} \Pi^{\alpha\beta} \right). $$

• When we impose $e'$ for the actual given defect distribution $\{y\}$, since in the presence of the defects, the crystal is an inhomogeneous medium.

• The total deformation can not be characterized so easily; it is necessary to follow the crystal step by step and to split the velocity gradient $v_j^i$ into the elastic part $v_j^{ii}$ and the plastic part $v_j^{ii}$: $v_j^i = v_j^{ii} + v_j^{iii}$ which gives $\partial \dot{x}^i / \partial X^K = v_j^i \partial x^j / \partial X^K$ and $\dot{x}_\alpha^i = v_j^{ii} x_\alpha^j$.

## 1.2. Physical model

Indeed, it is necessary to build a physical model to get the elastic and plastic parts of the velocity gradient, $v_j^{ii}$ and $v_j^{iii}$, explicitly. Here we only recall how we obtain the plastic part:

• We know that for metallic single crystals, the plastic strains result from slips on some crystallographic planes and along some crystallographic directions, these planes and directions being fully determined in the frame $\vec{e}_\alpha$ which represents the mean lattice.

• We know that if we denote by $\vec{h}$ ($h^i$), the unit vector along the slip direction, by $\vec{n}$ ($n_j$), the unit normal vector to a slip plane which contains atoms, and by $\gamma$ the relative velocity along $\vec{h}$ between two planes which are unit length apart, the velocity gradient is then given by:

$$ v_j^{ii} = h^i n_j \gamma $$

and that when there are several such active slip systems ($r$), the velocity gradient is expressed by:

$$ v_j^{iii} = \sum_{(r)} h^i(r) n_j^{(r)} \gamma^{(r)}. $$

• During this transformation, the frame $\vec{e}_\alpha$ remains fixed and the frame $\vec{F}_K$ are transformed (strain and rotation); the vector (Burgers vector) which links two atoms, will have constant components in the frame $\vec{e}_\alpha$,

$$ \vec{b} = b^\alpha \vec{e}_\alpha = b^\alpha x_\alpha^i \vec{g}_i \Rightarrow \vec{h} = \frac{\vec{b}}{b} = h^i \vec{g}_i, $$
the normal $\vec{N}$ to a plane which contains atoms, will have constant covariant components in the frame $\vec{e}^{\alpha}$,

$$ (x_{\alpha}^{i})^{-1} \vec{g}_{i}, \quad \vec{N} = N_{\alpha} \vec{e}^{\alpha} = N_{i} \vec{g}_{i} \Rightarrow \vec{n} = \frac{\vec{N}}{N}. $$

• But we need to know when a slip system (r) is active and what is the amplitude of the relative velocity $\gamma^{(r)}$ (and here is the starting point of all the various proposals, each person claiming that he has the best model as do the authors of this paper); usually, based on Schmid's law, symbolically, we introduce

$$ \tau^{(r)} \equiv h_{i}^{(r)} n_{j}^{(r)} \boldsymbol{\sigma}^{ij}; \quad f^{(r)} \equiv \tau^{(r)} - g^{(r)}(\{y\}) \geq 0 $$

and in view of the fact that $\gamma^{(r)}$ results from delayed motions of the defects, we assume $\gamma^{(r)} \equiv \gamma^{(r)}(\tau^{(k)}, \{y\})Y(f^{(r)})$ where Y is the Heaviside function; during these motions, some parameters among the $\{y\}$ family may evolve,

$$ \{\dot{y}\} \equiv \sum_{(r)} G^{(r)}(\tau^{(k)}, \{y\}) Y(f^{(r)}). $$

## 1.3. Face-centered cubic single crystal

### 1.3.1. Simplified model for room temperature

In Zarka (1968), we proposed a model where we took into account the dislocation distribution based on a schematic representation (Fig. 2) of the classical electron microscope images of the observed microscopic physical phenomena.

At room temperature and for quasi-static loadings, an idealized model was proposed as follows: We defined a mechanism (r) of dislocations by grouping the set of randomly distributed small segments of dislocations which have the same Burgers vector $\vec{b}^{(r)}$ and parallel slip planes with unit normal $\vec{n}^{(r)}$, $\mathcal{N}^{(k)}$, $l^{(k)}$ are respectively the mean number per unit volume of segments and their mean length (the product $\mathcal{N}^{(k)}l^{(k)}$ is the total length of such dislocations per unit volume which means a quantity which may be measured at any time).

<figure><img src="image_3.png" /><figcaption>Fig. 2. Crude schematization of the distribution of dislocations.</figcaption></figure>
Then by using simple hypothesis (as initiated by Taylor, Stroh, Saada, Cambell, ...), we were able:

• To compute the interactions between mechanisms of dislocations as:

$$ \mathcal{F}^{(r)} = \mathcal{F}_{PN}^{(r)} + \mu b^{(r)} \sum_{(k)} \beta_{(k)}^{(r)} (\mathcal{N}^{(k)} l^{(k)})^{1/2}, $$

where $\mu$ is the elastic shear modulus and $\beta_{(k)}^{(r)}$ was explicitly obtained after very long calculations, to be

$$ \beta_{(k)}^{(r)} = \begin{cases} \beta_{(k)}^{(k)} = 1/8, & \text{for self-hardening,} \\ 1/16, & \text{for } (r) \neq (k) \text{ but with parallel slip plane,} \\ 1/20, & \text{for } (r) \neq (k) \text{ but with same Burgers vector,} \\ 1/12, & \text{for all other mechanisms } (r) \neq (k). \end{cases} $$

(It was possible to provide a physical basis for this formula. The formula, however, may be taken as a phenomenological extension of the Taylor's model.) The formula shows that (what we call) the self-hardening is more important than the latent hardening i.e. it is easier, after the activation of the principal slip mechanism to move a dislocation which belongs to another slip system. Note that experiments (Nowacki and Zarka, 1971), seemed to confirm this point, although many other researchers have reported contrary results!

$F^{(r)} = b^{(r)}\tau^{(r)} \equiv b^{(r)}h_i^{(r)}n_j^{(r)}\sigma^{ij}$ is the applied force per unit length on the dislocation (Peach-Koelher formula); $F^{(r)}$ must be greater or equal $\mathcal{F}^{(r)}$ in order to move the dislocation, i.e., to activate the mechanism $(r)$.

• To compute some mean times of realization (such as, for example, the mean time to cross over another dislocation) which were written in the form $t_m = t_m(\langle F^{(r)} - \mathcal{F}^{(r)} \rangle, F^{(k)}, \mathcal{N}^{(k)}, l^{(k)})$ ($\langle \rangle$ being the Heaviside function).

• To give the full incremental evolution of the crystal deformation using the following hypothesis during the small increment of time $dt$,

(i) for each active mechanism $(r)$, some $(r)$ dislocations are sources and they emit $d\mathcal{N}^{(r)}$ loops of dislocations per unit volume; the loops with the mean length, $2\pi\mathcal{L}_{(r)}$, sweep the mean area $\pi\mathcal{L}_{(r)}^2$ by easy glide before being stopped on the dislocations which pierce their slip plane;

(ii) with the aid of some probability laws, it was possible to evaluate $\pi\mathcal{L}_{(r)}^2$, $d\mathcal{N}^{(r)}$ as functions of the dislocation distribution and the applied forces $F^{(k)}$:

$$ \gamma^{(r)} dt = b d\mathcal{N}^{(r)} \pi \mathcal{L}_{(r)}^2, $$

$$ d(\mathcal{N}^{(r)} l^{(r)}) = d\mathcal{N}^{(r)} 2\pi \mathcal{L}_{(r)}, $$

$$ \pi \mathcal{L}_{(r)}^2 \approx \frac{x^{(r)}}{\sum_{(k)} n^{(k)}} \equiv \frac{x^{(r)}}{\underline{\rho}^{(r)}}, $$

$$ d\mathcal{N}^{(r)} \approx w^{(r)} \mathcal{N}^{(r)} \frac{dt}{x^{(r)} t_{\text{GD}}^{(r)}}, $$

$$ t_{\text{GD}}^{(r)} = \frac{e_0^{(r)}}{\langle F^{(r)} - \mathcal{F}^{(r)} \rangle}, $$

where $\underline{\rho}_{(r)} = \sum_{(k)} n^{(k)}$ is the density per unit surface of all the dislocations which pierce the slip plane of $(r)$, $x^{(r)}$ is the number of jogs on the new loop, and $w^{(r)}$ is the fraction of dislocations which may be sources and may be taken as constant (there is no cross-slip or climb in this approximation).
### 1.3.2. Viscoplastic Potential (Rice, 1970)

By examining the above formulae, we obtain for the glide velocity relative to each mechanism:

$$ \gamma^{(r)} \approx b \frac{w^{(r)} \mathcal{N}^{(r)}}{C_0^{(r)} \underline{\rho}^{(r)}} \langle F^{(r)} - \mathcal{F}^{(r)} \rangle, $$

which expresses that this velocity is, by the intermediacy of $F^{(r)}$, only a function of the applied shear stress, $\tau^{(r)} = h_i^{(r)} n_j^{(r)} \sigma^{ij}$ on the mechanism $(r)$ (the mean pressure $\sigma_i^i$ does not take part when $\vec{h}$ and $\vec{n}$ are orthogonal). Thus, as shown by Rice (1970), by putting for each $(r)$,

$$ \mathcal{F}^{(r)}(\sigma, \mathcal{N}^{(k)}, l^{(k)}) \equiv \frac{w^{(r)} \mathcal{N}^{(r)}}{2 C_0^{(r)} \underline{\rho}^{(r)}} \langle F^{(r)} - \mathcal{F}^{(r)} \rangle^2 $$

and for the crystal

$$ \mathcal{F} = \sum_{(r)} \mathcal{F}^{(r)}, $$

the plastic part of the velocity gradient can be written

$$ v_j^{\prime i} = \left( \frac{\partial \mathcal{F}}{\partial \sigma^{ij}} \right)_{\{y\} = \text{const}} = \sum_{(r)} \left( \frac{\partial \mathcal{F}^{(r)}}{\partial \sigma^{ij}} \right) = \sum_{(r)} \left( \frac{\partial \mathcal{F}^{(r)}}{\partial F^r} \right) \left( \frac{\partial F^{(r)}}{\partial \sigma^{ij}} \right) = \sum_{(r)} h^{i(r)} n_j^{(r)} \gamma^{(r)}. $$

(The virtual variations of the Cauchy stress tensor $\sigma$ are taken in $R^9$, and not in $R^6$, although the real stresses are in $R^6$ due to the symmetry $\sigma^{ij} = \sigma^{ji}$.) $\mathcal{F}$ is the Rice's viscoplastic potential of the single crystal.

In the virtual generalized stress space ($R^9$), the equipotential surfaces $\mathcal{F} = \text{const} > 0$ are regular (without any singular point) and the surface $\mathcal{F} = 0$ corresponds to the elastic surface since we have $F^{(r)} \le \mathcal{F}^{(r)}$. In this space, the vector $v_j^{\prime i}$ is always normal to the equipotential surface which is defined by the actual value of $\sigma$ (Fig. 3).

### 1.3.3. Multiple plastic potential

When all the ratios $w^{(r)}/C_0^{(r)}$ are very large and the strain rates are rather small, the differences $\langle F^{(r)} - \mathcal{F}^{(r)} \rangle$ for the active systems are thus always very small. This means that, at any time, the actual stress state $\sigma$ is very near the elastic limit surface and that the viscosity may be neglected (classical plasticity). For each active system $(r)$ between $t$ and $t+dt$, we have: $F^{(r)}(t) \approx \mathcal{F}^{(r)}(t)$ for any active $(r)$ or $F^{(r)}(t) \le \mathcal{F}^{(r)}(t)$

<figure><img src="image_4.png" /><figcaption>Fig. 3. Schematization of the viscoplastic potential surfaces.</figcaption></figure>
for all (r) then, at time t, $F^{(r)} \simeq \mathcal{F}^{(r)}(t)$ for potentially active (r) and $F^{(r)}(t+dt) \le \mathcal{F}^{(r)}(t+dt)$ for all (r) with $(\mathcal{N}^{(r)} l^{(r)}) \ge 0$.

This may be explicitly written in the form

$$F^{(r)} = b h_i^{(r)} n_j^{(r)} \sigma^{ij} = b \tau^{(r)} \approx \mathcal{F}_{PN}^{(r)} + \mu b^2 \sum_{\text{at all } (k)} D_{(k)}^{(r)} (\mathcal{N}^{(k)} l^{(k)})^{1/2} \frac{\overline{(\mathcal{N}^{(r)} l^{(r)})}}{2(\mathcal{N}^{(k)} l^{(k)})^{1/2}}$$

by using $(\mathcal{N}^{(k)} l^{(k)}) = 2\gamma^{(k)}/bL^{(k)}$ we obtain

$$\dot{\tau}^{(r)} \approx \sum_{\text{at active } (k)} H_{(k)}^{(r)} \gamma^{(k)} \quad (\text{with } \gamma_{(k)} \ge 0).$$

So, indeed, we find that

$$v_j^{ni} = \sum_{\text{at active } (k)} \frac{\partial \tau^{(k)}}{\partial \sigma^{ij}} \gamma^{(k)} = \sum_{\text{at active } (k)} h_i^{(k)} n_j^{(k)} \gamma^{(k)}$$

and that the increase of the limit for any system (r) is a linear function of $\gamma^{(k)}$ ($\gamma^{(k)} \ge 0$) of the active systems (k).

But the classical interaction matrix: $H_{(k)}^{(r)} \equiv \mu D_{(k)}^{(r)} / (\mathcal{N}^{(k)} l^{(k)})^{1/2} L^{(k)}$ as introduced by Taylor (1934), Koïter, Mandel (1965) in the multiple plastic potential theory, and used thereafter by others, is here found explicitly. It is clear that this matrix is not constant and is not symmetric. This is linked to the fact that large changes of the dislocation distribution are occurring during the deformation.

It is evident that, although $D_{(r)}^{(r)} \ge D_{(k)}^{(r)}$ for $r \ne k$, since $L^{(r)} \gg L^{(k)}$ for $r \ne k$ (the density of dislocations $\rho(r)$ which pierce the actual principal slip plane stays constant, while the other densities $\rho(k)$ have considerably increased, we obtain $H_{(r)}^{(r)} \le H_{(k)}^{(r)}$ (this is called by most of the researchers, except the authors of this work, the self- and latent hardening).

Then, there is the very important problem of how to determine the active systems (k) and to calculate $\gamma^{(k)}$ when we impose the objective stress rate, or $\dot{\tau}^{(r)}$, on each system (r). By using as the new intermediate variables

$$\dot{X}^{(r)} = \frac{\overline{(\mathcal{N}^{(r)} l^{(r)})}}{2(\mathcal{N}^{(k)} l^{(k)})^{1/2}}$$

it is easy to show that, since the matrix $D(D_{(k)}^{(r)})$ is symmetric and is strictly positive for any number of systems lower than 8, we have a quadratic minimization problem with the linear constraints $\dot{X}^{(r)} \ge 0$ which has one unique solution. This may then be solved easily with the help of many classical mathematical programming algorithms.

### 1.3.4. Simulation of tensile experiments

Generally, only experiments with tension or compression are performed on single crystals. These experiments are very difficult to realize and to interpret, chiefly owing to the large changes in the geometry.

During the test, usually, a constant ramp velocity is applied and only the applied force $Q$ and the crystal length $D$ are measured, but the tensile curves are often given in the plane $(\tau, \Gamma)$ where $\tau$ is the applied shear stress and $\Gamma$ is the glide on the principal slip system; this is a strict mathematical fiction, except when there is only one active slip mechanism.

We built two different programs: the first one was based on the direct explicit integration of the evolution viscoplastic equations; the Treanor algorithm (which is similar to the order 4 Runge-Kutta) with an adaptive
<figure><img src="image_5.png" /><figcaption>Fig. 4. Geometrical representation of a tensile test on the crystal.</figcaption></figure>

time increment, was found to be very efficient; the second one was based on the multiple plastic potential and then we used an implicit scheme of integration and the quadratic programming algorithms. We integrated for various initial conditions of the tensile axes and the initial states. We were able to interpret not only the two important stages of work-hardening which are stage I and stage II, but also the transition stages (since here, cross-slip and climb are neglected in this approximation, stage III is not reached). We were able to get the dependence on the initial orientation of the tensile axes without adding any phenomenological constant, and to understand what happens during the straining, explicitly. We compared the theoretical curves obtained based on the model to the experimental results; the tendency was very well represented (Figs. 4 and 5).

Although we still consider our model as the only real physical model which truly preserves the rational 3D hypothesis (and which differs from Havner (1972), Asaro (1983), ...), we think that its main interest is at the fundamental level and it can not be directly used for engineering applications: it concerns the single crystal; too many parameters have to be taken into account and moreover a better experimental fitting (less than a few percents) needs to be reached.

<figure><img src="image_6.png" /><figcaption>Fig. 5. Global tensile responses with the orientation of the crystal.</figcaption></figure>
## 2. Physical modeling of metallic polycrystal

### 2.1. Averaging methods

A polycrystal which is considered as a continuum point, contains several crystals at a microscopic scale. According to the stress gradient and the rate of loading, different representative volume elements may be introduced for which global behavior is defined based on the local behavior of the crystals and their relative interactions.

In Nemat-Nasser and Hori (1993), the reader may find a perfect comprehensive review of the state-of-the-art in the averaging homogenization techniques, bounds, theorems and simplified methods for the elastic or inelastic behavior of the RVE which may contain cracks. This state results from all the works by several hundred authors.

Under some regularity hypotheses, particularly for homogeneous uniform boundary tractions or linear boundary displacements, it was shown that the macroscopic or global stress, the global strain and the global stress-work for the RVE, are direct volume averages of the respective microscopic quantity:

$$ \Sigma = \frac{1}{V} \int_V \sigma \, dv = \langle \sigma \rangle = \bar{\sigma}, $$

$$ E = \frac{1}{V} \int_V \varepsilon \, dv = \langle \varepsilon \rangle = \bar{\varepsilon}, $$

$$ \Sigma E = \langle \sigma \varepsilon \rangle = \bar{\sigma} \bar{\varepsilon}. $$

But, even for small deformations, if we split the local strain $\varepsilon$ into its elastic $\varepsilon' \equiv M\sigma$ (or $\sigma = L\varepsilon$) and plastic $\varepsilon'' \equiv \varepsilon - \varepsilon'$ component, and we perform the decomposition for the global strain $E \equiv E' + E''$ with $E' \equiv \bar{M}\Sigma$ (or $\Sigma = \bar{L}\varepsilon$), more complex formulae are found by defining the concentration tensors with $\sigma \equiv A\Sigma$ or $\varepsilon = BE$, it may be shown that:

$$ E \equiv \frac{1}{V} \int_V \varepsilon A \, dv, \quad E' \equiv \frac{1}{V} \int_V \varepsilon' A \, dv, \quad E'' \equiv \frac{1}{V} \int_V \varepsilon'' A \, dv. $$

The same concentration tensors allow us to express the global elastic tensors as

$$ \bar{M} = \frac{1}{V} \int_V M A \, dv \quad \text{or} \quad \bar{L} = \frac{1}{V} \int_V L B \, dv. $$

One of the main problems for the global analysis of the polycrystal is the determination of the concentration tensors. Several simplified models (Sachs, Voigt, 1889; Taylor, 1938; Lin, 1971) and particularly the self-consistent model (Kröner, 1958; Budiansky and Wu, 1962) were proposed to solve it.

In the case of local elastic isotropy, when the shape of the crystals are approximated by a sphere, the self-consistent model gives a simple formula to compute the interaction between crystals:

$$ \sigma = \Sigma + 2\mu(1 - \beta)(E'' - \varepsilon'') $$

where $\beta$ is a constant = $2(4 - 5\nu)/15(1 - \nu)$, $\mu$ is the shear modulus and $\nu$ is the Poisson ratio.

The model was extended to anisotropic elasticity by various authors which drives to a modified interaction formula (Bui, 1969)

$$ \sigma = A[\Sigma + 2\mu(1 - \beta)(E'' - \varepsilon'')] $$
or to take into account the elastoviscoplastic stress relaxation:

$$ \sigma = \Sigma + 2\mu(1 - \beta)(E'' - \varepsilon'')\lambda(\Sigma, E'', P) $$

with explicit A and $\lambda(\Sigma, E'', P)$ (Zaoui, 1970).

## 2.2. Physically based model

Since many researchers are still working in this area, we only give here the results that we obtained for the global behavior of FCC polycrystals for small deformation, based particularly on some unpublished work with Compere (1973) and Engel (1975). In fact, we underline the use of the single crystal model, the relation in a RVE with the self-consistent model, simplified methods, or the more sophisticated finite element method; we underline also the use of the direct explicit time integration of the evolution viscoplastic equations (Treanor' algorithm) or the implicit scheme of time integration of the plastic equations based on quadratic programming algorithms.

### 2.2.1. Uniaxial tension

We consider a polycrystal RVE which has an infinite size and with an infinite number of crystals with initial random orientation characterized by 3 orthogonal vectors $\vec{X}$, $\vec{Y}$, $\vec{Z}$. In a fixed frame of reference (following the notations of Zaoui (1970)), $\vec{x}$, $\vec{y}$, $\vec{z}$, with $\vec{u}$ being the axis along the intersection between the two planes $(x, y)$ and $(X, Y)$, we define the angles $\Psi = (\vec{x}, \vec{u})$ rotation with $\vec{z}$ axis $\theta = (\vec{z}, \vec{Z})$ rotation with $\vec{u}$ axis and $\varphi = (\vec{u}, \vec{X})$ rotation with $\vec{Z}$ axis.

To cover all the orientations, it is just sufficient to limit the variations of the 3 angles $(\Psi, \theta, \varphi)$ to:

$$ \sin \varphi \sin \theta \geq 0, \quad \cos \theta \leq \sin \varphi \sin \theta, \quad \sin \theta \cos \varphi \geq 0, \quad \cos \theta \geq 0, \quad \cos \theta \leq \cos \varphi \sin \theta $$

which define in the space $(\varphi, \cos \theta, \Psi)$ a cylinder parallel to $\Psi$, bounded by the 2 planes $\Psi = 0$ and $\Psi = 2\pi$, as indicated in Fig. 6. In a tensile test, the angle $\Psi$ is not taking part.

We have then addressed the following questions:

(1) How many different orientations of crystals are necessary to reach the description of an initially isotropic polycrystal? We have used our fcc single crystal model for the local behavior. We have also used the simple

<figure><img src="image_7.png" /><figcaption>Fig. 6. Limits on the angles for the isotropic polycrystal.</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>Fig. 7. Differences with the number of crystals for the global response.</figcaption></figure>

self-consistent model. We have selected 3, 24, 36, 54, 117 orientations uniformly distributed in the plane $(\varphi, \cos \theta)$.

We did not notice great differences between them (Fig. 7). Then, we concluded that 3 well selected orientations were sufficient for the simulation.

(2) Are the simplified models, such as the self-consistent model, sufficient to represent the interactions between crystals? We have taken the previous model of the polycrystal but in the same time, we performed the finite element numerical analysis. At that time, it was necessary to develop a program for rate-dependent or rate-independent crystal plasticity (the program was also written for large transformations).

We took a polycrystalline specimen containing 216 tetrahedral elements crystals with distributed orientation as in Fig. 8.

<figure><img src="image_9.png" /><figcaption>Fig. 8. Finite element meshes for the numerical simulation.</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>Fig. 9. Differences between local and global simulations.</figcaption></figure>

Notation for the orientations of the crystals

<table><thead><tr><th></th><th>A</th><th>B</th><th>C</th><th>D</th><th>E</th><th>F</th></tr></thead><tbody><tr><th>φ</th><td>22.5°</td><td>37.5°</td><td>37.5°</td><td>52.5°</td><td>52.5°</td><td>67.5°</td></tr><tr><th>θ</th><td>82.5°</td><td>82.5°</td><td>67.5°</td><td>82.5°</td><td>67.5°</td><td>82.5°</td></tr><tr><th>ψ</th><td>0°</td><td>60°</td><td>120°</td><td>180°</td><td>240°</td><td>300°</td></tr></tbody></table>

Repartition per layer of 9 cubes (A1 meaning φ = 22.5°, θ = 82.5° and ψ = 0)

<table><tr><td>A1</td><td>E2</td><td>B3</td><td>C6</td><td>D1</td><td>A2</td><td>D3</td><td>C2</td><td>F1</td><td>F2</td><td>B1</td><td>E6</td></tr><tr><td>D6</td><td>C5</td><td>F4 → (1 to 9),</td><td>B5</td><td>F6</td><td>C1 → (10 to 18),</td><td>A4</td><td>E5</td><td>B6 → (13 to 27),</td><td>E1</td><td>A6</td><td>D5 → (28 to 36).</td></tr><tr><td>F5</td><td>B4</td><td>E3</td><td>E4</td><td>A3</td><td>D2</td><td>C3</td><td>D4</td><td>A5</td><td>B2</td><td>F3</td><td>C4</td></tr></table>

We did not find great differences between the two kinds of analysis and we concluded that the simple self-consistent model, without any other improvement, was sufficient for the analysis (Fig. 9).

(3) Is it necessary to modify the values of the coefficients of the interaction matrix $D_{(k)}^{(r)}$ between slip systems which were explicitly computed with the physical model of the crystals?

<figure><img src="image_11.png" /><figcaption>Fig. 10. Influence of the absolute value of the interaction matrix.</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>Fig. 11. Influence of the relative value of the latent hardening.</figcaption></figure>

In Figs. 10 and 11, we have performed a number of different simulations using various assumed values for the interaction coefficients (as from a phenomenological point of view, these coefficients may take different values) and even that artificially, the latent hardening could be more important than the self-hardening.

This may allow a greater flexibility to reach a better agreement with the experimental results (we were able to represent rather well the tensile behavior of copper).

In Fig. 12, we have also simulated the case where the density of dislocations for any slip system reaches a limiting value following the idea that the mean distance D between 2 dislocations can not be smaller than a critical value.

### 2.2.2. Multiaxial loadings

We have also analyzed the global behavior of the polycrystal during biaxial loading along the z and y axis, particularly, when there is a viscoplastic potential $\Omega$ when the crystals admit viscoplastic potentials $\mathcal{F}$ as shown by Rice (1970),

$$ \Omega = \frac{1}{V} \int_V \mathcal{F} \, dv. $$

We have simulated the evolution of the elastic yield surface, $\Omega = 0$ and the potential surface, $\Omega = 10^{-9}$ of

<figure><img src="image_13.png" /><figcaption>Fig. 12. Phenomenological introduction of a saturation state.</figcaption></figure>
<figure><img src="image_14.png" /><figcaption>Fig. 13. Elastic yield surface with the axial strain.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>Evolution of the equipotential with the strain (isotropic elasticity)</figcaption></figure>

Fig. 14. Viscoplastic potential surface with the axial strain.

<figure><img src="image_16.png" /><figcaption>Evolution of the equipotential with the strain (anisotropic elasticity)</figcaption></figure>

Fig. 15. Viscoplastic potential surface with the axial strain.
the polycrystal, during a loading along the z-axis. When the initial yield surface looks like the Tresca criterium, it quickly changes with its size is considerably reduced. On the other hand, when the viscoplastic surface initially looks like the Mises criterium, its subsequent evolution consists of a rotation, a translation and a distortion, which are functions of the amplitude of the actual axial strain Figs. 13 and 14.

Finally, in Fig. 15, we also examined the effect of the elastic anisotropy of the crystals and their plastic rotation, using the modified self-consistent formula.

### 2.2.3. Final remark

We have shown that, to represent the global physical behavior of the polycrystal, it is sufficient to use only 3 orientations of the crystals for both rate-dependent or rate-independent crystal plasticity and the simple self-consistent model. Although we have obtained a reasonable fit to the experimental data, we do not think that the method can be used in industrial applications; the calculations require at least 72 parameters!

## 3. Quasi-physical modeling of aggregates

A more useful method was developed by retaining the main features of the microscopic approach, but using simple and elementary analysis (inter and intra crystalline couplings and the construction of explicit relations). This led to a new approach of analysis of modelling materials which ensures a much better fit to the tests and chiefly drives to elementary analysis of structures where only elastic computations give practical bounds of its response to many complex loadings. We refer the reader to this book for its complete description. Here, we shall just underline the main points.

### 3.1. Hypothesis

Our aim is to consider the simplest model which is still able to represent the experimental results instead of the exact physical model.

Based on some works by Mandel (1971) and Nguyen (1973) or Halphen and Nguyen (1975), we assume that, at a local scale, there are some inelastic elementary mechanisms which are the sources of producing

* instantaneous plastic strains with a threshold (as in the glide of a dislocation);

* viscous strains without a threshold (as in the thermally activated diffusion of a point defect or precipitate);

* visco-plastic strains with a threshold (as in the thermally activated crossing between two dislocations).

### 3.2. Basic equations

The basic equations associated with this quasi-physical modeling involve

* some internal parameters say,

$$x^T = \begin{bmatrix} \alpha \\ \beta \\ \gamma \end{bmatrix}$$

which are linked to the local inelastic strains of each of the mechanisms, local sources of inelasticity.

* The local stresses associated with on each of these mechanisms are expressed as

$$\begin{bmatrix} \sigma_{\alpha} \\ \sigma_{\beta} \\ \sigma_{\gamma} \end{bmatrix} = \begin{bmatrix} A_{\alpha} \\ A_{\beta} \\ A_{\gamma} \end{bmatrix} \Sigma - \begin{bmatrix} \hat{\alpha} \\ \hat{\beta} \\ \hat{\gamma} \end{bmatrix} = A \Sigma - y,$$
where $\Sigma$ is the global applied stress tensor, $\sigma_\alpha$, $\sigma_\beta$, $\sigma_\gamma$, are the respective local stresses, $A_\alpha$, $A_\beta$, $A_\gamma$ are the elastic stress localization tensors and $y$ stand for the transformed internal parameters.

• These transformed internal parameters (similar to eigenstresses) are linked to the internal parameters by

$$y = \begin{bmatrix} \hat{\alpha} \\ \hat{\beta} \\ \hat{\gamma} \end{bmatrix} = B \begin{bmatrix} \alpha \\ \beta \\ \gamma \end{bmatrix} = Bx$$

where $B$ is a non-negative symmetrical matrix with only the following two possible properties:

(i) $B$ is strictly positive; there is then a one-to-one relation between the transformed internal parameters and the internal parameters;

(ii) $B$ is singular; the transformed internal parameters $y$ must belong to a subspace of compatibility.

• The evolution laws for each of these inelastic mechanisms are written as

(i) $\dot{\alpha} \in \partial\Psi_C(\sigma_\alpha)$ for the instantaneous irreversible strains with a threshold;

(ii) $\dot{\beta} = f(\sigma_\beta)$ for the viscous strains where $f$ may be taken as a quadratic function;

(iii) $\dot{\gamma} = \partial\Omega(\sigma\gamma)/\partial\sigma\gamma$ for the delayed irreversible strains with a threshold and $\Omega$ is the viscoplastic potential.

• The global inelastic strain rate is taken as

$$\dot{E}'' = A_{\alpha}^{T} \dot{\alpha} + A_{\beta}^{T} \dot{\beta} + A_{\gamma}^{T} \dot{\gamma} = A^{T} \dot{x}$$

and the global elastic strain rate as

$$\dot{E}' = \overline{M} \dot{\Sigma}$$

### 3.3. Remark

We have shown that, using the above outlined modelling approach, it is possible to represent all the observed experimental results during monotonic or cyclic multiaxial loadings with a very limited number of internal parameters. Moreover, we have shown that, when using these relationships in numerical simulations of complex structures, it is easy to obtain evaluation of the response with only a very limited number of linear elastic analysis which may reduce the execution time of the computer by a factor of 20 to 1000.

However, in our view, this was not yet sufficient for real practical industrial applications.

## 4. Intelligent modeling of elastic aggregates

### 4.1. Principles

Recently, we have developed an approach based on automatic learning and optimization techniques, which can use the entire body of the existing knowledge of the experts, the experimental results, computational simulations, and produce practical results for engineering applications.

For this approach, it is necessary:

(i) To build a data base, i.e., to obtain some experimental, real or simulated results, where the experts of each particular problem identify all variables or descriptors which they think should be essential for the solution of the considered problem. This is, at first, done with some primitive descriptors $x$ which are usually limited in number. Then, the experts transform the data into a (generally much larger) number of intelligent descriptors XX. These descriptors represent the actual existing knowledge and all the relevant (but often insufficient)
theories. Each descriptor may be a number, boolean, alphanumeric, name of files to provide access to data bases, or curve, signal, or image. The results or conclusions may be classes (good, not good, ...), or numbers (Young modulus, cost, weight, life time, ...). Usually, the data base may contain 20 to 50 examples with 10 to 1000 descriptors and 1 to 20 conclusions.

(ii) To generate the rules with any automatic learning tool. The intelligent descriptors help these learning algorithms. Each conclusion is explained as a function or a set of rules of some among the input intelligent descriptors with known reliability or accuracy. If this reliability is too low, it implies that there is insufficient data or there are bad, missing descriptors, or that the problem was not well described.

(iii) To optimize at two levels (inverse problems):

• Considering the intelligent descriptors as independent; it is possible to obtain the optimal solution satisfying some specially required properties and allowing the discovery of new mechanisms.

• Considering the intelligent descriptors linked to primitive descriptors; it is possible to obtain the optimal solution which is technologically feasible.

So, not only a practical optimal solution is obtained (since good interpolations and extrapolations for new cases are made) but also the experts of the problem may learn the missing parts, may build models or theories based only on the retained intelligent descriptors and guided by the structure of the rules or relationships.

## 4.2. Application to an elastic aggregate

We refer particularly to the chapter in the book by Nemat-Nasser and Hori (1993) on elastic solids with inclusions. Several models for approximated values or bounds of the global elastic tensors were proposed. Almost all these classical models are based on only the volume fraction $v_i = V_i/V$ of the inclusions.

Here we consider a special RVE which contains randomly distributed aligned fibers with the same rectangular section (plane strain problem). The fibers are along the z-direction. Since this RVE has a finite size $V$, we have a special distribution which implies a general orthotropic material instead of a transversely isotropic

<figure><img src="image_17.png" /><figcaption>Fig. 16. Inclusions/fibers in the representative volume element RVE: finite element mesh.</figcaption></figure>
one; see Fig. 16. We want to obtain the values of the effective moduli of a particular finite sized RVE, using random distributions of the fibers. For each case, a special distribution needs to be defined. Some mathematical morphological tools are available for this purpose.

Here, let us assume that this distribution may just be defined by the following parameters: $v_i$, the volume fraction; $n_i$, the number per unit volume or cross-section; $l_i/h_i$, the relative aspect ratio of the cross-section of fibers; $C_{ix}$, $C_{iy}$, the relative coordinates of the center of inertia $C_i$ of the fibers; $R_{max}$, the radius of the circle with center at $C_i$ which contains all the cross-sections of the fibers; $d_{min}$, the minimum distance between two centers of cross-section; dis1, the mean distance between two centers of cross-section; and dis2, the mean distance between a center and the center of inertia $C_i$.

### 4.2.1. Building the data base

Let us assume that we have a data base of examples where the distribution of the fibers, and the local and global elastic properties are given. These global elastic coefficients may have been measured during tests or, as a less interesting case, they may have been computed numerically using for example a finite element program.

Since our illustrative example has only a pedagogical interest here, we have created the data base numerically in plane strain. We applied an homogeneous boundary displacement field:

$$ \vec{U} = \begin{cases} u = ax + by \\ v = cy \\ w = 0 \end{cases} $$

which implies the global strain tensor

$$ \Rightarrow \mathbf{E} = \begin{pmatrix} E_{xx} = a & E_{xy} = b/2 & 0 \\ S & E_{yy} = c & 0 \\ S & S & 0 \end{pmatrix} $$

or equivalently the global strain vector

$$ \Rightarrow \vec{E} = \begin{cases} E_1 = E_{xx} = a \\ E_2 = E_{yy} = c \\ E_3 = 0 \\ E_4 = 0 \\ E_5 = 0 \\ E_6 = 2 E_{xy} = b \end{cases} $$

With any linear finite element program (Navidi and Azhdari used respectively CADSAP from Algor and NIKE-2D from Laurence Livermore Laboratory), we compute the local stress field $\sigma$ and by directly averaging, we define the global stress tensor

$$ \Sigma = \frac{1}{V} \int_V \sigma \, dv $$
or equivalently its vector form

$$ \vec{\Sigma} = \begin{pmatrix} \Sigma_1 = \Sigma_{xx} \\ \Sigma_2 = \Sigma_{yy} \\ \Sigma_3 = \Sigma_{zz} \\ \Sigma_4 = \Sigma_{yz} = 0 \\ \Sigma_5 = \Sigma_{xz} = 0 \\ \Sigma_6 = \Sigma_{xy} \end{pmatrix}. $$

Then the effective global elastic coefficients are defined by

$$ \bar{L} = \begin{pmatrix} \bar{L}_{11} & \bar{L}_{12} & \bar{L}_{13} & 0 & 0 & 0 \\ S & \bar{L}_{22} & \bar{L}_{23} & 0 & 0 & 0 \\ S & S & \bar{L}_{33} & 0 & 0 & 0 \\ S & S & S & \bar{L}_{44} & 0 & 0 \\ S & S & S & S & \bar{L}_{55} & 0 \\ S & S & S & S & S & \bar{L}_{66} \end{pmatrix}. $$

Indeed due to our hypothesis of plane strain, only 6 of these coefficients are identified: $\bar{L}_{11}$, $\bar{L}_{12}$, $\bar{L}_{22}$, $\bar{L}_{66}$, $\bar{L}_{13}$, $\bar{L}_{23}$.

The description of each example contains:

(i) The classical primitive descriptors $E_m$, $v_m$ (elastic constants of matrix), $E_i$, $v_i$ (elastic constants of inclusion), $v_i$ (concentration of inclusions or volume fraction) from which we compute for the matrix: Lamé coefficients $\lambda_m = E_m v_m / (1 + v_m)(1 - 2v_m)$, $\mu_m = E_m / 2(1 + v_m)$ and $K_m = E_m / 3(1 - 2v_m)$ and the same coefficients for the inclusion: $\lambda_i$, $\mu_i$ and $K_i$.

(ii) The intelligent descriptors of the inclusions distribution within the particular volume as described earlier: $n_i$, $l_i/h$, $C_{ix}$, $C_{iy}$ ($C_i$ coordinates), $R_{max}$, $d_{min}$, dis1, dis2.

(iii) The intelligent descriptors which represent, in principle, the best actual knowledge of the experts (for simplification reasons, here, we only take the predicted values given by Behrens (1967) but one of the Hashin–Shtrikman bounds or exact Nemat–Nasser’s bounds or improved Willis bounds, could have been taken, keeping in mind, however, that the simplest one is usually the best one; we do not want expensive descriptors which may need cumbersome computations or complex testings).

The values are only given for the transversely isotropic medium and only for the concentration of inclusions or volume fraction, with

$$ A = \mu_i + \lambda_i + \mu_m, \quad B = \mu_i + \lambda_i - \mu_m - \lambda_m, \quad C = \mu_i - \mu_m, \quad q = \frac{2\mu_m + \lambda_m}{\mu_m}, \quad A' = \mu_i + \lambda_i + \lambda_m $$

$$ \bar{L}_{11}^P = \bar{L}_{22}^P = (2\mu_m + \lambda_m) \left[ \frac{A}{A - Bv_i} + \frac{Cv_i}{A' - C(v_i)^q} \right] $$

$$ \bar{L}_{13}^P = \bar{L}_{23}^P = \lambda_m + \frac{(\lambda_i - \lambda_m)(2\mu_m + \lambda_m)v_i}{A - Bv_i} $$
$$ \bar{L}_{33}^{P} = 2\mu_{m} + \lambda_{m} + (2\mu_{i} + \lambda_{i} - 2\mu_{m} - \lambda_{m})v_{i} - \frac{(\lambda_{i} - \lambda_{m})^{2}v_{i}(1 - v_{i})}{A - B} $$

$$ \bar{L}_{66}^{P} = \mu_{m} \left[ 1 + \frac{Cv_{i}}{\mu_{i} - C(v_{i})^{1/q}} \right] $$

$$ \bar{L}_{12}^{P} = \lambda_{m} + (2\mu_{m} + \lambda_{m})v_{i} \left[ \frac{B}{A - Bv_{i}} - \frac{C}{A' - C(v_{i})^{q}} \right] $$

The data base has then the following structure:

* 18 input descriptors; filename (only as reference), y_inc, p_inc, ymat, pmat, volfrac, nb_incl, aspect, C,y, C.z, R, dis1, dis2, dmin, and the 4 Berhens descriptors L11P, L31P, L12P, L66P.

* 6 output descriptors or conclusions; L11, L22, L12, L66, L13, L23.

### 4.2.2. Learning (Fig. 17)

Using the numerical option of LES from LMS of $X^2$ we can learn each of the six conclusions.

This automatic learning when using all the input descriptors and special learning control parameters (then taking the Berhens intelligent descriptors) may give for the first conclusion L11:@# L11 = 4.715796561017032e-01*y.mat + 4.557336454329526e-01*L11P - 7.575868226106983e-03*y_inc + 4.760595153407522e-01*L31P + 6.069191287762840e-01*y_inc*volfrac - 4.600487272203139e-01*volfrac*L11P + 1.262632291054333e-06*p_inc*L12P**2 - 1.879706062534128e-06*y.mat**2*volfrac + 4.151452553084349e-07*C.y*L12P**2 - 1.219946907404129e-17*L31P**2*L12P**2 with the errors: standard-deviation: $1.33 \times 10^6$, skewness: 0.895992, kurtosis: 2.78957.

The direct automatic learning (without the introduction of the Berhens descriptors) gives for the same conclusion: L11 = 3.88E-01*y.mat - 9.21E + 03*aspect + 3.12E - 06*y_inc*y.mat - 7.82E - 01*y_inc*p.mat + 1.76E + 00*y_inc*volfrac + 8.46E - 01*p_inc*y.mat + 2.14E + 00*y.mat*p.mat - 1.41E + 00*y.mat*volfrac + 2.24E + 05*p.mat**2 - 8.45E + 04*volfrac**2, but the errors are much higher: standard-deviation: $9.85 \times 10^6$, skewness: 1.46374, kurtosis: 5.89781.

We can also use other automatic learning techniques particularly those which are based on neural networks such as NEUROSHELL<fnref n="3" />.

So, we can predict for any mixture and any distribution of inclusions the effective properties with a reasonable accuracy, i.e., we can conclude.

However, what is more important in our personal view, is that we may try to find the special mixture and distribution which will satisfy a priori given requirements using for example, the functionality OPTIMIZE of LES of X or the program GENEHUNTER of Ward Systems, which are both based on genetic algorithms. Here, we may find the composite material with the lowest weight (the objective function), which will have however some required bounds for its elastic moduli knowing that the materials for the matrix and the inclusions belong

<footnote n="2">² The LES system was developed at the Laboratoire de Mécanique des Solides de l'Ecole Polytechnique. The software is given to the participants of the intensive workshop Intelligent Optimal Design of Materials and Structures organized regularly at Ecole Polytechnique (France) and at University of California San Diego (USA). It is also distributed by CADLM, 9 Rue Raoul Dautry, 91190 Gif/Yvette, France. Tel.: +33-1-69072922; fax: +1-33-1-69072809; e-mail: cadlm@wanadoo.fr. The system has a symbolic option (by Dr. M. Sebag and M. Schonauer) and a numeric option (by Mr. M. Terrien).</footnote>

<footnote n="3">³ This system is developed and distributed by WARDS Systems Group, Inc., Executive Park West, 5 Hillcrest Drive, Frederick, MD 21703, USA. Tel.: +1-301-6627950; fax: +1-301-6625666; e-mail: wardsystems@msn.com.</footnote>
<figure><img src="image_18.png" /><figcaption>Fig. 17. Learning of the 6 conclusions and coupling with optimization within EXCEL.</figcaption></figure>

to some classes with bounded elastic properties. But any new objective may be given and the results are obtained instantaneously (Fig. 17).

## Acknowledgements

Most of this paper was written during a visit to UCSD's Center of Excellence for Advanced Materials, under ONR contract N00014-96-1-0631 to the University of California San Diego. The authors want to thank Dr. M. Beizaie for his efficient assistance during the writing of this paper and Dr. A. Azhdari for the numerical simulations of the elastic aggregates to obtain the effective elastic coefficients. They want to thank Mr. M. Terrien for his advises in the use of his Learning Expert System generator, and at last, Professor J. Willis for his review of the paper and Professor S. Nemat-Nasser for his friendly discussions who initiated the paper and for his intensive review of the text.

## References

* Asaro, R.J., 1983. Micromechanics of crystals and polycrystals. Adv. Solid Mech. 23, 1-115.

* Budiansky, B., Wu, T.K., 1962. Theoretical prediction of plastic strains of polycrystals. Proc. 4th US Natl. Cong. Appl. Mech., 1175-1185.

* Behrens, E., 1967. Elastic constants of filamentary composites with rectangular symmetry. J. Acoustical Soc. America 42, 367-386.

* Bui, H.D., 1969. Thesis, Paris.

* Compere, P., 1973. These de Docteur-Ingenieur, Universite Paris XI.

* Engel, J.J., 1975. Thesis, Paris.
Halphen, B., Nguyen, Q.S., 1975. Sur les materiaux standards generalises, J. Mec. 14-28.

Havner, K.S., 1972. An analytical model of large deformation effects in crystalline aggregates. In: Sawczuk, A. (Ed.), International Symposium on Foundations of Plasticity, Varsovie. Noordhof-International Publishing.

Lin, T.H., 1971. Physical theory of plasticity, in: Advances in Applied Mechanics, vol. 2. Academic Press.

Mandel, J., 1965. Generalisation de la theorie de plasticite de W.T. Kolter. Int. J. Sol. Struct. 1, 273-296.

Mandel, J., 1971. Courses and Lectures, in: Plasticite Classique et Viscoplasticite, CIDM, Udine, Springer, New York.

Nemat-Nasser, S., Hori, M., 1993. Micromechanics, Elsevier.

Nowacki, W.K., Zarka, J., 1971. Etude de la frontiere elastique des monocristaux d'aluminium. Int. J. Solids Structure 7, 1277-1287.

Nguyen, Q.S., 1973. Thesis, Paris.

Rice, J.R., 1970. On the structure of stress-strain relations for time-dependent plastic deformation in metals. J. Appl. Mech. 37, 728-737.

Taylor, G.I., 1934. The mechanism of plastic deformation of crystals; Part I theoretical. Proc. Roy. Soc. London A 145, 362-387.

Zaoui, A., 1970. Thesis, Paris.

Zarka, J., 1968. Thesis, Paris. In: Sur la Viscoplasticite des Metaux. Mem. Art. Fr., 2, 223-292.