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

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.


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 phasechange 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 -1 .K -1 ), 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), 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.

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).

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.

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 (φ =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 Taverage(k) is associated with the equivalent convective exchange coefficient heq.The coefficient heq 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 hfluid corresponding to the plate number k.The heat transfer fluid exchange coefficient has been calculated from the Colburn correlation [9].
Fig. 5 2D model in the plane of a composite material plate.

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].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.
Fig. 9 Evolution of solidified fraction according to temperature for the three cooling speeds.
A kinetic study was performed on the sample by imposing three heating and cooling rates (5, 10, 30 K.min -1 ).
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 KOzawa(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(KOzawa(T)). ( We have found that n=1.3.And a validation of the KOzawa values obtained has been carried out by comparing the solidified fraction obtained using the DSC and the solidified fraction calculated, using the KOzawa values, with the formula (1).Fig. 10 shows a comparison between the calculated and measured solidified fraction.The curves show a good agreement.
Fig. 10 Comparison between the calculated and measured solidified fraction.Finally, the curve of KNakamura coefficient according to the temperature [7] has been calculated from the formula: Fig. 11 shows the evolution of KNakamura according to temperature.These values has been implemented in Comsol to calculate the source term in the energy equation: 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.8W. m -1 .K -1 and λtransverse = 6.12 W. m -1 .K -1 .The density value of the composite material has been measured to ρ= 878 Kg. m -3 [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).

Fig. 12
Algorithm illustrates the 2D simulation in Comsol coupled with an analytical model in Matlab. ( where J is the iteration number, k is the number of plates, Taverage (t) is the average temperature of heat transfer fluid on the plate thickness, Tin (t) and Tout (t) are respectively the inlet temperature and the outlet temperature of the heat transfer fluid, heq is the equivalent exchange coefficient, ϕ is the heat flux, ṁfluid (t) is the measured mass flow rate, Cp (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: (Taverage(j+1)-Taverage(j)) ² >10 -4 .This criterion means that as long as the average temperature variation between two iterations is greater than 1x10 -4 , 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%.

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.

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: (4) where ∆T inlet/outlet is the temperature difference of the heat transfer fluid between the inlet and outlet of the heat exchanger.Cpfluid (t) and ṁfluid (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: (5) where hfluid is the convective heat transfer coefficient in the heat transfer fluid, calculated by the Colburn correlation [9].Rtc is the thermal contact resistance between copper tube and composite material fixed to 5.10 -4 m 2 .K -1 .W -1 [8].Theat transfer fluid is the average temperature of the heat transfer fluid calculated numerically in Matlab.Tmaterial 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 Cp 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 Cp method, enthalpy method, and source term of crystallization kinetic method during the solidification phase.The apparent Cp method shows a convergence problem related to the presence of the two peaks of the apparent Cp 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.

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 Cp 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 Cp 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.

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 Cp method, enthalpy method during the melting phase.The limits of the apparent Cp method appeared in the case of a non-linearity of the Cp, 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.

Fig. 3
Fig. 3 Position of thermocouples in the plates, depending on the height of storage.

Fig. 4
Fig. 4 Composite material plate crossed by five tubes.

Fig. 7
Fig. 7 Specific heat of RT70 paraffin according to temperature (solidification phase).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.

Fig. 13
Fig. 13 Temperature evolution in the fifteen plate and at the inlet and outlet during melting phase.

Fig. 17
Fig. 17 Evolution of temperature in the fifteenth plate obtained experimentally and numerically (melting phase).

Fig. 18
Fig. 18 Evolution of temperature in the fifteenth plate obtained experimentally and numerically (solidification phase).