Numerical Study of a Latent Heat Storage System’s Performance as a Function of the Phase Change Material’s Thermal Conductivity

Featured Application

Guidance in the development of phase change materials (PCMs), as well as knowledge towards the design and development of latent heat thermal energy storage (LHTES) systems.

Abstract

The thermal conductivities of most commonly used phase change materials (PCMs) are typically fairly low (in the range of 0.2 to 0.4 W/m·K) and are an important consideration when designing latent heat energy storage systems (LHESSs). Because of that, material scientists have been asking the following question: “by how much does the thermal conductivity of a PCM needs to be increased to positively impact the design and performance of a LHESS?” The answer to this question is not straightforward as the performance of a LHESS depends on the PCM’s thermal conductivity, other PCM thermophysical properties, the type of heat exchange system geometry used, the mode of operation, and the targeted power/energy storage of the LHESS. This paper presents work related to this question through a numerical study based on a simplified 2D model of an experimental setup studied previously in the authors’ laboratory. A model created in COMSOL Multiphysics, based on conduction and accounting for the solid-liquid phase change process, was initially validated against experimental results and then used to study the impact of the PCM’s thermal conductivity (dodecanoic acid) on the discharging power of the LHESS. The results show that even increasing the thermal conductivity of the PCM by a factor of 50 only leads to a maximum instantaneous power increase by a factor of 2 or 3 depending on the LHESS configurations.

1. Introduction

With the ever-increasing emergence of renewable energy technologies in the overall worldwide energy production [1], it is now clear that energy storage solutions are needed to deal with natural fluctuations in the availability of renewables [2]. Knowing that a large percentage of energy end-use is in the form of thermal energy (space heating and cooling, domestic hot water, power generation, thermally driven industrial processes, etc.), it appears logical that energy could be stored thermally in these instances, especially when it has been demonstrated that thermal energy storage (TES) is a less expensive option compared to electrical storage [3]. For example, storing heat in water (or other sensible storage media like rocks) is very inexpensive compared to electrical storage in Li-ion batteries. In contrast, it is very clear that the quality of energy stored in a Li-ion battery is much greater than that stored in a given mass of water (try charging your smartphone from a hot water tank).

That is to say that TES solutions still face many challenges, of which the main one is designing a system that matches the need of a given application at a competitive market price [4]. Latent heat energy storage systems (LHESSs) use phase change materials (PCMs) to store heat with a much greater energy density than sensible heat, through their solid-liquid phase transition [5]. LHESSs also offer the advantage that the phase transition of the PCM happens over a small and well-defined temperature range, allowing for more efficient design of components linked to the LHESS [6,7].

Both the selection of thermophysical properties of the PCM [8,9] and the design of the heat exchange system [4,10,11,12] are central parts in the design of a LHESS. The solid-liquid phase transition temperature often plays the most important role in the selection of a specific PCM for a given application, depending on the operating temperature range of that application. The latent heat and, to a certain extent, specific heat of a PCM are central to its ability to store thermal energy [13].

However, for the design of a LHESS based on transient operation, the (typically) very small thermal conductivity of PCMs (apart from metallic ones and exceptions [14]) are in the order of 0.2 to 0.4 W/m·K [15], leading to limitations on how rapidly heat can be transferred in and out of the storage system, i.e., the “rate problem” [16]. The rate problem is currently more often addressed by designing the heat exchange system used in a LHESS to increase the heat transfer rate to a desired value, for example through the use of additional heat transfer fluid (HTF) passes [5,17], the addition of fins [18,19,20], increasing the contact area [21,22,23], or by favoring close-contact melting during the charging phase [24,25,26]. These approaches aim to optimize the efficiency and performance of LHESSs by overcoming limitations posed by the PCM’s thermal conductivity.

Increasing the PCM’s thermal conductivity could also directly affect that material’s ability to move heat to and from the exchange media, increasing the system’s heat transfer rate [27]. But by how much does the thermal conductivity of a PCM need to be increased to positively impact the design and performance of a LHESS? This question was asked by material scientists from the International Energy Agency (IEA) Annex 40/Annex 67 expert group on Compact Thermal Energy Storage. Optimizing PCM thermal conductivity holds the potential for cost benefits, offering a cheaper solution compared to relying on modifications to heat exchange systems.

This increase in thermal conductivity has been attempted through the addition of conductive nanoparticles [28,29,30]. Numerous issues are impeding the widespread use of those nano-enhanced PCMs (NePCMs), including thermal limitations still stemming from the minimal thermal conductivity increase [31], the cyclability of those systems, the economics of higher NePCM costs, etc. [28]. These same material scientists have tools in their arsenal to design new PCMs with higher thermal conductivities. However, determining the optimal increase in thermal conductivity remains a complex challenge due to the interplay of numerous parameters and properties affecting PCM behavior throughout cyclic operation.

This is a very complex question to answer since the melting and solidification of a PCM during the cyclic operation of a LHESS depend on numerous interrelated parameters and properties, with the importance of the thermal conductivity of the PCM varying throughout the process and for different geometries and operating parameters. The work presented in this paper, therefore, aims at answering this question for a given LHESS that was studied and characterized experimentally in the author’s Laboratory of Applied Multiphase Thermal Engineering (LAMTE) at Dalhousie University [18]. A simplified numerical model, accounting for conductive heat transfer in the PCM as well as the solid-liquid phase change process, was validated using experimentally measured heat transfer rates during the discharging of the system. From this validated model, a numerical parametric study was performed to determine the system’s behavior, in terms of overall heat transfer rates, for increasing thermal conductivity values of a hypothetical PCM, having all the same other thermophysical properties as the PCM used (dodecanoic acid) in the experimental work.

2. Physical Model

2.1. Experimental Setup and Geometry

The experimental apparatus, used in Herbinger and Groulx [18], was made with multiples of the same aluminum finned-tube geometry, as seen in Figure 1a,b. Once assembled, the finned-tube heat exchanger was incorporated in a box containing 10 kg of PCM and measuring 30 cm × 30 cm × 15.25 cm (Figure 1c,d). The fin dimensions are shown in Figure 2.

The experimental apparatus was arranged in 4, 8, and 12 (evenly spaced) fin configurations which were tested at HTF flowrates (water) of 6, 9, and 12 LPM. The charging and discharging of the LHESS were tested at initial and final temperatures varying from 21 to 65 °C. At the start of each experiment, the inlet HTF temperature was ramped up or down from the initial PCM temperature to the desired final PCM temperature, and the system was left to fully charge or discharge to steady state from that point.

2.2. Material Properties

The PCM used in this study was dodecanoic acid, for which all the required thermophysical data are presented in

Table 1. Thermophysical properties of dodecanoic acid [32].

Property	Value
cp,s @ 20 °C	1.95 ± 0.03 kJ/kg∙K
cp,l @ 45 °C	2.4 ± 0.2 kJ/kg∙K
ks @ 30°C	0.160 ± 0.004 W/m∙K
kl @ 50 °C	0.150 ± 0.004 W/m∙K
ρs @ 22 °C	930 ± 20 kg/m3
ρl @ 50 °C	885 ± 20 kg/m3
Tm	43.3 ± 1.5 °C
L	184 ± 9 kJ/kg

. The materials making up the fin were pure copper and 1350 aluminum alloy, with these having defined material profiles in COMSOL Multiphysics.

3. Numerical Model

3.1. Geometry

For the purpose of this numerical study, a 2D sub-domain representing half of a finned tube (and half of the distance between two finned tubes) was modelled in SolidWorks, as shown in Figure 3. The 2D geometry was then imported into COMSOL Multiphysics 6.0 for simulation with appropriate symmetry and insulated boundary conditions.

For the simulation to be representative of the experimental setup and experimental results, particularly in terms of the total mass of PCM in the system, it was assumed that the domain would be 30 cm deep (in the third dimension into the page), while the PCM domain was calculated to be 12.8525 cm high when having a fin spacing of 2.5 cm for 12 fins in the system for the system containing 10 kg of PCM.

Specific to the finned-tube geometry, shown in Figure 2, the thickness of the thin extending aluminum fin was 0.5 mm. Overall, the fin length and PCM depth were kept constant, while the fin spacing was varied to keep the mass of PCM constant in the model. This required calculations to slightly adjust the fin spacing, changing from the theoretical 2.5 cm, 3.75 cm, and 7.5 cm to 2.5 cm, 3.691 cm, and 7.263 cm for 12, 8, and 4 fins, respectively. This change yielded a negligible change in the behavior of the system when compared.

Therefore, the geometry used is shown with the boundary definitions in Figure 3. In the COMSOL model, the entire rectangular domain is meant to be rotated 90 degrees clockwise so that the insulation boundaries (red) are on the top and the bottom; it is presented in this orientation in this paper for space constraints. The symmetry line of the fin and the edge of the PCM domain are defined with a symmetry condition (green), reducing the simulation domain size, and hence the simulation time and file size. The HTF surface (blue) is defined through a convective heat flux using Newton’s law of cooling. The convection coefficient is determined through empirical correlations detailed in Section 3.4, and the HTF temperature is a function of time coming from the experimental inlet temperature data that were ramped up or down at the start of experimental runs, as mentioned previously.

3.2. Heat Transfer Physics

As demonstrated experimentally [33], conduction is the dominant mode of heat transfer during the solidification of PCM (natural convection is present in the cooling liquid PCM during this process but can mostly be neglected [34]). Therefore, COMSOL Multiphysics 6.0 was used to perform the numerical modelling using the heat transfer in solids module to solve the transient conduction equation with temperature (and time) as dependent properties; this equation is pre-defined within the classical transient heat conduction module in COMSOL Multiphysics.

The heat transfer rate from the solidifying PCM to the HTF for the half-fin geometry modelled was calculated using COMSOL’s built in Line Integration function with the expression ht.ntflux, which finds the normal total heat flux crossing the fin surface at every time step (in units of W/m). This heat transfer rate was then multiplied by 2 to obtain the power per fin, then multiplied by the system’s number of fins and the system depth (30 cm) to find the system equivalent power (in W) presented as results in the rest of this paper.

3.3. Energy Conservation and Phase Change Modelling

Energy conservation has been successfully modelled with COMSOL Multiphysics in previous studies by modifying the specific heat capacity of the PCM to simulate the latent heat of fusion ($L$) for phase change once the melting temperature ($T_m$) is reached [35,36]. The phase change effect is based around a piecewise melt fraction function:

$$\varphi \left(T\right)=\left\{\begin{array}{c}0,T<T_m-\Delta T/2\\ \frac{T-(T_m-\mathsf{\Delta }T/2)}{\Delta T},T_m-\Delta T/2<T<T_m+\Delta T/2\\ 1,T>T_m+\Delta T/2\end{array}\right.$$

(1)

where this piecewise function models the state of the PCM as a solid with a value of 0 and liquid with a value of 1 and varies linearly between 0 and 1 over the mushy zone, defined over an interval ΔT around the melting temperature T_m. This function was used to define the change in PCM properties from solid to liquid as can be seen in Equation (3). The latent heat of fusion component of the modified specific heat capacity was created by applying a Gaussian function centered about $T_m$ [37].

$$D\left(T\right)=e^{-\frac{{\left(T-T_m\right)}^2}{{\left(\Delta T/4\right)}^2}}/\sqrt{\pi {\left(\Delta T/4\right)}^2}$$

(2)

The Gaussian function has a value of zero except within the interval of T_m − ΔT/2 to T_m + ΔT/2. The function integrated over the range of all temperatures is equal to 1, which allows D(T) to be multiplied by L to produce an energy balance through the phase change, centered about T_m. Using D(T) combined with the melt fraction function ϕ(T) creates the modified specific heat capacity function:

$$C_p\left(T\right)=c_{p,s}+\varphi \left(T\right)\left(c_{p,l}-c_{p,s}\right)+L·D\left(T\right)$$

(3)

A similar function (excluding the latent heat of fusion) was created for the thermal conductivity to simulate the gradual change in k as the PCM transitions between the solid and liquid phase in the mushy zone:

$$k\left(T\right)=k_s+\varphi \left(T\right)(k_l-k_s)$$

(4)

The density of the PCM was considered constant in the simulation to ensure the conservation of mass and energy, since the volume changes of the PCM were not simulated. The selected density value was the average between the solid and liquid density values: $\rho$ = 907.5 kg/m³.

3.4. Convective Boundary Modelling

In order to more accurately simulate the overall thermal and phase change processes within the PCM, it was imperative to capture as closely as possible the order of magnitude of the convection coefficient related to the HTF on the inside surface of the finned tube. For the various experimental trials used as the basis for this numerical study, the internal flow in the tube was in the transition region, with Reynolds numbers varying from 1000 to 4600, based mainly on flow rates and the number of finned tubes in the parallel configuration and calculated at the final HTF temperature. Therefore, using correlations from Chapter 5 of Rohsenow, et al. [38] defined for 2100 ≤ Re ≤ 10⁶ and the assumption of uniform wall temperature:

$$\mathrm{N}{\mathrm{u}}_t=\mathrm{N}{\mathrm{u}}_0+\frac{\left(0.079{\left(f/2\right)}^{0.5}\mathrm{R}\mathrm{e}\mathrm{P}\mathrm{r}\right)}{{\left(1+{\mathrm{P}\mathrm{r}}^{4/5}\right)}^{5/6}}$$

(5)

with $\mathrm{N}{\mathrm{u}}_0=4.8,\mathrm{R}{\mathrm{e}}_D=\frac{\rho U{D}_h}{\mu },$ and $D_h=\frac{4A_c}{P}$. The required friction factor ($f$) for this Nusselt number correlation is defined as:

$$f=A+\frac{B}{\mathrm{R}{\mathrm{e}}^{\frac{1}{m}}}\left\{\begin{array}{c}\mathrm{Re}<2100;A=0,B=16,m=1\\ 2100<\mathrm{Re}\le 4000;A=0.0054,B=2.3\times {10}^{-8},m=-2/3\\ \mathrm{Re}>4000;A=1.28\times {10}^{-3},B=0.1143,m=3.2154\end{array}\right.$$

(6)

Considering the effect of entrance length on the convection coefficient, it was determined that a weighted average approach was the most physically meaningful and accurate. Equation (7) presents the entrance length correction to be applied to Nu_t (=Nu_∞) determined from Equation (5). The constants C and n were selected for the 90° round bend case presented in Rohsenow, et al. [38].

$$\frac{{\mathrm{N}\mathrm{u}}_m}{\mathrm{N}{\mathrm{u}}_{\infty }}=1+\frac{C}{{\left(x/D_h\right)}^n}\mathrm{w}\mathrm{i}\mathrm{t}\mathrm{h}C=1.0517\mathrm{a}\mathrm{n}\mathrm{d}n=0.629$$

(7)

Finally, since this entrance region correction is dependent on the distance from the 90° bend, a weighted average (${\overline{\mathrm{N}\mathrm{u}}}_m$) was taken. Local corrected Nu_m values were calculated at every centimeter over the 30 cm length of the finned tube, integrated over the last 29 cm using Riemann’s sums, then divided by 29 slices. The average convection coefficient used in a given simulation was calculated from the standard definition of the Nusselt number.

$$\overline{h}={\overline{\mathrm{N}\mathrm{u}}}_m\left(k_{water}/D_h\right)$$

(8)

As mentioned earlier, values of Re for the experimental trials varied from 1000 to 4600, with the correlation in Equation (5) being valid for Re ≥ 2100. It is important to mention that convection coefficients coming from Equation (5) were found to be more accurate and therefore led to simulation results matching experimental validation data better than the purely laminar correlation from Re < 2100. Therefore, the correlation of Equation (5) was used for all numerical studies. The calculated convection coefficients for each validated experiment will be presented in Section 3.7.

3.5. Mesh Independence Study

A mesh independence study was performed to determine the effect of mesh quality on numerical results. Only free-triangular elements were used to create the mesh for these numerical simulations, as shown in Figure 4. The meshing size steps were set differently for the fin and the PCM domains, with the fin being much smaller and not experiencing phase change. Five mesh sizes were used to study the effect of varying the average element sizes in both the fin and PCM domains, as shown in

Table 2. Statistics of the meshes used in the mesh independence study.

Mesh	Domain Elements	Boundary Elements	Total Elements	Average Element Size	Simulation Time
1	5135	888	6023	0.00090 m	17 s
2	8152	985	9137	0.00057 m	16 s
3	15,419	1140	16,559	0.00046 m	23 s
4	30,907	1805	32,712	0.00030 m	46 s
5	130,768	3617	134,385	0.00013 m	194 s

.

Every mesh study simulation had a time step of 30 s and simulation time of 7200 s, using the model with 12 finned tubes, a flow rate of 9 LPM, a discharging temperature range of 65 to 21 °C, and a mushy-zone temperature range of ΔT = 4 K. A defined temperature on the inner side of the copper HTF interface was used at this stage instead of a convective boundary condition and was ramped down using the function made from the experimental data as mentioned previously.

Figure 5 shows the system’s equivalent power; the difference in the studied mesh made a negligible difference. Though negligible, a mesh was selected, as outlined in, Section 3.6 following the mushy-zone range study.

3.6. Impact of Mushy-Zone Temperature Range

The mushy-zone temperature range (ΔT) defines how rapidly the PCM properties change from the liquid to solid values as a function of temperature, and over which range of temperature latent heat is accounted for during the phase change. As presented by C. Kheirabadi and Groulx [37], the mushy-zone temperature range (ΔT) also has an impact on the simulated results when the solid–liquid phase change is accounted for. Previous studies have also shown that the smallest ΔT values lead to more physically accurate results at the expense of increased computational cost [36]. Therefore, simulations were performed to look at the impact of ΔT on the local results and on the extracted power results.

The effect of modifying ΔT was studied regarding its impact on the conductive heat transfer and phase change behavior through a parametric mushy-zone ΔT study. The simulations used mesh 3, as defined in Section 3.5, and were otherwise identical to the simulations ran in Section 3.5 for the mesh independence study.

Figure 6 presents the position of the solidification interface after 7020 s for the four values of ΔT used (1, 4, 8, and 12 K). The results are shown at such a long time after the start of solidification because the impact of the mushy-zone temperature range (ΔT) becomes more noticeable then. It can be observed that the position of the interface, therefore the amount of PCM solidified, varies slightly with values of ΔT. The position of the interface is mostly the same for ΔT of 1 and 4 K, with the difference increasing with an increase in ΔT.

Figure 7 shows the effect of changing ΔT on the system equivalent power curve. In reality, there is no mushy range for the PCM used. Rather, the material goes through its phase transition at a well-defined melting/solidification temperature, and the impact of the latent heat of fusion is felt solely at this transition temperature. Therefore, numerically, increasing ΔT yields a smoother, more continuous power curve with a slightly higher value for the system equivalent power. When ΔT = 1 K, the power curve closer mimics the discontinuous behavior seen in the experimental data at approximately 1250 s.

Also, as ΔT decreases, the modified specific heat capacity function becomes much steeper and taller due to the latent heat correction having to be accounted for within a smaller temperature range. Numerically, this makes the simulation more time consuming and causes convergence issues when varying thermal conductivities individually. Therefore, a ΔT = 2 K was chosen to perform the simulations for this study.

Finally, for the final mesh selection, mesh 3 (from Table 2) was selected, composed of approximately 16,000 elements (for 12 finned tubes). This mesh also eliminated some convergence issues encountered when using coarser meshes with the selected mushy-zone temperature range of ΔT = 2 K.

3.7. Model Validation

The numerical model was validated using experimental data from the trials using 4, 8, and 12 finned tubes, a flow rate of 9 LPM, and discharging temperature ranges of 65 to 21 °C and 55 to 31 °C. The HTF temperature was also ramped down using trial-specific functions made from the experimental data in Excel to mimic each trial as closely as possible. The simulation time was set to 7200 s with a strict maximum time step of 30 s. Figure 8 shows how the numerical results from this simplified geometry and model compare to the experimental results from Herbinger and Groulx [18].

Based on the experimental conditions of each trial, the HTF convection coefficients inside each finned tube were calculated following the approach presented in Section 3.4. These convection coefficient values are presented in

Table 3. Validation study and convection coefficients used.

Validation Study	# of Fins	$\dot{\mathit{V}}$ (L/min)	${\mathit{T}}_{\mathit{r}\mathit{a}\mathit{n}\mathit{g}\mathit{e}}(°\mathbf{C})$	$\overline{\mathit{h}}(\mathbf{W}/{\mathbf{m}}^2\mathbf{K})$
1	12	9	65 to 21	1590.32
2	12	9	55 to 31	1624.05
3	8	9	65 to 21	1953.44
4	8	9	55 to 31	2285.66
5	4	9	65 to 21	3809.42
6	4	9	55 to 31	4180.27

.

From the validation results, some cases showed great agreement; for example, the validation study 1 (see Table 3 and Figure 8a) for 12 fins was able to capture the maximum power and the representative power decline over time. Others captured general orders of magnitude for the system equivalent power without fully capturing the maximum power. It appears clear that this simplified 2D model is not as accurate or predictive for cases using only four finned tubes. However, for the purpose of this study, the validation results from Figure 8 prove that the simplified 2D models capture enough physics from the process, providing general trends. Therefore, the models were found to be valuable for parametric studies on the impact of the variation in the PCM thermal conductivity in this case.

For reference, during this work, simulations were performed using an Intel Core i7-10700k CPU (3.8 GHz with 5.1 GHz Turbo Boost) and 32 GB of RAM.

4. Results and Discussion

The parametric study performed in this paper looked at varying the thermal conductivity of dodecanoic acid by multiplying both the solid and liquid thermal conductivities by factors of 1, 1.05, 1.1, 1.2, 1.5, 2, 5, 10, 20, and 50, and determining the resulting system equivalent power over time. Simulations were performed for 12-, 8-, and 4-finned-tube geometries, with a system HTF flow rate of 9 LPM and over two sets of initial and final system temperatures: T_i = 65 °C to T_f = 21 °C and T_i = 55 °C to T_f = 31 °C.

Firstly, Figure 9 shows the temperature distribution, as well as the solidification interface, at six specific times (t = 1200, 1800, 2400, 3000, 4500, and 6000 s) for the original dodecanoic acid, and PCMs with thermal conductivities of 2k, 10k, and 50k obtained in the eight-finned-tube geometry at a starting temperature of 65 °C and discharged with an HTF at 21 °C. In the figure, the solidification interface is indicated by a white line in each image at the phase change temperature of dodecanoic acid (T_m = 43.3 °C = 316.45 K). From the temperature legend, purple and dark-pink colors correspond to temperatures below the phase change temperature and correspond to the PCM being solid; orange and yellow colors correspond to the PCM being liquid.

As expected, solidification started on the surface of the finned tube and progressed away as time went on, with faster solidification around the tube section of the fin compared to solidification from the fins. Therefore, solidification ended at the top and bottom of the system around the symmetry line.

The progression of the solidification interface for the PCM with a thermal conductivity of 2k was very similar to the original dodecanoic acid, although slightly faster, and a solid bridge was created by the time the simulation reached 6000 s. Increasing the thermal conductivity to 10k, the PCM was fully melted by 6000 s. Further increasing the thermal conductivity to 50k, the PCM was fully melted by 4500 s, but it can also be observed that with this increased thermal conductivity, the fins play a lesser role in conducting heat away from the solidifying PCM; the solidification interface is much more centered on the tube. This points to the fact that an increase in the thermal conductivity of the PCM could allow a change in geometry (the removing of fins in this case), reducing the necessary heat exchange material and potentially reducing the initial cost of the TES infrastructure.

Figure 10 presents the system equivalent power obtained for each finned-tube configuration and operating temperatures for the various PCM thermal conductivity values used in this parametric study. Figure 10a,c,e show the results when the systems were at 65 °C initially and were discharged with an HTF at 21 °C for the 12-, 8-, and 4-finned-tube configurations, respectively. The system with 12 fins always had a larger discharging power compared to the other two. For the original dodecanoic acid (k), the maximum system equivalent power was approximately 500 W for 12 fins, 400 W for 8 fins, and 200 W for 4 fins. Increasing the thermal conductivity to 2k led to a maximum power of approximately 600 W for 12 fins, 450 W for 8 fins, and 260 W for 4 fins. The increase in absolute power was only by approximately 12 to 30%, with a doubling of the PCM thermal conductivity. Finally, increasing the thermal conductivity by a factor of 50 (50k) led to a maximum power of approximately 1075 W for 12 fins, 1050 W for 8 fins, and nearly 800 W for 4 fins. As expected, the increase in thermal conductivity at this point made the number of fins less impactful. An increase in thermal conductivity by a factor of 50 only led to an increase in maximum power by factors of approximately 2.1 to 4. The increase in thermal conductivity was more impactful with fewer finned tubes in the system.

The thermal conductivity increase led to a reduction in the total time needed to fully discharge the system. However, it took an increase in conductivity by a factor of 5 or 10 to bring the discharge time below 2 h, and by a factor of 5 in the case for the four-finned-tube geometry.

Figure 10b,d,f show similar results when the systems were at 55 °C initially and discharged with an HTF at 31 °C. It has been demonstrated that the temperature differential between the HTF driving temperature and the PCM melting temperature plays a crucial role in power extraction [18]. In this case, the HTF temperature being higher led to smaller system equivalent power and a higher difference between power extracted with 12, 8, and 4 fins. For the original dodecanoic acid (k), the maximum power was approximately 375 W for 12 fins, 275 W for 8 fins, and 160 W for 4 fins. For 2k, it was 450 W for 12 fins, 325 W for 8 fins, and 200 W for 4 fins; absolute power increased by approximately 20% to 25% when compared to a standard k. When increased by 50k, the power was 975 W for 12 fins, 825 W for 8 fins, and 600 W for 4 fins.

In this case, the increase in thermal conductivity by a factor of 50 only led to an increase in system equivalent maximum power by factors of approximately 2.6 to 3.8, again with the largest impact of the increase in thermal conductivity being felt with the smaller number of finned tubes.

Finally, the numerical model allows studies where just the solid or liquid thermal conductivity of the PCM is changed (with the thermal conductivity of the other phase staying at the original dodecanoic acid value of k). Figure 11a shows the system equivalent power when the liquid thermal conductivity remains at k and only the solid thermal conductivity is increased.

The power remains the same for each scenario until solidification starts (just before reaching the 0.4 h mark), which is consistent with the liquid thermal conductivity being the same for each scenario. In this scenario, the maximum power extraction reaches a value of 800 W when the solid PCM thermal conductivity is 50k_s, compared to a maximum value of 1075 W when the conductivity of both phases is 50k, as per Figure 10a.

Figure 11b shows the system equivalent power when the solid thermal conductivity remains at k and only the liquid thermal conductivity is increased. This figure shows the impact of the liquid thermal conductivity when the PCM is entirely liquid prior to solidification (again just before reaching the 0.4 h mark). It was also observed that the liquid thermal conductivity increase on its own plays a very minor role in the overall power extracted from the system. From the information obtained through these studies and the results from Figure 11, it is clear that an increase in the solid PCM thermal conductivity plays the most important role overall, and any increase in the thermal conductivity of the liquid PCM then affirms this increase by further increasing the overall system equivalent power.

5. Conclusions

A numerical study, using a simplified and validated 2D numerical model, looking at the impact of the PCM thermal conductivity on the system equivalent power is presented in this work. The results show that with a substantial thermal conductivity increase of 100% (doubling the thermal conductivity), in the range of some of the best results presented when using nanoparticles, for example, the overall system power only increases by approximately 30% at best. From the power vs. time results, it can be concluded that this increase in the PCM thermal conductivity barely affects the overall system discharge history. Increases in the PCM thermal conductivity of at least 500% to 1000% are required to substantially change the discharging power history, with an ultimate power increase up to 400% for a 5000% increase in thermal conductivity.

It is important to keep in mind that the numerical model used for this analysis is simplified, using a 2D approach to capture the important behaviors of a more complex 3D system. Although the validation shows the similarities in the captured behaviors, it is not perfect, so additional processes might play different roles as the PCM thermal conductivity is increased. It is equally important to keep in mind that this analysis is based on a given system, with a given geometry and operating conditions; the impact of the PCM conductivity could be different for other systems. As mentioned in the introduction, the answer to the simple question related to the impact of the PCM thermal conductivity requires complex work because of the dependence of the heat transfer rates during phase change on much more than that one PCM properties. Additional analysis from the community should look at the impact of the PCM thermal conductivity on other systems, as well as the impact that modifying the thermal conductivity might have on other PCM thermophysical properties like the latent heat, density, and specific heat.