© R. Freile and M. Kimber, published by EDP Sciences, 2019

1 Introduction and background

During the last 50 years there has been a rising interest on the use of molten salts as a heat transfer fluid. In particular, they have been considered among the best candidates for the advanced reactor designs of the Generation IV reactor concepts [1].

Molten salts’ primary function in a nuclear reactor is to act as a coolant, extracting the heat resulting from the nuclear fission. These designs may have variations according to how the fuel is arranged in the core. For example, the fuel may be present in the form of a ceramic fuel in a prism or pebble bed core design (Liquid-salt-very-high-temperature reactor) or rather dissolved in the salt itself (i.e. liquid fuel). One very promising design is the liquid-fueled molten salt reactor, which is the only design of the Generation IV concepts which employs a liquid fuel. Fluorides of fissile and/or fertile elements such as UF4, PuF3 and/or ThF4 are combined with carrier salts to form fluids. The most famous carrier salt candidates are FLiNaK (46.5 LiF-11.5 NaF-42 KF mol%) and FLiBe (66 LiF-34 BeF2 mol%).

In a liquid fuel molten salt reactor (MSR) a key factor towards its modeling is the strong coupling between the physics describing the neutronics and the thermal-hydraulics. The velocity field affects the neutron precursors transport and the temperature field affects the neutron population through Doppler effects and density changes in the coolant. Hence, it is of paramount importance to be able to describe the fluid flow and the thermal characteristics of the molten salt fluid. In order to do so, it is required that the turbulence models used are valid for the molten salt flow and that the FLiNaK material properties are accurately described.

Historically, the investigation of molten salts began in the late 1950s at Oak Ridge National Laboratory, aimed at the development of a molten-salt nuclear reactor to power an airplane. In 1954 the first molten salt reactor for the Aircraft Reactor Experiment was built at Oak Ridge to investigate the use of molten fluoride fuels for an aircraft propulsion reactor.

During the same time period, forced convection experiments were performed to examine heat transfer properties of molten salts, in particular FLiNaK. In 1954 Grele and Gedeon [2] and Hoffman and Lones [3] performed experimental measurements using FLiNaK flowing through a heated test section with different piping materials. Nusselt numbers were calculated and contrasted with the Dittus-Boelter correlation, depicting an underprediction of about 50% in value. Final results revealed good agreement when using Ni and SS316, but as much as 50% of underprediction using Inconel. The results obtained by these experimentalists is shown in Figure 1.

Later inthe 1970’s, Vriesema performed measurements on forced convection in a vertical pipe [4]. The results of Nusselt number were 15% off the Dittus Boelter prediction, as shown in Figure 2.

Ignat’ ev et al. [5] also performed heat transfer studies using FLiNaK in a circular tube of iron-based steel. The results were compared with the Petukhov-Kirillov formula, yielding a maximum error of 10% with respect to their experimental data.

Throughout all these experiments appearing in the literature, they all had the common problem of not knowing precisely the value of the thermophysical properties of the molten salt. In particular, the wide range of values adopted by these researchers for thermal conductivity is worth noting. Grele and Gedeon and Hoffman and Loans used for their Nusselt calculations a thermal conductivity of λ = 4.5 W/mK, while Vriesema used λ = 1.3 W/mK. In the Greleand Gedeon report it was mentioned that the reason for the huge discrepancy between their measurements and Dittus Boelter was due to the formation of a film associated with the chromium that was corroded from the Inconel piping, resulting in an additional resistance to heat transfer in the Inconel tube. Studies done by [8] suggest thatanother possible reason for the existent large discrepancy was the adopted value of Flinak thermal conductivity.

Experimentally determining thermal conductivity for molten salts is a challenging task, as it is necessary to differentiate the contributions from the radiation, natural convection and conduction heat transfer modes in a given experiment [6, 8]. Because the thermal conductivity values in the past overestimated by a factor of 4 the most widely accepted values nowadays [7], it can be argued that all the heat transfer modes were not taken into account during the thermal conductivity measurements in that era.

In 1962, Ewing [6] was the first to realize that the thermal radiation in molten salts could significantly affect the experimental measurements for the conductivity. He described the measured value as an effective thermal conductivity, accounting for a molecular conduction plus a thermal radiation conduction. Ewing concluded that the effective thermal conductivity values were very sensitive to the conditions of the experiment. In his measurements, the FLiNaK conductivity values ranged from 0.6 W/m K to 5 W/m K. The former was constant during each run equal to λ = 0.6 W/m K.

Another set of experiments found in literature for measuring the thermal conductivity of FLiNaK was done by Smirnov [7] in 1987. He used the method of coaxial cylinders made of platinum, with radiation heat transfer taken into account. He proposed an empirical correlation for the thermal conductivity equal to λ = 0.36 + 0.00056.T(W/m K) over the range of 517 °C to 817°C.

In 2009, Ambrosek [8] reanalyzed each of the previously mentioned experiments using the most widely accepted thermophysical properties of FLiNaK. In particular, he used Smirnov’s correlation, and found thermal conductivity values of k = 0.81 − 0.93 W/m K over the temperature range of validity. The final result of this analysis agreed within 15% with respect to Dittus Boelter correlation for Inconel piping. The reanalyzed data by Ambrosek is shown in Figure 3. Even though in this figure, the results for Ni and SS316 pipes was not shown, it must be said that they presented a large discrepancy when compared with the Dittus Boelter correlation. No explanation could be found for these results.

After this work done by Ambrosek, it was concluded that the main reason why the first experimentalists had observed results in great discordance with the known correlations was the misprediction in the thermal conductivity of FLiNaK used for their calculations.

In 2013, Idaho National Laboratory published an Engineering database of liquid salt thermophysical properties [9]. It consists of acompilation of the different commonly accepted correlations for the temperature dependent FLiNaK properties such as specific heat, density, viscosity and thermal conductivity. In Table 1 some of the most widely accepted correlations for these properties are shown. Note that for the thermal conductivity there is only one expression in the table, which is the most widely accepted temperature dependent expression. However, according to studies done by Williams et al. [12], the thermal conductivity value should range between 0.6 W/m K and 1.0 W/m K at 973 K.

As noted earlier, it is unclear which set of expressions would be best to use. With this in mind, the authors in [9] determined the effect of uncertainty in the resulting Nusselt number, but assumed complete validity of the Dittus Boelter equation. For the worst case scenario, the maximum reported error due to thermophysical uncertainties in the Nusselt number was + ∕−8%. However, those authors suggested a more detailed analysis should be conducted, one that does not rely on Nusselt number correlations or the assumption of constant properties.

In the scientific community other efforts were made to perform sensitivity analysis of heat transfer in molten salts. Xu et al. [10] conducted a sensitivity study of heat transfer mechanisms in a packed-bed molten salt thermocline thermal storage system using a two-phase model. Different interstitial heat transfer correlations, which are functions of the Reynolds number and the porosity, were used to compare the thermal performance of the storage system. Sabharwall et al. [11] performed hand-made calculations of the sensitivity coefficients for the different thermophysical properties in the VHTR (very high temperature reactor). In their calculations, an analytical equation for the heat transfer coefficient was obtained by using the Dittus Boelter equation and using a fixed pressure drop with friction coefficients for the velocity estimation, which was made using the Blasius formula. Partial derivatives of the heat transfer coefficient with respect to each individual thermophysical property were taken and dimensionless sensitivity coefficients were obtained. Therefore, in the calculation of the sensitivity coefficients, they assumed that Dittus Boelter is valid and that the bulk temperature of the salt does not strongly depend on the molten salt properties.

To the best of our knowledge, there is no definitive expression that best defines molten salt material properties, and more importantly,a thorough investigation is needed to truly quantify the impact uncertainty in material properties has on the pertinent thermal hydraulics. Several expressions for each material property are found in the literature, with noticeable differences in the values one would predict. In the present work, computational fluid dynamics (CFD) was used to obtain sensitivity coefficients for each material property. With this information, the effect of the uncertainties in the material properties on the heat transfer coefficient was accounted for. The CFD simulations using RANS models were compared to a published experiment on molten salt done by Grele and Gedeon [2] and the adequacy of different Nusselt number correlations was assessed.

Fig. 1 Experimental results obtained by Grele and Gedeon [2] and Hoffman and Loans [3] compared with Dittus Boelter correlation.

Fig. 2 Experimental results obtained by Vriesema [4] compared with Dittus Boelter.

Fig. 3 Grele and Gedeon, Hoffman and Loans and Vriesema results for Inconel piping reanalyzed by Ambrosek et al. [8].

2 Model

2.1 Hydrodynamics

In this work, equations of flow motion were solved using the Reynolds Averaged Navier Stokes models. Using this approach, the instantaneous velocity is expressed as the sum of a mean velocity , and a fluctuating part, , such that: (1)

where

Replacing the instantaneous velocity in the Navier Stokes equations, and time averaging both yields the Reynolds averaged equations ofmotion: (2) (3)

It must be clarified that the compressibility effect (change in the density of the fluid) can be caused either by temperature or by pressure. In the case of molten salt flow in a pipe, from the standpoint of the pressure influence on the density, it can be said that it will remain constant. The way to justify this statement is by analyzing the Mach number, which in this particular case is much lower than one. On the other hand, in heated flows, there will be a density variation due to its dependence on temperature. Although the dilatation term was modeled, it can be said that its value will still be negligible compared with the rate of strain tensor.

The last term on the right hand side of equation (3) represents the effects of turbulence and it is commonly known as the Reynolds stress tensor. Note that for this Reynolds stresses an overbar is used to denote the time averaged quantities of those variables. In this work, the Boussinesq approximation was used. The expression for the Reynolds tensor using this approximation is shown in the next equation: (4)

where k is defined as the turbulent kinetic energy and μ_t is the turbulent viscosity, whose expression varies accordingly to the model used.

In this work a two-equation model was used. In these types of models the turbulent viscosity μ_t is taken as a function of the turbulent kinetic energy and the turbulent dissipation . Particularly for this work, the k-ω SST model was used, where is called specific dissipation or turbulent frequency.

The shear-stress transport (SST) k-ω model was developed by Menter [15] with the objective of combining two of the most popular two-equation models, k-ω and k-ɛ, through the use of blending functions. Both of the aforementioned models are multiplied by a blending function and both models are added together. The blending function is generally a hyperbolic tangent whose range is restricted between zero and one. It is designed in such a way that near the wall, the standard k-ω model is predominant, and zero away from the surface, which activates the k-ɛ model in the free stream region. By doing this, the k-ω SST model is directly usable all the way down to the wall through the viscous sub-layer, hence being able to resolve turbulent parameters up to the wall region. In addition, in the free stream region it avoids the common k-ω problem of being overly sensitive to the inlet free stream turbulence properties and takes advantage of the k-ɛ model’s advantage for free-shear flows.

There are several works which have been done on the CFD RANS modeling of molten salts. Ferng et al. [20] proposed a CFD methodology for investigating thermal-hydraulic characteristics of FLiNaK in a pipe geometry using k − ɛ model. Results were compared to existing correlations. Chen et al. [19] performed calculations of Nusselt numbers for Hitec molten salt using four different RANS models and compared them with a present experiment of Hitec and several other experiments with different salts.

This turbulence model was chosen because of the good performance it has shown in other simulations related to the geometry of interest (i.e., pipe flow or channel flow). Kim et al. [16] investigated the effects of non-uniformity of fluid properties on forced convection heat transfer using different turbulence models by comparing against experiments. In this study, the effect of the non-uniformity of fluid properties was accounted for by applying a factor to the value of the Nusselt number for constant properties. The prediction with the k − ω model formulation was in good agreement with the experimental results. Menter et al. [17] compared the performance of the k − ω SST model in heat transfer applications with other turbulence models such as Low-Re models and k-epsilon models against experiments for different geometries. The best overall performance was achieved with the k − ω SST model. Given the conclusions of these two studies, the k − ω SST model is employed in the current work as well.

2.2 Heat transfer

Similarly to what was shown in equation 1, the energy E and enthalpy H are also decomposed into a mean and a fluctuating value. In order to compute the temperature distribution in the fluid (Eq. (5)) and in the solid (Eq. (6)), the equations of energy conservation were used: (5) (6)

where , and is defined as a volumetric heat source in the solid. The term is commonly known as the turbulent heat flux.

In order to model this last term, a similar concept to the one used in Boussinesq approximation (Eq. (4)) is used: (7)

where α_t is the eddy-diffusivity and Pr_t is the turbulent Prandtl number. This number was given a generally used constant value of 0.85 for the simulations [15].

Regarding radiation effects, results provided by Ambrosek et al. [8] proved that, because of the transparency of FLiNaK [6], the amount of energy transferred by radiation can be significant in applications involving high temperatures (T = 1123 K) and laminar flow conditions (Re < 500) in pipes with a diameter of 1 cm or greater. In the current experimental configuration, the pipe diameter is too small to take into account radiation heat transfer, as the convective heat transfer coefficient prevails over the radiative one. This was also proven by Ethan S. Chaleff et al. [22].

3 Methodology

3.1 Computational domain

RANS simulations were carried out in this study using ANSYS Fluent v19 commercial package. For this work the k-ω SST model was used for reasons previously noted. CFD codes require generating a computational domain and mesh to simulate flow behaviors in a specific geometry. The computational domain modeling the experimental test section in Grele and Gedeon experiment [2] was completed using ICEM CFD 19.0. The computational domain is shown in Figure 4.

The model consisted of a solid volume modeling the Inconel X pipe wall and a fluid volume to represent the molten salt flow. The diameters, the pipe wall thickness and the test section’s length dimensions were in agreement with the test section used in [2]. Specifically, it had an outside diameter of 3/8 inches and a wall thickness of 0.065 inches. The total length was 24 inches.

A mesh convergence study using Richardson extrapolation [18] was made to assure that the model was producing a mathematically accurate solution. The number of elements was gradually increased until the wall temperature error for the highest Reynolds number casewas below 1E-2 K. The final total number of elements in the mesh was 1,054,585. The minimum orthogonal quality was 5.14e−1 and the maximum aspect ratio was 50.4 in the cells which are the nearest to the wall. As seen in Figure 4, in the radial direction a 30 × 30 uniform square cartesian grid was used at the center of the pipe. The vertices of the square were equidistant from the origin at a distance of 1.8 mm. Thirty additional radial layers of cells exist between the bounds of the square and the inner diameter of the pipe. A growth ratio of 1.005 from the wall was imposed in this section, and the first mesh division was set at 2E-6 m. This refinement in the fluid region near the wall was made so that the model could resolve the equations for the temperatures and velocities up to the viscous sub-layer. In internal flows, the distance from the wall is generally represented by using the y⁺ coordinate, which is defined as , where y is the distance from the wall and τ_w is the wall shear stress. For the highest Reynolds case, the mesh contained at least one grid point below normalized wall distance, i.e y⁺ =1, and at least five grid points in the viscous sublayer region, i.e. y⁺ = 5. In the solid region, 8 equidistant divisions were imposed radially. In the axial direction, two hundred uniformly spaced divisions were imposed.

The numerical solution was obtained using the SIMPLE algorithm [21]. The gradients were computed by using the Green-Gauss Cell-Based method, which is taken from an average of the values at the neighboring cell centers in a certain direction. For the convection terms, face values of different variables are needed. For the momentum, turbulent kinetic energy, dissipation and energy equation a second order upwind scheme was used, which takes information from the cell-center value and gradient value in the upstream cell. The simulations were considered converged once the residuals levels became lower than 1E-8.

Fig. 4 Perpendicular section of the computational domain used in the simulations.

3.2 Boundary conditions

The boundary conditions applied to the numerical model include standard hydrodynamic conditions for pipe flow as well as thermal boundary conditions in both the fluid and solid regions of the model. The hydrodynamic boundary conditions are provided in Table 2, and include a prescribed hydrodynamically-developed velocity at the inlet, pressure outlet and no slip conditions at the wall. The inlet velocity is reported in [2], and was experimentally determined for each trial by measuring the amount of time it took for the fluid level to rise from one level to another in a volume-measuring tank locate d within the loop. Those authors used a pair of contact points and an electric stopwatch to make this calculation. The turbulence inlet boundary conditions were inputted by using standard expressions for a fully developed internal flow, defining turbulence intensity and a reference length.

Solving the energy equation also required appropriate boundary conditions to be defined and are needed for both the fluid and the solid within the domain. These are shown in Tables 3 and 4, respectively. The fluid inlet temperature for each trial is taken from [2] and the outlet is considered adiabatic. At the surface between the molten salt and the solid tube wall, an interface condition is applied which serves to thermally couple the two element types (i.e., solid and fluid). For the solid, this interface condition is also applied, while the remaining solid surfaces are all adiabatic. The coupled boundary condition was set in the interface. The solution for the temperatures in the interface are solved for by using the information of the temperatures of the interface adjacent cells, iterating until the energy balance is satisfied at the interface.

Enhanced wall treatment was used during the simulations, which means that for the dimensionless velocity, a blending function is used to smoothly merge both the log layer and the viscous layer. For the thermal formulation, an elliptic blending is also used for merging the laminar and logarithmic profiles accordingly.

In Grele and Gedeon experiment, electrical heaters were fixed around the pipe circumference which provided the required heating. To model this, cell zone conditions were imposed in the solid region in order to add the volumetric heat source due to electrical heating. The constant volumetric heat source may have discrepancies with the actual experiment because of several factors such as non-uniformity of the electrical current and because of the variation of the pipe temperature. However, because of the lack of data regarding the heat flux profiles, it was decided that the constant volumetric heat source was the best choice for reproducing the experiment as it represents the simplest solution for applying the thermal load.

In each one of the cases, the total heat supplied to the fluid was computed using the temperatures at the inlet and the outlet together with the choice of material properties (Tab. 1), particularly density and specific heat. The formula used was the following: (8)

where superscripts x and y refer to the different possibilities of thermophysical properties (Tab. 1), T_b is the bulk temperature calculated as and A is the cross- sectional area of the pipe. After the total heat supplied to the fluid was computed, using the assumption of uniformly distributed heat, the volumetric heat source was inputted as a cell zone condition in the solid region.

4 Results

In order to assess the RANS models for molten salt flows (FLiNaK), different cases in the Grele and Gedeon experiment were simulated and compared to their experimental results. In the experiment, the external average wall temperature was reported in each run, which was found by simply averaging the 12 thermocouples located on the pipe wall at different axial positions. The wall temperature in the external wall of the modeled solid pipe at each one of the thermocouples’ positions was averaged in order to directly compare with the reported external wall temperature from the experiment.

The thermophysical properties of FLiNaK were considered as temperature dependent in the simulations, implemented as user defined functions in Fluent. Due to the lack of trustworthy material properties, for each case, different combinations of the material properties of FLiNaK presented in Table 1 were used. The thermal conductivity used was given by Smirnov’s correlation, as the majority of the published literature consider this correlation as the most accurate for FLiNaK. In regards to the solid region, the Inconel X pipe thermophysical properties were taken from [23].

From the 52 cases presented by Grele and Gedeon, 10 specific cases were simulated, and encompass the range of operating conditions encountered by the experimentalists. The key parameters taken from the experimental data were:

• for the fluid: the velocity inlet, the temperature inlet;

• for the solid: using equation (8), taking the inlet and outlet temperature, the uniform volumetric heat flux was calculated and imposed as a cell zone condition.

The 10 cases selected for the simulations are shown in Table 5.

For clarification, the range for the Reynolds and Prandtl numbers observed in column 6 and 7 of Table 5 corresponded to the minimum and maximum values using the aforementioned different combination of thermophysical properties.

Taking into account the two correlations for density, the two correlations for specific heat and the four correlations for viscosity, the total amount of 16 simulations were run per case, yielding a total of 160 simulations. Thus, in order to proceed in a practical way, a loop which went through all of the possible combinations was programmed in Scheme language, a language seamlessly interpreted by Fluent. In this script, the case was initialized, the properties were changed accordingly, the simulation was run andfinally results were exported. Using Matlab, a table containing all the pertinent results was created automatically. In particular, this code calculated the heat transfer coefficient and compared it with Dittus-Boelter and Sieder Tate correlation, both of which useReynolds and Prandtl numbers as inputs, and therefore a spread exists in that data as well due to dependence on thermophysical properties. Unsurprisingly, the Nusselt number spread using the Sieder Tate correlation was smaller in comparison to that using the Dittus Boelter correlation. On average, the simulation results showed a 1.4% difference when compared to the Sieder Tate correlation and a 4.9% of the Dittus Boelter correlation. This showed that the variation of the viscosity values at the wall surface with respect to the bulk viscosity was significant and could not be neglected.

The results for the averaged outer wall temperature are shown in Figure 5. On the y-axis, the difference between the experimental results and the simulation results for the averaged wall temperature was plotted. The x-axis denotes the specific case number as previously noted.

The most noticeable conclusion from Figure 5 is that the variations on the predicted thermophysical properties can lead to a wall temperature spread of almost 15 ℃. The constant line with a zero value represents the match between the simulation and the experimental results.

The initial thought with running these experiments was to identify the set of material properties that best describe the experimental results. The properties’ correlations that best describe the experiment were not consistent between any two experiments, meaning that the best set of material properties could not be found. It should be noted that even though for some cases (e.g., 30 and 46), results suggest that there is no overlap between experimental results and simulation results, the reader should bear in mind that experimental uncertainties were not reported and cannot accurately be accounted for at this stage.

Figure 6 shows the axial temperature profile of the CFD simulations compared to the measurement of the thermocouples in the external wall for case number 43. It is worth noting that thermocouples readings were shown for only 3 cases in [2]. As can be seen in the figure, the temperature distribution suggests the heat source is not uniform as it drops in the middle of the test section. However, it is not possible to recreate the experimental conditions given the lack of further details. Figure 6 shows the best combination of thermophysical properties along with the two combinations that represent the bounding cases. It can be seen that the uniform volumetric heat source approach is valid, especially with a lack of rationale for choosing a less standard approach to applying the thermal energy.

Due to the aforementioned spread in the wall temperature values, it was observed that the uncertainty in the material thermophysical properties could lead to a significant variance in the heat transfer properties of FLiNaK. This uncertainty coming from the material properties, added to the lack of reliable estimates of experimental uncertainties, impeded an appropriate validation of RANS models for molten salt flows. From this point onwards, it was assessed that a sensitivity analysis on how these properties’ individual variations affected heat transfer properties was needed.

Fig. 5 Difference between the wall temperature measured in the experiment and the wall temperature resulting from the simulations for each of the ten cases. The various data points in each case correspond to different selections of material properties of FLiNaK.

Fig. 6 Axial temperature on the outer wall of the test section measured with thermocuples in [2] and obtained in the simulations for case 43.

5 Sensitivity analysis

A widely used parameter to measure the sensitivity of a simulation result S(X) to changes in a parameter X_i is termed asthe sensitivity coefficient, whose definition is the following: (9)

where X_i is one element of the vector , which includes all dependencies of the variable S.

Using a Taylor series approximation, the uncertainty of a function S(X) with n uncorrelated parameters can be accounted for by using the following expression: (10)

For simple cases of an algebraic model, these sensitivity coefficients may be calculated analytically. However, the most common case is that the model is a complex numerical simulation where no algebraic model can be used. Using this analysis as an example, a change in one of the material properties may produce a change in the velocity profile of the fluid and also in the heat transfer properties between the liquid and the solid. Because of the existence of these mentioned coupled effects, the sensitivity coefficient was calculated with data from computational fluid dynamics simulations.

The procedure used in this section was to run the simulation with nominal values of the parameter vector . For the next simulations, perturbed values for the input parameter X_i are used. Finally, using a forward finite difference approximation in the parameter space, the sensitivity coefficient is calculated as follows [24]: (11)

For this particular case, the simulation result of interest S(X) was the heat transfer coefficient between the molten salt and the solid pipe. The vector of parameters consisted of the four thermophysical properties of FLiNaK.

It is worth pointing out that the One Factor At a Time Approach (OAT) [25] was used in the following analysis. In this approach, only one parameter changes its value between consecutive simulations, and so, in a deterministic model, the analyst can determine exactly what effect is caused by changing the parameter. In order to obtain good results, the model must be linear in the sense that Gaussian distributions are assumed between input and output values. In addition, there should not be a significant covariance between the input parameters of the sensitivity analysis. Even though it is known that the thermophysical properties are related to similar underlying molecular mechanisms [26] and that the model may have non-linear effects, the current analysis presents a first approach to sensitivity analysis in molten salt flows. Therefore, for the sake of simplification, the model response is assumed linear and that the covariance is not significant. Additional studies should be conducted in the future to validate this claim.

Taking into account every correlation presented in Table 1, an average temperature dependent function was calculated for each property. In order to do this, the temperature dependent expression for each property was evaluated at each temperature in the range of 800–1100K. Next, these values were averaged, yielding one average value of the property for each temperature. Finally, a fitting analysis was made to obtain the average function which better described the behaviour of the material property in question (see Fig. 7). From this point onwards, in each one of the simulations, only one of the properties was varied while keeping this same functional shape obtained in the fitting. The variations for the analysis were constant steps and the range of these variations was in agreement with respect to maximum and minimum values of the property of interest (as shown in Fig. 7). The remaining average temperature dependent properties were kept invariant for each analysis. As an example, the input viscosities for the sensitivity analysis are shown in Figure 8.

In Figure 7 the four models for viscosity in addition to the average viscosity in the temperature range of interest are plotted. From Figure 8, the average viscosity function and the simulation inputs for the sensitivity analysis canbe observed. An analogous process of obtaining an average function and adding constant steps to build the simulation inputs was done for the density, the thermal conductivity and the specific heat. Particularly for the thermal conductivity case Smirnov’s correlation [7] (shown in Tab. 1) was used, and the variations were done according to the studies done by Williams et al. [12], which states that the values should range between 0.6 W/m K and 1.0 W/m K at 973 K.

As it was mentioned before, the sensitivity analysis in this section was targeted at obtaining sensitivity coefficients for the heat transfer coefficient variable. Thus, each one of the eight black curves labeled as Simulation Inputs in Figure 8 were used inthe simulations in order to obtain the corresponding heat transfer coefficients as a function of the property value evaluated atthe bulk temperature. The simulations were done in three different temperature ranges, determined by different temperatureinlet boundary conditions. Hence, some of the typical temperatures of operation present in a MSR were covered.

The resulting heat transfer coefficients obtained in the simulations as a function of viscosity, density, specific heat and thermal conductivities are plotted in Figure 9. Taking as an example the viscosity sensitivty analysis, from Figure 9d it can be observed that there are eight different values of heat transfer coefficients for each temperature range. This matches the number of simulation inputs described in Figure 8. From the data shown in Figure 9, the results of the calculation of the sensitivity coefficients for different temperature ranges are shown in Table 6.

Note that the sensitivity coefficient is positive for three of the four properties considered with viscosity alone having a negative value. This is not surprising since an increase in viscosity would decrease the Reynolds number, thereby diminishing the amount of turbulence as well as the transfer of thermal energy.

Nusselt number calculations using Dittus Boelter and Sieder Tate correlations were performed in order to compare and assess the validity of using both of these correlations. The comparison was done using the results for the viscosity sensitivity analysis, which was arepresentative case (see Fig. 10).

As it can be seen from Figure 10, Sieder Tate predictions agree much better with CFD simulation results compared to Dittus Boelter. This is not surprising in that the Sieder Tate formula accounts for the variability in fluid properties near the wall, which can influence the heat transfer performance. The slope of the curves in Figure 10 also reveal a good agreement between Sieder Tate and our CFD simulations suggesting similar sensitivity from the viscosity on the heat transfer coefficient could be predicted based on the Sieder Tate correlation. It is important to note that the CFD simulations allow a direct calculation of the surface temperature, and therefore provide the evidence of the validity of the Sieder Tate formula in this scenario. Additional work should be performed with other turbulence models and higher fidelity simulations to more completely establish this validity.

Apart from the sensitivity coefficients presented above in equation (11), a dimensionless sensitivity parameter is an alternative way to provide information about how a certain parameter affects the variable of interest. Consequently, the dimensionless sensitivity coefficient is defined as: (12)

Equation (12) suggests a rather simple interpretation for S^*, namely a fractional variation in the calculated observable quantity (in this case the heat transfer coefficient) that is produced by a specific fractional change in a selected input parameter (in this case material properties). In equation (12), the observable quantity h is defined as the heat transfer coefficient values obtained by using the average values of the material properties x. For the viscosity sensitivity analysis, the heat transfer coefficient h appearing in equation (12) would be calculated with the average viscosity red curve in Figure 8 and the remaining property average values. The amount of variation of x is δx while the resulting change in h is δh. Using these definitions, the dimensionless sensitivity coefficients were calculated and are presented in Table 7. Also, the maximum variations in percentage of a material property according to the different correlations can be observed.

According to the results shown in Table 7, viscosity, thermal conductivity and density have the highest dimensionless sensitivity coefficient. Nevertheless, the percentage change in density is very small compared to the thermal conductivityand viscosity, making it less relevant in terms of how it might affect the prediction of the heat transfer coefficient. This is due to the fact that the determination of the value of density is a less complex task than the measurement of the other two properties. In contrast, the thermal conductivity values reported by [6] and [12], may vary between 0.6 and 1.0 W/m K approximately, which accounts for a 26% spread.

Therefore, to be able to producing high quality experimental data for future validation purposes, a more confident understanding of material properties is fundamental. More specifically to future molten salt efforts, thermal conductivity and viscosity need to be accurately described since they are the most sensitive parameter.

Fig. 7 Temperature dependent correlations for viscosity and the average viscosity function obtained by exponential fitting.

Fig. 8 Viscosity inputs for the sensitivity analysis CFD simulations.

Fig. 9 Results of the sensitivity analysis 34.

Fig. 10 Heattransfer coefficient versus thermal conductivity values obtained from the sensitivity analysis. Comparisons with Dittus Boelter and Sieder Tate weremade.

6 Conclusions

Reynolds averaged Navier Stokes using the k − ω SST model was used to perform computational fluid dynamics simulations attempting to reproduce the experiment results obtained by Grele and Gedeon [2] in 1954. Lack of detail and uncertainty estimates of the experimental data prevented a full validation effort. Despite the inability to conducting full validation studies, the influence of the variation of the thermophysical properties was gauged as they pertained to external wall temperatures and heat transfer coefficients. A large spread of values was observed for the external wall temperature of the simulations.

In light of the uncertainty of FLiNaK’s thermophyisical properties, a sensitivity analysis using a one-factor at a time approach (OAT) was done to account for the influence of these variations on the heat transfer coefficient for each one of the properties. The analysis proved that viscosity and thermal conductivity are the most crucial properties for conjugate heat transfer simulations using FLiNaK. The simulations also demonstrated that Sieder Tate correlation can be used with a reasonable level of trust since the results agreed well with that expression. It is not surprising that the Dittus Boelter expression did not agree as well, especially considering it cannot account for different viscosity at the wall compared to the bulk. This considered analysis can be extended to other molten salts due to the similarity in the dimensionless numbers (Prandtl and Reynolds) used in the study.

This work provides evidence on the importance of conducting further research in the material properties for FLiNaK, as they have a significant effect in the prediction of the heat transfer coefficient. Additional turbulence models should also be run to further validate the use of the Sieder Tate correlation and gain further insight into the flow physics.

For future studies, a sensitivity analysis using Large Eddy Simulations (LES) and DNS is suggested with the aim of reducing the impact arising from the error in the turbulence model used and the assumptions built into those models. In addition, it is important to note that this entire analysis was based on main effects of the material properties. Further studies to quantify any interactions between the variations of more than one property at the same time are suggested. This potential analysis would serve to quantify the covariance between the uncertainty of two properties. In this way a probability density functions of inputs and outputs in the current sensitivity study could be characterized in order to obtain the non-linearity effects of the thermophysical properties’ effect on the heat transfer coefficient.

Author contribution statement

Ramiro Freile has written the article. Mark Kimber has contributed to this work by providing financial support through a federal grant and advising on the different topics discussed.

References

All Tables

All Figures

	Fig. 1 Experimental results obtained by Grele and Gedeon [2] and Hoffman and Loans [3] compared with Dittus Boelter correlation.
In the text

	Fig. 2 Experimental results obtained by Vriesema [4] compared with Dittus Boelter.
In the text

	Fig. 3 Grele and Gedeon, Hoffman and Loans and Vriesema results for Inconel piping reanalyzed by Ambrosek et al. [8].
In the text

	Fig. 4 Perpendicular section of the computational domain used in the simulations.
In the text

	Fig. 5 Difference between the wall temperature measured in the experiment and the wall temperature resulting from the simulations for each of the ten cases. The various data points in each case correspond to different selections of material properties of FLiNaK.
In the text

	Fig. 6 Axial temperature on the outer wall of the test section measured with thermocuples in [2] and obtained in the simulations for case 43.
In the text

	Fig. 7 Temperature dependent correlations for viscosity and the average viscosity function obtained by exponential fitting.
In the text

	Fig. 8 Viscosity inputs for the sensitivity analysis CFD simulations.
In the text

	Fig. 9 Results of the sensitivity analysis 34.
In the text

	Fig. 10 Heattransfer coefficient versus thermal conductivity values obtained from the sensitivity analysis. Comparisons with Dittus Boelter and Sieder Tate weremade.
In the text