Development and Verification of CUPID-MSR Using the de Vahl Davis Natural-Convection Benchmark

1. Introduction

In recent years, molten-salt reactors (MSRs) have attracted growing attention because they employ liquid fuel and offer enhanced safety behavior. In these systems, the fissile material is mixed directly into molten salt and circulated through the primary circuit. This contrasts with conventional light-water reactors (LWRs), in which the fuel remains solid. The molten salt operates at temperatures approaching $700°\mathrm{C}$, enabling high thermal efficiency and potentially improved economic performance. A key advantage of MSRs is their strong passive safety potential. In accident situations, the fuel salt can be drained into dedicated tanks, either automatically or through operator action, where the resulting configuration becomes subcritical and the progression of the event is halted.

Historically, the first nuclear reactor to use molten salts was the Aircraft Reactor Experiment (ARE), built in 1954 at Oak Ridge National Laboratory (ORNL), USA, and designed to power a nuclear-propelled bomber [1]. Although the aircraft nuclear propulsion program was eventually abandoned, ORNL later achieved first criticality with the Molten-Salt Reactor Experiment (MSRE) in 1965 [2]. This $8{\mathrm{MW}}_{\mathrm{th}}$ test reactor, fueled with LiF–${\mathrm{BeF}}_2$–${\mathrm{ZrF}}_4$–${\mathrm{UF}}_4$, operated successfully for 4.5 years before shutting down in 1969. The MSRE was widely regarded as a success, demonstrating the viability of molten salt fuels despite corrosion challenges. However, MSR development in the United States was halted as attention shifted toward other reactor concepts.

Currently, no MSR is in operation worldwide, but MSRs remain a promising alternative to water-cooled reactors, as demonstrated by numerous international projects, including Terrestrial Energy’s IMSR-400 in Canada [3], Kairos Power’s KP-FHR in the United States [4], and Saltfoss Energy’s Compact Molten Salt Reactor in Denmark [5]. South Korea is also advancing MSR development through the Korea Molten Salt Reactor (K-MSR), a $100{\mathrm{MW}}_{\mathrm{th}}$ reactor intended for marine propulsion [6,7]. In 2024, the Korea Atomic Energy Research Institute (KAERI) also signed a Memorandum of Understanding with Saltfoss Energy to collaborate on MSR technology development [8].

Reactor concepts that rely on coolants other than water introduce unconventional design features and require adapted safety strategies. Because their thermal-hydraulic behavior differs significantly from that of water-cooled reactors, extensive safety assessments are necessary, and these must rely on verified numerical tools. Most existing power plants employ light water as the primary coolant; therefore, commonly used system analysis codes such as RELAP5 [9], CATHARE [10], and MARS [11] were developed specifically for LWR applications. As a result, these codes are not well suited for analyzing MSRs, leaving a gap in available computational tools. To address this need, KAERI developed the system code GAMMA for reactor concepts using non-water coolants, including various MSR designs [12,13,14]. Such tools are commonly referred to as system codes because they represent the overall reactor layout and major components using simplified one-dimensional models.

In LWR cores, coolant channels are arranged axially and the flow proceeds primarily in the vertical direction. Under these conditions, one-dimensional models can provide acceptable accuracy. In MSRs, however, particularly in fast-spectrum designs, the core typically lacks internal flow-guiding structures such as fuel assemblies, control rods, and graphite moderators. In the absence of solid structures that constrain the flow path, radial circulation develops in addition to axial motion, leading to considerable turbulent mixing within the reactor vessel. These three-dimensional flow effects make system-scale one-dimensional modeling insufficient for accurately capturing MSR core thermal-hydraulics.

For this reason, CUPID-MSR was developed by extending the CUPID code with correlations and models appropriate for molten salt coolants [15,16]. CUPID has been developed at KAERI since 2011 and is a component-scale thermal-hydraulics solver for transient two-phase flow based on a two-fluid formulation. Depending on the selected computational resolution, it can be applied in component-scale, porous-media, or full CFD modes. In recent years, significant development efforts have focused on extending CUPID’s applicability to coolant fluids beyond water, including molten salts.

At present, CUPID-MSR includes two molten salt mixtures: potassium chloride–uranium(III) chloride (KCl–${\mathrm{UCl}}_3$) and sodium chloride–magnesium chloride–transuranic(III) chloride (NaCl–${\mathrm{MgCl}}_2$–${\mathrm{TRUCl}}_3$). In molten salt reactors, the circulating fuel serves simultaneously as the fissile medium and the coolant, resulting in a strong coupling between neutron behavior and thermal-hydraulic processes. This interaction involves reactivity feedback mechanisms as well as the transport of delayed neutron precursors and fission products within the flowing salt. Consequently, MSR simulations must employ coupled neutron-kinetics and thermal-hydraulics models within a multiphysics framework. Although CUPID provides multiphysics coupling with a three-dimensional neutron kinetics code [17], its current application is limited to LWRs, as suitable three-dimensional neutron kinetics tools for MSRs are not yet available. With future coupling to an MSR neutronics code, CUPID-MSR is expected to support the design and safety analysis of advanced MSR systems.

To assess the reliability of CUPID-MSR for molten salt applications, its predictions are compared with results from a well-established benchmark problem. For this purpose, the natural convection cavity problem proposed by de Vahl Davis [18] is employed, and results are evaluated for both molten salt mixtures.

In addition to code verification, this study examines the applicability of the Boussinesq approximation. Although widely used to simplify buoyancy-driven flow simulations, its validity over large temperature differences has not been systematically assessed for fluids with strongly temperature-dependent properties. By comparing the full variable-density formulation of CUPID-MSR with the Boussinesq model over a wide range of density variations, a practical applicability criterion is identified in terms of the nondimensional parameter $\beta \Delta T=\Delta \rho /{\rho }_{\mathrm{ref}}$. Because this parameter represents the relative density change, the results are relevant not only for molten salts but also for natural convection flows in a broad class of Newtonian fluids.

Section 2 outlines the implementation of molten salt thermophysical properties in CUPID-MSR. Section 3 presents the code verification using the de Vahl Davis benchmark. Section 4 introduces the assessment of the Boussinesq approximation and discusses its general applicability. The conclusions are summarized in Section 5.

2. Implementation of Molten-Salt Thermophysical Properties

The CUPID-MSR code was developed to extend the CUPID thermal-hydraulic solver for applications involving molten salt reactor systems, in which the working fluid is a molten salt rather than light water. Among the various candidate fuel salts, two chloride-based mixtures were implemented: KCl–${\mathrm{UCl}}_3$ and NaCl–${\mathrm{MgCl}}_2$–${\mathrm{TRUCl}}_3$. The eutectic temperature of the NaCl–${\mathrm{MgCl}}_2$–${\mathrm{TRUCl}}_3$ mixture is lower than that of KCl–${\mathrm{UCl}}_3$, making it particularly attractive from an operational standpoint.

Each salt mixture is represented using temperature-dependent correlations for the main thermophysical properties, namely density, specific heat capacity, thermal conductivity, and dynamic viscosity. These correlations are implemented in CUPID-MSR through a dedicated material property table file, allowing the code to capture the strong coupling between temperature and fluid behavior that is characteristic of molten salts.

2.1. Density Correlation

The density of each salt mixture was calculated using an empirical linear model for the density of its pure components i, expressed as

$${\rho }_i=A_i+B_{i}T,$$

(1)

where $A_i$ and $B_i$ are empirical constants and T is the temperature. The molar volume contribution of each component i in the mixture is given by

$$x_i{V}_i={\displaystyle \frac{x_i{M}_i}{{\rho }_i}},$$

(2)

where $x_i$ and $M_i$ denote the mole fraction and molar mass of component i, respectively, and the product corresponds to the mass contribution of that component. The final expression for the density of the entire salt mixture is given by

$${\rho }_{\mathrm{mix}}={\displaystyle \frac{{\sum }_i{x}_i{M}_i}{{\sum }_i{x}_i{V}_i}}.$$

(3)

From the linear correlations, the following relations were derived for the two implemented salts:

$${\rho }_{{\mathrm{KCl}-\mathrm{UCl}}_3}=3.69\times {10}^3-0.39T,$$

(4)

and

$${\rho }_{{\mathrm{NaCl}-\mathrm{MgCl}}_2{-\mathrm{TRUCl}}_3}=3.93\times {10}^3-0.33T.$$

(5)

2.2. Heat Capacity Correlation

The specific heat capacity determines the amount of energy required to raise the fluid temperature by one kelvin and was calculated using a polynomial correlation for each component:

$$C_{p,i}=A_i+B_{i}T+C_i{T}^2+D_i{T}^3+E_i{T}^{-2},$$

(6)

where $A_i$, $B_i$, $C_i$, $D_i$, and $E_i$ are empirical constants. The final expression for the specific heat capacity of the mixture is given by

$$C_{p,\mathrm{mix}}={\displaystyle \frac{{\sum }_i{x}_i{C}_{p,i}}{{\sum }_i{x}_i{M}_i}},$$

(7)

where $x_i{C}_{p,i}$ represents the contribution of component i, with $x_i$ denoting its mole fraction and $C_{p,i}$ its molar heat capacity. In CUPID-MSR, the KCl–${\mathrm{UCl}}_3$ mixture exhibits an approximately constant $C_p$ over the investigated temperature range, whereas NaCl–${\mathrm{MgCl}}_2$–${\mathrm{TRUCl}}_3$ shows stronger temperature dependence between 500 K and 1000 K before stabilizing at higher temperatures.

2.3. Thermal Conductivity Correlation

Thermal conductivity $\lambda$ characterizes the ability of the salt to conduct heat and was modeled using temperature-dependent linear relations for each component i, expressed as

$${\lambda }_i=A_i-B_i(T-T_{\mathrm{m},i}),$$

(8)

where T is the local temperature and $T_{\mathrm{m},i}$ is the melting temperature of component i. The effective thermal conductivity of the salt mixture is computed as the weighted sum of the individual component contributions:

$${\lambda }_{\mathrm{mix}}=\sum _i{x}_i{\lambda }_i.$$

(9)

As shown in Figure 1(c), both salt mixtures exhibit decreasing thermal conductivity with increasing temperature, while NaCl–${\mathrm{MgCl}}_2$–${\mathrm{TRUCl}}_3$ consistently maintains higher values. This behavior is advantageous for achieving a more uniform temperature field during reactor operation.

2.4. Viscosity Correlation

In CUPID-MSR, the dynamic viscosity of the mixture is calculated from the viscosities of the individual components using an Arrhenius-type expression:

$${\mu }_i=A_{i}exp\left({\displaystyle \frac{B_i}{RT}}\right),$$

(10)

where R is the universal gas constant, equal to $8.314\mathrm{J}{\mathrm{mol}}^{-1}{\mathrm{K}}^{-1}$. The dynamic viscosity of the mixture is evaluated using a logarithmic mixing rule:

$${\mu }_{\mathrm{mix}}=exp\left(\sum _i{x}_{i}ln{\mu }_i\right).$$

(11)

Both salt mixtures exhibit an exponential decrease in viscosity with increasing temperature, as shown in Figure 1(d). This negative feedback mechanism is beneficial for reactor safety: as local temperature increases, viscosity decreases, buoyancy-driven flow intensifies, and convective heat removal is enhanced, thereby contributing to passive system stabilization.

2.5. Summary of Property Behavior

Figure 1 summarizes the implemented thermophysical property correlations. The observed trends confirm that both chloride-based salts exhibit high density, moderate thermal conductivity, and rapidly decreasing viscosity with increasing temperature. These characteristics are favorable for establishing efficient and stable natural-circulation flow in molten salt reactor systems. The temperature-dependent correlations were implemented in CUPID-MSR to enable realistic simulation of molten salt reactor operating conditions.

3. Verification Against the de Vahl Davis Benchmark

3.1. Problem Description

The molten salt property implementation in CUPID-MSR was verified using a two-dimensional thermal cavity benchmark problem originally proposed by de Vahl Davis [18]. Figure 2 presents the schematic of the thermal cavity configuration and the computational mesh employed. The square cavity ($H=W$) is oriented such that gravity acts perpendicular to its surface. Constant wall temperature boundary conditions are imposed on the left and right walls, with the left wall maintained at a higher temperature than the right wall. The horizontal boundaries at the top and bottom are treated as adiabatic. Under these conditions, buoyancy-driven natural convection is expected to develop within the cavity.

In this thermal cavity problem, the formation of natural circulation is governed by the Rayleigh number ($Ra$), defined as

$$Ra={\displaystyle \frac{{\rho }^{2}g\beta (T_{\mathrm{hot}}-T_{\mathrm{cold}})H^3}{{\mu }^2}}Pr,$$

(12)

where $\rho$, $\beta$, and g denote the fluid density, thermal expansion coefficient, and gravitational acceleration, respectively. The Prandtl number ($Pr$) is defined as $Pr=\left(C_p\mu \right)/\lambda$, where $C_p$ is the specific heat capacity and $\lambda$ is the thermal conductivity. As the Rayleigh number increases, the flow regime gradually transitions from laminar toward turbulent.

The calculations were performed for ${10}^3\le Ra\le {10}^6$, corresponding to laminar flow conditions in the cavity. The hot and cold wall temperatures were set to 900.5 K and 899.5 K, respectively, which are representative of the typical operating temperature range of proposed molten salt reactor designs.

3.2. Nusselt Number Evaluation Method

In this benchmark, the Nusselt number ($Nu$) is used as the primary metric for evaluating the verification results. The Nusselt number is a dimensionless quantity that characterizes the intensity of convective heat transfer relative to pure conduction. A value of $Nu=1$ corresponds to purely conductive heat transfer, while values greater than unity indicate enhancement due to convection. As the Rayleigh number increases, buoyancy forces intensify, leading to higher fluid velocities and an increased heat transfer coefficient, and consequently a larger Nusselt number.

In this work, the local heat transfer near the hot wall is quantified using the local Nusselt number, defined as

$$Nu=-{\displaystyle \frac{H}{\Delta T}}{\left({\displaystyle \frac{\partial T}{\partial x}}\right)}_{\mathrm{hot}},$$

(13)

where $\Delta T=T_{\mathrm{hot}}-T_{\mathrm{cold}}$. When a nondimensional temperature formulation is employed, the Nusselt number can be expressed as

$$Nu={\left({\displaystyle \frac{\partial T^{*}}{\partial x^{*}}}\right)}_{\mathrm{hot}},$$

(14)

with the non-dimensional variables defined as

$$T^{*}={\displaystyle \frac{T-T_{\mathrm{cold}}}{T_{\mathrm{hot}}-T_{\mathrm{cold}}}},x^{*}={\displaystyle \frac{x}{H}}.$$

(15)

The average Nusselt number along the hot wall is obtained by integrating the local Nusselt number over the vertical direction:

$$\overline{Nu}={\displaystyle \frac{1}{H}}{\int }_0^{H}Nudy.$$

(16)

For numerical evaluation, this integral is approximated using Simpson’s rule:

$$\overline{Nu}={\displaystyle \frac{1}{6}}\left[Nu\left(0\right)+4Nu\left({\displaystyle \frac{H}{2}}\right)+Nu\left(H\right)\right].$$

(17)

Accordingly, the local Nusselt number is evaluated at three vertical locations along the hot wall: the bottom, mid-height, and top of the cavity.

3.3. Mesh Sensitivity Assessment

A mesh sensitivity study was conducted to ensure that the predicted Nusselt numbers are independent of grid resolution. The two-dimensional cavity was discretized using structured meshes ranging from $20\times 20$ to $120\times 120$ cells for both KCl–${\mathrm{UCl}}_3$ and NaCl–${\mathrm{MgCl}}_2$–${\mathrm{TRUCl}}_3$ salt mixtures. All simulations were performed at a reference temperature of 900 K. For each Rayleigh number, calculations were carried out on all six mesh resolutions to evaluate the influence of grid refinement on the numerical solution.

As an initial assessment, Figure 3 presents the column-averaged temperature profile $\overline{T}(x/H)$ across the cavity, obtained by averaging the temperature along the vertical direction at each non-dimensional horizontal position $x/H$. This global metric reflects the overall thermal distribution within the cavity and provides a convenient measure of convergence without focusing on localized quantities. For clarity, results obtained using four meshes ($20\times 20$, $60\times 60$, $80\times 80$, and $120\times 120$) are shown.

Across all Rayleigh numbers and for both salt mixtures, temperature profiles obtained on the finer meshes ($60\times 60$, $80\times 80$, and $120\times 120$) collapse onto a single curve, indicating that the global temperature distribution becomes essentially mesh-independent once a moderately fine grid is used. The coarsest mesh ($20\times 20$) exhibits only minor deviations, which become noticeable primarily at the highest Rayleigh numbers. This global agreement provides qualitative evidence of mesh independence prior to examining more sensitive near-wall quantities.

A more sensitive indicator of mesh dependence is the temperature in the first column of cells adjacent to the hot wall, as shown in Figure 4. Because this region contains the steepest thermal gradients, it is particularly sensitive to grid resolution and is therefore well suited for identifying the mesh density required to accurately resolve the thermal boundary layer, especially at higher Rayleigh numbers. As the mesh is refined, the near-wall temperature monotonically approaches the imposed hot-wall temperature of 900.5 K, demonstrating consistent convergence behavior.

Both molten salt mixtures exhibit nearly identical convergence trends, with differences between them remaining negligible across all Rayleigh numbers. At higher Rayleigh numbers, coarse meshes may lead to distorted solutions, requiring finer grids compared with low-$Ra$ cases. The mesh selection criterion adopted in this study is the coarsest grid for which further refinement produces only marginal changes in the near-wall average temperature, specifically when the relative deviation with respect to the $120\times 120$ reference solution falls below ${10}^{-3}$ (less than 0.1%). Based on this criterion, the selected meshes are $40\times 40$ for $Ra={10}^3$, $60\times 60$ for $Ra={10}^4$ and ${10}^5$, and $80\times 80$ for $Ra={10}^6$. A summary of the selected mesh resolutions and their relative deviations is provided in Table 1.

3.4. Comparison with Benchmark Solution

Figures 4 and 5 illustrate the temperature contour fields for KCl–${\mathrm{UCl}}_3$ and NaCl–${\mathrm{MgCl}}_2$–${\mathrm{TRUCl}}_3$ at various Rayleigh numbers. The corresponding streamline patterns are shown in Figure 6 and 7. In all cases, buoyancy-driven natural convection develops clockwise within the cavity. As the Rayleigh number increases, the buoyant flow intensifies, leading to thinner thermal boundary layers and the appearance of additional corner vortices. At $Ra={10}^3$ and $Ra={10}^4$, the flow is dominated by a single recirculating cell; at $Ra={10}^5$, two vortices are observed; and at $Ra={10}^6$, three distinct vortical structures emerge. The progressive reduction of boundary-layer thickness with increasing $Ra$ reflects enhanced heat transfer at both the hot and cold walls.

A qualitative comparison of the temperature contours (Figure 5 and 6) and streamline patterns (Figure 7 and 8) with the reference de Vahl Davis benchmark confirms that CUPID-MSR successfully reproduces the characteristic temperature distributions and flow-field structures. The close agreement in vortex formation and isotherm topology provides strong qualitative validation of the implemented buoyancy and thermophysical property models and confirms the correct computation of the Rayleigh number in the code.

Figure 5. KCl–${\mathrm{UCl}}_3$ temperature contours for different Ra.

Figure 5. KCl–${\mathrm{UCl}}_3$ temperature contours for different Ra.

Figure 6. NaCl–${\mathrm{MgCl}}_2$–${\mathrm{TRUCl}}_3$ temperature contours for different Ra.

Figure 6. NaCl–${\mathrm{MgCl}}_2$–${\mathrm{TRUCl}}_3$ temperature contours for different Ra.

Figure 7. KCl–${\mathrm{UCl}}_3$ flow streamlines for different Ra.

Figure 7. KCl–${\mathrm{UCl}}_3$ flow streamlines for different Ra.

Figure 8. NaCl–${\mathrm{MgCl}}_2$–${\mathrm{TRUCl}}_3$ flow streamlines for different Ra.

Figure 8. NaCl–${\mathrm{MgCl}}_2$–${\mathrm{TRUCl}}_3$ flow streamlines for different Ra.

Quantitatively, the average Nusselt numbers predicted by CUPID-MSR are compared with benchmark data in Table 2 for both molten salt mixtures. The relative deviation of the average Nusselt number is defined as the percentage difference between the computed value and the benchmark reference at the same Rayleigh and Prandtl numbers:

$${\epsilon }_{Nu}={\displaystyle \frac{\left|{\overline{Nu}}_{\mathrm{CUPID}-\mathrm{MSR}}-{\overline{Nu}}_{\mathrm{DVD}}\right|}{{\overline{Nu}}_{\mathrm{DVD}}}}\times 100\%.$$

(18)

The deviations from the reference solution range from 0.66% to 3.86% for KCl–${\mathrm{UCl}}_3$ and from 0.38% to 1.78% for NaCl–${\mathrm{MgCl}}_2$–${\mathrm{TRUCl}}_3$. A direct comparison between the two salt mixtures shows that both follow the same heat transfer trend, with only minor differences in the predicted Nusselt numbers. The maximum difference between the two solutions remains within 1.04% for $Ra\le {10}^3$, 2.24% for $Ra\le {10}^4$, 1.27% for $Ra\le {10}^5$, and 2.04% for $Ra\le {10}^6$.

In both molten salt cases, the predicted Nusselt numbers increase monotonically with Rayleigh number, consistent with the benchmark trend and reflecting progressively stronger convective heat transport. The discrepancy between the numerical predictions and the reference solution increases slightly at higher Rayleigh numbers, as the flow approaches transitional or weakly turbulent regimes.

3.5. Temperature-Variation Sensitivity

To further assess the temperature dependence of the implemented thermophysical property correlations, an additional sensitivity analysis was performed by varying the reference temperature of the thermal cavity. The objective of this study was to quantify how changes in the molten salt fluid temperature influence the accuracy of the predicted average Nusselt number in CUPID-MSR relative to the benchmark solution.

Simulations were carried out for both molten salt mixtures over a range of Rayleigh numbers, with fluid temperatures spanning from 800 K to 2000 K and appropriately selected temperature increments, while maintaining the same boundary-condition configuration as in the reference case. Because molten salt properties such as density, viscosity, thermal conductivity, and heat capacity exhibit strong temperature dependence, even modest temperature variations can significantly alter the resulting Rayleigh and Prandtl numbers and, consequently, the heat transfer behavior.

Figures 9 and 10 present the relative deviations of the average Nusselt numbers from the de Vahl Davis benchmark for the KCl–${\mathrm{UCl}}_3$ and NaCl–${\mathrm{MgCl}}_2$–${\mathrm{TRUCl}}_3$ mixtures, respectively. For KCl–${\mathrm{UCl}}_3$, the deviation remains below 4.3%, while for NaCl–${\mathrm{MgCl}}_2$–${\mathrm{TRUCl}}_3$ it consistently remains below 3.0%. These results confirm that the implemented property correlations provide stable and reliable predictions over a wide temperature range.

Figure 9. Relative error compared to the de Vahl Davis solution for the KCl–${\mathrm{UCl}}_3$ case.

Figure 9. Relative error compared to the de Vahl Davis solution for the KCl–${\mathrm{UCl}}_3$ case.

Figure 10. Relative error compared to the de Vahl Davis solution for the NaCl–${MgCl}_2$–${TRUCl}_3$ case.

Figure 10. Relative error compared to the de Vahl Davis solution for the NaCl–${MgCl}_2$–${TRUCl}_3$ case.

The results also indicate that the deviation in the Nusselt number generally increases with increasing Rayleigh number. This behavior can be attributed to stronger buoyancy-driven convection at higher $Ra$ values. At low Rayleigh numbers, heat transfer is dominated by conduction and the velocity field is weak, resulting in smooth temperature gradients and close agreement with the benchmark solution. As $Ra$ increases, buoyancy forces intensify and more pronounced circulation develops, producing steeper velocity and temperature gradients near the hot and cold walls. These gradients are more challenging to resolve numerically, leading to slightly larger deviations from the reference data.

Overall, this sensitivity analysis demonstrates that both molten salt mixtures maintain good agreement with the benchmark solution across the investigated temperature range, verifying that CUPID-MSR can reliably simulate molten salt reactor conditions under typical operating temperatures as well as during transients involving moderate temperature increases.

4. Assessment of the Boussinesq Approximation

Natural convection in molten salt systems is strongly influenced by temperature-dependent density variations. In many CFD and system codes, buoyancy effects are simplified using the Boussinesq approximation, which assumes constant density in all terms of the momentum equation except for the buoyancy force, where a linear dependence on temperature is retained. Although this approximation is widely applied to water and air flows, its applicability to molten salts—characterized by high density, strong thermophysical gradients, and potentially nonlinear $\rho \left(T\right)$ behavior—has not been systematically assessed.

Because CUPID-MSR relies on temperature-dependent molten salt properties to accurately model buoyancy-driven flow, it is important to determine the temperature range over which the Boussinesq approximation remains valid for chloride-based MSR fluids. Therefore, this section evaluates the accuracy of the Boussinesq approximation by comparing it with the full variable-density formulation of CUPID-MSR using the de Vahl Davis natural-convection benchmark. Although molten salt thermophysical properties are employed, the applicability of the Boussinesq approximation is governed by the nondimensional parameter $\beta \Delta T$, which represents the relative density variation and is therefore independent of the specific working fluid. Consequently, the results obtained here are relevant not only for molten salts but also for a broad class of buoyancy-driven flows.

4.1. Full Variable-Density and Boussinesq Models

The Boussinesq approximation assumes that the fluid density remains constant throughout the computational domain, except in the buoyancy term of the momentum conservation equation, where a linear dependence on temperature is retained. This simplification reduces the computational cost of solving the Navier–Stokes equations while still capturing the dominant driving mechanism of natural convection.

In the full momentum conservation equation, density is treated as a temperature-dependent variable and appears in all inertial and gravitational terms. The governing equation can be expressed as

$$\rho \left(T\right)\left({\displaystyle \frac{\partial \overrightarrow{u}}{\partial t}}+\overrightarrow{u}·\nabla \overrightarrow{u}\right)=-\nabla p+\mu {\nabla }^2\overrightarrow{u}+\rho \left(T\right)\overrightarrow{g},$$

(19)

where $\rho \left(T\right)$ denotes the temperature-dependent fluid density.

When the Boussinesq approximation is applied, the momentum equation is written as

$${\rho }_{\mathrm{ref}}\left({\displaystyle \frac{\partial \overrightarrow{u}}{\partial t}}+\overrightarrow{u}·\nabla \overrightarrow{u}\right)=-\nabla p+\mu {\nabla }^2\overrightarrow{u}+{\rho }_{\mathrm{ref}}\beta (T-T_{\mathrm{ref}})\overrightarrow{g},$$

(20)

where T is the local fluid temperature, $T_{\mathrm{ref}}$ is the reference temperature, and ${\rho }_{\mathrm{ref}}$ is the density evaluated at $T_{\mathrm{ref}}$. The coefficient $\beta$ denotes the thermal expansion coefficient evaluated at the reference temperature.

These two formulations highlight the fundamental difference between the full variable-density model and the Boussinesq approximation. In the full formulation, $\rho \left(T\right)$ affects both inertial and buoyancy terms, whereas under the Boussinesq approximation density is treated as constant in all terms except for the buoyancy force, where its linear temperature dependence is preserved. This approach retains the correct physical mechanism driving natural convection while significantly simplifying the governing equations.

4.2. Numerical Methodology

The applicability of the Boussinesq approximation was assessed using the same de Vahl Davis cavity configuration described in Section 3. All simulations in this study were performed for the KCl–${\mathrm{UCl}}_3$ salt mixture at a reference temperature of 900 K and at a fixed Rayleigh number of ${10}^3$. All geometric parameters, boundary conditions, and solver settings were kept identical to ensure that differences in the predicted Nusselt numbers originate solely from the treatment of density.

The temperature difference between the hot and cold walls was varied from 0.2 K to 400 K. For each case, the thermal expansion coefficient $\beta$ was obtained from the implemented density correlations, including the O1 linear reference model and two O3 nonlinear variants (O3–Variant A and O3–Variant B). For each density model, two types of simulations were conducted:

• the full variable-density formulation, and
• the Boussinesq formulation, in which density is treated as constant in all momentum terms except for buoyancy.

The latter was implemented through a minor modification of the CUPID-MSR source code to ensure that density remained fixed in all non-buoyancy momentum contributions. Deviations from the benchmark were evaluated using the same relative error definition introduced in Section 3.4. Presenting the results as a function of the nondimensional parameter $\beta \Delta T$ enables a direct assessment of the temperature range over which the Boussinesq approximation remains accurate for molten salt natural convection. Because $\beta \Delta T$ represents the relative density variation, the same methodology is applicable to any Newtonian fluid.

4.3. Linearization of Density

Figure 11 compares the three density–temperature correlations examined in this study. The O1 linear correlation corresponds to the original KCl–${\mathrm{UCl}}_3$ density model introduced in Eq. (4). This linear model is equivalent to retaining only the first-order term of the Taylor expansion of the density field about the reference temperature $T_{\mathrm{ref}}=900\mathrm{K}$:

$$\rho \left(T\right)=\rho \left(T_{\mathrm{ref}}\right)+{\left.{\displaystyle \frac{\partial \rho }{\partial T}}\right|}_{T_{\mathrm{ref}}}(T-T_{\mathrm{ref}})+{\displaystyle \frac{1}{2}}{\left.{\displaystyle \frac{{\partial }^2\rho }{\partial T^2}}\right|}_{T_{\mathrm{ref}}}{(T-T_{\mathrm{ref}})}^2+\dots ,$$

(21)

Truncating this expansion after the linear term yields the familiar form

$$\rho \left(T\right)\approx {\rho }_{\mathrm{ref}}\left[1-\beta (T-T_{\mathrm{ref}})\right],$$

(22)

which is exactly the density model employed under the Boussinesq approximation. Because the Boussinesq formulation preserves only the linear term of the density expansion, it is expected that the Boussinesq and full variable-density formulations yield identical results when the O1 linear density correlation is used.

To investigate the influence of higher-order density effects, two nonlinear cubic O3 variants were constructed by modifying the curvature of the density profile while keeping $\rho \left(T_{\mathrm{ref}}\right)$ and $\beta \left(T_{\mathrm{ref}}\right)$ unchanged. O3–Variant A introduces negative curvature, causing density to decrease more rapidly with temperature on both sides of $T_{\mathrm{ref}}$. In contrast, O3–Variant B introduces positive curvature, resulting in a slower density decrease near 900 K and higher density values at elevated temperatures. As shown in Figure 11, both nonlinear variants deviate from the O1 linear correlation by approximately $\pm 5.6\%$ at 700 K and $\pm 7.6\%$ at 1100 K.

These controlled modifications isolate the influence of density curvature on buoyancy-driven flow behavior without altering the reference density at 900 K.

4.4. Results and Discussion

Figures 12 and 13 present the deviations of the predicted average Nusselt number from the de Vahl Davis benchmark as a function of the nondimensional parameter $\beta \Delta T$ for the O1 linear density correlation and the two O3 nonlinear variants.

For the linear density correlation (Figure 12), the Boussinesq formulation and the full variable-density model yield nearly identical deviations over the entire range of $\beta \Delta T$. This behavior is expected, as discussed in Section 4.1: the Boussinesq approximation retains only the first-order term of the Taylor expansion of $\rho \left(T\right)$, which corresponds exactly to the O1 linear density model. As a result, both formulations produce identical buoyancy forces and therefore the same natural-convection behavior.

For the nonlinear O3 density correlations (Figure 13), the variable-density results progressively deviate from the linear trend as $\beta \Delta T$ increases, reflecting the growing influence of higher-order curvature in the density–temperature relationship. Variant A, which introduces negative curvature, causes density to decrease more rapidly with temperature than in the linear model. This leads to stronger buoyancy forces, enhanced convection, and higher Nusselt numbers. Consequently, the Boussinesq approximation underestimates convective heat transfer for this variant. In contrast, Variant B exhibits positive curvature, resulting in a slower decrease in density and therefore weaker buoyancy forces. In this case, the Boussinesq approximation overestimates the convective heat transfer.

Despite these opposite tendencies, the Boussinesq formulation yields identical predictions for both nonlinear variants because it retains only the linear term $\beta (T-T_{\mathrm{ref}})$ of the density expansion and is therefore insensitive to higher-order density variations.

A key outcome of this analysis is the establishment of a practical applicability criterion for the Boussinesq approximation. Because the approximation relies on linearizing the density field, the relevant nondimensional parameter is the relative density change,

$${\displaystyle \frac{\Delta \rho }{{\rho }_{\mathrm{ref}}}}=\beta \Delta T,$$

(23)

which corresponds to the first-order term of the Taylor expansion of $\rho \left(T\right)$. As shown in Figure 13, noticeable deviations from the full variable-density solution emerge once $\beta \Delta T\gtrsim 0.1$, corresponding to a relative density variation of approximately 10%. This threshold is consistent with common engineering practice, in which density variations exceeding about 10% are generally considered beyond the valid range of the Boussinesq approximation. Beyond this limit, higher-order density effects become increasingly important and a fully variable-density formulation is required.

The results indicate that the Boussinesq approximation remains accurate for $\beta \Delta T\lesssim 0.1$, corresponding to relative density variations below approximately 10%. Although this threshold was identified using chloride-based molten salt properties, the criterion itself is fundamentally fluid-independent. Therefore, the applicability limit obtained in this study is relevant not only for molten salts but also for other fluids undergoing buoyancy-driven natural convection.

5. Conclusions

This work presented the development and verification of CUPID-MSR, an extended version of the CUPID thermal-hydraulics code incorporating the thermophysical properties of KCl–${\mathrm{UCl}}_3$ and NaCl–${\mathrm{MgCl}}_2$–${\mathrm{TRUCl}}_3$ molten salt mixtures. Temperature-dependent correlations for density, specific heat capacity, thermal conductivity, and viscosity were implemented to represent realistic chloride-based molten salt behavior.

The code was verified against the de Vahl Davis natural-convection benchmark over Rayleigh numbers ranging from ${10}^3$ to ${10}^6$. The predicted average Nusselt numbers agreed with the reference benchmark within 0.4–3.9%, demonstrating the correctness of the numerical formulation, material property implementation, and boundary-condition treatment. The benchmark simulations also reproduced the characteristic flow structures reported in the literature, further confirming the physical fidelity of the model. Additional sensitivity studies showed that variations in reference temperature and temperature-dependent thermophysical properties do not degrade numerical stability or accuracy, indicating that CUPID-MSR can reliably capture natural-convection behavior over a wide thermal range.

In addition to code verification, this study evaluated the applicability of the Boussinesq approximation by comparing it with the full variable-density formulation over a wide range of temperature differences. Using the nondimensional parameter $\beta \Delta T$, which corresponds to the relative density change $\Delta \rho /{\rho }_{\mathrm{ref}}$, a practical validity limit was identified. The results indicate that the Boussinesq approximation remains accurate for $\beta \Delta T\lesssim 0.1$, while larger density variations lead to noticeable deviations that require a fully variable-density treatment. Although demonstrated using molten salt correlations, this criterion depends only on the normalized density variation and is therefore applicable to buoyancy-driven natural-convection flows in general Newtonian fluids.

Overall, this study verifies the accuracy and robustness of CUPID-MSR for modeling buoyancy-driven flows and establishes a general, fluid-independent guideline for determining when the Boussinesq approximation may be used in place of a full variable-density formulation.