# Hybrid analytical and finite element simulation of a latent heat storage exchanger: Comparison of heat transfer models for enhanced thermal conductivity phase change material

Mariam Jadal, Didier Delaunay, Lingai Luo, Jérôme Soto, Nicolas Boyard

► To cite this version:

Mariam Jadal, Didier Delaunay, Lingai Luo, Jérôme Soto, Nicolas Boyard. Hybrid analytical and finite element simulation of a latent heat storage exchanger: Comparison of heat transfer models for enhanced thermal conductivity phase change material. 16th International Heat Transfer Conference, IHTC 2018, Aug 2018, Pékin, China. 10.1615/IHTC16.ecs.023118. hal-02376167

HAL Id: hal-02376167

https://hal.science/hal-02376167

Submitted on 14 Jun 2021

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
# Hybrid analytical and finite element formulation of a latent heat storage exchanger simulation: Comparison of heat transfer models for enhanced thermal conductivity phase change material

Mariam JADAL¹, ²*, Didier DELAUNAY¹, Lingai LUO¹, Jérôme SOTO², Nicolas BOYARD¹

¹Laboratoire de Thermique et Energie de Nantes, Polytech’Nantes – UMR 6607
La Chantrerie, rue Christian Pauc – 44306 Nantes Cedex 3

²Institut Catholique d’Arts et Métiers de Nantes
35 avenue du Champ de Manoeuvres, 44470 Carquefou

## ABSTRACT

In order to cut down the computing time, a simplified method for heat transfer analysis within a multi-tube and multi-plate latent heat storage exchanger is proposed. This method associates a 2D Finite Element method to solve energy equation within each plate and an analytical model to link the plates between them. Different formulations of the energy equation have been tested and compared with experimental results on a simple configuration. Enthalpy formulation, temperature dependent specific heat method, and source term of crystallization associated to kinetic equation method are the different options studied. The limitations of each methods is then exhibited.

**KEY WORDS:** Energy efficiency, Heat exchanger, Thermal storage, Heat transfer enhancement.

## 1. INTRODUCTION

In order to reduce industry waste heat, the energy has to be reused. For sufficiently high temperature heat, the energy could be converted into electricity. For low temperature heat, the energy is mainly used for heating the premises of a factory. Alternatively, the waste heat could be used within an industrial process where the energy could be stored during the cooling phase and reused during the heating phase.

Latent heat storage is a clever way to store energy. Compared to sensible storage, where the energy is stored by temperature rising, the latent heat storage stores the energy at a nearly constant temperature during the phase change of the material. This property allows to store a large amount of energy and to obtain energy of better quality. In addition, the latent heat storage provides a better energy density.

The materials used for this type of storage are called "Phase Change Materials (PCM)". The use of phase-change materials as storage materials has the disadvantage of low thermal conductivity. It turns out that most PCMs generally has a low thermal conductivity (0.1 - 0.4 W·m⁻¹·K⁻¹), which limits the thermal exchange between the heat transfer fluid and the storage environment.

Therefore, the design of an efficient storage system requires the development of thermal transfer improvement techniques. There are many intensification methods in the literature, such as improving thermal conductivity, or increasing the exchange surface with PCM [1] [2] [3]. Based on numerous previous studies [4], a concept of storage by latent heat with enhanced thermal conductivity has been developed. The storage is based on a technology incorporating a conductive structure of Expanded Natural Graphite (ENG),

*mariam.jadal@univ-nantes.fr
impregnated with a PCM. The ENG has an orthotropic structure, leading to a concept of storage where the
heat transfer fluid crosses orthogonally the stack of plates in which transfers in the plane are preferred. This
system presents the most promising results [5] for short thermal cycle industrial applications.

Based on this concept, a 6kWh – 100kW prototype has been developed and tested for a sterilization
application. The experimental results show a good accordance compared to theoretical prediction. The
second step of the study is the optimization of the heat exchanger topology including the thermal properties
of the Composite Phase Change Material (CPCM). In order to carry out the optimization study, a numerical
tool has to be developed. The complexity of the geometry to optimize leads to large time consuming 3D
calculations inappropriate for a parametric study. Therefore, a method has been developed where a 2D Finite
Element method is associated with an analytical model for the third direction. Also, different formulations of
the energy conservation equation have been chosen: enthalpy formulation, temperature dependent specific
heat method and crystallization kinetic method. The last method was done in order to take into account the
kinetic of the crystallization during the cooling phase. This choice is interesting if the PCM has a
crystallization kinetic that depends on the thermal history (cooling rate). The difference between the three
methods will be discussed later on. All the numerical results have been confronted to the experimental
results obtained on a simple configuration [1]. The limitations of each methods is then exhibited.

## 2. Numerical and experimental study of the latent heat storage exchanger

### 2.1 Description of the experimental device

The storage module studied is a multi-tube heat exchanger where five copper tubes with a diameter of 12mm pass through 20 plates of composite material measuring 100x100x20mm (Fig. 1). To anticipate thermal expansion of the composite material, plates are spaced by a few millimeters in the thickness direction of the plates. The composite material is composed of 25% ENG and 75% of RT70 paraffin from RUBITHERM®. The heat exchanger storage is tested in a hydraulic bench where the fluid temperature and flow rate are imposed with a thermo-regulator and a cooling system (Fig. 2).

<figure><img src="image_2.png" /><figcaption>Fig. 1 Picture of the heat exchanger.</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>Fig. 2 Experimental set-up diagram.</figcaption></figure>

### 2.2 Instrumentation

Temperature distribution within the composite material has been obtained using 10 K-type thermocouples positioned in different composite material plates according to the fluid flow (Fig. 3). These thermocouples, in location 2 (Fig. 4), have been used to study the temperature gradient according to the height of the heat exchanger. PT100 probes were used in the exchanger to evaluate the inlet and outlet temperatures of the heat transfer fluid. The flow rate has been measured by a Coriolis type flowmeter.

<figure><img src="image_4.png" /><figcaption>Fig. 3 Position of thermocouples in the plates, depending on the height of storage.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Fig. 4 Composite material plate crossed by five tubes.</figcaption></figure>
## 2.3 Description of the numerical modeling method

The energy equation is solved in a x-y 2D plane using Finite Element methods under Comsol Multiphysics®. Each plate is assumed to be isothermal in the z-direction, i.e. there is no temperature gradient within the plate thickness. And, the plates are associated with an analytical model to couple the plates in the heat transfer fluid direction under Matlab®.

### Boundaries conditions

For symmetry reasons, each plate geometry has been modelised with one-quarter of the plate with five tubes (Fig. 5). Then, all the plate walls are assumed to be adiabatic ($\varphi = 0$) because of symmetry or efficient insolation, at the exception of the surfaces of the tubes where a convection conditions is imposed. At this boundary, the average temperature of the heat transfer fluid $T_{average}(k)$ is associated with the equivalent convective exchange coefficient $h_{eq}$. The coefficient $h_{eq}$ includes the wall thermal resistance, the thermal contact resistance between the composite material and the tube identified previously [8], and the convective exchange coefficient of heat transfer fluid $h_{fluid}$ corresponding to the plate number k. The heat transfer fluid exchange coefficient has been calculated from the Colburn correlation [9].

<figure><img src="image_6.png" /><figcaption>**Fig. 5** 2D model in the plane of a composite material plate.</figcaption></figure>

### Thermal properties of the composite material

A DSC analysis has been performed at different cooling and heating rates on the RT70 paraffin from Rubitherm®. The evolution of specific heat according to temperature has been calculated and implemented in Comsol as polynomial fit by temperature intervals. This specific heat corresponds to the apparent heat capacity. The melting range of RT70 paraffin is 66 to 74°C (Fig. 6). The solidification range is between 58 and 69°C (Fig. 7). There are two exothermic peaks. The first peak represents the liquid/solid phase change, while the second peak represents the solid/solid transition. This phenomenon is known and observed for this type of PCM [11].

<figure><img src="image_7.png" /><figcaption>**Fig. 6** Specific heat of RT70 paraffin according to temperature (melting phase).</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Fig. 7** Specific heat of RT70 paraffin according to temperature (solidification phase).</figcaption></figure>

The specific heat of ENG obtained by DSC, has been implemented into the Comsol model by the trend curve equation. A mixing law has been used to calculate the specific heat of the composite material PCM-ENG.
The DSC flux integration allows to obtain the total enthalpy variation of the material in the melting and solidification phase (Fig. 8). This paraffin has a latent heat of 270 KJ/Kg. There is no difference in enthalpy variation between the melting and solidification phases, indicating good thermal stability of the material.

<figure><img src="image_9.png" /><figcaption>**Fig. 8** RT70 paraffin Enthalpy according to temperature.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Fig. 9** Evolution of solidified fraction according to temperature for the three cooling speeds.</figcaption></figure>

A kinetic study was performed on the sample by imposing three heating and cooling rates (5, 10, 30 K.min⁻¹). Based on the heat flow results according to the temperature obtained by DSC, the transformed fraction is calculated on the liquid/solid transition peak. The choice of the kinetic study on a single peak has been made for simplicity reasons, in order to test, initially, the feasibility of the 2D finite element simulations methods, associated with a coupling with Matlab. Nevertheless, using the crystallization kinetics method in Comsol, the enthalpy value takes into account the two solidification peaks of the RT70 paraffin. A deconvolution of the peaks has been carried out. The solidified fraction has been classically identified by the ratio between the partial area under the peak between the beginning of the solidification and the current temperature T, and the total area under the peak (Fig. 9) [7]. This variable without unit ranges from 0 (solid sample) to 1 (liquid sample). We note on Fig. 9 that as the cooling rate increases, the initial crystallization temperature shifts to lower temperatures.

Crystallization of RT70 paraffin at constant cooling rate has been exploited using Ozawa theory [7]. The study has been carried out between 65 and 72° C. This choice has been done in order to take into account the evolution of the solidified fraction according to the three cooling rates for each temperature. Ozawa theory gives the evolution of the solidified fraction α(t) by equation (1) where V is the cooling rate and n a coefficient. The coefficients n and K<sub>Ozawa</sub>(T) have been identified by plotting ln(-ln(1-α)) according to ln(V). For each temperature, we obtain a portion of straight slope (-n) and intercept ln(K<sub>Ozawa</sub>(T)).

$$ \alpha(t) = 1 - \exp\left(\frac{K_{oz}(T)}{Vn}\right) \qquad (1) $$

We have found that n=1.3. And a validation of the K<sub>Ozawa</sub> values obtained has been carried out by comparing the solidified fraction obtained using the DSC and the solidified fraction calculated, using the K<sub>Ozawa</sub> values, with the formula (1).

Fig. 10 shows a comparison between the calculated and measured solidified fraction. The curves show a good agreement.

<figure><img src="image_11.png" /><figcaption>**Fig. 10** Comparison between the calculated and measured solidified fraction.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Fig. 11** Evolution of K<sub>Nakamura</sub> according to temperature.</figcaption></figure>
Finally, the curve of K_Nakamura coefficient according to the temperature [7] has been calculated from the formula:

$$K_{Nak} = \left( - \frac{d(K_{oz}^{n}(T))}{dT} \right)^{\frac{1}{n}} \qquad (2)$$

Fig. 11 shows the evolution of K_Nakamura according to temperature. These values has been implemented in Comsol to calculate the source term in the energy equation:

$$\frac{d\alpha}{dt} = n \times K_{Nak}(T) \times (1 - \alpha) \times \left[ \ln\left(\frac{1}{1 - \alpha}\right) \right]^{\frac{n-1}{n}} \qquad (3)$$

The value of plane thermal conductivity has been calculated by using an experimental device developed in our laboratory [10] while the transverse conductivity has been measuring using classical guarded hot plate test apparatus. The values of plane and transverse thermal conductivity of the composite material, are respectively λ-plane = 19.8 W. m⁻¹. K⁻¹ and λ-transverse = 6.12 W. m⁻¹. K⁻¹. The density value of the composite material has been measured to ρ= 878 Kg. m⁻³[4].

## Numerical modeling by finite element method

The energy equation with different phase change models is solved within the plate. The models are the following:

**Temperature dependent specific heat method:** This method consists in taking into account the amount of heat absorbed or released during the phase change in the specific heat coefficient. The evolution of the Cp during the transformation obtained directly using the DSC measurement is called “Apparent specific heat”.

**Enthalpy method:** The enthalpy method solves the heat equation in terms of enthalpy variation where the enthalpy of the phase change material has been calculated from the values of Cp.

**Crystallization kinetics method:** The crystallization kinetics model has been developed in our LTeN [6] laboratory. This numerical model is used to simulate the crystallization of phase change material with a crystallization kinetic. The numerical model implemented in Comsol allows to take into account the evolution of the Nakamura function according to temperature [7].

## Calculation methodology

The resolution algorithm is described in the following diagram. It illustrates the coupling between finite element 2D simulation in Comsol and the link performed by the routine in Matlab (Fig. 12).

<figure><img src="image_13.png" /><figcaption>Fig. 12 Algorithm illustrates the 2D simulation in Comsol coupled with an analytical model in Matlab.</figcaption></figure>
where J is the iteration number, k is the number of plates, T<sub>average</sub> (t) is the average temperature of heat transfer fluid on the plate thickness, T<sub>in</sub> (t) and T<sub>out</sub> (t) are respectively the inlet temperature and the outlet temperature of the heat transfer fluid, h<sub>eq</sub> is the equivalent exchange coefficient, φ is the heat flux, ḍ<sub>fluid</sub> (t) is the measured mass flow rate, C<sub>p</sub> (t) is specific heat. The thermal properties of the heat transfer fluid were calculated according to the experimental temperatures at each moment.

A quadratic criterion has been set for the average temperature of the heat transfer fluid in Matlab, by: (T<sub>average</sub>(j+1)- T<sub>average</sub>(j)) 2 >10<sup>-4</sup>. This criterion means that as long as the average temperature variation between two iterations is greater than 1x10<sup>-4</sup>, the software will continue to iterate. Convergence was achieved each time after four iterations.

A validation of energy conservation in the numerical model has been carried out. For that, a comparison between the enthalpy variation between the beginning and the end of the storage, and the energy absorbed or released by the heat transfer fluid in all the storage module has been done. This comparison has been performed using the three numerical methods. The differences obtained vary between 0.1 to 8%.

## 2.4 Numerical and experimental results comparison

### Experimental temperature variation within the storage module

Fig. 13 and 14 show the temperature evolution of the composite material in the fifteen plate at thermocouple 2 location. The initial composite material temperature at the beginning of the cycle was maintained at 59°C. Then a temperature step of 15K was imposed at the inlet temperature of the heat transfer fluid. At the early times, a particular behavior can be observed with a large difference between the inlet and outlet temperatures of the heat transfer fluid in the heat storage exchanger. This can be explained by the transitional period of water evacuation remaining in the pipes following the tests carried out previously. This is not taken into account in the simulation.

<figure><img src="image_14.png" /><figcaption>**Fig. 13** Temperature evolution in the fifteen plate and at the inlet and outlet during melting phase.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>**Fig. 14** Temperature evolution in the fifteen plate and at the inlet and outlet during solidification phase.</figcaption></figure>

### Comparison between numerical and experimental heat flow

The evolution of the heat flow delivered by the heat transfer fluid during the heating and cooling phase is obtained by the following formula:

$$ \Phi(t)_{\text{experimental}} = \dot{m}_{\text{fluid}}(t) \times C_{p_{\text{fluid}}}(t) \times \Delta T_{\substack{\text{inlet} \\ \text{outlet}}} (t) \quad (4) $$

where ΔT<sub>inlet/outlet</sub> is the temperature difference of the heat transfer fluid between the inlet and outlet of the heat exchanger. C<sub>pfluid</sub> (t) and ḍ<sub>fluid</sub> (t) are respectively the specific heat and mass flow rate of the fluid. The thermal properties of the heat transfer fluid (water) were used to calculate these parameters at all times.
The convective heat flow of the heat transfer fluid passing through the composite material has been numerically evaluated in Comsol by the following formula:

$$ \Phi(t)_{\text{Numerical}} = \frac{(T_{\text{heat transfer fluid}}(t) - T_{\text{material surface}}(t))}{(R_{tc} + \frac{1}{h_{\text{fluid}}(t)})} \quad (5) $$

where $h_{\text{fluid}}$ is the convective heat transfer coefficient in the heat transfer fluid, calculated by the Colburn correlation [9]. $R_{tc}$ is the thermal contact resistance between copper tube and composite material fixed to $5.10^{-4}$ m².K⁻¹.W⁻¹[8]. $T_{\text{heat transfer fluid}}$ is the average temperature of the heat transfer fluid calculated numerically in Matlab. $T_{\text{material surface}}$ is the surface temperature of the composite material calculated with Comsol.

Fig. 15 shows a comparison between the heat flow measured experimentally and the heat flow calculated numerically using apparent $C_p$ method and enthalpy method, during the melting phase. We can notice that the curves show a behavior similar and give values of the same order of magnitude. These values vary between 500 and 1800W throughout the entire storage module, according to the time. There is a higher heat flow at the beginning of the cycle. This is due to the highest temperature difference at the beginning of the cycle between the heat transfer fluid and the composite material. Experimental heat flow oscillations are due to regulation.

Fig. 16 shows a comparison between the heat flow measured experimentally and the heat flow calculated numerically using apparent $C_p$ method, enthalpy method, and source term of crystallization kinetic method during the solidification phase. The apparent $C_p$ method shows a convergence problem related to the presence of the two peaks of the apparent $C_p$ according to the temperature where the phase change takes place. Nevertheless, using the enthalpy method, we can observe that the evolution of numerical and experimental heat flow show a good fit, because enthalpy takes into account the phase change without presenting as strong non-linearity. Comparing between the enthalpy method and the crystallization kinetics method, we can notice that the curves are very close and present a good agreement with the experimental result.

<figure><img src="image_16.png" /><figcaption>**Fig. 15** Heat flow measured experimentally and calculated numerically (melting phase).</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>**Fig. 16** Heat flow measured experimentally and calculated numerically (solidification phase).</figcaption></figure>

Temperature comparison in the stock

Fig. 17 shows the temperature variation of thermocouple 2 in the fifteenth plate from the experimental values, compared to the temperature evolution obtained numerically using apparent $C_p$ method and enthalpy method during the melting phase. The curves obtained show a good fit with a slight difference between the experimental and numerical methods, especially for the apparent $C_p$ method. This can be explained by the difference of energy conservation achieved previously for each method studied.

Fig. 18 shows the evolution of the temperature in the fifteenth plate obtained experimentally and numerically using enthalpy method and crystallization kinetic method during the solidification phase. The agreement is good.
<figure><img src="image_18.png" /><figcaption>**Fig. 17** Evolution of temperature in the fifteenth plate obtained experimentally and numerically (melting phase).</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>**Fig. 18** Evolution of temperature in the fifteenth plate obtained experimentally and numerically (solidification phase).</figcaption></figure>

### 3. CONCLUSION

A numerical modeling of a multi-tubes and multi-plates storage exchanger was performed under Comsol using three different numerical methods. The calculated and measured heat flows were found very close, using the apparent $C_p$ method, enthalpy method during the melting phase. The limits of the apparent $C_p$ method appeared in the case of a non-linearity of the $C_p$, in particular for RT70 paraffin during the solidification phase. However, using enthalpy method, the comparison with the heat flow measured experimentally and calculated numerically showed a good fit with very close value. The results using the enthalpy method and the crystallization kinetic method, were appeared to be very close. In future works we will verify that this result is also true by taking into account the two peaks that appear in the DSC during cooling phase.

### ACKNOWLEDGEMENTS

The authors would acknowledge the European Commission for the financial support from H2020 MSCA programme and ADEME and Pays de la Loire for PhD thesis funding.

### REFERENCES

[1] F. Agyenim, P. Eames, and M. Smyth, "Heat transfer enhancement in medium temperature thermal energy storage system using a multitube heat transfer array," Renew. Energy, vol. 35, pp. 198-207, 2010.

[2] M. Martinelli, "Stockage d'énergie thermique par changement de phase - Application aux Réseaux de Chaleur", Thèse de doctorat en Mécanique des fluides, Energétique et Procédés, Université de Grenoble Alpes, 2016.

[3] A. Soupart, P. Marty, J. Fourmigué, and R. Couturier, "Le stockage de chaleur dans les Matériaux à Changement de Phase : étude expérimentale des transferts thermiques dans les échangeurs de stockage latent," in Congrès Français de Thermique, 2014.

[4] K. Merlin, "Caractérisation thermique d'un matériau à changement de phase dans une structure conductrice", Thèse de doctorat en Energétique, Thermique, Combustion, Université de Nantes, 2016.

[5] K. Merlin, D. Delaunay, J. Soto, L. Traonvouez, "Heat transfer enhancement in latent heat thermal storage systems: comparative study of different solutions and thermal contact investigation between the exchanger and the PCM," Applied Energy, Vol 166, pp.107-116, 2016.

[6] A. Levy, "A Novel Physics Node for Nakamura Crystallization Kinetics," International Journal of Theoretical and Applied Mathematics, 2017.

[7] E. KOSCHER, Thèse de doctorat en Matériaux polymeres et composites, Université CLAUDE BERNARD, Lyon 1, 2002.

[8] P. Garcia, M. Olcese, and S. Rougé, "Experimental and Numerical Investigation of a Pilot Scale Latent Heat Thermal Energy Storage for CSP Power Plant," Energy Procedia, vol. 69, pp. 842-849, May 2015.

[9] M. Martinelli, F. Bentivoglio, A. Caron-Soupart, R. Couturier, J. Fourmigue, P. Marty, "Experimental study of a phase change thermal energy storage with copper foam", Applied Thermal Engineering, vol 101, pp. 247-261, 2016.

[10] M. Jadal, K. Merlin, D. Delaunay, L. LUO, J. Soto, "Etude d'un dispositif de caractérisation thermique de matériaux composites : application à un matériau à changement de phase dans une structure conductrice", Compte rendu du Congrès de la Société Française de la Thermique SFT, Marseille 2017.

[11] C. Vélez, M. Khayet, J.M. Ortiz de Zárate, "Temperature-dependent thermal properties of solid/liquid phase change even-numbered n-alkanes: n-Hexadecane, n-octadecane and n-eicosane", Applied Energy, Vol.143, pp. 383-394, 2015.