# Simulation du Stockage de l'Energie Thermique dans un Lit Fixe de Sphères Contenant un Matériau à Changement de Phase

Abdeldjellil Benmansour¹ et Mohammad Hamdan²

¹ Département de Chimie, Faculté des Sciences, Université des Sciences et de la Technologie Mohamed Boudiaf d'Oran B.P. 1505 El M'NAOUAR ORAN, ALGERIE

² Mech. Eng. Dept., Faculty of Engineering and Technology, Jordan University P.O. Box 13240 AMMAN, JORDAN

**Résumé** – Ce travail présente une étude numérique du stockage de l'énergie thermique par chaleur latente à basses températures (0 °C à 100 °C) dans un lit cylindrique rempli de sphères uniformes, disposées au hasard et contenant chacune un matériau à changement de phase (MCP), traversé par un flux d'air. Un modèle numérique mono dimensionnel à deux phases séparées est appliqué. Il a permis de prédire la distribution axiale de la température du fluide et du matériau fusible le long du lit ainsi que les performances de ce lit fixe dans les deux modes de stockage et de récupération de la chaleur pour une température du fluide à l'entrée constante et un nombre de Reynolds variable. Le modèle développé s'applique à la fois pour un processus de changement de phase isotherme et non isotherme. L'évolution de la température à l'intérieur du lit a permis de déterminer les domaines de stockage les plus efficaces et de définir par conséquent les dimensions optimales de l'unité de stockage.

**Abstract** – This work presents a numerical study of a latent thermal energy storage for temperature 0 °C to 100 °C in a cylindrical randomly packed bed of uniform spheres containing a phase change material (PCM) subject to air flow. A one dimensional separate phase model is applied. It allows to predict the bad performance as well as fluid and PCM temperatures distribution within the packed bed in the axial direction for charging and recovery modes at constant inlet fluid temperature and variable Reynolds number during both isothermal and nonisothermal phase change process. Temperature evolution of a packed bed allows to determine the efficient storage domain and to define an optimization for storage units.

**Mots clés:** Stockage d'énergie - Chaleur latente - Changement de phase - Sphère.

## 1. INTRODUCTION

La nature intermittente de l'énergie solaire pose des problèmes de stockage. Il est nécessaire de développer des dispositifs de stockage permettant le transfert d'énergie des périodes excédentaires aux périodes déficitaires et adapter ainsi sa production à sa consommation. Ce stockage peut être sous forme de chaleur sensible ou de chaleur latente. Différents types d'unités de stockage ont été conçus et étudiés [1,2]. Ils diffèrent tous par la nature du matériau servant de support de stockage et la géométrie du lit.

Le stockage thermique par chaleur latente dans les liquides ou les solides peut être un moyen très intéressant en utilisant des matériaux à changement de phase. Ces matériaux présentent l'avantage de stocker une grande quantité d'énergie avec un bon rapport masse/volume.

Ainsi, et depuis quelques années, l'analyse et la modélisation du phénomène de stockage de l'énergie thermique par la chaleur latente ont fait l'objet d'un certain nombre de travaux que résume Adebiyi dans son article [3].

Dans des travaux antérieures [4,5], nous avons procédé à une analyse de ce type stockage et nous avons mis en évidence l'influence des paramètres dynamiques et thermiques du fluide sur l'évolution de la température du MCP le long du lit, dans les deux modes de stockage et
de récupération de la chaleur, dans le cas d'un processus de changement de phase à la fois isotherme et non isotherme. Dans ce travail, notre contribution consiste à évaluer les performances (durée de stockage, efficacité) d'un lit fixe ainsi que la distribution axiale de la température du matériau fusible dans ce lit de stockage constitué de sphères uniformes, disposées au hasard et contenant chacune un MCP, dans les deux modes de stockage et de récupération de la chaleur. Pour cela, on utilise un modèle mathématique à deux phases séparées (où les températures instantanées du fluide et du matériau sont supposées différentes) et unidimensionnel suivant l'axe parallèle à la vitesse d'écoulement du fluide. Après avoir introduit certaines hypothèses simplificatrices, les équations aux dérivées partielles régissant les phénomènes de transfert de chaleur sont résolues par une méthode de différences finies.

## 2. FORMULATION MATHEMATIQUE

### 2.1. Hypothèses

Le lit de stockage est représenté par le schéma suivant:

<figure><img src="image_1.png" /><figcaption>Fig. 1: Schéma du lit de stockage</figcaption></figure>

Le modèle étudié est basé sur les hypothèses suivantes:

* les parois du lit sont adiabatiques,

* l'influence de la paroi du cylindre est négligée,

* le phénomène est unidirectionnel,

* la vitesse moyenne d'écoulement du fluide est constante dans tout le lit,

* la conduction entre les particules sphériques est ignorée,

* les gradients de températures sont faibles à l'intérieur des particules sphériques.

### 2.2. Equations simplifiées

Le système proposé est un cylindre orienté verticalement contenant des particules sphériques remplies chacune de MCP et traversé par un flux de chaleur. Les équations d'énergie qui régissent le phénomène de stockage / déstockage de la chaleur pour un modèle à deux phases séparées s'écrivent[2]:

- pour le fluide:

$$ \frac{\partial T}{\partial t} + u \frac{\partial T}{\partial x} = \frac{k}{\rho_F c_F} \frac{\partial^2 T}{\partial x^2} + \frac{ha}{\rho_F c_F \varepsilon} (\Theta - T) \quad (1) $$
- pour le MCP:

$$ \rho_s c_s (1-\varepsilon) \frac{\partial \Theta}{\partial t} = h a (T - \Theta) \quad \text{phase solide} \quad (\Phi=0) \quad (2) $$

$$ \rho_s h_l (1-\varepsilon) \frac{\partial \Phi}{\partial t} = h a (T - \Theta_s) \quad \text{deux phases} \quad (0 < \Phi < 1) \quad (3) $$

$$ \rho_s c_l (1-\varepsilon) \frac{\partial \Theta}{\partial t} = h a (T - \Theta) \quad \text{phase liquide} \quad (\Phi=1) \quad (4) $$

Ces équations ne tiennent compte que de la dispersion thermique axiale.

A ce système d'équations aux dérivées partielles, nous devons associer des conditions initiales et aux limites

## 2.3. Les conditions initiales et aux limites

On suppose qu'au début du stockage (ou du déstockage), les températures du fluide et du MCP sont en équilibre avec le milieu ambiant. Ceci se traduit par:

$$ T(x, t=0) = \Theta(x, t=0) = T_a $$

A l'entrée du lit les températures du fluide et du MCP sont en équilibre avec la température de l'air. Donc :

$$ T(x=0, t) = \Theta(x=0, t) = T_e $$

A la sortie du lit de stockage, on suppose que les températures pour $X \ge H$ restent constantes. On aura donc:

$$ \frac{\partial T}{\partial x} (x = H, t) = \frac{\partial \Theta}{\partial x} (x = H, t) = 0 $$

## 2.4. Les performances du lit

Les performances caractérisant un lit fixe de stockage de l'énergie thermique peuvent être représentées par la durée du stockage et l'efficacité du lit.

a - la durée du stockage $\tau$ est définie comme le temps nécessaire au lit pour atteindre l'état de saturation thermique.

b - l'efficacité du lit $\pi$ est définie par :

$$ \pi_c (t) = \Theta - T_a / T_e - T_a \quad \text{en mode charge} $$

$$ \pi_r (t) = \Theta - T_a / T_s - T_a \quad \text{en mode récupération} $$

## 3. METHODE DE RESOLUTION

Pour la résolution des équations régissant le transfert de chaleur à l'intérieur du lit, nous avons utilisé la méthode des différences finies. La discrétisation des équations à partir d'un maillage spatial et temporel de pas respectivement égaux à $\Delta x$ et $\Delta t$ permet d'écrire en utilisant la méthode purement implicite des différences finies [6], pour le MCP:

$$ \Theta_i^{j+1} = (\Theta_i^j + \omega \Delta t T_i^{j+1}) / (1 + \omega \Delta t) \quad (5) $$
où:

$$ \omega = ha / \rho_s c_s (1 - \varepsilon) \quad \text{phase solide} $$

$$ \omega = ha / \rho_s c_1 (1 - \varepsilon) \quad \text{phase liquide} $$

Dans la région des deux phases:

- avec un processus de changement de phase isotherme:

$$ \omega = 0, \quad \Theta = \Theta_s $$

et

$$ \Phi_{i}^{j+1} = \omega' \Delta t (T_{i}^{j} - \Theta_{s}) + \Phi_{i}^{j} \quad (5') $$

où:

$$ \omega' = ha / \rho_s h_1 (1 - \varepsilon) $$

- avec un processus de changement de phase non isotherme:

$$ \omega = ha / \rho_s c_p (1 - \varepsilon) $$

En substituant l'équation (5) dans l'équation (1) écrite sous forme de différences finies, on obtient:

$$ -(2A+N)T_{i-1}^{j+1} + ET_{i}^{j+1} - NT_{i+1}^{j+1} = T_{i}^{j} + C\Theta_{i}^{j} \quad (6) $$

où on a posé:

$$ A = u \Delta t / 2 \Delta x $$

$$ B = ha \Delta t / \rho_F c_F $$

$$ C = B / (1 + \omega \Delta t) $$

$$ D = B \omega \Delta t / (1 + \omega \Delta t) $$

$$ E = 1 + 2A + B + 2N - D $$

$$ N = k \Delta t / \rho_F c_F (\Delta x)^2 $$

A partir de ce système, on calcule $T_{i}^{j+1}$ connaissant $T_{i}^{j}$ et $\Theta_{i}^{j}$ et ceci pour tout i. Puis on calcule $\Theta_{i}^{j+1}$ à partir de l'équation (5) pour tout i.

## 4. RESULTATS ET DISCUSSION

L'étude numérique a été faite dans le cas d'un lit cylindrique de 0.90 m de haut et 0.20 m de diamètre, rempli de sphères en plastiques de diamètre égal à 0.038 m contenant chacune de la cire de paraffine (domaine de fusion : 58°C - 60°C) et traversé par de l'air préalablement chauffé et maintenu à une température d'entrée constante et supérieure à la température de fusion du matériau choisi.

Nous avons représenté sur les Figures 2 (a), 2 (b) et 2 (c) l'efficacité du lit durant la période de stockage pour un rapport hauteur/diamètre du lit (H/D) variant de 1 à 4 et un débit d'air (Re) variable. La tendance générale de toutes ces Figures est que l'efficacité du lit durant le stockage diminue avec H/D. C'est dû au fait que le volume de stockage augmente avec H/D et par conséquent plus de chaleur est nécessaire au lit pour augmenter sa température et son efficacité de stockage.

Aussi, on peut remarquer que pour un temps donné et une même valeur de H/D, l'efficacité augmente avec Re ( débit d'air ). De plus ces figures indiquent que pour toutes les
valeurs de H/D, la durée du stockage (temps nécessaire au lit pour atteindre sa température maximum, correspondant à l'état de saturation thermique) diminue avec Re ( débit d'air ).

<figure><img src="image_2.png" /><figcaption>Fig. 2: Variation de l'efficacité du lit en fonction du temps, durant le stockage, pour différents Re et H/D.</figcaption></figure>

C'est dû au fait que comme plus d'air chaud est introduit dans le lit, plus de chaleur est stockée dans le système et ce qui par conséquent fait croître le coefficient de transfert de chaleur entre le fluide en mouvement et les particules sphériques.

La stratification thermique initiale du lit durant la période de stockage est illustrée sur la Figure 3, qui montre la distribution horaire de la température le long du lit durant le stockage de la chaleur à l'intérieur du lit. Du fait d'un coefficient de transfert de chaleur important entre l'air et le matériau constituant le lit, la zone du lit située près de l'entrée est chauffée, alors que la température de la partie du lit située près de la sortie demeure inchangée et la température de l'air à la sortie reste par conséquent très proche de la température initiale du lit. Au fur et à mesure que le temps augmente, on remarque que le front de température se déplace vers les couches les plus profondes du lit pour atteindre, dans le cas choisi, après trois heures de fonctionnement la saturation thermique et le lit est dit alors chargé.
<figure><img src="image_3.png" /><figcaption>Fig. 3: Distribution de la température dans le lit pour une température de l'air à l'entrée constante, durant le stockage.</figcaption></figure>

Pour ces mêmes raisons, nous constatons sur les Figures 4.(a) et 4 (b) où sont représentées l'évolution de la température du MCP en fonction du temps pour deux positions axiales et un débit d'air constant, dans le cas d'un processus de changement de phase respectivement isotherme et non isotherme, que plus on s'éloigne de l'entrée du lit plus le palier correspondant à la fusion du matériau est grand.

<figure><img src="image_4.png" /><figcaption>Fig. 4: Evolution de la température du MCP pour deux positions axiales, en fonction du temps, durant le stockage.</figcaption></figure>

Nous avons représenté sur la Figure 5 l'évolution de la température du MCP en fonction du temps et de la profondeur du lit durant la période de déstockage. Cette figure montre que plus la profondeur augmente plus la température du matériau met du temps pour atteindre la température d'entrée de l'air. Cette température tend par ailleurs à s'uniformiser puis à se stratifier.

La Figure 6 représente l'évolution de la température moyenne du lit en fonction du temps,
au cours d'un cycle (stockage suivi du déstockage), obtenue pour une vitesse d'écoulement et
une température d'entrée d'air fixes.

Des Figures comparables ont été obtenues par Adebiyi [3].

<figure><img src="image_5.png" /><figcaption>Fig. 5: Evolution de la température du MCP pour une distance x, en fonction du temps, durant le déstockage.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Fig. 6: Evolution de la température moyenne du lit en fonction du temps au cours d'un cycle (stockage et déstockage).</figcaption></figure>

## 5. CONCLUSION

Cette étude a porté sur une analyse du stockage de l'énergie thermique par la chaleur latente à basses températures dans un lit fixe rempli de sphères uniformes contenant chacune un matériau à changement de phase traversé par un flux d'air.

Un modèle théorique unidimensionnel à deux phases séparées est appliqué pour les deux modes (stockage, déstockage) et résolu numériquement lors d'un changement de phase isotherme et non isotherme. Ce modèle a permis de mettre en évidence l'influence des paramètres dynamiques et thermiques de l'air à l'entrée du lit et de montrer que la vitesse d'écoulement de l'air ainsi que le rapport hauteur/diamètre du lit avaient une influence notable sur le processus de stockage.
De plus il a permis prédire les performances qui caractérisent un tel système de stockage telle que la durée du stockage et l'efficacité du lit.

## NOMENCLATURE

a taux de surface de contact, $6 (1 e)/d_s$, m⁻¹

c<sub>F</sub> chaleur spécifique du fluide, J/kg °C

c<sub>S</sub> chaleur spécifique du MCP solide, J/kg °C

c<sub>l</sub> chaleur spécifique du MCP liquide, J/kg°C

D diamètre du lit, m

d<sub>S</sub> diamètre des sphères, m

h coefficient d'échange thermique, W/m²°C

h<sub>l</sub> chaleur latente de fusion, J/kg

H hauteur du lit, m

k conductivité thermique, W/m °C

T température du fluide, °C

T<sub>a</sub> température ambiante, °C

T<sub>e</sub> température du fluide à l'entrée, °C

T<sub>s</sub> température du fluide à la sortie, °C

u vitesse interstitielle du fluide, m/s

ε porosité

ρ<sub>F</sub> masse volumique du fluide, m³/kg

ρ<sub>S</sub> masse volumique moyenne du MCP, m³/kg

Θ température du MCP, °C

Θ<sub>S</sub> température de saturation du MCP, °C

Φ titre du mélange des deux phases du MCP (0<Φ<1)

## REFERENCES

[1] J. Duffie. and W. Beckmann, "Solar Engineering of Thermal Processes", John Wiley and Sons.(1980) New York.

[2] D. E. Beasley et C. Ramanarayan, "Thermal Response of Packed Bed of Spheres Containing a Phase Change Material", International Journal of Energy Research, 13(1989) 253 - 265

[3] G. A. Adebiyi, "A Second-Law Study on Packed Bed Energy Storage Systems Utilising Phase Change Materials", ASME journal of solar Energy Engineering, Vol.113(1991) pp146 - 156

[4] A. Benmansour et M. A. Hamdan, "Stockage de l'Energie Thermique dans un Matériau à Changement de Phase", CHEMSS'2000. Alger, 13 - 15 mai 2000, pp 58 - 63.
[5] A. Benmansour et M. A. Hamdan, "Analyse du Stockage Energie Thermique par Chaleur Latente", CIMASI'2000 Casablanca, Maroc, 23-25 octobre 2000.

[6] S. V. Patankar, "Numerical Heat Transfer and Fluid Flow", Hemisphere Publishing Corporation. Mc Graw Hill Company 1980.

# ANNEXE

Cette annexe explique l'origine des équations (2)-(4)

Considérons un lit de capsules sphériques contenant chacune un matériau fusible comme un milieu poreux. Lorsqu'il est traversé par un courant d'air chaud à une température supérieure à la température de fusion (charge du lit), il se produit un échange de chaleur entre la phase fluide et la phase solide constituant le lit. A l'entrée du lit il y'a apparition d'un front de fusion qui progresse en son sein dans le sens de l'écoulement de l'air jusqu'à ce que le matériau contenu dans toutes les capsules soit fondu. L'effet inverse se produit en déchargeant le lit par un courant d'air froid.

Les équations traduisant les transferts d'énergie ( conservation de l'enthalpie ) s'écrivent dans une approche macroscopique en tenant compte des hypothèses suivantes:

- le milieu considéré est isotrope,

- se plaçant dans le cas d'une convection forcée dans la phase fluide, nous négligeons l'action de la pesanteur,

- le couplage gradient de concentration-gradient de température est négligé.

- pour la phase fluide:

$$ \rho_F c_F \varepsilon \left( \frac{\partial T}{\partial t} + \bar{u} \cdot \vec{\nabla} T \right) = \varepsilon \vec{\nabla} (k_F \vec{\nabla} T) + h \alpha (\theta - T) \quad (A1) $$

- pour le matériau ( phase solide ):

$$ \rho_S c_S (1 - \varepsilon) \frac{\partial \theta}{\partial t} = (1 - \varepsilon) \vec{\nabla} (k_S \vec{\nabla} \theta) + h \alpha (T - \theta) \quad (A2) $$

En considérant le problème monodimensionnel et la conduction entre les particules sphériques négligeable, les d'équations (A1) et (A2) se réduisent à un système de deux équations simplifiées:

- pour la phase fluide:

$$ \rho_F c_F \varepsilon \left( \frac{\partial T}{\partial t} + u \frac{\partial T}{\partial x} \right) = \varepsilon \left( k_F \frac{\partial^2 T}{\partial x^2} \right) + h \alpha (\theta - T) \quad (A1') $$

- pour le matériau:

- Phase solide:

$$ \rho_S c_S (1 - \varepsilon) \frac{\partial \theta}{\partial t} = h \alpha (T - \theta) \quad (A2') $$

- Phase liquide:

On suppose que les propriétés physiques du matériau fusible à l'état liquide différent de leurs valeurs correspondant à l'état solide. L'équation (A2') devient:

$$ \rho_S c_L (1 - \varepsilon) \frac{\partial \theta}{\partial t} = h \alpha (T - \theta) \quad (A2'') $$

- Deux phases:
On détermine la position du front de fusion en exprimant en cet endroit que la chaleur cédée par le fluide est égale à la chaleur utilisée pour la fusion ( la chaleur transmise par conduction dans le matériau fusible n'étant pas prise en compte). Si $h_l$ est la chaleur latente massique de fusion, $\Theta_s$ la température de fusion et $m_l$ la masse fondue par unité de volume du lit, il vient:

$$ \frac{\partial(m_l h_l)}{\partial t} = h a (T - \Theta_s) \qquad (A3) $$

En introduisant la fraction massique fondue (le titre) $\Phi = \frac{m_l}{m}$, cette équation s'écrit:

$$ \frac{\partial \Phi}{\partial t} = \frac{h a}{m h_l} (T - \Theta_s) \qquad (A3') $$

où m est la masse de matériau fusible contenue dans l'unité de volume du lit.

En introduisant la fraction volumique de phase $\varepsilon = \frac{V_a}{V}$

et par suite:

$$ \frac{V_s}{V} = \frac{V - V_a}{V} = (1 - \varepsilon) $$

où V, $V_a$ et $V_s$ sont respectivement le volume du lit, du fluide et du matériau fusible, l'équation (A3') devient:

$$ \rho_s h_l (1 - \varepsilon) \frac{\partial \Phi}{\partial t} = h a (T - \Theta_s) \qquad (A3'') $$

Les distributions des températures dans le fluide et le lit s'obtiennent alors en résolvant les équations (A1'), (A2'), (A2'') et (A3'').