Exploring the Limits of Empirical Correlations for the Design of Energy Systems With Complex Fluids: Liquid Sulfur Thermal Energy Storage as a Case Study

New energy technologies are increasingly reliant on the use of fluids with complex properties. Thermal energy storage (TES), in particular, is likely to make use of unconventional working fluids to allow the replacement of heat generated from fossil fuel sources with recycled waste heat from manufacturing processes or accumulated renewable energy that has exceeded existing demand. The integration of TES into energy networks may heighten the reliability and energy efficiency of industrial systems while reducing operational costs and greenhouse gas emissions [1–3]. TES has demonstrated the potential to facilitate more effective use of large-scale energy substitutions that are economically feasible and viable [4], making it particularly attractive as an addition to sustainable energy systems. Concentrating solar power paired with TES, for example, can have an attractive marginal economic value due to its ability to incorporate simple, efficient, and cost-effective energy storage [5–8].

The performance of energy storage systems depends primarily on their storage mechanism and material. Fundamental thermal storage mechanisms and their associated mediums can be classified as either latent storage, where portions of the energy are stored in phase-change materials; thermochemical storage, where portions of the energy are stored in chemical bonds; or sensible processes, where the energy is stored entirely in a fluid that does not undergo phase change or chemical reactions. Latent-based heat storage is a popular approach because the isothermal condition of the phase-change process offers improved thermal properties compared to single-phase heat storage approaches [9]. The thermal performance enhancement of latent storage systems has not yet been fully realized [10–12], as the storage capacity benefit from high latent heat in phase-change materials is currently outweighed by inefficient heat transfer due to low thermal conductivities [13,14]. Thermochemical energy storage offers similar advantages of high energy density, long-term storage, and reduced heat loss [15,16]. Nevertheless, the comparably basic nature of single-phase techniques has made sensible energy storage the most utilized mechanism in commercial TES systems.

In addition to performance characteristics, issues such as storage material availability, storage material cost, structural material degradation, and operational cost often arise in sensible, latent, and thermochemical storage applications [5,17]. Molten salt mixtures are the conventional material choice for sensible TES due to their stable chemistry, low vapor pressure, and high specific heat capacity [18–20]. However, ideal storage materials should also have high thermal stability and widespread accessibility at a feasible cost [21,22]. A major drawback of commonly used 60% NaNO₃/40% KNO₃ molten salts is their limited operational range due to a high freezing temperature of 240 °C and degradation at temperatures above 565 °C [23–25]. Additionally, the high price of molten salts and their lack of availability at the expected scale [17] elevate the cost of TES implementation, hindering its practicality across a broad range of applications. Given these restrictions, extensive research has been devoted to finding alternative materials for TES [26–28].

Recent studies have identified molten sulfur as a promising molten salt substitute [18,24,29–31]. Sulfur is a cheap and abundant material with a low melting temperature of 118.8 °C [18]. As an element, sulfur has negligible degradation with temperature [32]. While the average specific heat of liquid sulfur (1.2 kJ/kg-K [18]) is slightly lower compared to 60% NaNO₃/40% KNO₃ molten salt (1.5 kJ/kg-K [23]), the low cost of sulfur (40–90 $/ton [18]) provides a significant advantage over molten salt (800–1200 $/ton [20]). Furthermore, while molten sulfur has relatively low thermal conductivity, heat transfer is enhanced in sulfur-based TES due to the presence of buoyancy-driven flows [33]. These qualities make sulfur a promising candidate for reliable, efficient, and cost-effective energy storage that will accelerate the commercialization of sustainable energy technologies.

Sulfur’s highly variable viscosity reflects its unusual molecular behavior. Sulfur melts at about 120 °C (393.15 K), forming a low-viscosity liquid mostly consisting of eight-membered ring molecules [36]. When the fluid temperature increases to around 159 °C (∼430 K), sulfur experiences the unusual phenomenon of polymerization in the liquid state [34,36,37]. While almost all properties of liquid sulfur suffer a discontinuity at this temperature, the most striking effect is the sudden change in viscosity [37]. Chain formation causes a spike in viscosity that extends by several orders of magnitude until the effect is reversed by thermal bond dissociation at about 188 °C (∼460 K) [32,38]. Beyond this temperature, the viscosity of sulfur rapidly decreases. It has been found that small amounts of certain chemical substances, such as hydrogen sulfide or polysulfanes, can interfere with the polymerization of sulfur chains to lower the viscosity at temperatures above 158 °C [18,32,38]. Investigations of the effects of H₂S on the properties of sulfur [34,38,39] suggest that this common impurity in nearly all recovered sulfur exhibits uniquely tremendous potential for modification of sulfur viscosity in chemical processing applications [32,40]. Nevertheless, H₂S is highly toxic and flammable [41], making it of great risk to personnel and equipment when present in large enough quantities.

Previous studies of heat transfer in supercritical carbon dioxide suggest that dimensionless scaling of systems with unconventional fluid property behavior is challenging [42,43]. In this work, through the comparison of results from steady-state computational fluid dynamics (CFD) simulations to widely used heat transfer correlations, we determine whether such correlations can remain valid when modeling sulfur natural convection in the range of extreme property variations. We then go on to assess if any correlation can be constructed to model sulfur in this regime or if alternative approaches are needed for the design and analysis of sulfur TES. This assessment is likely to be relevant to a broad range of energy systems using complex fluids.

The model used in this work is a system of steady-state free convection in concentric cylinders, as shown in Fig. 2. The annular fluid is sulfur with 35 ppm H₂S impurities, operating in the temperature range of 400–600 K. The fluid region is bounded by an inner cylinder of radius r_i = 0.02 m and an outer cylinder of radius r_o = 0.1 m. The inner and outer cylinder walls are at constant temperatures T_i and T_o, with cellular flow generated by the temperature gradient across the annular space. The geometry resembles that of a TES system with an uninsulated outer wall maintained at the initial temperature of the fluid, which is not an existing or likely TES configuration. However, the existence of well-established heat transfer correlations for this geometry [44,45] makes it an ideal test case for the predictive ability of dimensionless models. While TES is an unsteady process, steady-state correlations are often useful in quasi-steady approaches, which have been explored and validated for computing design parameters in TES systems [46].

For these correlations, properties are generally evaluated at the film temperature T_f, which is an average of T_i and T_o. This correlation is expected to be accurate for systems of fluid that fall within the bounds of $0.7≲Pr≲6000$ and $Rac≲107$⁠. Studies of the charge and discharge behavior of elemental sulfur in high-temperature thermal energy storage systems have reported Rayleigh numbers on the order of 10⁷–10⁹ for pipe diameters greater than 0.05 m and temperatures greater than 473.15 K [18]. Coupled with the characteristic decrease in 35 ppm H₂S sulfur viscosity beyond a temperature of 480 K, it is expected that the Prandtl and Rayleigh numbers at some of the operating temperatures in this application will not meet the criteria of the standard correlation.

Computational fluid dynamics simulations for this free convection heat transfer model were executed using ansys fluent 2021 r1 [52] software. The 2D concentric cylinder mesh, depicted in Fig. 3, consisted of hexahedral structured cells with refinement around the inner cylinder to accurately capture the near wall effects and symmetrical splitting of the geometry in the XY plane to reduce computational cost. The grid resolution for the computational domain was verified using a series of simulations, where resolution was refined until changes to solutions were negligible. A mesh of ∼25,000 elements with a cell density of ∼170 cells/cm² was deemed sufficient and used to conduct further CFD simulations. The outer cylinder wall was assigned a constant temperature T_o, while the inner wall was assigned a warmer constant temperature T_i. The temperature-dependent thermophysical properties of sulfur with 35 ppm H₂S impurities [24,32,34,35] were implemented into the numerical setup based on the data presented in Fig. 1.

Each simulation imitated a charging scenario, with the initial temperature of the sulfur flow field equal to that of the cold outer wall. The PREssure STaggering Option (PRESTO!) scheme, which is appropriate for buoyancy-driven flows, was used to compute face pressure from cell center values [49]. Pressure-velocity coupling was achieved through the coupled algorithm due to its efficient single-phase implementation for steady-state transport. Each simulation was initially performed in the steady-state setting, with face values for the convective terms in the governing equations obtained using the third-order Monotonic Upstream-Centered Scheme for Conservation Laws (MUSCL) discretization scheme for momentum and the second-order upwind discretization scheme for energy. Residual convergence criteria were at least 10⁻⁶ for continuity and 10⁻⁹ for momentum and energy. If the convergence criteria could not be met with these settings, discretization schemes for momentum and energy were brought down to first order, and the pseudo-transient setting was applied. An additional check for complete convergence was performed by ensuring the energy imbalance between the two cylinders was <0.1%. In temperature regions where sulfur viscosity is especially nonlinear, first-order discretization schemes and employment of the pseudo-transient condition were still not enough to achieve convergence. For cases where this occurred, transient analysis was run using partially converged pseudo-transient data as input. Each transient simulation had a maximum time-step size of 0.1 s. Convergence criteria of 10⁻⁴ for continuity and 10⁻⁶ for momentum and energy were used to ensure numerical accuracy. The simulation was assumed to have reached a steady-state when the net heat transfer rate imbalance was <1%. The heat transfer rates from the inner and outer cylinders were calculated as an area-weighted average across the medium.

The computational model in this study was first verified by simulating the system with air in the annular space, a fluid whose thermophysical properties are standard and well-established. Studies of free convection heat transfer through horizontal cylinders with air or water as the fluid medium often list ±10% variation as the criteria for acceptable accuracy when comparing empirical literature correlations to experimental data or numerical solutions [50,53]. Figure 4 shows comparisons of heat transfer data from CFD simulations to both the Raithby and Hollands [45] and Kuehn and Goldstein [46] correlations for air in horizontal concentric cylinders when T_o = 400 K and 500 K.

When air with an ideal gas equation of state is used as the annular fluid, the models of Raithby and Hollands [45] and Kuehn and Goldstein [46] agree with results from simulations for most cases to within ±10%. Simulation and correlation heat transfer rates agree particularly well when the outer wall temperature is 500 K, with the difference remaining less than ±7.1% for Raithby and Hollands and less than ±5.5% for Kuehn and Goldstein. There is a greater disagreement between CFD results and correlation predictions when the outer wall temperature is 400 K, with differences of up to ±10.13%. Nevertheless, ±10% deviation from the simulations is only exceeded with the Kuehn and Goldstein model and only when the inner wall temperature surpasses 580 K (T* ≥ 0.9). Figure 4 demonstrates that the simulation heat transfer rates for both outer wall temperatures fit between the predictions made by each correlation. The disagreement with simulation results by the Raithby and Hollands and Kuehn and Goldstein models is also consistent with deviations found in comparisons of these correlations to experimental data for air in horizontal circular annuli [50]. Further validation of the computational model is found in Fig. 4(b), which shows the distinct heat transfer curves produced by CFD at each outer wall temperature collapsing into a single curve when plotted as $Nu′Di$ versus $RaDi$ [20]. The results provide sufficient evidence to prove that this computational model is a valid framework for predicting the characteristics of natural convection heat transfer between horizontal concentric cylinders.

Following CFD model verification, simulations were conducted for the concentric cylinder geometry with 35 ppm H₂S sulfur replacing air as the annular fluid. Equations (1)–(11) were performed for outer wall temperatures varying between 400 K and 580 K. Inner wall temperatures were assigned values greater than T_o, up to 600 K. Table 1 provides a full list of the 87 cases that were simulated. As shown in Fig. 1, the most extreme variations of liquid sulfur fluid properties occur in the lower temperature regime, below a temperature of ∼450 K. Comparison of Raithby and Hollands’ [45] and Kuehn and Goldstein’s [46] predictions to results from CFD simulations for outer wall temperatures ≥460 K are shown in Fig. 5.

The results in Fig. 5(a) show that the Kuehn and Goldstein correlation better predicts natural convection heat transfer in sulfur than the Raithby and Hollands correlation. However, neither model can make predictions with an acceptable level of accuracy for the entire temperature range. This is true even when all temperatures are above 460 K, which excludes the most extreme regions of sulfur property variations. Raithby and Hollands’ predictions exhibited a 20−30% overestimation of CFD-generated heat transfer rates for most of the data points at temperatures ≥460 K. Deviations using the Kuehn and Goldstein model from CFD results were within ±10% for 24 of the 32 simulated cases. Variations for almost all the remaining data points remained below ±20%, with the exception of two cases (T_o = 460 K/T_i = 470 K and T_o = 500 K/T_i = 510 K) whose heat transfer rates were both overpredicted by nearly 24%. Dimensionless curves produced by the correlation in Fig. 5(b) are continuous and appear to collapse together quite well across all four outer wall temperatures, indicating the ability of the Kuehn and Goldstein model to remain self-similar throughout this temperature regime. Self-similarity aside, large deviations of the Kuehn and Goldstein model from simulation results are a concern made especially evident through nondimensionalization. Two distinct $Nu′Di$ versus $RaDi$ curves are visible in Fig. 5(b), one associated with correlation data and the other with simulation data. The simulation curve has a steeper slope than the correlation curve, with the most appreciable deviations of the correlation from CFD heat transfer statistics occurring at lower Nusselt and Rayleigh numbers. The data points in this low $Nu′Di$ and $RaDi$ region correspond to cases with lower system temperatures, especially those with T_o = 460 K, 500 K and T_i < 530 K. According to Fig. 1(a), this is the temperature range where sulfur viscosity passes through a local maximum. Such findings suggest that large variations in liquid sulfur viscosity lead to faulty modeling of heat transfer dynamics by existing correlations at low temperatures. Comparisons of correlation predictions to simulation results of heat transfer rates and coefficients for systems incorporating the 400−440 K temperature regime are shown in Fig. 6.

The results demonstrate the extent to which extreme fluid property variations impact the ability of correlations to make accurate predictions of heat transfer. Figure 6(a) illustrates that when the outer wall temperature is dropped down to 440 K, the level of accuracy in heat transfer rates predicted by the correlations resembles that of the higher temperature regime. The Kuehn and Goldstein model underestimated CFD results by a maximum of 21%, with over a third of the deviations from data points remaining within 15%. Coincidentally, the Raithby and Hollands model produced extremely accurate predictions for this outer wall temperature, with a maximum variation of ±5%. Regardless, the lack of consistency in the Raithby and Hollands correlation across temperatures both above and below 440 K provides strong evidence of its invalidity in this application. While the most obvious fluctuations in the fluid properties of 35 ppm H₂S sulfur take place during a temperature increase from 430 K to 450 K, the results show that self-similar predictions of heat transfer can still be made by the literature correlations at temperatures as low as 440 K. Figures 6(b) and 6(c) suggest that it is not until the inner-to-outer wall temperature gradient encompasses the sulfur polymerization transition temperature (∼430 K) that both correlations become completely invalid.

Heat transfer rates for single-phase systems, having a directly proportional relationship to the size of the temperature gradient across the fluid medium, are expected to increase monotonically with ΔT. At outer wall temperatures <440 K, CFD simulations’ heat transfer rate curves increase monotonically but show plateaus across those data points whose temperature gradients largely encompass the sulfur polymerization viscosity transition. Related plots of UA versus T_i reveal drastic decreases in the ability of sulfur to transfer heat once the inner wall temperature surpasses 435 K. These decreases in UA counteract the effect of increasing temperature gradient on the heat transfer rate, resulting in the observed plateaus. Looking at the curves produced from the correlations, dips in UA are predicted to occur far more expeditiously. Consequently, the heat transfer rates predicted by the correlations do not proceed in a monotonic fashion, exhibiting a period of decrease with increasing temperature differences. The resulting peaks seen in the q versus T_i curves of Figs. 6(b) and 6(c) are analogous to the critical point phase transition of heat transfer during boiling, a process with various heat transfer regimes that are notoriously challenging to predict without an assortment of rigid and complex empirical models [54,55].

To further understand the changes in heat transfer and convective flow behavior associated with viscosity variations in the low- and high-temperature regimes, sulfur viscosity contours, temperature contours, and convection currents for outer wall temperatures of 500 K and 400 K are demonstrated in Figs. 7 and 8, respectively.

Empirical correlations for free convection in concentric cylinders operate on the assumption that if T_i > T_o, the cellular flow will follow a pattern of symmetric ascension along the inner cylinder and descension along the outer cylinder [47]. The results in Fig. 7 suggest that buoyancy-driven flow patterns of sulfur at temperatures between 500 K and 600 K exhibit this behavior and maintain a steady semicircular current regardless of the size of the temperature gradient across the fluid region. An increase in temperature difference between the cylinder walls generates more localized flow activity concentrated in the top half of the annular space, causing an increase in temperature nonuniformity throughout the sulfur medium. Still, the general flow path remains consistent and matches expectations of the correlation, thus substantiating the self-similarity of Kuehn and Goldstein’s predictions across higher temperatures. In contrast, Fig. 8 shows that streamlines for an outer wall temperature of 400 K are highly nonuniform and considerably impacted by the size of the temperature gradient across the fluid region. Various patterns of multi-cell formation with different degrees of convective current intensity are evident with increasing temperature gradient, specifically when the gradient enters a new heat transfer regime such as the plateau observed in Fig. 6(c). All the streamlines exhibit a lack of compatibility with the semi-circular flow behavior assumed by the correlations, with localized and limited fluid circulation being apparent even when the inner-to-outer wall temperature gradient is large.

The contours of viscosity in these two examples help to elucidate the heat transfer physics observed in Figs. 5 and 6 and further confirm that property variations are the source of predictive model inaccuracy. As shown in Fig. 8, sulfur flow patterns when T_o = 400 K directly correspond to the distribution of low- and high-viscosity fluid that develops throughout the system geometry. Since the sulfur viscosity initially decreases with increasing temperature, low-viscosity sulfur circulates through the system first. The temperature distribution is such that viscosity becomes lower in the top and bottom of the geometry but remains slightly higher in the middle. When the inner wall temperature is raised above 430 K, the transition from decreasing to increasing viscosity occurs, as well as a spike in specific heat capacity. The sulfur in contact with or near the inner wall now has a higher viscosity than the sulfur that occupies the upper region of the geometry, which has a lower viscosity than the sulfur that occupies the middle region of the geometry. Free convection circulation cells that appear at each temperature gradient directly overlap these transitional viscosity boundaries. The distribution of higher temperatures from the inner wall appears to be hindered, likely in part by the very high specific heat capacity of sulfur at around 430 K. More importantly, as the temperature and viscosity of sulfur lining the inner wall increases, streamlines near the inner wall dissipate. This clearly demonstrates that increases in viscosity have a huge impact on the transfer of heat through the system, resulting in the abnormal heat transfer dynamics observed in Fig. 6(c). When the inner wall temperature exceeds 500 K and the viscosity of sulfur lining the wall begins to decrease, fluid flow near the inner wall is restored. This corresponds to the temperature gradient in Fig. 6(c) at which heat transfer physics returns to increasing proportionally with the temperature gradient.

At an outer wall temperature of 500 K, the viscosity and specific heat capacity of sulfur are past their peaks. As shown in Fig. 7, sulfur at this temperature has better circulation ability, with improved heat transfer capabilities allowing a broader distribution of heated sulfur throughout the system prior to equilibrium. Sulfur viscosity does become less uniformly distributed as it steadily decreases with increasing inner wall temperature, explaining the previously discussed changes in flow behavior and temperature distribution at larger temperature gradients. Still, viscosity transitions are visibly smoother, preventing the random formation of cells that derail convective flow patterns at outer wall temperatures <430 K. These correspondences between convection currents and viscosity contours seen in Figs. 7 and 8 provide clear evidence of the extreme impact of both the low- and high-temperature viscosity transitions on the heat transfer dynamics of the fluid.

Because literature models rely on dimensionless parameters for making heat transfer predictions, which are evaluated using the average fluid properties over a given temperature gradient, large degrees of nonuniformity in property distributions throughout a system can be detrimental to any predictive model’s accuracy. The findings of Figs. 7 and 8 clarify that the nonmonotonic heat transfer rates predicted by the correlations in Figs. 6(b) and 6(c), and the large deviations of correlations from CFD results at both low and high liquid-phase temperatures, are a consequence of the failure of averaged properties to adequately capture the true dynamics of these systems. A more detailed analysis of the correlation predictions in Fig. 6 revealed that the Raithby and Hollands and Kuehn and Goldstein models do a similarly poor job of predicting heat transfer physics in the low-temperature regime. Kuehn and Goldstein was the better model for an outer wall temperature of 400 K, with multiple variations from simulation data points remaining below ±20%. Be that as it may, the viscosity transition at T_i = 430 K caused both correlations to become invalid, and over a third of the data points predicted by Kuehn and Goldstein were underestimated by >30%. For an outer wall temperature of 420 K, Kuehn and Goldstein only performed better than Raithby and Hollands when the inner wall temperature was <520 K. Following the initial correlation break at the viscosity transition temperature, property variations at T_i = 520 K generated multiple roots in the Kuehn and Goldstein solving criteria. The corresponding discontinuity shown in Fig. 6(b) gave rise to underpredictions of simulation results reaching up to 55%. While the Raithby and Hollands model did not experience a discontinuity, nearly three-quarters of its variations from simulation data points exceeded ±20%. The most conclusive evidence of the failure of these correlations can be seen in Fig. 9, which shows the dimensionless representations of heat transfer calculated from simulation results and predictive models at outer wall temperatures ≤440 K.

Figure 9 confirms that the ability of correlations to accurately model heat transfer physics becomes critically hindered by the presence of complex property variations in liquid sulfur. Figure 9(a) shows that at an outer wall temperature of 440 K, while certain deviations of the literature models from simulation results were found to be significant, dimensionless representations of the correlation predictions still exhibit the same general trend as those generated from CFD results. As illustrated in Figs. 9(b) and 9(c), this is not the case when the temperature gradient includes the viscosity transition at 430 K. In both examples, when the inner wall temperature is <440 K and viscosity changes are still very slight, both simulation and correlation curves follow a trend where the Nusselt and Rayleigh numbers increase with increasing temperature gradient. Figures 9(b) and 9(c) show that once the viscosity transition temperature becomes present in the gradient, corresponding to the maximum observed Nusselt numbers of ∼20–25 and the maximum observed Rayleigh numbers of ∼10⁶–10⁷, the curve switches direction and exhibits a decrease in $Nu′Di$ and $RaDi$ with increasing temperature difference. The progression of this curve is opposite to that of the higher temperature cases, which can be explained by the plunging heat transfer coefficients observed in Figs. 6(b) and 6(c). As the temperature gradients continue to increase, the heat transfer dynamics predicted by the correlations share less and less of a resemblance to those produced from CFD data. In particular, the different positive slopes of q and UA seen in Figs. 6(b) and 6(c) for T_i > 500 K transform into an utter loss of consistency between correlation predictions and simulation results upon nondimensionalization. A comparison of the dimensionless curves for the three outer wall temperatures shows the additional loss of self-similarity across overlapping temperature regimes that exist for correlations and simulations alike. Given the effect of complex viscosity and associated fluid property variations on the adequacy of dimensionless parameter predictions in this temperature range, these findings demonstrate the difficulty of using any literature correlation to achieve accurate modeling of liquid sulfur heat transfer.

This study investigated the validity of literature correlations when modeling sulfur heat transfer in the temperature range of extreme property variations, as is applicable to the system design and optimization of sulfur thermal energy storage. A strong understanding of the influence of orders of magnitude changes in viscosity on the heat transfer dynamics of 35 ppm H₂S liquid-phase sulfur was developed and demonstrated using CFD simulations and empirical models. Simulations of sulfur natural convention were performed using a 2D concentric cylinder model with temperature conditions ranging from 400 K to 600 K. Computational heat transfer results were subsequently compared to predictions of well-established literature correlations calculated for each temperature case. At outer wall temperatures above 440 K, the results show inconsistencies in the degree of deviation between CFD results and correlations. Despite a handful of cases exhibiting acceptable agreement between simulations and literature models, large deviations at scattered temperature gradients and differing slopes in q versus T and Nu versus Ra plots give the first signs of correlation failure. At outer wall temperatures below 440 K, the heat transfer ability of sulfur begins to show unconventional fluctuations with respect to temperature gradient. CFD contours and streamlines reveal that in these cases, rapid shifts in sulfur viscosity with increasing temperature cause regions of sulfur with varying thermophysical properties to be introduced into the system. Heat transfer physics and convective flow behaviors of the fluid are significantly altered by these variations, resulting in an unpredictable and highly nonuniform distribution of fluid properties throughout the geometry. The inability of averaged property values to adequately capture the true dynamics of these systems ultimately leads to extreme deviations of correlation predictions from CFD simulation results and a loss of self-similarity in the Nusselt and Rayleigh numbers across systems with different initial conditions.

This analysis shows that literature correlations are of limited utility in modeling liquid sulfur natural convection, as nondimensionalization of heat transfer physics cannot be adequately performed in the presence of complex property variations. While dimensionless models are the most computationally inexpensive way of correlating results from CFD simulations to scalable industrial systems, the data shows that standard definitions of the Nusselt and Rayleigh numbers cannot produce accurate representations of heat transfer for an elementary steady-state case, let alone the transient sulfur TES process. We conclude that achieving rapid optimization of molten sulfur TES, or of any system using a working fluid with highly nonlinear properties, will require a more detailed approach for predicting heat transfer, whether that be utilizing alternative methods to construct case-specific dimensionless correlations or abandoning the use of dimensionless correlations in the design process altogether. While the analysis and conclusions in this study are specific to heat transfer within liquid-phase sulfur, similar challenges can be expected for other fluids with complex thermophysical properties, which are becoming increasingly common in advanced energy systems.

This work was authored by the National Renewable Energy Laboratory, operated by Alliance for Sustainable Energy, LLC, for the U.S. Department of Energy (DOE) under Contract No. DE-AC36-08GO28308. Funding was provided by the Advanced Manufacturing Office’s HPC4EI (High Performance Computing for Energy Innovation) program through the project “Optimization of Sulfur Thermal Energy Storage.” Ms. Oliver’s participation was made possible by a fellowship from the EERE High Performance Computing Summer Internship program. The research was performed using computational resources sponsored by the Department of Energy’s Office of Energy Efficiency and Renewable Energy and located at the National Renewable Energy Laboratory. The views expressed in the article do not necessarily represent the views of the DOE or the U.S. Government.

There are no conflicts of interest.

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