# Conductivité thermique des composites à fort contraste de phase : approche théorique et numérique

Jean-Pierre Vassal, Laurent Orgéas, Steven Le Corre, Denis Favier, Jean-Louis Auriault

► To cite this version:

Jean-Pierre Vassal, Laurent Orgéas, Steven Le Corre, Denis Favier, Jean-Louis Auriault. Conductivité thermique des composites à fort contraste de phase : approche théorique et numérique. MATERIAUX 2006: "Nano Mat - Tech: des NanoMatériaux aux NanoSystèmes.", Nov 2006, Dijon, France. hal-01008016

HAL Id: hal-01008016

https://hal.science/hal-01008016

Submitted on 9 May 2023

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
Conductivité thermique des composites à fort contraste de phase
Approche théorique et numérique

J.-P. Vassal<sup>a</sup>, L. Orgéas<sup>a</sup>, S. Le Corre<sup>b</sup>, D. Favier<sup>a</sup>, J.-L. Auriault<sup>a</sup>

<sup>a</sup> Laboratoire Sols-Solides-Structures (3S)
CNRS – Universités de Grenoble (INPG – UJF)
BP 53, 38041 Grenoble cedex 9, France
Jean-Pierre.Vassal@hmg.inpg.fr, Laurent.Orgeas@hmg.inpg.fr, Denis.Favier@hmg.inpg.fr,
Jean-Louis.Auriault@hmg.inpg.fr

GeM – Institut de Recherche en Génie Civil et Mécanique
CNRS – Ecole Centrale de Nantes
BP 92101, 44321 Nantes cedex 3, France
steven.le-corre@ec-nantes.fr

RESUME

Cette étude propose une modélisation des phénomènes de conduction thermique dans les composites à forte concentration de renforts fibreux et présentant une conductivité de la matrice négligeable devant celle des renforts. La méthode d'homogénéisation des structures périodiques par développements asymptotiques est utilisée pour déterminer analytiquement la forme des équations de transfert thermique en régime transitoire à l'échelle macroscopique. La nature élancée des renforts permet l'obtention d'une formulation discrète du problème. On génère ensuite numériquement des milieux fibreux pour étudier l'influence de différents paramètres microstructuraux à l'aide d'un code de type éléments discrets. Les résultats obtenus montrent une assez bonne adéquation avec des mesures expérimentales de conductivité issues de la littérature.

MOTS-CLES : conductivité, thermique, composite, fibres, homogénéisation, éléments discrets

I. INTRODUCTION

Augmenter la conductivité thermique des composites à matrice polymère permet d'envisager des applications industrielles intéressantes telles que des dissipateurs thermiques injectés, des composants électroniques s'échauffant moins, des systèmes de frein plus légers... En effet, charger une matrice polymère de particules conductrices telles que du graphite [1], du noir de carbone [2], des poudres métalliques [3], des fibres d'aluminium [4], ou de carbone [5-6] sont autant de solutions qui permettent d'accroître notablement la conductivité effective du composite. De telles propriétés thermiques s'associent aux qualités usuelles des matériaux composites à matrice polymère : faible masse, bonnes propriétés mécaniques, résistance à la corrosion, mise en œuvre aisée... Néanmoins la difficulté pour prédire la conductivité thermique de ces matériaux limite leur développement.

De très nombreux modèles ont déjà été proposés dans la littérature pour modéliser les propriétés thermiques des matériaux composites [7-16]. Ces derniers sont en règle générale peu adaptés aux matériaux à renforts fibreux considérés ici, constitués de fibres très élancées en très forte concentration, avec de possibles résistances de contact aux interfaces fibre-fibre ou fibre-matrice. En outre, ces modèles sont très souvent développés pour des situations de régime permanent (conduction pure), dans le cas linéaire (propriétés thermiques locales indépendantes du champ thermique), et en présupposant le milieu continu équivalent comme monophasique. Dans ce travail, on propose une modélisation des phénomènes de conduction transitoires et non linéaires dans les composites conducteurs thermiques à renforts fibreux présentant d'une part une conductivité thermique de la matrice négligeable devant celle des renforts et d'autre part des résistances de contact aux interfaces fibre-fibre. Pour cela, la méthode d'homogénéisation des structures périodiques par développements asymptotiques est utilisée [17-20]. La première partie (§II) rappelle les résultats théoriques obtenus par cette modélisation : trois modèles sont ainsi distingués en fonction de la qualité des contacts fibre-fibre (§II.A). Les deux premiers sont simplifiés pour obtenir une formulation discrète du problème thermique (§II.B), alors que le troisième est déjà discret sans simplification.
La validité des modèles discrets est ensuite étudiée (§III). Ces modèles, introduits dans un code de type éléments discrets, permettent d'étudier sur des milieux fibreux à fort enchevêtrement générés numériquement (§IV.A), l'influence des paramètres microstructuraux sur les propriétés thermiques effectives (§IV.B). Enfin, les résultats obtenus par simulation numérique sont comparés avec des résultats expérimentaux (§VI) issus de la littérature [5-6].

## II. FORMULATION DU PROBLEME THERMIQUE

### A. Description du problème local

Le composite thermo-conducteur à renforts fibreux est vu comme une répétition périodique d'un Volume Elémentaire Représentatif (VER) de volume $\Omega$ et de longueur caractéristique $l_c$. Le VER est constitué d'un ensemble $\mathcal{P}$ de P fibres $p_\alpha$ de volume $\Omega_\alpha$, de longueurs $l_\alpha$ telles que $l_\alpha < l_c$, et de sections droites $S_\alpha$ telles que $S_\alpha \ll l_\alpha^2$. Elles ont une conductivité thermique caractérisée par un tenseur de conductivité thermique $\Lambda_\alpha$, dont la valeur caractéristique des conductivités principales est notée $\Lambda_c$, et une capacité thermique volumique $c_\alpha$ de valeur caractéristique $c_c$. Elles sont immergées dans une matrice ayant une conductivité thermique considérée comme nulle au regard de celle des fibres. En outre, on suppose que le milieu est suffisamment concentré de sorte que chaque fibre est en contact avec au moins deux autres fibres. L'ensemble des $C$ contacts fibre-fibre du VER est noté $\mathcal{C}$ et l'ensemble des $C_\alpha$ contacts appartenant à la fibre $p_\alpha$ est noté $\mathcal{C}_\alpha$. La surface extérieure des fibres est divisée en $\Gamma_\alpha$, surface d'échange entre les fibres et la matrice, et $\Gamma_{i\alpha\beta}$ (l'indice i représente le $i^{ème}$ contact dans le VER), surfaces de contact de la fibre $p_\alpha$ avec les fibres $p_\beta$ (dimension caractéristique $\Gamma_c$). La physique retenue à l'échelle locale suppose que l'équilibre thermique est régi par une équation de la chaleur qui en tout point d'une fibre $p_\alpha$ s'écrit :

$$c_{\alpha} \dot{T}_{\alpha} = -\operatorname{div} \mathbf{q}_{\alpha} + r_{\alpha}, \qquad (1)$$

où $T_\alpha(\mathbf{X}, t)$ est la température au temps $t$ et au point considéré de position $\mathbf{X}$, $r_\alpha(\mathbf{X}, t)$ représente une source de chaleur volumique, et où le flux de chaleur $\mathbf{q}_\alpha$ suit une loi de Fourier :

$$\mathbf{q}_{\alpha} = -\Lambda_{\alpha} \operatorname{grad} T_{\alpha}, \qquad (2)$$

On suppose que le flux est nul sur $\Gamma_\alpha$, car la matrice est supposée isolante. En revanche un flux de chaleur proportionnel au coefficient de transfert thermique $h_{i\alpha\beta}$ de valeur caractéristique $h_c$ existe sur les surfaces de contact $\Gamma_{i\alpha\beta}$ :

$$\left\{ \begin{array}{ll} \mathbf{q}_{\alpha} \cdot \mathbf{n}_{\alpha} = 0 & \forall M \in \Gamma_{\alpha} \\ \mathbf{q}_{\alpha} \cdot \mathbf{n}_{i\alpha\beta} = \mathbf{q}_{\beta} \cdot \mathbf{n}_{i\alpha\beta} & \forall M \in \Gamma_{i\alpha\beta}, \\ \mathbf{q}_{\alpha} \cdot \mathbf{n}_{i\alpha\beta} = -h_{i\alpha\beta}(T_{\beta} - T_{\alpha}) = -h_{i\alpha\beta}\Delta T_{\alpha\beta} & \forall M \in \Gamma_{i\alpha\beta} \end{array} \right. \qquad (3)$$

où les vecteurs $\mathbf{n}$ représentent les vecteurs unitaires normaux aux surfaces considérées. A l'instant initial $t_0$, les températures en tout point sont identiques. Les températures en tout point sont considérées du même ordre de grandeur. Les valeurs de $\Lambda_\alpha$, $c_\alpha$ et de $h_{i\alpha\beta}$ peuvent dépendre de la variable d'espace $\mathbf{X}$ mais aussi de la température $T_\alpha$. Enfin, on suppose qu'il existe une bonne séparation d'échelle entre la longueur caractéristique « macroscopique » $L_c$ (dimension géométrique caractéristique de l'échantillon et/ou longueur d'onde de la sollicitation) et la longueur caractéristique « microscopique » $l_c$, ce qui permet d'introduire le paramètre de séparation d'échelles $\epsilon = l_c / L_c \ll 1$. L'ensemble des hypothèses et équations précédentes forment le problème local.

### B. Changement d'échelle

Le comportement macroscopique équivalent au problème thermique hétérogène précédent a été étudié par la méthode d'homogénéisation des structures périodiques par développements asymptotiques à échelles multiples [18-19] en suivant la méthode proposée en [20] qui se scinde en quatre étapes. La première étape consiste à adimensionnaliser le problème local afin de quantifier l'importance des différents phénomènes physiques. Dans notre cas [17], on fait ainsi apparaître trois nombres adimensionnels, i.e. un nombre de Biot $B$ qui quantifie l'importance des échanges aux contacts par rapport à la conduction dans les fibres, un nombre de Fourier $F$ qui quantifie l'importance des phénomènes capacitifs dans les fibres par rapport aux phénomènes de conduction, et un nombre $R$ qui quantifie l'importance des effets capacitifs vis à vis des sources de chaleur extérieures:

$$B = \frac{h_c l_c}{\Lambda_c}, \quad F = \frac{t_c \Lambda_c}{c_c l_c^2}, \quad R = \frac{t_c r_c}{c_c}. \qquad (4)$$

Il est possible de montrer que les situations homogénéisables correspondent à des nombres $F$ et $R$ d'ordre $\epsilon^{-2}$ [17]. De plus, trois comportements distincts apparaissent en fonction du nombre de Biot $B$. Si $B$ est
d'ordre $\epsilon^{-1}$, la physique est pilotée par la conduction dans les fibres (modèle I) ; s'il est d'ordre $\epsilon^0$, les contacts fibres-fibres et la conduction dans les fibres ont la même importance sur la conductivité globale du milieu (modèle II) ; s'il est d'ordre $\epsilon^1$, la physique est pilotée par les contacts (modèle III). La seconde étape consiste à développer les champs de température $T_\alpha$ et les sources de chaleur $r_\alpha$ sous formes de développements asymptotiques en puissances de $\epsilon$:

$$ 
\begin{aligned}
T_{\alpha}(\mathbf{X},t) &= T_{\alpha}^{[0]}(\mathbf{x},\mathbf{y},t) + \epsilon T_{\alpha}^{[1]}(\mathbf{x},\mathbf{y},t) + \epsilon^2 T_{\alpha}^{[2]}(\mathbf{x},\mathbf{y},t) + \dots \\
r_{\alpha}(\mathbf{X},t) &= r_{\alpha}^{[0]}(\mathbf{x},\mathbf{y},t) + \epsilon r_{\alpha}^{[1]}(\mathbf{x},\mathbf{y},t) + \epsilon^2 r_{\alpha}^{[2]}(\mathbf{x},\mathbf{y},t) + \dots
\end{aligned} 
\quad (5) $$

où $\mathbf{x} = \mathbf{X}/L_c$ est la variable d'espace macroscopique, $\mathbf{y} = \mathbf{X}/l_c$ la variable d'espace microscopique, et où les champs de température $T_\alpha^{[i]}$ (et les sources $r_\alpha^{[i]}$) sont du même ordre de grandeur et périodiques en $\mathbf{y}$ sur $\Omega$. La troisième étape consiste à introduire le développement (5) dans le problème local (3), puis à identifier les termes de même puissance de $\epsilon$. Enfin, la quatrième étape consiste à résoudre les différents problèmes aux limites qui apparaissent aux différents ordres de $\epsilon$. Cette méthode permet de montrer [17] : (i) que $T^{[0]}$ est seulement fonction de la variable macroscopique $\mathbf{x}$ (i.e. $T^{[0]}$ est constante sur le VER) ainsi le milieu macroscopique équivalent est toujours monophasique, (ii) que son équilibre est régi par une équation de la chaleur standard :

$$ \langle c_{\alpha}(T^{[0]}) \rangle \dot{T}^{[0]} = -\operatorname{div} \mathbf{q}_{\alpha}^{\text{mac}} + \langle r_{\alpha}^{[0]} \rangle, \quad (6) $$

où $\langle \cdot \rangle$ représente la moyenne volumique sur le VER de la quantité $a$, (iii) que le flux de chaleur macroscopique $\mathbf{q}_\alpha^{\text{mac}}$ suit toujours une loi de Fourier de la forme

$$ \mathbf{q}_{\alpha}^{\text{mac}} = -\Lambda^{\text{eff}}(T^{[0]}) \cdot \text{grad} \, T^{[0]}, \quad (7) $$

où $\Lambda^{\text{eff}}$ est le tenseur de conductivité macroscopique, (iv) qu'en fonction du nombre de Biot $B$, trois expressions différentes de $\Lambda^{\text{eff}}$ apparaissent :

* Modèle I, physique dominée par la conduction dans les fibres $B=\mathcal{O}(\epsilon^{-1})$:

$$ \Lambda^{\text{eff}}(T^{[0]}) = \langle \Lambda_{\alpha}(T^{[0]}) \rangle + \langle \Lambda_{\alpha}(T^{[0]}) \cdot \text{grad} \, \epsilon \, \boldsymbol{\theta}_{\alpha}^{[1]} \rangle \quad (8) $$

où les composantes des vecteurs $\epsilon \boldsymbol{\theta}_\alpha^{[1]}$ correspondent aux fluctuations de températures $\epsilon T_\alpha^{[1]}$ lorsque le gradient de température $\text{grad} \, T^{[0]} = \mathbf{e}_k$ ($k=1,2,3$). Les fluctuations de température sont d'une manière générale solutions du problème de localisation suivant, posé sur le VER:

$$ 
\left\{
\begin{array}{ll}
\text{div} \, \mathbf{q}_\alpha^{[1]} = 0 & \text{dans } \Omega_\alpha \\
\mathbf{q}_\alpha^{[1]} = -\Lambda_\alpha(T^{[0]}) \cdot (\text{grad} \, \epsilon T_\alpha^{[1]} + \text{grad} \, T^{[0]}) & \text{dans } \Omega_\alpha \\
\mathbf{q}_\alpha^{[1]} \cdot \mathbf{n}_\alpha = 0 & \text{sur } \Gamma_\alpha, \\
\epsilon T_\alpha^{[1]} = \epsilon T_\beta^{[1]} & \text{sur } \Gamma_{i\alpha\beta} \\
\mathbf{q}_\alpha^{[1]} \cdot \mathbf{n}_{i\alpha\beta} = \mathbf{q}_\beta^{[1]} \cdot \mathbf{n}_{i\alpha\beta} & \text{sur } \Gamma_{i\alpha\beta}
\end{array}
\right.
\quad (9) $$

où la température $T^{[0]}$ et son gradient $\text{grad} \, T^{[0]}$ sont donnés. Ainsi, pour ce modèle, les fluctuations de température et par conséquent la conductivité thermique effective ne dépend pas des coefficients de transfert thermique $h_{i\alpha\beta}$.

* Modèle II, physique mixte – La forme du tenseur de conductivité thermique effective est similaire à celle du modèle I. Cependant, l'équation (9d) est remplacée par :

$$ \mathbf{q}_{\alpha}^{[1]} \cdot \mathbf{n}_{i\alpha\beta} = -h_{i\alpha\beta}(T^{[0]}) \Delta\epsilon T_{\alpha\beta}^{[1]} \quad \text{sur } \Gamma_{i\alpha\beta}, \quad (10) $$

Ainsi les fluctuations de température et la conductivité thermique effective dépendent aussi des coefficients de transfert thermique $h_{i\alpha\beta}$.

* Modèle III, physique dominée par les échanges aux contacts – Le tenseur de conductivité thermique effective s'écrit sous une forme discrète :

$$ \Lambda^{\text{eff}}(T^{[0]}) = \frac{h_{i\alpha\beta}(T^{[0]})}{\Omega} \sum_{\mathbf{C}} \Gamma_{i\alpha\beta} \mathbf{X}_{\alpha\beta} \otimes (\mathbf{X}_{\alpha\beta} + \Delta\epsilon \bar{\boldsymbol{\theta}}_{\alpha\beta}^{[1]}), \quad (11) $$

où $\mathbf{e}_k$ vecteur reliant le centre de gravité des fibres $p_\alpha$ et $p_\beta$ est noté $\mathbf{X}_{\alpha\beta} = \mathbf{G}_\alpha \mathbf{G}_\beta$ et où les vecteurs $\Delta\boldsymbol{\theta}_{\alpha\beta}^{[1]}$ sont solutions du système

$$ \sum_{\mathbf{C}_\alpha} h_{i\alpha\beta}(T^{[0]}) \Gamma_{i\alpha\beta} \left( \Delta\epsilon \bar{\boldsymbol{\theta}}_{\alpha\beta}^{[1]} + \mathbf{X}_{\alpha\beta} \right) = \mathbf{0}. \quad (12) $$

Ici, le tenseur de conductivité thermique effectif ne dépend pas des tenseurs de conductivité $\Lambda_\alpha$.
Il est important de noter que si la physique locale envisagée est transitoire et non linéaire, les problèmes de localisation nécessaires aux calculs des tenseurs de conductivité effective sont linéaires et statiques. On peut également remarquer que dans les modèles I et II, la détermination du tenseur de conductivité thermique effective nécessite la détermination des fluctuations de température en tout point du VER par la résolution du système d'équations différentielles (9). En revanche, pour le modèle III, la détermination de ce même tenseur ne nécessite que la résolution d'un système discret de P équations linéaires (12). On peut enfin noter que quel que soit le modèle, la capacité thermique volumique macroscopique et la source de chaleur macroscopique sont de simples moyennes en volume (6), leur détermination étant complètement découplée du calcul de la conductivité thermique effective.

### C. Modèles I et II discrets

Pour simplifier les modèles I et II, on peut tenir compte de la nature élancée des fibres (cf. figure 1a). Pour simplifier, nous nous limitons ici au cas des fibres droites de vecteurs directeur $\mathbf{e}_{\alpha}$ de section droite $S$, de longueur $l$, de conductivité longitudinale $\Lambda$ et qui ont en leurs points de contact un coefficient de transfert $h$. On peut également faire l'hypothèse que loin des points de contact les flux de chaleur dans les fibres sont (i) constants dans chaque portion limitée par les points de contact et (ii) portés par l'axe des fibres. En utilisant l'équation (9b), il est possible d'approcher l'expression du flux de chaleur entre deux points de contact $M_{i\alpha}$ et $M_{j\alpha}$ appartenant à la même fibre $p_\alpha$ selon l'expression discrète suivante :

$$ \mathbf{q}_{i\alpha} \approx -\Lambda \left( \frac{\Delta\varepsilon T_{i\alpha}^{[1]}}{l_{i\alpha}} + \mathbf{e}_{\alpha} \cdot \mathbf{grad} \, T^{[0]} \right) \mathbf{e}_{\alpha}, \quad (13) $$

où $\Delta\varepsilon T_{i\alpha}^{[1]}$ est la différence de fluctuations de température entre $M_{i\alpha}$ et $M_{j\alpha}$ et où $l_{i\alpha}$ est la distance entre ces deux points (cf. figure 1a). Pour chaque contact entre deux fibres $p_\alpha$ et $p_\beta$, le flux de chaleur est porté par la perpendiculaire commune aux deux fibres de vecteur unitaire $\mathbf{e}_{i\alpha\beta} = \mathbf{e}_\alpha \wedge \mathbf{e}_\beta$. En utilisant l'équation (10), on peut écrire :

$$ \mathbf{q}_{i\alpha\beta} \approx -h\Delta\varepsilon T_{i\alpha\beta}^{[1]} \mathbf{e}_{i\alpha\beta}, \quad (14) $$

où $\Delta\varepsilon T_{i\alpha\beta}^{[1]}$ est la différence de fluctuations de température entre le point $M_{i\alpha}$ situé sur l'axe de la fibre $p_\alpha$ en face du contact et le point $M_{j\beta}$ sur l'axe de la fibre $p_\beta$ situé en face de ce même contact (la distance entre ces points est notée $d$). On peut montrer [17] que ces hypothèses permettent de simplifier la résolution du système d'équations différentielles (13) à la résolution du système linéaire à 2C équations :

$$ \Lambda S \left( \frac{\Delta\varepsilon T_{i\alpha}^{[1]}}{l_{i\alpha}} + \frac{\Delta\varepsilon T_{i\alpha}^{[1]}}{l_{i\alpha}} + \delta \mathbf{e}_{\alpha} \cdot \mathbf{grad} \, T^{[0]} \right) + h \Gamma_{i\alpha\beta} \Delta T_{i\alpha\beta}^{[1]} = 0, \quad (15) $$

où les points $M_{i\alpha}$ et $M_{j\alpha}$ sont les points de contact sur la fibre $p_\alpha$ situés « à gauche » et « à droite » de $M_{i\alpha}$ (cf. figure 1a), et où le paramètre $\delta$ vaut 1 pour le premier contact sur la fibre $p_\alpha$ -1 pour le dernier et 0 pour les autres. Les fluctuations $\varepsilon T_{i\alpha}^{[1]}$ sont les 2C inconnues à déterminer pour trois gradients de température macroscopiques $\mathbf{grad} \, T^{[0]} = \mathbf{e}_k$ ($k = 1,2,3$). Elles forment les composantes des vecteurs que l'on note $\varepsilon \mathbf{\Theta}_{i\alpha}^{[1]}$. Le modèle discret ainsi obtenu ne dépend plus de la fluctuation de température en tout point du VER mais seulement de la fluctuation de température proche des points de contact. Les hypothèses précédentes permettent de réécrire le tenseur de conductivité thermique selon [17] :

$$ \Lambda^{\text{eff}} = \Lambda^{\text{fibre}} + \Lambda^{\text{contact}}, \quad (16) $$

avec

$$ \Lambda^{\text{fibre}} = f \Lambda \left( 1 - \frac{1}{\bar{C}_{\alpha} + 1} \right) \frac{1}{P} \sum_{\mathbf{P}} k_{\alpha} \mathbf{e}_{\alpha} \otimes \left( \mathbf{e}_{\alpha} + \frac{\Delta\varepsilon \mathbf{\Theta}_{i\alpha}^{[1]} + j^{+} \alpha}{k_{\alpha} l^{-}} \right) \quad \text{et} \quad \Lambda^{\text{contact}} = 2\chi d S h \frac{1}{C} \sum_{\mathbf{C}} k_{i\alpha\beta} \mathbf{e}_{i\alpha\beta} \otimes \Delta\varepsilon \mathbf{\Theta}_{i\alpha\beta}^{[1]} \quad (17) $$

où l'indice $i$ repère le premier contact sur la fibre $p_\alpha$ et $j$ le dernier, $l$ est la moyenne des longueurs $l_\alpha$ entre le premier et le dernier contact sur les fibres, $k_\alpha = l_\alpha^{-} / l^{-}$, $\bar{C}_\alpha$ le nombre moyen de contacts par fibre, $f = \text{PSI} / \Omega$ la fraction volumique de fibres, $\chi = C / \Omega$ le nombre de connexions par unité de volume, et où $k_{i\alpha\beta} = \Gamma_{i\alpha\beta} / S$. Il est à noter que dans le cas du modèle I, les vecteurs $\Delta\varepsilon \mathbf{\Theta}_{i\alpha\beta}^{[1]}$ sont nuls si bien que $\Lambda^{\text{contact}} = 0$. Il est à noter aussi que pour de très mauvais contacts thermiques, l'expression (16) de $\Lambda^{\text{eff}}$ tend vers celle établie pour le modèle III, i.e. (11).

### III. MICROSTRUCTURES FIBREUSES ELEMENTAIRES

Dans ce paragraphe, on montre sur des microstructures fibreuses élémentaires les limites de validité des modèles discrets établis auparavant. Pour cela, on choisit des milieux fibreux très simples dont le VER est représenté sur la figure 1b. Il est constitué d'une boîte prismatique de dimensions $l_1, l_2$ et $l_3$ et de volume $\Omega$ contenant une matrice isolante et deux fibres droites $p_\alpha$ et $p_\beta$ de section carrée de surface $S$, de propriétés
thermiques isotropes et de conductivité $\Lambda = 1 \text{ Wm}^{-1}\text{K}^{-1}$. Ces deux fibres sont en contact en deux zones $\Gamma_{\alpha\beta}$ de surfaces $S_{AB}$ dont le coefficient de transfert thermique vaut $h$. La périodicité est assurée au niveau des surfaces $\Gamma_{\alpha\gamma}$. Les parties « actives » des fibres jouant un rôle dans la conductivité thermique du milieu (i.e. les parties des fibres situées entre le premier et le dernier point de contact) ont pour longueur $I_{\alpha}^{-} = 3 \text{ mm}$ et $I_{\beta}^{-} = 5 \text{ mm}$. Le tenseur de conductivité macroscopique du milieu continu équivalent a été calculé de deux façons :

(i) un calcul dit « continu » pour le modèle II, consistant à résoudre par éléments finis le problème de localisation (9). Pour cela, le code éléments finis FemLab 3.1 a été utilisé (cf. maillage tétraédrique sur la figue 1b contenant environ 60000 éléments quadratiques).

(ii) un calcul dit « discret » pour les modèles I, II et III, consistant à résoudre analytiquement les systèmes d'équations linéaires (12) ou (15).

<figure><img src="image_2.png" /><figcaption>Figure 1 : (a) Schéma d'une microstructure fibreuse à fibres droites, (b) VER élémentaire, (c) Evolution de la conductivité principale macroscopique en fonction du coefficient de transfert thermique pour les différents modèles proposés.</figcaption></figure>

## A. Continuité des modèles

La figure 1c donne à titre d'exemple l'évolution de la seule composante principale non nulle du tenseur de conductivité effectif pour différentes valeurs de coefficient d'échange thermique $h$, pour des fibres de section $S=1 \text{ mm}^2$ et pour $\Omega=5*5*4 \text{ mm}^3$. Lorsque $h$ croît, la conductivité thermique effective augmente plutôt linéairement puis tend progressivement vers une asymptote horizontale. On peut remarquer qu'il y a continuité des trois modèles et que le modèle II, plus riche, est applicable quelle que soit la valeur du coefficient de transfert thermique. Le modèle III reproduit fidèlement l'évolution linéaire pour les faibles valeurs de $h$. Le modèle I reproduit quant à lui fidèlement l'asymptote horizontale pour les fortes valeurs de $h$. Enfin, on constate un faible décalage (environ 10% d'erreur dans cet exemple) entre les valeurs données par le modèle II discret et le modèle II continu correspondant à l'erreur due aux hypothèses de simplification à la base du modèle discret (cf. § suivant).

## B. Limites de validité des modèles discrets

Les hypothèses formulées pour réduire les modèles continus I et II en modèles discrets sont en principe d'autant plus valides que les portions de fibres entre deux contacts sont élancées. Pour étudier quantitativement les limites des approches discrètes, des VER avec un élancement $el = I_{\beta} / \sqrt{S}$ variant de 2.5 à 100 ont été sollicités. Pour chacun d'eux, six valeurs de coefficient d'échange ont été étudiées variant de $h = 10^0 \text{ Wm}^{-2}\text{K}^{-1}$ à $h \to \infty$. Les résultats du modèle continu II et du modèle discret II sont alors comparés en évaluant l'erreur liée au modèle discret :

$$ \text{Érreur} = \frac{\Lambda_{I continu}^{\text{eff}} - \Lambda_{I discret}^{\text{eff}}}{\Lambda_{I continu}^{\text{eff}}}. \qquad (18) $$

Les résultats sont présentés figure 2. L'erreur due au modèle II discret est raisonnable. Elle est comprise entre $3.10^{-4}\%$ et 25% en fonction des élancements et des qualités des contacts étudiés. Elle décroît quasi linéairement lorsque l'élancement des fibres entre contacts augmente (figure 2a). Plus les fibres sont élancées, plus le flux est colinéaire à la fibre neutre des fibres et plus l'expression (13) est valide. En
revanche, l'erreur croît avec le coefficient d'échange thermique (figure 2b), dans un premier temps linéairement avant d'infléchir et d'atteindre une asymptote horizontale lorsque $h > 10^6 \text{Wm}^{-2}\text{K}^{-1}$. Cette asymptote correspond à des contacts assimilables à des contacts parfaits. Plus le coefficient de transfert thermique est grand plus l'apport de chaleur au niveau des surfaces de contact est important. Ce flux étant orthogonal à la fibre neutre, ce phénomène altère la validité de l'hypothèse (13).

<figure><img src="image_3.png" /><figcaption>Figure 2 : Evolution de l'erreur relative du modèle discret par rapport au modèle continu EF en fonction (a) de l'élancement des fibres entre deux contacts, et (b) du coefficient d'échange thermique $h$.</figcaption></figure>

## IV. MICROSTRUCTURES FIBREUSES « REALISTES »

L'objectif de cette partie est d'étudier en détail l'influence des différents paramètres microstructuraux sur la conductivité thermique effective en étudiant la réponse de réseaux fibreux discrets proches de ceux observés sur les composites conducteurs à fibres courtes.

### A. Génération de milieux fibreux

Les microstructures des matériaux fibreux étant très complexes et souvent difficilement observables, il est en général délicat d'obtenir une description précise de la microstructure. Ainsi, nous avons essayé de générer numériquement des microstructures fibreuses 3D simples les plus représentatives, comme le montrent les exemples de la figure 3. Ces réseaux fibreux contiennent P fibres droites cylindriques à section circulaire de rayon $a$ et de longueur $l$, le nombre $P$ étant piloté par la fraction volumique de fibres $f$ souhaitée dans une boîte de volume $\Omega$.

<figure><img src="image_4.png" /><figcaption>Figure 3 : Exemples de VER (a) isotrope, (b) isotrope transverse et (c) orientée suivant $e_1$.</figcaption></figure>

La position du centre d'inertie de chacune des fibres est ici aléatoire, ainsi que les composantes de leurs vecteurs directeurs $e_\alpha$. Les microstructures ainsi obtenues sont isotropes (figure 3a). Afin d'étudier l'influence de l'orientation des renforts, les vecteurs $e_\alpha$ sont ensuite modifiés pour permettre d'obtenir par exemple des microstructures isotropes transverses (fibres plutôt alignées dans le plan $(e_1, e_2)$, figure 3b), ou bien encore alignées dans une des directions du repère de référence (fibres plutôt alignées selon la direction $e_1$, figure 3c). Les distributions des orientations des fibres étant choisies sous forme de gaussiennes, on peut décrire en première approximation l'orientation des réseaux fibreux à l'aide du tenseur d'orientation du second ordre [21] (cf. figure 3):
$$
\mathbf{A} = \frac{1}{P} \sum_{\mathbf{P}} \mathbf{e}_{\alpha} \otimes \mathbf{e}_{\alpha} \qquad (19)
$$

L'orientation du réseau fibreux a été modulée à l'aide de deux critères : (i) l'anisotropie dans le plan (e₁,e₂) que l'on peut caractériser par le rapport A₁₁/A₂₂, et (ii) l'anisotropie plan (e₁,e₂) / épaisseur e₃ dont un descripteur possible est donné par le rapport ⟨A₁₁,A₂₂⟩/A₃₃, où ⟨A₁₁,A₂₂⟩ représente la moyenne des termes A₁₁ et A₂₂ : ⟨A₁₁,A₂₂⟩ = (A₁₁+A₂₂)/2.

La connectivité des réseaux, c'est-à-dire la position des *C* points de contact fibre-fibre présents dans le VER est calculée en utilisant une méthode déterministe inspirée du modèle de tube [22-24] : dans le cas présent, toute fibre *p*<sub>β</sub> qui intercepte la fibre *p*<sub>α</sub> est comptée dans les ensembles *A*<sub>α</sub> et *B*<sub>β</sub>. La surface de contact fibre-fibre est arbitrairement choisie comme étant égale à l'intersection de la surface projetée des deux fibres dans le plan de la liaison, ce dernier étant normal à la perpendiculaire commune aux deux fibres en contact (*e*<sub>iaβ</sub> = *e*<sub>α</sub> ∧ *e*<sub>β</sub>).

B. Calcul du tenseur de conductivité – Code éléments discrets

Connaissant la microstructure du VER, des expériences numériques peuvent être menées. Elles consistent
à (i) imposer un gradient de température macroscopique unitaire grad T[0] = e_k dans une des directions de
référence (= directions principales du VER), (ii) résoudre numériquement le système linéaire (15) pour
déterminer les fluctuations de température en chaque point de contact (méthode de Gauss), (iii) déterminer
le tenseur de conductivité macroscopique à l'aide de (16-17), (iv) reproduire ces étapes pour les deux autres
directions principales du VER. L'ensemble des ces calculs a été réalisé en modifiant un code de type
éléments discrets programmé sous MatLab, initialement développé pour simuler la rhéologie des
suspensions de fibres concentrées [25].

C. Influence des paramètres microstructuraux

Des VER constitués d'une boîte cubique de coté 25 mm et contenant des fibres droites de longueur 25 mm et de diamètre 0.5 mm ont été générés. Nous avons vérifié en amont que la taille des VER ainsi générés était suffisante pour obtenir des propriétés effectives indépendantes de cette taille. La conductivité thermique des fibres a été choisie égale à 1 Wm⁻¹K⁻¹. Nous avons étudié sur ces microstructures l'influence (i) de la fraction volumique de fibres, (ii) du coefficient de transfert thermique, (iii) de la longueur des fibres, (iv) de l'orientation du réseau de fibres, sur le tenseur de conductivité effectif. Remarquons que pour une orientation et un taux de fibres donnés, chaque calcul a été répété dix fois avec dix microstructures générées dans les mêmes conditions, pour tenir compte de la variabilité induite par le mode de génération des réseaux fibreux.

<figure><img src="image_5.png" /><figcaption>Figure 4 : Evolution de la conductivité principale macroscopique en fonction (a) de la fraction volumique de fibres *f* pour différentes qualités de contacts *h* (microstructures isotropes), (b) de *h* pour différentes orientations de microstructures (*f*= 0.15), (c) de la longueur des fibres *l* pour différentes orientations de microstructures (*f*= 0.15).</figcaption></figure>

1. Influence de la fraction volumique de fibres

Nous avons fait évoluer la fraction volumique de fibres de 0.025 à 0.3 pour des microstructures isotropes et
pour différents coefficients de transfert thermique. On constate sur la figure 4a que la conductivité thermique
augmente fortement avec la fraction volumique de fibres et ce, quel que soit le coefficient d'échange h. La
figure 4a met une nouvelle fois en évidence que plus la qualité des contacts est bonne (h élevé), plus la conductivité thermique est élevée. Les résultats ont été ajustés par des lois puissances ($A_{11} = k f^n$) représentées par des droites sur la figure 4a, de pentes respectives $n = 1$ ($h \to \infty$), $n = 1.9$ ($h = 1 \text{ Wm}^{-2}\text{K}^{-1}$) et $n = 2$ ($h = 10^{-3} \text{ Wm}^{-2}\text{K}^{-1}$). On constate que pour des contacts très isolants ($h = 10^{-3} \text{ Wm}^{-2}\text{K}^{-1}$), le meilleur exposant pour la loi puissance est de 2, alors que pour des contacts parfaits ($h \to \infty$), il est plutôt proche de 1.

## 2. Influence du coefficient de transfert thermique

L'influence du coefficient d'échange $h$ a été étudiée sur des microstructures isotropes ($A_{11} = A_{22} = A_{33} = 1/3$), des microstructures plutôt planes (⟨$A_{11}$, $A_{22}$⟩/$A_{33} = 50$) et isotropes dans le plan ($e_1$, $e_2$) ($A_{11} = A_{22}$), des microstructures planes et plutôt orientées dans le plan ($A_{11}/A_{22} = 3$). La fraction volumique de fibres est fixée à 0.15. Les résultats sont tracés figure 4b. Quelle que soit l'orientation du réseau fibreux et comme déjà noté sur les figures 1c et 2b, la conductivité thermique augmente d'abord rapidement, s'infléchit avant de se stabiliser pour des valeurs de coefficients de transfert thermique supérieurs à 1000 $\text{Wm}^{-2}\text{K}^{-1}$.

## 3. Influence de la longueur des fibres

Pour une même section droite de fibres (rayon $a = 0.25 \text{ mm}$), nous avons fait évoluer la longueur des fibres entre 5 mm et 25 mm en conservant une fraction volumique de fibres de 0.15 et des contacts parfaits ($h \to \infty$). Différentes orientations de réseaux fibreux ont été testées. Les résultats sont tracés figure 4c. On constate que la conductivité thermique effective augmente avec la longueur des fibres. Cette augmentation d'abord assez marquée se réduit sensiblement pour des longueurs supérieures à 15 mm, correspondant à des élancements de fibres supérieurs à 30. La même évolution est constatée quelle que soit l'orientation du réseau fibreux.

## 4. Influence de l'orientation du réseau

Deux cas de contacts ont été étudiés : des contacts très isolants ($h = 10^{-3} \text{ Wm}^{-2}\text{K}^{-1}$) et des contacts très conducteurs ($h \to \infty$). La fraction volumique de fibres est fixée à 0.15. Les résultats sont tracés figure 5. Nous avons étudié l'impact de l'anisotropie microstructurale sur celle du tenseur de conductivité thermique macroscopique dans le plan ($e_1$, $e_2$). On constate sur cette figure que lorsque l'anisotropie géométrique dans le plan augmente, l'anisotropie thermique croît. Elle semble égale à l'anisotropie géométrique lorsque les contacts sont parfaits (figure 5a) et ce, quelle que soit l'anisotropie géométrique plan/épaisseur. Lorsque les contacts sont très résistifs (figure 5b), cette remarque semble vérifiée pour une anisotropie géométrique plan/épaisseur de 1. Dans les autres cas, l'anisotropie thermique a tendance à être supérieure à l'anisotropie géométrique, ce phénomène s'amplifiant lorsque l'anisotropie géométrique dans le plan, mais aussi plan/épaisseur augmente.

<figure><img src="image_6.png" /><figcaption>Figure 5 : Evolution de l'anisotropie du tenseur de conductivité thermique dans le plan en fonction de l'anisotropie microstructurale dans le plan pour des contacts parfaits (a) et très résistifs (b) ($f = 0.15$).</figcaption></figure>

# V. COMPARAISON AVEC DES EXPERIENCES

Nous avons essayé de modéliser les mesures de conductivité thermique réalisées par Keith et al. [5] sur des composites à matrice Nylon 6.6 (masse volumique $1140 \text{ kg m}^{-3}$, conductivité thermique proche de $0.25 \text{ W m}^{-1} \text{ K}^{-1}$) chargés de fibres de carbone ThermalGraph DKD X fabriquées par Cytec (longueur initiale
200 µm, diamètre 10 µm, section circulaire, masse volumique 2200 kg m⁻³, conductivité thermique comprise entre 400 et 700 W m⁻¹ K⁻¹). Ce composite thermo-conducteur est produit sous forme de granules par un procédé d'extrusion. Un soin particulier a été apporté pour ne pas trop détériorer les fibres durant l'élaboration du composite. Les granules ainsi obtenus sont séchés puis fondus et injectés dans des moules sous forme de disques (diamètre 63.5 mm et d'épaisseur 3.18 mm, point d'injection sur le côté du disque). La conductivité thermique est caractérisée par une méthode de type source plane transitoire. Celle-ci permet à l'aide d'un disque chauffant (diamètre 6.4 mm) formé de spires concentriques, de caractériser simultanément la conductivité dans l'épaisseur et dans le plan principal d'un même disque. Le dispositif standard utilisé est un Mathis Instrument Hot Disk Thermal Constant Analyser.

Les informations concernant les microstructures des échantillons caractérisés sont assez partielles. Par exemple, au cours du processus d'élaboration des échantillons, la longueur moyenne des fibres chute de 200 µm à environ 100 µm [5], mais rien n'est dit sur la distribution des longueurs. De même, aucune information concernant l'éventuelle hétérogénéité de la microstructure dans le plan ou l'épaisseur des disques (distribution du taux et de l'orientation des réseaux fibreux) n'est donnée.

Dans la comparaison qui suit et en l'absence de ce type d'informations microstructurales, nous supposerons arbitrairement que taux, orientation et longueur de fibres sont homogènes dans les VERs et nous fixerons la longueur de fibres à 100 µm. Ainsi, pour modéliser ces essais, il est nécessaire de ne déterminer que deux paramètres du modèle discret : l'anisotropie du tenseur d'orientation et le coefficient de transfert thermique. En regard de la géométrie et du mode d'obtention des échantillons, nous avons considéré que les structures fibreuses dans ces derniers étaient plutôt isotropes dans le plan mais présentaient une anisotropie plan/épaisseur. D'après les résultats numériques obtenus dans le paragraphe précédent, celle-ci a été choisie égale à l'anisotropie moyenne du tenseur de conductivité thermique, i.e. ≈ 8 en moyenne pour les taux de fibres testés. En outre, les mesures de conductivité effectuées dans [5] recouvrent une très grande plage de fraction volumique de fibres, y compris dans des zones où les fibres risquent de ne pas être toutes connectées entre elles. Les modèles discrets développés dans ce travail ne sont pas adaptés à ces régimes de faibles concentrations de fibres. Pour cela, les expressions des conductivités macroscopiques établies précédemment ont été modifiées par l'ajout d'une contribution adaptée à la modélisation de la conductivité thermique dans les composites à fibres très conductrices dans les régimes diluées [24] :

$$ \Lambda^{\text{composite}} = \Lambda^{\text{eff}} + \Lambda_m \left( \mathbf{I} + f \left( \frac{l}{a} \right)^2 \frac{1}{6 \ln(2 l / a)} \mathbf{A} \right), \quad (20) $$

où I est le tenseur identité et où Λm représente la conductivité de la matrice.

<figure><img src="image_7.png" /><figcaption>Figure 6 : Evolution de la conductivité thermique (a) dans le plan et (b) transverse en fonction de la fraction de fibres. Influence des différentes contributions du modèle. Comparaison avec les résultats expérimentaux obtenus par Keith et al. [5].</figcaption></figure>

Les résultats obtenus à l'aide de ce modèle sont présentés sur les deux graphes de la figure 6. Le coefficient de transfert thermique h, supposé constant quelle que soit la fraction volumique de fibres a été évalué à 3.10⁵ W m⁻² K⁻¹. Les résultats sont donnés figure 6a pour la conductivité dans le plan et figure 6b pour la conductivité transverse. Pour chaque fraction volumique de fibres deux calculs sont réalisés, correspondant aux limites de variation de la conductivité thermique des fibres. Pour chaque point de calcul de Λeff, vingt microstructures sont générées et calculées pour tenir compte du mode de génération statistique des VER. Cette variation est représentée par les barres d'erreur. Les évolutions correspondant aux deux contributions apparaissant dans (20) ont également été portées dans les graphes de la figure 6.
Les résultats obtenus figure 6 montrent que la correction apportée par le modèle de suspension diluée joue
un rôle prépondérant dans les faibles taux de fibres. Son domaine d'influence est plus grand dans le sens
transverse (0 < *f* < 0.23) que dans le plan des disques (0 < *f* < 0.13). Au delà de ces zones, le modèle à
forte concentration développé ici prend le dessus et permet de bien rendre compte de la forte augmentation
de la conductivité thermique dans cette zone. Cette augmentation est plutôt quadratique et tend à montrer
que les contacts fibre-fibre pilotent la thermique (figure 4a et §IV.C.1). Dans l'ensemble, l'expression (20)
reproduit assez bien les tendances observées expérimentalement, sauf pour la conductivité effective dans le
plan dans les faibles fractions de fibres, où les résultats du modèle sous-estiment largement les données
expérimentales. Cette différence pourrait être due à la non uniformité de la microstructure dans l'épaisseur
des échantillons formant éventuellement des couches plus conductrices, ainsi qu'à l'éventuelle orientation
des fibres dans une des directions du plan. Une connaissance plus fine de la microstructure des composites
testés permettrait d'(in)valider cette hypothèse.

VI. CONCLUSION

En utilisant la méthode d'homogénéisation des structures périodiques par développements asymptotiques, une modélisation des phénomènes de conduction thermique transitoires et non linéaires dans les composites à renforts fibreux présentant une conductivité de la matrice négligeable devant celle des renforts ainsi que des résistances de contact entre fibres a été proposée. Sans aucun pré-requis sur la structure des équations de comportement ou de bilan à l'échelle macroscopique, les développements théoriques montrent que les continus équivalents sont des milieux monophasiques obéissant à une équation de la chaleur standard, et que les calculs de capacité thermique, de source de chaleur et de conductivité macroscopiques sont totalement découplés. Les deux premiers consistent à réaliser de « simples » moyennes volumiques sur le VER, le troisième consiste à résoudre sur le VER un problème de localisation linéaire et statique. Trois milieux continus équivalents sont possibles, selon la qualité des contacts fibres-fibres. Celui dédié aux mauvais contacts est d'emblée discret. Les deux autres peuvent être simplifiés et « discrétisés » en tenant compte de la géométrie particulière des fibres.

La qualité des prédictions des modèles discrets dépend fortement de l'élancement des fibres et de la qualité des contacts fibre-fibre : elle est excellente pour des fibres très élancées et pour de mauvais contacts, elle peut induire une erreur de l'ordre de 25% si l'élancement des fibres entre deux contacts est de l'ordre de 2.5 et si les contacts entre fibres sont excellents.

Les modèles discrets ont ensuite été implantés dans un code éléments discrets. En générant numériquement des réseaux de fibres enchevêtrées, cela a permis de souligner l'influence de paramètres microstructuraux tels que le taux, l'orientation et l'élancement des fibres, sur le tenseur de conductivité du milieu continu équivalent. Ainsi, nous avons pu montrer (i) que l'évolution de la conductivité thermique est linéaire en fonction de la fraction volumique de fibres lorsque les contacts entre fibres sont de bonne qualité et quadratique lorsque ils sont isolants, (ii) que pour des contacts conducteurs, une augmentation de la longueur des fibres génère une augmentation de la conductivité thermique effective et (iii) que l'anisotropie du tenseur de conductivité thermique est pratiquement égale à celle du tenseur d'orientation du second ordre.

Enfin, combinés à des modèles de conductivité dédiés aux régimes dilués, nous avons montré que les modèles développés dans ce travail, valables dans les régimes très concentrés, permettaient une assez bonne reproduction des résultats expérimentaux obtenus sur des composites composés de fibres de carbone baignées par une matrice de nylon 6.6.

References

[1] Y. Agari and T. Uno. *Thermal conductivity of polymer filled with carbon materials: Effect of conductive particle chains on thermal conductivity*. J. Applied Polymer Sc., 30:2225-2235,1985.

[2] T. Dupuis F. Danes, B. Garnier. *Predicting, measuring and tailoring the transverse thermal conductivity of composites from polymer matrix and metal filler*. Int. J. of Thermophysics, 24(3):771-784, 2003.

[3] Y. Agari and T. Uno. *Estimation on thermal conductivities of filled polymers*. J. Applied Polymer Sc., 2:5705-5712, 1986.

[4] T. Dupuis P. Lerendu T.-P. Nguyen F. Danès, B. Garnier. *Non-uniformity of the filler concentration and of the transverse thermal and electrical conductivities of filled polymer plates*. Composites Sc. and Tech., 5:945-951, 2005.

[5] J.M. Keith, C.D. Hingst, M.G. Miller, J.A. King, R.A. Hauser, *Measuring and predicting in-plane thermal conductivity of carbon-filled nylon 6,6 polymer composites*, Polymer Composites, 10:1001-1007, 2006
[6] E.H. Weber *Development and Modelling of thermally conductive polymer/Carbon composites*, PhD Dissertation, Michigan Technological University, 1999

[7] J C Maxwell. *A treatise on electricity and magnetism*. 1, 1891.

[8] M. Kaviany. *Principles of Heat Transfer in Porous Media*. Springer, New York, 1995.

[9] H.-J. Ott, *Thermal Conductivity of Composite Materials*. Plastic and Rubber Processing and Applications, Vol. 1, pp. 9-24, 1981.

[10] Eric S.G. Shaqfeh, *A nonlocal theory for heat transport in composites containing highly conducting fibrous inclusions* Phys. Fluids 31 (9), 1988

[11] M B Mackaplow, E S G Shaqfeh, and R L Schiek. *A numerical study of heat and mass transport in fibre suspensions*. Proc R Soc Lond A, 447:77-110, 1994.

[12] V. V. Calmidi and R. L. Mahajan. *The effective thermal conductivity of high porosity fibrous metal foams*. Journal of Heat Transfer, 121:466-471, 1999.

[13] L. E. Nielsen, *The Thermal and Electrical Conductivity of Two-Phase Systems*. I&EC Fundamentals, Vol. 13, No. 1, 17-20, 1974.

[14] L. E. Nielsen, and R. F. Landel, *Mechanical Properties of Polymers and Composites*. 2nd ed., Marcel Dekker, New York, 1994.

[15] D. J. Tanner-A. C. Cleland, J. K. Carson, S. J. Lovatt. *Thermal conductivity bounds for isotropic, porous media*. Int. J. of Heat and Mass Transfer, 48:2150-2158, 2005.

[16] X. Cheng, A. M. Sastry, and B. E. Layton. *Transport in stochastic fibrous networks*. J. Eng. Mater. Technol., 123:12-19, 2001.

[17] J.-P. Vassal, L. Orgéas, D. Favier, J.-L. Auriault, S. Le Corre *Micro-Macro models for transient diffusion in granular and entangled fibrous media*. Theoretical aspects, submitted in Phys. Rev. E, 2006

[18] E. Sanchez-Palencia, *Non-Homogeneous Media and Vibration Theory*, Lecture Notes in Physics Springer-Verlag, Berlin, Germany, 1980

[19] A. Bensoussan, J.-L. Lions, and G. Papanicolaou. *Asymptotic Analysis for periodic structures*. North Holland, Amsterdam, 1978.

[20] J.-L. Auriault. *Heterogeneous medium. is an equivalent description possible?* Int. J. Engng. Sci., 29:785-795, 1991.

[21] Advani, S. G. & Tucker, C. L., *The use of tensors to describe and predict fiber orientation in short fiber composites*, J. Rheol. 3(8), 751-784, 1987.

[22] Doi, M., and S. F. Edwards, *Dynamics of rod-like macromolecules in concentrated solution*, J. Chem. Soc. Faraday Trans. 2(74), 560-570, 1978.

[23] Ranganathan, S., and S. G. Advani, *Fiber-fiber interactions in homogeneous flows of non-dilute suspensions*, J. Rheol. 35, 1499-1522, 1991.

[24] Toll, S., *Note: On the tube model for fiber suspensions*, J. Rheol. 37, 123-125, 1993.

[25] A. Rocha, A. Acrivos, Q. J. Mech. Math. 26 (1973) 441, 1974.

[26] Le Corre, S., Dumont, P., Orgéas, L. & Favier, D., *Rheology of highly concentrated fiber suspensions*, J. Rheol. 49, 1029-1058, 2005.