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

HAL
open science

# Numerical experiments of the homogenization method for operators with periodic coefficients

Jean-François Bourgat

► To cite this version:

Jean-François Bourgat. Numerical experiments of the homogenization method for operators with periodic coefficients. [Research Report] IRIA-RR-277, IRIA. 1978, pp.29. hal-04716545

HAL Id: hal-04716545

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

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

Distributed under a Creative Commons Attribution 4.0 International License
RIA

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

NUMERICAL EXPERIMENTS OF
THE HOMOGENIZATION METHOD
FOR OPERATORS
WITH PERIODIC COEFFICIENTS

Jean-François BOURGAT

Rapport de Recherche N° 277

Janvier 1978

29p.
NUMERICAL EXPERIMENTS OF THE HOMOGENIZATION METHOD
FOR OPERATORS WITH PERIODIC COEFFICIENTS

Jean-François Bourgat

Résumé :

La résolution numérique d'un problème aux limites dans un milieu hétérogène est impossible 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». Nous présentons ici, dans le cas particulier d'une distribution périodique d'hétérogénéités, quelques résultats numériques sur le calcul de l'opérateur homogénéisé, l'approximation de la solution par la solution homogénéisée et sur l'efficacité de «correcteurs» pour approcher le comportement périodique de la solution.

Abstract :

A boundary value problem in an heterogeneous medium is numerically impossible to solve if the number of heterogeneities is very large. The homogenization consists in replacing 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. We present here some numerical experiments on the computation of the homogenized operator, the approximation of the real solution by the homogenized one and the efficiency of «correctors» to approximate the periodic behaviour of the solution.
<figure><img src="image_3.png" /></figure>
INTRODUCTION

Let us consider a domain Ω, bounded open set of $\mathbb{R}^n$ with $n=1,2,3$, which consists in a composite material with a periodic structure. Let us denote by $\varepsilon$ the period of the structure (we take here the same length $\varepsilon$ of the period in all directions for the sake of clearness but this need not be the case in general). We suppose that $\varepsilon$ is "small" that is to say the length of the period is small with respect to the size of $\Omega$. In each cell we have two regions in which the physical constants are different (see fig. 1).

Let us make the scaling $\frac{1}{\varepsilon}$ on a given cell, we obtain the "unit cell" $Y = ]0,1[\mathbb{R}^n$ and in $Y$ we have two regions $Y_1$ and $Y_0$ (see fig. 1).

If the material which corresponds to $Y_0$ (resp. $Y_1$) is modelled by the operator

$$ - \sum_{i,j=1}^{n} a_{ijo} \frac{\partial^2}{\partial y_i \partial y_j} \quad (\text{resp.} \ - \sum_{i,j=1}^{n} a_{ij1} \frac{\partial^2}{\partial y_i \partial y_j}) , $$

<figure><img src="image_4.png" /><figcaption>Figure 1 - A cell structure and the representative cell</figcaption></figure>
we are led to introduce

$$ a_{ij}(y) = a_{ijo} \text{ in } Y_0, a_{ij1} \text{ in } Y_1 $$

and then the state of the system is given by the solution of

$$ - \sum_{i,j=1}^{n} \frac{\partial}{\partial x_i} \left( a_{ij} \left( \frac{x}{\epsilon} \right) \frac{\partial u_{\epsilon}}{\partial x_j} \right) = f \quad \text{in } \Omega, $$

where f is a given function and where $u_{\epsilon}$ is subject to some boundary condition on $\Gamma = \partial\Omega$.

We now have the following questions :

- How to compute an operator with constant coefficients which models the macroscopic behaviour of the material.

- How to compute $u_{\epsilon}$ when the number of cells is very large.

- How to obtain an approximation of the microperiodic behaviour of $u_{\epsilon}$.

In this paper we present some numerical results obtained in these directions.

Recent results on the mathematical aspects of the homogenization have been given by BABUSKA [2], BENSOUSSAN-LIONS-PAPANICOLAU [1], [2], DE GIORGI-SPAGNOLO [1], SANCHEZ PALENCIA [1] and TARTAR [1].

For this report we use theoretical results of BENSOUSSAN-LIONS-PAPANICOLAU [2].
# 1. - THE HOMOGENIZED OPERATOR

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

## 1.1. Statement of the basic result.

We consider functions $a_{ij}(y)$ defined in $\mathbb{R}^n$, with real values, such that

(1) $a_{ij} \in L^\infty(\mathbb{R}^n)$, $a_{ij}$ is $Y$-periodic, (i.e. takes equal values on opposite sides of $Y$).

and

(2) $a_{ij}(y) \xi_i \xi_j \ge \alpha \xi_i \xi_j, \alpha > 0, \forall \xi \in \mathbb{R}^n$.

We define, for $\varepsilon > 0$, the operator $A^\varepsilon$ by

$$ (3) \qquad A^\varepsilon \phi = - \frac{\partial}{\partial x_i} \left( a_{ij} \left( \frac{x}{\varepsilon} \right) \frac{\partial \phi}{\partial x_j} \right) $$

and we consider the equation

$$ (4) \qquad \begin{cases} A^\varepsilon u_\varepsilon = f \text{ in } \Omega, f \in L^2(\Omega) \\ u_\varepsilon = 0 \text{ on } \Gamma. \end{cases} $$

The operators $A^\varepsilon$ are "uniformly elliptic in $\varepsilon$" in the following sense. For $u,v$ in $H_0^1(\Omega)$ we set

$$ (5) \qquad a^\varepsilon(u,v) = \int_\Omega a_{ij} \left(\frac{x}{\varepsilon}\right) \frac{\partial u}{\partial x_j} \frac{\partial v}{\partial x_i} \, dx, $$

then

$$ (6) \qquad a^\varepsilon(v,v) \ge \alpha ||v||^2_{H_0^1(\Omega)} $$

So, as $\varepsilon \to 0$, one has

$$ (7) \qquad \|u_\varepsilon\|_{H_0^1(\Omega)} \le \frac{\|f\|_{L^2(\Omega)}}{\alpha} $$

and therefore we can extract a subsequence, still denoted by $u_\varepsilon$, such that

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

Then one has the following result (see BENSOUSSAN-LIONS-PAPANICOLAU [2], TARTAR [1])
**THEOREM 1** : It exists a second order elliptic operator $A$ with constant coefficients
$q_{ij}$ such that $u$ is the solution of the equation

$$
(9) \qquad \left\{
\begin{array}{l}
-q_{ij} \displaystyle\frac{\partial^2 u}{\partial x_i \partial x_j} = f \quad \text{in } \Omega \\
u = 0 \quad \text{on } \Gamma.
\end{array}
\right.
$$

$\mathcal{A}$ is the "homogenized operator", $q_{ij}$ the "homogenized coefficients", (9) the "homogenized problem" and $u$ the "homogenized solution".

Construction of $\mathcal{A}$

We introduce the space W defined by

$$
(10) \quad W = \{ \psi \in H^1(Y) \mid \psi \text{ is Y-periodic (i.e. takes equal values on opposite sides of Y)} \}
$$

and the "reduced" operator $A_Y$ defined by

$$
(11) \quad A_Y \phi = - \frac{\partial}{\partial y_i} \left( a_{ij}(y) \frac{\partial \phi}{\partial y_j} \right).
$$

The bilinear form associated to $A_Y$ is

$$
(12) \quad a_Y(\phi, \psi) = \int_Y a_{ij}(y) \frac{\partial \phi}{\partial y_j} \frac{\partial \psi}{\partial y_i} \, dy
$$

and we consider the function $\chi^i \in W$ defined by

$$
(13) \qquad a_Y(\chi^i - y_i, \psi) = 0 \quad \forall \psi \in W.
$$

We remark that (13) defines $\chi^i$ up to an additive constant. We choose $\chi^i$ such that the mean value of $\chi^i$ is zero (i.e. $\int_Y \chi^i(y) \, dy = 0$). This choice will be justified in section 3.

Then the "homogenized coefficients" $q_{ij}$ are obtained by

$$
(14) \qquad q_{ij} = \frac{1}{|Y|} a_Y(x_j^i - y_j, x_i^i - y_i),
$$

where $|Y|$ denotes the measure of Y.

1.2. Asymptotic expansion

Let us briefly recall the multiple scale method (see BENSOUSSAN-LIONS-PAPANICOLAU
[2]) which gives formulas (14).

We set
$$
(15) \quad y_i = x_i / \epsilon \quad i=1, \dots, n.
$$

Applied to a function of x and y the operator $\frac{\partial}{\partial x_i}$ becomes $\frac{\partial}{\partial x_i} + \frac{1}{\epsilon} \frac{\partial}{\partial y_i}$ and we can
write

$$
(16) \quad A^{\epsilon} = \epsilon^{-2}A_{1} + \epsilon^{-1}A_{2} + A_{3},
$$

with

$$
(17) \quad A_1 = - \frac{\partial}{\partial y_i} \left( a_{ij}(y) \frac{\partial}{\partial y_j} \right),
$$

$$
(18) \quad A_2 = - \frac{\partial}{\partial y_i} \left( a_{ij}(y) \frac{\partial}{\partial x_j} \right) - \frac{\partial}{\partial x_i} \left( a_{ij}(y) \frac{\partial}{\partial y_j} \right),
$$

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

We look for $u_{\epsilon}$ in the form

$$
(20) \quad u_{\epsilon}(x) = w_{0}(x,y) + \epsilon w_{1}(x,y) + \epsilon^{2} w_{2}(x,y) + \dots
$$

where

$w_j(x,y)$ is $Y$-periodic in $y$.

Then we can write (4) in the form

$$
(\epsilon^{-2}A_1 + \epsilon^{-1}A_2 + A_3)(w_o + \epsilon w_1 + \epsilon^2 w_2 + \dots) = f
$$

and if we identify in $\epsilon$ we obtain

$$
(21) \quad A_1 w_o = 0,
$$

$$
(22) \quad A_1 w_1 + A_2 w_o = 0,
$$

$$
(23) \quad \mathring{A}_1 w_2 + A_2 w_1 + A_3 w_o = f .
$$

If we notice that the only solution in Y of $A_1\theta = 0$, $\theta$ periodic, is $\theta = \text{constant}$,
we see that

$$
(24) \quad w_o(x,y) = u(x).
$$

Then equation (22) becomes

$$
(25) \quad A_1 w_1 = \frac{\partial}{\partial y_i} (a_{ij}(y) \frac{\partial u}{\partial x_j}(x)).
$$
Therefore if we introduce $\chi^j$ by (13), (i.e. $A_1\chi^j = - \frac{\partial}{\partial y_i} (a_{ij}(y))$, we obtain

$$ (26) \qquad w_1(x,y) = - \chi^j(y) \frac{\partial u}{\partial x_j}(x) + \tilde{w}_1(x). $$

Using the Fredholm's alternative (i.e. a necessary and sufficient condition, for $A_1\theta = F$, $\theta$ periodic, to have a solution is $\int_Y F(y)dy = 0$) equation (23) in $w_2$ has a solution iff

$$ (27) \qquad \int_Y (f - A_2 w_1 - A_3 w_0) dy = 0 \quad \forall x. $$

Using (18) and (19) we obtain

$$ - \frac{\partial}{\partial x_i} \int_Y a_{ij}(y) \frac{\partial w_1}{\partial y_i} dy - \left[ \int_Y a_{ij}(y) dy \right] \frac{\partial^2 u}{\partial x_i \partial x_j} = |Y| f $$

and replacing $w_1$ by its value, we obtain

$$ (28) \qquad \left( - \int_Y a_{ij}(y) dy + \int_Y a_{ik}(y) \frac{\partial \chi^j}{\partial y_k}(y) dy \right) \frac{\partial^2 u}{\partial x_i \partial x_j} = |Y| f. $$

Then the coefficient of $\frac{\partial^2 u}{\partial x_i \partial x_j}$ equals

$$ (29) \qquad q_{ij} = a_Y(\chi^j - y_j, -y_i) $$

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

$$ (30) \qquad a_Y(\chi^j - y_j, \chi^i) = 0 $$

and therefore we obtain (14).

### 1.3. Computation of the homogenized coefficients

The homogenized operator coefficients are computed only on the unit cell Y. For solving the problem (13) which gives the functions $\chi^i$ we use the finite element method with Lagrange finite elements of degree 1 on triangles. This computation need not a priori a fine mesh however, as the homogenized coefficients $q_{ij}$ are obtained from the first derivatives of the functions $\chi^i$, we have to use a sufficiently fine mesh of Y.

Thus, if we consider the unit cell representation as in Fig. 2, with

$Y_1$ = disk of area 0.25

$a_{ijo}(y) = \delta_{ij}$ (the Kronecker's index)

$a_{ij1}(y) = 10 \delta_{ij}$

<figure><img src="image_5.png" /><figcaption>Figure 2.</figcaption></figure>
and a triangulation such as the discontinuities of the coefficients $a_{ij}$ coincide with sides of triangles, we obtain the results shown in table 1.

<table><caption>Table 1 - Accuracy of the homogenized coefficients computations.</caption><thead><tr><td>Number of triangles</td><td>$q_{11}$</td></tr></thead><tbody><tr><td>32</td><td>1.494</td></tr><tr><td>128</td><td>1.512</td></tr><tr><td>512</td><td>1.516</td></tr><tr><td>2048</td><td>1.516</td></tr></tbody></table>

For solving the linear system associated to the discretized problem, it is very difficult to use direct methods because the periodicity conditions disturb the band structure of the matrix. So we have used the overrelaxation method with an optimal parameter.

In several examples studied here, the coefficients $a_{ij}$ have some properties of symmetry, as it is often the case in physical examples. Then the question is: does the homogenized coefficients have the same properties?

One has (see BOURGAT-DERVIEUX [1]) the

### PROPOSITION 1

If the coefficients $a_{ij}$ are such as

(31) $a_{ii}$ is symmetric with respect to the middle hyperplans $\Delta_i$ ($i=1,...,n$) of the representative cell (see fig. 3),

Then

(32) $a_{ij} = 0$ for $i \neq j \Rightarrow q_{ij} = 0$ for $i \neq j$.

<figure><img src="image_6.png" /><figcaption>Figure 3</figcaption></figure>
If we consider the representative cell
as in fig. 4, with

$$
a_{ijo} = \delta_{ij}, a_{ij1} = 10 \delta_{ij},
$$

we have $a_{ij} = 0$ if $i \neq j$ but the assumption (31) is not satisfied. The homogenized coefficients obtained

$$
q_{11} = q_{22} = 1.915, \quad q_{12} = q_{21} = -0.101,
$$

<figure><img src="image_7.png" /><figcaption>Figure 4</figcaption></figure>

show that $q_{ij} \neq 0$ for $i \neq j$. Therefore an heterogeneous medium made of two isotropic materials can be homogenized into an effective non isotropic medium.

We present now various examples showing some properties of the homogenized operator.

**Example 1 : Homogenized coefficients and mean values**

Let us consider the representative cell
as in fig. 5, with

$$
(33) \quad \begin{cases} a_{ijo} = \delta_{ij} \\ a_{ij1} = \lambda \delta_{ij}, \quad \lambda > 0. \end{cases}
$$

In this case one has

$$
q_{11} = q_{22}, q_{12} = q_{21} = 0.
$$

<figure><img src="image_8.png" /><figcaption>Figure 5</figcaption></figure>

The values of $q_{11}$ as $\lambda$ varies from 0 to $+\infty$ are represented on fig. 6 and compared to

$$
(34) \quad \mathcal{M}(a_{11}) = \frac{1}{|Y|} \int_Y a_{11}(y) dy, \text{ the mean value of } a_{11}(y)
$$

and also to

$$
(35) \quad \left( \mathcal{M} \left( \frac{1}{a_{11}} \right) \right)^{-1} = \left[ \frac{1}{|Y|} \int \frac{1}{a_{11}(y)} \right]^{-1}
$$

which is the homogenized coefficient if $a_{11}(x)$ is periodic in the direction $x_1$ only.

We remark in the figure 6 that $q_{11}$ is very different from $\mathcal{M}(a_{11})$ or $(\mathcal{M}(\frac{1}{a_{11}}))^{-1}$
and that $q_{11}$ converges as $\lambda \to 0$ or as $\lambda \to \infty$. We remark also that one has (see TARTAR
[1])

$$
(36) \quad \left[ \mathcal{M} \left( \frac{1}{a_{11}} \right) \right]^{-1} \le q_{11} \le \mathcal{M}(a_{11}) \quad \forall \lambda > 0.
$$
<figure><img src="image_9.png" /><figcaption>Figure 6 - Comparison of $q_{11}$, $m(a_{11})$, $[m(\frac{1}{a_{11}})]^{-1}$.</figcaption></figure>

**Example 2 : Torsion of fiber-reinforced bars** (see BOURGAT-LANCHON [1])

The equation (4) models the elastic torsion of a cylindrical fiber-reinforced bar. $\Omega$ is the cross section of the bar, $Y_1$ and $Y_0$ are the image in the scaling $\frac{1}{\epsilon}$ of the fiber and of the matrix respectively. $u_{\epsilon}$ is called the "stress function" and gives the components of the stress tensor.

The operator $A^{\epsilon}$ is modelled by

$$a_{ijo} = \delta_{ij}, \quad a_{ij1} = \lambda \delta_{ij},$$

with $\lambda = \frac{\mu_1}{\mu_o}$, where $\mu_1$ and $\mu_o$ are the shear modulus of the fiber and of the matrix respectively. The representative cell is chosen as in fig. 5.

For a "glass-epoxy" composite (glass in $Y_1$) with

$$(37) \qquad \mu_o = 2.2 \ 10^5 \text{ psi}, \mu_1 = 4 \ 10^6 \text{ psi} \ (\text{i.e. } \lambda = \frac{1}{18})$$

we obtain the homogenized shear modulus $\bar{\mu} = 2.68 \ 10^5$ psi while the law of mixture would have given $2.46 \ 10^5$ psi.

For a "boron-epoxy" composite (boron in $Y_1$) with

$$(38) \qquad \mu_o = 2.2 \ 10^5 \text{ psi}, \mu_1 = 2.5 \ 10^7 \text{ psi} \ (\text{i.e. } \lambda = \frac{1}{114})$$
we obtain $\bar{u} = 2.75 \ 10^5$ psi in place of $2.47 \ 10^5$ by the law of mixtures.

To be able to compare our results with the experimental ones given in GARG-SVATBONAS-GURTMAN [1], we make the same computation as above for a "glass epoxy" composite with $\overline{\text{meas}(Y_1)} = 0.7$. We obtain $\bar{u} = 0.965 \ 10^6$ psi in place of $0.65 \ 10^6$ psi by the mixture law and $0.95 \ 10^6$ psi in GARG-SVATBONAS-GURTMAN [1].

### Example 3 :

Let us consider the three representative cells of the figure 7. The area of $Y_1$ has the value $1/4$ in each cell. One can notice that each cell satisfies the symmetry property (31).

<figure><img src="image_10.png" /><figcaption>Figure 7</figcaption></figure>

For the following choice of $a_{ij}$

$$a_{ijo} = \delta_{ij}, \quad a_{ij1} = 10 \delta_{ij}$$

we have obtained with a very fine triangulation (2048 triangles)

$$q_{11} = \begin{cases} 1.516 & \text{for the disk} \\ 1.548 & \text{for the square} \\ 1.573 & \text{for the lozenge} \end{cases}$$

Therefore, in an application, if the shape of $Y_1$ is not well defined, the results of the homogenization computations depend strongly on the chosen shape.

### Example 4 :

Let us consider the representative cell $Y = ]0,a[ \times ]0,1[)$ as in fig. 8, with

$$a_{ijo} = \delta_{ij}, \quad a_{ij1} = 10 \delta_{ij}$$

The two materials are in the same pro-

<figure><img src="image_11.png" /><figcaption>Figure 8</figcaption></figure>
portion for every value of a.

Figure 9 shows how the homogenized coefficient $q_{11}$ varies as a varies from 0 to $\infty$.

<figure><img src="image_12.png" /><figcaption>Figure 9 - Variation of $q_{11}$ in function of the side a of Y = ]0,a[×]0,1[</figcaption></figure>

Let us introduce

$$ \mathcal{M}_{yi}(\phi) = \frac{1}{a_i} \int_{0}^{a_i} \phi(y_1, y_2) dy_i, \quad a_1 = a, a_2 = 1. $$

and

$$ \mathcal{M}_{yi}^{-1}(\phi) = \left( \frac{1}{a_i} \int_{0}^{a_i} \frac{1}{\phi(y_1, y_2)} dy_i \right)^{-1}. $$

We obtain

$$ \mathcal{M}_{y_2}^{-1} \circ \mathcal{M}_{y_1}(a_{11}) = 1.333, $$

$$ \mathcal{M}_{y_2} \circ \mathcal{M}_{y_1}^{-1}(a_{11}) = 1.143. $$

If we compare these numbers with the bounds of $q_{11}$ (see. fig. 9) we notice that

$$ q_{11} + \mathcal{M}_{y_2}^{-1} \circ \mathcal{M}_{y_1}^{-1}(a_{11}) \text{ as } a \to 0 $$

and

$$ q_{11} + \mathcal{M}_{y_2} \circ \mathcal{M}_{y_1}^{-1}(a_{11}) \text{ as } a \to \infty. $$

Thus, when the ratio of the sides of the representative cell Y is small or large we
can approximate $q_{11}$ by two successive homogenizations in one direction. This can be proved if we define the periodicity of $a_{ij}$ by $a_{ij} = \left(\frac{x_1}{\epsilon}, \frac{x_2}{\epsilon}\right)$ (see LIONS [1]).

## 2. - COMPARISON OF THE "TRUE SOLUTION" WITH THE "HOMOGENIZED SOLUTION"

We study here the approximation of $u_{\epsilon}$, solution of the equation (4), by the homogenized solution $u$, solution of equation (9).

Let us consider the unit cell Y as in fig. 5, with

$$ (39) \qquad a_{ij0} = \delta_{ij}, a_{ij1} = \lambda\delta_{ij}. $$

We take

$$ \Omega = ]0,1[ \times ]0,1[ , $$

$$ f = 10. $$

$$ \epsilon = 1/2, 1/4, 1/8. $$

Figure 10 shows the periodicity of the coefficients $a_{ij}$, in the domain $\Omega$, for $\epsilon = \frac{1}{8}$.

### 2.1. Numerical techniques

For solving equations (4) and (9), which give $u_{\epsilon}$ and $u$, we use the finite element method with Lagrangian elements of degree one on the triangles.

The homogenized solution which is the solution of a second-order elliptic equation with constant coefficients, is obtained without any difficulty. On the other hand the computation of $u_{\epsilon}$ needs a triangulation which must be very fine as $\epsilon$ becomes small.

Figure 10 shows a triangulation of $\Omega$ for $\epsilon = 1/8$, which is such that the discontinuities of the coefficients $a_{ij}$ coincide with sides of the triangles. The number of triangles is 4608. In each zone which corresponds to $Y_1$ (hachured zones), the number of triangles is only 8, this is a minimum to approach $u_{\epsilon}$ in these regions.

In the case $\lambda \ll 1$, 32 triangles are needed in this zone. However by using the symmetries of the problem, we can reduce it to a problem on a quarter of domain divided in 4608 triangles, which corresponds to a global triangulation with 18432 triangles.

For $\lambda = \frac{1}{114}$ and $\epsilon = \frac{1}{8}$ we obtain, in the middle M of a zone which corresponds to $Y_1$, the values of the approximated solution $u_{\epsilon h}$ given in table 2 (h is the step size).
<figure><img src="image_13.png" /><figcaption>**Figure 10** - Triangulation of Ω for ε = 1/8 - 4608 triangles - Hachured zones correspond to Y₁.</figcaption></figure>

<table><caption>**Table 2** - Convergence of u<sub>eh</sub>, as h → 0, for λ = 1/114, ε = 1/8.</caption><thead><tr><td>Number of triangles</td><td>Values of u<sub>eh</sub> at M</td></tr></thead><tbody><tr><td>1152</td><td>0.839</td></tr><tr><td>4608</td><td>0.996</td></tr><tr><td>18432</td><td>1.029</td></tr></tbody></table>

One can see that it is very difficult to compute the solution of a two dimensional problem in an heterogeneous medium as the number of heterogeneities is greater than one hundred.

For theoretical results on the convergence of problems with interfaces and singularities see BABUSKA [1].
## 2.2. Convergence of $u_{\varepsilon}$ towards $u$.

We have obtained in section 1 the following convergence result

$$ (40) \qquad u_{\varepsilon} \to u \text{ in } H_o^1(\Omega) \text{ weakly} $$

and, if coefficients $a_{ij}$ are sufficiently regular (see LIONS [1]), one has

$$ (41) \qquad \|u_{\varepsilon}-u\|_{C^0(\bar{\Omega})} \le C\varepsilon. $$

In this section we show numerical results on the convergence of $u_{\varepsilon}$ towards $u$ and the figures presented are representations of the solutions $u_{\varepsilon}$ and $u$ restricted along a diagonal of $\Omega$ as indicated on figure 11.

Figures 12, 13, 14 show $u_{\varepsilon}$ and $u$ respectively for $\lambda = 10, \frac{1}{18}, \frac{1}{114}$ and in each one $u_{\varepsilon}$ is represented for $\varepsilon = 1/2, 1/4, 1/8$. Values $\lambda = 1/18$ and $\lambda = 1/114$ correspond to the torsion of a fiber-reinforced bar presented in Section 1.3.

<figure><img src="image_14.png" /><figcaption>Figure 11</figcaption></figure>

One can see in figures 12, 13, 14 that $u_{\varepsilon}$ converges towards $u$ as $\varepsilon \to 0$ but also that the difference between $u_{\varepsilon}$ and $u$ is greater for $\lambda < 1$.

More precisely, we have indicated in table 3, for $\lambda = 10$ and $\lambda = \frac{1}{114}$, the error $u_{\varepsilon}-u$ in the norms of spaces $L^{\infty}(\Omega)$, $L^2(\Omega)$, $H_o^1(\Omega)$.

<table><caption>Table 3 - Convergence of the error $u_{\varepsilon}-u$ for $\lambda =10, \lambda = \frac{1}{114}$</caption><thead><tr><th rowspan="2">ε</th><th colspan="3">λ = 10</th><th colspan="3">λ = 1/114</th></tr><tr><th>Sup|u<sub>ε</sub>-u|</th><th>||u<sub>ε</sub>-u||<sub>L<sup>2</sup></sub><br>||u||<sub>L<sup>2</sup></sub></th><th>||u<sub>ε</sub>-u||<sub>H<sub>o</sub><sup>1</sup></sub><br>||u||<sub>H<sub>o</sub><sup>1</sup></sub></th><th>Sup|u<sub>ε</sub>-u|</th><th>||u<sub>ε</sub>-u||<sub>L<sup>2</sup></sub><br>||u||<sub>L<sup>2</sup></sub></th><th>||u<sub>ε</sub>-u||<sub>H<sub>o</sub><sup>1</sup></sub><br>||u||<sub>H<sub>o</sub><sup>1</sup></sub></th></tr></thead><tbody><tr><td>1/2</td><td>0.128</td><td>0.1</td><td>0.310</td><td>2.24</td><td>0.56</td><td>0.973</td></tr><tr><td>1/4</td><td>0.072</td><td>0.054</td><td>0.357</td><td>0.554</td><td>0.187</td><td>0.92</td></tr><tr><td>1/8</td><td>0.040</td><td>0.028</td><td>0.376</td><td>0.133</td><td>0.049</td><td>0.77</td></tr></tbody></table>
<figure><img src="image_15.png" /></figure>

<figure><img src="image_16.png" /><figcaption>Figure 12 - Comparison of $u_{\epsilon}$ with $u$ for $\lambda = 10$.</figcaption></figure>
<figure><img src="image_17.png" /><figcaption>Figure 13 - Comparison of $u_{\epsilon}$ with $u$ for $\lambda = 1/18$.</figcaption></figure>
<figure><img src="image_18.png" /></figure>

<figure><img src="image_19.png" /><figcaption>Figure 14 - Comparison of $u_{\epsilon}$ with $u$ for $\lambda = 1/114$.</figcaption></figure>
We remark in table 3 that $||u_\epsilon-u||_{H_0^1(\Omega)}$ does not converge but $||u_\epsilon-u||_{L^2(\Omega)}$ and $||u_\epsilon-u||_{L^\infty(\Omega)}$ converge as we can see in fig. 12, 13, 14. We also remark that the values of $\sup |u_\epsilon-u|$ have a linear decrease as $\epsilon \to 0$, though the coefficients $a_{ij}$ be discontinuous.

### 2.3. An example with a variational inequality

Results of convergence of $u_\epsilon$ towards $u$ have been obtained for simple variational inequalities (see LIONS [1]).

We present here numerical results for the following problem. We introduce the convex set $K$ defined by

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

and we consider $u_\epsilon$ solution of the variational inequality

$$a_\epsilon(u_\epsilon, v-u_\epsilon) \ge (f, v-u_\epsilon) \quad \forall v \in K,$$

where $a_\epsilon$ is the bilinear form associated to $A^\epsilon$ defined in (5). Then the homogenized problem is: find $u \in K$ solution of

$$a(u, v-u) \ge (f, v-u) \quad \forall v \in K,$$

where $a$ is the bilinear form with constant coefficients associated to the homogenized operator $\mathcal{L}$ defined in section 1.1, and we have (see LIONS [1])

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

For the computations we consider the unit cell $Y$ as in fig. 5, with

$$a_{ijo} = \delta_{ij}, \quad a_{ij1} = 5 \delta_{ij}.$$

We take

$$\Omega = ]0,1[ \times ]0,1[$$

$$f = 360 x_1 x_2 (1-x_1) (1-x_2) - 5x_1 - 10.$$

$$\epsilon = 1/2, 1/4, 1/8.$$

Figure 15 shows sections of the functions $u_\epsilon$ and $u$ as indicated in fig. 11. One can see, in this simple case, the good approximation of $u_\epsilon$ by $u$ for $\epsilon = 1/4$ or $\epsilon = 1/8$. One can notice that the free boundary $S_\epsilon$ of the true solutions are very close to the free boundary of the homogenized solution.

On the numerical techniques for solving variational inequalities see GLOWINSKI-LIONS
<figure><img src="image_20.png" /><figcaption>**Figure 15** - Comparison of $u_{\epsilon}$ with $u$ in the case of a variational inequality</figcaption></figure>
TREMOLIERES [1].

### 3. - APPROXIMATION OF THE PERIODIC BEHAVIOUR OF THE SOLUTIONS $u_{\epsilon}$ . CORRECTORS.

One can notice in fig. 12, 13, 14 that solutions $u_{\epsilon}$ have important periodic variations. Thus for $\lambda < 1$ $u_{\epsilon}$ has "peaks" and for $\lambda > 1$ $u_{\epsilon}$ has "flat regions". In this section we look for an approximation of this periodic behaviour with terms of order 1 or 2 in $\epsilon$, provided with the asymptotic expansion.

#### 3.1. Asymptotic expansion - Error estimations.

We consider here the Dirichlet problem (4) introduced in section 1. Let us briefly recall some results presented in BOURGAT-DERVIEUX [1]. The notations are those of section 1. We look for $u_{\epsilon}$ in the form

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

and we have shown in section 1 that

$$(42) \qquad w_0(x,y) = u(x)$$

$$(43) \qquad w_1(x,y) = - \chi^i(y) \frac{\partial u}{\partial x_i}(x) + \tilde{w}_1(x).$$

If we continue the computations we obtain

$$(44) \qquad w_2(x,y) = - \chi^i(y) \frac{\partial \tilde{w}_1}{\partial x_i}(x) + n_{ij}(y) \frac{\partial^2 u}{\partial x_i \partial x_j}(x) + \tilde{w}_2(x),$$

with $n_{ij} \in W$ solution of the equation

$$(45) \qquad a_Y(n_{ij}, \psi) = \int_Y g_{ij}(y) \psi(y) dy, \quad \forall \psi \in W, \quad \mathfrak{M}(n_{ij}) = 0,$$

where $g_{ij}$ is defined by

$$(46) \qquad g_{ij}(y) = - \frac{\partial}{\partial y_k} (a_{ki} \chi^j) - a_{ik} \frac{\partial \chi^j}{\partial y_k} + a_{ij} - q_{ij}.$$

Let us introduce $\theta_{\epsilon} \in H^1(\Omega)$ solution of the equation

$$(47) \qquad \begin{cases} A^{\epsilon} \theta_{\epsilon} = 0 \\ \theta_{\epsilon} = \chi^{i} \left(\frac{x}{\epsilon}\right) \frac{\partial u}{\partial x_{i}}(x) \text{ on } \Gamma. \end{cases}$$
If $a_{ij} \in C^{1+\alpha}(\bar{\Omega})$, $\alpha > 0$, $f \in C^{5+\alpha}(\bar{\Omega})$ and the boundary $\Gamma$ is of class $C^5$, then one has the following estimations.

(48) $$ \| \theta_{\varepsilon} \|_{C^0(\bar{\Omega})} < K $$

(49) $$ \| u_{\varepsilon} - u - \varepsilon \left(-\chi^{i}\left(\frac{X}{\varepsilon}\right) \frac{\partial u}{\partial x_{i}}(x) + \tilde{w}_{1}(x) + \theta_{\varepsilon}\right) \|_{C^{0}(\bar{\Omega})} < K \varepsilon^{2} $$

(50) $$ \| u_{\varepsilon} - u - \varepsilon \left(-\chi^{i}\left(\frac{X}{\varepsilon}\right) \frac{\partial u}{\partial x_{i}}(x) + \theta_{\varepsilon}\right) \|_{H^{1}(\Omega)} \leq C \varepsilon $$

One can also find results of this type in BENSOUSSAN-LIONS-PAPANICOLAOU [2], BABUSKA [2].

### Symmetric context.

We consider here the particular case where the coefficients $a_{ij}$ have some properties of symmetry (see section 1.3).

One has the following result (see BOURGAT-DERVIEUX [1]),

**PROPOSITION 2** : If the coefficients $a_{ij}$ satisfy the symmetry properties (31) then for the Dirichlet problem (4) one has if $a_{ij}=a_{ji}$

(51) $$ \tilde{w}_1 = 0 $$

If, in addition, for $\varepsilon$ given, $\Omega$ is the reunion of entire cells one has

(52) $$ \theta_{\varepsilon} = 0 $$

**Remark 1** : Assumption $\mathcal{M}(\chi^i) = \frac{1}{|Y|} \int_Y \chi^i(y) dy = 0$ is used in the proof of (51). This justifies the choice made in section 1.1 to determine $\chi^i$.

If $u \in C^{6+\alpha}(\bar{\Omega})$, then using results of proposition 2, we obtain in a symmetric context the following estimations.

(53) $$ \| u_{\varepsilon} - u + \varepsilon \chi^{i} \left(\frac{X}{\varepsilon}\right) \frac{\partial u}{\partial x_{i}}(x) \|_{C^{0}(\bar{\Omega})} \leq K \varepsilon^{2} $$

(54) $$ \| u_{\varepsilon} - u + \varepsilon \chi^{i} \left(\frac{X}{\varepsilon}\right) \frac{\partial u}{\partial x_{i}}(x) \|_{H^{1}(\Omega)} \leq C \varepsilon $$

(55) $$ \| u_{\varepsilon} - u + \varepsilon \chi^{i} \left(\frac{X}{\varepsilon}\right) \frac{\partial u}{\partial x_{i}}(x) - \varepsilon^{2} n_{ij} \left(\frac{X}{\varepsilon}\right) \frac{\partial^{2} u}{\partial x_{i} \partial x_{j}}(x) \|_{H^{1}(\Omega)} \leq K \varepsilon^{3/2}. $$

One can see that in a symmetric context, only the periodic terms of order 1 and 2 in $\varepsilon$ give a better approximation of the solution $u_{\varepsilon}$. They are respectively called "first order periodic corrector" and "second order periodic corrector".
**Remark 2** : We have shown in BOURGAT-DERVIEUX that it is the second order periodic corrector $\eta_{ij}(y) \frac{\partial^2 u}{\partial x_i \partial x_j}$, which gives the correction of the peaks of $u_\epsilon$ in the case $\lambda < 1$. J.L. LIONS, using a "stiff expansion" method, obtains in LIONS [2] a similar correction.

### 3.2. Numerical techniques for the computation of the correctors.

Let us denote by

$u_1$ the corrected solution obtained with the first order periodic corrector

$$ (56) \qquad u_1(x) = u(x) - \epsilon \chi^i \left(\frac{x}{\epsilon}\right) \frac{\partial u}{\partial x_i}(x), $$

$u_2$ the corrected solution obtained with the first and second order periodic correctors,

$$ (57) \qquad u_2(x) = u(x) - \epsilon \chi^i \left(\frac{x}{\epsilon}\right) \frac{\partial u}{\partial x_i}(x) + \epsilon^2 \eta_{ij} \left(\frac{x}{\epsilon}\right) \frac{\partial^2 u}{\partial x_i \partial x_j}(x). $$

We notice that equations (45) and (13) giving respectively $\eta_{ij}$ and $\chi^i$ differ only in their right hand side. So we compute $\eta_{ij}$ as $\chi^i$. The only difficulty is the discontinuity of coefficients $a_{ki}$ in the term

$$ \int_Y - \frac{\partial}{\partial y_k} (a_{ki} \chi^j) \psi \, dy $$

which appears in the right hand side of equation (45). Using the Green's formula we obtain :

$$ - \int_Y \frac{\partial}{\partial y_k} (a_{ki} \chi^j) \psi \, dy = \int_Y a_{ki} \chi^j \frac{\partial \psi}{\partial y_k} \, dy - \int_{\partial Y} a_{ki} \chi^j \psi \cos(\vec{n}, \vec{y}_k) \, d\gamma. $$

But $a_{ki}, \chi^j, \psi$ are Y-periodic functions and $\cos(\vec{n}, \vec{y}_k)$ takes opposite values on the opposite sides of Y so

$$ - \int_Y \frac{\partial}{\partial y_k} (a_{ki} \chi^j) \psi \, dy = \int_Y a_{ki} \chi^j \frac{\partial \psi}{\partial y_k} \, dy. $$

For computing the derivatives $\frac{\partial u}{\partial x_i}$ and $\frac{\partial^2 u}{\partial x_i \partial x_j}$, which appear in the calculus of corrected solutions $u_1$ and $u_2$, a fine discretization of $\Omega$ is needed. We use here 1152 triangles. If the triangles are not locally too different and if the right hand side f is regular we can approximate $\frac{\partial u}{\partial x_i}$ at each node M by

$$ (58) \qquad \frac{\partial u}{\partial x_i}(M) = \frac{1}{N_v(M)} \sum_{T \in V(M)} \left[ \frac{\partial u}{\partial x_i} \right]_T, $$

where
v(M) is the set of triangles of which M is a vertex,
$N_v(M)$ is the number of triangles of v(M),
$[\frac{\partial u}{\partial x_i}]_T$ is the value of the derivative $\frac{\partial u}{\partial x_i}$ on the triangle T.
Using formula (58) twice one case obtain an approximation of $\frac{\partial^2 u}{\partial x_i \partial x_j}$

3.3. First order periodic corrector

Let us consider the case $\lambda > 1$. We take $\lambda = 10$. Figure 16 shows sections of $u_\epsilon$, $u$, $u_1$ along a diagonal of $\Omega$ for $\epsilon = 1/2$.

<figure><img src="image_21.png" /><figcaption>**Figure 16** - Diagonal sections of $u_\epsilon$, $u$, $u_1$ for $\lambda=10$, $\epsilon=1/2$.</figcaption></figure>
One can see in fig. 16 that the first order periodic corrector gives a good approxi-
mation of the slopes of $u_{\epsilon}$ (i.e. of its gradient).

Let us consider now the case $\lambda < 1$. We take $\lambda = 1/114$.

Figure 17 shows diagonal sections of $u_{\epsilon}$, $u$, $u_1$ for $\epsilon = 1/4$.

<figure><img src="image_22.png" /><figcaption>Figure 17 - Diagonal sections of $u_{\epsilon}$, $u$, $u_1$ for $\lambda = 1/114$, $\epsilon = 1/4$.</figcaption></figure>

One can see in fig. 17 that the first order periodic corrector does not give any
correction of the peaks. More precisely we have indicated in table 4 the error $u_{\epsilon}-u_1$
in the norms of spaces $L^{\infty}(\Omega)$, $L^2(\Omega)$, $H^1(\Omega)$, (compare to table 3).
<table><caption>Table 4 - Convergence of the error $u_{\epsilon}-u_{1}$ for $\lambda = 10, \lambda = 1/114$.</caption><thead><tr><th rowspan="2">ε</th><th colspan="3">λ = 10</th><th colspan="3">λ = 1/114</th></tr><tr><th>Sup|u<sub>ε</sub>-u<sub>1</sub>|</th><th>||u<sub>ε</sub>-u<sub>1</sub>||<sub>L</sub><sup>2</sup></th><th>||u<sub>ε</sub>-u<sub>1</sub>||<sub>H</sub><sup>1</sup></th><th>Sup|u<sub>ε</sub>-u<sub>1</sub>|</th><th>||u<sub>ε</sub>-u<sub>1</sub>||<sub>L</sub><sup>2</sup></th><th>||u<sub>ε</sub>-u<sub>1</sub>||<sub>H</sub><sup>1</sup></th></tr></thead><tbody><tr><td>1/2</td><td>0.020</td><td>0.032</td><td>0.085</td><td>2.24</td><td>0.56</td><td>0.952</td></tr><tr><td>1/4</td><td>0.0097</td><td>0.012</td><td>0.050</td><td>0.554</td><td>0.183</td><td>0.889</td></tr><tr><td>1/8</td><td>0.0033</td><td>0.004</td><td>0.033</td><td>0.133</td><td>0.042</td><td>0.707</td></tr></tbody></table>

### 3.4. Second order periodic corrector

Figure 18 shows, for $\lambda = 1/114$ and $\epsilon = 1/4$, diagonal sections of $u_{\epsilon}$, $u$, $u_{2}$. If we compare figures 17 and 18 we verify remark 3.2, that is to say that the second order periodic corrector gives the complete correction of the peaks.

We indicate in table 5 the error $u_{\epsilon}-u_{2}$ and one can compare these results to those of table 4.

<table><caption>Table 5 - Convergence of the error $u_{\epsilon}-u_{2}$ for $\lambda = 1/114$.</caption><thead><tr><th>ε</th><th>Sup|u-u<sub>2</sub>|</th><th>||u<sub>ε</sub>-u<sub>2</sub>||<sub>L</sub><sup>2</sup></th><th>||u<sub>ε</sub>-u<sub>2</sub>||<sub>H</sub><sup>1</sup></th></tr></thead><tbody><tr><td>1/2</td><td>0.181</td><td>0.031</td><td>0.181</td></tr><tr><td>1/4</td><td>0.011</td><td>0.012</td><td>0.012</td></tr><tr><td>1/8</td><td>0.013</td><td>0.012</td><td>0.034</td></tr></tbody></table>

We notice in table 5 that the results obtained for $\epsilon = 1/8$ are worst than those obtained for $\epsilon = 1/4$. This is due to the error made in the computation of $u_{\epsilon}$ (see section 2.1). Therefore, if $\lambda < 1$ and $\epsilon$ small, even if the direct computation of $u_{\epsilon}$ is possible, homogenization gives better results in a symmetric context.
<figure><img src="image_23.png" /><figcaption>Figure 18 - Diagonal sections of $u_\epsilon$, $u$, $u_2$ for $\lambda = 1/114$, $\epsilon = 1/4$.</figcaption></figure>

## CONCLUSION

It follows from the previous results that the homogenization method is an effective method for solving boundary value problems with periodic coefficients as soon as the number of cells is large. The homogenized solution can be computed easily by solving equations in a representative cell and then by solving an equation with constant coefficients on the whole domain. Furthermore the use of periodic correctors gives, when the coefficients $a_{ij}$ have some properties of symmetry, a good approximation of the periodic behaviour of the solution and this is useful even if the number of cells is not very small.
# ACKNOWLEDGEMENTS

I am thankful to A. DERVIEUX and H. LANCHON for their help at various stages of this report.

# REFERENCES

## I. BABUSKA

[1] "Solution of problems with interfaces and singularities", Tech. Note BN-789 (April 1974), Univ. of Maryland.

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

## A. BENSOUSSAN-J.L. LIONS-G. PAPANICOLAUOU

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

[2] *Asymptotic analysis of partial derivative equations with highly oscillating coefficients*. North Holland- 1978 - Amsterdam.

## J.F. BOURGAT-A. DERVIEUX

[1] "Méthode d'homogenéisation des opérateurs à coefficients périodiques. Etude des correcteurs provenant du développement asymptotique". IRIA-LABORIA, Rapport de recherche à paraître.

## J.F. BOURGAT-H. LANCHON

[1] "Application of the homogenization method to composite materials with periodic structure". IRIA-LABORIA, Rapport de Recherche N° 208, Dec. 76.

## E. DE GIORGI-S. SPAGNOLO

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

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

[1] *Analysis of structural composite materials*. Marcel Dekker Inc., New-York, 1973.

## R. GLOWINSKI-J.L. LIONS-R. TREMOLIERES

[1] *Analyse Numérique des Inéquations variationnelles*, Paris, Dunod, (1976)

## J.L. LIONS

[1] Cours au Collège de France 1975-1976

[2] Remarques sur les aspects numériques de la méthode d'homogenéisation dans les matériaux composites. Colloque IRIA-Novosibirsk, Juin 1976, Dunod (à paraître)

## E. SANCHEZ PALENCIA

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

## L. TARTAR

[1] Cours Peccot, Collège de France, 1977.