Imed Adji Diop

Ben Talouba
email: imed.bentalouba@univ-rouen.fr

Laurent Balland

Nordine Mouhab

Diop Adji

Ben Talouba Imed

Balland Laurent

Mouhab Nordine

Thermal characterization of a biodiesel nitration: Bio-additive's synthesis by calorimetric methods

Keywords: methyl oleate, cetane number, reaction calorimeter, DSC, MTSR

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.

Introduction

The industrialization and motorization in recent years come along with a high demand for petroleum products and an increase in air pollution [1]. It has been established by the Interprofessional Center for Technical Studies of Atmospheric Pollution (CITEPA), in France, that the automotive sector (mainly diesel engines) plays a leading role in air pollution (CO2, CO, NOX, SO2).

To reduce the pollution related to diesel engines, several solutions have been tested including the setting up of catalytic oxidation pots, particles filters [START_REF] Courtils | La dépollution des gaz d'échappement diesel au moyen de pots catalytiques d'oxydation[END_REF] and the use of additives, mainly cetane number improvers [START_REF] Song | Chemistry of Diesel Fuels[END_REF].

Cetane improver additives are commonly used to increase engine performances. According to Beatrice et al 1996 [4], an increase in the cetane number reduces NOx and particle emissions.

There are a lot of cetane number improvers [START_REF] Suppes | Cetane-Improver Analysis and Impact of Activation Energy on the Relative Performance of 2-Ethylhexyl Nitrate and Tetraethylene Glycol Dinitrate[END_REF] [6] but the 2-EthylHexyl Nitrate (ENH) is the main additive used for diesel. The EHN is an energy molecule obtained by nitration of isooctanol (fossil fuels) with mixed nitric/sulfuric acid (1:1.5 w/w).

Several studies showed that the iso-octanol nitration is a strongly exothermic reaction which takes place in a heterogeneous liquid-liquid mixture [START_REF] Quadros | Benzene nitration: validation of heterogeneous reaction models[END_REF] [START_REF] Chen | Thermal hazard evaluation for iso-octanol nitration with mixed acid[END_REF]. EHN additive is thermally unstable and sensitive to heat. Its thermal instability is mainly due to energetic functional groups like nitro (C-NO2) and nitrate (C-O-NO2). Those functional groups are composed of easily broken C-N covalent bonds. In fact, according to Yang et al 2014 [START_REF] Yang | Thermal stability of 2-ethylhexyl nitrate with acid[END_REF], several industrial accidents in China are due to the loss of control during the EHN synthesis. The thermal behavior of this reaction need to be carefully investigated to assess the associated risks [START_REF] Wang | Research on the thermal history's influence on the thermal stability of EHN and NC[END_REF].

With the scarcity of fossil fuels, it is important to reduce the dependency on petroleum products. The goal of the work described in this article is to synthesize a substitute to the EHN by nitrating a trans-esterified vegetable oil (biodiesel). Many studies on the use of biodiesel as biofuels [START_REF] Schmidt | BIODIESEL Cultivating Alternative Fuels[END_REF] [START_REF] Johnson | Biodiesel: Awareness, Use and Perceptions of Students at Four U. S. Universities[END_REF] or as bio-additive [START_REF] Canoira | Nitration of biodiesel of waste oil: Nitrated biodiesel as a cetane number enhancer[END_REF] have been realized over the years. The biodiesel's particularity is that it contains less carbon, less sulfide, less water and more oxygen atoms than fossil diesel. Therefore, it is less dangerous for humans [START_REF] Madden | Comparative toxicity and mutagenicity of soy-biodiesel and petroleumdiesel emissions: overview of studies from the U.S. EPA, Research Triangle Park, NC[END_REF] and its amount of oxygen atoms improves combustion [START_REF] Riadi | Synthesis of biodiesel from waste cooking oil by two steps process transesterification and ozonation[END_REF].

In this article, the biodiesel obtained from canola oil was used. According to Abdullah et al 2010 [START_REF] Abdullah | PRODUCTION OF CETANE IMPROVER FROM Jathropa curcas OIL[END_REF], the presence of double bonds favors the nitration reaction. And the canola oil has a lot of double bounds compared to the other vegetable oils. The classic mixed sulfuric/nitric acid [START_REF] Topchiev | Nitration of Hydrocarbons and Other Organic Compounds -1st Edition[END_REF] (1.5:1 w/w) was used as the nitrating agent.

First, a study about the identification and the characterization of the chemical compounds that results from nitration by various analytical techniques (infrared spectrometer, GCMS, electro spray ionization mass spectrometer, and iodine index measurement) was conducted. Then, a complete study on the thermal stability of the reaction mixture (reagents and products) was realized. This study uses the thermal information given by a reaction calorimeter RC1 and a DSC type microcalorimeter. The quality of the synthesized bio-additive was determined by measuring its cetane number according to the ASTM standard D 613.

Materials and products

Products

The biodiesel (Diester®) used was obtained by trans-esterification of canola oil and was provided by Saipol company (Grand-Quevilly-France). Nitric (65%w/w), sulfuric (95% w/w) acids used as nitrating agent and solvents (potassium iodide 1% volumetric, sodium thiosulfate (Na2S2O3) 0.1M, bicarbonate solution for analysis (>99.7%) and Wijs reactive 0.1 mol•L -1 ) are provided by Fisher Chemical. The methyl nonadecanoate (> 98%) for GC analysis was supplied by Sigma Aldrich.

Materials

Calorimetric reactor RC1-RTCal

The nitration was carried out in a calorimetric reactor RC1-RTCal (Mettler®) operating in semi-batch mode. It is equipped with an AP00-0.5-RTCal reactor vessel with 0.5L capacity.

Mixing system is made of an impeller with one baffle. The experimental device is presented in Fig. 1. The reactor has a temperature control system (between -20°C and +165°C) and an RTCal measurement technology based on thermal flow sensors located on the reactor wall.

The temperature regulation includes a thermostatic group which takes care of the various regulation modes thanks to a fluid (Rodorsil 47 V20). This device provides real-time monitoring of the heat released by the nitration reaction which is highly exothermic. The reactive and products stability study was based on thermograms analysis obtained with the DSC Q20. The DSC was provided by TA Instruments with purge using extra pure nitrogen (50 mL•min -1 ). It can work in isothermal and non-isothermal conditions. DSC measures temperature with accuracy of ± 0.05 K. Its specific sensitivity is 1mW•kg -1 [START_REF] Talouba | Kinetic and safety parameters of decomposition of neat Tert-Butyl (2-Ethylhexyl) monoperoxy Carbonate and in organic solvents[END_REF].The device detects temperature changes between the sample and the reference to measure heat flow over time and temperature. Gold-plated high pressure crucibles (M20 crucibles from Swiss Institute for the Promotion of Safety and Security) were used for measurement [START_REF] Shimkin | Optimization of DSC calibration procedure[END_REF].

Chemical analysis materials

The FTIR IR Affinity-1S from Shimadzu (400 to 4000 cm -1 ), GC-MS from Agilent technology 7890A GC/5975C VL MSD, and the mass spectrometry with electro spray ionization from Waters LCT Premier were used to identify and quantify chemical species. GC-MS analyses were performed in a GC6850 gas chromatography coupled with a 5975C MSD Quad Instrument with electrons impact ionization at 70 eV. A HP5-MS 5% Phenyl Methyl Siloxane capillary column (Agilent 19091S-433E) with a 0.25 mm i.d., 30 m nominal length and 0.25 µm phase film diameter was used. The carrier gas was helium at a flow rate of 15 mL•min -1 .

Experiment

Considering the chemical instability of the mixed acid (nitrating agent), it was prepared just before the synthesis reaction (nitration). 86.4 g of sulfuric acid was introduced with a flow rate of 4g•min -1 into the reaction calorimeter filled with 57.6 g of nitric acid. The mixture was stirred at 200 rpm and kept at 5°C during the addition. At the end, the mixture, composed of HNO3 27 %, H2SO4 58%, H2O 15% (w/w), was kept in the fridge at the same temperature.

To carry out the synthesis reaction, 144g of mixed acid was dosed with a flow rate of 4g•min -1 into the reactor initially filled with 150 g of biodiesel composed of 92% of methyl oleate (determined by GC-MS Table 1). These operating conditions are chosen after preliminary study on the effect of parameters such as a flow rate, stirring speed on thermal behavior in order to ensure the safe operation of the reactor. The synthesis reaction was carried out at 600 rpm stirring speed, in the isothermal mode. This experiment was realized at seven different temperatures (10-15-20-25-30-40 and 50°C). The Mettler RC1 operating software allows to record the thermal signal (power profile) and to monitor the reaction progress. The final product was then transferred to a separating funnel and separated in two phases. Only the organic phase (which contains the nitrated products) was analyzed. Before the analysis, the organic phase was, first, washed with a bicarbonate solution at 2% w/w, and then with distilled water until a neutral pH was obtained. Methyl stearate 0.5

Nitration Product analysis

Infrared analysis and electro-spray ionization mass spectrometry were used to identify and quantify the different chemical compounds produced during the biodiesel nitration.

Infrared analysis

Fig. 2 shows the FTIR spectra of biodiesel and nitrated products at different temperatures.

These spectra were obtained after defining the baseline corresponding to the air spectrum with 50 scans. An absorbance difference is noted between 800 and 1700 cm -1 . It allows to identify new functional groups present in the nitrated product (Table 2). The value of the nitrated products absorbance depends on the synthesis temperature.

The characteristic ester frequency C=O bonds at 1750 cm -1 , alkane frequency C-H bonds between 2900-3000 cm -1 and alkene frequency C=C stretch at 1650 cm -1 with weak intensity, obviously appear in all spectra. The spectrum of the nitrated product shows the presence of three new peaks at 1363, 1550 and at 1635 cm -1 . They respectively indicate the stretching of symmetric NO2, that of assymetric NO2 and the absorption of ONO2 (fig. 3).

The absorption spectrum at 1363cm -1 (fig. 3) indicating a symmetric vibration of R-NO2 for nitrated products, can also be observed in biodiesel but is less intense. This peak in the biodiesel spectrum corresponds to the vibration of CH2 groups which is supposed to appear at 1370 cm -1 . To complete this IR analysis which only identifies the functional groups, electro spray ionization mass spectrometry was used to identify the synthesized molecules.

Electro spray ionization mass spectrometry

The electro spray ionization mass spectrometry allows to avoid the degradation of the nitrated products during the heating in the chromatography column. The samples (biodiesel and nitro product) were diluted with acetonitrile to a concentration of 20 μg/mL, before being introduced into a capillary tube with a flow rate of 15 μL•min -1 . Peaks of the spectrum in Fig. 4a correspond to methyl oleate, methyl linolenate and the methyl linoleate. These molecules are in a ionic adduct form [M + Na] + at m / z = 319, m / z = 315 and m / z = 317. For the nitrated biodiesel (Fig. 4b), intense peaks at m/z 359, 377 and 422 are observed. These peaks correspond to [M+NH4] + ion. Table 3 shows the different identified compounds. This study reveals that the nitration reaction can be performed on the carbon atoms carrying the double bond as well as those in the beta position. This result contradicts Abdullah et al (2010) [START_REF] Abdullah | PRODUCTION OF CETANE IMPROVER FROM Jathropa curcas OIL[END_REF] work which indicates that during the nitration process the NO2 or NO3 groups are attached only on the unsaturated carbon atoms.

Three compounds were identified A (Methyl-8-nitro oleate), B (Methyl-10-nitrate stearate) and C (Methyl-10-nitrate-9-nitrostearate).

Study of the reaction selectivity

GC-MS and iodine index measurement were used to determine the molar proportion of each compound A, B and C.

Measurement of the conversion rate

The nitration reaction's conversion rate was determined by using a GC-MS with addition of an internal standard (methyl nonadecanoate) according to the 2011 EN14103 standard. Each solution was dissolved in 10 mL of toluene (GC purity ≥99.5%). 2µL of the sample was injected in the column previously heated at 60°C. This column was heated to 200°C with a 10°C/min heating rate then at 5 °C•min -1 to 240°C. Finally, the column was maintained at this temperature during 7 min. The runtime was 30 min. Splitless mode was used. Table 4 shows that the reaction conversion rate (the production of all molecules A, B and C) varies with temperature and is maximal at 40 °C. The selectivity of B and C was determined by using the iodine index method.

Results are presented in Table 4

Determination of the molecules selectivity

The iodine number of oils and fats measures their unsaturation degree. The Wijs method [START_REF] Sahoré | Quelques méthodes d'analyse biochimique de produits alimentaires[END_REF],

commonly used, requires 30 min absorption time. Sodium thiosulfate solution 0.1M was used as a titration reactant. The measurement of the biodiesel iodine index before and after the nitration reaction allowed to calculate the conversion of unsaturations (the loss of the double bonds). This value was used to determine the proportion of molecules B and C according to the equation 3-1.

=

3-1

The proportion of A was obtained by making a difference between the conversion rate of the reaction (Table 4) and the unsaturations conversion rate. Results are gathered in Table 5 which shows that at high synthesis temperature, the selectivity of B and C decreases. The results obtained in this study show that the reaction mechanism of biodiesel nitration (aliphatic chain) is more complex than that of the aromatic compounds. Indeed, three molecules with different selectivity had been identified (Table 5).

During the synthesis at low temperature (below 20°C), the final product is a stable micro emulsion preventing phase separation. This stability may be due to the hydrophilic behavior of the saturated nitrated products (B and C) at low temperatures.

Calorimetric study

Determination of the reaction heat

The integration of the experimental power profiles (equation 4-1) obtained according to the conditions presented in section 3, shows that the heat released by the reaction varies with the temperature (Table 6). This change in the thermal behavior has also been observed by Chen et al 2012 [START_REF] Chen | Thermal hazard evaluation for iso-octanol nitration with mixed acid[END_REF] in their work about the nitration of iso-octanol but no interpretation has been given. Power profiles (qr) are shown in Fig. 5. Sudden variation of the thermal power with time at the beginning of the reaction corresponds to the nitration of the small amount of the products as Methyl linolenate and Methyl linoleate. Such thermal behavior may present a risk of thermal runaway. The total heat released by the reaction (table 6) represents all the thermal effects linked to the products synthesis, to the mixing phenomena and to the possible decomposition in case of unstable compound. In section 3.3.2, a modification of the reaction mechanism according to the synthesis temperature was noticed. In addition to other phenomena mentioned above, the modification of the reaction mechanism is likely to contribute to the change in the total heat released by the reaction.

= !

Heat related to the mixing phenomena and the decomposition of unstable compounds are presented in sections 4.2 and 4.3.

Heat of mixing measurement

To measure the heat of mixing, the power profile obtained without nitration reaction have to be integrated. The operating conditions were kept identical to those in the synthesis reaction (section 3.1) except for the biodiesel reactant. In fact, to avoid the nitration reaction, the biodiesel was replaced by a molecule close to its main molecule (methyl oleate): methyl palmitate.

The absence of the reaction was confirmed by the palmitate infrared spectrum after the reaction (fig. 6: no nitro group was detected). The non-reactivity of the methyl palmitate is due to the absence of unsaturated carbons. The heat of mixing below 40°C could not be measured because the melting point of methyl palmitate is higher than 35°C.

Measurement of the decomposition heat

The measurement of the heat released during the decomposition was done in two steps. Pure products: nitrated products, nitric acid, sulfuric acid and mixed acid (sulfo-nitric acid) were analyzed first (fig. 8), then samples of mixed acid with different percentages (w/w) of nitrated products (fig. 9). These analyses were realized using a DSC in a non-isothermal mode with Gold-plated high-pressure crucibles and a heating rate of 2°C•min -1 . According to Liu et al [START_REF] Liu | Effects of thermal runaway hazard for three organic peroxides conducted by acids and alkalines with DSC, VSP2, and TAM III[END_REF], the advantage of DSC test under non-isothermal condition is its swiftness and high efficiency in a short time.

For pure products (fig. 8) no signal of decomposition is observed in our range of synthesis temperature (<50°C). The only thermal effects are detected beyond synthesis temperatures.

Nitric acid and mixed acid's DSC curves show an endothermic effect at Tonset = 181.3°C which is attributed to the decomposition of the nitric acid [START_REF] Kazakov | Kinetics and mechanism of thermal decomposition of nitric acid in sulfuric acid solutions[END_REF]. This temperature was determined by the tangent method from zero curvature. Regarding to the nitrated product, it exhibits an exothermic peak at Tonset = 150.9°C which correspond to its decomposition. Its thermal instability was mentioned in other works [START_REF] Yang | Thermal stability of 2-ethylhexyl nitrate with acid[END_REF] [24] [START_REF] Chen | Thermal hazard evaluation for iso-octanol nitration with mixed acid[END_REF]. The Tonset of our additive decomposition is closed to that of fossil nitrated product (EHN) [START_REF] Chen | Thermal hazard evaluation for iso-octanol nitration with mixed acid[END_REF]. 7). This enthalpy (∆H) reported to the total sample mass was determined by integrating the area under the exothermic curve. This study allows to demonstrate that the decomposition of the sulfo-nitric mixed acid and nitrated products also contributes to the modification of the total heat released during the synthesis reaction (Fig. 5).

Study of the safety parameters

In the light of the strong exothermicity of the nitration reaction and the mixed acid decomposition in the presence of the nitrated products, the determination of the Maximum Temperature of the Synthesis Reaction (MTSR) is more than necessary to avoid the possible risk of thermal runaway.

The Maximum Temperature of Synthesis Reaction (MSTR) is the temperature reached by the reaction mixture if the energy potential of the main reactions is released under adiabatic conditions [START_REF] Lerena | Assessment of hazards linked to accumulation in semi-batch reactors[END_REF].

This decomposition releases an excess of energy causing a thermal runaway. To evaluate this temperature, it is necessary to determine the ΔTad (the temperature variation in adiabatic mode) according to equation 4-2.

∆# $! = % & '() * + , ../ 4-2
with Cp the heat capacity of reaction mixture obtained with reaction calorimeter, Mr the mass of the reaction mixture and X'r the full conversion of the biodiesel.

The calculation of the MTSR was done in semi-batch mode according to equation 4-3 [START_REF] Zichao | Insights into maximum temperature of synthesis reactions in isothermal homogeneous semibatch reactors[END_REF] [25]. The results are given in Table 8. They show a considerable increase in MTSR from an isothermal of 40 °C. This can be explained by the fact that the value of the adiabatic temperature increase ∆0 12 used for the calculation of MTSR, considers the enthalpy of the nitration reaction in addition to the decomposition enthalpy.

MTSR = TP + ∆# $! . X ac, max 4-3
with X ac,max the maximum accumulation of the unreacted reactant and Tp the synthesis temperature.

In our case, the nitric acid is in stoichiometric excess over biodiesel. Thus X acc, max is obtained at the end of the addition. during the nitration process of biodiesel, NO3 or NO2 groups are grafted to the double bonds of carbon atoms. Based on this statement, biodiesel containing more double bonds, once nitrated, tend to be more energetic than others. That could explain the big different between these 2 results.

However, the difference between our additive and the EHN is not so high. Indeed a 2.5% and a 5% addition of EHN on BN20 (diesel with 20% of methyl ester) leads respectively to an increase in cetane number of 2.5 and 6.5 points [START_REF] Sathiyamoorthi | Panimalar Engineering College, poonamallee, chennai, Tamil Nadu. Evaluation of Neem oil biodiesel with 2-Ethyl Hexyl Nitrate 2-EHN as Cetane Improver on Performance and Emission characteristics of a DI Diesel Engine[END_REF].

Conclusion

The results of chemical analysis of the nitrated products showed the presence of 3 different compounds (methyl-8-nitro oleate, methyl-10-nitrate stearate and methyl-10-nitrate-9nitrostearate), whose selectivity varies according to the synthesis temperature. Some of these compounds, particularly the methyl-8-nitro oleate and the methyl-10-nitrate-9-nitrostearate, were also identified by Canoira et al. [START_REF] Canoira | Nitration of biodiesel of waste oil: Nitrated biodiesel as a cetane number enhancer[END_REF] and Malins et al. [START_REF] Malins | Nitration of methyl oleate with acetyl nitrate: a synthesis of methyl aminostearate[END_REF].

The heat released by the chemical transformation obtained with the reaction calorimeter RC1-RTCal varies with the synthesis temperature. This difference in the thermal behavior is due to the complexity of the reaction mechanism as well as the instability of the nitrated products and sulfo-nitric acid mixture. The study of the thermal stability shows that the Tonset of the decomposition is around 40°C. Thus, at high synthesis temperature (beyond the Tonset), the heat measured by the RC1 corresponds to the sum of the heat released by the nitration reaction and that of the mixture decomposition. This thermal and chemical behavior shows the complexity of the aliphatic compound nitration compared to the aromatic compounds one.

The modification of thermal behavior according to the synthesis temperature increases the risk of thermal runaway. The Maximum Temperature of the Synthesis Reaction (MTSR) obtained is higher than the sulfo-nitric acid temperature of decomposition, which leads to important risks of thermal runway. To avoid thermal runaway while having a high conversion rate, the synthesis temperature, the flow rate of sulfo-nitric acid as well as the temperature control of the cooling fluid are the main parameters on which we can lean on.

Fig. 1 .

 1 Fig. 1. Reaction calorimeter RC1-RTCal (Mettler®) Differential Scanning Calorimeter DSC Q20

Fig. 2 .Fig. 3 .RNO2

 23 Fig. 2. IR spectra of biodiesel and nitrated products at different temperatures

Fig. 4 .

 4 Fig. 4. Biodiesel and nitrated biodiesel mass spectrum

Fig. 5 .

 5 Fig.5. Power profiles at different temperatures

Fig. 6 :

 6 Fig. 6: Methyl palmitate IR spectrum before and after the reaction

Fig. 7 .Fig. 7 .

 77 Fig.7. Power profiles obtained with the palmitate Fig.7. shows that the power profiles are identical regardless of the temperature. The average heat is about 8J. All power profiles return to the baseline at the end of the addition, demonstrating the absence of any chemical transformations.

Fig. 8 .

 8 Fig.8. Pure sulfo-nitric acid and pure nitrated product decomposition Fig.9shows the DSC curves of sulfo-nitric acid mixed with nitrated products. The initial decomposition temperatures Tonset are 48.83°C, 70.39°C and 77.96°C respectively for 50 %, 70% and 90% of mixed acid. Obviously, when the quantity of nitrated products added into the mixed acid increased, lower onset temperatures are observed. It is also noted that the specific enthalpy decreases with the amount of nitrated products (table7). This enthalpy (∆H) reported

Fig. 9 .

 9 Fig.9. Study of the decomposition mixed sulfo-nitric and nitrated products

Table 1 :

 1 Fatty Acid Methyl Esters (FAME) from canola composition

	Compounds	% w/w
	Methyl palmitate
		2.2
	Methyl oleate
		92.0
	Methyl linoleate
		2.5
	Methyl linoleneate
		2.8

Table 3 :

 3 Molecules after the nitration reaction

Table 4 :

 4 Total conversion of the nitration reaction

	Temperature (°C)	Conversion rate (%)
	25	94.6
	30	92.7
	40	97.0
	50	89.2

Table 5 :

 5 Determination of the molecule's selectivity

	Temperature °C	Total conversion rate	Selectivity of B and C	Selectivity of A
		(% )	(%)	( %)
	25	94.6	95.2	4.8
	30	92.7	68.0	32.0
	40	97.0	42.9	57.1
	50	89.2	41.6	58.4

Table 6 :

 6 Heat released by the synthesis reaction at different temperatures

	Temperature °C	Heat (kJ)	Conversion rate (%)
	25	112.5	94.6
	30	106.3	92.8
	40	161.9	97.0
	50	171.9	89.2

Table 7 :

 7 Results of DSC analysis

	Component of sample (w/w)	Tonset (°C)	ΔH (J/g)
	90% of mixed acid and 10% nitrated products	77.96	502.0
	70% of mixed acid 30% nitrated products	70.39	474.1
	50% of mixed acid 50% nitrated products	48.83	332.1

Table 8 :

 8 Safety parameters of nitration reactionThe cetane number is the measurement of the diesel fuel ignition performance. It is obtained by making some comparison with reference fuels in a standardized test. Higher cetane number indicates that the ignition properties of the fuel are better[START_REF] Li | Effects of cetane number improvers on the performance of diesel engine fuelled with methanol/biodiesel blend[END_REF]. This index is obtained according to ASTM D 613. The equipment used was a Cooperative Fuel Research (CFR) engine. A standard diesel of cetane number 50 was used to prepare two samples of 1% and 5% concentration in bio-additive. For each mixture, the cetane number was measured at two synthesis temperatures 25 and 40 °C. The results of the cetane number analysis are listed in

	Temperature (°C)	Conversion (%)	Cp (J•g -1 •°C -1 )	∆0 12 (°C)	X ac,max (%)	MTSR (°C)
	25	94.6	2.1	194.4	30.3	83.8
	30	92.8	2.0	194.8	36.5	101.2
	40	97.0	2.2	262.8	49.5	171.1
	50	89.2	2.1	313.6	50.5	208.5
	If we compare results in table 7 with those in table 8, it appears that the lowest MTSR
	(83.8°C) is higher than the maximum Tonset obtained in decomposition study (77.96°C). In
	case of loss of temperature control, this means that the decomposition reaction would be
	triggered before completion of the synthesis reaction, which would further lead to very severe
	consequences. According to Stoessel's criticality index used for the risk assessment of a
	chemical reaction [27], this implies that the runaway is too fast for a safety barrier to be
	efficient. These results correspond to a Stoessel's criticality class of 5 if the nitration reaction
	conducted in these operated conditions.

Table 9 .

 9 They show that the bio-additive synthesized in this study allows to increase the cetane number of the reference diesel. A 5% addition of the bio-additive to the reference diesel improves the cetane number by 3 points.

Table 9 :

 9 Cetane number with the CFR engine

	Synthesis			CN fuel
	temperature (°C)	CN fuel	CN fuel +1% additive	+ 5% additive
	25	50.0	50.9	52.9
	40	50.0	50.9	53.2
	Canoira et al 2012 [13] obtained a 5 points increase in the fuel cetane number by mixing their
	nitrated biodiesel to fossil diesel in a concentration of 1000 mg.L -1 . With our bio-additive a
	5% addition (concentration of 49g.L

-1 ) to the reference diesel improves the cetane number by 3 points. It seems clear that our bio-additive is much less efficient than Canoira et al one. It may be due to the fact that Canoira et al used a biodiesel containing 33.0 %, 54.20% and 16% of methyl linoleate while ours contains only 2.5%. According to Abdullah et al 2010

[START_REF] Abdullah | PRODUCTION OF CETANE IMPROVER FROM Jathropa curcas OIL[END_REF]

,

Acknowledgments

The authors thank AMED project. The AMED project has been funded with the support of the European Union with the European Regional Development Fund (ERDF) and the Regional Council of Normandy.

The authors particularly thank Arnaud BARBIER for his help on chemical analyses.


Aina Rakotondrandisa
email: aina.rakotondrandisa@etu.univ-rouen.fr

Ionut Danaila
email: ionut.danaila@univ-rouen.fr

Luminita Danaila
email: luminita.danaila@coria.fr

Numerical modelling of a melting-solidification cycle of a phase-change material with complete or partial melting

Keywords: phase change, PCM, melting, solidification, Navier-Stokes, Boussinesq, FreeFem++, Newton method, finite element, mesh adaptivity

niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Hot, cold f Fusion

Introduction

The fundamental operational mode of latent thermal energy storage (LTES) systems based on phase-change materials (PCM) is made of alternate melting and solidification cycles that are not necessarily periodic. Partial melting and/or solidification of the PCM are often observed in applications and, in particular, in applications for buildings [START_REF] Zhu | Dynamic characteristics and energy performance of buildings using phase change materials: a review[END_REF][START_REF] Ascione | Energy refurbishment of existing buildings through the use of phase change materials: Energy savings and indoor comfort in the cooling season[END_REF]. Using the total latent heat storage potential offered by the PCM in energy storage requires a complex design process. This could benefit from accurate numerical simulations of such incomplete charging/discharging cycles. A wide range of recent applications is concerned by these modelling issues, including thermal energy storage (e. g. for solar power generation) and passive temperature control (e. g. for modern portable electronics) devices. For a review of various applications of PCMs with different melting temperatures in thermal energy storage systems, see recent reviews by [START_REF] Agyenim | A review of materials, heat transfer and phase change problem formulation for latent heat thermal energy storage systems (lhtess)[END_REF] and [START_REF] Kalnaes | Phase change materials and products for building applications: a state-of-theart review and future research opportunities[END_REF].

Actual challenges in the mathematical and numerical description of a melting-solidification cycle include i) the derivation of a veracious theoretical framework, or transport equations for different involved quantities (velocity, temperature, viscosity, density) and ii) the design of robust, accurate and efficient numerical methods for solving these equations, for different initial and boundary conditions.

As far as point i) is concerned, the solution of the Stefan problem has been provided by [START_REF] Rubinstein | On the solution of stefan's problem[END_REF]. Later, other important physical phenomena have been accounted for (gravity effects, convection in the liquid phase, the presence of a mushy region containing both solid and liquid parcels at the interface between the two phases, etc.). For a comprehensive review of these approaches, see [START_REF] Kowalewski | Phase change with convection: modelling and validation[END_REF] and [START_REF] Faghri | Transport Phenomena in Multiphase Systems[END_REF]. A historical review of the role played by convection in phase-change problems is provided by e.g. [START_REF] Yao | Melting and freezing[END_REF]. In particular, the natural convection in the liquid was proved to play an important role in the heat transfer between phases and in the propagation of the melting/solidification front [START_REF] Morgan | A numerical analysis of freezing and melting with convection[END_REF][START_REF] Voller | A fixed grid numerical modelling methodology for convection-diffusion mushy region phase-change problems[END_REF][START_REF] Jany | Scaling theory of melting with natural convection in an enclosure[END_REF][START_REF] Evans | Development of a 2-d algorithm to simulate convection and phase transition efficiently[END_REF][START_REF] Vidalain | An enhanced thermal conduction model for the prediction of convection dominated solid-liquid phase change[END_REF]Wang et al., 2010a). As a consequence, modern simulations of phase-change systems are dealing with the Navier-Stokes equations for the liquid phase, with the Boussinesq approximation for thermal effects.

Single domain approaches are very convenient for numerical implementation, since the same system of Navier-Stokes-Boussinesq equations is solved inside both liquid and solid phases. Two ingredients are necessary to make possible the use of the single domain model. First, the velocity inside the solid phase has to be set to zero. This is achieved by directly setting the velocity to zero in finite-volume methods (e. g. Wang et al. (2010a,b)) or by using penalty models in finite-elements methods, based on viscosity [START_REF] Danaila | A Newton method with adaptive finite elements for solving phase-change problems with natural convection[END_REF] or Carman-Kozeny terms [START_REF] Belhamadia | An enhanced mathematical model for phase change problems with natural convection[END_REF][START_REF] Zhang | Numerical study of melted pcm inside a horizontal annulus with threads in a three-dimensional model[END_REF][START_REF] Mencinger | Numerical simulation of melting in two-dimensional cavity using adaptive grid[END_REF]. Second, the energy equation is written using enthalpy method (Voller et al., 1987;[START_REF] Cao | A numerical analysis of Stefan problems for generalized multidimensional phase-change structures using the enthalpy transforming model[END_REF][START_REF] Cao | A numerical analysis of phase change problem including natural convection[END_REF]. The feature of the enthalpy method is its capability to deal with both mushy and single point phase changes. Indeed, in case of non-isothermal phase change, a mushy zone between the liquid and the solid phases characterize the system. However, in case of pure materials, the phase change occurs at a fixed temperature, and a numerical "regularization-zone" is introduced between the solid and the liquid parts, in order to regularize the enthalpy and other discontinuous parameters.

The challenge for numerical systems solving the single-domain Navier-Stokes-Boussinesq model is to accurately capture the moving solid-liquid interface. The problem is even more challenging when several melting-solidification fronts, with distorted shapes are present (e. g. the solidification after a partial melting). When fixed uniform meshes are used, which is the case of a great majority of existing finite-volume codes, the grid density has to be considerably increased in the entire domain, making the simulation very expensive. When the trade-off between accuracy and computational cost is sought, the fixed grid approach allows to place only a few computational cells inside the regularization region.

Dynamical mesh adaptivity becomes in this context a valuable tool to concentrate the grid refinement effort only in regions displaying high gradients of the computed variables (melting-solidification fronts, thermal or viscous boundary layers). For the classical twophase Stefan problem, Belhamadia et al. (2004a) suggested an anisotropic mesh adaptation algorithm based on solution-dependent metrics. The authors extended their algorithm for the three-dimensional simulation of the same problem (Belhamadia et al., 2004b) and showed that the use of locally adaptive meshes with strong anisotropy proved to be very effective in reducing the number of computational nodes for such phase-change systems without convection. To simulate melting or solidification problems with convection, [START_REF] Danaila | A Newton method with adaptive finite elements for solving phase-change problems with natural convection[END_REF] recently suggested a dynamical mesh adaptation algorithm based on metrics control and implemented with the FreeFem++ software [START_REF] Hecht | FreeFem++ (manual)[END_REF][START_REF] Hecht | New developments in Freefem++[END_REF]. The advantage of this adaptive finite-element method, that will be also used in the present study, is to make possible, with reasonable computational cost, the re-meshing of the computational domain at each time step. A very refined discretization of the regularization zone between solid and liquid phases is thus obtained, while regions with low gradients are de-refined in order to balance the overall computational effort.

The previously mentioned modern numerical approaches were mostly applied to simulate separately melting or solidification problems and only recently for alternate melting and solidification complete cycles (Wang et al., 2010b). However, cyclic, or periodic, melting and solidification problems have attracted considerable attention in the literature. [START_REF] Ho | Periodic melting within a square enclosure with an oscillatory surface temperature[END_REF] and [START_REF] Voller | Cyclic phase change with fluid flow[END_REF] studied numerically, periodic melting in a square enclosure. Recently, [START_REF] Hosseini | Experimental and computational evolution of a shell and tube heat exchanger as a pcm thermal storage system[END_REF] presented experimental studies for the melting and the solidification of a cylindrical PCM during charging and discharging process and [START_REF] Chabot | Solid-liquid phase change around a tube with periodic heating and cooling: Scale analysis, numerical simulations and correlations[END_REF] studied analytically the effect of an alternate heating and cooling in a cylindrical PCM, with periodical boundary conditions.

The present contribution is scoped to offer an accurate numerical description of the alternate melting and solidification of a typical PCM. We use a finite-element numerical system with second-order accuracy in time and space to solve the single-domain model based on Navier-Stokes equations with Boussinesq approximation. The main advantage of our method is that the mesh adaptivity algorithm could be applied each time step. The mesh is thus dynamically refined with respect to velocity and temperature variables, allowing to accurately capture the interface between solid and liquid phases, the boundarylayer structure at the walls and the details of the unsteady convection cells in the liquid.

We simulate a typical PCM configuration represented by a differentially heated square cavity filled with an octadecane paraffin. This is a well-established benchmark documented experimentally by [START_REF] Okada | Analysis of heat transfer during melting from a vertical wall[END_REF] and extensively used to validate numerical codes (Wang et al., 2010a;[START_REF] Jany | Scaling theory of melting with natural convection in an enclosure[END_REF][START_REF] Mencinger | Numerical simulation of melting in two-dimensional cavity using adaptive grid[END_REF]. First, we simulate the melting phase and use this case to validate our numerical system against experimental and previously reported numerical data. Second, we consider two operating cases for the solidification process. In the first study case the solidification starts after a complete melting of the PCM (liquid fraction of 95%), while in the second case after a partial melting (liquid fraction of 50%). All cases are analysed in detail by providing temporal evolution of solid-liquid interface, liquid fraction, Nusselt number and accumulated heat input. Different heat transfer regimes are identified and explained using scaling correlation theory. Several practical implications for the two operating modes are finally drawn.

The paper is organized as follows. Section 2 introduces the governing equations. Section 3 presents the numerical system. The final section 4 is devoted to extensive analysis of the results for the two operating cases.

Governing equations

We consider a solid-liquid system placed in a two-dimensional square cavity of height H. In the following, subscripts s and l will refer to the solid and liquid phases, respectively.

For the numerical implementation, it is convenient to adopt a single-domain approach to describe both phases using the same system of equations. The model is based on the Navier-Stokes equations with Boussinesq approximation, which is the natural description of the fluid flow with natural convection. A penalty term is then added to the momentum equations to bring the velocity to zero inside the solid region. For the energy conservation equation, an enthalpy method is used to model the phase change process. The singledomain model is described in detail in the following sections. 4

Enthalpy method

The phase change process is modelled using an enthalpy method (Voller et al., 1987;[START_REF] Cao | A numerical analysis of Stefan problems for generalized multidimensional phase-change structures using the enthalpy transforming model[END_REF][START_REF] Cao | A numerical analysis of phase change problem including natural convection[END_REF] with temperature-based formulation. We start from the energy equation:

∂(ρh) ∂t ϕ + ∇ • (ρhũ) -∇ • (k∇T ) = 0, ( 1
)
where t ϕ is the physical time, h the enthalpy, ρ the density, ũ the velocity vector, T the temperature and k the thermal conductivity. To make Equation (1) valid for the entire domain containing both liquid and solid phases, the total enthalpy h is regarded as the sum of the sensible heat and the latent heat:

h = c(T + s(T )), (2)
with c the local specific heat. The function s(T ) is introduced to model the jump of the enthalpy due to the phase change and is theoretically a Heaviside step function depending on the temperature: it takes the zero value in the solid region and a large value in the liquid, equal to h sl /c, with h sl the latent heat of fusion. Linear (Voller et al., 1987;Wang et al., 2010a) or smoother functions [START_REF] Danaila | A Newton method with adaptive finite elements for solving phase-change problems with natural convection[END_REF] can be used to regularize s(T ) and also the jump of material properties (from solid to liquid). In this paper we use a regularization of all step-type functions by a continuous and differentiable hyperbolictangent function suggested by [START_REF] Danaila | A Newton method with adaptive finite elements for solving phase-change problems with natural convection[END_REF] (see below). We assume moreover the undercooling problem to be negligible (see also Wang et al. (2010b); [START_REF] Kowalewski | Phase change with convection: modelling and validation[END_REF]). Equation (1) can be further simplified by considering the following assumptions: (i) the density difference between solid and liquid phases is negligible, i. e. ρ l = ρ s = ρ is constant; (ii) the regularization zone is narrow and the velocity inside this zone is very low. Consequently, the final expression of the energy equation is obtained by combining (2) and (1) and neglecting the convection term ∇ • (csũ)1 :

∂ (cT ) ∂t ϕ + ∇ • (cT ũ) -∇ • k ρ ∇T + ∂ (cs) ∂t ϕ = 0. (3)
Furthermore, the essential feature of the current approach is that the phase change front is not tracked explicitly but is instead recovered a posteriori from the computed temperature field.

Navier-Stokes equations with Boussinesq approximation

The natural convection in the liquid part of the system is modelled using the incompressible Navier-Stokes equations, with Boussinesq approximation for buoyancy effects. To make this model valid for both liquid and solid phases, the momentum equation is modified as follows:

∂ ũ ∂t ϕ + (ũ • ∇)ũ + 1 ρ ∇p - µ l ρ ∇ 2 ũ -f B (T )e y = A(T )ũ, (4)
where p denotes the pressure, µ l the dynamic viscosity of the liquid (assumed to be constant) and f B (T ) the Boussinesq force. The penalty term A(T )ũ is artificially introduced in (4) to extend this equation in the solid phase, where the velocity, pressure, viscosity and Boussinesq force are meaningless. Consequently, A(T ) is modelled to vanish in the liquid, where the Navier-Stokes-Boussinesq momentum equation is recovered. A large value of A(T ) is imposed in the solid, reducing the momentum equation (4) to A(T )ũ = 0, equivalent to ũ = 0. Exact expressions for f B and A will be given in the next section.

Finally, the conservation of mass in the liquid phase is expressed by the continuity equation:

∇ • ũ = 0.
(5)

Final system of equations for the single-domain approach

It is convenient to numerically solve a dimensionless form of the previous equations. Using the cavity height H as length scale and a reference state (ρ, V ref , T f ), we can define the following scaling for the space, velocity, temperature and time variables:

x = x H , u = ũ V ref , θ = T -T f δT , t = V ref H t ϕ , (6)
Temperatures T h (hot) and T c (cold) will be used to set isothermal walls of the cavity.

The difference δT = T h -T f , with T f the temperature of fusion, is considered as the representative temperature scale for the natural convection onset in the liquid region.

As far as the solidification process is concerned, a distinct discussion will be provided in section 4.3. Thus δT is used to define the Rayleigh number of the flow:

Ra = gβH 3 δT ν l α l , ( 7
)
where α = k/(ρc) is the thermal diffusivity and β the thermal expansion coefficient. Note that the reference temperature in this scaling is T f , resulting in θ f = 0. This simplifies the identification of the regularization zone, defined for -ε ≤ θ ≤ ε.

Finally, the dimensionless system of equations to be solved in both liquid and solid regions can be written as:

∇ • u = 0, (8) ∂u ∂t + (u • ∇)u + ∇p - 1 Re ∇ 2 u -f B (θ) e y -A(θ)u = 0, (9)
∂ (Cθ) ∂t + ∇ • (Cθu) -∇ • K Re Pr ∇θ + ∂ (CS) ∂t = 0, ( 10
)
where the linearised (Boussinesq) buoyancy force (f B ), the Reynolds (Re) and Prandtl (Pr) numbers are defined as:

f B (θ) = Ra Pr Re 2 θ, Re = ρV ref H µ l = V ref H ν l , Pr = ν l α l . ( 11
)
Non-dimensional conductivity and specific heat are functions of the temperature θ,

K(θ) = k k l , C(θ) = c c l , ( 12
)
and have to take into account the variation of material properties between the solid and the liquid regions.

In energy equation ( 10), the non-dimensional function S = s/δT , introduced by the enthalpy model, is regularized across the regularization region using a hyperbolic-tangent function [START_REF] Danaila | A Newton method with adaptive finite elements for solving phase-change problems with natural convection[END_REF]:

S(θ) = S l + S s -S l 2 1 + tanh θ r -θ R r , ( 13
)
where θ r is the central value around which we regularize (typically θ r = θ f = 0) and R r the smoothing radius (typically R r = ε). Note that S s = 0 and

S l = h sl /c l δT = 1 Ste , ( 14
)
with Ste the Stefan number. Regularizations similar to ( 13) can be used to model the variation inside the regularization region of functions ( 12) defining material properties.

Finally, the penalty term in the momentum equation ( 9) takes the form [START_REF] Belhamadia | An enhanced mathematical model for phase change problems with natural convection[END_REF][START_REF] Kheirabadi | The effect of the mushy-zone constant on simulated phase change heat transfer[END_REF]:

A(θ) = - C CK (1 -L F (θ)) 2 L F (θ) 3 + b , (15)
where L F (θ) is the liquid fraction, which is 1 in the fluid region and 0 in the solid. Inside the regularization region, L F (θ) is regularized using a hyperbolic-tangent similar to (13).

The constant C CK is set to a large value (as discussed below) and the constant b = 10 -6 is introduced to avoid division by zero.

Numerical method

To solve the system of equations ( 8)-( 10) we use a finite-element method that was implemented using the open-source software FreeFem++ [START_REF] Hecht | FreeFem++ (manual)[END_REF][START_REF] Hecht | New developments in Freefem++[END_REF], using a large variety of triangular finite elements to solve partial differential equations.

FreeFem++ is an integrated product with its own high level programming language and a syntax close to mathematical formulations, making the implementation of numerical algorithms very easy. Among the numerous numerical tools offered by FreeFem++, the use of the powerful mesh adaptivity function proved mandatory in this study to obtain accurate results within reasonable computational time. The numerical code was optimized to afford the mesh refinement every time step: the mesh density was increased around the phase change interfaces, offering an optimal resolution of the large gradients of all regularized functions (S, K, C, L F ), while the mesh was de-refined (larger triangles) in the solid part, where a coarser mesh could be used. A simulation using a globally refined mesh would require a prohibitive computational time for an equivalent accuracy of the melting front resolution. Similar algorithms based on FreeFem++ were successfully used for solving different systems of equations with locally sharp variation of the solution, such as Gross-Pitaevskii equation [START_REF] Danaila | A finite element method with mesh adaptivity for computing vortex states in fast-rotating Bose-Einstein condensates[END_REF][START_REF] Vergez | A finite-element toolbox for the stationary Gross-Pitaevskii equation with rotation[END_REF] or Laplace equations with nonlinear source terms [START_REF] Zhang | Existence and numerical modelling of vortex rings with elliptic boundaries[END_REF].

The space discretization is based on Taylor-Hood finite elements, approximating the velocity with P 2 Lagrange finite elements (piecewise quadratic), and the the pressure with the P 1 finite elements (piecewise linear). The temperature and the enthalpy are discretized using P 1 finite elements. The weak formulation of ( 8)-( 10), necessary for the finite-element implementation is described in detail in [START_REF] Danaila | A Newton method with adaptive finite elements for solving phase-change problems with natural convection[END_REF]. There are two novelties in the present numerical approach, when compared to [START_REF] Danaila | A Newton method with adaptive finite elements for solving phase-change problems with natural convection[END_REF]: (i) we use an approach based on the Carman-Kozeny model to bring the velocity to zero inside the solid phase, as described in the previous section, instead of the viscosity penalty method (imposing a large value of the viscosity in the solid); (ii) we increase the time accuracy of the algorithm by replacing the first-order Euler scheme with the second-order Gear (BDF2) scheme (see also [START_REF] Belhamadia | An enhanced mathematical model for phase change problems with natural convection[END_REF],

dφ dt 3φ n+1 -4φ n + φ n-1 2δt , (16)
computing the solution φ n+1 at time t n+1 = (n + 1)δt by using two previous states (φ n , φ n-1 ). We use this scheme to advance in time both velocity (φ = u) and temperature fields (φ = θ). The other terms in equations ( 8)-( 10) are treated implicitly (i. e. taken at time t n+1 ). The resulting non-linear equations are solved using a Newton algorithm.

Mesh adaptivity

Mesh adaptivity by metric control is a standard function offered by FreeFem++ Hecht (2012). The key idea for the mesh adaptivity (see also [START_REF] Castro-Diaz | Anisotropic grid adaptation for inviscid and viscous flows simulations[END_REF]; [START_REF] Hecht | Mesh adaptation by metric control for multi-scale phenomena and turbulence[END_REF]; [START_REF] George | Delaunay triangulation and meshing[END_REF]) is to modify the scalar product used in an automatic mesh generator to evaluate distance and volume, in order to construct equilateral elements according to a new adequate metric. The scalar product is based on the evaluation of the Hessian H of the variables of the problem. Indeed, for a P 1 discretization of a variable χ, the interpolation error is bounded by:

E = |χ -Π h χ| 0 ≤ c sup T ∈T h sup x,y,z∈T |H(x)|(y -z, y -z) (17)
where Π h χ is the P 1 interpolate of χ, |H(x)| is the Hessian of χ at point x after being made positive definite. We can infer that, if we generate, using a Delaunay procedure (e.g. [START_REF] George | Delaunay triangulation and meshing[END_REF]), a mesh with edges close to the unit length in the metric M = |H| (cE) , the interpolation error E will be equally distributed over the edges a i of the mesh. More precisely, we have

1 cE a T i Ma i ≤ 1. ( 18
)
The previous approach could be generalized for a vector variable χ = [χ 1 , χ 2 ]. After computing the metrics M 1 and M 2 for each variable, we define a metric intersection M = M 1 ∩ M 2 , such that the unit ball of M is included in the intersection of the two unit balls of metrics M 2 and M 1 (for details, see the procedure defined in [START_REF] Frey | Maillages. Hermès[END_REF]).

For the cases considered in this study, we used five metrics intersections to adapt the mesh, based on S n+1 , S n , T n+1 , T n , u n+1 . Note that to reduce the impact of the interpolation on the global accuracy for time-depending problems, for the same variable used for adaptivity, we consider the metrics computed from actual (at t n+1 ) and previous (at t n ) values (see also Belhamadia et al. (2004a)). The anisotropy of the mesh is a parameter of the algorithm and it was set to values close to 1. This is an inevitable limitation since we also impose the minimum edge-length of triangles to avoid too large meshes.

Some examples of the adapted mesh during computations are given in section 4.3 during the solidification phase. The method is able to accurately capture the liquid-solid interface during the melting process and the two solidification fronts appearing during the solidification of the PCM. Mesh adaptivity is performed at each time step and offers a refined discretization of the regularization region where sharp gradients have to be accurately captured. The number of triangles for the melting case is N t = 12, 000 and N t = 17, 000 for the solidification phase. Non-adapted grids offering the same spatial resolution everywhere inside the computational domain would have resulted for the two cases in N t = 9.94 • 10 10 and N t = 10 • 10 10 triangles, respectively. Consequently, mesh adaptivity also greatly helps in reducing the computational time.

Results

The basic configuration considered in this study is the experimental study of [START_REF] Okada | Analysis of heat transfer during melting from a vertical wall[END_REF]. It consists of a differentially heated square cavity (see Fig. 1a), filled with an octadecane paraffin. The physical (non-dimensional) parameters are: Ra = 3.27 • 10 5 , Pr = 56.2 and Ste = 0.045. We numerically investigate the following cases: (i) First we proceed to the validation of our numerical method on the basis of the melting of the PCM. The material is initially solid (θ 0 = -0.01) and melts progressively starting from the left boundary, maintained at the hot temperature θ h = 1. The right boundary is also isothermal with cold temperature θ c = -0.01. This case is a well established benchmark used to validate numerical codes for phase-change systems. Three simulations are carried out to demonstrate the accuracy and the robustness of our method. The first validation case reproduces the experimental study of [START_REF] Okada | Analysis of heat transfer during melting from a vertical wall[END_REF], the second corresponds to the experimental investigation of [START_REF] Gong | Numerical simulation of melting with natural convection based on lattice boltzmann method and performed with cuda enabled gpu[END_REF] and the last considers the numerical simulations presented in [START_REF] Bertrand | Melting driven by natural convection a comparison exercise: first results[END_REF]. For the last cases, the melting front, the Nusselt number and the liquid fraction provided by our code are compared with experimental and numerical results.

(ii) The melting of the PCM (Fig. 1b). Having validated our code, we pay a closer attention to other physical characteristics of the melting, such as the temporal evolution of the temperature distribution in the cavity, the Nusselt number, etc. Ẇe discuss and favorably compare our results with those previously published: experimental [START_REF] Okada | Analysis of heat transfer during melting from a vertical wall[END_REF] and numerical [START_REF] Okada | Analysis of heat transfer during melting from a vertical wall[END_REF]Wang et al., 2010a;[START_REF] Ma | Solid velocity correction schemes for a temperature transforming model for convection phase change[END_REF][START_REF] Danaila | A Newton method with adaptive finite elements for solving phase-change problems with natural convection[END_REF]. The scaling formulae suggested by [START_REF] Jany | Scaling theory of melting with natural convection in an enclosure[END_REF] nicely collapse with the Nusselt number of our simulation. This is described in Section 4.2.

(iii) After the complete melting of the PCM (i. e. the melting front is very close to the right wall and the liquid fraction is 0.95), we trigger the solidification process by changing the left-wall temperature to θ c = -0.01 whilst the right wall is still kept at θ = -0.01. The solid phase will propagate into the cavity from both left and right sides (Fig. 1c). This case is computationally challenging, since two melting/solidification fronts have to be accurately followed during the simulation. The process is simulated up to the complete solidification of the PCM and the non-trivial evolution of the liquid phase is depicted in detail. This case is described in Section 4.3.1.

(iv) This case is similar to the previous one, but the solidification starts after a partial melting of the PCM (i. e. the melting front is located approximatively at half distance between the two vertical walls and the liquid fraction is 0.5). The analysis of this case attempts to provide answers about the effectiveness of different possible functioning cycles of the PCM. This case is described in Section 4.3.2.

Numerical validation

We validate our numerical method against experimental and numerical studies of the melting of the octadecane PCM available in the literature. Three cases are investigated. The first consists of an experimental study of the melting of the PCM in a differentially heated square cavity of height H = 1.5 cm by [START_REF] Okada | Analysis of heat transfer during melting from a vertical wall[END_REF]. The second reproduces the melting of the PCM included in a transparent building brick of height H = 15.2 cm, investigated experimentally and numerically by [START_REF] Gong | Numerical simulation of melting with natural convection based on lattice boltzmann method and performed with cuda enabled gpu[END_REF]. The last compares our results with various numerical methods, presented by [START_REF] Bertrand | Melting driven by natural convection a comparison exercise: first results[END_REF], simulating the melting of octadecane, considering a higher value of the Rayleigh number. For the simulation of the melting process (Fig. 1b), we use the following choice for the scaling introduced in §2.3, equations ( 6) and ( 11):

V ref = ν l H ⇒ t = t ϕ ν l H 2 ⇒ Re = 1. (19)
Moreover, a second dimensionless time τ is introduced in order to assess our results with respect to the numerical data of [START_REF] Bertrand | Melting driven by natural convection a comparison exercise: first results[END_REF] and the analytical correlation of [START_REF] Jany | Scaling theory of melting with natural convection in an enclosure[END_REF]:

τ = Ste • F o = Ste • αt ϕ H 2 = Ste • t P r , ( 20
)
where F o is the Fourier number.

We first examine the location of the interface obtained with our simulations. The comparison with the experimental results of [START_REF] Okada | Analysis of heat transfer during melting from a vertical wall[END_REF] and [START_REF] Gong | Numerical simulation of melting with natural convection based on lattice boltzmann method and performed with cuda enabled gpu[END_REF] is presented in Figure 2. The experimental study of [START_REF] Okada | Analysis of heat transfer during melting from a vertical wall[END_REF] in Figure 2(a) consists of a differentially heated square cavity of dimensions 1.5 cm × 1.5 cm, filled with an octadecane paraffin. The non-dimensional parameters are: Ra = 3.27 • 10 5 , Pr = 56.2 and Ste = 0.045.

For two particular time instants (τ = 0.032 and τ = 0.063), we could compare to available experimental [START_REF] Okada | Analysis of heat transfer during melting from a vertical wall[END_REF] and numerical [START_REF] Okada | Analysis of heat transfer during melting from a vertical wall[END_REF]Wang et al., 2010a;[START_REF] Danaila | A Newton method with adaptive finite elements for solving phase-change problems with natural convection[END_REF] data. In the experimental set up of [START_REF] Okada | Analysis of heat transfer during melting from a vertical wall[END_REF], the author have reported that the top of the PCM was not perfectly insulated and consequently the growth of the experimental upper melting front was delayed. In Figure 2(a), for the two time instants τ = 0.032 and τ = 0.063, the current work agrees well with the experimental results of [START_REF] Okada | Analysis of heat transfer during melting from a vertical wall[END_REF] at the bottom part of the melting front. However, our results overestimate the location of the front in the top part of the cavity, which is in total agreement with the experimental heat loss mentioned by the author.

Moreover, our results are qualitatively in a better agreement with experimental data than previously published numerical results. This is a direct consequence of the precise tracking of the melting front achieved by the mesh adaptivity performed at each time step. This assessment also allowed us to finely tune the value of the constants used in the model ( 15). Even though it is generally assumed that a large value for C CK must be set, the exact value of this constant could influence the accuracy of the results [START_REF] Kheirabadi | The effect of the mushy-zone constant on simulated phase change heat transfer[END_REF][START_REF] Kumar | Influence of mushy zone constant on thermohydraulics of a pcm[END_REF]. This choice of the value of this constant is a still open problem. Very good agreement with the experimental result of [START_REF] Okada | Analysis of heat transfer during melting from a vertical wall[END_REF] is obtained for C CK varying in the range [10 6 , 10 8 ]. Nevertheless, imposing a too large value C CK = 10 10 results in artificially slowing the propagation of the melting front. We set for all subsequent simulations C CK = 10 6 . Figure 2(b) illustrates the interface location in the experiment and simulations of [START_REF] Gong | Numerical simulation of melting with natural convection based on lattice boltzmann method and performed with cuda enabled gpu[END_REF], who studied the melting of an octadecane PCM inside a transparent building brick of dimensions 15.2cm×3cm. Their numerical simulation has been performed using a Lattice Boltzmann method. The non-dimensional parameters were: Ra = 2.48 • 10 8 , Pr = 50 and Ste = 0.072. The difficulty here compared to the first validation case is the presence of a stronger natural convection flow in the fluid due to the high value of the Rayleigh number. The location of the interface is compared for five particular time instants: τ = 0.0002, 0.00050, 0.00067, 0.00125 and 0.00252. We can notice thus a very good agreement with the numerical and the experimental data of [START_REF] Gong | Numerical simulation of melting with natural convection based on lattice boltzmann method and performed with cuda enabled gpu[END_REF].

A last validation case is also investigated to test the robustness of the method. The physical parameters are: Ra = 10 8 , Pr = 50 and Ste = 0.1. [START_REF] Bertrand | Melting driven by natural convection a comparison exercise: first results[END_REF] compiled results provided by five different authors (Lacroix, Le Quéré, Gobin-Vieira, Delannoy and Binnet-Lacroix). Results provided by these investigators will be hereafter referred to as (say) 'Lacroix, from [START_REF] Bertrand | Melting driven by natural convection a comparison exercise: first results[END_REF]'. They have attempted a first comparison by taking several numerical methods to compute the basic configuration presented in this section. Two investigators among the five failed to predict the process and showed unrealistic behaviors in Figures 3 and4: Lacroix and Delannoy seem to be insufficiently converged as shown by Figure 3, and Binet-Lacroix overestimates the average Nusselt number by more than 30% (Figure 4). Hence, this collection of results allows us to compare our numerical method with the existant code and check whether or not realistic results are obtained for complex physical configurations. We further inspect the melting front, the temporal evolution of the liquid fraction L f and the Nusselt number Nu at the left wall (x = 0), for each of the five methods presented by [START_REF] Bertrand | Melting driven by natural convection a comparison exercise: first results[END_REF]. For the liquid fraction, the initial solid state corresponds to L f = 0, while L f = 1 indicates the complete melting of the PCM. The Nusselt number Nu is defined as follows:

Nu = 1 0 ∂θ ∂x x=0 dy. ( 21
)
The phase-change interface for four time steps, τ = 5 • 10 -4 , τ = 2 • 10 -3 , τ = 6 • 10 -3 and τ = 1 • 10 -2 is represented in Figure 3. Our results are for each case in fairly good agreement with those of Gobin and Le Quéré. Details on the method and the numerical procedures are listed in [START_REF] Gobin | Melting from an isothermal vertical wall. synthesis of numerical comparison exercise[END_REF]. Gobin uses a front-tracking method using a coordinate transformation with a finite volume method in a 62 × 42 grids. Le Quéré solves a single domain method using a second order scheme with a finite volume method in a 192 × 192 grid points.

The time evolution of the Nusselt number and the liquid fraction are presented in Figure 4. A very good agreement is still obtained with Gobin and Le Quéré. A relative difference, less than 2% is noticed for the Nusselt number, and a dispersion smaller than 4% for the melted fraction.

The high Rayleigh number Ra = 10 8 is a very demanding numerical test. The high velocity, inducing a very narrow thermal boundary layer can lead to unrealistic results and some numerical methods have failed. The interest of the mesh adaptation is clearly evidenced since we initially use only 40 × 40 grid points.

Melting of the PCM

After the validation of the code, we now pay a closer attention to the temporal evolution of different physical parameter of the system, during the melting phase. We recall that the reference temperature is the fusion temperature and thus θ f = 0. The regularization range is defined for -ε 1 ≤ θ ≤ ε 2 , with ε 1 = ε 2 = 0.01. For this case, physical properties of the material are identical in both liquid and solid phases, and, consequently, we obtain from (12) that C(θ) = 1 and K(θ) = 1. This choice of the scaling was made to have the same set of parameters as in previous numerical simulations of this case [START_REF] Okada | Analysis of heat transfer during melting from a vertical wall[END_REF]Wang et al., 2010a;[START_REF] Ma | Solid velocity correction schemes for a temperature transforming model for convection phase change[END_REF][START_REF] Danaila | A Newton method with adaptive finite elements for solving phase-change problems with natural convection[END_REF]. The physical parameters correspond to the basic configuration presented in Figure 1: Ra = 3.27 • 10 5 , Pr = 56.2 and Ste = 0.045.

Time evolution

We start by analysing the time evolution of the melting process. At τ = 0, the material is solid and the initial temperature is set to θ 0 = -0.01 everywhere inside the cavity. Then, the temperature of the left wall is suddenly increased to θ h = 1, while the right wall is maintained at the same cold temperature θ c = -0.01. The material starts to melt, with a melting (identified by the iso-line θ = θ f = 0) propagating from the left to the right side of the domain. The time evolution of the phase-change system is depicted in Fig. 5 for representative time instants, also reported in previous studies.

From Fig. 5, we can easily identify three different regimes describing the time evolution of the melting process.

• From τ = 0 to τ = 0.004 (Fig. 5a), we note the vertical shape of the melting front, well predicted by the classical conduction model of [START_REF] Stefan | Über die theorie der eisbildung, insbesondere über die eisbildung im polarmeere[END_REF]. This indicates that, at this stage, heat transfer is dominated solely by conduction.

• Between τ = 0.016 to τ = 0.032 (Fig. 5b), the natural convection in the fluid phase starts to alter the shape of the melting front. A mixed conduction and convection regimes rule the heat transfer. Convection mainly affects the upper part of the fluid motion, while conduction is still dominating in the lower part. As the volume thermal expansion coefficient β is positive, we expect a clockwise circulation of the liquid inside the convection cell, as noted by [START_REF] Jany | Scaling theory of melting with natural convection in an enclosure[END_REF]. This also makes the liquid-solid interface to move faster at the top of the cavity, explaining the deformed shape of the melting front, which is a signature of the convection effects (see also [START_REF] Kowalewski | Phase change with convection: modelling and validation[END_REF]).

• After τ = 0.032 (Fig. 5c-d), natural convection dominates the heat transfer process and impacts radically the solid-liquid interface shape and motion. The melting front line exhibits four distinct regions characterized by different slopes with respect to the vertical axis. The largest slope is observed at the top of the cavity and is related to the particular shape of the convection cell. Note that top and bottom parts of the interface are normal to the cavity boundaries because of the imposed adiabatic boundary conditions. • After τ = 0.08 the melting front is nearly touching the right wall of the cavity, firstly at the top (Fig. 5e) of the cavity. The melting process continues and the fluid progressively fills the cavity, with a melting front deforming to a vertical line.

The simulation of the melting process is stopped at τ = 0.2 (Fig. 5f), when it is numerically difficult to separate the melting front from the right wall boundary. At this time instant, the fluid fraction reaches the value of 0.95 and the melting of the PCM is considered to be complete, even though a small region of solid PCM remains at the lower right bottom of the cavity. Note from Fig. 5f the existence in the fluid of two recirculating zones instead of a single one observed during previous stages.

Scaling analysis

The melting of the octadecane was theoretically studied by [START_REF] Jany | Scaling theory of melting with natural convection in an enclosure[END_REF]. Combining scaling theory and numerical modelling, they suggested closed-form correlations for the temporal evolution of the average Nusselt number Nu defined at the hot boundary (x = 0), under the form:

Nu(τ ) = 1 √ 2τ + c 1 Ra 1/4 - 1 √ 2τ 1 + c 2 Ra 3/4 τ 3/2 n 1/n . ( 22
)
The values of the constants were fitted from numerical data: c 1 = 0.27, c 2 = 0.0275, and n = -2.

In Fig. 6 we compare the time evolution of the Nusselt number obtained from our numerical data (see Figure 1) to experimental results of [START_REF] Okada | Analysis of heat transfer during melting from a vertical wall[END_REF] and predictions obtained from the correlation (22). Our results perfectly fit the theoretical prediction of [START_REF] Jany | Scaling theory of melting with natural convection in an enclosure[END_REF]. They are also in good agreement with experimental data, suggesting, however, that very accurate measurements and numerical simulations are needed to validate theoretical scaling analysis.

The time evolution of the Nusselt number can be correlated with the different heat transfer regimes analysed in the previous section:

1. A pure conduction regime for τ 0 (corresponding to Fig. 5a), characterized by the law Nu ∼ (2τ ) -1/2 . Since the temperature gradient has initially huge values because of the sudden increase of the temperature of the left wall, the Nusselt number rapidly decreases during the first stage of the flow evolution. The evolution law Nu ∼ τ -1/2 can be also obtained from the Neumann exact solution of [START_REF] Grigull | Heat conduction[END_REF]. The signature of this conduction regime is the slow heat transfer characterized by a monotonic decrease of the Nusselt number, down to a minimum value obtained for τ ∼ Ra -1/2 = 0.02. 2. A mixed conduction-convection regime for 0.02 ≤ τ ≤ 0.05 (illustrated in Fig. 5b).

The influence of the Rayleigh number in ( 22) starts to be important and a good approximation for this regime is: Nu ∼ τ -1/2 + Ra τ 3/2 . 3. A convection dominated regime for τ > Ra -1/2 (corresponding to Figs. 5c-e). In the asymptotic limit of large τ , the simplified law Nu ∼ Ra 1/4 is obtained. The plateau at the value of Ra 1/4 corresponds to the pure convective transfer and is observed in Fig. 6 for 0.05 ≤ τ ≤ 0.1. Numerical results show a slight decrease of the Nu in the final stage (τ ≥ 0.1), when the melting front starts to touch the right wall of the cavity (see Figs. 5e-f). The correlation model of [START_REF] Jany | Scaling theory of melting with natural convection in an enclosure[END_REF] is not valid for this late evolution of the melting process.

Another important basic quantity describing the melting process is the liquid fraction L f . The time evolution of the liquid fraction (Fig. 7a) displays three regimes during the melting process. L f initially grows as τ 1/2 , which is a typical law for a conductiondominated heat transfer. Then, a linear temporal evolution is observed, until the melting front reaches the right wall. This linear regime corresponds to the quasi-steady state observed in the evolution of the Nusselt number (Fig. 6).

Using the asymptotic limits of Eq. ( 22) for τ → 0 (pure conduction) and τ → ∞ (pure convection), [START_REF] Jany | Scaling theory of melting with natural convection in an enclosure[END_REF] suggested the following correlation law for the time evolution of the liquid fraction:

L f (τ ) = √ 2τ 5 + c 1 Ra 1/4 τ 5 1/5 , (23)
where c 1 = 0.27 is the same constant as in ( 22). We compare in Fig. 7b our numerical results with the predictions based on (23) within the validity domain of the analysis, i. e. before the melting front reaches the right wall of the cavity. A very good agreement is found with theoretical predictions and also with previously published numerical results (Wang et al., 2010a).

Influence of the Rayleigh number

To assess the influence of the Rayleigh number on the evolution of the melting process, we performed two other simulations by multiplying the initial value of Ra = 3.27 • 10 5 by a factor of 5 and 10, respectively. The exact values are: Ra = 1.62 • 10 6 and Ra = 3.27 • 10 6 . First, we increase the height H of the cavity by a factor of 3 √ 5 and 3 √ 10 and consider the same δT . Thus the Stefan number Ste is kept constant. Second, we increase the temperature difference parameter δT by keeping H constant. It corresponds of an increased value of the Stefan number by a factor of 5 and 10: Ste = 0.223 and Ste = 0.45. Figures 8 and9 show the temporal evolution of the liquid fraction L f (a), and the average Nusselt number defined at the hot wall, (b). The same heat transfer regimes described previously are observed for each case: conduction, mixed conduction-convection and convection.

Figure 8(a) indicates that increasing the Rayleigh number by keeping δT constant induces a slower melting rate. This is the expected behaviour since the size of the PCM is increased by a factor of 2, and the velocity u is hence decreasing in order to satisfy the condition Re = 1. We note however a non-monotonic variation of time necessary to melt a fixed value of fluid. For instance, to achieve L f = 0.5 (50% of the volume is melted), an increase of Ra by a factor of 10 leads to a growth of the time by a factor of 1.7. Nonetheless, when Ra is 5 times larger, the necessary time only increases by a factory of 2. This is most likely due to the non-linear intricacies of the problem and requires further investigation. Furthermore, the Nusselt number reported in Figure 8(b) shows that the higher the Rayleigh number, the higher the Nusselt number. This is consistent, since the temperature gradient is integrated along a greater heated wall. Figure 9(a) shows that by increasing the value of δT and consequently increasing the Rayleigh number and the Stefan number, the PCM melts faster. We note that the heigh H of the cavity is kept constant, hence the natural convection flow in the melted PCM is enhanced when the Rayleigh number keep increasing. As a consequence, the convection-dominated regime is reached earlier, as shown by the shift of the minimum of the Nu to lower values of t ϕ in Figure 9(b). This evolution is also observed for the liquid fraction. As expected, an increase of the Rayleigh number and the Stefan number is followed by an enhancement of the heat transfer during the melting, and consequently an improved efficiency of the PCM.

Solidification of the PCM

After the melting of the PCM, we simulate the solidification process. We consider two cases: -case CM: solidification after a Complete Melting of the material (L f = 0.95, Fig. 5f) and -case PM: solidification after a Partial Melting (L f = 0.5, Fig. 5d).

As emphasized previously, the natural convection occurring in the melting PCM is driven by the temperature difference δT = T h -T f . The dimensionless number that decipher the ratio between the forces creating, and those refraining motion, is the Rayleigh number, which appears in the dimensionless form of the Navier-Stokes equations with Boussinesq approximation (section 2). The higher its value, the more intense the heat transfer.

Conversely, during the solidification, the phase-change is handled by the discharged temperature T c , thus the relevant temperature difference in the solid phase of the PCM is T f -T c and the dimensionless temperature in the solid phase should be defined with respect to T f -T c . It is then obvious, see Eqs. ( 7) and ( 11), that the Rayleigh number should be defined using the same temperature difference. However, because the Rayleigh number, as emphasized earlier, amounts for the motion created by hot temperature difference, we choose to keep the same definition for the Rayleigh number as for the melting case, still relevant for the melted core of the flow, where motion still persists, acting as a boundary condition for the solidification process. Under these conditions, in regard with the solidification process, we introduce a new parameter, r δ = (T f -T c )/δT , the normalised temperature is with respect to T f -T c and the relevant Rayleigh number will then r δ × Ra. In the following, we will describe the process of solidification using three different values of r δ . A new scaling is moreover introduced:

V ref = α H ⇒ t = t ϕ ν l H 2 Pr ⇒ Re = 1 Pr . ( 24
)
The solidification stage is indeed a slower process compared to the melting, therefore the use of an adapted scaling is more relevant. This leads to a different time scaling for each cycle.

The simulation of the solidification process starts by imposing at the left-wall a constant (cold) temperature.

The solid phase will propagate into the cavity from both left and right sides, which makes this case computationally challenging. The mesh adaptivity capabilities of our numerical code made possible to accurately track the two solidification fronts identified by the iso-line θ = 0. In the discussion below, the solidification process starts at physical time t ϕ = 185 min (corresponds to τ = 0.2) for case CM and at t ϕ = 59 min (τ = 0.06) for case PM.

Solidification after a complete melting. Case CM

The simulation continues from the state corresponding to Fig. 5 at t ϕ = 185 min (τ = 0.063) and follows after a complete melting. Figure 10 shows the evolution of the PCM during the solidification process. At t ϕ = 185 min (Fig. 10a), the liquid fraction is L f = 0.95 and the melting/solidification front is close to the right wall of the cavity. Setting a low temperature θ c = -1 at the left wall, while the right wall is maintained at a constant temperature (θ right = -0.01 ≤ θ f ) triggers the formation of a second solidification front, propagating from the left side of the domain. Figures 10b and10c illustrate the left part of the cavity solidifying at a faster rate because of the very low temperature imposed at the left wall, inducing a non symmetric evolution of the solid-liquid interfaces. The solid part is represented in blue and corresponds to the region of temperature θ ≤ 0. The signature of the conductive heat transfer is characterized by the vertical shape of the left front. Inside the liquid, the initial convection cells facilitate the heat transfer from the boundaries, resulting in a very rapid decrease of the fluid temperature. Temperature gradients being smoothed out during this first stage, the influence of the convection inside the liquid region is considerably reduced. As a result, the velocity inside the liquid is reduced to very low values. From t ϕ = 430 min (Fig. 10d), the shape of both interfaces is almost symmetrical. This is a signature of a conduction dominated process. At t ϕ = 510 min (Fig. 10e) the liquid region starts to shrink at the bottom sides of the cavity. This process is accelerated and finally the liquid is trapped in a thin pocket and disappears completely from the top of the cavity (Fig. 10e). The complete solidification ends at t ϕ = 530 min, i. e. the liquid fraction is L f = 0. The adapted mesh, refined along the two solidification fronts, at t ϕ = 300 min is reported in Fig. 10f, illustrating the efficiency of the adaptive mesh tool.

Solidification after a partial melting. Case PM

In this case, the solidification starts from the state corresponding to Fig. 11a at t ϕ = 59 min (τ = 0.032), when the liquid fraction is L f = 0.5. The temperature of the left wall is suddenly lower at θ c = -1 as in the previous solidification simulation. The time evolution of the process is illustrated in Figs. 11a-e, while the adapted mesh corresponding to t ϕ = 90 min is plotted in Fig. 11f. As in the previous case, a second solidification front starts to propagate from the left side of the cavity. The straight shape of the left solid front is always observed while the right solid front is impacted by the convection cell present in the middle liquid region (Fig. 11b). The stronger convective effect is most likely due to the huge temperature difference that occurs over smaller space distance (almost half of the volume is occupied by the solid state). This leads to stronger temperature gradients in the liquid region, and then to a stronger transfer. The two fronts merge to form a pocket of fluid which is connected to the top of the cavity (Figure 11c-e). It is interesting to note that, as in the previous solidification case, the left part is solidifying at a faster rate, hence the pocket of melted PCM disappears completely from the right at the top side of the cavity (Figs. 11c-e).

Analysis of the solidification cycle from two different initial conditions: complete and partial melting. Cases CM and PM

The aim of this subsection is to investigate the temporal evolution of some physical properties of the solidification process, from two different initial conditions: i) completely melted volume (case CM) and ii) partially melted volume (50% of the fluid is melted, case PM).

Figure 12 represents the temporal evolution of the liquid fraction, the Nusselt number calculated at the cooled wall defined similarly to (21) but it can be negative in this case and the accumulated heat input Q 0 , for the two investigated cases. Q 0 is defined as follows: where t is a dummy variable.

Q 0 = tϕ 0 Nu dt , (25)
Simulations for three values r δ = 1, r δ = 5 and r δ = 10 are carried out.

Figure 12(a) illustrates the temporal evolution of the liquid fraction L f for CM case. Complete melting occurs for t ϕ = 185 min, after which solidification starts, with a continuous decrease of L f till complete solidification is achieved. For the lowest value of r δ , corresponding to r δ × Ra = 3.27 • 10 5 , the solidification process ends at t ϕ = 530 min. Then, the higher value of r δ , the faster the discharge process, with final times t ϕ = 260 min and t ϕ = 230 min for cases r δ × Ra = 1.62 • 10 6 and r δ × Ra = 3.27 • 10 6 , i.e. a drop of the cold boundary temperature by a factor of 5 and 10 respectively. The solidification speed, quantified by dL f /dt ϕ is nearly constant during almost the whole process for each case. This uniformity of the process indicates that the natural convection flow vanishes during the solidification, and conduction remains the only heat transfer mode.

Figure 12(b) plots the temporal evolution of L f for PM case. As previously discussed, 50% of the volume is melted, at time t ϕ = 59 min, then solidification starts. Furthermore, noticeable is the fact that, despite the solidification process is launched, L F continues to increase slightly at the very beginning of the discharge stage, and then decreases monotonically towards 0 at t ϕ = 240 min. The heat stored in the melted PCM continues to melt the remaining solid PCM until the convection becomes negligible. It is worth noticing that this behavior is not observed in the complete melting case because of the imposed cold temperature in the right wall.

Let us now pay attention to the transfers occurring at the left wall, suddenly submitted to a lower temperature. This is done through the temporal evolution of the Nusselt number, and the temporal-integrated values of the Nusselt number, or the accumulated heat input. Panels (c) and (d) of the Figure 12 illustrate the Nusselt number for the CM and PM. The three investigated Rayleigh numbers are shown, with clear differences between them. This difference corroborates with that already reported for the melting case, over shorter times scales. This indicates that the heat transfer during the solidification process is fundamentally different from the melting one.

For the CM case, for r δ × Ra = 3.27 • 10 5 , the Nusselt number first decreases sharply, for t ϕ ≤ 18 min, then it reaches a plateau at Nu = 7 during the complete melting. At t ϕ = 185 min, solidification starts and Nu suddenly decreases over very short times, reaching negative values (Nu ≈ -15). It follows an increase of Nu with time, up to reaching an asymptotic value very nearly equal to 0 (zero temperature gradients, i.e. uniform temperature at the left wall). The same mechanism is observed over a shorter time steps when r δ × Ra is increased.

For the PM case, the Nusselt number also decreases sharply to a negative values when the solidification starts. However, the convection flow remaining in the melted influences the heat transfer at the very beginning of the solidification process. The hot fluid in the middle of the melted PCM is advected by the natural convection flow to the boundaries and induces a temperature gradient at the left wall, resulting to an oscilating trend of the Nusselt number before reaching the asymptotic value. This is in total agreement with the previous comment about the continuing melting in the right part of the cavity, despite the solidification has started, and the slight increase of the liquid fraction at the very first time steps of the discharging process.

Both charge and discharge cycles are better illustrated in the time evolution of the accumulated heat Q 0 defined in (25), as it is shown in panels (d) and (e) of Fig. 12.

Heat is first stored during the melting stage, corresponding to t ϕ ≤ 185 min for CM (Figure 12(d)) and τ ≤ 59 min for PM (Figure 12(e)), and is then restored during the solidification stage. The CM case indicates higher value of Q 0 (Q 0 = 1400, for r δ × Ra = 3.27 • 10 5 ) compared to the PM case (Q 0 = 500), meaning that PCM is more efficient in terms of heat storage. However, PM case exhibits well balanced characteristic times between the solidification and the melting stages for r δ × Ra = 3.27 • 10 5 . Besides, when the Ra number increases, the stored heat is discharged faster. Moreover, the temperature and the velocity profiles drop sharply during the first step of the cooling process and become almost equal to zero very early in the whole domain. This means that conduction dominates the solidification process, and the convection is rapidly negligible. As a consequence, the melting fronts are vertical and have a symmetric position.

Final discussion and conclusions

The n-octadecane PCM we simulated in this paper is generally used for buildings purposes, due to its phase change temperature, of 28 o C. [START_REF] Zhu | Dynamic characteristics and energy performance of buildings using phase change materials: a review[END_REF] and [START_REF] Kalnaes | Phase change materials and products for building applications: a state-of-theart review and future research opportunities[END_REF] listed various applications, starting from free cooling, peak load shifting, passive building systems and solar energy storage. For each case, PCM is assumed to melt during daytime and to solidify during the nighttime.

We have developed and validated a numerical tool based on adaptive finite elements methods to simulate melting and solidification processes of the PCM. Once our tool was validated against available results, we used it to obtain new physical behaviours of PCM and therefore to predict their ecological use.

It was noticed that, when the same Rayleigh number is used to compare the melting and solidification cycles, meaning that the charge and the discharge modes occur other the same δT , the solidification is always slower than the melting. This behavior is linked to the heat transfer mode leading each of the cycles. Convective heat transfer dominates indeed the melting process, enhancing thus the heat transfer, while conduction is the main heat transfer mode during the solidification, resulting to a slower operating process. However, when the discharge temperature is decreased by a factor of 5, i.e both r δ × Ra and r δ × Ste are increased by a factor of 5, the solidification and the melting occur over similar times.

A first issue that has been brought up by [START_REF] Ascione | Energy refurbishment of existing buildings through the use of phase change materials: Energy savings and indoor comfort in the cooling season[END_REF] is the difficulty of the PCM systems to completely discharge during night-time. Though, if the PCM does not solidify entirely, the effectiveness of the system may be considerably reduced. In this case, it is not advised to melt the PCM completely in order to have a shorter cooling period.

However, for solar energy storage applications, full melting of the PCM is needed to utilize its latent heat storage capacity. Hence, a partial melting is not optimal. For other applications, where shorter discharge time is needed, the use of external cooling techniques is needed in order to ensure a colder discharge temperature.

Figure 1 :

 1 Figure 1: Sketch of the computational domain and boundary conditions. General configuration (panel a) with isothermal (θ = cst.) vertical (x = 0 and x = 1) walls and adiabatic (∂θ/∂n = 0) top and bottom walls. Configuration for the melting phase (panel b) with a hot left wall (θ = θ h > 0) and a cold right wall (θ = θc < 0), followed by a solidification phase (panel c), when the temperature of the left wall is lowered to θ = θc < 0.

Figure 2 :

 2 Figure 2: Location of the interface during the melting of the PCM. (a) Comparison with experimental data of[START_REF] Okada | Analysis of heat transfer during melting from a vertical wall[END_REF] and numerical results of[START_REF] Danaila | A Newton method with adaptive finite elements for solving phase-change problems with natural convection[END_REF] andWang et al. (2010a) for two time instants (τ = 0.032 and 0.063). Ra = 3.27 • 10 5 , Pr = 56.2 and Ste = 0.045 (b) Comparison with both experiment and simulation of[START_REF] Gong | Numerical simulation of melting with natural convection based on lattice boltzmann method and performed with cuda enabled gpu[END_REF] for five time instants (τ = 0.0002, 0.00050, 0.00067, 0.00125, 0.00252). Ra = 2.48 • 10 8 , Pr = 50 and Ste = 0.072.

Figure 3 :Figure 4 :

 34 Figure3: Location of the solid-liquid interface at dimensionless time (panels a to d) τ = 0.0005, τ = 0.002, τ = 0.006 and τ = 0.01, compared with five simulations presented by[START_REF] Bertrand | Melting driven by natural convection a comparison exercise: first results[END_REF]. Ra = 2 • 10 8 , Pr = 50 and Ste = 0.1.

Figure 5 :

 5 Figure 5: Complete melting of the PCM. Temperature iso-lines and streamlines in the fluid phase. The solid part is represented in blue and corresponds to the region of temperature θ ≤ θ f = 0. Time instants (panels a to f): τ = 0.004; 0.016; 0.032; 0.063; 0.08; 0.2. Ra = 3.27 • 10 5 , Pr = 56.2 and Ste = 0.045.

Figure 6 :

 6 Figure 6: Complete melting of the PCM. Time evolution of the average Nusselt number defined at the hot (left) wall (cf. eq. 21) (solid line). Comparison with the experimental results of Okada (1984) (dashed line) and the predictions using the correlation (22) suggested by Jany and Bejan (1988) (dash-dot line). Ra = 3.27 • 10 5 , Pr = 56.2 and Ste = 0.045.

Figure 7 :

 7 Figure 7: Complete melting of the PCM. (a) Time evolution of the liquid fraction for the complete melting of the PCM. (b) Comparison of results (solid line) with the numerical results of Wang et al. (2010a) (dashed line) and the predictions using the correlation (23) suggested by Jany and Bejan (1988) (dash-dot line).

Figure 8 :Figure 9 :

 89 Figure 8: Complete melting of the PCM. Influence of the value of the Rayleigh number (Ra) on the time evolution of the average Nusselt number defined at the hot (left) wall (a) and liquid fraction (b). The reference case (Ra = 3.27 • 10 5 ) is represented by red continuous lines. The value of the Ra was increased by a factor of 5 and 10, respectively while the Stefan number Ste is kept constant.

Figure 10 :

 10 Figure 10: Solidification of the PCM after a complete melting. Temperature iso-lines in the fluid phase. The solid part is represented in blue and corresponds to the region of temperature θ ≤ θ f = 0. Time instants (panels a to e): tϕ = 185 min, tϕ = 231 min, tϕ = 300 min, tϕ = 430 min and tϕ = 510 min. The adapted mesh corresponding to tϕ = 300 min is plotted in panel (f). r δ × Ra = 3.27 • 10 5 .

Figure 11 :

 11 Figure 11: Solidification of the PCM after a partial melting. Temperature iso-lines in the fluid phase. The solid part is represented in blue and corresponds to the region of temperature θ ≤ θ f = 0. Time instants (panels a to e): tϕ = 59 min, tϕ = 70 min, tϕ = 90 min, tϕ = 131 min and tϕ = 200 min. The adapted mesh corresponding to tϕ = 90 min is plotted in panel (f). r δ × Ra = 3.27 • 10 5 .

Figure 12 :

 12 Figure 12: Temporal evolution of the liquid fraction (L f ), the Nusselt number Nu, and the accumulated heat input Q 0 during the entire melting-solidification cycle. Case CM, (left); Case PM, (right).

In the liquid phase, ∇ • (csũ) = h sl ∇ • ũ = 0; in the solid phase, s = 0; in the regularization region, ũ = 0.

Acknowledgements

This project was co-financed by the European Union with the European regional development fund (ERDF, HN0002081) and by the Normandy Regional Council via the M2NUM project. The authors are grateful to CRIANN (Centre Régional Informatique et d'Applications Numériques de Normandie, France) for providing computational resources (project 2017005).