# Prédiction par transfert inverse d'un champ de conductance thermique de contact dans un mur de réacteur métallurgique

CLÉMENT ROUSSEAUª, MARCEL LACROIX ET MARTIN DÉSILETS

Université de Sherbrooke 2500, Boul. de l'Université, Sherbrooke (Québec), J1K 2R1, Canada

Reçu le 14 avril 2011, accepté le 13 mars 2014

**Résumé** – Cet article porte sur une méthode d'estimation, novatrice et non intrusive, de la conductance thermique de contact dans un mur de réacteur métallurgique. Dans le cas d'un réacteur métallurgique, il est essentiel de déterminer la conductance thermique de contact dans l'assemblage des parois. Cela permet d'identifier et de corriger les défauts de contact avant la mise en fonctionnement du réacteur métallurgique. Cette étude recourt à la méthode inverse du gradient conjugué avec un problème adjoint, pour réaliser cette estimation. La méthode a été validée à l'aide de tests numériques représentant un diagnostic du contact thermique entre une paroi réfractaire et une paroi d'acier du réacteur métallurgique avec différentes conductances thermiques de contact. Dans un second temps, des tests supplémentaires ont été réalisés pour observer l'évolution de l'erreur d'estimation en fonction de deux nombres sans dimension, le rapport de conductivité thermique et le nombre de Biot. Cela a permis de conclure qu'il faut avoir simultanément un rapport de conductivité thermique supérieur à 1 et un nombre de Biot supérieur à 0,005 pour obtenir une estimation précise dans les cas étudiés.

**Mots clés** : Méthode inverse / conductance thermique de contact / gradient conjugué avec problème adjoint / réacteur métallurgique

**Abstract** – Identification by inverse method of a thermal contact conductance field in a metallurgical reactor. This work provides an inverse estimation method for the prediction of the thermal contact conductance field between two materials layers inside the wall of a metallurgical reactor. In this case, it is very important to determine the thermal contact conductance in the wall because if it's importance on the heat balance. This technique allows identifying and correcting contact defaults before the reactor is started. This study uses an inverse method based on the conjugate gradient method with adjoint problem, to achieve this identification. Firstly, this method is validated for three test cases with different thermal contact conductances at the interface between two materials usually found on the side wall of a metallurgical reactor: a refractory layer in contact with a steel plate. Secondly, additional tests are carried out to study the evolution of the error estimation as a function of two adimensional parameters, the Biot number and the thermal conductivity ratio $k^*$. It is shown that in order to have a good estimation, the Biot number must be larger than 0.005 while $k^*$ needs to be larger to 1. This method of estimation of the thermal contact conductance field is innovative and non intrusive.

**Key words:** Inverse method / thermal contact conductance / conjugate gradient method with adjoint problem / metallurgical reactor

## 1 Introduction

Les murs de fours industriels à haute température sont faits de couches de matériaux superposées (Fig. 1). Or, lors de l'assemblage, il est rare que les différentes parois des murs s'épousent parfaitement. Des interstices peuvent subsister entre les parois. Lorsque le four est exploité à

haute température, ces interstices agissent comme un isolant et créent des points chauds. En ces endroits, la chaaleur est mal évacuée par les murs du four. À long terme, ces défauts d'assemblage peuvent influer sur le procédé industriel et même menacer l'intégrité de l'installation.

Il est toutefois très difficile de localiser ces défauts de contact. En effet, il n'y a aucune méthode de measure non intrusive qui permette de connaître la conductance thermique de contact dans un assemblage sans

<footnote>ª Auteur pour correspondance : rousseau.cl@gmail.com</footnote>
Nomenclature

<table><thead><tr><th>Symbole</th><th>Définition</th></tr></thead><tbody><tr><td colspan="2">Lettres romaines</td></tr><tr><td>Biot</td><td>Nombre de Biot</td></tr><tr><td>De</td><td>Direction de descente</td></tr><tr><td>E</td><td>Erreur d'estimation générale (%)</td></tr><tr><td>Ep</td><td>Erreur d'estimation ponctuelle (%)</td></tr><tr><td>h</td><td>Coefficient de convection extérieure (W/m².K)</td></tr><tr><td>H(y)</td><td>Conductance thermique de contact estimée (W/m².K)</td></tr><tr><td>H<sub>réel</sub></td><td>Conductance thermique de contact recherchée (W/m².K)</td></tr><tr><td>Jq</td><td>Gradient trouvé à l'aide du problème adjoint</td></tr><tr><td>k*</td><td>Rapport des conductivités thermiques</td></tr><tr><td>k</td><td>Conductivité thermique (W/m.K)</td></tr><tr><td>L</td><td>Largeur de l'assemblage (m)</td></tr><tr><td>Ly</td><td>Hauteur de l'assemblage (m)</td></tr><tr><td>M</td><td>Résultat numérique dans la méthode inverse</td></tr><tr><td>n<sub>y</sub></td><td>Nombre de points de mesure</td></tr><tr><td>P</td><td>Paramètre dans la méthode inverse</td></tr><tr><td>R</td><td>Résistance thermique K.m²/W</td></tr><tr><td>S</td><td>Écart quadratique</td></tr><tr><td>T<sub>c</sub></td><td>Température imposée en x = 0 (K)</td></tr><tr><td>Tm(y)</td><td>Température mesurée en x = L (K)</td></tr><tr><td>T<sub>ext</sub></td><td>Température extérieure (K)</td></tr><tr><td>T(i)</td><td>Température au point de mesure i (K)</td></tr><tr><td>U</td><td>Source dans la méthode inverse</td></tr><tr><td>X<sub>0</sub></td><td>État initial dans la méthode inverse</td></tr><tr><td>Y</td><td>Résultat expérimental dans la méthode inverse</td></tr><tr><td colspan="2">Lettres grecques</td></tr><tr><td>β</td><td>Coefficient utilisé pour calculer De</td></tr><tr><td>ΔT(i)</td><td>Température donnée par le problème en variations au point de mesure i (K)</td></tr><tr><td>ε</td><td>Critère d'arrêt</td></tr><tr><td>ρ</td><td>Paramètre de descente</td></tr><tr><td>θ</td><td>Température adimensionnée</td></tr><tr><td>θ<sub>c</sub></td><td>Température adimensionnée imposée en x* = 0</td></tr><tr><td>θ<sub>ext</sub></td><td>Température adimensionnée extérieure</td></tr><tr><td colspan="2">Exposant</td></tr><tr><td>*</td><td>Valeur adimensionnée</td></tr><tr><td>k</td><td>Itération</td></tr><tr><td colspan="2">Indices</td></tr><tr><td>1</td><td>Zone 1</td></tr><tr><td>2</td><td>Zone 2</td></tr><tr><td>1c</td><td>Conduction dans la zone 1</td></tr><tr><td>2c</td><td>Conduction dans la zone 2</td></tr><tr><td>c</td><td>Zone de contact</td></tr><tr><td>ext</td><td>Zone extérieure</td></tr></tbody></table>

endommager celui-ci [1]. Dans cet article, une méthode novatrice est développée reposant sur les méthodes inverses pour réaliser l'estimation de la conductance thermique de contact. Cette méthode inverse utilise les températures mesurées sur la surface extérieure $T_m(y)$ en $x = L$, voir figure 1, pour estimer $H(y)$, la conductance thermique de contact. Le principe s'appuie sur la

mesure de l'effet, la distribution de température sur la surface extérieure, pour estimer la cause, les défauts de contact. La méthode utilisée est le gradient conjugué avec un problème adjoint. Plusieurs techniques d'estimation peuvent être utilisées parmi lesquelles, l'estimation de paramètre et l'estimation de fonction. Une estimation de paramètres est utilisée lorsque le nombre de paramètres
<figure><img src="image_1.png" /><figcaption>Fig. 1. Montage d'un mur de réacteur métallurgique.</figcaption></figure>

recherchés est faible. Une estimation de fonction permet plutôt d'identifier une fonction inconnue à l'aide de nombreux points.

Plusieurs auteurs utilisent les méthodes inverses pour estimer les conductances thermiques de contact [2-4], mais en considérant cette conductance comme uniforme sur toute la surface de contact. Gill et al. [5] ont réussi à développer une méthode de mesure de la conductance de contact variant dans l'espace, mais en utilisant une prise de mesure intrusive, c'est à dire en relevant des températures internes. Feuillet [6], quant à lui, a réussi à mettre au point une méthode permettant l'estimation d'un champ de résistance thermique de contact entre deux plaques à l'aide d'une technique de modélisation des résistances en différences finies. L'estimation réalisée dans cette étude est une estimation de la résistance thermique de contact selon certaines zones définies. L'auteur a recherché une conductance thermique de contact constante par zone. L'ensemble de ces zones forme une grille sur la surface. Il s'agit d'une estimation de paramètres moins précise qu'une estimation de fonction.

L'objectif de ce travail est de mettre au point une méthode d'estimation de fonction adaptable à différentes géométries et conditions aux limites, avec une prise de mesure non intrusive. Les éléments suivants sont pris en compte dans cette étude :

- la conductance thermique de contact n'étant pas uniforme sur la surface, il ne sera pas possible de se baser sur les études [2-4];

- l'estimation se fera pendant le diagnostic d'un réacteur métallurgique. La prise de mesure se fera donc d'une manière non intrusive, à la surface extérieure ($T_m(y)$ à $x = L$), pour ne pas endommager

le four. L'origine du repère est le coin inférieur gauche du mur du réacteur métallurgique;

- l'estimation de fonction est utilisée pour permettre une meilleure précision et une meilleure adaptabilité du procédé.

Avec une méthode de mesure de la conductance thermique de contact, l'industriel peut réaliser un diagnostic pendant la conception des fours pour vérifier la qualité du contact thermique, et donc corriger les défauts avant la mise en fonctionnement du four.

Dans ce qui suit, la méthode directe permettant de modéliser le transfert thermique dans l'assemblage sera d'abord présentée. Puis le principe des méthodes inverses sera décrit ainsi que les équations du problème en variation et celles du problème adjoint. Ensuite, les résultats de tests numériques d'estimation seront présentés pour valider la méthodologie. Enfin, une étude sur l'influence de paramètres adimensionnés permettra d'observer la plage de validité de la méthode.

# 2 Méthode directe

## 2.1 Hypothèses

Dans cette étude, certaines hypothèses générales ont été prises en compte :

- régime stationnaire;

- propriétés thermiques constantes;

- deux dimensions;

- la chaleur est transmise par conduction à travers la paroi réfractaire et la paroi d'acier, $T_c$ est constante en $x = 0$.
Cette étude est réalisée en considérant le régime stationnaire, un état qui sera atteint si la phase de chauffe du four dure assez longtemps. Cette hypothèse permet d'augmenter la sensibilité de la prise de mesure par rapport à la zone de contact. En effet, l'impact des points chauds sur le profil de température $T_m(y)$ est plus important lorsque le régime stationnaire est établi.

La deuxième hypothèse peut être expliquée par le fait que la plage de variation de la température demeure faible (environ 40 K dans le cas présenté ci-dessous), donc la variations de conductivité thermique (liées au change-ment de température) de la paroi réfractaire et du caisson d'acier sont négligeables.

Le contact est considéré comme uniforme selon la troisième dimension. En effet, les gradients thermiques dans les directions $x$ et $y$ sont habituellement nettement supérieurs à ceux dans la direction $z$. Les tests sont donc réalisés en deux dimensions. Une étude en 3 dimensions est très lourde en besoin de puissance de calcul, une étude en 2 dimensions est la première étape avant une étude plus complète.

Enfin, pour avoir un transfert thermique dans l'assemblage, une température est imposée en $x = 0$, $T_c$ (voir Fig. 1). Cette condition est nécessaire pour réaliser une estimation. Dans un test réel, il sera nécessaire d'imposer une température constante sur la surface intérieure du four pour réaliser le diagnostic.

## 2.2 Équations du problème

Le problème est défini avec les équations de diffusion de la chaleur suivantes [7, 8]. $T_1(x, y)$ est la distribution de la température dans la zone 1 et $T_2(x, y)$ est la distribution de la température dans la zone 2.

$$0 = \left( \frac{\partial^2 T_1}{\partial x^2} + \frac{\partial^2 T_1}{\partial y^2} \right) \qquad (1)$$

$$0 = \left( \frac{\partial^2 T_2}{\partial x^2} + \frac{\partial^2 T_2}{\partial y^2} \right) \qquad (2)$$

L'équation (3) permet de définir le transfert thermique au niveau du contact.

$$\begin{aligned} -k_1 \left. \frac{\partial T_1}{\partial x} \right|_{x=L_1} &= -k_2 \left. \frac{\partial T_2}{\partial x} \right|_{x=L_1} \\ &= H(y) [T_1(L_1, y) - T_2(L_1, y)] \end{aligned} \qquad (3)$$

Enfin, les conditions aux limites du système sont définies par les équations (4) à (6) :

$$T_1(0, y) = T_c \qquad (4)$$

$$-k_2 \left. \frac{\partial T_2}{\partial x} \right|_{x=L} = h [T_2(L, y) - T_{ext}] \qquad (5)$$

$$\left. \frac{\partial T_1}{\partial y} \right|_{y=0} = \left. \frac{\partial T_1}{\partial y} \right|_{y=L_y} = \left. \frac{\partial T_2}{\partial y} \right|_{y=0} = \left. \frac{\partial T_2}{\partial y} \right|_{y=L_y} = 0 \qquad (6)$$

$h$ est le coefficient de convection à l'extérieur et $T_{ext}$ est la température extérieure.

Pour permettre de traiter le cas d'une manière générale, les équations (1) à (6) ont été adimensionnées avec :

$$k^* = \frac{k_1}{k_2} \qquad (7)$$

$$x^* = \frac{x}{L_2} \qquad (8)$$

$$y^* = \frac{y}{L_2} \qquad (9)$$

$$H^* = \frac{H \cdot L_2}{k_2} \qquad (10)$$

$$Biot = \frac{h \cdot L_2}{k_2} \qquad (11)$$

$$\theta = \frac{T - T_{ext}}{T_c - T_{ext}} \qquad (12)$$

Le nombre de Biot, exprime un rapport entre les résistances au transfert thermique à l'intérieur et à la surface d'un corps. Le rapport $H^*$ exprime le rapport de la résistance au transfert thermique à l'intérieur de la zone 2 sur la résistance thermique de contact.

$$0 = \left( \frac{\partial^2 \theta_1}{\partial x^{*2}} + \frac{\partial^2 \theta_1}{\partial y^{*2}} \right) \qquad (13)$$

$$0 = \left( \frac{\partial^2 \theta_2}{\partial x^{*2}} + \frac{\partial^2 \theta_2}{\partial y^{*2}} \right) \qquad (14)$$

$$\begin{aligned} -k^* \left. \frac{\partial \theta_1}{\partial x^*} \right|_{x^*=L_1^*} &= - \left. \frac{\partial \theta_2}{\partial x^*} \right|_{x^*=L_1^*} \\ &= H^*(y^*) [\theta_1(L_1^*, y^*) - \theta_2(L_1^*, y^*)] \end{aligned} \qquad (15)$$

$$\theta_1(0, y^*) = 1 \qquad (16)$$

$$- \left. \frac{\partial \theta_2}{\partial x^*} \right|_{x^*=L^*} = Biot [\theta_2(L^*, y^*)] \qquad (17)$$

$$\begin{aligned} \left. \frac{\partial \theta_1}{\partial y^*} \right|_{y^*=0} &= \left. \frac{\partial \theta_1}{\partial y^*} \right|_{y^*=L y^*} = \left. \frac{\partial \theta_2}{\partial y^*} \right|_{y^*=0} \\ &= \left. \frac{\partial \theta_2}{\partial y^*} \right|_{y^*=L y^*} = 0 \end{aligned} \qquad (18)$$

Le milieu 2 a été pris comme référence pour l'adimensionalnement. En effet, ce sont les propriétés thermophysiiques du milieu 2 qui influencent le plus la sensibilité du champ de température mesuré en $x^* = L^*$ par rapport à la conductance thermique de contact.

Pour valider ce modèle numérique, une comparaison a été réalisée, au tableau 1, entre les températures adimensionalnées en $x^* = L^*$ de la solution exacte, équation (29), et celles du modèle numérique. Dans cette comparaison, la conductance thermique de contact est considéré constante sur $y$. La solution exacte, l'équation (29), a été trouvée à l'aide d'une approche basée sur le calcul des résistances thermiques [9] de la manière suivante, voir figure 2.
<table><caption>**Tableau 1.** Comparaison solution exacte et numérique.</caption><thead><tr><th>H<sup>*</sup></th><th>&theta;(x<sup>*</sup> = L<sup>*</sup>) Solution exacte</th><th>&theta;(x<sup>*</sup> = L<sup>*</sup>) Modèle numérique</th></tr></thead><tbody><tr><td>0,25</td><td>0,7619</td><td>0,7619</td></tr><tr><td>2,50</td><td>0,7890</td><td>0,7890</td></tr><tr><td>25,00</td><td>0,7918</td><td>0,7918</td></tr><tr><td>250,00</td><td>0,7920</td><td>0,7920</td></tr></tbody></table>

<figure><img src="image_2.png" /><figcaption>**Fig. 2.** Circuit thermique.</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>**Fig. 3.** Schéma du problème direct.</figcaption></figure>

D'après la loi de conservation des flux thermiques :

$$R_{1c} = \frac{L_1^*}{k^* L_y^*} \qquad (19)$$

$$R_{ct} = \frac{1}{H^* L_y^*} \qquad (20)$$

$$R_{2c} = \frac{1}{L_y^*} \qquad (21)$$

$$R_{ext} = \frac{1}{Biot L_y^*} \qquad (22)$$

$$R_{totale} = R_{1c} + R_{ct} + R_{2c} + R_{ext} \qquad (23)$$

$$\phi = \frac{\theta_c - \theta_{ext}}{R_{totale}} \qquad (24)$$

$$\phi = \frac{\theta_c - \theta_{x^*=L^*}}{R_{1c} + R_{ct} + R_{2c}} \qquad (25)$$

$$\theta_c = 1 \qquad (26)$$

$$\theta_{ext} = 0 \qquad (27)$$

$$\frac{1}{\frac{L_1^*}{k^* L_y^*} + \frac{1}{Biot L_y^*} + \frac{1}{H^* L_y^*} + \frac{1}{L_y^*}} = \frac{1 - \theta_{x^*=L^*}}{\frac{L_1^*}{k^* L_y^*} + \frac{1}{H^* L_y^*} + \frac{1}{L_y^*}} \qquad (28)$$

Cela permet d'en déduire l'équation (29) :

$$\theta_{x^*=L^*} = \theta_{2(x^*=L_2^*)} = 1 - \frac{\frac{L_1^*}{k^* L_y^*} + \frac{1}{L_y^* H^*} + \frac{1}{L_y^*}}{\frac{L_1^*}{k^* L_y^*} + \frac{1}{Biot L_y^*} + \frac{1}{H^* L_y^*} + \frac{1}{L_y^*}} \qquad (29)$$

Le modèle numérique a été résolu avec la méthode des éléments finis avec le logiciel commercial COMSOL

MULTIPHYSIC 3.5a. Le maillage utilisé dans ce modèle, est composé d'éléments triangulaires avec comme contrainte de taille maximale 0.1 m dans la zone 1 et 0.001 m dans la zone 2 (pour améliorer la précision des résultats), la taille maximale étant fixée automatiquement par le logiciel. Le modèle est composé de 21 396 éléments quadratiques.

Les températures adimensionnées du modèle numérique et de la solution exacte sont identiques pour toutes les valeurs de $H^*$. Le modèle numérique est donc validé.

## 3 Méthode inverse

### 3.1 Principe

Les méthodes directes [10, 11] sont les méthodes utilisées classiquement pour la modélisation numérique. Elles consistent à construire un modèle, voir figure 3, avec l'état initial, les paramètres et l'entrée du système pour en déduire les résultats du modèle, ici $M(P)$, l'inconnue du système. Les effets sont déduits des causes.

Afin de référer aux variables de la présente étude :

- $U : T_c$ température imposée,

- $P : H(y)$ conductance thermique de contact,

- $X_0 : \text{température initiale du système,}$

- $M(P) : T(i)$ prises de températures.

Dans cette étude, la température initiale n'est pas nécessaire pour résoudre le problème étant donné que le régime est stationnaire.
<figure><img src="image_4.png" /><figcaption>Fig. 4. Schéma du problème inverse.</figcaption></figure>

Dans les méthodes inverses [12, 13], les effets sont connus et les causes sont recherchées.

Le système est donc de la forme (Fig. 4) :

Les résultats $M(P)$ sont construits comme dans le cas des méthodes directes avec les $U, P$ et $X_0$. La différence est qu'ici les résultats du modèle $M(P)$ sont recherchés pour être égaux aux mesures $Y$. C'est donc $U, X_0$ et $P$ qui sont les inconnues dans les méthodes inverses. L'inconnue du système est donc recherchée pour faire concorder $M(P)$ à $Y$.

Il existe différentes méthodes inverses, les plus connues et utilisées étant les méthodes de descente. Ces dernières sont basées sur la minimisation de l'écart quadratique $S$, équation (30), qui est la différence entre les prédicitions du modèle $M(P)$ et les mesures expérimentales $Y$.

$$S^k = \frac{1}{2}[M(P) - Y]^2 \quad (30)$$

Dans le cas présent, les prédicitions $M(P)$ sont les prises de températures $T(i)$, i variant de 1 à $n_y$ nombre de mesures, avec la conductance thermique de contact $H$ comme paramètre $P$ recherché.

Le problème inverse devient donc un problème de minimisation de cette fonction, équation (30), appelée fonctionnelle ou fonction objectif. La minimisation est arrêtée quand l'écart quadratique atteint une valeur seuil appelée critère d'arrêt, $\epsilon$, voir l'équation (31).

$$S^k \le \epsilon \quad (31)$$

Dans cette étude, le modèle mathématique direct sera utilisé pour créer un champ de température à l'intérieur du domaine de calcul. Ces prédicitions vont bien sûr varier en fonction de la conductance de contact de base, $H_{réel}$, qui aura été choisie pour représenter la paroi industrielle. Ces prédicitions du modèle direct seront considérées comme des mesures expérimentales notées $Y$.

Dans cet article la méthode de descente utilisée est basée sur la méthode du gradient conjugué avec un problème adjoint, GCPA [13, 14].

## 3.2 Méthode GCPA

La méthode du GCPA est une méthode d'estimation de fonction, permettant la minimisation de la fonctionnelle, l'équation (30). Ici la fonction recherchée est la conductance thermique de contact $H^*(y)$. Cette méthode est basée sur l'algorithme suivant, voir l'équation (32) :

$$p^{k+1} = p^k + \rho^k De^k \quad (32)$$

- $p$: paramètre calculé à l'itération $k$;

- $De$: la direction de descente;

- $\rho$: le paramètre de descente.

Dans le cas d'une estimation de fonction, les paramètres recherchés sont les points de la fonction recherchée. Dans cette étude, les paramètres $p$ sont donc les points de la fonction $H^*(y)$.

L'algorithme de la méthode du GCPA se rapproche fortement de celui de la méthode du gradient conjugué, la différence étant le calcul du paramètre de descente et de la direction de descente [14].

La direction de descente $De$ de la méthode correspond au sens, c'est à dire à la diminution ou à l'augmentation, associé à la modification des paramètres recherchés. Elle dépend du gradient, $J_q$, calculé par le problème adjoint (voir Sect. 3.3) et des directions de descente précédentes, l'équation (33) :

$$De^k = -J_q^k + \beta De^{k-1} \quad (33)$$

Cette variable est initialisée à 0 dans le cas de la première itération ($k=1$).

Dans cette étude, la formule de Fletcher Reeve [14] a été utilisée pour calculer $\beta$, l'équation (34) :

$$\beta = \frac{\int_{0}^{tf} J_q(k)^2 dt}{\int_{0}^{tf} J_q(k-1)^2 dt} \quad (34)$$

Le paramètre de descente est un paramètre qui dépend du problème en variation (voir Sect. 3.3) et de la différence entre les valeurs mesurées et les valeurs expérimentales.
La formule utilisée pour le paramètre de descente est donnée par l'équation (35) [14] :

$$ \rho = \frac{\sum_{i=1}^{n_y} \int_{0}^{tf} [T(i) - Y] \Delta T(i) dt}{\sum_{i=1}^{n_y} \int_{0}^{tf} \Delta T(i)^2 dt} \quad (35) $$

Avec

- $\Delta T(i)$ : température donnée par le problème en variation au point de mesure $i$,

- $T(i)$ : température donnée par le problème direct au point de mesure $i$,

- $Y$ : température mesurée au point $i$,

- $n_y$ : nombre de points de mesure.

### 3.3 Équations de la méthode inverse

Un problème en variation est nécessaire pour la méthode du GCPA. Le problème en variation correspond à la réponse du système à une variation de la fonction recherchée. Les inconnues de ce problème sont donc les variations de température $\Delta\theta$. Ce problème permet de calculer le paramètre de descente de la méthode.

Voici les équations (13) à (18) modifiées pour le problème en variation adimensionné. Pour plus de détails au sujet de la construction du problème en variation, le lecteur peut consulter les références [14, 15].

$$ 0 = \left( \frac{\partial^2 \Delta\theta_1}{\partial x^{*2}} + \frac{\partial^2 \Delta\theta_1}{\partial y^{*2}} \right) \quad (36) $$

$$ 0 = \left( \frac{\partial^2 \Delta\theta_2}{\partial x^{*2}} + \frac{\partial^2 \Delta\theta_2}{\partial y^{*2}} \right) \quad (37) $$

$$ \begin{aligned} -k^* \frac{\partial \Delta \theta_1}{\partial x^*} \Big|_{x^*=L_1^*} &= - \frac{\partial \Delta \theta_2}{\partial x^*} \Big|_{x^*=L_1^*} \\ -k^* \frac{\partial \Delta \theta_1}{\partial x^*} \Big|_{x^*=L_1^*} &= H^*(y^*) [\Delta\theta_1(L_1^*, y^*) - \Delta\theta_2(L_1^*, y^*)] \\ &\quad + \Delta H^*(y^*) [\theta_1(L_1^*, y^*) - \theta_2(L_1^*, y^*)] \end{aligned} \quad (38) $$

$$ \Delta\theta_1(0, y^*) = 0 \quad (39) $$

$$ \left. \frac{\partial \Delta \theta_2}{\partial x^*} \right|_{x^*=L^*} = \operatorname{Biot} \Delta \theta_2(L^*, y^*) \quad (40) $$

$$ \begin{aligned} \left. \frac{\partial \Delta \theta_1}{\partial y^*} \right|_{y^*=0} &= \left. \frac{\partial \Delta \theta_1}{\partial y^*} \right|_{y^*=L_y^*} = \left. \frac{\partial \Delta \theta_2}{\partial y^*} \right|_{y^*=0} \\ &= \left. \frac{\partial \Delta \theta_2}{\partial y^*} \right|_{y^*=L_y^*} = 0. \end{aligned} \quad (41) $$

Pour minimiser la fonctionnelle, équation (30), il est nécessaire de faire intervenir des multiplicateurs de Lagrange. En effet, la température adimensionnée $\theta$ qui intervient dans l'équation (30) dépend de la méthode directe, il y a donc une contrainte à respecter. Un problème adjoint est donc ajouté au problème en variation pour obtenir le gradient $J_q$, utilisé pour le calcul de la direction

de descente. Pour plus de détails au sujet du problème adjoint, le lecteur peut consulter les références [14, 16].

$$ \begin{aligned} 0 = k^* \left( \frac{\partial^2 P_1}{\partial x^2} + \frac{\partial^2 P_1}{\partial y^{*2}} \right) &+ \sum_{i=1}^{ny} \left[ (Y_i^* - \theta_1(x^*, y^*)) \right. \\ & \left. \times \delta(x^* - x_i^*)(y^* - y_i^*) \right] \end{aligned} \quad (42) $$

$$ \begin{aligned} 0 = & \left( \frac{\partial^2 P_2}{\partial x^{*2}} + \frac{\partial^2 P_2}{\partial y^{*2}} \right) + \sum_{i=1}^{ny} \left[ (Y_i^* - \theta_2(x^*, y^*)) \right. \\ & \left. \times \delta(x^* - x_i^*)(y^* - y_i^*) \right] \end{aligned} \quad (43) $$

$$ -k^* \left. \frac{\partial P_1}{\partial x^*} \right|_{x^*=L_1} = - \left. \frac{\partial P_2}{\partial x^*} \right|_{x=L_1^*} \quad (44) $$

$$ -k^* \left. \frac{\partial P_1}{\partial x^*} \right|_{x=L_1} = H^*(y^*) [P_1(L_1^*, y^*) - P_2(L_1^*, y^*)] \quad (45) $$

$$ P_1(0, y^*) = \left. \frac{\partial P_1}{\partial x^*} \right|_{x^*=0} = 0 \quad (46) $$

$$ k^* \left. \frac{\partial P_2}{\partial x^*} \right|_{x^*=L^*} = \operatorname{Biot} \times P_2(L^*, y^*) \quad (47) $$

$$ \begin{aligned} \left. \frac{\partial P_1}{\partial y^*} \right|_{y^*=0} &= \left. \frac{\partial P_1}{\partial y^*} \right|_{y^*=L y^*} = \left. \frac{\partial P_2}{\partial y^*} \right|_{y^*=0} \\ &= \left. \frac{\partial P_2}{\partial y^*} \right|_{y^*=L y^*} = 0 \end{aligned} \quad (48) $$

Le gradient $J_q$ est déduit de ce problème, équation (49) :

$$ J_q = (\theta_1|_{x^*=L_1^*} - \theta_2|_{x^*=L_1^*})(P_1|_{x^*=L_1^*} - P_2|_{x^*=L_1^*}) \quad (49) $$

Celui-ci sera utilisé dans l'équation (33).

La fonction de Dirac $\delta$ dans les équations (42) et (43) permet d'imposer un flux thermique en $x^* = x_i^*$ et $y^* = y_i^*$, c'est à dire au niveau des prises de mesure. Ce flux thermique est égal à la différence entre les températures du modèle et les mesures expérimentales.

### 3.4 Résolution inverse

Une méthode d'éléments finis construite avec le logiciel Comsol Multiphysic 3.5a a été utilisée pour la résolution des équations des problèmes directs, en variation et adjoint.

L'algorithme du GCPA, qui a été programmé est présenté à la figure 5.

- L'étape 1 consiste à initialiser les variables pour la création des modèles, propriétés géométriques et thermiques, ainsi que les valeurs utilisées dans la méthode inverse, comme le critère d'arrêt et la valeur initiale des points de la fonction $H^*$.

- L'étape 2 construit les valeurs $Y^*$, considérées comme des mesures expérimentales, à l'aide de la conductance thermique de contact $H_{r\acute{e}el}^*$ imposée par l'utilisateur.

- À l'étape 3, la résolution du modèle direct est réalisée avec la valeur de $H^*$ à l'itération $k$. Dans le cas de
<figure><img src="image_5.png" /><figcaption>Fig. 5. Méthodologie de la méthode d'estimation.</figcaption></figure>

la première itération, cette valeur est la valeur initiale choisie à l'étape 1. Les valeurs $\theta^{*k}$ sont construites à cette étape.

- L'étape 4 illustre la résolution du modèle adjoint qui permet la construction du gradient $J_q$ utilisé dans le calcul de la direction de descente, l'équation (33).

- Le problème en variation est résolu pendant l'étape 5. Elle permet de construire les températures en variation utilisées pour le calcul du paramètre de descente.

- La nouvelle fonction $H^*$ est calculée dans l'étape 6 ainsi que le nouveau champ de température $\theta^{*k+1}$.

- L'étape 7 permet de calculer la fonctionnelle, l'équation (30). Si cette fonctionnelle est inférieure au critère d'arrêt de l'équation (31), alors l'estimation est finie. Sinon, la méthode retourne à l'étape 4 avec $k = k + 1$.

Cet algorithme a été programmé à l'aide du logiciel Matlab R2008a. L'avantagee d'utiliser Matlab R2008a est la possibilité de faire appel aux modélisations Comsol Multiphysic 3.5a et de pouvoir récupérer les résultats facilement.

# 4 Résultats

## 4.1 Validation de la méthode

La méthode sera éprouvée à l'aide de trois cas, notés A, B et C, possédant la même géométrie et les mêmes

<table><caption>Tableau 2. Paramètres géométriques et propriétés thermophysiques des cas A, B et C.</caption><tr><td>$k_1$</td><td>30 W/m.K</td></tr><tr><td>$k_2$</td><td>40 W/m.K</td></tr><tr><td>$L_1$</td><td>0,15 m</td></tr><tr><td>$L_2$</td><td>0,01 m</td></tr><tr><td>$L_y$</td><td>1 m</td></tr><tr><td>$T_c$</td><td>450 K</td></tr><tr><td>$T_{ext}$</td><td>293 K</td></tr><tr><td>$h$</td><td>50 W/m².K</td></tr></table>

<table><caption>Tableau 3. Paramètres géométriques et propriétés thermophysiques adimensionnés des cas A, B et C</caption><tr><td>$k_1^*$</td><td>0,75</td></tr><tr><td>$k_2^*$</td><td>1</td></tr><tr><td>$L_y^*$</td><td>100</td></tr><tr><td>$L_1^*$</td><td>15</td></tr><tr><td>$L_2^*$</td><td>1</td></tr><tr><td>Biot</td><td>0,025</td></tr></table>

paramètres présentés dans les Tableaux 2 et 3, mais avec une conductance thermique de contact $H_{réel}^*$ différente.

La configuration proposée ci-dessus, représente le cas d'un diagnostic d'un mur de réacteur métallurgique, avec dans la zone 1 une dalle de carbone et dans la zone 2 une plaque d'acier.
<figure><img src="image_6.png" /><figcaption>Fig. 6. Conductance thermique de contact, cas A.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Fig. 7. Conductance thermique de contact, cas B.</figcaption></figure>

Ces tests numériques permettent de montrer l'efficacité de la méthode avec différents contacts thermiques.

Le cas A représente le cas avec un défaut moyen dans l'assemblage, situé au centre de la géométrie (Fig. 6).

Le cas B représente un cas avec deux défauts situés en $y^* = 30$ et $y^* = 70$ (Fig. 7).

Le cas C représente un cas avec une conductance thermique de contact plus complexe, un cas plus réaliste (Fig. 8).

Tous ces tests ont été réalisés avec comme critère d'arrêt $\epsilon = 1e-5$ sur un ordinateur équipé d'un processeur Intel(R) Pentium(R) M de 1,70 GHz avec 1 Go de mémoire vive.

Pour vérifier la qualité des estimations, une erreur d'estimation générale a été calculée pour chaque test (Tab. 4). Cette erreur permet de quantifier la précision de l'estimation. Elle est calculée avec l'équation (50) :

$$E = \frac{\sum_{i=1}^{ny} (H_{réel}^{*} - H^{*})^2}{\sum_{i=1}^{ny} (H_{réel}^{*})^2} \quad (50)$$

<figure><img src="image_8.png" /><figcaption>Fig. 8. Conductance thermique de contact, cas C.</figcaption></figure>

<table><caption>Tableau 4. Résumé des tests d'évaluations A, B et C.</caption><thead><tr><th>Cas</th><th>A</th><th>B</th><th>C</th></tr></thead><tbody><tr><td>Nombre itération</td><td>103</td><td>85</td><td>103</td></tr><tr><td>Érreur d'évaluation générale en %</td><td>0,9</td><td>1,6</td><td>1,7</td></tr><tr><td>Temps de calcul (s)</td><td>2103</td><td>1985</td><td>1873</td></tr></tbody></table>

Une étude préliminaire a permis de définir une limite d'erreur d'estimation générale maximale de 5 %. En effet, avec une erreur supérieure à 5 %, la conductance thermique de contact estimée est trop éloignée de la conductance thermique recherchée. Cette limite permet de réaliser une étude sur la qualité de l'estimation.

Une erreur d'estimation ponctuelle sera aussi calculée à l'aide de l'équation (51), pour observer la précision de l'estimation au niveau des défauts de contact thermique.

$$Ep(i) = \frac{(H_{réel}^{*}(i) - H^{*}(i))^2}{(H_{réel}^{*}(i))^2} \quad (51)$$

La figure 9 représente la conductance de contact réelle et celle retrouvée suite à l'estimation du cas A, présentant un défaut moyen dans l'assemblage. La forme générale est donc retrouvée avec précision, le défaut est identifié à la bonne position et à la bonne valeur. En effet, d'après la figure 9, l'erreur d'estimation ponctuelle reste inférieure à 3 % sur l'ensemble du contact thermique. Au niveau du défaut, entre $y^* = 40$ et $y^* = 70$, elle reste inférieure à 1 %. Ce cas est donc fonctionnel, car l'erreur d'estimation générale est inférieure à 1 %, voir tableau 4, ce qui est bien inférieur à la limite de 5 % fixée ci-dessus.

La figure 10 représente les résultats de l'estimation du cas B, avec deux défauts. Comme précédemment la forme générale, bien que plus complexe, a été retrouvée. Les deux défauts sont positionnés au même endroit que sur la conductance thermique de contact réelle et leurs valeurs sont bien estimées. Dans ce test, l'erreur d'estimation générale est supérieure au cas précédent, 1,7 % au lieu de 0,9 % (voir Tab. 4) mais reste inférieure à la limite d'erreur maximale tolérée dans cette étude, 5 %.
<figure><img src="image_9.png" /><figcaption>Fig. 9. Estimation de la conductance thermique de contact, cas A.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Fig. 10. Estimation de la conductance thermique de contact, cas B.</figcaption></figure>

La figure 11 représente les résultats de l'estimation de conductance thermique de contact du cas C, cas avec une conductance thermique de contact de forme stochastique. L'estimation a fonctionné malgré la complexité, les défauts étant localisés et estimés avec précision. L'erreur d'estimation générale reste inférieure à 2 % (voir Tab. 4).

Les figures 9–11 montrent que les erreurs d'estimation les plus importantes sont dans les zones où la conductance thermique de contact varie rapidement, les défauts sont quant à eux bien localisés. Une conductance thermique de contact ayant des variations importantes sera donc plus difficile à estimer.

Enfin, il est intéressant de constater que, malgré une complexité différente, ces trois tests possèdent des

<figure><img src="image_11.png" /><figcaption>Fig. 11. Estimation de la conductance thermique de contact, cas C.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Fig. 12. Variation de l'erreur d'estimation globale en fonction de Biot et de $k^*$ avec le cas A.</figcaption></figure>

temps de résolution quasi-identiques, entre 2103 et 1873 s (Tab. 4).

## 4.2 Influence des paramètres

Suite à la validation de la méthode, des tests ont été réalisés pour observer l'influence de deux paramètres adimensionnés, $k^*$ et Biot, sur l'estimation. La géométrie est fixée, $L_1^* = 150$ et $L_y^* = 100$. La forme de la conductance de contact adimensionnée, $H_{réel}^*$, est celle du cas A avec un défaut, présentée à la figure 6. La figure 12 représente l'erreur d'estimation générale en fonction du nombre de Biot et de $k^*$. Elle montre que, pour avoir une estimation la plus précise possible, il faut avoir une combinaison d'un nombre de Biot et $k^*$ élevés.
<figure><img src="image_13.png" /><figcaption>Fig. 13. Variation de l'erreur d'estimation globale en fonction de Biot et de $k^*$ avec le cas C.</figcaption></figure>

Ces résultats s'expliquent par le fait que, dans cette situation, la sensibilité du profil de température mesuré en $x^* = L^*$ par rapport au profil de conductance de contact thermique est grande. En effet, avec un rapport $k^*$ grand, la diffusion de chaleur de $x^* = 0$ à $x^* = L^*$ sera importante, l'estimation est donc plus facilement réalisable. Dans le cas contraire, avec un $k^*$ petit, la conduction de chaleur dans la zone 1 est faible, voir équation (15), et donc l'estimation est plus difficile à réaliser. Ceci explique les erreurs importantes pour $k^* = 0,025$.

De même, si le refroidissement dû à la convection extérieure est important, c'est à dire pour un Biot élevé, alors les points chauds seront plus visibles et donc plus facilement identifiables. Le lien entre le profil de température mesuré et la nature du contact thermique est alors plus direct.

La figure 12 permet aussi de constater qu'avec un nombre de Biot inférieur à 0,0025, même avec un nombre $k^*$ élevé, l'erreur d'estimation générale est supérieure à la limite d'erreur d'estimation. Le cas d'une convection naturelle est donc à éviter lors d'une estimation. En effet, dans ce cas particulier le nombre de Biot est faible, de l'ordre de 0,0025, car le transfert thermique par convection est faible par rapport au transfert par conduction.

Enfin, pour obtenir une erreur inférieure à 5 %, erreur maximale tolérée, $k^*$ doit être supérieur à 0,25 et Biot doit être supérieur ou égal à 0,05 dans ce cas.

Une étude similaire a été réalisée (voir la Fig. 13), avec la conductance thermique de contact du cas C possédant une forme plus complexe, tel que présenté à la figure 8.

La figure 13 permet d'observer que, dans un cas plus complexe, il est nécessaire d'avoir $k^*$ supérieur à 1 et Biot supérieur à 0,005.

Ces tests permettent de conclure que, plus la complexité du test augmente, plus le nombre de Biot et le rapport $k^*$ doivent être grands pour diminuer le plus possible l'erreur d'estimation.

## 5 Conclusion

Une méthode totalement novatrice utilisant les méthodes inverses réalisant l'estimation de la conductance thermique de contact a été développée. Les tests numériques ont permis de montrer l'intérêt de cette méthode de calcul dans le cas de l'estimation d'un contact thermique dans un réacteur métallurgique.

Des tests numériques ont été réalisés avec plusieurs formes de conductance thermique de contact. L'influence des paramètres $k^*$ et Biot a également été présentée avec une géométrie donnée pour deux formes de conductance thermique de contact.

Cette étude permet de conclure que l'estimation d'une conductance thermique de contact de forme complexe est réalisable dans le cas d'un réacteur métallurgique. Elle a aussi permis de conclure que, pour obtenir une estimation précise dans le cas présenté, il faut avoir simultanément un $k^*$ supérieur à 1 et un nombre de Biot supérieur à 0,005. Cette méthode peut être utilisée dans d'autres applications industrielles telles que l'estimation du contact thermique lors de l'assemblage de puces électroniques.

*Remerciements.* Les auteurs remercient le Conseil de recherches en sciences naturelles et en génie du Canada pour son support financier.

## Références

[1] A.S. Marchand, Étude de la résistance thermique de contact à l'interface de solides déformables en frottement. Application aux procédés de forgeage, Ph.D. thesis, Institut National des Sciences Appliquées de Lyon, 1998

[2] K. Atchonouglo, M. Banna, C. Vallee, J.C. Dupre, "Inverse transient heat conduction problems and identification of thermal parameters," Heat and Mass Transfer/Waerme- und Stoffuebertragung 45 (2008) 23–29

[3] C. Fieberg, R. Kneer, "Determination of thermal contact resistance from transient temperature measurements," Int. J. Heat Mass Transfer 51 (2008) 1017–1023

[4] M.H. Shojaeefard, K. Goudarzi, M.S. Mazidi, "Inverse heat transfer problem of thermal contact conductance estimation in periodically contacting surfaces," J. Thermal Sci. 18 (2009) 150–159

[5] J. Gill, E. Divo, A.J. Kassab, "Estimating thermal contact resistance using sensitivity analysis and regularization," Eng. Anal. Bound. Elem. 3 (2009) 54–62

[6] V. Feuillet, Développement d'outils d'analyse thermique pour la conception de composant électronique de puisance, Ph.D. thesis, Université de Nantes, 2006

[7] M.G. Cooper, B.B. Mikic, M.M. Yovanovich, "Thermal contact conductance," Int. J. Heat Mass Transfer 12 (1969) 279–300

[8] R. Perez, "Heat transfer; le transfert thermique," Electricité de France, Bulletin de la Direction des Etudes et Recherches, Serie A : Nucleaire, Hydraulique, Thermique 1-2 (1973) 113–188

[9] W.M. Rohsenow, J.P. Hartnett, Y.I., Handbook of heat transfer, 1998
[10] J.V. Beck, "Determination of optimum, transient experiments for thermal contact conductance," Int. J. Heat Mass Transfer 12 (1969) 621-33

[11] J.V. Beck, Y. Jarny, D.G. Inverse Engineering Handbook, 2003

[12] D. Petit, D. Maillet, "Techniques inverses et estimation de paramètres", Technique de l'ingenieur, 1995

[13] M.J. Colaco, H.R.B. Orlande, "Inverse problems in science and engineering : Foreword", Inverse Problems in Science and Engineering 17 (2009) 1

[14] M.N. Ozisik, H.R.B. Orlande, "Inverse heat transfert", 1999

[15] Y. Jarny, M.N. Ozisik, J.P. Bardon, "A general optimization method using adjoint equation for solving multi-dimensional inverse heat conduction", Int. J. Heat Mass Transfer 34 (1991) 2911-2919

[16] R. Fletcher, C.M. Reeves, "Function minimization by conjugate gradients", Comput. J. 7 (1964) 149-154