# Application of the homogenization method to composite materials with periodic structure

Jean-François Bourgat, Hélène Lanchon

► To cite this version:

Jean-François Bourgat, Hélène Lanchon. Application of the homogenization method to composite materials with periodic structure. [Research Report] IRIA-RR-208, IRIA. 1976, pp.35. hal-04716344

HAL Id: hal-04716344

https://inria.hal.science/hal-04716344v1

Submitted on 1 Oct 2024

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 Attribution 4.0 International License
Inv 202535 RR 208

IRIA

laboria

Institut de Recherche
d'Informatique
et d'Automatique

Domaine de Voluceau
Rocquencourt
B. P. 105 78150 - Le Chesnay
France
Tél.: 954 90 20

laboratoire de recherche
en informatique
et automatique

ARCHIVES

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

APPLICATION OF THE
HOMOGENIZATION METHOD
TO COMPOSITE MATERIALS
WITH PERIODIC STRUCTURE

115

Jean François BOURGAT
Hélène LANCHON

Rapport de Recherche N° 208
Décembre 1976

35p
APPLICATION OF THE HOMOGENIZATION METHOD TO
COMPOSITE MATERIALS WITH PERIODIC STRUCTURE

J.F. Bourgat*, H. Lanchon**

ABSTRACT

A boundary value problem in an heterogenous medium is very difficult, and even impossible, to solve if the number of heterogeneities is very large. The homogenization consists to replace the heterogeneous medium by an "equivalent" homogeneous one. In the particular case of a periodical distribution of heterogeneities, results on the mathematical aspects of the homogenization have been obtained recently. An application to the problem of torsion for a cylindrical fiber-reinforced bar is presented here. The homogenized solutions are compared with the real solutions for an increasing number of fibers.

RÉSUMÉ

Dans les milieux physiques hétérogènes les problèmes aux limites sont difficiles à résoudre numé-
riquement, parfois impossibles, si le nombre des hétérogénéités est très grand. L'homogénéisation
consiste à remplacer le milieu hétérogène par un milieu homogène "équivalent". Dans le cas parti-
culier de distributions périodiques d'hétérogénéités, des résultats sur les aspects mathématiques de
l'homogénéisation ont été obtenus récemment. Nous présentons ici une application au problème
de la torsion d'une barre cylindrique renforcée par des fibres. Les solutions obtenues par la mé-
thode d'homogénéisation sont comparées aux solutions réelles pour des nombres de fibres croissants.

<footnote>* IRIA-LABORIA, Laboratoire de Recherche en Informatique et Automatique, B.P. 105 - 78150 Le Chesnay (France)</footnote>

<footnote>** LEMTA, Laboratoire d'Energétique et de Mécanique Théorique et Appliquée, B.P. 850 - 54011 Nancy (France)</footnote>
# INTRODUCTION.

We are here concerned with a method of homogenization. The origin of this word is related to the question of replacement of the heterogeneous medium by an "equivalent" homogeneous one, either to obtain global properties or to determine some behaviours of real medium impossible to reach directly.

If we have to solve a boundary value problem, homogenization means replacement of the corresponding differential operator with variable coefficients by an "equivalent" one with constant coefficients. The constant coefficients can then be considered as a new kind of "effective modulus" because they are corresponding to the same formulation of the same boundary value problem.

In the particular case of a periodical distribution of heterogeneities, recent results on the mathematical aspects of the homogenization have been given by I. BABUSKA [1], A. BENSOUSSAN-J.L. LIONS-G. PAPANICOLAOU [2],[3],[4],[5], E. GIORGI-S. SPAGNOLO [6], and L. TARTAR [13].

The method has been applied by G. DUVAUT [7], [8] to composite plates and also to three dimensional elasticity.

Using some asymptotic approach, with the help of physical interpretation, E. SANCHEZ PALENCIA [12] has obtained the same homogenized coefficients.

In this paper we recall the mathematical model and the construction of the homogenized coefficients, we apply the method to the torsion of a cylindrical fiber-reinforced bar and we present numerical results comparing the homogenized solutions to the real solutions. Some elements about the mathematical and the numerical techniques are added.

# 1. - DEFINITION OF A PERIODIC HETEROGENEOUS MEDIUM.

Let consider an heterogeneous medium made of a regular distribution of several components. We can idealize this medium, admitting a kind of space periodicity, in this way : if $C(x)$ is a variable coefficient characterizing the (mechanical, thermical, electromagnetical...) comportment of the material, and if $\Omega$ is the bounded domain of $\mathbb{R}^N$ ($N = 1,2,...,3$) occupied by it, the space periodicity is defined by :

$$ (1) \qquad C(x+k\epsilon) = C(x) $$
∀k, N-integer and ∀x ∈ Ω such that x+kε ∈ Ω,

with

ε = (ε₁, ..., εₙ), εᵢ being the period in the xᵢ direction,

kε = (k₁ε₁, ..., kₙεₙ).

Then we can, for instance, introduce

$$ (2) \qquad \epsilon = \max \epsilon_i, \quad y_i = \frac{x_i}{\epsilon} \quad \text{and} \quad p_i = \frac{\epsilon_i}{\epsilon}, \quad 1 \le i \le N, $$

in such a way that :

$$ C(x) = C(\epsilon y) = a(y) $$

and

$$ a(y+kp) = a(y). $$

a, defined on the rectangular parallelepiped

$$ P = [0, p_1] \times \dots \times [0, p_N] \subset \mathbb{R}^N $$

is said P-periodic. P is the image of the representative cell of the medium.

*Example*: for a composite material with 2 components we shall have

$$ a(y) = \begin{cases} C_r & \text{on } p_r \text{ (image of the reinforced part of the cell)} \\ C_m & \text{on } p_m \text{ (image of the matrix part of the cell)} \end{cases} $$

The constants C_r and C_m being respectively the values of the considered coefficient for the reinforcement and the matrix. (See Fig. 1).

It is important to point out that the basic cell of the medium is not necessarily simple. In the case of an hexagonal array, the representative cell has the shape shown in Fig. 2.
<figure><img src="image_3.png" /><figcaption>Fig. 1 : A cell structure.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>Fig. 2 : Cell structure with hexagonal reinforcement.</figcaption></figure>
## 2. BOUNDARY VALUE PROBLEM.

Suppose now that we formulate a boundary value problem for a periodic heterogeneous medium. We have a partial differential equation with variable coefficients and, some boundary (eventually initial) conditions. Even if we can prove "existence and uniqueness of the solution, we are generally unable to obtain it by numerical method as soon as $ε$ is "small" because the rapid variations of the coefficients. We are then looking for an approached solution and that is the topic of the present theory.

To be clearer, let us choose, as example, a simple mathematical model obtained by the formulation of several physical situations. After having given a quick description of this model, we shall explicit a mechanical corresponding problem.

We shall use the convention of summation on the repeated index.

### 2.1. Mathematical model :

Let us consider the differential operator

$$ (3) \qquad A^{\epsilon} = - \frac{\partial}{\partial x_i} \left[ a_{ij} \left( \frac{x_j}{\epsilon} \right) \frac{\partial}{\partial x_j} \right] $$

where $a_{ij}$ are P-periodic bounded coefficients such that

$$ (4) \qquad \exists \beta > 0 ; a_{ij}(y) \xi_i \xi_j \ge \beta |\xi|^2, \forall \xi \in \mathbb{R}^N \text{ and } \forall y \in P, $$

and are not necessarily symmetrical in i and j.

Suppose now that we have to solve the following Dirichlet problem :

Problem ($P_ε$) : Find $u_ε$ such that :

$$ (5) \qquad \left\{ \begin{array}{ll} A^{\epsilon} u_{\epsilon} = f & \text{in } \Omega \\ u_{\epsilon} = 0 & \text{on } \partial\Omega \end{array} \right. $$

where $\Omega$ is an open bounded domain in $\mathbb{R}^N$ and $f$ is a given function in $L^2(\Omega)$.

We know by the general theory of partial differential equations that, thanks to the property of coercivity, there exists $u_ε$, unique solution of (5).
If $\epsilon$ is small by comparison with the size of $\Omega$, it is impossible to compute $u_{\epsilon}$. So it is natural to wonder if $u_{\epsilon}$ tends toward a limit $u$ when $\epsilon$ tends toward 0 and then if $u$ can be considered as a good approximation of $u_{\epsilon}$ as soon as $\epsilon$ is sufficiently small.

## 2.2. Torsion of fiber-reinforced bars :

If we formulate the standard problem of torsion for a cylindrical fiber-reinforced bar, it is known that :

<figure><img src="image_5.png" /><figcaption>Fig. 3
Fiber-reinforced bar.</figcaption></figure>

The stress tensor is such that

$$ (6) \quad \left\{ \begin{array}{l} \sigma_{ij} = 0 \text{ except} \\ \sigma_{13} = \frac{\partial \theta}{\partial x_2} = \sigma_{31} \text{ and} \\ \sigma_{23} = - \frac{\partial \theta}{\partial x_1} = \sigma_{32} \end{array} \right. $$

$\theta$, called "stress function" being defined on the cross section $\Omega$ of the bar and verifying the following boundary value problem :

$$ (7) \quad \left\{ \begin{array}{l} - \frac{\partial}{\partial x_i} \left[ \frac{1}{\mu(x)} \frac{\partial \theta}{\partial x_i} \right] = 2\alpha \quad \text{on } \Omega \\ \theta = 0 \quad \quad \quad \quad \quad \quad \quad \quad \quad \text{on } \partial\Omega \end{array} \right. $$

Where : $\mu$ is the shear modulus, ($\mu_f$ on the fibers, $\mu_m$ on the matrix) $\alpha$ the torsion angle imposed at the terminal section.

We shall suppose here that the distribution of fibers is in such a way that the period $\epsilon$ is the same in the two directions $\alpha x_1, \alpha x_2$.

Then, with the same change of coordinates as above :

$$ y_i = \frac{x_i}{\epsilon} $$

We have :

$$ \mu(x) = a(y) = \begin{cases} \mu_f & \text{on } P_f \\ \mu_m & \text{on } P_m \end{cases} $$
This problem is actually a ($P_\varepsilon$) one with :

$$a_{ij} = \frac{\delta_{ij}}{a(y)} \quad (\delta_{ij}, \text{ kronecker's index}),$$

$$f = 2\alpha, \text{ constant,}$$

$$\beta = \frac{1}{\sup_{y \in \Omega} a(y)} = \frac{1}{\mu_f} \quad (\text{if, as it is realistic } \mu_f > \mu_m)$$

**Remark 1** : What does mean here the fact to make $\varepsilon$ tend toward 0 ?

It implies that we consider a sequence of cylinders of same shape, same components, with an increasing number of fibers but a constant degree of reinforcement, that is to say : if $V_f$ is the total volume of fibers and $V$ that of the cylinder,

$$\frac{V_f}{V} = \text{constant.}$$

This condition is very logic on a physical point of view and very comfortable for the mathematics, because it ensures the fixity of $P_f$ and $P_m$ when $\varepsilon$ is decreasing.

**Remark 2** : A lot of other problems are entering in the ($P_\varepsilon$) scope, for example, the diffusion equation in steady cases.

### 3. - THE HOMOGENIZATION METHOD.

3.1. **Theorem** : With the hypothesis mentioned in section 2.1. about $a_{ij}$ and $f$, if $u_\varepsilon$ is the unique solution of

$$ (P_\varepsilon) \quad \left\{ \begin{array}{l} A^\varepsilon u_\varepsilon = f \text{ in } \Omega \\ u_\varepsilon = 0 \text{ on } \partial\Omega \end{array} \right. $$

Then

$$u_\varepsilon \to u \text{ weakly in } H_0^1(\Omega),$$

where $u$ is the unique solution of

$$ (P) \quad \left\{ \begin{array}{l} Au = f \text{ in } \Omega \\ u = 0 \text{ on } \partial\Omega \end{array} \right. $$
and $\mathcal{A}$ the differential operator defined by

$$ (8) \qquad \mathcal{A} = - q_{ij} \frac{\partial^2}{\partial x_i \partial x_j} , $$

$q_{ij}$ being here constant coefficients defined by the way mentioned in the following section.

For the proof of this theorem see BENSOUSSAN-LIONS-PAPANICOLAU [4], [5] and TARTAR [13].

## 3.2. Construction of $\mathcal{A}$.

Let us remind the definition of the Sobolev space $H^1(\Omega)$.

$$ H^1(\Omega) = \{ \phi | \phi \in L^2(\Omega), \frac{\partial \phi}{\partial x_i} \in L^2(\Omega) \text{ i=1,...,N} \} , $$

where

$$ L^2(\Omega) = \{ \phi | \int_{\Omega} |\phi|^2 dx < \infty \} . $$

We introduce the space $W$ defined by

$$ (9) \qquad W = \{ \phi | \phi \in H^1(\Omega), \phi \text{ is P-periodic, i.e. takes equal values on} \\ \qquad \text{opposite sides of P} \} , $$

and the "reduced" operator $A_p$ defined by

$$ A_p \phi = - \frac{\partial}{\partial y_i} (a_{ij}(y) \frac{\partial \phi}{\partial y_j}) . $$

The bilinear form associated to $A_p$ is:

$$ (10) \qquad a_p(\phi, \psi) = \int_p a_{ij}(y) \frac{\partial \phi}{\partial y_j} \frac{\partial \psi}{\partial y_i} dy . $$

We define the functions $\chi^i \in W$ by

$$ (11) \qquad a_p(\chi^i, \psi) = a_p(y_i, \psi) \quad \forall \psi \in W, \chi^i \in W. $$

The equation (11) defined $\chi^i$ up to an additive constant. Then $q_{ij}$ is obtained by

$$ (12) \qquad q_{ij} = \frac{1}{V(P)} a_p(\chi^j - y_j, \chi^i - y_i), $$
where $V(p)$ is the volume of the parallelepiped $P$.

**Remark 3** : It is clear, right now, that the $q_{ij}$ can be interpreted like some "homogenized" or "effective" coefficients for the considered (mechanical, thermal, electrical,...) problem. These constant coefficients are in fact independent of the shape of $\Omega$ and of the boundary conditions, they are essentially dependent of :

• The representative cell of the heterogeneous medium.

• The differential operator.

**Remark 4** : In fact, $u$ is the first term of an asymptotic expansion :

$$ (13) \qquad u_{\varepsilon}(x) = u(x) - \varepsilon \left[ \chi^{i} \left( \frac{x}{\varepsilon} \right) \frac{\partial u}{\partial x_{i}} + \tilde{w}(x) \right] + O(\varepsilon^2), $$

whose the first part of the first order term can be computed without too much difficulties. One can found the following result in BENSOUSSAN-LIONS-PAPANICOLAUO [4]

$$ \| u_{\varepsilon} - u + \varepsilon \left[ \chi^{i} \left( \frac{x}{\varepsilon} \right) \frac{\partial u}{\partial x_{i}} \right] \| \to 0 \text{ strongly in } H^1(\Omega). $$

**Remark 5** : It is instructive to have a look on what are $\chi_i$ and $q_{ij}$ in the above mentioned problem of torsion :

• $\chi_1$ is the solution of the following transmission problem :

$$ \Delta \chi_1 = 0 \text{ in } P_f $$

$$ \Delta \chi_1 = 0 \text{ in } P_m $$

and, if $[\phi]_{\partial P_f}$ means "jump of the quantity $\phi$" across the surface $\partial P_f$ and $n$ is the unitary normal on $\partial P_f$ oriented toward the interior of $P_f$, we have

$$ [\chi_1]_{\partial P_f} = 0 $$

$$ \left[ \frac{1}{\mu} \frac{\partial x_1}{\partial n} \right]_{\partial P_f} = \left[ \frac{1}{\mu} \right]_{\partial P_f} \cos(\vec{n}, \vec{y}_1) . $$

• By virtue of (11) and (12), the "homogenized coefficients" are given by :

$$ q_{ij} = \frac{1}{V} a_p(y_j, y_i) - \frac{1}{V} a_p(\chi_j, y_i) . $$
If we calculate the first part of $q_{ij}$ we obtain

$$ \frac{1}{V} a_p(y_j, y_i) = \frac{1}{\mu_f} \frac{V_f}{V} + \frac{1}{\mu_m} \frac{V-V_f}{V}, $$

and we notice that it represents exactly the law of mixtures.

### 3.3. Advantages of the homogenization :

They are almost obvious

a) While the solution $u_{\epsilon}$ of $(P_{\epsilon})$ is generally impossible to compute as soon as $\epsilon$ is small, $u$ can be obtained by solving simple problems independent of $\epsilon$ :

The $\chi^i$ are solutions of elliptic equations on the very simple domain P.

The $q_{ij}$ are obtained by integration on the same domain P.

Lastly, $u$ is given by a Dirichlet problem on $\Omega$ but now with constant coefficients.

b) The $q_{ij}$ have to be computed only once for a whole family of problems corresponding to the same differential operator.

## 4. - NUMERICAL RESULTS.

### 4.1. Description of the worked out computations :

We suppose that we are looking for the stresses in a cylindrical fibers reinforced bar B, subjected to a torque, the data for this problem being chosen as following :

* $\Omega$, square cross section of B with sides of unit length.

* Fibers of glass ($\mu_f=4.10^6$ psi) or Boron ($\mu_f=2,5.10^7$ psi), with square cross section such that

$$ \frac{V_f}{V} = \frac{1}{9} $$

* matrix epoxy : $\mu_m = 2,2.10^5$ psi.

The representative cell is the unit cell shown on Fig. 4.

We have computed :

* The real solution $\theta_{\epsilon}$ for $\epsilon = \frac{1}{2}$, $\epsilon = \frac{1}{4}$ and $\epsilon = \frac{1}{8}$.
<figure><img src="image_6.png" /><figcaption>Figure 4 : Representative cell.</figcaption></figure>

The solution θ of the associated (P) problem.

The corrected solutions :

$$ (14) \qquad \tilde{\Theta}_{\varepsilon}(x) = \theta(x) - \varepsilon \left\{ \chi_i \left(\frac{x}{\varepsilon}\right) \frac{\partial u}{\partial x_i} + \tilde{W}(x) \right\} $$

with $\tilde{W} \equiv 0$ here because of the symmetries of the basic cell P.

The following given curves are the respective sections of the surfaces $\theta_{\varepsilon}(x)$, $\theta(x)$ and $\tilde{\theta}_{\varepsilon}(x)$ by a diagonal plan (or a "parallel to the sides" one) perpendicular to the domain $\Omega$ (see Fig. 5).

## 4.2. Some comments about the following results :

### 4.2.1. For the "glass-epoxy" composite :

a) The homogenized coefficient obtained is $\bar{\mu} = 2,68 \cdot 10^5$ psi while the law of mixture would have given $\bar{\mu} = 2,46 \cdot 10^5$ psi.

b) For $\varepsilon = \frac{1}{2}$, which is not small, the difference between $\theta_{\varepsilon}$ and $\theta$ is relatively big (see Fig. 6 and 9), particularly in the neighbourhood of the fibers. However, the corrected term $\tilde{\theta}_{\varepsilon}$ is already much better if we remember that, in fact, the stress tensor is given by the components of grad $\theta_{\varepsilon}$.

c) When $\varepsilon$ is becoming smaller (cf. Figures 7, 8, 10, 11) the results are better and better, specially for the slopes of $\theta_{\varepsilon}$ and $\tilde{\theta}_{\varepsilon}$. We can see that the theorem of Section 3.1 is numerically verified although $\varepsilon$ is not very small.
$$ \epsilon = \frac{1}{2} $$

diagonal sections

$$ \epsilon = \frac{1}{4} $$

parallel sections

$$ \epsilon = \frac{1}{8} $$

<figure><img src="image_7.png" /><figcaption>Fig. 5 : Distribution of heterogeneities in the different treated cases :</figcaption></figure>
### 4.2.2. For the "Boron-Epoxy" composite :

a) The homogenized coefficient is $\mu = 2,75 \times 10^5$ in place of $\bar{\mu} = 2,47 \times 10^5$ by the law of mixtures.

b) We cannot present the curves for $\epsilon = \frac{1}{2}$ because the tops reached by the solution $\theta_{\epsilon}$ on the fibers are too much high, but, for $\epsilon = \frac{1}{4}$ (see Figures 12, 14) which is not so small, the solutions $\theta$ and $\tilde{\theta}_{\epsilon}$ give already very good approaches of $\theta_{\epsilon}$ and still better of grad $\theta_{\epsilon}$. For $\epsilon = \frac{1}{8}$, $\theta_{\epsilon}$ is very difficult to reach directly but the homogenized results seem to be satisfying.

**Remark 6** : The scale of the curves in ordinate is not significant because, in fact, $\theta_{\epsilon}$, $\theta$ and $\tilde{\theta}_{\epsilon}$ are proportional to $\alpha$ ; we have here chosen $\alpha\mu_m=5$ to obtain curves of sufficient size.

**Remark 7** : To be able to compare our results to the experimental ones given in the book [9] p. 68, we made again the same computing for $\frac{V_f}{V} = 0.7$ and we obtained

$\mu = 0.965 \times 10^6$ for the Glass-Epoxy in place of $\bar{\mu} = 0.65 \times 10^6$ by the mixture law and $\mu = 0.95 \times 10^6$ in [9].

**Remark 8** : The "parallel sections" (cf. Figures 9,10,11,14 and 15) allow to determine directly $\sigma_{32} = -\frac{\partial\theta}{\partial x_1}$ which is nothing else in each point, as the slope of the curves obtained. Because of the symmetries of the domain chosen, it is possible to deduct of that, $\sigma_{31} = \frac{\partial\theta}{\partial x_2}$ at the symmetric points with respect to the diagonales of $\Omega$.

**Remark 9** : The corrected solution obtained by formula (14) is not sufficient to approach the picks (see fig. 12 for example). For a better approximation of the picks see J.L. LIONS [11].
<figure><img src="image_8.png" /><figcaption>Figure 6 : $\epsilon = \frac{1}{2}$. COMPOSITE "GLASS-EPOXY". (diagonal section)</figcaption></figure>

$$ \left. \begin{array}{l} \mu_f = 4.10^6 \\ \\ \mu_m = 2,2.10^5 \end{array} \right\} \mu = 2,68.10^5 \text{ for } \frac{V_f}{V} = \frac{1}{9} $$
<figure><img src="image_9.png" /><figcaption>Figure 7 : $\epsilon = \frac{1}{4}$. "GLASS-EPOXY". (diagonal section)</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>Figure 8 : ε = 1/8 . "GLASS EPOXY" (diagonal section)</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>Figure 9 : $\epsilon = \frac{1}{2}$. "GLASS-EPOXY". (parallel section)</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>Figure 10 : $\epsilon = \frac{1}{4}$. "GLASS-EPOXY" (parallel section)</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>Figure 11 : $\epsilon = \frac{1}{8}$. "GLASS-EPOXY". (parallel section)</figcaption></figure>
<figure><img src="image_14.png" /><figcaption>Figure 12 : $\epsilon = \frac{1}{4}$. Composite "Boron Epoxy". (diagonal sections)</figcaption></figure>

$$ \left. \begin{array}{l} \mu_f = 2,5.10^7 \\ \\ \mu_m = 2,2.10^5 \end{array} \right\} \quad \mu = 2,75.10^5 \text{ for } \frac{V_f}{V} = \frac{1}{9} $$
<figure><img src="image_15.png" /><figcaption>Figure 13 : $\epsilon = \frac{1}{8}$. "BORON-EPOXY". (diagonal section)</figcaption></figure>
<figure><img src="image_16.png" /><figcaption>Figure 14 : $\epsilon = \frac{1}{4}$. "Boron-Epoxy". (parallel section)</figcaption></figure>
<figure><img src="image_17.png" /><figcaption>Figure 15 : $\epsilon = \frac{1}{8}$. "Boron-Epoxy". (parallel section)</figcaption></figure>
## 5. - SOME ELEMENTS ABOUT THE MATHEMATICAL TECHNIQUES USED TO OBTAIN THE PREVIOUS RESULTS :

The hypothesis and notations are those of sections 2.1 and 3.1.

### 5.1. Multiple scales method (see BENSOUSSAN-LIONS-PAPANICOLAU [2], [4], [5])

The solution $u_{\epsilon}$ has no reasons to be periodic itself because $\Omega$ is bounded and $\partial\Omega$ is not necessarily coinciding with the boundary of some basic cells. However, $u_{\epsilon}$ is obviously depending of the period of distribution of the heterogeneities ; for these reasons we are going to look for a solution of the form :

$$ (15) \qquad u_{\epsilon}(x) = W_{0}(x,y) + \epsilon W_{1}(x,y) + \epsilon^{2} W_{2}(x,y) + \dots $$

$$ \text{where} \qquad y_i = \frac{x_i}{\epsilon} \quad \text{and} \quad (x,y) \in \Omega \times P. $$

If we set generally : $\phi(x,y) = \phi^*(x)$, we have

$$ \frac{\partial \phi^*}{\partial x_i} = \frac{\partial \phi}{\partial x_i} + \frac{1}{\epsilon} \frac{\partial \phi}{\partial y_i} \quad \text{and then} $$

$$ A_{\epsilon}^{\alpha} u_{\epsilon} = f $$

can be written

$$ (16) \qquad (\epsilon^{-2} A_1 + \epsilon^{-1} A_2 + A_3) (W_0 + \epsilon W_1 + \epsilon^2 W_2 \dots) = f $$

with

$$ (17) \qquad A_1 = - \frac{\partial}{\partial y_i} \left[ a_{ij} \frac{\partial}{\partial y_j} \right] $$

$$ (18) \qquad A_2 = - \frac{\partial}{\partial y_i} \left[ a_{ij} \frac{\partial}{\partial x_j} \right] - a_{ij} \frac{\partial^2}{\partial x_i \partial y_j} $$

$$ (19) \qquad A_3 = - a_{ij} \frac{\partial^2}{\partial x_i \partial x_j} $$

By identification in (16) for the different orders of $\epsilon$, we obtain successively:

$$ (20) \qquad \epsilon^{-2} : \qquad A_1 W_0 = 0 $$

$$ (21) \qquad \epsilon^{-1} : \qquad A_1 W_1 + A_2 W_0 = 0 $$

$$ (22) \qquad \epsilon^0 : \qquad A_1 W_2 + A_2 W_1 + A_3 W_0 = f. $$
In (20) x plays the role of a parameter, consequently for all x, $W_o$ is a constant, i.e.

$$ (23) \qquad W_o(x,y) = u(x) . $$

Using (23), (21) reduces to

$$ A_l W_l = \frac{\partial}{\partial y_i} (a_{ij}(y)) \frac{\partial u}{\partial x_j} . $$

If we introduce $\chi^j$ by (11), i.e.

$$ A_l \chi^j = - \frac{\partial}{\partial y_i} (a_{ij}(y)) , $$

we obtain

$$ (24) \qquad W_l(x,y) = - \chi^j(y) \frac{\partial u}{\partial x_j}(x) + \tilde{W}_l(x) . $$

We have now to solve (22). Using the Fredholm's alternative, (22) has a solution iff

$$ (25) \qquad \int_p (f - A_2 W_l - A_3 W_o) dy = 0 . $$

Using (18), (19), (24), in (25) we obtain

$$ (26) \qquad \left[ - \int_p a_{ij}(y) dy - \int_p a_{ik} \frac{\partial \chi^j}{\partial y_k}(y) dy \right] \frac{\partial^2 u}{\partial x_i \partial x_j} = f . $$

The coefficient of $\frac{\partial^2 u}{\partial x_i \partial x_j}$ in (26) equals

$$ a_p(\chi^j - y_j, -y^i), $$

and since $\chi^i$ is periodic we have

$$ a_p(\chi^j - y_j, \chi^i - y_i) = 0 . $$

So we obtain (12).
## 5.2. Energy method. (See L. TARTAR [13])

We consider the variational formulation of the problem $(\rho_{\epsilon})$:

$$a^2(u_{\epsilon}, v) = (f, v) = (f, v) \quad \forall v \in H_0^1(\Omega),$$

where $a^\epsilon(u_\epsilon, v) = \int_{\Omega} a_{ij} \left(\frac{x}{\epsilon}\right) \frac{\partial u_\epsilon}{\partial x_j} \frac{\partial v}{\partial x_i} dx$,

and we introduce $\xi_i^\epsilon$ by

$$\xi_i^\epsilon = a_{ij} \left(\frac{x}{\epsilon}\right) \frac{\partial u_\epsilon}{\partial x_j}$$

by virtue of (4), $u_\epsilon$ is bounded in $H^1(\Omega)$ and $\xi_i^\epsilon$ is bounded in $L^2(\Omega)$, so we have

$$u_{\epsilon} \to u \text{ weakly in } H_0^1(\Omega)$$

$$\xi_i^\epsilon \to \xi_i \text{ weakly in } L^2(\Omega).$$

We can prove by a non trivial way that

$$\xi_i = q_{ij} \frac{\partial u}{\partial x_j}$$

and that $u$ is the solution of

$$(Au, v) = (f, v) \quad \forall v \in H_0^1(\Omega)$$

where

$$(Au, v) = q_{ij} \int_{\Omega} \frac{\partial u}{\partial x_j} \frac{\partial v}{\partial x_i} dx.$$

**Remark 10** : In order to fix the ideas we have considered the problem in $H_0^1(\Omega)$ but the results are the same in a space $V$ such as $H_0^1(\Omega) \subset V \subset H^1(\Omega)$.

**Remark 11** : The "multiple scales" method is practical because naturally constructive. It allows in particular to obtain the first order corrector, but this method needs regularity of the coefficients. On the other hand, the "energy method" needs only bounded coefficients but does not introduce naturally the correctors. So we can judge the very complementary character of those two mathematical techniques.
## 6. - ELEMENTS ABOUT NUMERICAL TECHNIQUES :

For solving the boundary value problems introduced in this paper, we use the finite element method.

The domains P and Ω are divided into triangles and the spaces H¹(Ω) and H¹(P) are approached by :

$$V_h = \{v | v \in C^0(\Omega), \text{ or } v \in C^0(P) \text{ polynomials of degree } 1 \text{ on each triangle}\}$$

The functions θ and θε are approached in the space

$$V_{oh} = \{v \in V_h | v = 0 \text{ on } \partial\Omega\}$$

and the function χᵢ in the space

$$W_h = \{\phi \in V_h | \phi \text{ "P-periodic," } \phi=0 \text{ at the four corners of P}\}$$

The triangulation of P, to compute χᵢ, and the triangulation of Ω to compute θε are such as the discontinuities of the coefficients aᵢⱼ coincide with sides of triangles.

Generally a few hundred triangles are sufficient to compute χᵢ because the geometry of P is very simple while, to compute θε several thousands are necessary in order to approximate the numerous discontinuities of aᵢⱼ.

For θ which is solution of an operator with constant coefficients, a few number of triangles is sufficient but for $\frac{\partial\theta}{\partial x_i}$, which appears in the first corrective term, we need about a thousand of triangles.

For the elastic torsion example mentioned above, we took

288 triangles for computing χᵢ
1152 triangles for computing θ
4608 triangles for computing θε (see fig. 16).

The linear systems obtained after approximation of the boundary value problems are solved by the overrelaxation method with optimal parameter.

The computing times obtained with an IBM 370-168 are

1.5 s. to compute the coefficients of the homogenized operator,
4 s. to compute $\theta$,
30 s. to compute $\theta_{\epsilon}$.

In order to compute the first corrective term we approximate $\frac{\partial\theta}{\partial x_i}$ by a function of $V_f$ defined on each node by

$$ \frac{\partial\theta}{\partial x_i}(M) = \frac{1}{N_v(M)} \sum_{T \in v(M)} \left[ \frac{\partial\theta}{\partial x_i} \right]_T $$

where

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

$M$ is a node of the triangulation,

$v(M)$ is the set of triangles whose $M$ is a vertex,

$N_v(M)$ is the number of triangles of $v(M)$,

$[\frac{\partial\theta}{\partial x_i}]_T$ is the value of the derivative of $\theta$ on the triangle $T$.

We can use this formula because $\theta$ is very regular as a solution of a problem with constant coefficients and constant second member.

## CONCLUSIONS.

It follows from the previous results that the method of homogenization is an effective method for solving problems in periodic heterogeneous medium as soon as the number of heterogeneities is large. Furthermore it can be extended to numerous directions :

*Neumann problem : for instance*

$$ \left\{ \begin{array}{l} A^{\epsilon} u_{\epsilon} + a_{0} \left(\frac{x}{\epsilon}\right) u_{\epsilon} = f \quad \text{on } \Omega \\[1em] \frac{\partial u_{\epsilon}}{\partial v} = 0 \quad \text{on } \partial\Omega \\[1em] A^{\epsilon} \frac{\partial u_{\epsilon}}{\partial x} = 0 \quad \text{on } \partial\Omega \end{array} \right. \quad (\text{see Figure 17}) $$

*Variational inequalities* : that is to say, when the solution has to be found in a convex subset of a vectorial space for instance.

$$ K = \{ \phi \mid \phi \in H_0^1(\Omega) , \phi \le 0 \text{ in } \Omega \} $$

It is interesting to note, on the Figure 18, that the free boundary (unknown of the problem) between
$$ \Omega_{\epsilon}^{0} = \{x | x \in \Omega , u_{\epsilon}(x) = 0\} \quad \text{and} $$

$$ \Omega_{\epsilon}^{-} = \{x | x \in \Omega , u_{\epsilon}(x) < 0\}, $$

is given with a good precision by the homogenized solution u.

### Operators of higher orders (cf. G. DUVAUT [7])

#### Evolution problems

#### Systems of partial differential equations: for instance in elasticity (cf. G. DUVAUT [8]) or electromagnetism.

#### Operators with coefficients depending more generally on $\epsilon$:

$$ a_{ij}(x, \frac{x}{\epsilon}); a_{ij}(x, \frac{x}{\epsilon}, \frac{x}{\eta}). $$

Several physical examples are corresponding to these cases (cf. A. BENSOUSSAN-LIONS-PAPANICOLAU [4]).

#### Porous media

#### Mechanics of fluid suspensions.
<figure><img src="image_19.png" /><figcaption>Fig. 16 : Triangulation of Ω. 4608 triangles.

ε = 1/8 . The fibers are hachured.</figcaption></figure>
<figure><img src="image_20.png" /><figcaption>Figure 17 : Neumann problem. (diagonal section)</figcaption></figure>
<figure><img src="image_21.png" /><figcaption>Figure 18 : Variational Inequalities. (diagonal section)</figcaption></figure>

$u_{\varepsilon} < 0$
REFERENCES

I. BABUSKA

[1] "Homogenization and its applications. Mathematical and computational problems". Tech. Note BN-821 (July 1975). Institute for Fluid Dynamics and Applied Mathematics, University of Maryland.

A. BENSOUSSAN-J.L. LIONS-G. PAPANICOLAUO.

[2] "Sur quelques phénomènes asymptotiques stationnaires". C.R. Acad. Sc. Paris, t. 281, (21.7.1975).

[3] "Sur quelques phénomènes asymptotiques d'évolution". C.R. Acad. Sc. Paris, t. 281, (1975), p. 317-322.

[4] "Homogénéisation, correcteurs et problèmes non-linéaires". C.R. Acad. Sc., Paris (1976).

[5] Book in preparation - North Holland.

E. DE GIORGI - S. SPAGNOLO

[6] "Sulla convergenza delli integrali dell'energia per operatozi ellittici del secondo ordine". Boll. Uni. Math. Ital. 8, 1973, p. 391-411.

G. DUVAUT

[7] "Homogénéisation d'une plaque mince en flexion à structure périodique". Will appear in C.R. Acad. Sc. Paris.

[8] Invited lecture at 14th International Congress of Theoretical and Applied Mechanics, Delft 1976. Will appear in the proceedings.

S.K. GARG - V. SVATBONAS - G.A. GURTMAN

[9] Analysis of Structural Composite Materials. Marcel Dekker Inc., New-York 1973.

J.L. LIONS

[10] Cours au Collège de France, Paris 1975-1976.

[11] "Remarks on some Numerical Aspects of the Homogenization Method in Composite Materials". Colloque Franco Soviétique, Novossibirsk 1976. (To appear in 1977 Dunod).
E. SANCHEZ PALENCIA

[12] "Comportement local et macroscopique d'un type de milieux physiques hétéro-gènes". Int. J. Eng. Sci., 1974, Vol. 12 p. 331-351.

L. TARTAR

[13] To appear.