Experimental and Numerical Investigation of Lauric Acid Melting at Suboptimal Inclines

Solidification and melting processes are of interest for numerous engineering applications, including thermal management of electronic equipment, manufacturing processes such as casting, and thermal energy storage (TES). The thermal-hydraulic behavior of solid-to-liquid phase change significantly influences the performance of these systems, and the focus on optimizing the charging and discharging of phase change materials (PCMs) has become prominent in support of TES research. There has been a surge of interest in the numerical modeling of PCMs, which are useful in many engineering applications, but face many difficulties, such as treatment of material properties through phase change, tracking solid–liquid interfaces, and conditions at those interfaces.

Accurately modeling solid-to-liquid phase change is extremely challenging. Despite a plethora of early work, progress has been slow. Research in the modeling of melting and solidification has been an important aspect of TES research for considerable time. In the 1970s, research focused on developing the enthalpy method [1,2] or the effective heat capacity model [3], which are lumped models that do not take into account the separate solid and liquid phases, essentially neglecting the natural convection in the melt. Around that time, investigators also began to realize the importance of the impact of natural convection in the liquid on the modeling of solid-to-liquid phase change processes [4–6]. In 1977, Sparrow et al. [7] were one of the first to perform a numerical solution taking into account both the solid and the fluid motion in the liquid phase. By using the Boussinesq approximation, they numerically simulated the melting of a PCM with a vertical, embedded heater tube. Unlike the prediction provided by the pure-conduction model (enthalpy method), the thickness of the melt layer along the vertical tube varied as a function of height. In their model, many assumptions were made, and experimental validation was not provided. In 1982, Sparrow and Broadbent [8] melted paraffin wax within a cylindrical container and found that by comparing their experimental results with numerical pure-conduction models, the fluid movement had a significant impact on the melt speed.

A number of early investigators used flow-visualization techniques, such as photography [9–14], shadowography [15], interferometry [16], and direct melt-front measurement [17], to study the effect of the natural convection on the phase change process. These studies identified the ideal experimental formats and modes of data reduction for investigating melting experiments, such as the simple cylindrical and rectangular melting and use of the Stefan, Fourier, and Grashof number correlations, and set the precedent for how flow-visualization solid-to-liquid phase change experiments are conducted today. In the pioneering work, Hale and Viskanta [9] optically observed the liquid–solid interface during the melting of n-octadecane from an isothermal vertical wall in a rectangular container. In follow-on work, Ho and Viskanta [15] used a shadowgraph technique to also investigate n-octadecane melting in a similar experimental setup. They also compared their numerical results with experimental data and found excellent agreement for early melt times. For later times, Ho and Viskanta opined that the reason for the discrepancy between the experimental and numerical data had to do with neglecting the volume of expansion during the solid-to-liquid phase change and the conduction along the bottom of the wall in the numerical model.

By the mid-1980s, a variety of geometries and PCMs had been experimentally examined, and different numerical approaches were applied. Transformed grids that track the fluid domain were frequently tested and were able to match experimental results well within certain problem conditions [10,15,18,19]. Methods that tracked the solid–liquid interface using an enthalpy function within a constant mesh were attractive due to their numerical simplicity. Potential methods for approaching these problems were tested and compared by Voller et al. [20]. Treatment of the solid–liquid interface as a porous media was found to be the most appealing due to its simplicity and apparent physicality. This method, which is now referred to as the enthalpy-porosity method, was refined and tested in a stand-alone study [21]. Both transformed grids and enthalpy-porosity were reviewed and compared by Lacroix and Voller [22], who found that the two different methods ultimately performed similarly to each other when comparing melt fronts. They discussed only a few niche situations where one may be favorable over the other, mostly focused on mesh generation and associated computational loads. One note regarding the porosity method was the apparent ambiguity about the porosity coefficient (now commonly referred to as the mushy zone constant) and its lack of a physical interpretation. Ultimately, the enthalpy-porosity model rose in popularity due to its ease of use and integration into commercial software.

One of the first applications of the enthalpy-porosity model was validated using a rectangular flow-visualization experiment based on gallium [17]. The study provided melt-front curves and volume fraction (now typically referred to as melt fraction) for gallium during melting and solidification. The melt curves were used to validate a numerical study by Brent et al. [23] and for the comparative study discussed earlier [22]. The use of gallium data is limiting for other materials due to its high thermal conductivity and viscosity. Lauric acid has been the subject of many modern PCM flow-visualization studies due to its potential as a TES material, similarity to other commercial PCMs, and the ease of visualizing the flow.

Shokouhmand and Kamkari [24] were one of the first to conduct experiments of lauric acid melting within a 50 mm-by-120 mm rectangular container that had an isothermal, vertical wall. The vertical wall was denoted as the $θ=90deg$ orientation (θ defined in Fig. 1). Melt-fraction and melt-front data as a function of time utilizing photographic data collection were provided. This initial study was later expanded upon to include $θ=0deg$ and $45deg$ isothermal melting orientations [25]. The experiments showed a large dependence on both the wall temperature and angle for the melt speed and melt-front shape of the PCM, with a wavy melt front observed due to the formation of multiple convection cells when the apparatus was in the $0deg$ orientation. The study, utilizing the same container, was again expanded through the addition of fins attached to the isothermal surface [26]. These initial studies provided experimental data covering a large range of potential melting configurations for lauric acid; however, the data were limited to only a single container aspect ratio, and configurations where the heated surface was situated above the PCM were excluded (i.e., $90deg<θ≤180deg$⁠).

Numerical simulation of these experiments [24,25] is a common exercise for those seeking to model PCMs due to the well-posed nature of the initial experiments, making them useful validation cases. The experiment for $θ=90deg$ was modeled by Kaheirabadi and Groulx [27] utilizing commercial finite-volume method (FVM) and finite-element method (FEM) softwares, ansys fluent and comsol, respectively, utilizing the enthalpy-porosity model formulated by Voller and Prakash [21]. There were considerable differences in the results between the two methods, and sensitivity of the numerical method to the parameter known as the mushy zone constant was noted. Variations in results due to the mushy zone constant were much larger with FVM compared to FEM. The same group expanded on the study, investigating 0 deg, 45 deg, and 90 deg with FEM software [28]. The additional angles, 0 deg and 45 deg, where complex convection currents form, were not predicted as well by the numerical method compared to the 90 deg case, where the convection currents form a relatively smooth melt front.

Another lauric acid study by Kamkari and Amlashi [29] implemented the enthalpy-porosity model into a FVM code, which presented results that closely matched the experimental data [24,25] for all three incline angles (θ = 0 deg, 45 deg, and 90 deg). The melt fraction over time agreed with the experimental data very well, and the melt front was well predicted; however, the simulation struggled when more complex convection cells were present. Correlations to predict the liquid fraction, energy storage, and time-averaged Nusselt number as a function of a modified Stefan number, Fourier number, Rayleigh number, and angle were developed from the data. The study was further expanded by Karami and Kamkari [30] to inspect the experiments that utilized the fins for heat transfer enhancement. The numerical results also matched the experiment well but struggled again to match the melt front created from the complex currents that formed around the tips of the fins. Many other investigators have studied various heat transfer enhancements for improving solid-to-liquid phase change. Exhaustive reviews on fins [31] and embedded metal foams [32] for enhancing solid-to-liquid phase change are found in the literature.

Additional experimental studies have taken into account different types of boundary conditions. One study, completed by Sathe and Dhoble [33], examined both lauric acid and paraffin wax in a similar rectangular container as shown in Refs. [24,25], melting under a heat flux boundary condition. Avci and Yazici [34] investigated n-eicosane melting at different inclines due to a heat flux condition with a large amount of air in the container, creating an additional fluid interface. Fadl and Eames [35] tested RubiTherm44HC, a commercial PCM, melting from a flux condition on two sides of a rectangular container. RubiTherm25 was also numerically investigated for use within a heat sink by Shatikian et al. [36,37] under both isothermal and constant heat flux boundary conditions.

Numerous studies have also been conducted to examine other container geometries. Ye [38] investigated different aspect ratios numerically, including the effects of the container wall in a conjugate model, which showed that the melting rate increased with the increasing aspect ratio. A more recent study by Duan et al. [39] investigated single-wall isothermally heated containers with various aspect ratios, and a set of correlations for melt fraction and Nusselt number were proposed. Shmueli et al. [40] melted RubiTherm27 within an open-to-air cylinder with isothermal walls. They used their experimental results to validate a numerical model, which involved both a comparison of pressure discretization schemes and a calibration study for the mushy zone constant. Longeon et al. [41] investigated an annular container in a vertical configuration with RubiTherm35 and validated a numerical model utilizing the results. The annular container with RubiTherm35 was also investigated by Siyabi et al. [42], with the addition of different melting inclines. This study also validated a numerical model and found that the 45 deg incline was best for increasing melt speed, and the flowrate of heat transfer fluid did not significantly affect the melt speed. Heat transfer enhancement has also been implemented within annular containers. Karami and Kamkari [43] investigated the addition of various fin configurations in vertical, annular containers melting lauric acid.

Recent experimental studies have employed a larger variety of methods toward the analysis of melting PCMs in rectangular containers. Particle image velocimetry was used to track flow behavior along the melt front of the commercial PolyFin Paraffin within a rectangular container [44]. A two-part study by Azad et al. [45,46] presented experimental results for the onset of convection from an isothermal cylinder for two different PCMs and developed correlations for onset conditions with dimensionless numbers. Thermography, using liquid crystals, has been used to gather more detailed temperature data and, in some cases, provided new information on opaque PCMs. Hu et al. [47,48] used liquid crystals for nano-enhanced PCMs to located the melt front. Marri and Balaji [49] added liquid crystals to visualize the temperature of tetracosane during melting in a rectangular container. The experimental results were used within the same study to investigate and compare different mushy zone constants.

Furthermore, experiments have been designed to directly inspect the mushy zone using simple one-dimensional conditions. Yang et al. [50] used a confocal optical microscope focused on the melt front of paraffin wax and computed an empirical formula for the mushy zone as a function of the local melt fraction and dynamic viscosity. Jiang et al. [51] performed a similar experiment for a binary nitrate salt, adding a row of thermocouples along the melt direction. The study determined a relationship between mushy zone thickness and the temperature of the isothermal wall. The mushy zone constant is a frequent subject of numerical investigations as it has been shown to have a large influence on melt behavior [52–55].

In the present study, the baseline lauric acid investigations [24,25] are extended to include cases where heating is above the PCM (configurations where melting is suboptimal). The work includes, for isothermal wall temperatures of T_w = 60 °C and 70 °C, (1) θ = 180 deg melting experiments, (2) θ = 135 deg melting experiments, (3) numerical modeling of θ = 135 deg melting, and (4) investigations of mushy zone constant and density method interactions. The choice of density method (e.g., setting the density as a function of temperature or employing the Boussinesq approximation) when conducting natural convection simulations can affect the final results, especially when large density gradients exist, as is potentially the case in the melt region with solid-to-liquid phase change. To the best knowledge of the authors, no investigation on the combined effect of the prescribed density method and choice of mushy zone constant when simulating solid-to-liquid phase change using the enthalpy-porosity method have been conducted.

A schematic diagram of the experimental setup is displayed in Fig. 2. The PCM test section, containing lauric acid with 99% purity (Acros Organics, Carlsbad, CA), is made of acrylic with 25 mm-thick walls and 14 mm-thick base, and has interior dimensions of 50 mm × 120 mm × 120 mm (Fig. 3). The container is sealed with a chemical-resistant Viton o-ring and bolted to a 25.4 mm-thick aluminum heat spreader to fully enclose the lauric acid. A flow loop driven by a constant-temperature bath (A40, Anova, Houston, TX) is attached to the heat spreader and designed to provide a controllable, isothermal condition on one wall of the test section. The heat spreader performance was verified by setting the temperature bath to the designated temperature and recording the temperature on the melting surface with five E-type thermocouples calibrated to ±0.1 °C. The thermocouples on the heat spreader showed that there was approximately a 1 °C drop from the bath to heat spreader, and the standard deviation between the five thermocouples was 0.22 °C, thus confirming an isothermal condition.

Lauric acid was melted and poured into the acrylic container. This ensured that the container was full while the PCM was in the liquid state and also limited potential air bubbles. Next, the lauric acid was allowed to solidify, and then the PCM container was bolted to the heat spreader. An initial melt-freeze cycle was completed where all of the PCM in the container was melted and solidified in the 0 deg orientation to ensure good contact with the heat spreader while isolating any voids on the opposite side of the container.

After an initial melt, the test section was oriented to the desired test angle and allowed to reach a steady-state temperature of 23 °C in equilibrium with the ambient. A DSLR CMOS camera (EOS Rebel T6, Canon, Tokyo, Japan) was positioned, and the macro lens was manually focused to ensure limited image variation. The constant-temperature bath was set to the designated setpoint temperature (value dependent on desired isothermal wall temperature) and turned on. The camera was set up using associated software to take images at a frequency of 5 min. When testing the 180 deg melting, the frequency was reduced to 15 min due to the longer length of time required for full melting. The software began to take images at the same time when the bath was turned on, and images taken in the predetermined warm-up period were removed from the data set during processing. After the melt was completed, the container was set to a 0 deg orientation to solidify, which ensured perfect thermal contact for the next experiment.

Experiments were performed for two different angles, 135 deg and 180 deg, at two different isothermal temperature conditions, 60 °C and 70 °C. The angle is defined as the degree from the base position when the isothermal surface is below the PCM, as shown in Fig. 1. The angle was set, with an accuracy of ±1 deg, using an adjustable experimental apparatus and a digital angle gauge. For the angles examined in this study, the heat spreader was above the PCM container. To ensure that there was limited influence from the outside environment, 51 mm-thick insulation was placed on the opposite surface of the heat spreader.

Although some prior experiments embedded an array of thermocouples into the PCM (as many as 32 in a 50 mm-by-120 mm area), they were omitted here to prevent external interference on the melting behavior due to the thermocouple wires.

The cropping and tuner for the conversion to black and white both require calibration for each data set, which presents some opportunity for human influence. To limit this influence, these settings were verified for three to five images before use with the entire data set. Variation in lighting during the experiment was minimized by isolating the experiment from outdoor lighting and keeping artificial lighting consistent.

The numerical solutions were conducted using ansys fluent 19.2. The solid-to-liquid phase change is modeled using the melting and solidification model in fluent, based on the enthalpy-porosity method developed by Voller and Prakash [21]. The pressure–velocity coupling was solved with the SIMPLE scheme, PRESTO! for pressure, and second-order upwind for both momentum and energy. Least squares cell based for gradient was used for spatial discretization. The transient formulation utilized first-order upwind. The residuals for the convergence criteria were set to 1 × 10⁻⁴ for continuity, 1 × 10⁻⁵ for velocity, and 1 × 10⁻⁹ for energy. The under-relaxation factors were set to 0.3 for pressure, 1 for density, 1 for body forces, 0.7 for momentum, 0.9 for local liquid fraction, and 1 for energy.

The numerical model represents an ideal version of the experiment, modeling the interior of the enclosure as a two-dimensional domain with one isothermal side and three insulated sides. A schematic diagram of the solution domain is shown in Fig. 5, where W = 50 mm and L = 120 mm. The domain was set to an initial temperature of T_i = 23 °C. The orientation is set with respect to the gravity vector, as shown in Fig. 1. Material properties of lauric acid were taken from previous studies [28–30] to maintain consistency within the set of investigations. The properties are provided in Table 1.

In the literature, several different methods of dealing with the density when simulating solid-to-liquid phase change exist. As will be shown later, the choice of density treatment for the enthalpy-porosity method can have a significant effect on the results. The preferred method for modeling the density for the present problem is the use of tabulated data dependent on temperature taken from either equations of state, other theoretical analysis, or experiments. However, since detailed tabulated data for the density as a function of temperature is not readily available, four different density approaches, inspired by techniques found in the literature, were examined in the present study and are labeled as follows: (1) standard-variable density, (2) liquid density-based Boussinesq, (3) solid density-based Boussinesq, and (4) standard-variable density with thermal expansion. A graphical representation of the effective density, which essentially replaces the density in the ρg term found in Eq. (3), as a function of temperature for the methods can be seen in Fig. 6. The treatment of the density term in the other parts of the governing equations is discussed alongside each method presented.

This approach was assumed to be used by Fadl and Eames [52,56]. Since, in this method, the density is a function of temperature, it applies to all density terms in the governing equations.

Here, ρ_ref is the reference density, T_o is the operating temperature, and β is the coefficient of thermal expansion. For all other densities appearing in the governing equations, the specified reference density is utilized. Here, the operating temperature is set to the melting temperature, T_liquidus. There are some variations in the literature on what reference density is chosen. Most frequently, the reference density is set to the value of the liquid density [27,28,30,57]. Alternatively, others calculate the reference density based on the average of the solid and liquid densities [58] or based on a solid density reference rather than the liquid density [42]. Although it may not make sense to use the solid-based density value in the Boussinesq method, it has been done before in the literature, and it also provides an opportunity to study the sensitivity of the density value on the results.

An initial mesh-sensitivity study was conducted for $θ=135deg$⁠, T_w = 70 °C, and A_mush = 1 × 10⁵ using the standard-variable density method. A structured grid was created using Pointwise 18.0R3. The mesh was generated for three different refinements of 15,360, 24,000, and 37,500 nodes. The model was run using all three mesh resolutions, and the resultant melt fraction (η) was compared at the same time-step for all three meshes, resulting in less than 1% variance. The grid size of 24,000 nodes was selected for the remaining cases.

In parallel, a time-step sensitivity study was completed, and time-steps of 0.1, 0.05, and 0.02 s were evaluated. The melt fraction at equivalent time-steps for each model were compared, and the variance was found to be less than 1%. The 0.05 s time-step was selected based on these results. Since the density treatment and the value of A_mush were varied, further mesh- and time-sensitivity studies were conducted as needed.

Validation of the model and the mushy zone constant is the subject of the current investigation. The numerical model is directly compared to the experimental melt fraction data in the results and discussed. The models did require some basic melting calibration that is required for most simulations utilizing the enthalpy-porosity model. Starting with the standard-variable density method, A_mush was set to 1 × 10⁵ for $θ=135deg$ and T_w = 60 °C. Initially, the model was simulated with a mushy zone constant of 1 × 10⁵ due to its frequent appearance in similar previous studies, shown in Table 2. It was found that the results with A_mush = 1 × 10⁵ closely matched the experimental data for the T_w = 60 °C case, but underpredicted the melt speed for the T_w = 70 °C case. A set of mushy zone constants from 1 × 10⁴ up to the original 1 × 10⁵ for T_w = 70 °C were evaluated, and the results are shown in Fig. 7. The value of 3 × 10⁴ was found to match the experimental data the closest for T_w = 70 °C. Furthermore, these results show a strong dependence on liquid-phase temperature for selection of the mushy zone constant.

Figure 8 shows the phase change images at several time-steps for the 180 deg experiments at two different wall temperatures (60 °C and 70 °C). The placement of the heated wall above the PCM creates no pathways for convection cells to form, thus the melting process is slow (compared to other inclination angles where natural convection plays a dominant role in the speed of melting) as it appears to be driven solely by conduction through the liquid phase. The melt interface can be generally described as horizontal with a uniform melt thickness, which is expected from the lack of convection cells to create irregular interfaces seen at other melt angles. However, there are noticeable curves near the walls of the container. The additional melting at the wall could be resultant for a few possible reasons. Most likely, there is heat conduction down the side of the container from the heat spreader that aids in the melting of the PCM. Additional experiments that record the wall temperature or that track the liquid-phase movement may be able to better inform the cause of this melt feature.

The experimental melt fraction (η) as a function of time for the 180 deg case is shown in Fig. 9. In general, the rate of melting decreases as the experiment time increases. According to experiments performed by Kamkari et al. [25], the most optimal angle for melting lauric acid in this container was 0 deg. According to the present results, the 180 deg incline melt time was 165% longer than 0 deg at the same 70 °C wall temperature.

To generally inspect the change in melt behavior for θ = 135 deg, the experimental data were analyzed with dimensionless numbers using the form applied by Ho and Viskanta [15], which unifies different Stefan and Rayleigh number cases in relation to the melt fraction. Figure 10 shows the melt fraction as a function of the dimensionless group SteFoRa^0.25 for the 135 deg incline experiment. This generalization is based on the same procedure used for n-octadecane melting within a rectangular container [15]. Recent studies [25,29] have made adjustments to this particular dimensionless group to better fit their unique data, and for the present study, a small adjustment was made to the Rayleigh number in Eq. (16) to include an angle dependence. The vertical line in Fig. 10 marks the shift in melt speed, which corresponds to the time when the solid PCM has fully separated from the isothermal wall. The complete separation of PCM from the isothermal wall marks the transition between regions I and II. In region I, the solid PCM attached to the wall causes a faster rate of melting, and once the PCM separates and region II starts, the melting is slowed due to natural convection heat transfer being the only source continuing to melt the PCM.

As mentioned earlier, the experimental data and numerical model are not always in full agreement. Figure 11 shows both the experimental data and the numerical simulations based on the standard-variable density formulation, utilizing “best-fit” mushy zone constants of 1 × 10⁵ for 60 °C and 3 × 10⁴ for 70 °C (see Fig. 7 for T_w = 70 °C). The first portion of the melt, or region I, is well modeled, and the transition time for slower melting is accurately predicted; but after the first region, the model consistently predicts slower melt speeds than what occurs experimentally. The percent differences in melt fractions between the experimental and numerical data at the transition point are 0.2% for 60 °C at 210 min and 2.1% for 70 °C at 120 min. The maximum percent differences, 8.5% for 60 °C and 11.5% for 70 °C, occur after the transition point.

The experimental melt fraction visualization, numerical temperature contours, and numerical velocity streamlines are shown in Figs. 12 and 13 for wall temperatures of 60 °C and 70 °C, respectively. The most apparent difference between the numerical and experimental results is seen in the melt curves. In prior studies [27–29,68], the numerical melt curves also deviated slightly from the experimental data. The variance from the experimental data in these cases appears to be due to underprediction of the melting near the bottom of the container, which may be attributed to a couple of reasons. First, in the numerical model, the nonisothermal walls are modeled as perfectly adiabatic; and second, it may be possible that the mushy zone modeling method prevents any liquid PCM from circulating into the corner, which can prevent the melt front from progressing into that area.

Region II starts after the solid phase completely separates from the isothermal wall, which causes the rate of melting in the experimental setup to slow. The numerical model shows similar behavior, but the melting slows significantly more than in the experiment. This can possibly be explained by the availability of a larger cold surface of the remaining solid PCM in the experiment to create a stronger natural convection circulation. The melt front, or cold surface, in the experiment has a greater length for a longer period of the melting process. The melt front shown by the numerical model shrinks away from the adiabatic top and the isothermal side walls toward the lowest corner of the container and, generally, is significantly shorter than the corresponding solid surface in the experiment. Interestingly, at 120 min in the 70 °C case (Fig. 13), a secondary recirculation zone near the top forms, rotating opposite to the primary current. In that case, the cold surface of the PCM is no longer able to entrain enough of the liquid PCM to take full advantage of heating from the isothermal wall.

All density methods were tested initially with the same base case, T_w = 70 °C with a mushy zone constant of 3 × 10⁴, which was calibrated for the standard-variable density method. Figure 14 shows the melt fraction versus time for the four different density treatments compared to the experimental data for θ = 135 deg and T_w = 70 °C. It can be seen in the figure that each method has similar melt trends, but the slopes vary slightly. The percent differences between the numerical data and the experiment at the transition point (120 min) are 2.8% for the standard-variable density with thermal expansion, 8.8% for liquid density-based Boussinesq, and 9.9% for solid density-based Boussinesq. After the region transition point, the slopes of the Boussinesq schemes (both liquid and solid based) show faster melting, matching the melting rate of the experimental data more closely. The standard-density method with thermal expansion shows a faster melting rate than the standard method after the region transition but not to the same degree as the Boussinesq methods, which suggests the possibility of a more complex relationship between the melting speed and density. The third variable in this problem, the mushy zone constant, further adjusts the results. It was seen previously in the validation process, in Fig. 7, that the mushy zone constant influences the melt speed. This makes it possible that all four density treatments, which may or may not reflect the true properties of the material, could bring the numerical model within acceptable ranges of experimental data by solely changing the value of the mushy zone constant. The influence of these effects is examined later.

There are additional differences when comparing the temperature contours with velocity vector results at a particular time. The melt fraction and velocity streamlines comparing different density methods for θ = 135 deg and T_w = 70 °C at a time of 60 min (before the transition in melting rate) are shown in Fig. 15. In the figure, the movement of the fluid phase was altered for the variable density with thermal expansion (Fig. 15(b)) and Boussinesq schemes (Figs. 15(c) and 15(d)) compared to the standard-density scheme (Fig. 15(a)). A small fluid circulation forms along the solid-phase wall, rather than the fluid circulating throughout the entire container. Major differences in the melt-front shapes shown in the figure can be attributed to the differences in melt speed.

The temperature contours, also shown in Fig. 15, further reveal a large difference between the standard-variable density method and the other methods. In this figure, the temperature of the liquid region shown is more uniform for the variable density with thermal expansion (Fig. 15(b)) and Boussinesq methods (Figs. 15(c) and 15(d)) compared to the standard-density scheme (Fig. 15(a)). These results reinforce that the density changes in the liquid region play an important role. Without experimental temperature data, investigating the temperature results as a method for discerning one approach as more accurate than another is not possible. However, based on previous experimental data, it can be extrapolated that a more uniform temperature in the liquid region is probable. Within the study by Shokouhmand and Kamkari [24], lauric acid was melted at a 90 deg orientation from a wall at 60 °C. The liquid region above the remaining solid PCM within the container had a uniform temperature matching that of the isothermal wall. This suggests that similar temperature trends could be expected at other orientations, but more experimental data are required to validate this hypothesis.

Figure 16 has been prepared to compare the temperature contours and velocity vectors at 180 min, after the transition point in the melting rate. The standard-variable density method, shown in Fig. 16(a), continues to show a large variation in temperature in the fluid domain. In addition, with a larger fluid region and smaller solid surface to induce the fluid movement, there are several points of recirculation that form. The upper portion of the container circulates up and away from the solid PCM as well as a small pocket at the lower left corner of the melt front, opposite of the isothermal wall. This recirculation of fluid away from the melt front could potentially explain the reduced melting rate of the standard-variable density method after the transition point, as this fluid movement is not contributing to further melting.

In the standard-variable density method with thermal expansion, Fig. 16(b), the fluid movement is now primarily a single circulation throughout the fluid domain and across the melt front. A small amount of opposing circulation can still be seen within the mushy zone. This circulation appears to be warmer fluid in the major current cooling and encountering the wall of the container, becoming a back-flow. This creates a small circulation along the melt front and the lower wall of the container that can be seen turning back into the larger current in the upper-left corner. The weaker back-flow currents along the melt surface and lower wall is where the basic Boussinesq methods, Figs. 16(c) and 16(d), differ from the standard-variable density method with thermal expansion. While there is a more noticeable small current that travels from the melt front along the lower wall of the container in both standard-variable density methods (with and without thermal expansion), this circulation is not as large in the Boussinesq-based methods. This is likely due to the weaker buoyant forces in those regions from the smaller magnitude in density change.

The interaction of the mushy zone constant and density method was investigated using the T_w = 70 °C case for three values of A_mush = 2 × 10⁴, 3 × 10⁴, and 4 × 10⁴. Figures 17–20, which display both the melt fraction and the specific latent heat storage rate (Eq. (15)) as functions of time, and the specific latent heat storage as a function of melt fraction, have been prepared. It can be seen that for each density method, the melting rate varies slightly for different mushy zone constants in the first 80% of the melt, and then reaches a critical point where the melting speed converges to a uniform slope. Interestingly, a melt fraction of 0.8 corresponds to when all of the solid above the horizontal plane has been melted (see Fig. 21) This is best seen in Figs. 17–20(c), which show the melting rate with respect to the melt fraction. For example, Fig. 17(c) shows each mushy zone constant starting at three different melting rates, 2332 W/kg for 2 × 10⁴, 1872 W/kg for 3 × 10⁴, and 1591 W/kg for 4 × 10⁴, and then all converging to around 170 W/kg at a melt fraction of 80%.

Table 3 presents the standard deviation of the initial melt rate between the three mushy zone constants for each density method. The standard-variable density method shows the largest difference in melt behavior for each mushy zone constant, followed by the second-largest, standard-variable density with thermal expansion. In Fig. 17, the latent heat storage rate is very mushy zone constant dependent before 150 min. After 150 min, when the melt fraction exceeds 80%, the latent heat storage rate becomes independent of the mushy zone constant. The experimental data in the figures show more variation than the numerical data, but has similar magnitudes and transition points between the regions.

The method utilizing standard-density with thermal expansion (Fig. 18) shows the fastest melting rates of any other method, and a similar interaction with the mushy zone constant, as seen with the standard-variable density method, is observed. However, the difference between the three mushy zone cases is less than the standard-variable density method. The smaller dependence on the mushy zone constant is likely due to additional convective forces from the density change occurring in the liquid outside the melt range and, therefore, outside of the mushy zone. The convergence of the melting rates occurs at an earlier time than the other three methods, coinciding with the faster melting rates, which leads to reaching the critical melt fraction sooner.

Both the liquid-based and solid-based Boussinesq methods have similar trends, seen in Figs. 19 and 20, including the melt fraction plots being nearly identical. Overall, for these methods, the mushy zone constant has an even smaller impact on the melting rate. The influence of the mushy zone constant appears to be tied to the magnitude of the density change within the melt range. Larger changes in density correspond to stronger convective forces, and if these occur within the mushy zone, then the mushy zone pressure drop has a larger influence on the fluid movement.

The most appropriate method is likely the standard-variable density with thermal expansion as it best approximates the actual physical behavior without detailed temperature-dependent density measurements. To improve on this method, detailed measurements of density over a range of temperatures are required, which may be difficult to complete without specialized equipment, especially over the phase change region. Lack of the required information may be overcome by making sacrifices in the material-property accuracy, which may be necessary with new materials. These changes should be clearly communicated so that the influence can be approximated.

This study presented here shows the restraint that suboptimal inclines place on the melt speed of PCM, which is heavily reliant on natural convection currents formed during heating. In addition, numerical models were shown to predict the initial linear region of the melt fraction well, but failed to predict after the separation from the isothermal wall. The predictions of the numerical model are also very sensitive to the material properties within the melting temperature range, and this influence is scaled with the applied mushy zone constant. The model also lacks in predicting accurate melt interface curves, but in most engineering applications, this is of less consequence. Potential reasons for the shortcomings of the numerical model are discussed and explored through additional numerical test cases.

It is interesting to note that when looking at the melt fraction data with comparison to only experimental data, it is easy to calibrate the mushy zone constant so that the results match well enough for the numerical simulation to be “valid,” even when there may be potential errors with certain aspects of the model, such as material properties or the problem conditions. It is conceivable that a numerical simulation with incorrect material properties (or other issues) could be found to be valid through a simple calibration of the mushy zone constant so that the melt fraction data matches. This is especially concerning when attempting to calibrate a model with experimental data that uses a different material than the subject of the study or without melt-curve data for comparison. Care should be taken to validate numerical models with similar materials and with comparison to different facets of experimental data.

The next steps of the study are to seek higher fidelity material properties measurements, especially within the melt-temperature range, complete experiments with alternative boundary conditions to further examine the numerical model weaknesses, and examine potential changes that can be made to improve upon the modeling methods. Focused studies that couple specific experimental data and numerical models will allow for a thorough understanding of the enthalpy-porosity model so that it can be more accurately applied for a wider array of engineering problems.

The authors gratefully acknowledge Embry-Riddle Aeronautical University for providing computing time on Vega Supercomputer.

There are no conflicts of interest. This article does not include research in which human participants were involved. Informed consent not applicable. This article does not include any research in which animal participants were involved.

The datasets generated and supporting the findings of this article are obtainable from the corresponding author upon reasonable request.

g =

gravity (m/s²)

h =

sensible enthalpy (J/kg)

k =

thermal conductivity (W/(m K))

p =

pressure (Pa)

t =

time (s)

v =

fluid velocity vector (m/s)

A =

area of domain (m²)

H =

enthalpy (J/kg)

J =

partial enthalpy of fusion (J/kg)

L =

height of container (m)

N =

number of pixels

S =

mushy zone source

T =

temperature (K)

W =

width of container (m)

I =

unit tensor

c_p =

specific heat at constant pressure (J/(kg K))

h_sl =

specific latent heat of fusion (J/kg)

q_sl =

specific latent energy storage rate (W/kg)

A_mush =

mushy zone constant

T_liquidus =

melting temperature (K)

T_solidus =

solidification temperature (K)

Fo =

Fourier number

Ra =

Rayleigh number

Ste =

Stefan number

α =

thermal diffusivity (m²/s)

β =

coefficient of thermal expansion (1/K)

γ =

local liquid fraction

ɛ =

numerical divisor

η =

melt fraction

θ =

incline angle (deg)

μ =

dynamic viscosity (kg/(m s))

ν =

kinematic viscosity (m²/s)

ρ =

density (kg/m³)

τ =

stress tensor (Pa)

i =

initial

l =

liquid or at liquidus

m =

melting temperature, averaged

o =

operating

s =

solid or at solidus

w =

wall

ref =

reference