Study on effective front region thickness of PCM in thermal energy storage using a novel semi-theoretical model

Thermal energy storage in mobile applications, particularly battery of electric vehicles, is currently gaining a lot of importance. In this paper, a semi-theoretical time-dependent mathematical model of the phase change in a double shell thermal energy storage module has been developed where the inner tube is a heat exchange surface. An effective front region thickness for the melting and solidification process has been studied. The proposed model is calibrated based on our experimental data. The purpose of such a model is to enable the optimization of the geometry of the energy storage modules in terms of the PCM to the TES container mass ratio and enhancement of phase change rate. In addition, results obtained for a single tube can be used in the bundle of tubes in shell and tube TES design. The results have shown that for the larger diameter of the module (smaller difference between the working tube and shell diameter) the optimal working time is around 2000 s.


Introduction
The importance of thermal energy storage in mobile applications, particularly in battery of electric vehicles, is significantly increasing. Many studies focus on the cooling of lithium-ion batteries as high temperatures can negatively impact the energy storage efficiency and cause the battery failure [1,2]. Our research, however, focuses on constructing a double shell or shell and tube heat exchanger to not only regulate the battery temperature, including preheating in winter, but also to decrease the energy consumption of HVAC systems [3]. This approach opens up the possibility of replacing part of the Li-Ion battery with heat storage using natural and biodegradable materials, which is motivated by both economic and environmental concerns to minimize the use of rare earth elements.
Mobile applications require heat storage to have high capacity while maintaining minimal size and weight. These parameters depend on both PCM (Phase Change Material) selection and heat storage design. With recent developments, many studies are concentrated on various kinds of PCMs as paraffin, esters, alcohols, and acid facts [4]. The cheapest material for thermal energy storage is water. Unfortunately, the relatively low temperature of phase change makes it very often useless in HVAC systems. Taking into account the appropriate phase change temperature the cheapest and most common are raw organic materials i. e., animal and plant-based fats [5]. In this field, coconut oil is easily available in the market in virgin and refined forms. Both of these oils have a good shelf life (2 years for virgin oil and 5 years for refined one) [6]. It is worth noting that in mobile applications this time is much longer than the life of e.g., engine oil. What is more important is that these substances are non-toxic and can be sustainable, and abundant in nature. It should be stressed that coconut palms can absorb and use saltwater during growth [7]. It is a big advantage of that plant due to the limited resources of sweet water in the world.
Due to the low thermal conductivity of PCMs, it's essential to improve the design of thermal storage containers and heat exchange surfaces. There have been various experimental and numerical studies conducted to enhance the efficiency of charging/discharging latent energy storage systems. Pourakabar and Darzi [8] focused on the shape of the module and for the investigation, nine separate cases of multi tubes TES with different fin configurations and same quantity of PCM are taken into consideration. Results indicate that the melting rates are highest and lowest for cases with two tubes arranged vertically and one tube, respectively. The cases with two tubes positioned along the x axis and the case with four tubes arranged in a diamond array had the longest and quickest solidification times, respectively. Rostami et al. [9] studied the use of nanoparticles to increase thermal conductivity in PCMs but the effects on living organisms are not fully understood [10,11]. Youssef et al. [12] suggested a heat exchanger with spiral-wired tubes arranged vertically to enhance natural convection in the container. Although the study showed an even buildup of liquid/solid phase during the melting/ Downloaded from mostwiedzy.pl solidification process, it did not analyze velocity distribution along the tubes. Coconut oil has the advantage of solidifying not only on the heat exchange surface, or the phase change front but also in the volume of the container which is enhanced by natural convection. Furthermore, its abundance in nature, easy availability in marketplace, nontoxic and mostly sustainable behavior make coconut oil as an ideal choice.However, the wide range of solidification temperatures and occurrence of subcooling make it difficult to model it. Many studies have been done on numerical modeling of its melting and solidification process, such as the work of Beust et al. [13] on the influence of the mushy zone constant and phase change temperature by using the enthalpy-porosity model. Which is also could be found in the works of Longeon et al. [14] and Andrassy and Szantho [15]. However, the standard enthalpy porosity model used in these studies does not allow for precise modeling of subcooling.
The specific heat in the phase change region is obtained from DSC (Differential Scanning Calorimetry). The DSC measurement of PCMs requires careful preparation and interpretation of results, as the phase change temperature range is often larger than in macroscale experiments. This offset is a property described in detail in various works [13,16,17] and depends on the heat rate. To reproduce the thermal characteristics of the material [16], the scanning rate and sample mass must be well-controlled, and the sample must be homogeneous. Obtaining material homogeneity in a sample of a mixture of components can be challenging, particularly for samples of just a few milligrams.
In the aspect of the CFD numerical modeling, Günther et al. [18] propose an improvement to model PCM solidification with subcooling. Namely, a new condition has been added to the basic balance equation set, which suppresses phase change in cells that are not in contact with other cells with some solid fraction or with the wall. Thanks to this modification the temperature of PCM can be above the solidification temperature if the cell is not in contact with phase changing front. However, conducted experiments show that in the case of coconut oil phase change may occur in the volume of the container, far from the heat exchange surface or the phase change front. The development of the CFD approach is not the only approach currently underway. Previously, authors presented an analytical model for modeling temperature and heat flux during phase change processes in an experimental TES module [19]. Mazzeo et al. [20] explored an analytic model to solve the Stefan Problem in a PCM layer and the results showed good agreement with experimental results in terms of temperature trends at different layers, surface heat flux trends, and stored energy. The proposed procedure enabled the determination of thermal conductivity and specific heat in the liquid and solid phases, as well as the latent heat. While TES for heating systems at high and intermediate temperatures has been well-researched, little attention has been given to TES for HVAC applications in ambient and sub-ambient temperature ranges. A modular experimental module and test loop were developed to study this issue. The experimental work presented in this paper focuses on obtaining heat transfer from PCM to the volume under constant wall temperature conditions. This paper aims to develop a semi-theoretical time-dependent Table 1 Thermal properties of coconut oil [21].  [21].

Parameter
Operating range Uncertainty mathematical model of the phase change process in coconut oil. Coconut oil was selected as a widely available, biodegradable phase change material with a phase change temperature appropriate for HVAC systems. The proposed model is calibrated on our experimental data. The purpose of such a model is to enable the optimization of the geometry of the energy storage module in terms of the PCM to the TES container mass ratio and enhancement of phase change rate. The double shell design of the TES is proposed as a special case of shell and tube heat exchanger, which is one of the most common designs. Moreover, results obtained for a single tube can be used in a bundle of tubes in shell and tube TES design.

Experimental TES module
In this paper, the mathematical modeling was verified by comparing it with reported experimental data [19]. The described experimental setup, the TES module design, and the procedure is given in detail in reference [19]. The container was made of a PMMA tube with an internal diameter of 40 mm and an external diameter of 50 mm and a height of 165 mm. The flanges were also made from PMMA. A copper tube with a 10 mm outer diameter and a 1 mm wall thickness was used as the heat transfer element. To monitor the temperature, the module was equipped with 15 T-type thermocouples with a 1 mm jacket diameter and 300 mm length, as well as 4 Pt100b class A (4-wire) sensors with a 1.5 mm jacket diameter and 150 mm length. The thermocouples were positioned along expected isothermal lines to reduce measurement errors and were arranged on four cylindrical planes (R7, R10, R13, R16) placed at 60 • intervals, with 5 sensors on each surface (S1 -S5). The sensors were covered at four different radii and the ends were positioned on five surfaces along the axis of the module. Fig. 1 illustrates the sensor arrangement.
Efforts were made to reduce heat loss to the surroundings by insulating all elements of the system and tubing. Only about a quarter of the module casing and a portion of the top flange were left uncovered to allow for visual monitoring of the process using thermal imaging and high-speed cameras.

Thermodynamic properties of PCM
The authors have verified the thermodynamic properties of coconut oil used in the experiments, as reported in the literature showing a wide range of values for the latent heat of coconut oil, from 72 to 249 kJ/kg   4. Dimensionless thickness of phase change with the share of the liquid phase β = 0.5 at surfaces S1 -S5 as function of time for solidification process. Points represent numerical modeling results, linessemi-theoretical approximation. [19]. High-quality cosmetic-grade coconut oil was utilized in the tests. The thermal conductivity and specific heat for the solid phase were obtained using a Tempos Thermal Properties Analyzer (TTPA) at 17.5 • C with an accuracy of ±10%. The melting temperature and latent heat were determined through temperature measurement and energy balance of the experimental module and confirmed with parallel measurements using a differential scanning calorimeter (DSC). The properties used in the present study are summarized in Table 1 [21].
It is noted that the results of measurements performed using DSC differ from those obtained from macroscale measurements. In the current study, the specific heat measured using Tempos Thermal Properties Analyzer (TTPA) was 3450 J/(kg*K) at 17.5 • C, while the DSC measurement showed 7170 J/(kg*K) at the same temperature, 3870 J/ (kg*K) at 10.0 • C, and 2070 J/(kg*K) for the liquid phase at 29.0 • C. The melting temperature determined through visual analysis of the phase distribution in the module and temperature measurement was in good agreement with the DSC measurement, with values of 25.5 • C and 25.0 • C, respectively. The specific heat and fraction of the liquid phase as a function of temperature are shown in Fig. 2 [21].
The results from the DSC measurement for the solidification process of coconut oil are different from the measurements taken on the macroscale. The DSC measurement showed the end of the solidification process to be under 10 • C with a latent heat that does not correspond exactly to the liquid fraction during the phase change process. This phenomenon is related to the intensity of the solidification process, where lower wall temperatures result in more dynamic solidification and larger subcooling. Additionally, coconut oil being a complex mixture, its components have varying solidification temperatures and latent heat, further contributing to the differences in results.
In conclusion, the use of a different measurement method (DSC vs. energy balance of TES module) leads to different results in terms of the thermal properties of coconut oil. The DSC measurement is a dynamic measurement on a small sample, while the value measured based on the energy balance of the TES module was obtained using a large TES module (207 cm 3 ) during a long-lasting steady-state process. This difference in measurement methods should be taken into consideration when interpreting the results and choosing the appropriate method for a given application.
Moreover, uncertainties of all measured parameters are also given Table 2 [19].

Determination of the effective front region thickness
The effective phase change front region thickness "δ" is crucial in determining the efficiency of PCM thermal energy storage systems in a shell-tube design. The thickness depends on various factors such as the PCM's thermal conductivity, the heat transfer coefficient between the PCM and its surroundings, and the rate of heat transfer. To determine the effective front region thickness, numerical/analytical techniques or experiments can be used. Optimizing the thickness leads to improved thermal performance and more efficient energy storage and release.
The effective front region thickness could be referred to the dimension of the PCM layer that experiences changes from a solid to liquid state or vice versa, considering a specified operating time and the unique design of the TES. For PCMs with relatively large phase change temperature range like coconut oil, the effective front region thickness varies in height of the heat transfer surface.
To find the equation for calculating the thickness of the phase transformation front, the following analytical determination was performed.
When the thickness of the PCM layer is small, it can be simplified as a flat partition for thermal resistance calculation purposes. This assumption is based on the dominance of heat flow perpendicular to the heat For this calculation, an alternative thermal conductivity coefficient (λ z ) was used, which is a combination of the average conductivity coefficient of the solid phase of the PCM (near the phase transition point) and a coefficient that depends on free convection.
After using the basics of average heat transfer energy and the definition of rate of heat conduction we can modify Eq. 2 as given below: In conducted investigation dimensionless parameters are employed. The characteristic diameter of the storage module is defined as length available for the PCM: where D s -diameter of the shell, D t -diameter of the inner tube (the heat exchange surface). The Fourier number characterizes transient heat conduction and is defined as: where ttime, α = λ ρ*Cp -thermal diffusivity. Thermal diffusivity was measured as 0.12 mm/s at 17.5 • C using TTPA. Values obtained with above relation are 0.17 mm/s for solid and 0.045 mm/s for liquid.
The Stefan number defines the ratio of sensible heat to latent heat: where ΔT -temperature difference between melting point T m and wall  temperature T w . The Rayleigh number describes the natural convection, in case of natural convection near a vertical wall it is defined as product of Grashof and Prandtl numbers: ν*α (9) where ggravity constant, β -coefficient of thermal expansion, νkinematic viscosity, αthermal diffusivity.
Ultimately, the following equation is obtained from Eq. 5 after some mathematical simplifications and modifications.
Determination of the effective phase change radius cannot be conducted using only the experimental model. Both, in solidification and melting process the radius is not visible. In case of melting, it is obvious. However, during the solidification process the front of phase change cannot be clearly determined due to wide phase change zone. Coconut oil is a complex mixture, for this reason, the hypothesis was assumed that the solid phase visible in the volume of the module, far from the phase change front, may be the effect of solidification of a mixture component of low energy significance. Above a semi theoretical formula for the related thickness is presented. The presented formula is modification of correlation presented previously by A. Bejan and Z. Zhang [22].
The dimensionless thickness δ/L relates to the surface of the isoline for which the share of liquid phase is 0.5. The layer thickness was measured at surfaces S1 -S5 (Fig. 1). The liquid phase equal to 0.5 was chosen because of the following reasons: (1) The phase change front for liquid share 0.0 is growing too slowly in relation to the experiment; (2) The energy share of the partially solidified phase is significant; (3) The isoline for liquid share 0.5 coincides with the fading speed isoline. The last reason is important, because when interpreting the experiment, we can consider the solid as one in which there are no convective movements.
Dimensionless thickness of the solidified surface during the solidification process is presented in Fig. 3. TES is presented horizontally, and x-axis is a heat transfer surface offset from the axis of symmetry by the radius of the inner tube with working fluid.
The function given with the eq. (12) was adjusted to experiment and numerical modeling results [19] by dividing it into two intervals. For the first interval, the coefficients C 1 = 0.14, m 1 = 0.60, n 1 = 0.15, and p 1 = − 0.12 were selected. For the second interval, they are C 2 = 0.35, m 2 = 1.5, n 2 = 0.05, and p 2 = − 0.20. Transition between intervals occurs when both functions give the same result for given time step: Eq. (13) presents a dimensionless cross-sectional area occupied by the solid phase during the solidification process. The formula is given with an integral of the eq. (13). The area under the curve A corresponds to the cross-sectional area occupied by the solid phase after a certain time during the solidification process related to the cross-section of the module A max given as H max D s /2 . The A and A max areas refer to the half of the TES module using the symmetry axis boundary. Also, in this situation if diameter D t of the internal tube increases, A max will remain constant.
Chart in Fig. 4 presents the dimensionless thickness of the solidified layer δ/L within function of time for surfaces S1 -S5. The points correspond to the CFD model calibrated on the basis of the experiment; [19]lines refer to the semi-theoretical model.
Values for layer S2 and S3 are the same for most of the points because of the resolution of the discretization grid, which is 0.5 mm outside the boundary layer. The semi-theoretical lines for those layers are close to each other but are not the same. Fig. 5 presents the visualization of the A/A max dimensionless area.
Eqs. (14)-(16) present a derivation of the formula for dimensionless volume occurred by solidified PCM. What is important, the maximum volume V max refers to the volume of whole module given as H max *π*D 2 s /4 . In this situation if diameter of the internal tube d increases, the volume of PCM decreases but the V max value is constant. Thanks to this, the change of the diameter of the module and the inner tube affects not only the time of the phase change but also the storage capacity.
Where (H) = δ + D t /2. To shorten the expression a variable B is introduced: The shorten formula is as follows: Fig. 6 presents a dimensionless area A/A max and dimensionless volume V/V max as a function of time. As can be seen after 9600 s 0.75 of the surfaces and 0.94 of the module's volume is solidified. This is the limit of the D t /D s = 10/40 configuration. Additionally, the effect of increasing the diameter of the inner tube is shown. In the case of larger inner tube diameter, the solidification process is faster. For inner tube diameter D t = 25 mm the available cross-section surface (0.38) and volume (0.61) is solidified in 2400 s.
Using the above scheme, the melting process was described. Fig. 7 presents phase change isolines using dimensionless thickness of the solidified layer δ/L with the share of the liquid phase β = 0.8 for specific time steps of melting process. The liquid phase share was selected using the same criteria as for solidification process. Fig. 8 presents exemplary phase change isolines for which the share of the solid phase β = 0.8 during the melting process.
The function given with the eq. (13) was adjusted to experiment and numerical modeling results by dividing it into two intervals. For the first interval, the coefficients C 1 = 0.50, m 1 = 0.30, n 1 = 0.13, and p 1 = 0.20 were selected. For the second interval, they are C 2 = 0.10, m 2 = 4.00, n 2 = 0.20, and p 2 = 3.00. The function is not exactly adjusted to points of the dimensionless thickness of the phase change front.
The main difference between solidification and melting process is that liquefied phase relatively fast fills the upper part of the container. This is due to the much larger temperature difference between the wall and the phase change temperature (T liquidus ). As a result, the process is faster. The forming funnel of the liquid phase accelerated the melting process due to developed convection, however it also increased the difference in results based on mean radius values. Fig. 9 presents visualization of A/A max dimensionless area at the cross-section of the TES module.
Chart in Fig. 10 presents a dimensionless area A/A max and dimensionless volume V/V max as a function of time for inner tubes diameters D t = 10 and 25 mm. The shape of curves is different from that of the solidification process. The A/A max and V/V max curves approach asymptotically to their maximum values for specified diameters D s /D t ratios. As can be seen in 10,000 s, 0.64 of the surfaces and 0.76 of the module's volume is melted for D t = 10 mm. For inner tube diameter D t = 25 mm 0.36 of the available cross-section surface and 0.58 of the volume is melted after 6000 s. As for the solidification process the maximum available area and volume ratios for inner tube diameter D t = 25 mm is 0.38 and 0.61 respectively. Those values are not reached in 10,000 s.
The results have shown that for larger diameter of module (smaller difference of working tube diameter and shell diameter) the optimal working time is around 2000s. After that time the melting process slows down. It is possible to melt almost maximum available volume of PCM, however, this requires twice the module runtime. In case of TES with D t = 10 mm the available volume of PCM (possible thermal energy to storage) is visible larger. The results have shown that for smaller diameter of module the optimal working time is around 6000 s. After this time the process slows down. It should be noted that for TES with D t = 10 mm it is possible to reach, like for solidification process, 0.94 of the module's volume. It seems that optimal relation between working tube diameter and the shell diameter is a function of available of charging/ discharging time. For example, if available time is 4000 s the better configuration is TES with D t = 25 mm (0.56 of V max will be melted, compared to 0.35 for TES with D t = 10 mm). If, however the available time will be longer (more than 6000 s) the better option will be TES with D t = 10 mm.

Conclusions
The study shows experimental and analytical research of organic PCM based thermal energy storage module in double shell design. Thermodynamic properties of coconut oil were experimentally obtained by authors.
A semi-theoretical time dependent mathematical model of the phase change in a double shell thermal energy storage has been developed where the inner tube is a heat exchange surface. It was reported that in case of D t = 10 mm, 0.75 of the surface and 0.94 of the module's volume was solidified after 9600 s. While in the case of larger inner tube diameter, the solidification process was faster. For inner tube diameter D t = 25 mm the available cross-section surface (0.38) and volume (0.61) is solidified in 2400 s. On the other hand, during melting process A/A max and V/V max curves approach asymptotically to their maximum values for specified diameters D s /D t ratios. It was concluded that the optimal working time for the module having larger diameter (smaller difference between the working tube and shell diameter) is around 2000s. However, it was further concluded that available charging/discharging time could also be the function for optimal relation between working tube diameter and shell diameter.
Results obtained for this type of TES can be used in shell and tube design. Such technical solutions are often found in heat exchangers due to low production costs and good efficiency. The proposed semitheoretical model can be used to determine the spacing of tubes in the bundle. The model is based on basic physical parameters of the working fluid, which greatly simplifies its application in practice. This module confirmed that optimal relation between of working tube diameter and the shell diameter it is a function of available of discharging/charging time, for the same temperature difference between heat transfer surface and the PCM.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Data availability
Data will be made available on request.