## Abstract

Sustainable energy technologies often use fluids with complex properties. As an example, sulfur is a promising fluid for use in thermal energy storage (TES) systems, with highly nonlinear thermophysical properties. The viscosity of liquid-phase sulfur varies by four orders of magnitude due to polymerization of sulfur rings between 400 K and 500 K, followed by depolymerization of long rigid chains, and a decrease in viscosity, as temperature increases. These properties may compromise the accuracy of long-established empirical correlations in the design of TES systems. This work uses computational fluid dynamics to compute steady-state free convection heat transfer coefficients of sulfur in concentric cylinders at temperatures between 400 K and 600 K. The results show that uneven distributions of high and low-viscosity sulfur in the system cause variations in flow patterns and highly nonlinear heat transfer coefficients as temperature gradients increase. As a result, existing empirical correlations for describing system performance become inaccurate. Comparisons of simulation results to predictions from well-established literature correlations show that deviations may surpass 50%. Nusselt versus Rayleigh number correlations for heat transfer are significantly affected by the loss of self-similarity. The analysis proves that existing correlations are not able to capture the complex properties of sulfur in this temperature range, suggesting that alternative modeling techniques are needed for the design and optimization of sulfur TES systems. These challenges are unlikely to be limited to sulfur as a working fluid or TES but will appear in a range of energy systems.

## Introduction

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_{3}/40% KNO_{3} 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_{3}/40% KNO_{3} 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.

*a*)

*is a nondimensionalized convective heat transfer coefficient*

_{L}*h*based on the length scale of the system

*L*and the fluid thermal conductivity

*k*

*b*)

*is a dimensionless number used to characterize buoyancy-driven flows based on the gravitational acceleration*

_{L}*g*, the thermal expansion coefficient of the fluid

*β*, the thermal diffusivity

*α*(defined as

*k*divided by the product of the density

*ρ*and the specific heat

*c*

_{p}), and the kinematic viscosity

*ν*(defined as the viscosity

*μ*divided by the density

*ρ*)

*c*)

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_{2}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_{2}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.

## Problem and Simulation Setup

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_{2}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].

### Correlation Details.

*a*)

*q*is the heat transfer per unit length. Computation of the effective thermal conductivity,

*k*

_{eff}, provides a correction to the thermal conductivity of a fictitious stationary fluid that yields the same amount of heat that would be transferred by an actual moving fluid [47]

*b*)

_{c}is a Rayleigh number based on the length scale

*c*)

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\u2272Pr\u22726000$ and $Rac\u2272107$. 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^{7}–10^{9} 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_{2}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.

*a*)

*b*)

*c*)

_{D}_{i,conv}represents natural convection heat transfer from a single horizontal cylinder, and Nu

_{D}_{o,conv}represents quasi-steady natural convection heat transfer to a fluid inside a horizontal cylinder [45]. Rayleigh numbers for the inner and outer cylinders per this method are evaluated at a bulk fluid temperature,

*T*

_{b}[49]

*a*)

*b*)

*a*) and (4

*b*) cannot be evaluated until the bulk fluid temperature is known.

*T*

_{b}is calculated based on the following relation:

*T*

_{b}, and convective Nusselt number Eqs. (3

*a*) and (3

*b*) are iteratively solved until the Eq. (5) condition is met. In this study,

*T*

_{b}values are determined using the bisection root-finding method with an error tolerance of 10

^{−7}. The overall Nusselt number found at any Rayleigh number, $Nu\u2032Di$, is then evaluated by combining the overall Nusselt number for natural convection between two cylinders with the Nusselt number for conduction between two cylinders

*a*)

*b*)

*c*)

*UA*, and heat transfer rate,

*q*

*f*at some temperature in the range of

*T*

_{1}–

*T*

_{2}where

*f*(

*T*) equals the average value of that property over the given temperature range. When using the Raithby and Hollands’ [44] correlation, each property is evaluated with the outer cylinder temperature as the lower bound and the inner cylinder temperature as the upper bound. The Kuehn and Goldstein [45] correlation requires two separate evaluations of each fluid property, one with Δ

*T*= (

*T*

_{i}−

*T*

_{b}) for the calculation of Ra

_{Di}(4

*a*) and one with Δ

*T*= (

*T*

_{b}−

*T*

_{o}) for the calculations of Ra

_{D}_{o}(4

*b*) and Pr in (3

*c*).

### Simulation Method.

^{9}at low Prandtl numbers and Rayleigh numbers <10

^{11}at high Prandtl numbers, the flow is expected to be laminar [18,29]. Correlation results suggest that Rayleigh numbers for the cases considered do not increase beyond an order of magnitude of 10

^{8}, so laminar conditions were employed in the 2D computational framework. The governing equations for these simulations include conservation of mass, momentum, and energy

*a*)

*b*)

*c*)

**V**is the fluid velocity vector,

*p*is the pressure,

*t*is the time,

*T*is the temperature, $\tau \xaf\xaf$ is the stress tensor, and

*S*

_{g}is the source term incorporating buoyancy force due to gravity. Determination of

*S*

_{g}is done using the Boussinesq approximation

*ρ*

_{avg}is evaluated at the film temperature (

*T*

_{i}

*+ T*

_{o})/2, and the reference temperature,

*T*

_{ref}, is that of the lower temperature outer wall. While sulfur’s viscosity and specific heat vary widely with temperature, the change in sulfur’s density is approximately linear, suggesting the Boussinesq approximation remains valid within the studied temperature range.

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^{2} 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_{2}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^{−6} for continuity and 10^{−9} 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^{−4} for continuity and 10^{−6} 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.

*Q*, W/m

^{2}) from the numerical simulations is converted to heat transfer rate (

*q*, W/m) and heat transfer coefficient (

*UA*, W/m-K) by the following equations:

*k*is the average thermal conductivity of sulfur for the given temperature gradient.

### Computational Model Verification.

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\u2032Di$ 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.

## Results

Following CFD model verification, simulations were conducted for the concentric cylinder geometry with 35 ppm H_{2}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.

Outer wall temperature T_{o} (K) | Inner wall temperature T_{i} (K) |
---|---|

580 | 590, 600 |

540 | 550, 560, 570, 580, 590, 600 |

500 | 510, 520, 530, 540, 550, 560, 570, 580, 590, 600 |

460 | 470, 480, 490, 500, 510, 520, 530, 540, 550, 560, 570, 580, 590, 600 |

440 | 450, 460, 470, 480, 490, 500, 510, 520, 530, 540, 550, 560, 570, 580, 590, 600 |

420 | 430, 435, 440, 450, 460, 470, 480, 490, 500, 510, 520, 530, 540, 550, 560, 570, 580, 590, 600 |

400 | 410, 420, 430, 440, 450, 460, 470, 480, 490, 500, 510, 520, 530, 540, 550, 560, 570, 580, 590, 600 |

Outer wall temperature T_{o} (K) | Inner wall temperature T_{i} (K) |
---|---|

580 | 590, 600 |

540 | 550, 560, 570, 580, 590, 600 |

500 | 510, 520, 530, 540, 550, 560, 570, 580, 590, 600 |

460 | 470, 480, 490, 500, 510, 520, 530, 540, 550, 560, 570, 580, 590, 600 |

440 | 450, 460, 470, 480, 490, 500, 510, 520, 530, 540, 550, 560, 570, 580, 590, 600 |

420 | 430, 435, 440, 450, 460, 470, 480, 490, 500, 510, 520, 530, 540, 550, 560, 570, 580, 590, 600 |

400 | 410, 420, 430, 440, 450, 460, 470, 480, 490, 500, 510, 520, 530, 540, 550, 560, 570, 580, 590, 600 |

Note: The simulations imitated a charging scenario, with the inner wall always set as the higher temperature boundary.

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\u2032Di$ 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\u2032Di$ 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_{2}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^{6}–10^{7}, the curve switches direction and exhibits a decrease in $Nu\u2032Di$ 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.

## Conclusion

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_{2}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.

## Acknowledgment

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.

## Conflict of Interest

There are no conflicts of interest.

## Data Availability Statement

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

## References

^{1}