MODELISATION DE LA CONDUCTION
THERMIQUE DANS UN MILIEU HETEROGENE -
NOTION DE CONDUCTIVITE THERMIQUE
EFFECTIVE

Franck Enguehard

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

To cite this version:

Franck Enguehard. MODELISATION DE LA CONDUCTION THERMIQUE DANS UN MILIEU HETEROGENE - NOTION DE CONDUCTIVITE THERMIQUE EFFECTIVE. CEA/DAM, Le Ripault. 2002. hal-05307817

HAL Id: hal-05307817

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

Submitted on 10 Oct 2025

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
MODELISATION DE LA CONDUCTION
THERMIQUE DANS UN MILIEU HETEROGENE
NOTION DE CONDUCTIVITE THERMIQUE
EFFECTIVE

Auteur :
Franck Enguehard

Rapport d'étude

CEA / Le Ripault

Janvier 2002
# MODELISATION DE LA CONDUCTION THERMIQUE DANS UN MILIEU HETEROGENE NOTION DE CONDUCTIVITE THERMIQUE EFFECTIVE

## RESUME

Dans ce document, nous abordons la modélisation du transfert conductif stationnaire dans les milieux hétérogènes par le biais d'un programme informatique de simulation thermique que nous présentons puis testons sur différentes architectures matérielles.

Après avoir posé dans la partie 1- le problème mathématique de "conductivimétrie thermique numérique" à traiter, nous décrivons dans la partie 2- la méthode de résolution que nous avons retenue : il s'agit d'une simple technique de discrétisation en différences finies, laquelle aboutit à un système linéaire qui est résolu au moyen de l'algorithme de Fletcher et Reeves.

La partie 3- est consacrée dans un premier temps à la validation de notre programme de calcul. Dans un deuxième temps, on soumet le programme à des géométries matérielles hétérogènes assez variées, bidimensionnelles pour la grande majorité, sensées constituer des images simplifiées mais réalistes approchant les microstructures de grandes familles de matériaux comme les composites {fibres + résine}, les isolants thermiques microporeux ou les mélanges binaires {matrice + charges}. L'étude des matériaux composites est abordée sous deux éclairages différents : celui de la périodicité spatiale et celui de l'anisotropie. Quant aux matériaux isolants thermiques microporeux, ceux-ci sont examinés aux deux échelles caractéristiques de leur microstructure, c'est-à-dire les échelles nano- et micrométrique. Enfin, les mélanges binaires sont évoqués en géométries bi- et tridimensionnelle et en situation de contraste thermique particulièrement fort.

Ce document ouvre de nombreux chantiers intéressant une grande variété de matériaux ; parmi ceux-ci, deux catégories particulières de grand intérêt pour nos applications : les matériaux composites et les isolants thermiques microporeux. Parallèlement, une multitude de difficultés nouvelles se font jour au fil des échantillons testés, difficultés qu'il faudra progressivement résoudre pour améliorer notre connaissance des phénomènes de transfert dans les matériaux hétérogènes.
détermination du champ de température

s'établissant au sein d'un milieu périodique...................page 38

3-1-5- Milieu à deux constituants et anisotrope...................................................page 39

3-1-5-1- Calcul par différences finies de la conductivité thermique effective Λ en fonction de l'angle α d'orientation du milieu hétérogène...................................................................................page 39

3-1-5-2- Recherche par voie analytique de la conductivité thermique effective Λ du milieu homogène anisotrope équivalent...................................................................................................................page 40

3-1-5-3- Peut-on faire coïncider les résultats des approches numérique et analytique ? ...................................................................................page 45

3-1-6- Quelques milieux particuliers à deux constituants...................................................page 47

3-1-6-1- Matériaux fissurés...................................................................................................................page 47

3-1-6-1-1- Fissure isolée...................................................................................................................page 47

3-1-6-1-2- Distribution de fissures ...................................................................................................................page 49

3-1-6-2- Matériau poreux ...................................................................................................................page 51

3-1-6-3- Matériau en grains et matrice microporeuse .................................................................page 51

3-1-6-3-1- Matériau en grains...................................................................................................................page 52

3-1-6-3-2- Matrice microporeuse ...................................................................................................................page 52

3-1-6-4- Dispersion dans une matrice ...................................................................................................................page 57

3-2- Un exemple en géométrie tridimensionnelle : estimation numérique des conductivités thermiques effectives de deux échantillons élastomères chargés ...................................................................................................................................................................................page 59

4- Conclusion - les suites à donner à ce travail...................................................................................................................page 61

Liste des figures...................................................................................................................................................................................page 64

Bibliographie ...................................................................................................................................................................................page 130
## Avertissement

Ce document contenant un nombre important de figures, celles-ci ont été regroupées à la suite du texte dans une section « liste des figures » commençant à la page 64.
# 1- Introduction - quel problème veut-on résoudre ?

Ce rapport vient en complément de deux documents précédents [1,2]. Ici nous nous limitons à l'étude du transfert thermique par conduction au sein d'un matériau hétérogène dont les différents constituants sont tous thermiquement isotropes ; cette étude aboutit à la détermination de la conductivité thermique effective de ce type de matériau dans une direction donnée.

L'évaluation des propriétés physiques effectives d'un milieu hétérogène est la finalité de la science de l'homogénéisation : il s'agit de partir des caractéristiques physiques (supposées connues) et de l'organisation spatiale (supposée connue) des différents constituants entrant dans la composition du milieu hétérogène et d'en déduire les caractéristiques physiques de son milieu homogène équivalent, c'est-à-dire du milieu homogène dont le comportement se rapproche le plus du sien. En conduction thermique, on va donc supposer connus les conductivités thermiques ainsi que l'arrangement dans l'espace des différents constituants, et on va essayer d'en déduire le tenseur des conductivités thermiques effectives du milieu hétérogène.

Dans le cadre de ce rapport, l'étape d'homogénéisation thermique par changement d'échelle qui vient d'être évoquée va être traitée de manière informatique au moyen d'une technique qu'on pourrait qualifier de technique des plaques gardées numérique en référence à la technique expérimentale classique des plaques gardées pour la mesure de conductivité thermique (voir par exemple la référence [3] pour plus de détails concernant cette technique expérimentale). Dans ses grandes lignes, le principe de conductivimétrie thermique numérique qui va être appliqué est le suivant. Nous allons donner une forme rectangulaire en 2-D, parallélépipédique en 3-D, à notre échantillon hétérogène. Puis nous imposerons sur deux faces opposées de la géométrie deux températures différentes, l'une chaude $T_c$ et l'autre froide $T_f$ (avec $T_c > T_f$), tandis que les autres faces de la géométrie seront supposées adiabatiques (la figure 1 illustre en 2-D l'ensemble des conditions aux limites qui seront imposées à la géométrie). La résolution du champ de température stationnaire $T(M)$ (T désignant la température et M la position) s'établissant au sein de l'échantillon hétérogène nous permettra alors d'évaluer la puissance thermique $\Phi$ traversant l'échantillon de la face chaude à la face froide ; de cette valeur de $\Phi$ découlera celle de la conductivité thermique effective $\Lambda$ de l'échantillon dans la direction normale aux deux faces à température imposée.

Le gros du travail va donc consister à résoudre le problème de conduction thermique stationnaire à travers l'échantillon hétérogène ; ce problème s'écrit sous la forme mathématique suivante :

$$ egin{cases} 
\vec{\nabla} \cdot \vec{\phi} = 0 \quad \text{avec } \vec{\phi} = -\lambda \vec{\nabla} T \quad \text{dans le volume de l'échantillon} \\ 
T = T_c \quad \text{sur la face chaude de l'échantillon} \\ 
T = T_f \quad \text{sur la face froide de l'échantillon} \\ 
\vec{\phi} \cdot \vec{n} = -\lambda \frac{\partial T}{\partial n} = 0 \quad \text{sur les autres faces de l'échantillon} 
\end{cases} 
\quad (1) $$

Dans cette formulation, $\vec{\nabla}$ est l'opérateur nabla ; $\vec{\phi}$ désigne le flux thermique conductif s'exprimant selon la loi de Fourier par $\vec{\phi} = -\lambda \vec{\nabla} T$ ; enfin, dans l'équation $\vec{\phi} \cdot \vec{n} = -\lambda \frac{\partial T}{\partial n} = 0$, le
vecteur $\vec{n}$ et la dérivée partielle par rapport à n réfèrent tous deux à la direction normale à la surface considérée.

Rappelons que la loi de Fourier $\bar{\phi} = -\lambda \vec{\nabla} T$ est une relation locale faisant intervenir la conductivité thermique $\lambda$ au point considéré. Lorsque le milieu d'étude est homogène et thermiquement isotrope, la quantité $\lambda$ ne dépend pas de la position ; l'équation de volume du système (1) s'en trouve grandement simplifiée puisqu'elle devient la célèbre équation de Poisson :

$$ \Delta T = 0 \tag{2} $$

dans laquelle $\Delta$ désigne l'opérateur laplacien. Dans ces circonstances, le système (1) est particulièrement simple à résoudre : le champ de température $T(M)$ est fonction de la seule variable $x$ définie comme étant la distance du point $M$ à la face chaude (voir les notations sur la figure 1), et ce champ $T(x)$ est une fonction affine de $x$ :

$$ T(x) = T_c - (T_c - T_f) \frac{x}{X} \tag{3} $$

$X$ désignant la distance entre la face chaude et la face froide.

Malheureusement, le type de milieu que nous étudions ici étant hétérogène, il s'ensuit que la conductivité thermique $\lambda$ n'est pas constante d'un point à un autre de l'échantillon et doit donc être considérée, au même titre que la température $T$, comme un champ. Dès lors, on ne peut plus extraire de manière simple le terme de conductivité thermique $\lambda$ de l'équation de volume $\vec{\nabla} \cdot (-\lambda \vec{\nabla} T) = 0$, et il ne reste plus que des méthodes numériques pour résoudre le problème (1).

Une fois cette difficulté contournée et le problème résolu, il restera à calculer la puissance thermique $\Phi$ et la conductivité thermique effective $\Lambda$. Il est aisé de constater que, les faces de l'échantillon autres que les deux faces à température imposée étant adiabatiques, le flux thermique conductif $\vec{\phi}$ se trouve canalisé dans l'échantillon, de sorte que la puissance thermique $\Phi$ qui s'exprime par :

$$ \Phi = \iint_S \vec{\phi} \cdot \vec{n} \, dS \tag{4} $$

ne dépend pas de la surface d'intégration $S$ pourvu que cette surface constitue une section de ce tuyau de flux thermique limité par les faces adiabatiques. [ Cette indépendance de $\Phi$ vis-à-vis de la section $S$ s'établit de la façon suivante : de l'équation de volume $\vec{\nabla} \cdot \vec{\phi} = 0$ on tire, en utilisant le théorème d'Ostrogradski, l'égalité :

$$ \iint_{\Sigma} \vec{\phi} \cdot \vec{n} \, dS = 0 \tag{5} $$
pour toute surface fermée Σ entièrement incluse dans l'échantillon. Prenant pour Σ une portion du tuyau de flux thermique limitée par deux sections distinctes S₁ et S₂, l'égalité (5) devient :

$$ \iint_{S_1} \bar{\phi} \cdot \bar{n} \, dS = \iint_{S_2} \bar{\phi} \cdot \bar{n} \, dS \qquad (6) $$

d'où le résultat recherché. On pourra donc par exemple calculer la puissance thermique Φ par intégration sur n'importe quelle section plane de l'échantillon S parallèle aux deux faces à température imposée par la formule :

$$ \Phi = \iint_S \bar{\phi} \cdot \bar{x} \, dS = \iint_S \phi_x \, dS \qquad (7) $$

x étant le vecteur unitaire de la figure 1 et se trouvant être par la même occasion le vecteur normal à la section S considérée. Mais si le volume occupé par notre milieu hétérogène d'étude est maintenant rempli par son milieu homogène effectif de conductivité thermique effective Λ dans la direction x, alors le principe général de l'homogénéisation thermique selon lequel une sollicitation thermique d'amplitude (T_c - T_f) donnée se traduit par deux puissances thermiques Φ identiques dans les deux milieux amène une nouvelle expression de Φ vue cette fois du milieu homogène effectif :

$$ \Phi = C(T_c - T_f) = \frac{\Lambda S}{X}(T_c - T_f) \qquad (8) $$

C désignant la conductance thermique dans la direction x du milieu homogène effectif. La combinaison des relations (7) et (8) aboutit à l'expression suivante de la conductivité thermique effective Λ :

$$ \Lambda = \frac{X}{S(T_c - T_f)} \iint_S \phi_x \, dS = \frac{X}{S(T_c - T_f)} \iint_S (-\lambda \frac{\partial T}{\partial x}) \, dS \qquad (9) $$

et c'est cette expression qui nous permettra, à partir de la connaissance du champ de température stationnaire T(M) s'établissant au sein de l'échantillon hétérogène, d'évaluer la conductivité thermique effective Λ.

Ce rapport sera structuré en deux grandes parties. Dans un premier temps (chapitre 2-), on décrira en détail la technique utilisée pour le calcul du champ de température stationnaire T(M), de la puissance thermique Φ et de la conductivité thermique effective Λ. Ensuite (chapitre 3-), on exposera plusieurs résultats de calcul illustrant les performances et limitations de l'outil numérique développé.
## 2- Description de la technique numérique de détermination du champ de température et de la conductivité thermique effective

Avant de commencer cette description, nous ferons deux remarques en guise de préambule :

1. Par souci de simplification, le milieu hétérogène d'étude sera considéré bidimensionnel. Cela dit, tous les développements 2-D qui vont suivre peuvent très facilement être étendus à un espace tridimensionnel.

2. Comme il a été dit plus haut, la chose à faire en premier lieu est de résoudre le système (1). Si dans ce système on remplace la variable température T par une nouvelle variable θ sans dimension, qu'on pourrait qualifier de température réduite et définie par :

$$ \theta = \frac{T - T_f}{T_c - T_f} \qquad (10) $$

alors le système (1) devient le système (11) suivant :

$$ \left\{ \begin{array}{l} \vec{\nabla} \cdot (-\lambda \vec{\nabla} \theta) = 0 \quad \text{dans le volume de l'échantillon} \\[1em] \theta = 1 \quad \text{sur la face chaude de l'échantillon} \\[1em] \theta = 0 \quad \text{sur la face froide de l'échantillon} \\[1em] -\lambda \frac{\partial \theta}{\partial n} = 0 \quad \text{sur les autres faces de l'échantillon} \end{array} \right. \qquad (11) $$

Le problème mathématique à résoudre est donc absolument le même, à ceci près que les températures imposées aux deux faces opposées de l'échantillon sont maintenant 1 sur la face chaude et 0 sur la face froide. Par ailleurs, les valeurs de la puissance thermique Φ et de la conductivité thermique effective Λ peuvent tout à fait se déduire de la connaissance du champ des températures réduites θ(M) ; en effet, les expressions (7) et (9) deviennent respectivement :

$$ \Phi = \iint_S \left(-\lambda \frac{\partial T}{\partial x}\right) dS = (T_c - T_f) \iint_S \left(-\lambda \frac{\partial \theta}{\partial x}\right) dS \qquad (12) $$

$$ \Lambda = \frac{X}{S} \iint_S \left(-\lambda \frac{\partial \theta}{\partial x}\right) dS \qquad (13) $$

Raisonner sur les températures réduites θ à la place des températures réelles T présente l'avantage suivant : le nouveau problème de conduction thermique (11) incluant des températures imposées valant 0 et 1, il s'ensuit que le champ des températures réduites θ(M) est forcément compris entre 0 et 1 en tout point de l'échantillon hétérogène. Comme
on le verra par la suite, cette propriété du champ $\theta(M)$ sera bien utile pour l'étude de la
convergence de processus itératifs de résolution numérique. D'autre part, l'expression (13)
de la conductivité thermique effective $\Lambda$ en fonction des températures réduites $\theta$ a ceci de
rassurant qu'elle ne contient aucune des deux températures $T_c$ et $T_f$ ; le champ $\theta(M)$ étant
indépendant de ces deux quantités (puisqu'elles sont absentes du système (11)), il apparaît
clairement dans l'expression (13) que la valeur de la conductivité thermique effective $\Lambda$ ne
dépend pas des valeurs des deux températures $T_c$ et $T_f$ imposées aux extrémités de
l'échantillon hétérogène. Pour ces deux raisons, nous nous attach erons dans ce qui suit à la
résolution du système (11), puis du champ $\theta(M)$ obtenu nous déduirons la puissance
thermique $\Phi$ et la conductivité thermique effective $\Lambda$ par les relations (12) et (13).

2-1- **_Discrétisation du problème en différences finies_**

La géométrie considérée est un rectangle dont les côtés X et Y sont découpés respectivement en I et en J petits segments d'égale longueur $\delta = X/I = Y/J$. Le rectangle se trouve ainsi fractionné en $I \times J$ petits carrés de côté $\delta$, et on suppose que ce fractionnement est suffisamment fin pour que la matière occupant chacun des petits carrés puisse être considérée comme homogène et isotrope.

Nos inconnues sont les températures $\theta_{ij}$, $1 \le i \le I$, $1 \le j \le J$, aux centres des $I \times J$ petits carrés du maillage. Certaines de ces $I \times J$ inconnues sont fixées par les deux conditions aux limites de température imposée : il s'agit en l'occurrence des températures $\theta_{1,j}$, $1 \le j \le J$, et des températures $\theta_{I,j}$, $1 \le i \le I-1$, $1 \le j \le J$. Quatre cas sont à étudier suivant la localisation du carré auquel l'on s'intéresse :

1. Dans un premier temps, nous considérons un carré (centré au noeud N) se trouvant dans le "volume" du maillage, c'est-à-dire remplissant les deux conditions suivantes :

1. Il a des voisins à sa gauche (carré centré au noeud G), à sa droite (carré centré au noeud D), à son haut (carré centré au noeud H) et à son bas (carré centré au noeud B),

2. Les températures $\theta_G$, $\theta_D$, $\theta_H$ et $\theta_B$ aux centres des quatre carrés voisins sont tout aussi inconnues que la température $\theta_N$ du carré d'intérêt.

La figure 2 résume les notations adoptées et montre l'ensemble des carrés concernés.
L'équation de volume $\vec{\nabla} \cdot \vec{\psi} = 0$ avec $\vec{\psi} = -\lambda \vec{\nabla} \theta$ intégrée sur la frontière du carré donne
$$ \sum \Psi_{\text{sortant}} = 0 $$ : la somme des puissances thermiques sortant du carré est égale à zéro. Le
bilan énergétique $$ \sum \Psi_{\text{sortant}} = 0 $$ peut encore s'écrire, avec les notations évidentes :

$$ \Psi_{N \to G} + \Psi_{N \to D} + \Psi_{N \to H} + \Psi_{N \to B} = 0 \quad (14) $$
La puissance thermique $\Psi_{N \to G}$ sortant du carré dans la direction $N \to G$ s'exprime simplement par $\Psi_{N \to G} = C_{NG}(\theta_N - \theta_G)$, $C_{NG}$ désignant la conductance thermique entre les noeuds N et G. Or, les carrés centrés en N et en G étant occupés de matières homogènes et isotropes de conductivités thermiques respectives $\lambda_N$ et $\lambda_G$, le conducteur thermique entre les noeuds N et G apparaît dès lors comme étant constitué de deux résistances thermiques en série, toutes deux de section $\delta$ et de longueur $\delta/2$, la première étant décrite par une conductivité $\lambda_N$ et la seconde par une conductivité $\lambda_G$ (voir figure 3). La conductivité thermique équivalente $\lambda_{NG}$ du conducteur thermique entre les noeuds N et G s'obtient donc par l'égalité :

$$ \frac{1}{\lambda_{NG}} \frac{\delta (\text{longueur})}{\delta (\text{section})} = \frac{1}{\lambda_N} \frac{\delta / 2 (\text{longueur})}{\delta (\text{section})} + \frac{1}{\lambda_G} \frac{\delta / 2 (\text{longueur})}{\delta (\text{section})} \quad (15) $$

soit :

$$ \frac{1}{\lambda_{NG}} = \frac{\frac{1}{\lambda_N} + \frac{1}{\lambda_G}}{2} \quad (16) $$

et il apparaît ainsi que la conductivité $\lambda_{NG}$ est la moyenne harmonique des deux conductivités $\lambda_N$ et $\lambda_G$. La conductance thermique $C_{NG}$ et la puissance thermique $\Psi_{N \to G}$ ont alors pour expressions :

$$ C_{NG} = \lambda_{NG} \frac{\delta (\text{section})}{\delta (\text{longueur})} = \lambda_{NG} \quad (17) $$

$$ \Psi_{N \to G} = C_{NG}(\theta_N - \theta_G) = \lambda_{NG}(\theta_N - \theta_G) = \frac{2\lambda_N \lambda_G}{\lambda_N + \lambda_G} (\theta_N - \theta_G) \quad (18) $$

et l'équation de volume (14) devient :

$$ \left( \frac{2\lambda_N \lambda_G}{\lambda_N + \lambda_G} + \frac{2\lambda_N \lambda_D}{\lambda_N + \lambda_D} + \frac{2\lambda_N \lambda_H}{\lambda_N + \lambda_H} + \frac{2\lambda_N \lambda_B}{\lambda_N + \lambda_B} \right) \theta_N \dots \\ \dots - \frac{2\lambda_N \lambda_G}{\lambda_N + \lambda_G} \theta_G - \frac{2\lambda_N \lambda_D}{\lambda_N + \lambda_D} \theta_D - \frac{2\lambda_N \lambda_H}{\lambda_N + \lambda_H} \theta_H - \frac{2\lambda_N \lambda_B}{\lambda_N + \lambda_B} \theta_B = 0 \quad (19) $$

L'équation ci-dessus peut être simplifiée par 2 puis réécrite à l'aide des indexations i et j introduites précédemment ; ainsi on obtient la forme discrétisée de l'équation de volume du problème (11) à résoudre :
$$
\begin{equation}
\begin{aligned}
& \left( \frac{\lambda_{i,j} \lambda_{i,j-1}}{\lambda_{i,j} + \lambda_{i,j-1}} + \frac{\lambda_{i,j} \lambda_{i,j+1}}{\lambda_{i,j} + \lambda_{i,j+1}} + \frac{\lambda_{i,j} \lambda_{i-1,j}}{\lambda_{i,j} + \lambda_{i-1,j}} + \frac{\lambda_{i,j} \lambda_{i+1,j}}{\lambda_{i,j} + \lambda_{i+1,j}} \right) \theta_{i,j} \dots \\
& \qquad - \frac{\lambda_{i,j} \lambda_{i,j-1}}{\lambda_{i,j} + \lambda_{i,j-1}} \theta_{i,j-1} - \frac{\lambda_{i,j} \lambda_{i,j+1}}{\lambda_{i,j} + \lambda_{i,j+1}} \theta_{i,j+1} \dots \\
& \qquad - \frac{\lambda_{i,j} \lambda_{i-1,j}}{\lambda_{i,j} + \lambda_{i-1,j}} \theta_{i-1,j} - \frac{\lambda_{i,j} \lambda_{i+1,j}}{\lambda_{i,j} + \lambda_{i+1,j}} \theta_{i+1,j} = 0, \\
& \forall (i,j) \in \{3 \dots I-2\} \times \{2 \dots J-1\}
\end{aligned}
\tag{20}
\end{equation}
$$

2. Nous considérons maintenant un carré se trouvant au voisinage des zones influencées par les deux conditions aux limites de température imposée ; la figure 4 montre l'ensemble des carrés concernés.

Lorsque le carré se trouve au voisinage de la condition "chaude", c'est-à-dire en $i = 2$, l'équation de volume discrétisée (20) ci-dessus subit une modification liée au fait que le carré centré au noeud H est situé dans la zone à température imposée $\theta = 1$ (voir figure 4). Par conséquent, cette équation de volume devient :

$$
\begin{equation}
\begin{aligned}
& \left( \frac{\lambda_{2,j} \lambda_{2,j-1}}{\lambda_{2,j} + \lambda_{2,j-1}} + \frac{\lambda_{2,j} \lambda_{2,j+1}}{\lambda_{2,j} + \lambda_{2,j+1}} + \frac{\lambda_{2,j} \lambda_{1,j}}{\lambda_{2,j} + \lambda_{1,j}} + \frac{\lambda_{2,j} \lambda_{3,j}}{\lambda_{2,j} + \lambda_{3,j}} \right) \theta_{2,j} \dots \\
& \qquad - \frac{\lambda_{2,j} \lambda_{2,j-1}}{\lambda_{2,j} + \lambda_{2,j-1}} \theta_{2,j-1} - \frac{\lambda_{2,j} \lambda_{2,j+1}}{\lambda_{2,j} + \lambda_{2,j+1}} \theta_{2,j+1} \dots \\
& \qquad - \frac{\lambda_{2,j} \lambda_{3,j}}{\lambda_{2,j} + \lambda_{3,j}} \theta_{3,j} = \frac{\lambda_{2,j} \lambda_{1,j}}{\lambda_{2,j} + \lambda_{1,j}}, \\
& \forall j \in \{2 \dots J-1\}
\end{aligned}
\tag{21}
\end{equation}
$$

De même, lorsque le carré se trouve au voisinage de la condition "froide", c'est-à-dire en $i = I - 1$, le carré centré au noeud B se situe dans la zone à température imposée $\theta = 0$ (voir figure 4), et par conséquent l'équation de volume discrétisée (20) devient :
$$
\begin{equation}
\left( \frac{\lambda_{I-1,j} \lambda_{I-1,j-1}}{\lambda_{I-1,j} + \lambda_{I-1,j-1}} + \frac{\lambda_{I-1,j} \lambda_{I-1,j+1}}{\lambda_{I-1,j} + \lambda_{I-1,j+1}} + \frac{\lambda_{I-1,j} \lambda_{I-2,j}}{\lambda_{I-1,j} + \lambda_{I-2,j}} + \frac{\lambda_{I-1,j} \lambda_{I,j}}{\lambda_{I-1,j} + \lambda_{I,j}} \right) \theta_{I-1,j} \dots
\end{equation}
\begin{equation*}
\begin{split}
& \dots - \frac{\lambda_{I-1,j} \lambda_{I-1,j-1}}{\lambda_{I-1,j} + \lambda_{I-1,j-1}} \theta_{I-1,j-1} - \frac{\lambda_{I-1,j} \lambda_{I-1,j+1}}{\lambda_{I-1,j} + \lambda_{I-1,j+1}} \theta_{I-1,j+1} \dots \\
& \dots - \frac{\lambda_{I-1,j} \lambda_{I-2,j}}{\lambda_{I-1,j} + \lambda_{I-2,j}} \theta_{I-2,j} = 0,
\end{split}
\end{equation*}
\begin{equation*}
\forall j \in \{2 \dots J-1\}
\end{equation*}
$$

3. Lorsque le carré d'intérêt est soumis à l'une des deux conditions aux limites d'adiabaticité (voir figure 5), l'équation de volume $\vec{\nabla} \cdot \vec{\psi} = 0$ avec $\vec{\psi} = -\lambda \vec{\nabla} \Theta$ doit être reprise et intégrée sur la frontière du carré en tenant compte de l'adiabaticité d'un de ses côtés.

Ainsi, si l'on considère un carré se trouvant à la frontière gauche du domaine d'étude, c'est-à-dire en $j = 1$, l'équation de volume $\vec{\nabla} \cdot \vec{\psi} = 0$ avec $\vec{\psi} = -\lambda \vec{\nabla} \Theta$ donne, après intégration (voir plus haut) :

$$ \Psi_{N \to D} + \Psi_{N \to H} + \Psi_{N \to B} = 0 \quad (23) $$

Reprenant ensuite les expressions (18) pour les puissances thermiques $\Psi_{N \to D}$, $\Psi_{N \to H}$ et $\Psi_{N \to B}$, on aboutit à l'expression discrétisée de l'équation de volume aux noeuds soumis à la condition d'adiabaticité sur la frontière gauche du domaine d'étude :

$$
\begin{equation}
\left( \frac{\lambda_{i,1} \lambda_{i,2}}{\lambda_{i,1} + \lambda_{i,2}} + \frac{\lambda_{i,1} \lambda_{i-1,1}}{\lambda_{i,1} + \lambda_{i-1,1}} + \frac{\lambda_{i,1} \lambda_{i+1,1}}{\lambda_{i,1} + \lambda_{i+1,1}} \right) \theta_{i,1} \dots \\
\qquad \dots - \frac{\lambda_{i,1} \lambda_{i,2}}{\lambda_{i,1} + \lambda_{i,2}} \theta_{i,2} - \frac{\lambda_{i,1} \lambda_{i-1,1}}{\lambda_{i,1} + \lambda_{i-1,1}} \theta_{i-1,1} - \frac{\lambda_{i,1} \lambda_{i+1,1}}{\lambda_{i,1} + \lambda_{i+1,1}} \theta_{i+1,1} = 0, \tag{24}
\end{equation}
\[
\forall i \in \{3 \dots I-2\}
\]
$$

Effectuant le même travail aux carrés situés cette fois à la frontière droite du domaine d'étude, c'est-à-dire en $j = J$, on aboutit à l'équation de volume intégrée :

$$ \Psi_{N \to G} + \Psi_{N \to H} + \Psi_{N \to B} = 0 \quad (25) $$

qui s'écrit encore :
$$
\begin{gather*}
\left( \frac{\lambda_{i,J} \lambda_{i,J-1}}{\lambda_{i,J} + \lambda_{i,J-1}} + \frac{\lambda_{i,J} \lambda_{i-1,J}}{\lambda_{i,J} + \lambda_{i-1,J}} + \frac{\lambda_{i,J} \lambda_{i+1,J}}{\lambda_{i,J} + \lambda_{i+1,J}} \right) \theta_{i,J} \dots \\
\dots - \frac{\lambda_{i,J} \lambda_{i,J-1}}{\lambda_{i,J} + \lambda_{i,J-1}} \theta_{i,J-1} - \frac{\lambda_{i,J} \lambda_{i-1,J}}{\lambda_{i,J} + \lambda_{i-1,J}} \theta_{i-1,J} - \frac{\lambda_{i,J} \lambda_{i+1,J}}{\lambda_{i,J} + \lambda_{i+1,J}} \theta_{i+1,J} = 0, \tag{26} \\
\forall i \in \{3 \dots I-2\}
\end{gather*}
$$

4. Il reste maintenant à traiter le cas des quatre carrés situés aux quatre angles du domaine d'étude (voir figure 6). Pour ces quatre carrés, la discrétisation de l'équation de volume fait intervenir une condition de température imposée et une condition d'adiabaticité.

Prenons l'exemple du carré situé en bas à gauche du domaine, donc correspondant au
couple (i = I - 1, j = 1) (voir figure 6). Le bilan des puissances thermiques sortant de ce
carré s'exprime par :

$$
\Psi_{N \to D} + \Psi_{N \to H} + \Psi_{N \to B} = 0 \tag{27}
$$

Comme $\theta_B = 0$, ce bilan devient :

$$
\begin{equation}
\begin{split}
& \left( \frac{\lambda_{I-1,1} \lambda_{I-1,2}}{\lambda_{I-1,1} + \lambda_{I-1,2}} + \frac{\lambda_{I-1,1} \lambda_{I-2,1}}{\lambda_{I-1,1} + \lambda_{I-2,1}} + \frac{\lambda_{I-1,1} \lambda_{I,1}}{\lambda_{I-1,1} + \lambda_{I,1}} \right) \theta_{I-1,1} \dots \\
& \qquad - \frac{\lambda_{I-1,1} \lambda_{I-1,2}}{\lambda_{I-1,1} + \lambda_{I-1,2}} \theta_{I-1,2} - \frac{\lambda_{I-1,1} \lambda_{I-2,1}}{\lambda_{I-1,1} + \lambda_{I-2,1}} \theta_{I-2,1} = 0
\end{split}
\tag{28}
\end{equation}
$$

Le même travail effectué aux trois autres carrés (i = 2, j = 1), (i = 2, j = J) et
(i = I - 1, j = J) aboutit aux trois équations discrétisées suivantes :

$$
\begin{equation}
\begin{split}
& \left( \frac{\lambda_{2,1} \lambda_{2,2}}{\lambda_{2,1} + \lambda_{2,2}} + \frac{\lambda_{2,1} \lambda_{1,1}}{\lambda_{2,1} + \lambda_{1,1}} + \frac{\lambda_{2,1} \lambda_{3,1}}{\lambda_{2,1} + \lambda_{3,1}} \right) \theta_{2,1} \dots \\
& \quad - \frac{\lambda_{2,1} \lambda_{2,2}}{\lambda_{2,1} + \lambda_{2,2}} \theta_{2,2} - \frac{\lambda_{2,1} \lambda_{3,1}}{\lambda_{2,1} + \lambda_{3,1}} \theta_{3,1} = \frac{\lambda_{2,1} \lambda_{1,1}}{\lambda_{2,1} + \lambda_{1,1}}
\end{split}
\tag{29}
\end{equation}
$$

$$
\begin{equation}
\begin{split}
& \left( \frac{\lambda_{2,J} \lambda_{2,J-1}}{\lambda_{2,J} + \lambda_{2,J-1}} + \frac{\lambda_{2,J} \lambda_{1,J}}{\lambda_{2,J} + \lambda_{1,J}} + \frac{\lambda_{2,J} \lambda_{3,J}}{\lambda_{2,J} + \lambda_{3,J}} \right) \theta_{2,J} \dots \\
& \quad - \frac{\lambda_{2,J} \lambda_{2,J-1}}{\lambda_{2,J} + \lambda_{2,J-1}} \theta_{2,J-1} - \frac{\lambda_{2,J} \lambda_{3,J}}{\lambda_{2,J} + \lambda_{3,J}} \theta_{3,J} = \frac{\lambda_{2,J} \lambda_{1,J}}{\lambda_{2,J} + \lambda_{1,J}}
\end{split}
\tag{30}
\end{equation}
$$
$$
\left( \frac{\lambda_{I-1,J} \lambda_{I-1,J-1}}{\lambda_{I-1,J} + \lambda_{I-1,J-1}} + \frac{\lambda_{I-1,J} \lambda_{I-2,J}}{\lambda_{I-1,J} + \lambda_{I-2,J}} + \frac{\lambda_{I-1,J} \lambda_{I,J}}{\lambda_{I-1,J} + \lambda_{I,J}} \right) \theta_{I-1,J} \cdots \\
\qquad \cdots - \frac{\lambda_{I-1,J} \lambda_{I-1,J-1}}{\lambda_{I-1,J} + \lambda_{I-1,J-1}} \theta_{I-1,J-1} - \frac{\lambda_{I-1,J} \lambda_{I-2,J}}{\lambda_{I-1,J} + \lambda_{I-2,J}} \theta_{I-2,J} = 0 \tag{31}
$$

Ici se termine l'étape de discrétisation du problème (11) que nous avons à résoudre : avec le système constitué des équations (20), (21), (22), (24), (26), (28), (29), (30), (31), nous disposons de la version discrétisée de l'équation de volume pour chaque carré de température inconnue du maillage. Le problème de la détermination des températures inconnues $\theta_{i,j}$, $2 \le i \le I-1$, $1 \le j \le J$, est donc formellement résolu.

2-2- Détermination des températures inconnues

Le système discrétisé {(20),(21),(22),(24),(26),(28),(29),(30),(31)} que nous avons à résoudre est un
système linéaire, conséquence directe de la linéarité par rapport au champ de température θ de
l'équation de volume et des conditions aux limites du problème "continu" (11). L'étape de
détermination des températures inconnues θᵢⱼ, 2 ≤ i ≤ I-1, 1 ≤ j ≤ J, va donc pouvoir être effectuée
par des techniques numériques classiques de résolution de systèmes linéaires.

Le système est réécrit en ordonnant les températures inconnues θᵢⱼ, 2 ≤ i ≤ I-1, 1 ≤ j ≤ J, dans un grand vecteur x̄ : les (I - 2) premières composantes de x̄ sont constituées de la première colonne de la matrice [θᵢⱼ]₂≤i≤I-1,₁≤j≤J, les (I - 2) composantes suivantes de x̄ sont constituées de la deuxième colonne de la matrice [θᵢⱼ]₂≤i≤I-1,₁≤j≤J, et ainsi de suite. Le vecteur x̄ ainsi construit est un vecteur à (I-2)×J composantes.

Le système linéaire {(20),(21),(22),(24),(26),(28),(29),(30),(31)}, lorsqu'exprimé non plus en
fonction de la matrice des inconnues $[\theta_{i,j}]_{2\le i \le l-1, 1 \le j \le s}$ mais en fonction du vecteur des inconnues
$\bar{x} = [x_k]_{1\le k \le (I-2)\times J}$, peut s'écrire sous la forme d'un produit matriciel $A\bar{x} = \bar{b}$, les quantités A et $\bar{b}$
étant respectivement une matrice à $(I-2) \times J$ lignes et $(I-2) \times J$ colonnes et un vecteur à
$(I-2) \times J$ composantes. Les $((I-2) \times J)^2$ composantes de la matrice A et les $(I-2) \times J$
composantes du vecteur $\bar{b}$ se déduisent aisément du système d'équations discrétisées
{(20),(21),(22),(24),(26),(28),(29),(30),(31)}; les allures de ces deux tableaux sont représentées
schématiquement à la figure 7.

Remarquons que les équations (20), (21), (22), (24), (26), (28), (29), (30) et (31) constituant le
système linéaire à résoudre n'ont pas été simplifiées au maximum : ainsi, l'équation (20) aurait pu
encore être simplifiée par $\lambda_{i,j}$, l'équation (21) par $\lambda_{2,j}$, etc ... Ceci n'est pas un oubli : il se trouve
qu'en conservant les différentes équations sous les formes écrites plus haut, la matrice A du système
linéaire $A\bar{x} = \bar{b}$ possède les propriétés remarquables d'être une matrice symétrique (i.e. $^tA = A$, la
notation $^t$ désignant l'opérateur de transposition) définie positive (i.e. $\forall \bar{x} \neq \bar{0}$, $^tA\bar{x} > 0$). Dès lors,

1. La matrice A est inversible, et par conséquent il existe une solution $\bar{x}$ unique au système linéaire $A\bar{x} = \bar{b}$ ;
2. On peut faire appel à toute une panoplie de techniques numériques adaptées à la résolution de systèmes linéaires à matrice symétrique définie positive.

Parmi les différentes techniques numériques que nous avons testées pour la résolution de notre problème, nous avons retenu l'algorithmme de Fletcher et Reeves ; ce choix sera motivé plus bas au chapitre 2-4. La résolution numérique du système linéaire $A\bar{x} = \bar{b}$ par le truchement de cet algorithme aboutit à l'évaluation du vecteur $\bar{x}$ et donc de la matrice des températures inconnues $[\theta_{i,j}]_{2\le i\le1-1, 1\le j\le J}$.

## 2-3- Evaluation de la conductivité thermique effective

Comme il a été expliqué au chapitre 1-, la conductivité thermique effective $\Lambda$, exprimée indifféremment par les relations (9) ou (13), peut être évaluée sur n'importe quelle section plane de l'échantillon S parallèle aux deux faces à température imposée. En géométrie bidimensionnelle, la relation (13) devient, avec les notations de la figure 1 :

$$ \Lambda = \frac{X}{Y} \int_{y=0}^{Y} \psi(x, y) \, dy \qquad (32) $$

et doit fournir une valeur numérique de $\Lambda$ indépendante de la variable x. En version discrétisée, cette nouvelle relation (32) s'écrit :

$$ \Lambda = \frac{2(I-1)}{J} \sum_{j=1}^{J} \frac{\lambda_{i,j} \lambda_{i+1,j}}{\lambda_{i,j} + \lambda_{i+1,j}} (\theta_{i,j} - \theta_{i+1,j}), \quad (33) \\ \forall i \in \{1 \dots I-1\} $$

et c'est cette expression (33) qui est utilisée pour évaluer la conductivité thermique effective $\Lambda$ à partir de la matrice des températures $[\theta_{i,j}]_{2\le i\le1-1, 1\le j\le J}$. En outre, la valeur numérique de $\Lambda$ étant théoriquement indépendante de l'indice i, le calcul des $(I-1)$ quantités exprimées à droite du signe égal dans la relation (33) permet de quantifier l'écart relatif maximal $\varepsilon_{\Lambda}$ entre ces $(I-1)$ quantités et leur valeur moyenne

$$ \bar{\Lambda} = \frac{1}{I-1} \sum_{i=1}^{I-1} \frac{2(I-1)}{J} \sum_{j=1}^{J} \frac{\lambda_{i,j} \lambda_{i+1,j}}{\lambda_{i,j} + \lambda_{i+1,j}} (\theta_{i,j} - \theta_{i+1,j}) \quad (34) $$

et par conséquent d'apprécier la qualité du résultat numérique fourni par l'algorithmme de résolution.
2-4- Quelques informations d'ordre algorithmique

Comme il a été indiqué plus haut, nous faisons appel à l’algorithme de Fletcher et Reeves pour la résolution de notre système linéaire. Deux raisons motivent le choix de cette procédure numérique : son caractère itératif et sa rapidité de convergence. L’algorithme de Fletcher et Reeves se base sur celui, plus général, du gradient conjugué : il s’agit de minimiser la fonction quadratique

$$
f(\bar{y}) = \frac{1}{2} \bar{y}^T A \bar{y} - \bar{b}^T \bar{y} \quad (35)
$$

La matrice A étant symétrique, on démontre aisément que, si x̄ désigne la solution du système linéaire A x̄ = b̄, alors :

$$
f(\bar{y}) - f(\bar{x}) = \frac{1}{2} (\bar{y} - \bar{x}) A (\bar{y} - \bar{x}) \quad (36)
$$

La matrice A étant de plus définie positive, il s'ensuit que f(\bar{y}) > f(\bar{x}) pour tout vecteur \bar{y} ≠ \bar{x} : le vecteur \bar{x} solution du système linéaire A\bar{x} = \bar{b} est donc bien le minimum absolu de la fonction quadratique f. Par ailleurs, le gradient de cette fonction quadratique f valant \nabla f|_{\bar{y}} = A\bar{y} - \bar{b} (car la matrice A est symétrique), rechercher le minimum \bar{x} de la fonction f est bien équivalent à résoudre l'équation A\bar{x} - \bar{b} = \bar{0}. L'algorithme proposé par Fletcher et Reeves pour converger vers le vecteur \bar{x} réalisant le minimum de la fonction quadratique f est le suivant :

1. On part d'un vecteur initial $\bar{x}_0$ [ Il est judicieux de choisir ce vecteur $\bar{x}_0$ le plus proche possible de la solution $\bar{x}$ recherchée. Le $\bar{x}_0$ que nous avons choisi pour initier nos calculs est le vecteur correspondant à un milieu homogène isotrope, c'est-à-dire défini par les températures réduites $\theta_{i,j} = (I-i)/(I-1), i \in \{1...I\}, j \in \{1...J\}$ ] et on calcule les deux vecteurs $\bar{g}_0$ et $\bar{h}_0$ exprimés par $\bar{g}_0 = \bar{h}_0 = A\bar{x}_0 - \bar{b}$.

2. Le passage du triplet de vecteurs ($\bar{x}_n, \bar{g}_n, \bar{h}_n$) à l'itération n au triplet de vecteurs ($\bar{x}_{n+1}, \bar{g}_{n+1}, \bar{h}_{n+1}$) à l'itération (n + 1) (avec n ≥ 0) se fait par la séquence suivante :
$$
\left\{
\begin{array}{ll}
1- & \text{on calcule } \alpha_n = \frac{\bar{\bar{g}}_n \bar{\bar{g}}_n}{\bar{\bar{h}}_n A \bar{\bar{h}}_n} \\[2ex]
2- & \text{on en déduit } \bar{x}_{n+1} = \bar{x}_n - \alpha_n \bar{h}_n \text{ et } \bar{\bar{g}}_{n+1} = \bar{\bar{g}}_n - \alpha_n A \bar{\bar{h}}_n \\[2ex]
3- & \text{on calcule } \beta_n = \frac{\bar{\bar{g}}_{n+1} \bar{\bar{g}}_{n+1}}{\bar{\bar{g}}_n \bar{\bar{g}}_n} \\[2ex]
4- & \text{on en déduit } \bar{h}_{n+1} = \bar{\bar{g}}_{n+1} + \beta_n \bar{h}_n
\end{array}
\right.
\qquad (37)
$$

3. L'itération du processus est arrêtée lorsque le vecteur $\bar{x}_{n+1}$ est jugé suffisamment proche de son vecteur prédécesseur $\bar{x}_n$. A ce sujet, rappelons que depuis le début de ce chapitre nous travaillons avec des températures réduites $\theta$ toutes comprises entre 0 et 1. Il s'ensuit que le vecteur $\bar{x}$ solution du système linéaire $A\bar{x} = \bar{b}$ ainsi que les différents vecteurs $\bar{x}_n$ calculés par l'algorithme de Fletcher et Reeves doivent avoir toutes leurs coordonnées elles aussi comprises entre 0 et 1, de sorte qu'un critère de convergence de l'algorithme peut être exprimé tout simplement par :

$$
\| \bar{x}_{n+1} - \bar{x}_n \|_{\infty} < \text{un certain } \varepsilon \quad (38)
$$

la notation $\|\dots\|_\infty$ désignant la norme infinie classique.

L'algorithme de Fletcher et Reeves converge assez rapidement. A titre d'illustration, nous avons choisi l'échantillon d'étude représenté à la figure 8 : cet échantillon bidimensionnel comporte 400 × 400 pixels et est composé de deux constituants, les pixels noirs de conductivité thermique choisie égale à 5 W/m/K et les pixels blancs de conductivité thermique choisie égale à 1 W/m/K (cet échantillon peut être considéré comme un échantillon composite unidirectionnel, les pixels noirs et les pixels blancs matérialisant respectivement des fibres et de la matrice). Nous avons calculé le champ de température θ (c'est-à-dire les 400 × 400 inconnues θᵢⱼ) s'établissant au sein de l'échantillon à l'aide de deux algorithmes : d'une part l'algorithme de Fletcher et Reeves qui a été présenté plus haut et d'autre part un algorithme très connu dit "de relaxation" [4]. Le graphe de la figure 9 montre l'évolution en fonction de l'indice n d'itération des deux normes $\|\bar{x}_{n+1} - \bar{x}_n\|_\infty$ caractérisant la convergence des deux processus numériques. Cette figure est éloquente : certes la convergence de l'algorithme de Fletcher et Reeves est-elle plus chaotique que celle de l'algorithme de relaxation, mais aussi tellement plus rapide ! Au bout de 2202 itérations produites en 9.7 minutes de calcul sur un micro-ordinateur classique (fréquence de l'horloge = 400 MHz - taille de mémoire vive = 256 Mo), l'algorithme de Fletcher et Reeves calcule une norme $\|\bar{x}_{n+1} - \bar{x}_n\|_\infty$ inférieure à $10^{-9}$ (valeur seuil retenue dans le cadre de cet exemple) et le programme s'arrête. En comparaison, après 25.1 minutes et 5000 itérations, la norme $\|\bar{x}_{n+1} - \bar{x}_n\|_\infty$ est toujours supérieure à $10^{-5}$ avec l'algorithme de relaxation.

Un article récemment paru [5], traitant de la modélisation de la conduction thermique dans le silicium poreux et abordant le calcul de la conductivité thermique effective de façon tout à fait
analogue à notre approche, fait appel à un algorithme itératif particulier de résolution de systèmes linéaires. Cet algorithme baptisé BASIS (pour "Block-correction Aided Strongly Implicit Solver") est présenté par ses concepteurs comme étant plus puissant et plus robuste que la plupart des algorithmes itératifs classiques (voir référence [6]). A l'avenir, si les tailles des systèmes linéaires et par conséquent les temps de calcul et les quantités de mémoire informatique augmentent de façon importante (on pense en particulier à des applications tridimensionnelles), il pourra donc être opportun d'implanter cet algorithme BASIS et de comparer ses performances à celles de l'algorithme de Fletcher et Reeves.

Le matériau composite unidirectionnel que nous avons choisi comme milieu d'étude dans ce chapitre sera l'objet d'études plus approfondies au chapitre 3-1-5.

## 3- Présentation de quelques résultats

Dans cette partie, nous allons tester notre programme de calcul sur des exemples de complexité croissante. Ces exemples seront pour la plupart en géométrie bidimensionnelle ; nous présenterons néanmoins au chapitre 3-2- quelques calculs 3-D qui ont été effectués dans le cadre d'une action de diffusion technologique.

### 3-1- Exemples en géométrie bidimensionnelle

#### 3-1-1- Milieu homogène - conductivité thermique fonction de la position - résistances thermiques en parallèle

Le problème le plus simple à résoudre, après celui du milieu homogène isotrope, est celui d'un milieu présentant certes un champ de conductivité thermique $\lambda$ non uniforme mais uniquement fonction de la variable y (se référer aux notations de la figure 1). Dans cette configuration, il est aisé de concevoir que le champ de température $\theta$ sera uniquement fonction de la variable x : en effet, le champ $\lambda$ étant uniquement fonction de y, l'échantillon se présente comme un groupement de résistances thermiques en parallèle.

Pour plus de clarté, nous allons prendre l'exemple d'un échantillon dont la conductivité thermique $\lambda$, uniquement fonction de la variable y, s'exprime par $\lambda(y) = 10 + 9 \cos(6\pi \frac{y}{\lambda})$ (en W/m/K) ; le champ de conductivité thermique dans l'échantillon est représenté sous forme d'image (200 x 200 points) à la figure 10. On conçoit aisément que ce genre de structure thermique n'aura aucune incidence particulière sur le champ de température par rapport à un milieu homogène isotrope : le champ de température restera une fonction affine de la seule variable x (voir expression (3) qui en version réduite devient $\theta(x) = 1 - \frac{x}{\lambda}$), le flux thermique restera orienté

selon l'axe (Ox) et se partagera entre les différentes résistances thermiques constituant l'échantillon. Ces différents points sont confirmés par notre programme de calcul : le champ de température qu'il prévoit est représenté sous forme d'image à la figure 11. Conformément à nos attentes, ce champ est bien fonction de la seule variable x, et l'analyse de son profil en x (voir figure 12) révèle bien un comportement affine.
Concernant la conductivité thermique effective $\Lambda$ de la structure thermique considérée, là encore cette quantité peut être anticipée. Le flux thermique étant en tout point de la structure orienté selon l'axe (Ox), l'échantillon apparaît dès lors comme constitué de résistances thermiques en parallèle. Or, des résistances thermiques en parallèle signifiant des conductances thermiques en série, il s'ensuit que la conductivité thermique effective $\Lambda$ de la structure s'exprime simplement par :

$$ \Lambda = \frac{1}{Y} \int_{y=0}^{Y} \lambda(y) dy = \frac{1}{Y} \int_{y=0}^{Y} \left( 10 + 9 \cos\left(6\pi \frac{x}{Y}\right) \right) dy = 10 \text{ W/m/K} \quad (39) $$

Cette valeur est confirmée par le calcul numérique : celui-ci aboutit à une conductivité thermique effective moyenne $\bar{\Lambda}$ (définie en (34)) valant exactement 10 W/m/K et un écart relatif maximal $\epsilon_{\Lambda}$ (défini au chapitre 2-3-) de $8.5 \times 10^{-15}$ ...

### 3-1-2- Milieu homogène - conductivité thermique fonction de la position - résistances thermiques en série

Nous allons maintenant compliquer quelque peu le problème en positionnant les résistances thermiques non plus en parallèle mais en série. Le champ de conductivité thermique $\lambda$ va donc rester fonction d'une seule coordonnée d'espace, mais ladite coordonnée d'espace ne sera plus la variable y mais la variable x. La fonction $\lambda(x)$ qu'on retiendra dans ce chapitre sera identique à celle adoptée au chapitre précédent : $\lambda(x) = 10 + 9 \cos(6\pi \frac{x}{X})$ (en W/m/K) ; le champ de conductivité thermique dans l'échantillon est représenté sous forme d'image (200 x 200 points) à la figure 13.

Cette fois, le champ de température s'établissant au sein de l'échantillon est un peu plus compliqué à déterminer : certes il reste fonction de la seule variable x, mais il faut passer par la résolution du système (11) pour obtenir son expression mathématique. Le champ $\theta$ étant indépendant de la variable y, le système (11) se simplifie et devient :

$$ \left\{ \begin{array}{l} \displaystyle \frac{d}{dx} \left( \lambda(x) \frac{d\theta}{dx} \right) = \frac{d}{dx} \left[ \left( 10 + 9 \cos \left( 6\pi \frac{x}{X} \right) \right) \frac{d\theta}{dx} \right] = 0 \quad \text{pour } 0 \le x \le X \\[1em] \theta = 1 \quad \text{en } x = 0 \\[1em] \theta = 0 \quad \text{en } x = X \end{array} \right. \quad (40) $$

La résolution de ce problème différentiel fournit la solution analytique suivante pour le champ $\theta(x)$ :
$$
\theta(x) = \begin{cases} 1 - \frac{1}{3\pi} \operatorname{arctg} \left( \frac{1}{\sqrt{19}} \operatorname{tg} \left( 3\pi \frac{x}{X} \right) \right) & \text{pour } 0 \le \frac{x}{X} \le \frac{1}{6} \\ \frac{2}{3} - \frac{1}{3\pi} \operatorname{arctg} \left( \frac{1}{\sqrt{19}} \operatorname{tg} \left( 3\pi \frac{x}{X} \right) \right) & \text{pour } \frac{1}{6} \le \frac{x}{X} \le \frac{1}{3} \\ \theta \left( \frac{x}{X} - \frac{1}{3} \right) - \frac{1}{3} & \text{pour } \frac{1}{3} \le \frac{x}{X} \le 1 \end{cases} \quad (41)
$$

Maintenant que nous connaissons la solution analytique du problème posé, observons le résultat
numérique fourni par notre programme de calcul. Le champ de température obtenu par calcul est
représenté sous forme d'image à la figure 14 : ce champ est bien fonction de la seule variable x ; par
ailleurs, son profil en x est comparé à la solution analytique (41) à la figure 15, et on constate un
excellent accord entre les résultats numérique et analytique.

La conductivité thermique effective Λ de la structure thermique considérée peut être évaluée par un
raisonnement analogue à celui présenté au chapitre précédent. Le flux thermique étant en tout point
de la structure orienté selon l'axe (Ox), l'échantillon apparaît dès lors comme constitué de
résistances thermiques en série. Or, des résistances thermiques en série signifiant des conductances
thermiques en parallèle, il s'ensuit que la conductivité thermique effective Λ de la structure
s'obtient par l'égalité :

$$
\frac{1}{\Lambda} = \frac{1}{X} \int_{x=0}^{x} \frac{dx}{\lambda(x)} = \frac{1}{X} \int_{x=0}^{x} \frac{dx}{10 + 9 \cos\left(6\pi \frac{x}{X}\right)} = \frac{1}{\sqrt{19}} \quad (42)
$$

soit $\Lambda = \sqrt{19} \approx 4.3589$ W/m/K. Cette évaluation analytique de $\Lambda$ est parfaitement reproduite par le calcul numérique, puisque celui-ci aboutit au couple $(\bar{\Lambda}, \varepsilon_{\Lambda}) = (4.3589 \text{ W/m/K}, 3.9 \times 10^{-8})$.

3-1-3- Milieu homogène - conductivité thermique fonction de la position - résistances thermiques "en produit série-parallèle"

Le cas que nous allons maintenant traiter constitue le prolongement direct des études présentées aux
deux chapitres précédents. Nous allons considérer une structure thermique dite "en produit série-
parallèle", c'est-à-dire caractérisée par un champ de conductivité thermique λ fonction certes des
deux variables x et y mais de façon séparée. En d'autres termes, le champ λ peut s'écrire sous la
forme d'un produit d'une fonction de la seule variable x par une fonction de la seule variable y :

$$
\lambda(x,y) = \lambda_x(x) \times \lambda_y(y) \tag{43}
$$
Pour cette étude, nous avons retenu les expressions suivantes pour les deux fonctions $\lambda_x(x)$ et $\lambda_y(y)$ :

$$ \lambda_x(x) = 10 + 9 \cos\left(6\pi \frac{x}{X}\right) \text{ et } \lambda_y(y) = 1 + 9 \left(\frac{y}{Y}\right)^2 $$

Le champ de conductivité thermique dans l'échantillon est représenté sous forme d'image (200 x 200 points) à la figure 16.

Au premier abord, on est tenté de se dire qu'avec une telle distribution de conductivité thermique, le champ de température s'établissant dans l'échantillon ne peut être obtenu que de façon numérique. Cependant, il s'avère que le champ $\theta$ calculé par notre programme (et présenté sous forme d'image à la figure 17) est uniquement fonction de la variable x et que de plus il est exactement égal à celui obtenu dans le chapitre précédent. Ce dernier point est confirmé par la comparaison, faite à la figure 18, du profil en x du champ $\theta$ à la solution analytique (41).

Ce résultat numérique, quoique surprenant, est juste. Le système à résoudre est le système (11) ; remplaçant dans ce système le champ $\lambda$ par le produit $\lambda_x(x) \times \lambda_y(y)$ d'une part et postulant un champ $\theta$ uniquement fonction de la variable x d'autre part, le système (11) devient :

$$ egin{cases} 
\frac{d}{dx} \left( \lambda_x(x) \frac{d\theta}{dx} \right) = 0 \quad \text{pour } 0 \le x \le X \\ 
\theta = 1 \quad \text{en } x = 0 \\ 
\theta = 0 \quad \text{en } x = X 
\end{cases} 
\quad (44) $$

et on a là un système d'une part complètement débarrassé de la composante en y du champ de conductivité thermique et d'autre part rigoureusement identique au système (40) du chapitre précédent. D'où ce résultat assez étonnant : lorsque le champ de conductivité thermique est à variables séparées, le champ de température s'établissant dans l'échantillon n'est fonction ni de y ni de la composante en y du champ de conductivité thermique ; le champ de température est uniquement le reflet de la composante en x du champ de conductivité thermique.

L'évaluation analytique de la conductivité thermique effective $\Lambda$ se fait de la façon suivante. Considérons une section de l'échantillon parallèle aux deux faces à température imposée et située à l'abscisse x ; on sait (se référer au chapitre 1-) que la puissance thermique $\Psi$ traversant cette section est indépendante de x et vaut $\frac{\Lambda Y}{X}$. Par ailleurs, $\Psi$ s'exprime aussi par :

$$ \Psi = \int_{y=0}^{Y} \psi_x dy = - \int_{y=0}^{Y} \lambda \frac{\partial \theta}{\partial x} dy = - \lambda_x(x) \frac{d\theta}{dx} \int_{y=0}^{Y} \lambda_y(y) dy \quad (45) $$

Le système (44) peut facilement être intégré, et aboutit à l'expression suivante du champ de température $\theta$ :
$$ \theta(x) = \frac{\int_{0}^{x} \frac{dx}{\lambda_{x}(x)}}{\int_{0}^{x} \frac{dx}{\lambda_{x}(x)}} \qquad (46) $$

De cette expression on tire :

$$ -\lambda_x(x) \frac{d\theta}{dx} = \frac{1}{\int_0^x \frac{dx}{\lambda_x(x)}} \qquad (47) $$

La formule (45) de la puissance thermique Ψ peut alors être complétée, et on aboutit finalement à l'expression analytique suivante de la conductivité thermique effective Λ :

$$ \Lambda = \frac{\frac{1}{Y} \int_{0}^{Y} \lambda_{y}(y) dy}{\frac{1}{X} \int_{0}^{X} \frac{dx}{\lambda_{x}(x)}} \qquad (48) $$

Cette formule (48) justifie l'appellation de résistances thermiques "en produit série-parallèle" que nous avons adoptée pour décrire notre milieu : en effet, elle apparaît comme étant le produit de deux conductivités thermiques effectives, la première (numérateur de l'expression (48)) étant celle de résistances thermiques en parallèle (voir chapitre 3-1-1- et formule (39)) et la seconde (inverse du dénominateur de l'expression (48)) étant celle de résistances thermiques en série (voir chapitre 3-1-2- et formule (42)). Remplaçant les fonctions $\lambda_x(x)$ et $\lambda_y(y)$ par leurs expressions respectives (en l'occurrence $10 + 9 \cos(6\pi \frac{x}{X})$ et $1 + 9(\frac{y}{Y})^2$) dans l'expression (48), on prévoit par calcul intégral une valeur de $4\sqrt{19} \approx 17.436 \text{ W/m/K}$ à la conductivité thermique effective $\Lambda$ ; encore une fois, cette valeur théorique est reproduite par notre programme de calcul avec une excellente précision (celui-ci aboutit au couple $(\bar{\Lambda}, \varepsilon_{\Lambda}) = (17.436 \text{ W/m/K}, 2.4 \times 10^{-9})$).

### 3-1-4- Milieu à deux constituants et périodique

Dans les trois chapitres précédents, nous avons testé - avec succès - notre programme de calcul sur des configurations d'école. Nous allons maintenant passer à une catégorie de matériaux plus complexes et donc beaucoup plus difficiles à analyser de façon analytique, mais éminemment intéressants car aux applications pratiques très répandues : les matériaux périodiques.
Dans le but d'illustrer les possibilités offertes par notre programme, nous avons choisi de travailler sur un matériau particulier de type nid d'abeille (c'est-à-dire à cellules hexagonales). Dans une première partie, nous allons présenter différents résultats de calcul fournis par le programme et réfléchir à la notion de volume élémentaire représentatif; ensuite, nous décrirons une autre approche plus analytique permettant de résoudre le problème thermique (11) par le biais de l'analyse spectrale.

### 3-1-4-1- Calcul de champs de température à l'aide de notre programme en vue de la définition d'un volume élémentaire représentatif

La maille bidimensionnelle que nous avons retenue pour engendrer nos échantillons périodiques est illustrée sous forme d'image à la figure 19(a) : cette maille compte 72 × 42 pixels et est constituée de deux matériaux, les pixels noirs de conductivité thermique égale à 5 W/m/K et les pixels blancs de conductivité thermique égale à 1 W/m/K.

A partir de cette maille, répétée par translations successives dans les directions horizontale et verticale du plan, nous avons engendré dix échantillons périodiques bornés de tailles croissantes, le plus petit d'entre eux étant la maille elle-même et les neuf autres étant constitués respectivement de 2 × 2 mailles, 3 × 3 mailles, 4 × 4 mailles et ainsi de suite jusqu'à 10 × 10 mailles. La figure 19 montre cinq des dix échantillons ainsi engendrés (en l'occurrence les échantillons "1 × 1" (c'est-à-dire la maille), "2 × 2", "3 × 3", "5 × 5" et "9 × 9").

A l'aide de notre programme de calcul, nous avons déterminé les champs de température θ s'établissant dans les échantillons puis les conductivités thermiques effectives Λ de ceux-ci. Les différents résultats de ces calculs sont consignés dans le tableau 1 ci-dessous.

<table><caption>Tableau 1</caption><thead><tr><td>échantillon</td><td>nombre de températures inconnues</td><td>temps de calcul (minutes)</td><td>Λ̄ (W/m/K)</td><td>ε<sub>Λ</sub></td></tr></thead><tbody><tr><td>"1 × 1"</td><td>72 × 42</td><td>0.0</td><td>2.0430</td><td>2.1 10<sup>-9</sup></td></tr><tr><td>"2 × 2"</td><td>144 × 84</td><td>0.2</td><td>2.0299</td><td>9.9 10<sup>-9</sup></td></tr><tr><td>"3 × 3"</td><td>216 × 126</td><td>0.7</td><td>2.0254</td><td>1.1 10<sup>-8</sup></td></tr><tr><td>"4 × 4"</td><td>288 × 168</td><td>1.5</td><td>2.0232</td><td>2.3 10<sup>-8</sup></td></tr><tr><td>"5 × 5"</td><td>360 × 210</td><td>2.9</td><td>2.0219</td><td>2.3 10<sup>-8</sup></td></tr><tr><td>"6 × 6"</td><td>432 × 252</td><td>4.8</td><td>2.0210</td><td>2.7 10<sup>-8</sup></td></tr><tr><td>"7 × 7"</td><td>504 × 294</td><td>7.3</td><td>2.0203</td><td>2.2 10<sup>-8</sup></td></tr><tr><td>"8 × 8"</td><td>576 × 336</td><td>10.6</td><td>2.0198</td><td>2.3 10<sup>-8</sup></td></tr><tr><td>"9 × 9"</td><td>648 × 378</td><td>14.6</td><td>2.0195</td><td>2.9 10<sup>-8</sup></td></tr><tr><td>"10 × 10"</td><td>720 × 420</td><td>20.2</td><td>2.0192</td><td>5.2 10<sup>-8</sup></td></tr></tbody></table>

Notons tout d'abord deux points dans ce tableau : d'une part la faiblesse des valeurs des écarts relatifs maximaux ε<sub>Λ</sub> (les dix valeurs sont inférieures à 10<sup>-7</sup>), synonyme de fiabilité des résultats
numériques, et d'autre part l'efficacité de l'algorithme de Fletcher et Reeves : cet algorithme ne met
que 20 minutes pour résoudre un système linéaire à 720 × 420 = 302400 inconnues.

Les champs de température s'établissant dans les cinq échantillons de la figure 19 sont représentés sous forme d'images et de faisceaux de courbes isothermes à la figure 20. L'observation des cinq images et plus précisément des cinq faisceaux d'isothermes mène à la constatation intéressante suivante : plus l'échantillon comporte de mailles, plus les isothermes tendent vers des droites horizontales régulièrement espacées ; en d'autres termes, plus l'échantillon comporte de mailles, plus le champ de température θ se rapproche du champ qui s'établit dans un milieu homogène isotrope et qui, rappelons-le, a pour expression θ(x) = 1 - x/X. Ce point est confirmé à la figure 21 qui expose sous forme d'images les champs dits "de déviation thermique" régnant dans les cinq échantillons de la figure 19, ces champs de déviation thermique notés δθ étant définis par les différences entre les champs de température θ évalués par notre programme de calcul et le champ θ(x) = 1 - x/X s'établissant dans un milieu homogène isotrope. On constate en effet sur la figure 21 que plus l'échantillon comporte de mailles, plus la palette de valeurs couverte par le champ δθ se resserre autour de 0. On peut d'ailleurs tracer la courbe de la norme infinie du champ δθ (définie par $\max_{i,j} |\delta\theta_{i,j}|$) en fonction du nombre de mailles dans chacune des directions horizontale et verticale de l'échantillon (voir courbe en trait continu de la figure 22), et cette courbe montre bien une tendance asymptotique vers la valeur 0 quand le nombre de mailles augmente.

Les valeurs des conductivités thermiques effectives moyennes $\bar{\Lambda}$ indiquées dans le tableau 1 sont intéressantes en ce sens qu'elles viennent compléter ce qui vient d'être dit concernant les champs de déviation thermique $\delta\theta$. La courbe montrant l'évolution de $\bar{\Lambda}$ en fonction du nombre de mailles dans chacune des directions horizontale et verticale de l'échantillon est tracée à la figure 22 (courbe en trait pointillé). Cette courbe, à l'instar de celle de la norme infinie du champ $\delta\theta$, présente un comportement asymptotique quand le nombre de mailles augmente. Sa valeur limite correspond à un échantillon contenant un nombre de mailles infini, donc à un échantillon dans lequel le champ de température $\theta$ qui s'établit est égal à celui du milieu homogène isotrope ; par conséquent, la valeur limite de la courbe $\bar{\Lambda}$ est la conductivité thermique effective "vraie", c'est-à-dire celle du milieu homogène véritablement équivalent à notre milieu hétérogène d'étude constitué de cellules hexagonales.

Pour clore ce chapitre, notons que l'échantillon "1 × 1" (c'est-à-dire l'échantillon constitué d'une seule maille) fournit une conductivité thermique effective déjà assez proche de la valeur asymptotique de la courbe $\bar{\Lambda}$, et que la figure 22 montre par la suite une convergence assez lente de la courbe $\bar{\Lambda}$ vers ladite valeur asymptotique. De ces deux remarques il faut tirer les deux conséquences suivantes : d'une part, pour une évaluation approximative de la conductivité thermique effective d'un milieu hétérogène, il n'est peut-être pas nécessaire de considérer un volume matériel trop important ; d'autre part, l'obtention d'un résultat véritablement homogénéisé, c'est-à-dire en l'occurrence d'un champ de température à l'image de celui qui s'établirait dans un milieu homogène, requiert en revanche la prise en compte d'un volume matériel de taille conséquente.

### 3-1-4-2- Application de l'analyse spectrale à la détermination du champ de température s'établissant au sein d'un milieu périodique
### 3-1-4-2-1- Enoncé du problème à résoudre

La figure 21, qui montre quelques-uns des dix champs de déviation thermique $\delta\theta$ évalués par notre programme de calcul, est assez évocatrice : elle semble indiquer qu’un champ de conductivité thermique $\lambda$ périodique se traduit par un champ $\delta\theta$ lui aussi périodique et présentant la même période que le champ $\lambda$. Cette tendance se trouve confirmée à la figure 23 qui présente sous forme d’images les spectres bidimensionnels d’une part de l’échantillon "10 x 10" et d’autre part du champ de déviation thermique $\delta\theta$ obtenu sur cet échantillon : le spectre du champ $\delta\theta$ est discret, et les positions de ses pics se superposent parfaitement à celles des pics du spectre de l’échantillon.

Il est donc tentant d’essayer d’appliquer la théorie des séries de Fourier au problème mathématique (11) que nous avons à résoudre et d’étudier ce que cette approche spectrale peut apporter de nouveau : en particulier, permet-elle de passer des amplitudes complexes des pics du spectre bidimensionnel de l’échantillon aux amplitudes complexes des pics du spectre bidimensionnel du champ de déviation thermique $\delta\theta$ ?

### 3-1-4-2-2- Détermination de la série de Fourier caractérisant le champ de déviation thermique $\delta\theta$

Le champ de conductivité thermique $\lambda(x,y)$ étant périodique et défini par une maille comptant $N_x$ pixels dans la direction x et $N_y$ pixels dans la direction y, ce champ a donc pour fréquences spatiales fondamentales dans les directions x et y respectivement les quantités $\sigma_x = \frac{1}{N_x\delta}$ et $\sigma_y = \frac{1}{N_y\delta}$ (on rappelle que $\delta$ désigne la distance dans les directions x et y entre deux pixels voisins) ; il peut par conséquent s’écrire sous la forme de la série de Fourier suivante :

$$ \lambda(x,y) = \sum_{m=-\infty}^{+\infty} \sum_{n=-\infty}^{+\infty} \lambda_{m,n} h_m^x(x) h_n^y(y) \quad (49) $$

Dans cette expression, les $\lambda_{m,n}$ sont les coefficients de la série de Fourier bidimensionnelle de la fonction $\lambda(x,y)$, et les fonctions $h_m^x(x)$ et $h_n^y(y)$ désignent les harmoniques définies par :

$$ h_m^x(x) = \exp(2i\pi m \sigma_x x) \quad \text{et} \quad h_n^y(y) = \exp(2i\pi n \sigma_y y) \quad (50) $$

De même, supposant le champ de déviation thermique $\delta\theta(x,y)$ lui aussi périodique et défini par la même maille, on peut écrire ce champ sous la forme :

$$ \delta\theta(x,y) = \sum_{m=-\infty}^{+\infty} \sum_{n=-\infty}^{+\infty} \delta\theta_{m,n} h_m^x(x) h_n^y(y) \quad (51) $$
les $\delta\theta_{m,n}$ étant les coefficients de la série de Fourier bidimensionnelle de la fonction $\delta\theta(x,y)$. Le but de cette étude est le suivant : connaissant la maille génératrice du champ $\lambda(x,y)$, c'est-à-dire les coefficients $\lambda_{m,n}$ de la série de Fourier bidimensionnelle de ce champ, peut-on déduire les coefficients $\delta\theta_{m,n}$ de la série de Fourier bidimensionnelle du champ $\delta\theta(x,y)$ et donc ce champ $\delta\theta(x,y)$ lui-même ?

Pour cela, il faut remonter à l'origine du problème à résoudre, c'est-à-dire à l'équation de volume du système (11). Cette équation de volume s'écrit sous la forme :

$$ \vec{\nabla} \cdot \vec{\psi} = \frac{\partial \psi_x}{\partial x} + \frac{\partial \psi_y}{\partial y} = 0 \qquad (52) $$

en introduisant le vecteur flux thermique $\vec{\psi}$ défini par $\vec{\psi} = -\lambda \vec{\nabla}\theta$ et ses composantes en x et y notées respectivement $\psi_x$ et $\psi_y$. Or le champ de température $\theta(x,y)$ vaut, par définition du champ $\delta\theta(x,y)$ :

$$ \theta(x,y) = \delta\theta(x,y) + \left(1 - \frac{x}{X}\right) = \sum_{m=-\infty}^{+\infty} \sum_{n=-\infty}^{+\infty} \delta\theta_{m,n} h_m^x(x) h_n^y(y) + \left(1 - \frac{x}{X}\right) \qquad (53) $$

Le vecteur $\vec{\nabla}\theta$ a donc pour coordonnées :

$$ \left\{ \begin{aligned} \frac{\partial \theta}{\partial x} &= 2i\pi\sigma_x \sum_{m=-\infty}^{+\infty} \sum_{n=-\infty}^{+\infty} m \delta\theta_{m,n} h_m^x(x) h_n^y(y) - \frac{1}{X} \\ \frac{\partial \theta}{\partial y} &= 2i\pi\sigma_y \sum_{m=-\infty}^{+\infty} \sum_{n=-\infty}^{+\infty} n \delta\theta_{m,n} h_m^x(x) h_n^y(y) \end{aligned} \right. \qquad (54) $$

La composante $\psi_x$ du flux thermique $\vec{\psi}$, qui est égale - au signe près - au produit du champ $\lambda$ (périodique) par le champ $\frac{\partial\theta}{\partial x}$ (lui aussi périodique), apparaît dès lors comme étant un champ périodique et peut donc s'écrire sous la forme :

$$ \psi_x(x,y) = \sum_{m=-\infty}^{+\infty} \sum_{n=-\infty}^{+\infty} \psi_{m,n}^x h_m^x(x) h_n^y(y) \qquad (55) $$

Le même argument appliqué à la composante $\psi_y$ du flux thermique $\vec{\psi}$ permet d'écrire ce champ sous une forme analogue :
$$ \psi_y(x,y) = \sum_{m=-\infty}^{+\infty} \sum_{n=-\infty}^{+\infty} \psi_{m,n}^y h_m^x(x) h_n^y(y) \quad (56) $$

Les deux composantes $\psi_x$ et $\psi_y$ étant liées par l'équation de volume (52), il est facile d'en déduire la relation suivante entre les coefficients de Fourier $\psi_{m,n}^x$ et $\psi_{m,n}^y$ de ces deux champs :

$$ \forall m \in \mathbb{Z}, \forall n \in \mathbb{Z}, m \sigma_x \psi_{m,n}^x + n \sigma_y \psi_{m,n}^y = 0 \quad (57) $$

Exprimons maintenant les coefficients de Fourier $\psi_{m,n}^x$ et $\psi_{m,n}^y$ du champ de flux thermique $\bar{\psi}$ en fonction de ceux $\delta\theta_{m,n}$ du champ de déviation thermique $\delta\theta$. Le développement des produits $\lambda \frac{\partial\theta}{\partial x}$ et $\lambda \frac{\partial\theta}{\partial y}$ fait apparaître des convolutions de la série de Fourier $(\lambda_{m,n})_{m \in \mathbb{Z}, n \in \mathbb{Z}}$ du champ $\lambda$ avec la série de Fourier $(\delta\theta_{m,n})_{m \in \mathbb{Z}, n \in \mathbb{Z}}$ du champ $\delta\theta$ ; les expressions obtenues sont les suivantes :

$$ \left\{ \begin{aligned} \psi_{m,n}^x &= \frac{\lambda_{m,n}}{X} - 2i\pi\sigma_x \sum_{p=-\infty}^{+\infty} \sum_{q=-\infty}^{+\infty} \lambda_{m-p,n-q} p \delta\theta_{p,q} \\ \psi_{m,n}^y &= -2i\pi\sigma_y \sum_{p=-\infty}^{+\infty} \sum_{q=-\infty}^{+\infty} \lambda_{m-p,n-q} q \delta\theta_{p,q} \end{aligned} \right. \quad (58) $$

L'injection de ces expressions (58) dans l'équation de volume (57) aboutit à l'ensemble d'équations suivant :

$$ \forall m \in \mathbb{Z}, \forall n \in \mathbb{Z}, \sum_{p=-\infty}^{+\infty} \sum_{q=-\infty}^{+\infty} (mp+nq^2) \lambda_{m-p,n-q} \delta\theta_{p,q} = \frac{m\lambda_{m,n}}{2i\pi\sigma_x X} \quad (59) $$

dans lequel on a introduit la quantité $\sigma = \frac{\sigma_y}{\sigma_x} = \frac{N_x}{N_y}$. Or, $\sigma_x$ étant l'inverse de la longueur de la maille selon la direction x, il s'ensuit que le produit $\sigma_x X$ est égal au nombre de mailles de l'échantillon considéré dans cette même direction. Baptisant $P_x$ ce nombre, et l'insérant dans (59), on aboutit à l'ensemble d'équations suivant :

$$ \forall m \in \mathbb{Z}, \forall n \in \mathbb{Z}, \sum_{p=-\infty}^{+\infty} \sum_{q=-\infty}^{+\infty} (mp+nq^2) \lambda_{m-p,n-q} \delta\theta_{p,q} = \frac{m\lambda_{m,n}}{2i\pi P_x} \quad (60) $$
dans lequel, rappelons-le, la série de Fourier $(\lambda_{m,n})_{m \in Z, n \in Z}$ du champ $\lambda$ est supposée connue (car la maille de l'échantillon périodique est donnée) et la série de Fourier $(\delta\theta_{m,n})_{m \in Z, n \in Z}$ du champ $\delta\theta$ constitue l'ensemble de nos inconnues. L'ensemble d'équations (60) apparaît dès lors comme étant un système linéaire des inconnues $(\delta\theta_{m,n})_{m \in Z, n \in Z}$ ; le problème est donc formellement résolu.

### 3-1-4-2-3- Le champ $\delta\theta$ est-il inversement proportionnel à $P_x$ ?

Il est intéressant de constater que les termes à droite du signe égal dans l'ensemble d'équations (60) sont inversement proportionnels au nombre de mailles $P_x$ dans l'épaisseur X de l'échantillon ; par linéarité, il s'ensuit que la série de Fourier $(\delta\theta_{m,n})_{m \in Z, n \in Z}$ et donc le champ $\delta\theta$ sont eux aussi inversement proportionnels à $P_x$. Cette remarque vient compléter les observations faites au chapitre 3-1-4-1- : lorsque le nombre de mailles $P_x$ dans l'épaisseur X de l'échantillon tend vers l'infini, le champ de déviation thermique $\delta\theta$ tend effectivement vers 0 car ce champ $\delta\theta$ est en $O(1/P_x)$.

Un corollaire de ce résultat est que lorsque $P_x$ tend vers l'infini, le champ $P_x \times \delta\theta$ tend vers un champ constant, propriété que l'on peut essayer de vérifier sur la base des dix champs de déviation thermique $\delta\theta$ évalués au chapitre précédent. La figure 24 présente, pour chacun des dix échantillons traités, le champ $\delta\theta$ moyenné sur toutes les mailles constitutives de l'échantillon puis multiplié par le nombre $P_x$ correspondant (ainsi, pour l'échantillon "4 x 4", le champ $\delta\theta$ a été moyenné sur les 16 mailles de l'échantillon puis multiplié par 4). Définissant l'écart relatif entre deux images $I_1$ et $I_2$ par

$$ \text{la quantité} \quad \frac{\max(|I_2 - I_1|)}{\min(\max(I_1) - \min(I_1), \max(I_2) - \min(I_2))} \quad (\text{on remarquera que le numérateur de cette}$$

expression n'est autre que la norme infinie de l'image $(I_2 - I_1)$ tandis que son dénominateur a pour effet de rapporter cette distance entre les deux images à leurs amplitudes respectives), on peut dès lors quantifier cet écart relatif pour deux images successives de la figure 24 ; les calculs aboutissent aux valeurs numériques suivantes :

<table><caption>Tableau 2</caption><thead><tr><th>couple d'images de la figure 24</th><th>écart relatif entre les deux images</th></tr></thead><tbody><tr><td>"1 x 1" et "2 x 2"</td><td>0.176</td></tr><tr><td>"2 x 2" et "3 x 3"</td><td>0.064</td></tr><tr><td>"3 x 3" et "4 x 4"</td><td>0.033</td></tr><tr><td>"4 x 4" et "5 x 5"</td><td>0.020</td></tr><tr><td>"5 x 5" et "6 x 6"</td><td>0.013</td></tr><tr><td>"6 x 6" et "7 x 7"</td><td>0.010</td></tr><tr><td>"7 x 7" et "8 x 8"</td><td>0.007</td></tr><tr><td>"8 x 8" et "9 x 9"</td><td>0.006</td></tr><tr><td>"9 x 9" et "10 x 10"</td><td>0.005</td></tr></tbody></table>
Au vu des valeurs numériques portées dans le tableau 2 ci-dessus, il apparaît que l'écart relatif entre deux images successives de la figure 24 diminue lorsque les tailles des deux images augmentent. Par ailleurs, lorsque les deux images considérées sont les images "9 × 9" et "10 × 10", cet écart relatif tombe à 5‰ : on peut donc raisonnablement considérer que ces deux images sont égales, et que par conséquent la proposition "$P_x \times \delta\theta = \text{un champ constant}$" semble être vérifiée par nos calculs.

### 3-1-4-2-4- Détermination de la conductivité thermique effective $\Lambda$ du milieu périodique

Une fois la série de Fourier bidimensionnelle $(\delta\theta_{m,n})_{m\in Z, n\in Z}$ évaluée, il est possible de déterminer la conductivité thermique effective $\Lambda$ du milieu périodique de la façon suivante. On part de l'expression générale (13) de $\Lambda$ que l'on réécrit sous la forme suivante :

$$ \Lambda = \frac{X}{S} \iint_S \psi_x \, dS \qquad (61) $$

On rappelle que l'intégrale double dans la formule ci-dessus peut être évaluée sur n'importe quelle surface S pourvu que celle-ci constitue une section plane de l'échantillon parallèle aux deux faces à température imposée ; la face chaude de l'échantillon $S_c$ (d'équation $x = 0$) vérifiant ces différents critères, on peut donc remplacer $S$ par $S_c$ dans la suite des calculs. Par ailleurs, notre échantillon périodique étant constitué de $P_x$ périodes dans la direction $x$ et de $P_y$ périodes dans la direction $y$ (soit donc au total de $P_x \times P_y$ périodes), la formule (61) devient donc :

$$ \Lambda = \frac{P_x N_x \delta}{P_y N_y \delta} \int_{y=0}^{P_y N_y \delta} \psi_x(x=0, y) \, dy = \frac{P_x N_x}{P_y N_y} \int_{y=0}^{P_y N_y \delta} \psi_x(x=0, y) \, dy \qquad (62) $$

La périodicité du champ $\psi_x$ et sa décomposition en série de Fourier bidimensionnelle (55) autorisent les développements mathématiques suivants :

$$ \begin{align} \Lambda &= \frac{P_x N_x}{P_y N_y} P_y \int_{y=0}^{N_y \delta} \psi_x(x=0, y) \, dy \notag \\ &= \frac{P_x N_x}{N_y} \sum_{m=-\infty}^{+\infty} \sum_{n=-\infty}^{+\infty} \psi_{m,n}^x h_m^x(x=0) \int_{y=0}^{N_y \delta} h_n^y(y) \, dy \tag{63} \\ &= \frac{P_x N_x}{N_y} \sum_{m=-\infty}^{+\infty} \sum_{n=-\infty}^{+\infty} \psi_{m,n}^x \int_{y=0}^{N_y \delta} h_n^y(y) \, dy \notag \end{align} $$
Enfin, remarquant que :

$$ \int_{y=0}^{N_y \delta} h_n^y(y) dy = \begin{cases} 0 & \text{pour } n \neq 0 \\ N_y \delta & \text{pour } n=0 \end{cases} \qquad (64) $$

on aboutit à l'expression suivante de $\Lambda$ :

$$ \Lambda = P_x N_x \delta \sum_{m=-\infty}^{+\infty} \psi_{m,0}^x \qquad (65) $$

Il ne reste plus maintenant qu'à remplacer dans (65) les termes $\psi_{m,0}^x$ par leurs formulations (58) ; n'oubliant pas que $\sigma_x = \frac{1}{N_x \delta}$ et que $X = P_x N_x \delta$, on atteint l'expression mathématique suivante de la conductivité thermique effective $\Lambda$ en fonction des deux séries de Fourier $(\lambda_{m,n})_{m \in Z, n \in Z}$ et $(\delta\theta_{m,n})_{m \in Z, n \in Z}$ :

$$ \Lambda = \sum_{m=-\infty}^{+\infty} \lambda_{m,0} - 2i\pi P_x \sum_{m=-\infty}^{+\infty} \sum_{p=-\infty}^{+\infty} \sum_{q=-\infty}^{+\infty} p \lambda_{m-p,-q} \delta\theta_{p,q} \qquad (66) $$

Or il est aisé de s'apercevoir, en repartant de la décomposition (49) du champ de conductivité thermique $\lambda$ en série de Fourier bidimensionnelle, que la première sommation dans (66) n'est autre que la valeur moyenne du champ $\lambda$ sur le segment d'équation $x = 0$ de la maille ; en d'autres termes, cette première sommation vaut exactement la quantité

$$ \bar{\lambda}_{x=0} = \frac{1}{N_y \delta} \int_{y=0}^{N_y \delta} \lambda(x=0, y) dy \qquad (67) $$

qui a une signification physique bien précise. Considérons en effet un milieu bidimensionnel périodique fictif dont la maille est définie de la façon suivante :

1. Elle présente les mêmes dimensions et le même profil de conductivité thermique sur le segment d'équation $x = 0$ que la maille de notre milieu bidimensionnel périodique réel ;

2. Elle présente un champ de conductivité thermique indépendant de la variable $x$.

A titre d'illustration, la figure 25 montre d'une part la maille génératrice du milieu de type nid d'abeille qui a été étudié aux chapitres précédents et d'autre part la maille génératrice du milieu
fictif correspondant. Il est clair qu'ainsi défini, ce milieu fictif est un milieu de type "résistances thermiques en parallèle" ; d'après les résultats énoncés au chapitre 3-1-1-, sa conductivité thermique effective, égale à la valeur moyenne de son champ de conductivité thermique selon la direction y, vaut donc exactement la quantité $\bar{\lambda}_{x=0}$ exprimée en (67). Par conséquent l'expression (66) de $\Lambda$, réécrite sous la forme :

$$ \Lambda = \bar{\lambda}_{x=0} - 2i\pi P_x \sum_{m=-\infty}^{+\infty} \sum_{p=-\infty}^{+\infty} \sum_{q=-\infty}^{+\infty} p \lambda_{m-p,-q} \delta\theta_{p,q} \quad (68) $$

fait apparaître la sommation triple dans la formule ci-dessus comme étant un terme correctif rendant compte de l'écart du milieu bidimensionnel périodique réel par rapport au milieu bidimensionnel périodique fictif associé constitué de résistances thermiques en parallèle (on pourra d'ailleurs remarquer que si le milieu bidimensionnel périodique d'étude est lui-même de type "résistances thermiques en parallèle", alors le champ de déviation thermique $\delta\theta$ est nul, sa série de Fourier bidimensionnelle $(\delta\theta_{m,n})_{m\in Z, n\in Z}$ aussi, et par conséquent l'égalité (68) se réduit à $\Lambda = \bar{\lambda}_{x=0}$, ce qui est effectivement la valeur de $\Lambda$ lorsque le milieu est de type "résistances thermiques en parallèle").

Dernier point à souligner : au premier coup d'œil sur l'expression (68), il peut paraître curieux que la valeur de $\Lambda$ soit fonction du nombre $P_x$ de périodes dans la direction x. En fait, il ne faut pas oublier (voir chapitre précédent) que la série $(\delta\theta_{m,n})_{m\in Z, n\in Z}$ est inversement proportionnelle à $P_x$ ; par conséquent, l'expression (68) ci-dessus fournit une valeur de la conductivité thermique effective $\Lambda$ qui est bien indépendante du nombre $P_x$.

### 3-1-4-2-5- Filtrage passe-bas des champs de conductivité thermique $\lambda$ et de déviation thermique $\delta\theta$

Revenons maintenant au système linéaire (60) permettant la détermination de la série de Fourier $(\delta\theta_{m,n})_{m\in Z, n\in Z}$ du champ $\delta\theta$ ; il s'agit bien d'un système linéaire des inconnues $(\delta\theta_{m,n})_{m\in Z, n\in Z}$

mais ... de dimension infinie. Pour contourner cette difficulté liée à la taille infinie du système, on a classiquement recours à une technique de filtrage passe-bas qui consiste à éliminer, dans tous les développements mathématiques qui viennent d'être présentés, les termes harmoniques $h_m^x(x)$ et $h_n^y(y)$ d'ordres m et n supérieurs en valeur absolue respectivement à des constantes M et N données. On a donc recours aux approximations suivantes :

$$ \left\{ \begin{aligned} \lambda(x,y) &= \sum_{m=-\infty}^{+\infty} \sum_{n=-\infty}^{+\infty} \lambda_{m,n} h_m^x(x) h_n^y(y) \approx \sum_{m=-M}^{+M} \sum_{n=-N}^{+N} \lambda_{m,n} h_m^x(x) h_n^y(y) \\ \delta\theta(x,y) &= \sum_{m=-\infty}^{+\infty} \sum_{n=-\infty}^{+\infty} \delta\theta_{m,n} h_m^x(x) h_n^y(y) \approx \sum_{m=-M}^{+M} \sum_{n=-N}^{+N} \delta\theta_{m,n} h_m^x(x) h_n^y(y) \end{aligned} \right. \quad (69) $$
et tous les termes $\lambda_{m,n}$ et $\delta\theta_{m,n}$ sont remplacés par 0 dès que $|m| > M$ ou $|n| > N$. Dès lors, le système linéaire (60) à une infinité d'inconnues est réduit au système suivant :

$$ \forall m \in \{-M \dots M\}, \forall n \in \{-N \dots N\}, $$

$$ \sum_{|p| \le M \text{ et } |m-p| \le M} \sum_{|q| \le N \text{ et } |n-q| \le N} (mp+nq\sigma^2) \lambda_{m-p,n-q} \delta\theta_{p,q} = \frac{m\lambda_{m,n}}{2i\pi P_x} \quad (70) $$

Ce nouveau système linéaire (70) comporte $(2M+1) \times (2N+1)$ équations et autant d'inconnues ; il peut donc être résolu par le biais des techniques classiques de résolution numérique des systèmes linéaires.

Une fois cette étape franchie, les deux séries filtrées $(\lambda_{m,n})_{|m| \le M, |n| \le N}$ et $(\delta\theta_{m,n})_{|m| \le M, |n| \le N}$ étant dès lors connues, la conductivité thermique effective $\Lambda$ du milieu périodique s'obtient par filtrage passé-bas de la relation (68) :

$$ \Lambda \approx \bar{\lambda}_{x=0} - 2i\pi P_x \sum_{|p| \le M \text{ et } |m-p| \le M} \sum_{|q| \le N} p \lambda_{m-p,-q} \delta\theta_{p,q} \quad (71) $$

### 3-1-4-2-6- Présentation de quelques résultats obtenus au moyen de l'analyse spectrale

Nous allons maintenant appliquer les différents résultats issus de l'analyse spectrale que nous avons énoncés dans les chapitres précédents à des géométries bidimensionnelles périodiques de complexité croissante. Ces différents tests de validation seront l'occasion de souligner les points forts mais aussi les faiblesses de cette nouvelle approche par rapport à l'approche plus classique des différences finies présentée au chapitre 2-1-.

#### 3-1-4-2-6-1- Première étape : milieux monodimensionnels périodiques

Nous avons vu au chapitre 3-1-2- que lorsque le champ de conductivité thermique $\lambda(x,y)$ est uniquement fonction de la variable $x$, il en est de même du champ de température $\theta(x,y)$ et donc du champ de déviation thermique $\delta\theta(x,y)$. Or, du point de vue de l'analyse spectrale, l'indépendance des champs $\lambda$ et $\delta\theta$ vis-à-vis de la variable $y$ se traduit par le fait que les termes $\lambda_{m,n}$ et $\delta\theta_{m,n}$ de leurs séries de Fourier bidimensionnelles respectives sont nuls dès que $n \neq 0$. Dans ces conditions, les deux relations (70) et (71) permettant dans un premier temps la détermination de la série de Fourier du champ $\delta\theta$ puis dans un deuxième temps l'évaluation de la conductivité thermique effective se simplifient et deviennent respectivement :
$$ \forall m \in \{-M...M\}, \quad \sum_{|p| \le M \text{ et } |m-p| \le M} p \lambda_{m-p,0} \delta\theta_{p,0} = \frac{\lambda_{m,0}}{2i\pi P_x} \quad (72) $$

$$ \Lambda \approx \lambda(x=0) - 2i\pi P_x \sum_{|p| \le M \text{ et } |m-p| \le M} p \lambda_{m-p,0} \delta\theta_{p,0} \quad (73) $$

Deux remarques seront faites quant aux deux formulations mathématiques ci-dessus :

1. Le système (72) est un système linéaire comportant (2M+1) équations et autant d'inconnues.

2. Attention ! Il est tentant de procéder à une injection rapide et irréfléchie de (72) dans (73) ... Il faut bien noter que la sommation double de l'expression (73) de $\Lambda$ comporte des couples (m,p) qui n'apparaissent à aucun moment dans le système (72) : ainsi par exemple le couple (m = M+1, p = 1), et de façon générale tous les couples (m,p) pour lesquels $|m| > M, |p| \le M \text{ et } |m-p| \le M$.

Nous allons maintenant tester les deux formules (72) et (73) sur le milieu étudié au chapitre 3-1-2- : pour cet échantillon, le paramètre $P_x$ vaut 3 et la période spatiale du champ de conductivité thermique dans la direction x est définie par la restriction de la fonction génératrice
$$ \lambda(x) = 10 + 9 \cos\left(6\pi \frac{x}{X}\right) \quad (\text{en W/m/K}) \text{ sur l'intervalle } 0 \le \frac{x}{X} \le \frac{1}{3}. $$

Les calculs sont extrêmement simples. La série de Fourier monodimensionnelle $(\lambda_{m,0})_{m \in \mathbb{Z}}$ du champ $\lambda$ est constituée d'uniquement trois termes non nuls : $\lambda_{\pm 1,0} = 9/2$ et $\lambda_{0,0} = 10$. Une valeur du niveau de filtrage passe-bas M étant fixée, la résolution du système linéaire (72) aboutit à la série de Fourier monodimensionnelle filtrée $(\delta\theta_{m,0})_{|m| \le M}$ du champ $\delta\theta$ ; finalement, une transformation de Fourier inverse numérique est opérée sur cette série, et ainsi on obtient le champ $\delta\theta$ sur la période spatiale $0 \le \frac{x}{X} \le \frac{1}{3}$.

En ce qui concerne la détermination du champ $\delta\theta$, l'analyse spectrale se révèle plutôt efficace. A titre d'illustration, nous nous sommes limités à un niveau de filtrage passe-bas M de 5, ce qui veut dire que seules les cinq premières harmoniques ont été retenues dans l'expression du champ $\delta\theta$. La détermination des onze inconnues $(\delta\theta_{m,0})_{|m| \le 5}$ a été instantanée, et le champ de température monodimensionnel $\theta(x)$ que nous avons obtenu est représenté à la figure 26 : comme on peut le constater, ce champ est très proche de l'expression analytique (41) déterminée au chapitre 3-1-2-.

La valeur numérique de la conductivité thermique effective $\Lambda$ fournie par la formule (73) nécessite quant à elle la prise en compte d'un plus grand nombre M d'harmoniques pour une précision acceptable. En effet, le graphe des variations de la conductivité $\Lambda$ en fonction du nombre M, tracé à la figure 27, présente un comportement oscillant et lentement convergent vers une valeur asymptotique qui se trouve être la valeur théorique de la conductivité thermique effective, soit $\sqrt{19}$. L'analyse de ce graphe montre qu'il faut atteindre un nombre M d'harmoniques de 16 - soit un nombre somme toute encore assez petit, le système linéaire (72) n'ayant que 33 inconnues - pour que l'écart relatif entre la valeur numérique de $\Lambda$ et sa valeur théorique passe en dessous de 1‰.
Áinsi donc, l'analyse spectrale s'est avérée très efficace et très peu gourmande en termes de moyens de calcul pour la résolution du problème thermique du chapitre 3-1-2. Mais il faut ici remarquer que le cas que nous venons de traiter était, du point de vue de l'analyse spectrale, particulièrement simple : en effet, le champ de conductivité thermique monodimensionnel choisi présentait un spectre de fréquences spatiales très étroit. Nous allons maintenant compliquer quelque peu le problème à résoudre en choisissant un champ de conductivité thermique certes toujours monodimensionnel et toujours fonction de la seule variable x mais cette fois présentant des variations beaucoup plus brutales en fonction de x.

Le champ de conductivité thermique $\lambda(x)$ retenu pour cette nouvelle étude est représenté à la figure 28 : il s'agit d'un champ discontinu oscillant entre les valeurs 1 W/m/K et 5 W/m/K, et périodique constitué de quatre périodes.

Là encore, les calculs sont très rapides et sans grande difficulté. La seule nouveauté apportée par cet exemple concerne l'obtention de la série de Fourier monodimensionnelle $(\lambda_{m,0})_{m \in Z}$ du champ $\lambda\square$; pour cette opération, une transformation de Fourier numérique est pratiquée sur la courbe $\lambda(x)$ restreinte à sa période spatiale. Pour le reste, le traitement est rigoureusement identique à celui décrit précédemment : une valeur du niveau de filtrage passe-bas M étant fixée, la résolution du système linéaire (72) aboutit à la série de Fourier monodimensionnelle filtrée $(\delta\theta_{m,0})_{|m| \le M}$ du champ $\delta\theta$ ; puis une transformation de Fourier inverse numérique est opérée sur cette série, laquelle transformation fournit le champ $\delta\theta$ sur la période spatiale $0 \le \frac{x}{X} \le \frac{1}{4}$.

La figure 29 compare le champ de température monodimensionnel $\theta(x)$ issu de l'analyse spectrale avec un niveau de filtrage passe-bas M de 10 à la solution analytique du problème thermique posé ; l'accord entre les deux courbes est très satisfaisant malgré un seuillage en fréquence relativement sévère : en effet, seules les dix premières harmoniques ont été retenues dans les expressions des deux champs $\lambda(x)$ et $\delta\theta(x)$.

A ce sujet, il est intéressant d'observer de plus près l'influence dudit seuillage passe-bas sur la validité du champ de déviation thermique $\delta\theta(x)$ calculé par analyse spectrale ; cette préoccupation prend un sens tout particulier dans le cadre de l'exemple étudié ici en ce sens que, comme il a été indiqué plus haut, le champ de conductivité thermique choisi présente un spectre de fréquences spatiales assez large. Le graphe de la figure 30 montre trois de ces champs de déviation thermique $\delta\theta(x)$ calculés par analyse spectrale pour trois valeurs différentes du niveau de filtrage passe-bas M ; la confrontation de ces trois courbes à la solution analytique fait clairement apparaître l'apport des harmoniques "hautes" pour une définition plus fine de la courbe $\delta\theta(x)$, ceci au prix cependant d'une quantité de calculs plus importante.

Concernant la valeur numérique de la conductivité thermique effective $\Lambda$ fournie par la formule (73), un travail identique à celui effectué plus haut, à savoir le traçage de la courbe de $\Lambda$ en fonction du niveau de filtrage passe-bas M, aboutit au graphe de la figure 31. Cette fois, la convergence de la courbe $\Lambda(M)$ vers la valeur théorique de $\Lambda\left(\text{qui est } \frac{1}{0.84} \approx 1.1905 \text{ W/m/K}\right)$ est extrêmement lente puisque pour un niveau de filtrage passe-bas M de 50 (ce qui signifie la prise en compte des 50 premières harmoniques pour chacun des deux champs $\lambda$ et $\delta\theta$), la courbe $\Lambda(M)$ présente une oscillation d'amplitude égale à $\pm 12\%$ de sa valeur asymptotique. On met ici le doigt sur une limitation importante de l'analyse spectrale par rapport à la technique classique par différences
finies, à savoir une certaine faiblesse dans la production d'une valeur correcte de conductivité thermique effective lorsque le champ de conductivité thermique présente un contenu spectral important.

### 3-1-4-2-6-2- Deuxième étape : milieu bidimensionnel périodique - champ de conductivité thermique à variations douces

Jusqu'à présent, nous avons testé l'analyse spectrale sur deux problèmes monodimensionnels ; ces deux tests nous ont permis de constater les performances intéressantes de la méthode notamment en ce qui a trait à la détermination des champs de déviation thermique $\delta\theta$. Nous allons maintenant nous compliquer quelque peu la tâche en ajoutant une dimension d'espace supplémentaire aux problèmes que nous allons tenter de résoudre. Par ailleurs, à la lumière de la conclusion énoncée au chapitre précédent quant aux difficultés éprouvées par l'analyse spectrale dans la gestion de signaux à contenus spectraux importants, nous nous limiterons dans le cadre de ce chapitre au traitement d'un exemple pour lequel le champ de conductivité thermique présente un contenu fréquentiel relativement restreint.

L'échantillon d'étude, décrit sous forme d'image à la figure 32, est un milieu bidimensionnel périodique constitué de 3 × 3 périodes spatiales. La maille génératrice de ce milieu compte 100 × 100 pixels, et la répartition des conductivités thermiques sur cette maille, également explicitée à la figure 32, montre des variations lentes et continues de la conductivité thermique entre ses valeurs minimale (1 W/m/K) et maximale (5 W/m/K).

Le champ de déviation thermique $\delta\theta$ calculé par différences finies est représenté sous forme d'image à la figure 33 ; sur cette figure est également reporté le champ $\delta\theta$ moyenné sur les neuf périodes de l'échantillon. Quant à la conductivité thermique effective $\Lambda$ fournie par notre programme de différences finies pour le milieu d'étude, les résultats du calcul sont les suivants :

$$ (\bar{\Lambda}, \varepsilon_{\Lambda}) = (1.4413 \text{ W/m/K}, 4.6 \ 10^{-8}). $$

A l'image de ce qui a été présenté plus haut lors des deux tests monodimensionnels, le traitement du problème par analyse spectrale se fait en quatre étapes :

1. Transformation de Fourier bidimensionnelle numérique du champ de conductivité thermique $\lambda$ restreint à la maille génératrice du milieu ; ainsi on obtient une partie de la série de Fourier bidimensionnelle $(\lambda_{m,n})_{m \in Z, n \in Z}$ du champ $\lambda$.

2. Choix des deux niveaux de filtrage passe-bas M et N, puis résolution du système linéaire (70) au terme de laquelle on atteint la série de Fourier bidimensionnelle filtrée $(\delta\theta_{m,n})_{|m| \le M, |n| \le N}$ du champ de déviation thermique $\delta\theta$.

3. Transformation de Fourier inverse bidimensionnelle numérique de la série filtrée $(\delta\theta_{m,n})_{|m| \le M, |n| \le N}$ ; ainsi on aboutit au champ de déviation thermique $\delta\theta$ filtré et restreint à la maille génératrice du milieu.

4. Calcul de la conductivité thermique effective $\Lambda$ à l'aide de la relation (71).

Dans le cadre de l'exemple étudié ici, la maille comporte 100 × 100 pixels ; par conséquent, l'étape 1. du calcul fournit 100 × 100 composantes spectrales correspondant aux couples de fréquences spatiales $(m\sigma_x, n\sigma_y)_{|m| \le 50, |n| \le 50}$ (on ne s'étendra pas ici sur la limitation bien connue de
l'algorithme de transformation de Fourier numérique concernant la gestion des fréquences de Nyquist). Nos calculs ont été effectués avec la valeur 10 pour les deux niveaux de filtrage passe-bas M et N, ce qui veut dire que des 100 × 100 composantes spectrales ($\lambda_{m,n}$) $|m| \le 50, |n| \le 50$, nous n'en avons retenu que 21 × 21, à savoir les $\lambda_{m,n}$ pour lesquels $|m| \le 10$ et $|n| \le 10$. L'étape 2. du calcul nous a donc fourni la série filtrée ($\delta\theta_{m,n}$) $|m| \le 10, |n| \le 10$, et nous avons complété cette série filtrée par des zéros pour construire une série ($\delta\theta_{m,n}$) $|m| \le 50, |n| \le 50$ qui soit la duale spectrale du champ de déviation thermique $\delta\theta$ restreint à la maille.

Malgré le filtrage passe-bas plutôt sévère que nous nous sommes choisis, les résultats obtenus par analyse spectrale sont en très bon accord avec ceux issus des différences finies, et ce pour une quantité de calculs beaucoup moindre (le système linéaire à résoudre passe de 300 × 300 = 90000 inconnues dans le cadre de l'approche par différences finies à seulement 21 × 21 = 441 inconnues dans le cadre de l'approche par analyse spectrale). Ainsi, la figure 34 compare les deux champs de déviation thermique $\delta\theta$ issus des deux approches ; reprenant la définition de l'écart relatif entre deux images $I_1$ et $I_2$ introduite au chapitre 3-1-4-2-3-, cet écart relatif entre les deux champs $\delta\theta$ est évalué à seulement 1.6%. Quant à la conductivité thermique effective $\Lambda$, son calcul par analyse spectrale - via la relation (71) - aboutit à la valeur de 1.4407 W/m/K, tandis que les différences finies prévoient une valeur de $\Lambda$ de 1.4413 W/m/K ...

### 3-1-4-2-6-3- Troisième étape : milieu bidimensionnel périodique - champ de conductivité thermique discontinu

Il est clair que l'exemple bidimensionnel que nous avons traité ci-dessus était parfaitement adapté à une analyse spectrale : en effet, les couvertures fréquentielles très étroites des champs de conductivité thermique $\lambda$ et de déviation thermique $\delta\theta$ autorisaient le choix de valeurs faibles pour les paramètres M et N sans perte d'information aux hautes fréquences spatiales. Nous allons mettre un terme à cette étude consacrée à l'analyse spectrale avec un cas plus délicat pour lequel le champ de conductivité thermique $\lambda$ comporte des composantes à hautes fréquences spatiales non négligeables.

Pour cela, nous reprenons l'exemple du milieu bidimensionnel périodique de type nid d'abeille déjà étudié par différences finies au chapitre 3-1-4-1-. Le champ de conductivité thermique $\lambda$ étant discontinu, son contenu spectral est large, et par conséquent il faut s'attendre à des difficultés dans le traitement de ce cas par le biais de l'analyse spectrale.

Les figures 35(a) et 35(b) rappellent, sous forme d'images, respectivement l'allure de la maille considérée et le champ de déviation thermique $\delta\theta$ moyen calculé par différences finies sur l'échantillon "10 × 10". Dans ce qui suit, nous allons essayer de reproduire ce champ $\delta\theta$ au moyen de l'analyse spectrale, d'abord avec des niveaux de filtrage passe-bas M et N fixés à 10 (figures 35(c1), 35(c2) et 35(c3)), puis avec les niveaux M et N poussés à 20 (figures 35(d1), 35(d2) et 35(d3)).

Les images 35(c1) et 35(d1) montrent les deux champs de conductivité thermique $\lambda$ résultant des deux filtrages passe-bas évoqués ci-dessus. La comparaison des deux champs $\lambda$ filtrés à la maille d'origine (image 35(a)) est éloquente : les filtrages se traduisent par l'apparition d'oscillations caractéristiques des fréquences spatiales de coupure. Ces oscillations s'accompagnent de variations importantes des champs $\lambda$ filtrés au niveau des discontinuités du champ $\lambda$ d'origine : ainsi, tandis qu'au niveau d'une discontinuité la conductivité thermique passe entre les deux valeurs 1 W/m/K et
5 W/m/K sur l'image 35(a), au même endroit la conductivité thermique passe certes de façon plus continue mais entre les deux valeurs 0.101 W/m/K et 5.567 W/m/K sur l'image 35(c1). De telles distorsions du champ λ d'origine, très classiques des opérations de filtrage passe-bas, indiquent de façon claire que les filtrages qui ont été pratiqués se sont accompagnés de pertes d'information spectrale non négligeables ; les divergences des deux champs λ filtrés par rapport au champ λ d'origine peuvent d'ailleurs se quantifier au moyen de la notion d'écart relatif entre deux images $I_1$ et $I_2$ introduite au chapitre 3-1-4-2-3- : cet écart relatif vaut 59.2% entre les images 35(a) et 35(c1) et 26.3% entre les images 35(a) et 35(d1).

Les deux champs de déviation thermique $δθ$ issus de l'analyse spectrale sont représentés sous forme d'images aux figures 35(c2) et 35(d2). Au premier coup d'œil, ces deux champs semblent relativement proches de celui calculé par différences finies (représenté sous forme d'image à la figure 35(b)) ; les figures 35(c3) et 35(d3) font certes apparaître des différences, mais le calcul des écarts relatifs entre les images 35(b) et 35(c2) d'une part et entre les images 35(b) et 35(d2) d'autre part fournit respectivement les valeurs de 9.9% et 6.5%. A ce stade, il peut paraître étonnant que, se basant sur un champ λ filtré jugé passablement distordu par rapport au champ λ d'origine, le calcul aboutisse néanmoins à un champ $δθ$ en assez bon accord avec celui calculé par différences finies ; en termes plus simples, comment l'analyse spectrale fait-elle pour construire un champ $δθ$ erroné de seulement 9.9% à partir d'un champ λ erroné à hauteur de 59.2% ? La réponse à cette question se trouve à la figure 23 : le spectre fréquentiel du champ de déviation thermique $δθ$ est beaucoup plus étroit que celui du champ de conductivité thermique λ, et ceci pour une raison toute simple liée au fait que contrairement au champ λ, le champ $δθ$ est toujours un champ continu. Par conséquent il faut retenir le résultat important suivant : même si le champ de conductivité thermique λ ne semble a priori pas favorable à un traitement par analyse spectrale, on peut quand même attendre de cette analyse spectrale le calcul d'un champ de déviation thermique $δθ$ assez voisin du résultat exact.

En revanche, le calcul des conductivités thermiques effectives $Λ$ à l'aide de la relation (71) aboutit aux deux valeurs suivantes : 2.2802 W/m/K lorsque $M = N = 10$ et 1.8943 W/m/K lorsque $M = N = 20$, soit des valeurs très éloignées de celles calculées par différences finies au chapitre 3-1-4-1-. On confirme ici la conclusion énoncée au chapitre 3-1-4-2-6-1- lors de l'étude d'un milieu monodimensionnel périodique à fonction de conductivité thermique $λ(x)$ discontinue, à savoir les grandes difficultés de l'analyse spectrale à évaluer correctement la conductivité thermique effective $Λ$ lorsque le champ de conductivité thermique λ présente un spectre fréquentiel étendu.

### 3-1-4-2-7- Ce qu'il faut retenir de cette partie consacrée à l'emploi de l'analyse spectrale pour la détermination du champ de température s'établissant au sein d'un milieu périodique

On peut résumer les principales conclusions de cette étude en quatre points :

1. Lorsque l'échantillon est périodique et comporte un grand nombre de mailles, le champ de déviation thermique $δθ$ devient un champ périodique de mêmes périodes spatiales que le champ de conductivité thermique $λ$. Il est alors tout naturel d'essayer de résoudre le problème thermique au moyen des outils de l'analyse spectrale, car cette analyse spectrale porte en elle les propriétés de périodicité spatiale des champs $λ$ et $δθ$.

2. L'analyse spectrale des champs $λ$ et $δθ$ n'opérant que sur une période spatiale, elle implique beaucoup moins de temps de calcul et de mémoire informatique que le calcul classique par différences finies.
3. Malgré la nécessité imposée par la technique d'opérer des filtrages passe-bas sur les deux champs de conductivité thermique $\lambda$ et de déviation thermique $\delta\theta$, et même lorsque le champ $\lambda$ présente un contenu fréquentiel très large, l'analyse spectrale calcule des champs $\delta\theta$ en très bon accord avec ceux évalués par différences finies. Cet heureux état de fait tient à ce que, le champ $\delta\theta$ étant toujours un champ continu, son contenu fréquentiel est par conséquent beaucoup moins large que celui du champ $\lambda$.

4. En revanche, le filtrage passe-bas du champ de conductivité thermique $\lambda$ se traduit par une certaine incapacité de l'analyse spectrale à évaluer correctement des conductivités thermiques effectives $\Lambda$, surtout lorsque le champ $\lambda$ présente un spectre fréquentiel étendu.

Dès lors, les points forts et faiblesses de l'analyse spectrale étant bien cernés, il est peut-être envisageable, lorsque l'on a affaire à des milieux périodiques, d'utiliser de façon complémentaire cette technique avec celle des différences finies en exploitant au maximum les atouts des deux approches. Ainsi, on pourrait commencer par calculer au moyen de l'analyse spectrale le champ de déviation thermique $\delta\theta$ ; le résultat de ce calcul très rapide serait ensuite injecté dans le programme de différences finies en tant que champ à l'itération 0 de l'algorithme de Fletcher et Reeves ; la convergence dudit algorithme serait alors probablement très rapide car le champ à l'itération 0 serait déjà très proche du champ solution ; quant au calcul de la conductivité thermique effective $\Lambda$, il serait effectué par différences finies. Cette procédure mixte couplant analyse spectrale et différences finies pourrait bien réduire sensiblement les temps de calcul tout en assurant une bonne qualité des résultats fournis.

### 3-1-5- Milieu à deux constituants et anisotrope

Dans ce chapitre, nous allons nous intéresser à un aspect qui va souvent de pair avec l'hétérogénéité, à savoir l'anisotropie. Cette anisotropie est généralement produite en organisant les hétérogénéités dans le but d'optimiser les propriétés physiques dans des directions particulières. Ainsi, à titre d'illustration, un matériau composite unidirectionnel constitué de fibres bonnes conductrices de la chaleur, orientées toutes dans la même direction et noyées dans une matrice à faible conductivité thermique présentera une conductivité thermique effective dans la direction des fibres nettement supérieure à celle dans n'importe quelle direction orthogonale.

C'est ce genre de "cas d'école" que nous allons étudier dans ce qui suit. Nous allons créer informatiquement un échantillon composite unidirectionnel ; puis nous ferons subir à cet échantillon des rotations d'angles $\alpha$ divers, et nous calculerons par différences finies la conductivité thermique effective $\Lambda$ de l'échantillon orienté résultant en fonction de l'angle $\alpha$. Nous ferons ensuite appel à la théorie classique de la conduction thermique dans les milieux anisotropes, et nous comparerons les indicatrices de conductivité thermique fournies par les deux approches.

#### 3-1-5-1- Calcul par différences finies de la conductivité thermique effective $\Lambda$ en fonction de l'angle $\alpha$ d'orientation du milieu hétérogène

Le milieu composite unidirectionnel d'origine est composé de deux constituants : des pixels noirs, matérialisant les fibres, de conductivité thermique choisie égale à 5 W/m/K, et des pixels blancs, matérialisant la matrice, de conductivité thermique choisie égale à 1 W/m/K. Les fibres ont une épaisseur de 3 pixels (noirs) et sont orientées dans la direction verticale ; par ailleurs, deux fibres voisines sont distantes de 7 pixels (blancs), de sorte que la fraction surfacique de fibres dans le milieu $\eta_r$ est égale à 0.3.
Nous avons fait subir au milieu d'origine des rotations d'angles 0°, 15°, 30°, 45°, 60°, 75° et 90°, et de chacun des sept milieux résultants nous avons extrait un sous-ensemble de 400 × 400 pixels. Les sept échantillons ainsi produits sont représentés sous forme d'images à l'échelle ½ à la figure 36 (seul l'échantillon de la figure 36(c) est reproduit en grandeur réelle à la figure 8).

Nous avons lancé notre programme de différences finies sur les sept géométries décrites ci-dessus ; les résultats des calculs sont consignés dans le tableau qui suit.

<table><caption>Tableau 3</caption><thead><tr><td>angle de rotation α</td><td>Λ̄ (W/m/K)</td><td>εΛ</td></tr></thead><tbody><tr><td>0°</td><td>2.2000</td><td>6.3 10<sup>-14</sup></td></tr><tr><td>15°</td><td>2.0232</td><td>1.5 10<sup>-8</sup></td></tr><tr><td>30°</td><td>1.8198</td><td>2.7 10<sup>-8</sup></td></tr><tr><td>45°</td><td>1.6340</td><td>1.8 10<sup>-7</sup></td></tr><tr><td>60°</td><td>1.4676</td><td>5.0 10<sup>-8</sup></td></tr><tr><td>75°</td><td>1.3642</td><td>4.5 10<sup>-8</sup></td></tr><tr><td>90°</td><td>1.3158</td><td>1.6 10<sup>-8</sup></td></tr></tbody></table>

La colonne de valeurs de $\bar{\Lambda}$ du tableau ci-dessus illustre bien la notion d'anisotropie introduite plus haut : d'une part, le milieu hétérogène apparaît le plus conducteur thermique lorsqu'il est sollicité thermiquement dans la direction de ses fibres, et d'autre part, c'est dans la direction orthogonale aux fibres qu'il apparaît le plus isolant thermique. Par ailleurs, les deux valeurs extrêmes de $\Lambda$ peuvent facilement être calculées de manière théorique : en effet, lorsque l'angle de rotation $\alpha$ vaut 0° et 90°, l'échantillon apparaît comme étant constitué de résistances thermiques disposées respectivement en parallèle et en série ; les valeurs théoriques de $\Lambda(\alpha = 0°)$ et de $\Lambda(\alpha = 90°)$ sont donc respectivement (avec les notations évidentes) $\Lambda(\alpha = 0°) = \eta_f \lambda_f + \eta_m \lambda_m = 2.2000 \text{ W/m/K}$ et $\Lambda(\alpha = 90°) = \left(\frac{\eta_f}{\lambda_f} + \frac{\eta_m}{\lambda_m}\right)^{-1} \approx 1.3158 \text{ W/m/K}$, soit exactement les deux valeurs déterminées par notre programme de différences finies.

### 3-1-5-2- Recherche par voie analytique de la conductivité thermique effective $\Lambda$ du milieu homogène anisotrope équivalent

Les deux valeurs extrêmes $\Lambda(\alpha = 0°)$ et $\Lambda(\alpha = 90°)$ sont donc faciles à évaluer analytiquement ; mais qu'en est-il du reste de la courbe de la conductivité thermique effective $\Lambda$ en fonction de l'angle de rotation $\alpha$ ? Pour répondre à cette question, nous allons étudier comment la courbe $\Lambda(\alpha)$ que nous venons de déterminer à l'aide de nos calculs par différences finies se compare à celle que l'on peut déduire d'un raisonnement simple d'homogénéisation thermique.

Le raisonnement est le suivant : si on veut remplacer le milieu hétérogène réel par un milieu homogène équivalent, alors on se doit de donner au milieu homogène équivalent des propriétés de conduction thermique anisotropes. Le milieu homogène équivalent doit donc être défini, pour ce qui
a trait à son anisotropie thermique, d'une part par deux directions principales xp et yp orthogonales et d'autre part par deux conductivités thermiques principales $\lambda_{xp}$ et $\lambda_{yp}$ correspondant aux conductivités thermiques effectives dans les deux directions principales. Or ces deux directions principales xp et yp sont définies par la propriété caractéristique suivante : si le gradient de température est porté par l'une de ces deux directions, alors le flux thermique aussi ; dès lors, les deux directions principales xp et yp s'imposent comme étant respectivement la direction des fibres et la direction orthogonale aux fibres. Par conséquent, le milieu homogène équivalent au milieu hétérogène réel en termes de propriétés de conduction thermique est défini par (i) deux directions principales xp et yp orthogonales, égales respectivement à la direction des fibres et à la direction orthogonale aux fibres, et (ii) deux conductivités thermiques principales $\lambda_{xp}$ et $\lambda_{yp}$ valant respectivement $\eta_f \lambda_f + \eta_m \lambda_m = 2.2000 \text{ W/m/K}$ et $(\frac{\eta_f}{\lambda_f} + \frac{\eta_m}{\lambda_m})^{-1} \approx 1.3158 \text{ W/m/K}$. Le tenseur de conductivité thermique du milieu homogène équivalent dans ses axes principaux est un tenseur diagonal valant $\begin{pmatrix} \lambda_{xp} & 0 \\ 0 & \lambda_{yp} \end{pmatrix}$.

Faisons maintenant subir au milieu homogène anisotrope équivalent une rotation d'angle $\alpha$ ; le tenseur de conductivité thermique dudit milieu dans le système d'axes (x,y) non principal (voir notations sur la figure 37) perd son caractère diagonal et devient du type $\begin{pmatrix} \lambda_{xx} & \lambda_{xy} \\ \lambda_{xy} & \lambda_{yy} \end{pmatrix}$. Un calcul très classique de changement de repère permet d'exprimer les trois conductivités thermiques $\lambda_{xx}$, $\lambda_{xy}$ et $\lambda_{yy}$ en fonction des deux conductivités thermiques principales $\lambda_{xp}$ et $\lambda_{yp}$ ; les formules de transformation sont les suivantes :

$$ egin{cases} 
\lambda_{xx} = \lambda_{xp} \cos^2 \alpha + \lambda_{yp} \sin^2 \alpha \\ 
\lambda_{xy} = (\lambda_{xp} - \lambda_{yp}) \cos \alpha \sin \alpha \\ 
\lambda_{yy} = \lambda_{xp} \sin^2 \alpha + \lambda_{yp} \cos^2 \alpha 
\end{cases} 
\quad (74) $$

Revenons maintenant au problème thermique (11) que nous avons à résoudre, et supposons dans un premier temps que notre échantillon homogène anisotrope équivalent est d'extension infinie selon y. Dans ce cas, les conditions d'adiabaticité des deux faces de l'échantillon autres que celles à température imposée disparaissent (car ces deux faces sont repoussées à l'infini), de sorte que le système (11) devient :

$$ egin{cases} 
\bar{\nabla} \cdot \left( - \begin{pmatrix} \lambda_{xx} & \lambda_{xy} \\ \lambda_{xy} & \lambda_{yy} \end{pmatrix} \bar{\nabla} \theta \right) = 0 \quad \Leftrightarrow \quad \lambda_{xx} \frac{\partial^2 \theta}{\partial x^2} + 2 \lambda_{xy} \frac{\partial^2 \theta}{\partial x \partial y} + \lambda_{yy} \frac{\partial^2 \theta}{\partial y^2} = 0 \\ 
\theta = 1 \text{ en } x = 0 \\ 
\theta = 0 \text{ en } x = X 
\end{cases} 
\quad \text{dans le volume de l'échantillon} 
\quad (75) $$
Ce nouveau système (75) est extrêmement simple à résoudre : son champ de température solution $\theta(x,y)$ est en fait fonction de la seule variable $x$ et vaut $\theta(x) = 1 - \frac{x}{X}$. Le calcul du flux thermique $\bar{\psi} = -\begin{pmatrix} \lambda_{xx} & \lambda_{xy} \\ \lambda_{xy} & \lambda_{yy} \end{pmatrix} \bar{\nabla} \theta$ est alors immédiat : ce flux thermique $\bar{\psi}$ est le même en tous les points de l'échantillon et vaut $\bar{\psi} = \frac{1}{X} (\lambda_{xx} \bar{x} + \lambda_{xy} \bar{y})$. Notons que le vecteur $\bar{\psi}$ n'est pas colinéaire au vecteur $\bar{\nabla} \theta$, conséquence directe de l'anisotropie thermique du milieu homogène équivalent et du fait que $\bar{\nabla} \theta$ n'est porté par aucune des deux directions principales xp et yp ; $\bar{\psi}$ est donc incliné par rapport à l'axe des $x$ d'un certain angle $\beta$ (voir figure 38) dont la tangente vaut, en introduisant le paramètre $L = \frac{\lambda_{xp}}{\lambda_{yp}}$ mesurant le degré d'anisotropie thermique du milieu homogène équivalent :

$$ \operatorname{tg}\beta = \frac{\lambda_{xy}}{\lambda_{xx}} = \frac{(\lambda_{xp} - \lambda_{yp}) \cos\alpha \sin\alpha}{\lambda_{xp} \cos^2\alpha + \lambda_{yp} \sin^2\alpha} = \frac{(L-1) \operatorname{tg}\alpha}{L + \operatorname{tg}^2\alpha} \quad (76) $$

Considérons maintenant un tuyau de flux thermique tel que celui schématisé à la figure 38, partant d'une surface S quelconque située sur la face chaude de l'échantillon et engendré par translation de cette surface S selon le vecteur $\bar{\psi}$. Invoquant des arguments en tout point similaires à ceux exposés dans l'introduction de ce document, on arrive très simplement à prouver que la puissance thermique $\Psi$ traversant une section S' du tuyau est indépendante de S' et vaut $\Psi = \iint_S \bar{\psi} \cdot \bar{x} \, dS = \psi_x S = \frac{\lambda_{xx} S}{X}$.

Mais ici s'arrête le parallèle que l'on peut établir avec le raisonnement tenu au chapitre 1- : en effet, le tuyau de flux thermique dans le cas présent n'est pas un cylindre droit, de sorte que l'égalité $\Psi = C(\theta_c - \theta_f) = \frac{\Lambda S}{X}$ écrite en (8) n'est pas valide ici. Il faut donc se définir deux nouvelles sections planes du tuyau normales au vecteur $\bar{\psi}$ et reprendre le raisonnement en prenant bien garde au fait que, ces deux nouvelles sections planes n'étant plus normales au vecteur $\bar{x}$, ce ne sont donc plus des surfaces isothermes. La figure 39 tente de rendre la situation plus explicite : elle montre, en géométrie bidimensionnelle, la surface S (qui devient un segment de longueur d), le tuyau de flux thermique, les deux sections planes du tuyau normales au vecteur $\bar{\psi}$ (qui deviennent les segments [AB] et [A'B']) et le découpage de la résistance thermique d'intérêt en résistances thermiques élémentaires qui est opéré pour les développements mathématiques qui suivent :
$$
\begin{align*}
\Psi = \frac{\lambda_{xx} \, d}{X} &= \int_{M \in [AB]} d\Psi(M) \\
&= \int_{M \in [AB]} dC(M) [\theta(M) - \theta(M')] \\
&= \int_{M \in [AB]} \frac{\Lambda(\beta) \, dz}{MM'} [\theta(M) - \theta(M')] \tag{77}
\end{align*}
$$

A ce stade il est important de souligner que l'expression de la conductance thermique élémentaire
dC(M) = Λ(β)dz / MM' écrite ci-dessus fait intervenir la conductivité thermique effective Λ du milieu
homogène anisotrope équivalent dans la direction β qui est la direction du transport
énergétique. Le développement mathématique (77) peut être prolongé ; remarquant que
θ(M) - θ(M') = 1/X M M' ⃗ x = MM' cosβ / X , on aboutit très facilement à :

$$
\Psi = \frac{\lambda_{xx} \, d}{X} = \frac{\Lambda(\beta) \, d \cos^2 \beta}{X} \qquad (78)
$$

soit :

$$
\begin{align*}
\Lambda(\beta) &= \frac{\lambda_{xx}}{\cos^2 \beta} \\
&= \lambda_{xx} (1 + \operatorname{tg}^2 \beta) \\
&= \lambda_{yp} (L \cos^2 \alpha + \sin^2 \alpha) \left[ 1 + \left( \frac{(L-1) \operatorname{tg}\alpha}{L + \operatorname{tg}^2 \alpha} \right)^2 \right] \tag{79} \\
&= \lambda_{yp} \frac{L^2 + \operatorname{tg}^2 \alpha}{L + \operatorname{tg}^2 \alpha}
\end{align*}
$$

Il importe maintenant, pour pouvoir exploiter la relation Λ(β) = λ_yp (L² + tg²α) / (L + tg²α) obtenue ci-dessus, d'étudier les propriétés de la fonction β ↦→ α définie implicitement par la relation (76). Or il s'avère que la fonction α ∈ [0; π/2] ↦→ β = arctg[(L - 1) tgα / (L + tg²α)] n'est pas bijective, comme le montre le tableau des variations de cette fonction :
<figure><img src="image_3.png" /></figure>

Ce tableau apporte les deux enseignements suivants :

1. Quelle que soit l'orientation α de l'anisotropie du milieu homogène équivalent, l'orientation β du transport énergétique au sein dudit milieu reste dans un intervalle qui peut être assez étroit : ainsi, pour le cas numérique qui nous intéresse, L vaut 1.672 et le calcul de la valeur maximale de β donne βmax ≈ 14.57°.

2. Deux valeurs différentes de l'angle α judicieusement choisies de part et d'autre de l'angle α produisant l'angle β maximal peuvent aboutir à deux valeurs identiques de l'angle β.

Le point 2. ci-dessus amène à considérer la courbe $β \mapsto Λ(β) = λ_{yp} \frac{L^2 + tg^2α}{L + tg^2α}$ comme une courbe paramétrée par l'angle α : on se choisit des valeurs de α entre 0° et 90°, on calcule les valeurs de β et de Λ correspondantes et ainsi on trace la courbe $β \mapsto Λ$. De façon très similaire, le tableau 3 qu'on a empli des résultats des calculs par différences finies au chapitre précédent doit être complété d'une colonne supplémentaire contenant les valeurs des angles de transport énergétique β :

<table><caption>Tableau 4</caption><thead><tr><td>angle de rotation α</td><td>angle de transport énergétique β</td><td>Λ̄ (W/m/K)</td><td>εΛ</td></tr></thead><tbody><tr><td>0°</td><td>0.00°</td><td>2.2000</td><td>6.3 10<sup>-14</sup></td></tr><tr><td>15°</td><td>5.90°</td><td>2.0232</td><td>1.5 10<sup>-8</sup></td></tr><tr><td>30°</td><td>10.95°</td><td>1.8198</td><td>2.7 10<sup>-8</sup></td></tr><tr><td>45°</td><td>14.12°</td><td>1.6340</td><td>1.8 10<sup>-7</sup></td></tr><tr><td>60°</td><td>13.99°</td><td>1.4676</td><td>5.0 10<sup>-8</sup></td></tr><tr><td>75°</td><td>9.13°</td><td>1.3642</td><td>4.5 10<sup>-8</sup></td></tr><tr><td>90°</td><td>0.00°</td><td>1.3158</td><td>1.6 10<sup>-8</sup></td></tr></tbody></table>

La figure 40 compare le faisceau des sept points $(β, Λ̄)$ "expérimentaux" issus du tableau 4 ci-dessus à la courbe théorique $β \mapsto Λ(β) = λ_{yp} \frac{L^2 + tg^2α}{L + tg^2α}$ prévue par le raisonnement d'homogénéisation que nous avons développé dans ce chapitre ; comme on peut le constater, l'accord n'est pas très bon ...
### 3-1-5-3- Peut-on faire coïncider les résultats des approches numérique et analytique ?

Pourquoi les deux ensembles de points de la figure 40 ne s'accordent-ils pas ?

On peut invoquer un argument lié à la nature bornée du domaine de calcul. En effet, les développements analytiques présentés au chapitre précédent ont été faits sur la base d'une hypothèse de départ pour le moins exigeante, à savoir celle d'un échantillon d'extension infinie selon la direction y. Ce n'est évidemment pas le cas des échantillons que nous traitons par différences finies au chapitre 3-1-5-1-. Si dès lors on abandonne cette hypothèse simplificatrice, le système différentiel (11) à résoudre n'est plus équivalent au système (75), le champ de température solution $\theta(x,y)$ n'est plus fonction de la seule variable x et ce champ ne vaut plus $\theta(x) = 1 - \frac{x}{X}$.

Afin de mieux mesurer l'impact des frontières verticales de l'échantillon sur le champ de température solution, nous nous sommes livrés à de petits calculs simples à l'aide du code par éléments finis Castem 2000. Nous avons "construit" dans ce code un échantillon bidimensionnel carré, homogène et anisotrope, caractérisé par (i) deux directions principales xp et yp orthogonales et tournées d'un angle $\alpha = 30°$ par rapport au référentiel du carré et (ii) deux conductivités thermiques principales $\lambda_{xp}$ et $\lambda_{yp}$ valant respectivement 2.2000 W/m/K et 1.3158 W/m/K ; nous avons donc reproduit la version homogène anisotrope équivalente à l'échantillon de la figure 36(c) selon la logique développée au chapitre 3-1-5-2- précédent. Sur ce domaine de calcul, Castem 2000 a résolu l'équation de la chaleur accompagnée (i) des deux conditions aux limites de température imposée sur les deux frontières horizontales et (ii) des deux conditions aux limites d'adiabaticité sur les deux frontières verticales. La figure 41 présente sous forme d'images le champ de température $\theta$ déterminé par Castem 2000 et le champ de déviation thermique $\delta\theta$ résultant.

Comme on peut le constater à la figure 41(a), les isothermes ne sont pas des droites horizontales ; elles arrivent sur les deux frontières verticales adiabatiques avec un angle d'inclinaison qui est le même pour toutes les isothermes et qui est tel que le flux thermique $\bar{\psi}$ sur ces deux frontières est uniquement porté par l'axe $\bar{x}$. Or, ce flux thermique $\bar{\psi}$ s'exprimant par $\bar{\psi} = - \begin{pmatrix} \lambda_{xx} & \lambda_{xy} \\ \lambda_{xy} & \lambda_{yy} \end{pmatrix} \bar{\nabla} \theta$, l'écriture de la condition $\psi_y = 0$ amène à la conclusion que le vecteur $\bar{\nabla} \theta$ est incliné d'un certain angle polaire $\gamma$ par rapport à l'axe $\bar{x}$ (voir figure 41(a)), la tangente de cet angle $\gamma$ valant :

$$ \operatorname{tg}\gamma = \frac{\frac{\partial \theta}{\partial y}}{\frac{\partial \theta}{\partial x}} = - \frac{\lambda_{xy}}{\lambda_{yy}} \qquad (80) $$

Le calcul de cette quantité dans le cas que nous traitons ici (c'est-à-dire avec $\lambda_{xp} = 2.2000$ W/m/K, $\lambda_{yp} = 1.3158$ W/m/K et $\alpha = 30°$) donne $\operatorname{tg} \gamma = -0.2491$, soit $\gamma \approx 166°$.

Les deux frontières verticales adiabatiques étant, dans le cas présent, assez proches l'une de l'autre, il s'ensuit que la perturbation imposée au champ de température $\theta$ par leur présence est relativement importante. La figure 41(b), qui montre le champ de déviation thermique $\delta\theta$ sous forme d'image, confirme cette affirmation : la zone verte de cette image située autour de l'isovaleur $\delta\theta = 0$, qui représente la partie du domaine de calcul à l'intérieur de laquelle le champ $\theta(x,y)$ vaut
approximativement $1 - \frac{x}{X}$, est très étroite car vite "grignotée" par les effets des bords. Signalons
deux points supplémentaires concernant cette figure 41(b) :

1. L'isovaleur $\delta\theta = 0$ est quasiment rectiligne et partage l'échantillon en deux parties symétriques ; son orientation au centre du carré est celle du transport thermique, donc celle du vecteur $\vec{\psi}$, au point le moins influencé par les effets des bords.

2. L'anisotropie de l'échantillon se traduit par un surrefroidissement au niveau de sa frontière verticale gauche et par un suréchauffement au niveau de sa frontière verticale droite. Ces deux effets découlent de l'orientation du vecteur flux thermique $\vec{\psi}$ dont on vient de parler au point 1. ci-dessus : ce vecteur $\vec{\psi}$ présentant une petite composante orientée vers la droite, il s'ensuit qu'un excédent d'énergie thermique se trouve apporté à la frontière verticale droite de l'échantillon et qu'un excédent d'énergie thermique se trouve arraché à sa frontière verticale gauche.

L'échantillon homogène anisotrope carré donnant lieu à un champ de température $\theta(x,y)$ trop éloigné du champ $\theta(x) = 1 - \frac{x}{X}$ recherché, nous avons lancé de nouveaux calculs avec Castem 2000 mais cette fois en prenant pour domaine de calcul un rectangle obtenu à partir du carré précédent en divisant par 2 sa hauteur et en multipliant par 2 sa largeur. La figure 42 présente sous forme d'images le champ de température $\theta$ déterminé par Castem 2000 et le champ de déviation thermique $\delta\theta$ résultant. Cette fois, les deux images de la figure 42 montrent clairement que le champ $\theta$ régnant dans l'échantillon rectangle est très voisin du champ monodimensionnel $\theta(x) = 1 - \frac{x}{X}$ sur une grande majorité du domaine de calcul ; les deux frontières verticales adiabatiques induisent ici une perturbation du champ $\theta$ uniquement dans leurs voisinages immédiats.

A la lumière de ce résultat, nous avons décidé de relancer un calcul par différences finies sur notre échantillon hétérogène anisotrope à deux constituants de la figure 36(c), mais en en modifiant au préalable la géométrie : l'échantillon d'origine a vu sa hauteur divisée par 2 et sa largeur multipliée par 2, de sorte qu'il est passé d'une forme carrée avec $400 \times 400$ pixels à une forme rectangulaire avec $200 \times 800$ pixels. La figure 43(a) présente sous forme d'image la géométrie du nouvel échantillon, et les figures 43(b) et 43(c) montrent respectivement le champ de température $\theta$ s'établissant dans l'échantillon et le champ de déviation thermique $\delta\theta$ résultant.

Ces deux dernières images ressemblent beaucoup à leurs homologues des figures 42(a) et 42(b) obtenues sur l'échantillon rectangulaire homogène anisotrope équivalent ; en particulier, il apparaît clairement à la figure 43(c) que le champ bidimensionnel $\theta(x,y)$ est très proche du champ monodimensionnel $\theta(x) = 1 - \frac{x}{X}$ recherché en tous les points du domaine de calcul sauf aux voisinages immédiats des deux frontières verticales adiabatiques. Il est donc maintenant possible de se définir un tuyau de flux thermique de section suffisamment grande (c'est-à-dire, dans le cas présent d'une géométrie bidimensionnelle, de section contenant suffisamment de pixels dans la direction $\vec{y}$), reliant la face chaude de l'échantillon à la face froide de celui-ci, et incliné de l'angle de transport énergétique $\beta$ (qui dans le cas présent vaut environ $11^\circ$ - voir formule (76) avec $L = 1.672$ et $\alpha = 30^\circ$). La figure 43(d) présente, superposé au champ $\delta\theta$, le tuyau de flux thermique que nous avons retenu. Sur ce nouveau domaine, qui couvre environ la moitié du domaine rectangulaire
de calcul, le champ de température θ(x,y) peut être considéré égal à $1 - \frac{x}{X}$ ; par conséquent, si l'équivalence entre le milieu hétérogène et le milieu homogène est avérée, alors ce tuyau de flux thermique doit permettre l'évaluation de la conductivité thermique $\Lambda(\beta)$ dans la direction de transport énergétique $\beta$. Malheureusement, les calculs aboutissent à une moyenne $\bar{\Lambda}$ de la conductivité thermique effective sur le tuyau de flux thermique valant 1.8511 W/m/K ; on est donc encore bien loin de la quantité $\Lambda(\beta)$ attendue puisque ladite quantité $\Lambda(\beta)$, dont l'expression est donnée en (79), vaut 2.0530 W/m/K ...

Cette persistance de désaccord entre les propriétés thermiques effectives du milieu hétérogène anisotrope et celles du milieu homogène anisotrope réputé équivalent force à la réflexion quant à cette supposée équivalence entre les deux milieux : en d'autres termes, un milieu hétérogène anisotrope tel que celui sur lequel nous nous sommes concentrés tout au long de ce chapitre 3-1-5- peut-il, en ce qui a trait au transfert thermique homogénéisé, être remplacé par un simple tenseur de conductivité thermique anisotrope équivalent ? Les résultats de l'étude que nous venons de présenter semblent indiquer que non. Au stade où nous en sommes maintenant, il serait très intéressant de poursuivre cette réflexion dans deux directions :

1. Par des expériences : élaborer un matériau hétérogène anisotrope du genre de celui que nous avons retenu pour nos calculs (par exemple un composite 2-D carbone/résine) puis mesurer l'indicatrice de conductivité thermique effective de ce matériau à partir d'échantillons taillés selon différents angles permettrait de confronter cette indicatrice expérimentale à celle prévue par notre programme de différences finies.

2. Par la comparaison de nos résultats numériques à ceux issus d'autres approches d'homogénéisation thermique comme par exemple celle développée actuellement par nos collègues au Laboratoire des Composites Thermo-Structuraux (LCTS, UMR 5801 CNRS / CEA / SNECMA / Université de Bordeaux I).

### 3-1-6- Quelques milieux particuliers à deux constituants

Dans cette partie, nous allons exposer différents résultats fournis par notre programme de calcul sur des matériaux bidimensionnels à deux constituants de géométries diverses. Le but de ce "catalogue de résultats" est double : montrer la capacité du programme à gérer une variété importante de microstructures, mais aussi profiter de cette variété importante de microstructures pour souligner les difficultés auxquelles l'on se heurte lorsque l'on veut faire une modélisation fine des transferts thermiques dans ces architectures hétérogènes.

### 3-1-6-1- Matériaux fissurés

Dans cette partie, nous étudierons successivement deux configurations : d'abord le cas d'une fissure isolée, et ensuite le cas d'une distribution de fissures.

#### 3-1-6-1-1- Fissure isolée
Considérons un échantillon carré, homogène et isotrope, contenant en son sein une fissure isolée de taille importante. Notre échantillon d'étude est représenté sous forme d'image à la figure 44(a) : il comporte 200 × 200 pixels, et le segment noir épais matérialise la fissure. Pour les besoins du calcul, nous avons attribué aux pixels blancs de l'image une conductivité thermique locale de 1 W/m/K (il s'agit donc ici de la conductivité thermique du matériau homogène isotrope) ; quant aux pixels noirs matérialisant la fissure, postulant que celle-ci était occupée par de l'air à pression et température ambiantes, nous leur avons attribué une conductivité thermique locale de 0.025 W/m/K, c'est-à-dire la conductivité thermique de l'air dans les conditions normales de température et de pression.

Les figures 44(b) et 44(c) montrent les champs de température $\theta$ et de déviation thermique $\delta\theta$ obtenus sur l'échantillon fissuré de la figure 44(a). Visiblement, le champ $\theta$ est très perturbé par la présence de la fissure : il est sensiblement éloigné du champ monodimensionnel $\theta(x) = 1 - \frac{x}{X}$ qui règnerait dans un échantillon homogène isotrope, et par conséquent le volume d'échantillon qui est pris en compte dans le cadre de ces calculs ne peut en aucun cas donner lieu à une homogénéisation thermique satisfaisante. On pourra néanmoins noter, à l'examen de la figure 44(c), un point intéressant qui n'est pas sans rappeler une remarque émise au paragraphe 3-1-5-3- précédent lors de l'étude du milieu homogène anisotrope : la présence de la fissure se traduit par un suréchauffement de l'échantillon au voisinage "supérieur" de la fissure (c'est-à-dire du côté de la fissure orienté vers la face chaude de l'échantillon) et par un surrefroidissement de l'échantillon au voisinage "inférieur" de la fissure (c'est-à-dire du côté de la fissure orienté vers la face froide de l'échantillon). Ce suréchauffement et ce surrefroidissement s'expliquent aisément : la fissure constitue une barrière thermique intense (car sa conductivité thermique est très faible par rapport à celle de l'échantillon) qui bouleverse considérablement l'écoulement normal de l'énergie calorifique.

La barrière thermique que nous venons d'évoquer perturbe donc sensiblement le champ $\theta$ ; quant à la localisation précise de cette barrière et la détermination de son extension, il faut aller chercher ces deux informations non pas dans le champ $\theta$ mais dans le champ de la norme de son gradient $\|\vec{\nabla}\theta\|_2$ (la notation $\|\dots\|_2$ désignant la norme euclidienne classique). En effet, le transfert thermique dans nos domaines de calcul se faisant à ordre de grandeur constant pour le vecteur flux thermique $\vec{\psi} = -\lambda \vec{\nabla}\theta$, le passage d'un point A de conductivité thermique locale $\lambda_A$ à un point B de conductivité thermique locale $\lambda_B \ll \lambda_A$ se traduit immédiatement par $(\vec{\nabla}\theta)_B \gg (\vec{\nabla}\theta)_A$. Autrement dit, le champ $\|\vec{\nabla}\theta\|_2$ devrait rendre compte de l'existence et de l'intensité d'éventuelles barrières thermiques au sein du domaine de calcul. Connaissant les températures $[\theta_{i,j}]_{1\le i\le 1, 1\le j\le 1}$ aux différents pixels de l'image, ce champ $\|\vec{\nabla}\theta\|_2$ se calcule aisément - au facteur multiplicatif $1/\delta$ près dans lequel $\delta$ désigne la distance entre deux pixels voisins - à l'aide des deux formules classiques de différences finies qui suivent :

$$ \forall (i,j) \in \{1 \dots I-1\} \times \{1 \dots J-1\}, \quad \begin{cases} \left(\frac{\partial \theta}{\partial x}\right)_{i+\frac{1}{2}, j+\frac{1}{2}} \approx \frac{\theta_{i+1,j+1} + \theta_{i+1,j} - \theta_{i,j+1} - \theta_{i,j}}{2\delta} \\ \\ \left(\frac{\partial \theta}{\partial y}\right)_{i+\frac{1}{2}, j+\frac{1}{2}} \approx \frac{\theta_{i+1,j+1} + \theta_{i,j+1} - \theta_{i+1,j} - \theta_{i,j}}{2\delta} \end{cases} \quad (81) $$
Nous avons calculé le champ $||\vec{\nabla}\Theta||_2$ régnant dans l'échantillon fissuré au moyen des formules (81) (dans lesquelles la distance $\delta$ a été remplacée par 1), puis nous avons présenté ce champ $||\vec{\nabla}\Theta||_2$ sous forme d'image à la figure 44(d). Comme prévu, il révèle de façon très efficace l'existence et l'extension de la résistance thermique constituée par la fissure ; l'outil "champ $||\vec{\nabla}\Theta||_2$" sera donc par la suite régulièrement mis à profit pour la détermination et la localisation de telles barrières dans des situations de transfert thermique plus complexes que celle étudiée ici.

### 3-1-6-1-2- Distribution de fissures

A présent, changeons de matériau : nous allons maintenant considérer un échantillon carré, homogène et isotrope, contenant en son sein une distribution de fissures. Notre échantillon d'étude est représenté sous forme d'image à la figure 45(a) : il comporte 400 x 400 pixels qui sont soit blancs (matériau homogène isotrope, donc conductivités thermiques locales valant 1 W/m/K) soit noirs (fissures occupées par de l'air, donc conductivités thermiques locales valant 0.025 W/m/K). L'échantillon comporte 100 fissures aléatoires qui ont été produites de la façon suivante : une fissure étant complètement déterminée par (i) les coordonnées $i \in \{1... 400\}$ et $j \in \{1... 400\}$ de son centre, (ii) son angle polaire $\alpha \in [0, \pi]$ et (iii) sa longueur totale $\ell \ge 0$, les 100 quadruplets $(i, j, \alpha, \ell)$ ont été obtenus par 100 réalisations du quadruplet de variables aléatoires $(\hat{i}, \hat{j}, \hat{\alpha}, \hat{\ell})$ dans lequel (i) $\hat{i}$ et $\hat{j}$ ont été définies comme étant toutes deux la loi uniforme entre 1 et 400, (ii) $\hat{\alpha}$ a été définie comme étant la loi uniforme entre 0 et $\pi$, et (iii) $\hat{\ell}$ a été définie comme étant la loi uniforme entre 0 et 30 (ce choix de la loi $\hat{\ell}$ est arbitraire et donc discutable).

Les figures 45(b) et 45(c) présentent sous forme d'images les champs $\Theta$ et $\delta\Theta$ résultant du calcul par différences finies. A l'examen de ces deux images, on constate que les fissures apportent certes de petites perturbations locales au champ $\Theta$, mais que celui-ci ne s'éloigne jamais beaucoup du champ monodimensionnel homogénéisé $\Theta(x) = 1 - \frac{x}{X}$ (les valeurs couvertes par le champ $\delta\Theta$ sont, en valeur absolue, inférieures à 3%). Le présent échantillon peuplé de sa distribution de fissures semble donc homogénéisable du point de vue de la conduction thermique, et dès lors on peut s'intéresser à la valeur de la conductivité thermique effective moyenne $\bar{\Lambda}$ fournie par notre programme de calcul : 0.9404 W/m/K, soit un peu moins que la conductivité thermique du matériau homogène isotrope non fissuré (qui, rappelons-le, a été prise égale à 1 W/m/K).

La figure 45(d) présente sous forme d'image le champ $||\vec{\nabla}\Theta||_2$ résultant du champ $\Theta$. Comme on l'a expliqué plus haut, ce champ est un détecteur des barrières thermiques s'opposant à l'écoulement normal de l'énergie calorifique. A l'examen de la figure 45(d), deux informations nouvelles apparaissent concernant le champ $||\vec{\nabla}\Theta||_2$: (i) ce champ ne "voit" pas les barrières thermiques lorsque celles-ci sont orientées parallèlement au vecteur flux thermique $\vec{\psi}$ ; (ii) par contre, il révèle avec une sensibilité maximale les barrières thermiques orientées perpendiculairement à $\vec{\psi}$. En fait, cette sensibilité du détecteur $||\vec{\nabla}\Theta||_2$ à l'orientation de la barrière thermique est tout à fait logique ; on peut même pousser plus loin ce résultat en annonçant que la capacité du champ $||\vec{\nabla}\Theta||_2$ à détecter une
barrière thermique donnée est proportionnelle au cosinus de l'angle $\langle \tilde{n}, \tilde{\psi} \rangle$ dans lequel $\tilde{n}$ désigne le vecteur unitaire normal à la barrière thermique.

Dans le cadre des calculs par différences finies présentés à la figure 45, le vecteur flux thermique $\tilde{\psi}$ était dirigé principalement selon l'axe $\tilde{x}$, de sorte que l'amplitude du signal de détection d'une fissure de normale $\tilde{n}$ était proportionnel à $\cos(\tilde{n}, \tilde{x})$. Pour que le champ $\|\tilde{\nabla}\tilde{\theta}\|_2$ révèle les fissures verticales de l'échantillon, il faudrait que le vecteur $\tilde{\psi}$ soit dirigé principalement selon l'axe $\tilde{y}$ ; l'amplitude du signal de détection d'une fissure de normale $\tilde{n}$ serait alors proportionnel à $\cos(\tilde{n}, \tilde{y}) = \sin(\tilde{n}, \tilde{x})$. On disposerait donc ainsi de deux champs $\|\tilde{\nabla}\tilde{\theta}\|_2$ : le premier, baptisé $\|\tilde{\nabla}\tilde{\theta}_x\|_2$, correspondant à un transfert thermique dirigé principalement dans la direction $\tilde{x}$, et le second, baptisé $\|\tilde{\nabla}\tilde{\theta}_y\|_2$, correspondant à un transfert thermique dirigé principalement dans la direction $\tilde{y}$ ; ces deux champs $\|\tilde{\nabla}\tilde{\theta}_x\|_2$ et $\|\tilde{\nabla}\tilde{\theta}_y\|_2$ présentant des capacités de détection proportionnelles à $\cos(\tilde{n}, \tilde{x})$ et $\sin(\tilde{n}, \tilde{x})$ respectivement, il s'ensuit que le nouveau champ de gradient thermique généralisé $\|\tilde{\nabla}\tilde{\theta}_{xy}\|_2$ défini par :

$$ \|\tilde{\nabla}\tilde{\theta}_{xy}\|_2 = \sqrt{\|\tilde{\nabla}\tilde{\theta}_x\|_2^2 + \|\tilde{\nabla}\tilde{\theta}_y\|_2^2} \quad (82) $$

devrait présenter une sensibilité de détection qui soit totalement indépendante de l'orientation de la fissure mais uniquement fonction de l'importance de sa résistance thermique.

Nous avons essayé de mettre en application cette idée. Le changement d'orientation principale du vecteur flux thermique $\tilde{\psi}$ s'est fait non pas en écrivant un nouveau programme de différences finies adéquat mais beaucoup plus simplement en tournant de 90° l'échantillon initial de la figure 45(a).

La figure 46(a) rappelle la géométrie de l'échantillon d'étude. La figure 46(b) montre sous forme d'image le champ de température $\theta$ s'établissant dans le domaine de calcul : le flux thermique $\tilde{\psi}$ est bien porté principalement par le vecteur $\tilde{y}$, le champ $\theta$ est bien voisin du champ monodimensionnel homogénéisé dans la direction $\tilde{y}$ $\theta(\tilde{y}) = 1 - \frac{\tilde{y}}{\tilde{Y}}$, et le programme de différences finies fournit une conductivité thermique effective moyenne $\bar{\Lambda}$ de 0.9308 W/m/K, valeur assez voisine de celle (0.9404 W/m/K) tirée de l'analyse thermique dans la direction $\tilde{x}$. La figure 46(c) montre sous forme d'image le champ de déviation thermique $\delta\theta$ résultant du champ $\theta$. Enfin, la figure 46(d) rend compte du champ $\|\tilde{\nabla}\tilde{\theta}_y\|_2$ régnant dans le domaine de calcul : ce champ faisant apparaître les fissures de l'échantillon d'autant plus clairement que ces fissures sont plus orientées dans la direction normale à $\tilde{\psi}$, il ne "voit" quasiment pas les fissures horizontales et "voit" les fissures verticales avec une sensibilité maximale.

Disposant des deux champs de gradient thermique $\|\tilde{\nabla}\tilde{\theta}_x\|_2$ (figure 45(d)) et $\|\tilde{\nabla}\tilde{\theta}_y\|_2$ (figure 46(d)), nous avons calculé le champ de gradient thermique généralisé $\|\tilde{\nabla}\tilde{\theta}_{xy}\|_2$ au moyen de la formule (82) ; ce champ $\|\tilde{\nabla}\tilde{\theta}_{xy}\|_2$ est présenté sous forme d'image à la figure 47(b), à côté de la figure 47(a) qui rappelle la géométrie de l'échantillon fissuré. Cette image 47(b) parle d'elle-même : à la
Dans tout ce chapitre, l'échantillon sur lequel nous allons nous concentrer sera celui représenté sous forme d'image à la figure 51(a). Cette image a été produite au moyen d'un logiciel de morphologie mathématique baptisé Micromorph (logiciel développé au Centre de Morphologie Mathématique de l'Ecole des Mines de Paris) de la façon suivante : on a commencé par disséminer de façon aléatoire 75 pixels noirs dans un carré blanc de 256 x 256 pixels ; puis Micromorph a recherché le squelette des zones d'influence de ces 75 pixels, ce qui a fourni notre image d'étude.

Cette image peut être appréhendée de deux façons différentes : on peut y voir soit la structure d'un matériau polycristallin constitué de grains agglutinés les uns aux autres, soit l'architecture d'une matrice microporeuse faite de nanoparticules coalescées. Nous allons donc dans ce qui suit essayer des modélisations thermiques dans les deux cas de figure envisagés.

### 3-1-6-3-1- Matériau en grains

Dans ce chapitre, nous supposons avoir affaire à un matériau polycristallin constitué de grains agglomérés. Pour les besoins du calcul, il nous faut attribuer des valeurs de conductivité thermique locale $\lambda$ aux pixels blancs et noirs de la figure 51(a) ; si le choix de $\lambda$ pour les pixels blancs ne pose a priori pas de difficulté particulière (il suffit de prendre la valeur du matériau massif), il n'en va pas du tout de même pour les pixels noirs matérialisant les joints entre grains contigus, comme nous le verrons par la suite.

Pour notre calcul, nous avons supposé que les différents grains étaient séparés par de l'air (cette hypothèse sera discutée plus loin) ; nous avons donc assigné aux pixels noirs une conductivité thermique locale de 0.025 W/m/K, tandis que les pixels blancs ont été maintenus à une valeur de $\lambda$ de 1 W/m/K.

Les figures 51(b) et 51(c) présentent, sous forme d'images, les champs de température $\theta$ et de gradient thermique directionnel $||\vec{\nabla}\theta_x||_2$ calculés sur l'échantillon polycristallin. Rien de fondamentalement nouveau à signaler concernant le champ $||\vec{\nabla}\theta_x||_2$ ; par contre, comme le montre très clairement la figure 51(b) et encore plus clairement la figure 52, le champ $\theta$ qui s'établit dans le domaine de calcul a ceci de particulier que chaque grain de l'échantillon se trouve à température approximativement uniforme : les chutes de température entre la face chaude et la face froide de l'échantillon ont lieu quasi-exclusivement au niveau des résistances thermiques que constituent les couches d'air situées aux joints entre les différents grains.

Nous mettons ici le doigt sur un gros défaut de notre simulation : nous avons donné aux joints de grain les propriétés de transfert thermique de l'air. En d'autres termes, nous avons calculé les propriétés de conduction thermique d'un amas de grains simplement posés les uns sur les autres sans aucune adhérence entre eux. Une modélisation plus réaliste du transfert thermique au sein d'un milieu polycristallin devrait prendre en compte cette adhérence entre les grains, laquelle se traduit par des résistances thermiques d'interface beaucoup moins intenses que celles que nous avons implantées en séparant les grains par de l'air. Du point de vue de nos calculs, une telle modification serait très facile à apporter (elle consisterait simplement à augmenter la conductivité thermique locale aux joints de grain, c'est-à-dire aux pixels noirs) ; mais la mesure quantitative de la résistance thermique d'interface à un joint de grain dans un milieu polycristallin constitue encore aujourd'hui un sujet de recherche expérimentale dont l'état d'avancement est loin d'être abouti.

### 3-1-6-3-2- Matrice microporeuse
L'échantillon de la figure 51(a) peut être considéré avec un œil autre que celui des milieux polycristallins : il peut aussi faire penser aux architectures que l'on rencontre au sein des isolants microporeux. Cette classe de matériaux est l'objet de nombreuses études dans des domaines d'application très variés pour l'extraordinaire capacité d'isolation thermique qu'elle apporte : tandis que l'air, généralement considéré comme un très bon isolant thermique, a une conductivité thermique de 25 mW/m/K à pression et température ambiantes, la conductivité thermique effective est de quelques mW/m/K pour un matériau microporeux placé sous vide primaire.

Cette très grande capacité d'isolation thermique des matériaux microporeux trouve son explication dans la microstructure de ces matériaux. Ils sont en général constitués de particules de taille nanométrique de silice (le choix de la silice se justifie par la faiblesse de la conductivité thermique de ce matériau) et de particules de taille micrométrique de matériaux dits "opacifiants" dont le rôle consiste à absorber le rayonnement thermique au cours de sa progression au sein du milieu microporeux. Les fractions volumiques des nanoparticules et des microparticules étant respectivement d'environ 10% et de 1 ou 2%, il s'ensuit que les matériaux microporeux sont très poreux : la porosité avoisine les 90% et est ouverte, autorisant ainsi la raréfaction du gaz occupant le matériau par mise sous vide de celui-ci.

La microstructure d'un matériau microporeux se décline selon deux échelles d'observation différentes. A l'échelle du micromètre (c'est-à-dire à l'échelle des particules opacifiantes), le matériau semble constitué d'une matrice d'apparence homogène dans laquelle sont dispersées les microparticules en très faible fraction volumique. Mais l'observation à l'échelle du nanomètre (c'est-à-dire à l'échelle des particules de silice) de la matrice révèle que cette matrice prétendument homogène est en fait une architecture tridimensionnelle de nanoparticules coalescées constituant un squelette dans lequel le gaz se trouve emprisonné par viscosité. La figure 53 fait la synthèse de la description microstructurale à deux échelles qui vient d'être faite du matériau microporeux.

La modélisation du transfert thermique par conduction au sein du milieu d'étude doit intégrer cette description à deux échelles de sa microstructure ; elle doit donc commencer par l'étude du transfert thermique conductif et la détermination de la conductivité thermique effective $\Lambda_m$ de la matrice. Or évaluer cette conductivité thermique effective $\Lambda_m$ nécessite la connaissance préalable de 4 paramètres : (i) les conductivités thermiques $\lambda_s$ et $\lambda_g$ des deux constituants de la matrice, c'est-à-dire les nanoparticules de silice et le gaz ; (ii) la résistance thermique moyenne $\bar{R}$ produite au niveau de la zone de coalescence entre deux nanoparticules contiguës ; et (iii) la géométrie de la matrice, c'est-à-dire les positions dans l'espace tridimensionnel des nanoparticules la matérialisant.

Les 4 paramètres évoqués ci-dessus ($\lambda_s$, $\lambda_g$, $\bar{R}$ et géométrie) sont tous difficiles à connaître avec précision :

1. Les nanoparticules de silice ont des diamètres de quelques nm à quelques dizaines de nm, soit du même ordre de grandeur que le libre parcours moyen des phonons dans la silice à température ambiante. Il s'ensuit que la conductivité thermique $\lambda_s$ d'une nanoparticule de silice dépend de la géométrie de cette nanoparticule et qu'elle ne peut être prise égale à la conductivité thermique de la silice massive car elle lui est inférieure. Mais dans quelle mesure lui est-elle inférieure ? Est-elle inférieure d'un ou plusieurs ordres de grandeur, ou bien au contraire le gain en isolation thermique apporté par l'effet de taille est-il somme toute modeste ?

2. Les inclusions du squelette de nanoparticules de silice dans lesquelles le gaz se trouve emprisonné par viscosité ont des diamètres de quelques centaines de nm, soit du même ordre de grandeur que le libre parcours moyen des molécules gazeuses. Se pose donc un problème très connexe de celui évoqué ci-dessus, à savoir l'évaluation de la conductivité thermique $\lambda_g$ du gaz qui se trouve dans un état confiné.
3. Les chaînes de nanoparticules coalescées méritent une étude attentive : que se passe-t-il au niveau corpusculaire à la coalescence de deux nanoparticules ? Y a-t-il continuité ou non du réseau atomique ? Les zones de coalescence jouent-elles le rôle de barrières thermiques ? En d'autres termes, doit-on leur associer des résistances thermiques équivalentes $\overline{\mathfrak{R}}$ importantes, de sorte que le pouvoir isolant thermique de la matrice serait finalement plus lié aux barrières thermiques constituées par les zones de coalescence qu'au faible niveau de conductivité thermique atteint par les nanoparticules de silice ?

4. Enfin, de par les multiples intersections et circonvolutions qu'elle impose au flux conductif, la géométrie tridimensionnelle squelettique constituée par les nanoparticules de silice contribue elle aussi très certainement à l'extraordinaire pouvoir isolant thermique de la matrice. Vu l'échelle spatiale de cette architecture, une tomographie de la matrice pour en décrire de façon géométrique la microstructure squelettique paraît difficilement envisageable ; une modélisation du transfert conductif rendant correctement compte de la complexité du « cheminement thermique » dans la matrice doit donc intégrer une réflexion sur le problème de la génération par voie informatique (sur la base d'algorithmes générateurs aléatoires) de squelettes tridimensionnels de nanoparticules représentatifs de la microstructure réelle de la matrice.

Dans l'état actuel de nos travaux, nous ne sommes pas en mesure de répondre aux difficultés soulignées ci-dessus concernant la quantification des paramètres $\lambda_g$ et $\overline{\mathfrak{R}}$ (points 1. et 3.). Par contre, pour ce qui est de la géométrie de la matrice (point 4.), on peut tout à fait voir dans l'image de la figure 51(a) un squelette microporeux dans lequel les nanoparticules seraient localisées aux pixels noirs tandis que les pixels blancs définiraient les zones occupées par le gaz. Enfin, en ce qui concerne la conductivité thermique $\lambda_g$ d'un gaz à l'état confiné (point 2.), nous avons trouvé dans la littérature [7] la corrélation suivante :

$$ \lambda_g(T, P, \delta) = \frac{\lambda_g^0(T)}{1 + 2\beta \frac{\overline{\ell}}{\delta}} \qquad (83) $$

Dans cette formule, T, P et $\delta$ désignent respectivement la température du gaz, la pression du gaz et la dimension caractéristique du confinement ; $\lambda_g$ et $\lambda_g^0$ sont respectivement la conductivité thermique du gaz à l'état confiné (fonction de T, P et $\delta$) et la conductivité thermique du même gaz à l'état libre (uniquement fonction de T) ; $\beta$ est un paramètre sans dimension caractéristique du gaz ; enfin, $\overline{\ell}$ désigne le libre parcours moyen des molécules gazeuses que l'on calcule au moyen de la relation classique

$$ \overline{\ell} = \frac{1}{\sqrt{2 n \pi d^2}} \qquad (84) $$

dans laquelle n et d sont respectivement le nombre de molécules par unité de volume (qui vaut $\frac{P N_A}{R T}$ dans l'hypothèse des gaz parfaits, avec $N_A$ = nombre d'Avogadro et R = constante universelle des gaz parfaits) et le diamètre des molécules.
Nous avons appliqué les formules (83) et (84) ci-dessus pour évaluer la conductivité thermique de l'air emprisonné dans la matrice microporeuse. Les données suivantes ont été retenues pour les calculs numériques : T = 300 K, donc $\lambda_g^0 = 25 \times 10^{-3}$ W/m/K ; P = 1 bar (matrice microporeuse sous air à pression atmosphérique) puis P = 1 mbar (matrice microporeuse sous vide primaire d'air) ; $\delta = 100$ nm (valeur estimée pour le diamètre caractéristique des inclusions retenant l'air dans la matrice microporeuse) ; $\beta = 2.429$ (valeur tirée de [7]) ; enfin, d = 3.15 Å (diamètre moléculaire de l'azote gazeux N₂, constituant majoritaire de l'air, tiré de la référence [8]). Les deux valeurs obtenues pour la conductivité thermique gazeuse "confinée" $\lambda_g$ sont d'une part $\lambda_g(P = 1$ bar) $= 4.5 \times 10^{-3}$ W/m/K et d'autre part $\lambda_g(P = 1$ mbar) $= 5.5 \times 10^{-6}$ W/m/K ; on fera les deux remarques suivantes concernant ces deux valeurs numériques :

1. A pression atmosphérique, l'air occupant les inclusions dans la matrice microporeuse se trouve déjà dans un état confiné : le libre parcours moyen $\bar{l}$ de ses molécules gazeuses vaut 94 nm, il est donc très voisin du diamètre caractéristique d des inclusions, et par conséquent la conductivité thermique $\lambda_g$ est sensiblement plus faible que son homologue "non confinée" $\lambda_g^0$ (elle lui est environ 5 fois inférieure).

2. Ceci est encore plus vrai lorsque la matrice microporeuse est placée sous vide primaire : en effet, le libre parcours moyen $\bar{l}$ des molécules étant inversement proportionnel à la pression P, diviser P par 1000 multiplie $\bar{l}$ par 1000 et par conséquent renforce l'état de confinement du gaz. Dans le cas traité ci-dessus de l'air à T = 300 K et P = 1 mbar emprisonné dans des inclusions de diamètre caractéristique d = 100 nm, cet air devient un isolant thermique parfait du point de vue du phénomène de conduction thermique puisque sa conductivité thermique "confinée" est de quelques $10^{-6}$ W/m/K.

Nous avons essayé de simuler le transfert conductif au travers d'une matrice microporeuse au moyen de notre programme de différences finies en procédant de la façon suivante : (i) ne maîtrisant pas les effets de taille sur la valeur de la conductivité thermique $\lambda_s$ d'une nanoparticule de silice, nous avons fixé ce paramètre $\lambda_s$ à 1 W/m/K ; (ii) nous avons donné à la conductivité thermique $\lambda_g$ du gaz les deux valeurs calculées plus haut, c'est-à-dire $4.5 \times 10^{-3}$ W/m/K (matrice sous air à pression atmosphérique) et $5.5 \times 10^{-6}$ W/m/K (matrice sous air à 1 mbar) ; (iii) faute de connaissance sur le sujet, nous n'avons pas pris en compte de résistances thermiques de contact moyennes $\bar{R}$ au niveau des zones de coalescence entre nanoparticules de silice voisines ; (iv) enfin, la géométrie que nous avons retenue pour décrire la matrice microporeuse est celle de la figure 51(a), les pixels noirs matérialisant les nanoparticules et les pixels blancs le gaz. Les résultats des calculs numériques nous ont paru quelque peu surprenants au premier abord car ils ont conclu à une conductivité thermique effective moyenne $\bar{\Lambda}$ de la matrice valant $5.5352 \times 10^{-3}$ W/m/K sous air à pression atmosphérique et aussi faible que $6.8097 \times 10^{-6}$ W/m/K sous air à 1 mbar ! Ces deux valeurs numériques sont à rapprocher des valeurs de $\lambda_g$ utilisées dans les calculs : il semble que notre programme de différences finies ne soit quasiment pas sensible à l'existence du "filet conducteur thermique" constitué par le squelette de nanoparticules puisque, sous pression atmosphérique d'air comme sous vide primaire d'air, il prévoit pour la matrice une conductivité thermique effective moyenne $\bar{\Lambda}$ tout juste légèrement supérieure à la conductivité thermique du gaz.

Cette insensibilité de notre programme au squelette de nanoparticules s'explique par la méthode que nous avons employée lors de la discrétisation du problème (11) en différences finies (voir chapitre 2-1-) : de façon générale, le bilan énergétique au nœud (i, j) ne fait intervenir que ce nœud et ses 4 voisins immédiats que sont le nœud "au-dessus" (i-1, j), le nœud "au-dessous" (i+1, j), le nœud "à gauche" (i, j-1) et le nœud "à droite" (i, j+1). Il s'ensuit que si autour d'un nœud (i, j) "noir" (c'est-à-dire en lequel se trouve une nanoparticule de silice) la configuration est par exemple la suivante :
<figure><img src="image_4.png" /></figure>

alors l'information de bonne conductance thermique entre les nœuds (i, j) et (i, j-1) d'une part, (i, j) et (i+1, j) d'autre part, sera correctement répercutée par l'algorithme ; en revanche, si autour du nœud (i, j) "noir" la configuration est par exemple la suivante :

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

alors l'information de bonne conductance thermique entre les nœuds (i, j) et (i-1, j-1) d'une part, (i, j) et (i+1, j-1) d'autre part, ne sera pas répercutée de façon suffisamment forte par l'algorithme. En d'autres termes, si deux nanoparticules voisines ne sont pas alignées horizontalement ou verticalement dans le maillage, alors le programme insérera une résistance thermique artificielle sur le segment reliant les deux nanoparticules et par là-même ignorera l'information de bonne conduction thermique entre elles.
Ce défaut de notre programme de différences finies est illustré à la figure 54. La figure 54(a) montre, sous forme d'image, la partie supérieure gauche de l'échantillon d'étude agrandie dans un rapport 4 (l'échantillon au complet est présenté sous forme d'image à la figure 51(a)) ; quant à la figure 54(b), elle montre, également sous forme d'image et avec la même échelle, le champ de gradient thermique généralisé $\|\vec{\nabla}\Theta_{xy}\|_2$ régnant dans cette partie supérieure gauche de l'échantillon.

Comme on peut le constater sur cette figure 54(b), le champ $\|\vec{\nabla}\Theta_{xy}\|_2$ détecte des résistances thermiques intenses (pixels jaunes à rouges) là où il y a discontinuité importante des chaînes de pixels noirs sur l'image discrétisée de l'échantillon (figure 54(a)) : en insérant de la sorte des résistances thermiques artificielles dans les chaînes constituées par les nanoparticules coalescées, notre programme annihile la capacité de conduction thermique de ces chaînes et donc sous-estime la conductivité thermique effective $\Lambda_m$ de la matrice. Par contre, là où les nanoparticules sont correctement alignées sur l'image de la figure 54(a), le champ $\|\vec{\nabla}\Theta_{xy}\|_2$ est minimal (zones en bleu foncé sur l'image de la figure 54(b)) : en ces pixels, l'information de bonne conductance thermique est correctement répercutée par le programme, et donc le gradient thermique est quasi-nul.

Que faire pour remédier à cette imperfection de notre programme de différences finies ? L'idéal serait de remplacer le présent algorithme par une autre routine plus sophistiquée intégrant les effets des 8 nœuds voisins (en géométrie bidimensionnelle) dans le calcul de la température en un nœud donné, ceci évidemment au prix de la rapidité du calcul. En attendant cette amélioration, on peut avancer quelque peu dans le domaine de la modélisation du transfert thermique au sein de la matrice microporeuse par une astuce toute simple qui consiste à "éroder" l'image de la matrice. Cette opération morphologique qu'est l'érosion a ici pour résultat d'épaissir les chaînes de pixels noirs, comme on peut le constater aux figures 55 et 56 : chaque zone blanche de la figure 55(a) voyant sa frontière érodée d'un pixel, il s'ensuit que l'épaisseur moyenne des chaînes de pixels noirs passe de 1 à 3 pixels. On comprend aisément, au vu de la figure 56(b), que les problèmes évoqués plus haut concernant notre programme de différences finies seront contournés avec l'image érodée : en effet, chaque pixel noir de cette image aura parmi ses 4 pixels influents voisins au moins un pixel noir, de sorte que la "continuité thermique" des chaînes de nanoparticules sera reconnue par le programme. Nous avons lancé un calcul de conductivité thermique effective sur la matrice érodée (image de la figure 55(b)) en donnant aux pixels noirs une conductivité thermique $\lambda_s$ de 1 W/m/K et aux pixels blancs la conductivité thermique de l'air confiné à pression atmosphérique $\lambda_g = 4.5 \times 10^{-3}$ W/m/K ; le programme a fourni une conductivité thermique effective moyenne $\bar{\Lambda}$ valant $9.3365 \times 10^{-2}$ W/m/K, soit - comme attendu - beaucoup plus que la valeur $\bar{\Lambda}$ issue du calcul effectué dans les mêmes conditions sur la matrice non érodée (rappelons que cette valeur $\bar{\Lambda}$ était de $5.5352 \times 10^{-3}$ W/m/K).

### 3-1-6-4- Dispersion dans une matrice

La modélisation du transfert thermique par conduction au sein d'un matériau microporeux, évoquée au paragraphe précédent pour ce qui concernait la quantification de la conductivité thermique effective $\Lambda_m$ de la matrice microporeuse, pose un autre problème auquel notre programme par différences finies est susceptible d'apporter des solutions : comment passer de la conductivité thermique effective $\Lambda_m$ de la matrice à celle $\Lambda$ du milieu microporeux dans sa globalité ? Pour répondre à cette question, il faut se remémorer la description microstructurale du matériau qui a été faite plus haut : à l'échelle micrométrique, celui-ci apparaît comme étant constitué en grande majorité de la matrice microporeuse dans laquelle sont dispersées de façon approximativement uniforme les microparticules opacifiantes (voir figure 53(a)). Cette dispersion de microparticules étant de très faible fraction volumique (1 ou 2%), on peut raisonnablement émettre une hypothèse
de type "non-interférence des microparticules les unes avec les autres en ce qui concerne le champ thermique": une telle hypothèse consiste à considérer qu'une microparticule est suffisamment isolée dans la matrice pour que la perturbation apportée au champ de température par sa présence n'atteigne pas les microparticules voisines. Dès lors, de nombreuses corrélations existent dans la littérature pour évaluer de façon simple la conductivité thermique effective $\Lambda$ du mélange {matrice microporeuse + microparticules opacifiantes dispersées de façon uniforme et en faible fraction volumique} à partir des conductivités thermiques de la matrice $\Lambda_m$ et des microparticules opacifiantes $\lambda_o$ et de la fraction volumique des microparticules opacifiantes $\alpha_o$; on pourra par exemple se référer à la publication [9] dans laquelle les auteurs passent en revue un certain nombre de formulations de la conductivité thermique effective $\Lambda$ proposées par différents chercheurs.

Si les corrélations évoquées ci-dessus permettent d'évaluer la conductivité thermique effective $\Lambda$ dans de nombreux cas de mélanges de matériaux, elles ont aussi leurs limites d'applicabilité. En ce qui concerne nos matériaux microporeux, une difficulté apparaît liée au contraste thermique important entre les deux constituants du "mélange" que sont d'une part la matrice et d'autre part les microparticules opacifiantes : en effet, tandis que les microparticules ont une conductivité thermique $\lambda_o$ de l'ordre du W/m/K (à la différence des nanoparticules de silice, les microparticules opacifiantes ont un diamètre nettement supérieur au libre parcours moyen des phonons, de sorte qu'on peut retenir pour les microparticules la conductivité thermique du matériau massif), la conductivité thermique de la matrice $\Lambda_m$ peut le cas échéant tomber à quelques $10^{-3}$ W/m/K. S'ensuit alors la réflexion suivante : les microparticules sont certes en minorité dans le mélange puisque leur fraction volumique n'est que de 1 ou 2%; en revanche, leur conductivité thermique peut être jusqu'à 100 fois supérieure à celle de la matrice ; alors lequel de ces deux arguments l'emporte : l'argument "faible fraction volumique" conduisant à une conductivité thermique effective du mélange peu différente de celle du constituant majoritaire, ou l'argument "forte conductivité thermique" conduisant au contraire à une conductivité thermique effective du mélange augmentant sensiblement par rapport à celle du constituant majoritaire ?

Notre programme de différences finies peut aider à répondre à cette question. A titre d'illustration, nous avons dispersé de façon aléatoire et selon une loi uniforme 3200 pixels noirs dans un carré de 400 x 400 pixels blancs ; nous avons ainsi obtenu un mélange ressemblant fort à la dispersion de microparticules opacifiantes dans la matrice microporeuse puisque présentant une fraction surfacique de pixels noirs égale à 2%. Nous avons ensuite fait calculer notre programme de différences finies sur cette image en attribuant à ses pixels blancs et noirs des conductivités thermiques locales valant respectivement 0.01 W/m/K et 1 W/m/K ; la géométrie de l'échantillon d'étude est présentée à la figure 57(a), et les résultats des calculs sont exposés aux figures 57(b) (champ de température $\theta$), 57(c) (champ de déviation thermique $\delta\theta$) et 57(d) (champ de gradient thermique directionnel $\|\vec{\nabla}\theta_x\|_2$).

Peu de choses à signaler concernant ces trois dernières figures. Le champ $\theta$ est très voisin du champ monodimensionnel homogénéisé $\theta(x) = 1 - \frac{x}{X}$ (les valeurs couvertes par le champ $\delta\theta$ sont, en valeur absolue, inférieures à 3‰) ; le présent échantillon peuplé de sa distribution de particules conductrices semble donc homogénéisable du point de vue de la conduction thermique, et dès lors on pourra par la suite s'intéresser à la valeur de la conductivité thermique effective moyenne $\bar{\Lambda}$ fournie par notre programme de calcul. Le champ $\|\vec{\nabla}\theta_x\|_2$ détecte des gradients faibles (pixels de couleur bleue foncée) là où plusieurs particules conductrices sont en contact, et des gradients forts (pixels de couleur rouge) aux voisinages de ces petits amas de particules ; mais de façon générale ce champ $\|\vec{\nabla}\theta_x\|_2$ est à peu près uniforme sur tout le domaine de calcul.
La valeur de $\bar{\Lambda}$ fournie par le programme est de $1.0297 \times 10^{-2}$ W/m/K, avec une erreur relative maximale $\varepsilon_{\Lambda}$ valant $6.2 \times 10^{-6}$. Il semble donc, dans le cas des présents calculs, que la dissémination de particules conductrices dans la matrice isolante produit un matériau hétérogène dont la conductivité thermique effective n'est pas sensiblement modifiée par rapport à celle de la matrice ; l'effet "faible fraction surfacique" de l'ensemencement semblerait donc ici l'emporter sur l'effet "forte conductivité thermique" de celui-ci ... Une telle conclusion est à reporter à plus tard, à l'issue d'études plus approfondies :

1. L'échantillon d'étude de la figure 57(a) doit être considéré comme une représentation particulière de la variable aléatoire "carré blanc + dissémination de pixels noirs avec une fraction surfacique de 2%", car ces pixels noirs sont distribués de façon aléatoire et selon une loi uniforme dans le carré. Dès lors, la valeur $\bar{\Lambda}$ issue du programme de calcul doit elle aussi être considérée comme une variable aléatoire ; il serait donc à l'avenir important de se pencher sur les propriétés statistiques de cette variable aléatoire $\bar{\Lambda}$ (en commençant par son espérance mathématique et son écart-type), par exemple par la génération aléatoire d'une population de "carrés disséminés" puis le calcul des valeurs de $\bar{\Lambda}$ obtenues sur cette population.

2. Pour les travaux présentés dans cette partie comme pour ceux qui ont été faits précédemment dans le cadre de la modélisation du transfert thermique au sein d'une matrice microporeuse (voir chapitre 3-1-6-3-2-), il devient maintenant nécessaire de produire des domaines de calcul et de résoudre des champs thermiques dans l'espace tridimensionnel. Ainsi, on peut imaginer qu'un cube homogène saupoudré d'une dispersion de particules avec une fraction volumique de 2% ne soit pas équivalent, du point de vue des phénomènes de transfert, à son pendant bidimensionnel, c'est-à-dire un carré homogène saupoudré d'une dispersion de particules avec une fraction surfacique de 2%. Le chapitre ci-dessous est consacré au seul travail que nous ayons eu l'occasion d'effectuer jusqu'à présent dans l'espace à 3 dimensions.

## 3-2- Un exemple en géométrie tridimensionnelle : estimation numérique des conductivités thermiques effectives de deux échantillons élastomères chargés

Dans le cadre d'un travail de diffusion technologique, il nous a été demandé d'étudier les propriétés de transfert thermique stationnaire de deux matériaux constitués pour l'essentiel d'une matrice élastomère et chargés de grains d'aluminium à deux fractions massiques différentes : 20% pour l'un et 40% pour l'autre.

Tout d'abord, il a fallu "construire informatiquement" deux échantillons tridimensionnels représentatifs des deux matériaux à caractériser. Les études de microstructure ayant montré que les grains d'aluminium présentaient des tailles assez peu variables (de l'ordre de la dizaine de µm) et étaient distribués de façon relativement uniforme au sein des deux matériaux, l'étape de génération informatique des deux échantillons a consisté simplement à disperser de manière aléatoire des billes d'aluminium dans des cubes remplis de matrice élastomère tout en respectant les contraintes de fraction volumique d'aluminium dans les deux matériaux à caractériser. Préalablement à ce travail informatique, il a donc été procédé à des mesures de masse volumique des deux matériaux qui ont permis de quantifier les deux fractions volumiques d'aluminium : 9.8% pour le matériau le moins chargé, et 22.1% pour le matériau le plus chargé.

Les figures 58 et 59 montrent les deux "échantillons informatiques" qui ont été engendrés pour les calculs : ils sont tous deux constitués de 50 × 50 × 50 cubes élémentaires dont des fractions de respectivement 9.8% (soit donc 12250 cubes) et 22.1% (soit donc 27625 cubes) sont occupées par
des billes d'aluminium, les cubes restants étant quant à eux remplis d'élastomère. Notons au passage que 50 × 50 × 50 cubes élémentaires signifient autant de températures $\theta_{ij,k}$ inconnues, c'est-à-dire autant de mémoire informatique qu'un champ bidimensionnel $\theta_{ij}$ sur un carré découpé en $50^{3/2} \approx 350$ éléments ; par contre, le temps de calcul d'un problème tridimensionnel à $50^3$ inconnues est notablement supérieur à celui d'un problème bidimensionnel à $350^2$ inconnues, car le nombre de voisins influents d'un nœud donné passe de 4 à 6 lorsque la dimension de l'espace de travail passe de 2 à 3.

Pour le calcul des conductivités thermiques effectives des deux échantillons évoqués ci-dessus, nous avons dû donner des valeurs numériques aux conductivités thermiques des deux constituants que sont l'aluminium et la matrice élastomère. La conductivité thermique de l'aluminium a été trouvée dans la littérature : 240 W/m/K. Quant à la conductivité thermique de la matrice, celle-ci a été déduite de mesures expérimentales de masse volumique, de chaleur massique et de diffusivité thermique effectuées sur un échantillon non chargé en aluminium, donc constitué à 100% d'élastomère ; lesdites mesures ont abouti à la valeur de 0.20 W/m/K pour la conductivité thermique de la matrice. Notons que cette valeur nous a paru tout à fait recevable : en effet, l'élastomère entrant dans la composition des deux échantillons "ressemblait" beaucoup à un matériau du type PVC, et la conductivité thermique du PVC à température ambiante vaut 0.16 W/m/K.

Les calculs ont abouti aux deux valeurs suivantes de la conductivité thermique effective : $\bar{\Lambda} = 0.24$ W/m/K pour le matériau le moins chargé (avec une erreur relative maximale $\varepsilon_\Lambda$ valant $5.7 \times 10^{-6}$), et $\bar{\Lambda} = 0.41$ W/m/K pour le matériau le plus chargé (avec une erreur relative maximale $\varepsilon_\Lambda$ valant $1.2 \times 10^{-5}$). On fera au sujet de cette étude les deux remarques suivantes :

1. Il est intéressant de constater que si la présence des grains d'aluminium dans la matrice élastomère se traduit de façon assez modeste dans la conductivité thermique du mélange le moins chargé, en revanche cette présence double la conductivité thermique dans le mélange le plus chargé : celle-ci passe en effet de 0.20 W/m/K pour la matrice pure à 0.41 W/m/K pour la matrice chargée à 40% en masse, soit à 22.1% en volume.

2. Le matériau hétérogène qui a été étudié dans le cadre de ce travail de diffusion technologique est tout à fait intéressant en ce sens que, de par sa microstructure et la nature de ses constituants, il présente plusieurs points communs avec les matériaux microporeux tels que nous les avons abordés dans la partie 3-1-6-4- précédente : qu'il s'agisse d'un matériau de la présente étude observé à une échelle de l'ordre de la dizaine de µm (c'est-à-dire de l'ordre du diamètre des grains d'aluminium) ou d'un matériau microporeux observé à une échelle de l'ordre du µm (c'est-à-dire de l'ordre du diamètre des particules opacifiantes), les deux matériaux apparaissent comme étant constitués d'une matrice dans laquelle sont dispersées des particules de façon approximativement uniforme ; de plus, les deux matériaux partagent cette propriété de renfermer chacun deux constituants de conductivités thermiques radicalement différentes, le constituant "dispersé" étant dans les deux cas énormément plus conducteur thermique que le constituant "matrice". Pour toutes ces raisons, nous sommes d'avis que les mélanges {matrice élastomère + grains d'aluminium} pourraient à l'avenir nous être bien utiles pour avancer dans la compréhension du transfert thermique conductif au sein des milieux microporeux : en effet, les mélanges {matrice élastomère + grains d'aluminium} se prêtant incomparablement mieux à la caractérisation thermique que les mélanges {matrice microporeuse + particules opacifiantes} d'une part, les conductivités thermiques de leurs constituants étant considérablement mieux maîtrisées d'autre part, ces matériaux seront sans aucun doute beaucoup plus faciles à appréhender dans le cadre d'études de base en homogénéisation thermique alliant mesures expérimentales et modélisation numérique.
## 4- Conclusion - les suites à donner à ce travail

Dans ce document, nous avons voulu aborder la modélisation du transfert conductif stationnaire dans les milieux hétérogènes ; pour ce faire, nous avons développé un programme informatique de simulation thermique que nous avons ensuite testé sur différentes architectures matérielles.

Après avoir, en guise d'introduction au chapitre 1-, posé le problème mathématique de "conductivimétrie thermique numérique" à traiter, nous avons décrit dans le chapitre 2- la méthode de résolution que nous avons retenue : il s'agit d'une simple technique de discrétisation en différences finies, laquelle aboutit à un système linéaire qui est résolu au moyen de l'algorithme de Fletcher et Reeves. La technique de discrétisation par différences finies que nous avons implantée mériterait d'être sophistiquée de telle sorte que le bilan énergétique intègre les contributions thermiques de plus de nœuds voisins qu'il ne le fait actuellement. Quant à l'algorithme de Fletcher et Reeves, deux exigences principales nous ont mené vers le choix de cet algorithme itératif : d'une part la rapidité de calcul, et d'autre part la faible consommation de mémoire vive d'ordinateur. Jusqu'à présent, l'algorithme a convergé de façon très satisfaisante même dans des situations numériques critiques de matériaux hétérogènes présentant de forts contrastes thermiques ; cela étant, il serait intéressant de comparer ses performances à celles du solveur BASIS, outil numérique qui a été employé récemment par un groupe de chercheurs oeuvrant dans un domaine très connexe de nos intérêts (voir références [5] et [6]).

La partie 3- de ce rapport a été consacrée dans un premier temps à la validation de notre programme de calcul : pour cela, on a fait appel à des géométries simples sur lesquelles les champs de température pouvaient être déterminés de manière analytique. Dans un deuxième temps, on a soumis notre programme à des géométries matérielles hétérogènes assez variées, bidimensionnelles pour la grande majorité, sensées constituer des images simplifiées mais réalistes approchant les microstructures de grandes familles de matériaux comme les composites {fibres + résine}, les isolants thermiques microporeux ou les mélanges binaires {matrice + charges}. De tous les calculs et études qui ont été effectués, on retiendra les principaux résultats qui suivent :

1. L'étude des matériaux composites a été abordée sous deux éclairages différents : d'abord celui de la périodicité spatiale, ensuite celui de l'anisotropie.

L'aspect "périodicité spatiale" a été traité d'une part au moyen de notre programme et d'autre part par le biais d'un outil mathématique qui vient naturellement à l'esprit lorsqu'on a affaire à des signaux périodiques, à savoir l'analyse spectrale. La comparaison des performances des deux techniques d'analyse a conclu à leur complémentarité : si l'analyse spectrale produit des champs de déviation thermique $δθ$ en très bon accord avec ceux prévus par notre programme moyennant une quantité de calcul beaucoup moindre, en revanche elle pèche par sa relative incapacité à quantifier correctement la conductivité thermique effective d'un milieu hétérogène périodique, en particulier lorsque celui-ci présente un champ de conductivité thermique à fortes discontinuités. Cependant, pour compléter le volet "analyse spectrale" et conclure définitivement quant à sa faiblesse principale, il serait bon d'améliorer quelque peu la technique de résolution du système linéaire (70) permettant la détermination de la série de Fourier multidimensionnelle du champ $δθ$ : ainsi, le calcul de toutes les harmoniques accessibles de ce champ $δθ$ permettrait peut-être de produire des valeurs de conductivité thermique effective plus satisfaisantes. Dans un même souci de fignolage des techniques, signalons que notre programme de différences finies, dans son état actuel, n'exploite pas du tout les propriétés de périodicité spatiale du champ de température dans les directions autres que celle du transfert thermique principal ; tirer parti de ces propriétés par l'implantation de conditions
aux limites de périodicité adaptées permettrait de diminuer de façon conséquente les consommations de mémoire vive et de temps de calcul de notre programme.

L'aspect "anisotropie" a été effleuré au moyen d'un "échantillon informatique" hétérogène orienté. Notre programme a été mis à contribution pour l'évaluation de conductivités thermiques effectives directionnelles correspondant à différentes inclinaisons de l'échantillon - et donc de son anisotropie - par rapport à la direction du transfert thermique principal. Parallèlement, un modèle très simple (voire simpliste ?) d'homogénéisation thermique nous a conduits dans un premier temps à définir les propriétés d'un échantillon homogène anisotrope réputé équivalent à notre échantillon hétérogène anisotrope d'étude, et dans un deuxième temps à tirer de la théorie classique de la conduction thermique anisotrope une expression analytique de la conductivité thermique directionnelle de l'échantillon homogène équivalent. Les indicatrices de conductivité thermique directionnelle des deux échantillons ne coïncidant pas, nous avons provisoirement conclu à une inadéquation du modèle d'homogénéisation thermique que nous avons développé : un échantillon hétérogène anisotrope périodique de maille non infiniment petite vis-à-vis de ses dimensions caractéristiques ne pourrait pas, en ce qui a trait au transfert thermique homogénéisé, être remplacé par un simple tenseur de conductivité thermique anisotrope équivalent. Cela étant, notre modèle numérique par différences finies fournit des indicatrices de conductivité thermique effective directionnelle qui sont peut-être exactes ; pour vérifier ce point, il faudrait prolonger l'étude dans deux directions : d'une part la voie expérimentale, par la collection d'indicatrices de conductivité au moyen des techniques classiques de caractérisation thermique, et d'autre part la voie de la modélisation, par la confrontation de nos résultats à ceux issus d'autres approches d'homogénéisation thermique. Dans un futur plus proche, il serait bon de poursuivre la réflexion en ce qui a trait à l'optimisation numérique : pour prolonger ce qui a été dit plus haut, un échantillon périodique incliné reste un échantillon périodique et par conséquent induit peut-être des symétries particulières dans les champs de température ou de déviation thermique évalués par le programme. La prise en compte de ces symétries lors de l'étape de calcul numérique permettrait de libérer beaucoup d'énergie informatique (sous forme de mémoire vive et de temps de calcul) que l'on pourrait alors affecter à une définition plus fine des géométries.

2. Les matériaux isolants thermiques microporeux ont été étudiés à leurs deux échelles caractéristiques : d'une part l'échelle du micromètre à laquelle ils se présentent comme étant constitués d'une matrice d'apparence homogène dans laquelle des microparticules sont dispersées en très faible fraction volumique, et d'autre part l'échelle du nanomètre à laquelle la matrice prétendument homogène apparaît comme étant en réalité une architecture squelettique tridimensionnelle de nanoparticules coalescées.

L'évaluation de la conductivité thermique effective de la matrice microporeuse est un problème assez compliqué pour plusieurs raisons : (i) la taille des nanoparticules ; (ii) l'état confiné du gaz occupant les pores ; (iii) les résistances thermiques éventuelles au niveau des zones de coalescence ; (iv) la géométrie squelettique tridimensionnelle constituée par les chaînes de nanoparticules (pour plus de détails concernant ces 4 sources de difficulté, se référer au chapitre 3-1-6-3-2-). Nous avons quelque peu avancé en ce qui concerne les points (ii) et (iv), d'une part grâce à une corrélation (trouvée dans la littérature) permettant de quantifier la conductivité thermique d'un gaz dont le libre parcours moyen des molécules est du même ordre de grandeur que la taille caractéristique des inclusions dans lesquelles il se trouve piégé, et d'autre part à l'aide d'un opérateur de morphologie mathématique (en l'occurrence le squelette des zones d'influence) qui nous a permis de générer une image bidimensionnelle pouvant évoquer entre autres une architecture squelettique plane de nanoparticules. Nous nous sommes alors heurtés, lors du calcul sur les structures fines de cette image, à une déficience fâcheuse de notre discrétisation qui a
déjà été évoquée plus haut : le schéma de différences finies que nous avons retenu ne gère pas les nœuds diagonaux. Il faudra donc à l'avenir remédier à ce défaut, et aussi réfléchir à l'extension de l'opérateur "squelette des zones d'influence" à la génération de structures dans l'espace tridimensionnel.

3. Les matériaux microporeux se présentant, à l'échelle du micromètre, comme étant constitués d'une matrice d'apparence homogène chargée d'une dispersion de microparticules en très faible fraction volumique, nous nous sommes intéressés au problème relativement classique de la détermination de la conductivité thermique effective des mélanges binaires {matrice + charges}. Certes, de nombreuses corrélations existent dans la littérature pour répondre à cette question ; cependant, le contraste thermique très fort entre les deux constituants d'un "mélange microporeux" (les microparticules peuvent être plusieurs centaines de fois plus conductrices que la matrice) incite à une grande circonspection lorsqu'il s'agit de l'application de ces corrélations à un tel mélange. Notre programme de différences finies a donc calculé les conductivités thermiques effectives d'échantillons bidimensionnels et tridimensionnels constitués de matrices chargées à faibles taux volumiques de grains beaucoup plus thermiquement conducteurs qu'elles. Ce travail est là encore à l'état d'ébauche : maintenant que l'outil numérique existe et fonctionne, il va falloir l'utiliser de façon intensive sur des populations importantes d'"échantillons informatiques aléatoires" de façon à pouvoir étudier les propriétés statistiques de la variable aléatoire "conductivité thermique effective" ; il sera alors intéressant de comparer l'espérance mathématique de cette variable aléatoire aux conductivités thermiques effectives fournies par les corrélations classiques. Par ailleurs, signalons que dans le cadre d'un sujet de diffusion technologique, nous avons été amenés à évaluer les conductivités thermiques effectives de deux mélanges binaires constitués d'une matrice élastomère chargée de grains d'aluminium ; pour plusieurs raisons liées (i) à la bonne connaissance des conductivités thermiques des deux constituants, (ii) au très fort contraste thermique qui existe entre eux et (iii) à la grande facilité de caractérisation thermique de ces mélanges comparativement aux mélanges microporeux, nous sommes d'avis que les mélanges {matrice élastomère + grains d'aluminium} pourraient à l'avenir nous servir de matériaux images beaucoup plus faciles à appréhender dans des études de base que les mélanges microporeux {matrice microporeuse + microparticules}.

4. Enfin, divers échantillons (présentant une fissure isolée ou une distribution de fissures, poreux ou polycristallin) ont fait l'objet d'études sommaires et ont amené leurs difficultés propres. Parmi ceux-ci, l'échantillon polycristallin constitué de grains agglomérés mérite une mention particulière : l'exécution de notre programme de différences finies ayant abouti à des grains approximativement isothermes tandis que les chutes de température se produisaient quasi-exclusivement au niveau des joints, nous avons été amenés à constater que le calcul que nous venions d'effectuer péchait par notre méconnaissance des valeurs des résistances thermiques d'interface aux joints de grain.

Pour conclure, on soulignera que beaucoup de chantiers ont été ouverts avec ce travail préliminaire ; ceux-ci intéressent une grande variété de matériaux parmi lesquels deux catégories particulières de grand intérêt pour nos applications : les matériaux composites et les isolants thermiques microporeux. Parallèlement, une multitude de difficultés nouvelles se sont fait jour dans des domaines très divers : analyse numérique, morphologie mathématique, physique des petites échelles pour ce qui concerne les isolants microporeux ou encore maîtrise des propriétés thermiques des constituants pour ce qui concerne les matériaux composites. Il faudra à l'avenir continuer à avancer dans toutes ces directions pour améliorer notre connaissance des phénomènes de transfert dans les matériaux complexes.
Liste des
figures
<figure><img src="image_6.png" /><figcaption>Figure 1 : la géométrie considérée pour nos calculs (cas 2-D) et les conditions aux limites qui lui sont imposées.</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>Figure 2 : discrétisation de l'équation de volume.</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>Figure 3 : notion de conductivité thermique équivalente $\lambda_{NG}$ décrivant le conducteur thermique entre les noeuds N et G.</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>Figure 4 : discrétisation de l'équation de volume avec prise en compte des deux conditions aux limites de température imposée.</figcaption></figure>

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

<figure><img src="image_11.png" /></figure>
<figure><img src="image_12.png" /><figcaption>Figure 5 : discrétisation de l'équation de volume avec prise en compte des deux conditions aux limites d'adiabaticité.</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>Figure 6 : discrétisation de l'équation de volume aux quatre angles du domaine d'étude.</figcaption></figure>
<figure><img src="image_14.png" /><figcaption>Figure 7 : allures de la matrice A et du vecteur $\bar{b}$ (les ronds blancs et noirs désignent respectivement les éléments nuls et les éléments non nuls des deux tableaux).</figcaption></figure>

<figure><img src="image_15.png" /></figure>
<figure><img src="image_16.png" /><figcaption>Figure 8 : échantillon choisi pour l'étude de la convergence de l'algorithm de Fletcher et Reeves.</figcaption></figure>
Figure 9 : évolution en fonction de l'indice n d'itération des normes $\|\bar{x}_{n+1} - \bar{x}_n\|_\infty$ caractérisant la convergence des algorithmes de Fletcher et Reeves (courbe en trait continu) et de relaxation (courbe en trait pointillé).

<figure><img src="image_17.png" /></figure>
<figure><img src="image_18.png" /><figcaption>Figure 10 : milieu bidimensionnel homogène, conductivité thermique fonction de la position, résistances thermiques en parallèle. Champ de conductivité thermique dans l'échantillon d'étude (échelle de couleurs : 1 W/m/K et 19 W/m/K).</figcaption></figure>
Figure 11 : milieu bidimensionnel homogène, conductivité thermique fonction de la position, résistances thermiques en parallèle. Champ de température obtenu par calcul dans l'échantillon d'étude (échelle de couleurs : 0 1).

<figure><img src="image_19.png" /></figure>
<figure><img src="image_20.png" /><figcaption>Figure 12 : profil du champ de température de la figure 11 en fonction de la variable x.</figcaption></figure>
<figure><img src="image_21.png" /><figcaption>Figure 13 : milieu bidimensionnel homogène, conductivité thermique fonction de la position, résistances thermiques en série. Champ de conductivité thermique dans l'échantillon d'étude (échelle de couleurs : 1 W/m/K et 19 W/m/K).</figcaption></figure>
Figure 14 : milieu bidimensionnel homogène, conductivité thermique fonction de la position, résistances thermiques en série. Champ de température obtenu par calcul dans l'échantillon d'étude (échelle de couleurs : 0 1).

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

<figure><img src="image_23.png" /></figure>
Figure 15 : profil du champ de température de la figure 14 en fonction de la variable x (courbe en trait continu, échelle de gauche) comparé à la courbe du champ de température théorique (courbe en trait pointillé, échelle de droite).

<figure><img src="image_24.png" /></figure>
<figure><img src="image_25.png" /><figcaption>Figure 16 : milieu bidimensionnel homogène, conductivité thermique fonction de la position, résistances thermiques "en produit série-parallèle". Champ de conductivité thermique dans l'échantillon d'étude (échelle de couleurs : 1 W/m/K 190 W/m/K).</figcaption></figure>
Figure 17 : milieu bidimensionnel homogène, conductivité thermique fonction de la position, résistances thermiques "en produit série-parallèle". Champ de température obtenu par calcul dans l'échantillon d'étude (échelle de couleurs : 0 à 1).

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

<figure><img src="image_27.png" /></figure>
<figure><img src="image_28.png" /><figcaption>Figure 18 : profil du champ de température de la figure 17 en fonction de la variable x (courbe en trait continu, échelle de gauche) comparé à l'expression analytique (41) (courbe en trait pointillé, échelle de droite).</figcaption></figure>
<figure><img src="image_29.png" /><figcaption>Figure 19 : milieu à deux constituants et périodique. Représentation sous forme d'images de quelques-uns des dix échantillons étudiés : (a) la maille génératrice du milieu périodique ; (b) échantillon constitué de 2 × 2 mailles ; (c) échantillon constitué de 3 × 3 mailles ; (d) échantillon constitué de 5 × 5 mailles ; (e) échantillon constitué de 9 × 9 mailles.</figcaption></figure>
<figure><img src="image_30.png" /><figcaption>Figure 20 : milieu à deux constituants et périodique. Représentation sous forme d'images des champs de température s'établissant dans les cinq échantillons de la figure 19 (échelle de couleurs : 0 à 1 ; les isothermes sont tracées par pas de 0.1).</figcaption></figure>
<figure><img src="image_31.png" /><figcaption>Figure 21 : milieu à deux constituants et périodique. Représentation sous forme d'images des champs de déviation thermique $\delta\theta$ dans les cinq échantillons de la figure 19 (échelles de couleurs : $\delta\theta_{min}$ et $\delta\theta_{max}$ ; les valeurs de $\delta\theta_{min}$ et $\delta\theta_{max}$ sont précisées pour chaque image).</figcaption></figure>
<figure><img src="image_32.png" /><figcaption>Figure 22 : variations, en fonction du nombre de mailles dans chacune des directions horizontale et verticale de l'échantillon, de la norme infinie du champ de déviation thermique $\delta\theta$ (courbe en trait continu, échelle de gauche) et de la conductivité thermique effective moyenne $\bar{\Lambda}$ (courbe en trait pointillé, échelle de droite).</figcaption></figure>
Figure 23 : milieu à deux constituants et périodique. Représentation sous forme d'images en teintes de gris des spectres bidimensionnels (a) de l'échantillon "10 x 10" et (b) du champ de déviation thermique δθ obtenu sur cet échantillon (pour rendre les autres pics visibles, nous avons, sur chacune des deux images présentées ci-dessous, atténué le pic central - correspondant aux fréquences spatiales nulles - de sorte que son amplitude égale celle des pics immédiatement inférieurs).

<figure><img src="image_33.png" /><figcaption>(a)</figcaption></figure>

<figure><img src="image_34.png" /><figcaption>(b)</figcaption></figure>
<figure><img src="image_35.png" /><figcaption>Figure 24 : milieu à deux constituants et périodique. Représentation sous forme d'images des champs de déviation thermique δθ moyennés sur toutes les mailles constitutives des dix échantillons étudiés puis multipliés par les nombres Pₓ correspondants (échelles de couleurs : min max ; les valeurs de min et max sont précisées pour chaque image).</figcaption></figure>
Figure 25 : (a) maille génératrice du milieu bidimensionnel périodique réel ; (b) maille génératrice du milieu bidimensionnel périodique fictif. Echelles de teintes de gris : 1 W/m/K 5 W/m/K pour les deux images.

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

<figure><img src="image_37.png" /></figure>
<figure><img src="image_38.png" /><figcaption>Figure 26 : profil du champ de température θ(x) issu de l'analyse spectrale (courbe en trait continu, échelle de gauche) comparé à l'expression analytique (41) (courbe en trait pointillé, échelle de droite).</figcaption></figure>
<figure><img src="image_39.png" /><figcaption>Figure 27 : variations, en fonction du niveau de filtrage passe-bas M, de la conductivité thermique effective $\Lambda$ déduite de l'analyse spectrale pour le milieu étudié au chapitre 3-1-2. La ligne horizontale en trait pointillé marque la valeur théorique de $\Lambda$ (qui est $\sqrt{19}$).</figcaption></figure>
<figure><img src="image_40.png" /><figcaption>Figure 28 : champ de conductivité thermique λ(x) retenu pour l'étude des performances de l'analyse spectrale sur les milieux monodimensionnels périodiques à variations brutales de propriétés. La partie de courbe en trait épais matérialise la période spatiale du champ λ(x).</figcaption></figure>
<figure><img src="image_41.png" /><figcaption>Figure 29 : champ de température s'établissant dans le milieu monodimensionnel périodique décrit à la figure 28. Profil du champ de température issu de l'analyse spectrale (courbe en trait continu, échelle de gauche) comparé à la solution analytique du problème (courbe en trait pointillé, échelle de droite).</figcaption></figure>
Figure 30 : différentes courbes de déviation thermique δθ calculées au moyen de l'analyse spectrale sur une période spatiale du milieu monodimensionnel périodique décrit à la figure 28 pour différentes valeurs du niveau de filtrage passe-bas M. La solution analytique du champ de déviation thermique δθ est également reportée sur ce graphe. Pour une plus grande clarté de la figure, les trois courbes correspondant à M = 5, M = 10 et M = 50 ont été décalées respectivement de +0.03, +0.02 et +0.01 dans la direction verticale.

<figure><img src="image_42.png" /></figure>
<figure><img src="image_43.png" /><figcaption>Figure 31 : variations, en fonction du niveau de filtrage passe-bas M, de la conductivité thermique effective $\Lambda$ déduite de l'analyse spectrale pour le milieu monodimensionnel périodique décrit à la figure 28. La ligne horizontale en trait pointillé marque la valeur théorique de $\Lambda$ (qui est $\frac{1}{0.84}$).</figcaption></figure>
<figure><img src="image_44.png" /><figcaption>Figure 32 : milieu bidimensionnel périodique, conductivité thermique fonction de la position. Champ de conductivité thermique dans l'échantillon d'étude (échelle de couleurs : 1 W/m/K 5 W/m/K).</figcaption></figure>

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

<figure><img src="image_46.png" /><figcaption>Profil de conductivité thermique sur le segment [AB] :</figcaption></figure>

<figure><img src="image_47.png" /></figure>
Figure 33 : milieu bidimensionnel périodique, conductivité thermique fonction de la position. Représentation sous forme d'images du champ de déviation thermique $\delta\theta$ obtenu par calcul dans l'échantillon décrit à la figure 32 : (a) champ $\delta\theta$ sur l'échantillon complet ; (b) champ $\delta\theta$ moyenné sur les 9 périodes de l'échantillon. Echelles de couleurs : $\delta\theta_{min}$ à $\delta\theta_{max}$ ; les valeurs de $\delta\theta_{min}$ et $\delta\theta_{max}$ sont précisées pour chaque image.

<figure><img src="image_48.png" /><figcaption>(a) ; $\delta\theta_{min} = -0.046$ et $\delta\theta_{max} = +0.044$</figcaption></figure>

<figure><img src="image_49.png" /><figcaption>(b) ; $\delta\theta_{min} = -0.045$ et $\delta\theta_{max} = +0.044$</figcaption></figure>
Figure 34 : milieu bidimensionnel périodique, conductivité thermique fonction de la position. Comparaison des champs de déviation thermique δθ calculés par différences finies et par analyse spectrale sur la période spatiale décrivant l'échantillon de la figure 32 : (a) champ δθ moyen calculé par différences finies ; (b) champ δθ calculé par analyse spectrale ; (c) différence entre les deux champs. Echelles de couleurs : δθmin à δθmax ; les valeurs de δθmin et δθmax sont précisées pour chaque image.

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

(a) ; δθmin = -0.045 et δθmax = +0.044

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

(b) ; δθmin = -0.045 et δθmax = +0.045

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

(c) ; δθmin = 0.000 et δθmax = +0.001
<figure><img src="image_53.png" /><figcaption>Figure 35 : milieu bidimensionnel périodique à deux constituants, géométrie de type nid d'abeille ; (a) rappel de l'allure de la maille génératrice du milieu ; (b) champ de déviation thermique $\delta\theta$ moyen calculé par différences finies sur l'échantillon "10 × 10" ; (c1) champ de conductivité thermique $\lambda$ résultant d'un filtrage passe-bas caractérisé par $M = N = 10$ ; (c2) champ de déviation thermique $\delta\theta$ calculé par analyse spectrale avec $M = N = 10$ ; (c3) différence entre les deux champs de déviation thermique $\delta\theta$ des figures (c2) et (b) ; (d1) champ de conductivité thermique $\lambda$ résultant d'un filtrage passe-bas caractérisé par $M = N = 20$ ; (d2) champ de déviation thermique $\delta\theta$ calculé par analyse spectrale avec $M = N = 20$ ; (d3) différence entre les deux champs de déviation thermique $\delta\theta$ des figures (d2) et (b). Pour une meilleure lisibilité des images, celles-ci ont été agrandies dans un rapport 2.</figcaption></figure>