# Etude numérique de la fusion des matériaux à changement de phase pour le stockage thermique dans une cavité rectangulaire inclinée

Mohamed Serradj, Abdeldjellil Benmansour * et Abdewaheb Benzerdjeb

Laboratoire de Mécanique Appliquée
Université des Sciences et de la Technologie d'Oran Mohamed Boudiaf
B.P. 1505 El-M'Naouer, Oran, Algérie

(reçu le 22 Août 2015 – accepté le 29 Septembre 2015)

**Résumé** – Dans ce travail, une étude numérique bidimensionnelle du processus de la fusion d'un matériau à changement de phase (MCP) placé à l'intérieur d'une enceinte rectangulaire inclinée, dont toutes les faces sont isolées thermiquement à l'exception de la face verticale avant, qui est chauffée de manière isotherme. On a considéré que la conduction était le mode de transfert de chaleur qui régit le processus au début de la fusion, suivie par une courte période de transition avant que la convection naturelle devient le mode de transfert de chaleur dominant dans la phase liquide. Les équations qui gouvernent ce genre de phénomène sont résolues et une solution numérique est développée. Les effets de l'angle d'inclinaison et l'orientation de la cavité sur le phénomène de la fusion sont examinés.

**Abstract** – In this work, a two-dimensional numerical analysis of the melting process of a phase change material (PCM) disposed within a rectangular inclined enclosure. All faces are thermally insulated except the vertical front face which is heated isothermally. It was considered that the conduction mode was the heat transfer governing the process at the beginning of the fusion, followed by a short transition period before natural convection mode becomes dominant heat transfer in the liquid phase. The equations that govern this kind of phenomenon are solved and a numerical solution is developed. The effects of the tilt angle and orientation of the cavity on the phenomena of fusion were examined.

**Mot clés:** Fusion – MCP – Cavité inclinée – Stockage thermique.

## 1. INTRODUCTION

Beaucoup de problèmes surviennent dans lesquels le transfert de chaleur est accompagné par un changement de phase, où l'énergie thermique est absorbée ou libérée. Les exemples les plus remarquables sont les processus de fusion et de solidification dans lesquels une interface mobile sépare deux régions de propriétés thermo physiques différentes, à savoir les phases solide et liquide.

De tels problèmes de transfert de chaleur transitoires sont communément appelés changement de phase, frontière mobile, ou quelque fois problème de Stefan.

Plusieurs travaux ont été menés afin d'étudier le processus de la fusion et la solidification des matériaux à changement de phase dans une cavité rectangulaire verticale [1, 2]. Ces travaux révèlent que jusqu'à un passé assez récent, la conduction pure était considérée comme étant le seul mode de transfert de chaleur dans les problèmes de changement de phase. Or des études expérimentales ont montrées que seulement au début du processus de la fusion que le transfert de chaleur dans la phase liquide est contrôlé par la conduction, suivie par une courte période de transition. Puis comme le volume fondu augmente, la convection naturelle devient prédominante sur la conduction dans le liquide jusqu'à la fin du processus de la fusion.

<footnote>* abd_benmansour@yahoo.fr</footnote>
Dans le présent travail, on s’intéresse aux effets de l’angle d’inclinaison et de l’orientation d’une cavité rectangulaire sur l’écoulement et le transfert de chaleur lors de la fusion d’un matériau à changement de phase contenu dans cette cavité.

## 2. MODELE PHYSIQUE

Le système physique considéré dans ce travail est identique à celui déjà étudié expérimentalement par Webb et al. [3]. Il s’agit d’une enceinte de forme rectangulaire de hauteur H = 25.4 cm et d’épaisseur E = 6.35 cm, remplie de N-octadécane comme matériau à changement de phase, dont toutes les faces sont isolées thermiquement à l’exception de la face avant qui est chauffée de manière isotherme, (figure 1).

<figure><img src="image_1.png" /><figcaption>Fig. 1: Modèle physique</figcaption></figure>

## 3. MODELE MATHEMATIQUE

Dans cette partie, nous allons définir le modèle mathématique qui est utilisé pour résoudre ce type de problèmes. C’est un modèle où tous les champs (u, v, p, T) sont déterminés.

### 3.1 Hypothèses simplificatrices

Nous considérerons que l’écoulement dans la phase liquide est laminaire, incompressible et newtonien. Nous négligerons la dissipation visqueuse. Les propriétés physiques du matériau seront considérées constantes et identiques dans les deux phases. La variation de volume, résultant du changement de phase, est négligée.

De plus, nous utiliserons l’approximation de Boussinesq. Celle-ci implique que la densité ρ est considérée constante hormis dans le terme source de gravité, qui déclenche la convection naturelle dans ‘équation de quantité de mouvement.’

### 3.2 Equations de conservation

En tenant compte des hypothèses posées précédemment, nous pouvons écrire le système d’équation suivant:

Conservation de la masse-

$$ \frac{\partial u}{\partial x} + \frac{\partial u}{\partial y} = 0 \qquad (1) $$

Conservation de la quantité de mouvement suivant x -
$$ \rho \cdot \left( \frac{\partial u}{\partial t} + u \cdot \frac{\partial u}{\partial x} + v \cdot \frac{\partial u}{\partial y} \right) = - \frac{\partial P}{\partial x} + \mu \left( \frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2} \right) - \rho \cdot g \cdot \beta (T - T_{\infty}) \sin \varphi + \frac{(1 - f_1)^2}{(f_1^3 - \xi)} \cdot Au \quad (2) $$

Conservation de la quantité de mouvement suivant y -

$$ \rho \cdot \left( \frac{\partial v}{\partial t} + u \cdot \frac{\partial v}{\partial x} + v \cdot \frac{\partial v}{\partial y} \right) = - \frac{\partial P}{\partial y} + \mu \left( \frac{\partial^2 v}{\partial x^2} + \frac{\partial^2 v}{\partial y^2} \right) - \rho \cdot g \cdot \beta (T - T_{\infty}) \cos \varphi + \frac{(1 - f_1)^2}{(f_1^3 - \xi)} \cdot Av \quad (3) $$

Conservation de l'énergie-

$$ \frac{\partial(\rho h)}{\partial t} + u \cdot \frac{\partial(\rho h)}{\partial t} \frac{\partial v}{\partial y} + v \cdot \frac{\partial(\rho h)}{\partial y} = k \cdot \left( \frac{\partial^2 T}{\partial x^2} + \frac{\partial^2 T}{\partial y^2} \right) \frac{\partial v}{\partial y} \quad (4) $$

A est la constante de la zone molle qui quantifie la perte de vitesse dans la zone où a lieu le changement de phase. La valeur de A varie, selon les cas traités, entre 10⁴ et 10⁷ [4, 5]. En se référant à ces travaux sur un problème similaire, nous avons fixé sa valeur dans ce travail à 1.6×10⁶.

ζ est un nombre réel de faible valeur (~10⁻⁴) permettant d'éviter la division par zéro. Par ailleurs, nous allons utiliser la méthode 'enthalpie-porosité' pour prendre en compte le changement de phase de N-octadécane. C'est-à-dire que nous définirons l'enthalpie totale h comme étant la somme de l'enthalpie sensible et de la chaleur latente Lf.

$$ h = h_s + L_f = \int C_p d t + L_f \quad (5) $$

La fraction liquide locale ($f_1$) correspond à la chaleur latente absorbée localement ($\lambda$). Avec,

$$ f_1 = \begin{cases} 0 & (T < T_s) \\ \frac{T - T_s}{T_1 - T_s} & (T_s \le T \le T_1) \\ 0 & (T > T_1) \end{cases} \quad (6) $$

On utilise les équations (5) et (6), pour réécrire l'équation (4),

$$ \rho \cdot C_p \left( \frac{\partial T}{\partial t} + u \frac{\partial T}{\partial x} + v \frac{\partial T}{\partial y} \right) = k \cdot \left( \frac{\partial^2 T}{\partial x^2} + \frac{\partial^2 T}{\partial y^2} \right) - \lambda \cdot \rho \cdot \left( \frac{\partial f_1}{\partial t} + u \frac{\partial f_1}{\partial x} + v \frac{\partial f_1}{\partial y} \right) \quad (7) $$

### 3.3 Modèle numérique (CFD)

Nous allons utiliser le logiciel FLUENT® pour créer une base de données. Ce code commercial est basé sur la méthode des volumes finis qui consiste à discrétiser le domaine de l'écoulement, après avoir effectué sur chaque volume de contrôle les bilans de masse, de quantité de mouvement et de l'énergie.

Les résultats seront présentés comme un ensemble d'équations algébriques (matrices) que l'on résout de manière itérative.

A partir des schémas d'interpolations, choisis dans les travaux de [4-6], nous paramétrons nos modèles en utilisant le schéma d'interpolation Quick pour les termes convectifs des équations de conservation.

Le calcul de la pression se fera à l'aide du schéma Presto qui est bien adapté au gradient de chute de pression due au mouvement de rotation de la phase liquide. Le
couplage entre la vitesse et la pression se fera à l'aide de l'algorithme SIMPLE développé par Patankar et al. [4].

Toujours en nous référant aux travaux de [5, 6], nous prendrons les valeurs suivantes: les résidus de vitesse en x et y seront fixés à $10^{-3}$, ceux du balancement d'énergie seront fixés à $10^{-6}$. Enfin la masse locale sera de $10^{-3}$ par rapport à la masse totale de la cavité.

# 4. RESULTATS

## 4.1 Géométrie et maillage du système

On va étudier la fusion d'un MCP dans une cavité rectangulaire. La cavité a une hauteur H et une épaisseur E avec un rapport $H/E = 4$. Nous utilisons un maillage de type structuré (51 × 81) qui nous permet d'étudier l'effet de l'inclinaison sur la fusion. Toutes les mailles sont de même taille. Le maillage du domaine est représenté sur la figure 2.

<figure><img src="image_2.png" /><figcaption>Fig. 2: Géométrie et maillage</figcaption></figure>

## 4.2 Conditions aux limites

<table><caption>Tableau 1: Conditions aux limites</caption><thead><tr><th>Paroi</th><th>Conditions aux limites</th></tr></thead><tbody><tr><td>Paroi gauche</td><td>$T = T_c > T_f$, $u = v = 0$</td></tr><tr><td>Autres parois</td><td>$\frac{\partial T}{\partial y} = 0$, $u = v = 0$</td></tr><tr><td>Intérieur de la cavité</td><td>$T = T_i < T_m$</td></tr></tbody></table>

Où, $T_f = 301 \text{ K}$, $T_i = 292 \text{ K}$, $T_c = 312 \text{ K}$, $Pr=51$, $Ste=0.105$, $Sc=0.08$, $Ra = 3.87 \times 10^9$.

## 4.3 Effet de l'inclinaison sur la fusion

Pour étudier l'effet de l'inclinaison de cavité sur la fusion du MCP, on fait varier l'angle d'inclinaison $\phi$ de 0° à 90° et visualiser la forme et la position pour le même intervalle de temps (figure 3).
<figure><img src="image_3.png" /><figcaption>Fig. 3: Evolution de l'interface solide-liquide en fonction du temps pour différentes valeurs de φ</figcaption></figure>

Ces figures montrent clairement que le taux de fusion augmente avec l'angle d'inclinaison de la cavité. L'inclinaison d'un angle de 90° représente la limite et un cas particulier. Dans ce cas, l'évolution du processus de la fusion paraît très lente à cause de l'absence de la convection naturelle pendant longtemps où la forme de l'interface solide-liquide reste parallèle à la paroi chaude.

<figure><img src="image_4.png" /><figcaption>Fig. 4: Evolution de la fraction liquide en fonction du temps pour différentes valeurs de φ</figcaption></figure>

La figure 4 illustre l'effet de l'angle d'inclinaison sur la fraction liquide. Il apparaît clairement que pour une inclinaison supérieure à 45°, la fraction liquide augmente avec l'angle d'inclinaison. Par contre pour un φ inférieur à 45°, elle reste insignifiante.

<figure><img src="image_5.png" /><figcaption>Fig. 5: Evolution de la fraction liquide en fonction du temps pour une cavité verticale et une cavité inclinée</figcaption></figure>
## 4.4 Effet de l'orientation de cavité sur la fusion

Pour étudier l'effet de l'orientation de la cavité inclinée sur le phénomène de la fusion, on a incliné la cavité d'un angle de 45° et on a commencé par la chauffer d'abord du coté gauche où le MCP liquide se trouve au dessus du MCP encore à l'état solide, puis du côté de la paroi de droite de la cavité où le MCP liquide se trouve cette fois-ci au dessous du MCP à l'état solide. Nous avons gardé les mêmes dimensions et les mêmes conditions de températures.

Les figures 5a- et 5b- illustrent l'évolution de la fraction liquide en fonction du temps pour une cavité verticale et une cavité inclinée et chauffée d'abord par le haut puis par le bas.

Ces résultats montrent que l'inclinaison de la cavité peut, selon le cas soit augmenter ou diminuer le taux de fusion, donc améliorer ou retarder la fusion. Afin d'augmenter la partie fondue, le sens de l'inclinaison doit être choisi de telle façon que le MCP à l'état liquide se trouve au-dessous du MCP encore à l'état solide. Aussi pour une cavité inclinée et chauffée par le bas, la vitesse du liquide est plus élevée par rapport au cas où une cavité inclinée est chauffée par le haut, surtout lorsque la fusion atteint un stade très avancé.

## 5. CONCLUSION

Les résultats de ce travail montrent que l'inclinaison de la cavité joue un rôle très important, car elle influe directement sur la morphologie de l'écoulement du matériau à changement de phase liquide à l'intérieur de la cavité où la convection naturelle est le mode de transfert de chaleur qui régit le processus de la fusion. On constate que plus on n'incline la cavité, plus la fusion s'accélère et le taux de fusion augmente.

D'autre part, l'inclinaison elle-même peut retarder la fusion et diminuer le taux de la fusion. Cette inclinaison peut donc avoir un rôle inverse. De même que le choix du sens de l'orientation de la cavité peut soit accélérer la fusion ou la retarder.

## NOMENCLATURE

<table><tr><td>A, Constante de la zone molle</td><td>C<sub>p</sub>, Chaleur massique, J/kgK</td></tr><tr><td>E, Epaisseur de cavité, m</td><td>F<sub>l</sub>, Fraction liquide locale</td></tr><tr><td>H, Hauteur de la cavité, m</td><td>h, Enthalpie totale, J/kg</td></tr><tr><td>k, Conductivité thermique, W/mK</td><td>h<sub>s</sub>, Enthalpie sensible, J/kg</td></tr><tr><td>L, Chaleur latente, kJ/kg</td><td>MCP, Matériau à changement de phase</td></tr><tr><td>m, Masse, kg</td><td>P, Pression, Pa</td></tr><tr><td>T<sub>f</sub>, Température de fusion, K</td><td>T<sub>c</sub>, Température de paroi chaude, K</td></tr><tr><td>T<sub>s</sub>, Température de solide, K</td><td>T<sub>l</sub>, Température de liquide, K</td></tr><tr><td>u, Vitesse suivant x, m/s</td><td>v, Vitesse suivant y, m/s</td></tr><tr><td>Ra, Nombre de Rayleigh</td><td>Ste, Nombre de Stefan</td></tr><tr><td>V, Volume, m<sup>3</sup></td><td>Sc, Nombre de sous refroidissement</td></tr><tr><td>φ, Densité de flux de chaleur, W/m<sup>2</sup></td><td>ρ, Masse volumique, kg/m<sup>3</sup></td></tr><tr><td>λ, Ratio - Chaleur latente absorbée<br>localement</td><td>β, Coefficient d'expansion, 1/K</td></tr></table>
REFERENCES

[1] M. Serradj, ‘Contribution à l'Etude du Phénomène de la Fusion des Matériaux à Changement de Phase', Mémoire de Magister, USTO M.B, Oran, 2013.

[2] M.A Hamdan and F.A Elwerr, ‘Thermal Energy Storage Using a Phase Change Material', Solar Energy, Vol. 56, N°2, pp. 183 - 189, 1996.

[3] B.W. Webb and R. Viskanta, ‘Natural-Convection-Dominated Melting Heat Transfer in an Inclined Rectangular Enclosure', International Journal of Heat and Mass Transfer, Vol. 29, N°2, pp. 183 - 192, 1986.

[4] G. Vidalain, ‘Modélisation des Phénomènes Convectifs lors du Changement de Phase Solide Liquide par Utilisation de l'Equation de Diffusion de la Chaleur et d'une Forme Modifiée de la Conductivité', Thèse Maître Es Sciences, Université de Laval, Quebec, 2007.

[5] A.D. Brent, V.R. Voller and K.J. Reid, ‘The Entalpy-Porosity Technique for Modeling Convection Diffusion Phase Change: Application to the Melting of Pure Metal', Numerical Heat Transfer, Part B fundamentals, Vol. 13, N°3, pp. 297 - 318, 1992.

[6] R. Ankri et T. Zebbiche, ‘Résolution d'un Problème de Changement de Phase en Présence de la Convection Naturelle par la Méthode Enthalpique: Application au Gallium', Revue des Energies Renouvelables, CISM'08 Oum El Bouaghi, pp. 25 – 36. 2008.