© H. Mistry et al., Published by EDP Sciences, 2026

1. Introduction

As the global demand for sustainable and high-efficiency energy sources increases, advanced nuclear reactor designs are receiving renewed interest. Among these, Liquid Metal Fast Reactors (LMFRs)–such as the Sodium Fast Reactor (SFR) and the Lead Fast Reactor (LFR)–stand out as promising Generation IV nuclear systems due to their excellent neutron economy, closed fuel cycle potential, and superior thermal properties [1, 2]. The use of liquid metals as primary coolants in these reactors offers several advantages, including high thermal conductivity, low neutron moderation, and high boiling points, enabling efficient heat removal at near-atmospheric pressure [3].

However, despite these thermophysical advantages, the prediction of the turbulent heat transfer in LMFRs remains a significant challenge, primarily due to the very low-Prandtl numbers of liquid metals–typically ranging between 0.001 and 0.025 under reactor operating conditions [4]. In such regimes, the disparity between thermal and momentum diffusivities leads to thick thermal boundary layers, anisotropic turbulence, and heat transfer mechanisms that are not well captured by conventional models designed for unity-Pr fluids like water or air [5]. Due to the high thermal conductivity and low viscosity of lead, sodium, etc., the near-wall temperature shear layers differ significantly from the momentum boundary layers. This can lead to inadequate heat flux, temperature and velocity profiles in analyses with heat transfer from or to the wall. This is especially critical in transition regimes between convection-dominated and conduction flow, often encountered in reactor cores and internal components.

Conventional Computational Fluid Dynamics (CFD) simulations of turbulent heat transfer commonly rely on the Simple Gradient Diffusion Hypothesis (SGDH), which assumes a linear relationship between the turbulent heat flux and the mean temperature gradient via a turbulent Prandtl number [6]. While this approach yields acceptable predictions for engineering applications involving unity-Prandtl fluids, its limitations become increasingly apparent in low Prandtl number flows, where the assumptions of isotropy and flux-gradient alignment no longer hold [7–9]. These limitations are particularly significant in mixed and natural convection regimes, as noted in the OECD/NEA benchmark studies [4] and highlighted in [10].

Generally, the two equations model is based on eddy viscosity, eddy diffusivity and the Reynolds analogy [11]. In such models, the turbulent momentum transport and the turbulent heat transport are assumed to be proportional [12], which is not true for liquid metal flows, which exhibit higher thermal diffusivity compared to the molecular viscosity [13, 14]. To overcome the limitations of eddy diffusivity-based models, Algebraic Heat Flux Models (AHFMs) have emerged as an alternative that provides more physically accurate results. AHFMs are derived from second-order moment closures and account for anisotropic turbulence effects, mean velocity gradients, buoyancy, and other non-local interactions [15, 16]. One notable formulation was proposed by Kenjereš [17], incorporating these terms in an algebraic frame-work suitable for engineering CFD. However, the original model coefficients were calibrated primarily for flows with Pr ≈ 1, limiting their application to liquid metal scenarios.

Further implementation and validation of AHFM model is shown by Shams [18] for commercial software STAR-CCM+. An important step forward was taken by Shams [9], who recalibrated the AHFM coefficients for forced convection at low Prandtl numbers. Their formulation introduced a logarithmic dependency on the global Reynolds and Prandtl numbers. While this extension broadened the model’s theoretical applicability, it also introduced a significant limitation: the need for defining a global Reynolds number and characteristic length scale, which may not be well-defined or spatially consistent in complex geometries such as reactor cores, heat exchangers, or flow around fuel assemblies. Moreover, the requirement for a priori knowledge of bulk flow quantities poses challenges in transient simulations, where these quantities evolve over time [9].

Another attempt to improve the simulation of low-Prandtl number fluids with the help of AHFM can be found in the work of Habiyaremye [19]. There, a new AHFM formulation is presented, in which a local turbulent Reynolds number and the fluid’s molecular Prandtl number are used to calculate the model’s coefficient that is associated with forced convection. The coefficient was calibrated in two different forced convection cases by minimizing the error in the temperature profiles compared to reference DNS data at low Prandtl numbers (i.e. Pr < 0.2). Two additional forced convection cases were used to validate the new coefficient formulation. It was shown that in all geometries, the new formulation performs better than the commonly used SGDH at low Prandtl numbers.

To improve the modelling of the turbulent heat flux in liquid metal flows, Manservisi [20] decided to use a four parameters k − ε − k_θ − ε_θ turbulence model, which was then tested for the simulation of basic cases. The investigation of the heat transfer models based on variable turbulent Prandtl number for liquid metal flows is performed by using both momentum and energy transport turbulence equations.

Further attempt to deal with turbulent heat flux in liquid metal flows was done at KIT, Germany. It is based on look-up tables and the determination of a local variable turbulent Prandtl number from DNS data. The implementation of this approach depends on the local Prandtl number, the global Reynolds number, and the y+ value, which is estimated from a pre-calculation. Although, this approach is robust and performed reasonably well for pipe and rod bundle flows, buoyant, natural convection flows could not be predicted well, and the intermittency effects were not considered. The advantages of the look-up tables are robustness and the fact that additional equations are not required. The drawbacks are mainly related to the applied isotropic Reynolds analogy and the dependency of the implementation on the wall spacing.

In order to be able to accurately map the behaviour of fluids with a low Prandtl number with OpenFOAM, first implementation at GRS of an AHFM turbulence model in OpenFOAM was done within the RS1547 [21] research project. The available OpenFOAM incompressible non-linear turbulence model LienCubicKE [22, 23] was extended with an AHFM implementation called LienCubicKEAHFM and released within OpenFOAM Foundation Version 6. Also, a special solver (buoyantBoussinesqPimpleFoamAHFM) was developed for simulations, in which the calculation of the buoyancy is based on the Boussinesq approximation. However, this specific solver has not been further pursued after OpenFOAM Version 7, which was released in 2019 [24]. Its successor, the buoyantPimpleFoam solver, does not solve the transport equation for the temperature, but rather for the energy. In addition to that, OpenFOAM Foundation has made massive changes to the code structure, which lead to a high implementation effort of the AHFM in the new version, as well as duplication of some Open-FOAM internal libraries. Therefore, another solution was necessary.

In the present work, the developed turbulent heat flux model is extended on the basis of the existing material libraries and made usable with the basic solvers (fluid and incompressible Fluid) [25]. The turbulence model is formulated to be used for compressible and incompressible flows. In addition, a new turbulence model based on Non-Linear Eddy Viscosity model (NLEVM) has been implemented. The developments were implemented with the kind support of HZDR in their NUSAR-RCS code that is based on OpenFOAM Foundation version 12-s.1 [26]. The new development is verified and then validated against DNS data for natural convection in a slender geometry.

2. AHFM implementation in OpenFOAM

This section provides short details on the turbulent heat flux theory as well as on the AHFM implementation strategy.

2.1. Turbulent heat flux

The Reynolds-averaged momentum and energy equations with turbulent flows for most industrial application are defined as [27]:

$$\begin{array}{c}\hfill \frac{D{U}_i}{\mathit{Dt}}=F_i-\frac{1}{\rho }\frac{\partial P}{\partial x_i}+\frac{\partial }{\partial x_j}(\nu \frac{\partial U_i}{\partial x_j}-\overline{u_i{u}_j})\end{array}$$(1)

$$\begin{array}{c}\hfill \frac{\mathit{DT}}{\mathit{Dt}}=\frac{q_i}{\rho C_p}+\frac{\partial }{\partial x_j}(\frac{\nu }{{\mathit{Pr}}_t}\frac{\partial T}{\partial x_j}-\overline{\theta u_j})\end{array}$$(2)

where, $\frac{D}{\mathit{Dt}}=\frac{\partial }{\partial t}+U_k\frac{\partial }{\partial x_k}$ represents the material derivative, F_i is the body force per unit mass, q_i is the internal energy source, and C_p is heat capacity of fluid. In these equations, $\overline{u_i{u}_j}$ denotes the turbulent momentum flux, while $\overline{{\theta u}_j}$ is turbulent heat flux correlation with fluctuating velocity–temperature. To obtain an adequate closure of the mathematical system, models are added to close the nonlinear terms.

The traditional definition of the turbulent heat flux in the energy equation is based on the Boussinesq assumption, which assumes that the heat flux is aligned with the temperature gradient and scaled by the turbulent eddy viscosity. This approximation fails in the case of a flow regime that deviates from the energy equilibrium, such as a buoyancy-dominant flow. Consequently, the primary objective of a heat flux model should be to deliver accurate and reliable results across all types of flow regimes [27]. With such a limitation, sophisticated turbulent heat flux modelling approaches have to be employed.

2.2. Algebraic heat flux model

Several AHFMs have been developed within the European THINS [28] and SESAME projects [29] such as an explicit and implicit formulation. In the explicit formulation, the turbulent heat flux is determined using the gradient hypothesis combined with the eddy diffusivity concept. This involves the definition of suitable turbulence timescales, which results in the development of four-parameter turbulence models [20]. Due to its isotropic nature, this approach has some limitations and requires further validation. In the implicit approach, the turbulent heat flux is obtained by directly solving a nonlinear algebraic equation derived from the transport equation for turbulent heat flux, using an algebraic solution method [17]. This formulation is a function of Reynolds stress anisotropy and temperature variance, for which an additional transport equation is solved. Accordingly, this implicit formulation is implemented in this study as described in [17].

The performance of the AHFM model is closely related to the correct approximation of the turbulent behaviour near the wall, and therefore, requires a low Reynolds number turbulence model. In general, the unknown Reynolds stress tensor or turbulent momentum flux $\overline{u_i{u}_j}$ is assumed to obey a linear constitutive relation based on the mean strain rate [15, 27].

$$\begin{array}{c}\hfill -\overline{u_i{u}_j}=-\frac{2}{3}{\delta }_{\mathit{ij}}k+2{\nu }_t{S}_{\mathit{ij}}\end{array}$$(3)

The turbulent heat flux $\overline{{\theta u}_i}$ is commonly represented by the simple eddy-diffusivity model, which is based on the Reynolds analogy defined as:

$$\begin{array}{c}\hfill \overline{\theta u_i}=-{\alpha }_t\frac{\mathit{dT}}{d{x}_i}=-\frac{{\nu }_t}{{\mathit{Pr}}_t}\frac{\mathit{dT}}{d{x}_i}\end{array}$$(4)

where, k is the turbulent kinetic energy, ν_t is the eddy viscosity, δ_ij is the Kronecker delta, α_t is the turbulent thermal diffusivity, and Pr_t is the turbulent Prandtl number (constant value of 0.9) [27].

In AHFM, the algebraic equation of turbulent heat flux proposed by Kenjereš [17] is implicitly solved using an algebraic method,

$$\begin{array}{cc}\hfill \overline{\theta u_i}& =-C_{t0}\frac{k}{\epsilon }(C_{t1}\overline{u_i{u}_j}\frac{\partial T}{\partial x_j}+C_{t2}\overline{\theta u_j}\frac{\partial U_i}{\partial x_j}+C_{t3}\beta g_i\overline{{\theta }^2})\hfill \\ \hfill & +C_{t4}{a}_{\mathit{ij}}\overline{\theta u_j}\hfill \end{array}$$(5)

where, ε is turbulent kinetic energy dissipation rate, a_ij is the Reynolds stress anisotropy tensor defined as $a_{\mathit{ij}}=\frac{\overline{u_i{u}_j}}{k}-\frac{2}{3}{\delta }_{\mathit{ij}}$. The terms appearing in equation (5) with C_t1 and C_t2 model coefficients describe the production of turbulent heat flux due to non-uniform mean temperature and velocity field, respectively. The term with C_t3 model coefficient is amplification of turbulence fluctuations due to the effect of buoyancy, and the last term with C_t4 model coefficient is redistribution of turbulent heat flux due to anisotropic turbulence [15, 30]. On the other hand, the gravitational acceleration g_i, the thermal expansion coefficient β, and the turbulent fluctuations in temperature as temperature variance $\overline{{\theta }^2}$ [9, 18], which are modelled with a separate field equation, are used to solve the turbulent heat flow $\overline{{\theta u}_i}$.

$$\begin{array}{cc}& \frac{\partial (\rho \overline{{\theta }^2})}{\partial t}+U_k\frac{\partial (\rho \overline{{\theta }^2})}{\partial x_k}\hfill \\ \hfill & =2P_t-2\rho {\epsilon }_{\theta }+\frac{\partial }{\partial x_i}\left((\frac{\kappa }{C_p}+\frac{{\mu }_t}{{\sigma }_{\overline{{\theta }^2}}})\frac{\partial \overline{{\theta }^2}}{\partial x_i}\right)\hfill \end{array}$$(6)

where, κ is thermal conductivity of fluid, μ_t is the dynamic turbulent viscosity, $P_t=-\rho \overline{\theta u_i}\frac{\partial T}{\partial x_i}$ is the production of temperature variance, and ε_θ is its dissipation rate. The most generic approach for solving the above equation is when the dissipation rate ε_θ is solved with an additional transport equation. Velocity and the thermal fields in liquid metal flows are characterized by different length and time-scales. In the thermal boundary layer, the time-scales differ up to two orders of magnitude [27]. As a result, this type of modelling contains more than two free parameters including turbulent time-scales. Thus, most models stick with a simpler approach and assume a constant thermal $({\tau }_{\theta }=\overline{{\theta }^2}/2{\epsilon }_{\theta })$ to mechanical time-scale (τ_m = k/ε) ratio R:

$$\begin{array}{c}\hfill R=\frac{{\tau }_{\theta }}{{\tau }_m}=0.5\end{array}$$(7)

The complete AHFM model is expressed in three-equation terms of k, ε, and $\overline{{\theta }^2}$. The k and ε turbulence equations will be discussed in following section.

2.3. Turbulence modelling

As mentioned earlier, the Boussinesq linear stress-strain relation is unable to accurately simulate anisotropic turbulence in real-world flow scenarios, because it ignores variations in the normal stresses. By taking into account non-linear correlations between Reynolds stresses and strain rates, nonlinear models attempt to overcome this limitation. These NLEVMs belong to an intermediate class of turbulence models. In these, the Reynolds stress is explicitly expressed algebraically in terms of vorticity tensors and the strain rate. Depending on the accuracy level, theoretically functional term forms as quadratic or cubic model. In OpenFOAM, these terms have a simplified general stress-strain formulation following [23]:

$$\begin{array}{cc}\hfill a_{\mathit{ij}}& =-2{C_{\mu }{f}_{\mu }S}_{\mathit{ij}}+[{\beta }_1(S_{\mathit{ik}}{S}_{\mathit{kj}}-\frac{1}{3}{S}_{\mathit{lk}}{S}_{\mathit{lk}}{\delta }_{\mathit{ij}})\hfill \\ \hfill & +{\beta }_2(S_{\mathit{ik}}{\mathrm{\Omega }}_{\mathit{kj}}-{\mathrm{\Omega }}_{\mathit{ki}}{S}_{\mathit{jk}})\hfill \\ \hfill & +{\beta }_3(\frac{1}{2}({\mathrm{\Omega }}_{\mathit{ik}}{\mathrm{\Omega }}_{\mathit{kj}}+{\mathrm{\Omega }}_{\mathit{jk}}{\mathrm{\Omega }}_{\mathit{ki}})-\frac{1}{3}{\mathrm{\Omega }}_{\mathit{kl}}{\mathrm{\Omega }}_{\mathit{lk}}{\delta }_{\mathit{ij}})\hfill \\ \hfill & -({\gamma }_1{S}_{\mathit{kl}}{S}_{\mathit{lk}}+{{\gamma }_2\mathrm{\Omega }}_{\mathit{kl}}{\mathrm{\Omega }}_{\mathit{lk}})S_{\mathit{ij}}-{\gamma }_3({\mathrm{\Omega }}_{\mathit{ik}}{\mathrm{\Omega }}_{\mathit{kl}}{S}_{\mathit{lj}}\hfill \\ \hfill & +S_{\mathit{ik}}{\mathrm{\Omega }}_{\mathit{kl}}{\mathrm{\Omega }}_{\mathit{lj}}-{\mathrm{\Omega }}_{\mathit{kl}}{\mathrm{\Omega }}_{\mathit{lk}}{S}_{\mathit{ij}}-\frac{2}{3}{{\delta }_{\mathit{ij}}\mathrm{\Omega }}_{\mathit{km}}{S}_{\mathit{mn}}{\mathrm{\Omega }}_{\mathit{nk}})\hfill \\ \hfill & -{\gamma }_4(S_{\mathit{im}}{S}_{\mathit{mk}}{\mathrm{\Omega }}_{\mathit{kj}}-{\mathrm{\Omega }}_{\mathit{ik}}{S}_{\mathit{km}}{S}_{\mathit{mj}})]\hfill \end{array}$$(8)

In the above equation, the first line shows the linear form, and adding a term with (β₁, β₂, β₃) coefficients leads to quadratic products of mean strain and vorticity rate, S_ij and Ω_ij. Likewise, adding terms with γ₁, γ₂, γ₃, γ₄ coefficients generate cubic products of the strain and vorticity rate. The damping coefficient fμ, and other coefficients, i.e. (β and γ) support wall functions to model the flow with different mesh resolutions of the boundary layer at the wall. This formulation provides a better description of the anisotropic turbulent structures of streamline curvature in numerical flow predictions

$$\begin{array}{c}\hfill S_{\mathit{ij}}=\frac{1}{2}(\frac{\partial U_i}{\partial x_j}+\frac{\partial U_j}{\partial x_i}),{\mathrm{\Omega }}_{\mathit{ij}}=\frac{1}{2}(\frac{\partial U_i}{\partial x_j}-\frac{\partial U_j}{\partial x_i})\end{array}$$(9)

Furthermore, another important feature of the cubic models is the dependence of the eddy viscosity coefficient, C_μ. Usually, a constant C_μ (i.e. 0.09), which is used in conventional linear eddy viscosity models like the standard k − ε model, causes an excessive development of turbulence energy levels in impingement and stagnation zones. However, the cubic nonlinear model calculates lower turbulent kinetic energy in areas with large irrotational strains due to a functional dependency of C_μ based on strain and vorticity invariant [23], which is defined in Table 1. ν_t is generally defined as follows:

$$\begin{array}{c}\hfill {\nu }_t={f_{\mu }C}_{\mu }\frac{k^2}{\epsilon }\end{array}$$(10)

In this paper, four different turbulence models, implemented by GRS in OpenFOAM, are studied: LienCubicKEBase [22], SugaCubicKE, as well as the new developments LienCubicKEAHFM and SugaCubicKEAHFM. SugaCubicKE is originally developed by Suga [31], but GRS used the formulation of [32] to develop in the present work. All models are extended and implemented at GRS as general classes and can be used for compressible and incompressible flows.

2.3.1. Non-linear turbulence model

The standard LienCubicKE, which is based on an incompressible turbulence model, is formulated in this study to deal with compressible and incompressible fluids. This extended turbulence model is called LienCubicKEBase. For this, k and ε transport equations are multiplied by a non-uniform volume fraction field (alphaField) and a density field (rhoField), similar to other RAS (Reynolds Averaged Simulation) turbulence models in OpenFOAM.

$$\begin{array}{cc}& \frac{\partial ({\alpha }_f\rho k)}{\partial t}+U_k\frac{\partial \left({\alpha }_f\rho k\right)}{\partial x_k}-\frac{\partial }{\partial x_i}{\alpha }_f\left[(\frac{{\mu }_t}{{\sigma }_k}+\mu )\frac{\partial k}{\partial x_i}\right]\hfill \\ \hfill & ={\alpha }_f\rho G-{\alpha }_f\rho \frac{\epsilon }{k}k+S_k+S_{(fv)}\hfill \end{array}$$(11)

$$\begin{array}{cc}& \frac{\partial ({\alpha }_f\rho \epsilon )}{\partial t}+U_k\frac{\partial \left({\alpha }_f\rho \epsilon \right)}{\partial x_k}-\frac{\partial }{\partial x_i}{\alpha }_f\left[(\frac{{\mu }_t}{{\sigma }_{\epsilon }}+\mu )\frac{\partial \epsilon }{\partial x_i}\right]\hfill \\ \hfill & ={{\alpha }_f\rho C}_{\epsilon 1}G\frac{\epsilon }{k}-{{\alpha }_f\rho C}_{\epsilon 2}{f}_2\frac{{\epsilon }^2}{k}+{\alpha }_f\rho E+S_{\epsilon }+S_{(fv)}\hfill \end{array}$$(12)

where, f₂ is a damping function to allow high and low Reynolds number flow calculation, and G is turbulence production. Additionally, internal source terms (S_k and S_ε) are defined, and S_(fv) is the source term defined to support the correction of momentum transport (fvModels and fvConstraints). These source terms are modelled in both turbulence models.

SugaCubicKE model has optimized the coefficients appearing in the constitutive equation of the Suga model over a range of flows including simple shear, impinging, curved, and swirling flows [32]. Based on this, the transport equation for k is defined as follows [25]:

$$\begin{array}{cc}& \frac{\partial ({\alpha }_f\rho k)}{\partial t}+U_k\frac{\partial \left({\alpha }_f\rho k\right)}{\partial x_k}-\frac{\partial }{\partial x_i}{\alpha }_f\left[(\frac{{\mu }_t}{{\sigma }_k}+\mu )\frac{\partial k}{\partial x_i}\right]\hfill \\ \hfill & ={\alpha }_f\rho G-{\alpha }_f\rho \frac{(\epsilon +D)}{k}k+S_k+S_{(fv)}\hfill \end{array}$$(13)

This Suga model has a non-zero D term, which distinguishes homogeneous and inhomogeneous dissipation rates. In many references ε is written as $\tilde{\epsilon }$, but not in [23]. The ε transport equation with damping function is defined as follows:

$$\begin{array}{cc}& \frac{\partial ({\alpha }_f\rho \epsilon )}{\partial t}+U_k\frac{\partial \left({\alpha }_f\rho \epsilon \right)}{\partial x_k}-\frac{\partial }{\partial x_i}{\alpha }_f\left[(\frac{{\mu }_t}{{\sigma }_{\epsilon }}+\mu )\frac{\partial \epsilon }{\partial x_i}\right]\hfill \\ \hfill & ={{\alpha }_f\rho C}_{ϵ1}G\frac{\epsilon }{k}-{{\alpha }_f\rho C}_{ϵ2}{f}_2\frac{{\epsilon }^2}{k}+{\alpha }_f\rho (E+Y_c)+S_{ϵ}+S_{(fv)}\hfill \end{array}$$(14)

The additional near-wall source terms E and Y_c [33] are the limiters for the correct near-wall behaviour and turbulent length scale, respectively. In this model, the length scale is obtained from solving a dissipation rate equation that can define a local turbulent Reynolds number R_t without reference to position [34], which vanishes at solid boundaries. This approach leads to an improved prediction of laminarization, mixed convection and diffusion when using two-equation models. Here, we examine the differences between LienCubicKEBase and SugaCubicKE, see Table 1. The same model constants (listed in Tab. 2) are used in both turbulence models.

After implementation of both base classes (i.e., LienCubicKEBase and SugaCubicKE), LienCubicKEAHFM and SugaCubicKEAHFM turbulence models are defined as new subclasses in OpenFOAM, which are calling the constructor of their base class. In this new subclass, the contribution of turbulent production due to buoyancy is defined using internal kSource (S_k) and epsilonSource (S_ε) following [35], which are used in the transport equations implemented in the base classes:

$$\begin{array}{c}\hfill S_k=-{\alpha }_f\rho \beta \left(g_i\overline{\theta u_i}\right)\end{array}$$(15)

$$\begin{array}{c}\hfill S_{\epsilon }=C_{\epsilon 1}{C}_{\epsilon 3}{S}_k\left(\frac{\epsilon }{k}\right)\end{array}$$(16)

Additionally, in these new classes the algebraic equation (5) and the transport equation, for the temperature variance (6) are implemented. To solve these equations, some coefficients from the literature [18] are listed in Table 3. Using the algebraic method, equation (5) for $\overline{{\theta u}_i}$ is solved using $\overline{{\theta }^2}$. The values of $\overline{{\theta u}_i}$ are then used for the source term in the energy equation using the fvModels based mechanism. As discussed below, the energy equation (17) is equivalent to the temperature equation (2). When the AHFM is called, the transport equations k, ε and $\overline{{\theta }^2}$ are first solved, and consequently $\overline{{\theta u}_i}$.

This set of constants, as suggested in [18], was used for all cases presented in this work. It is based on extended experience of other scientists in this field. These might not always deliver the best agreement with experimental data, based on different geometries and boundary conditions. More research is necessary to find a set that can be of general purpose. In general, having multiple constants in the model is definitely a drawback. One possible way to get rid of this dependency is to introduce correlations for these constants. Such correlations might be based on dimensionless numbers like Pr, Re, and Ra, which might help to represent the local flow phenomena. The objective of the current study is not to optimize the constants in the turbulence model, but rather to explore the potential of AHFM to improve the simulation of liquid metal flows through advanced modelling of the turbulent heat flux.

2.3.2. Modelling of the energy correction term

In OpenFOAM, the transport equation is not solved for temperature, but rather for energy, which is implemented on the basis of total energy. This energy equation is based on enthalpy (h) or internal energy (e) for some solvers (e.g., compressible flow, heat transfer, multiphase flow), but without consideration of the mechanical source term. The simple energy equation is:

$$\begin{array}{cc}& \frac{\partial \rho (h\Vert e)}{\partial t}+\mathrm{\nabla }·\left(\rho U_i(h\Vert e)\right)\hfill \\ \hfill & +\frac{\partial \rho K}{\partial t}+\mathrm{\nabla }·\left(\rho U_{i}K\right)+[\mathrm{\nabla }·\left(U_{i}p\right)\Vert (-\frac{\partial p}{\partial t})]\hfill \\ \hfill & =-\mathrm{\nabla }·q+S_{\mathtt{fvModels}}(\mathtt{h}\Vert \mathtt{e})+\rho U_i{g}_i\hfill \end{array}$$(17)

where, heat flux vector ∇ ⋅ q = −α_eff∇e and α_eff is effective thermal diffusivity. A new general energyAHFM class in fvModels is developed in order to have the contribution of turbulence to the heat transport $\overline{{\theta u}_i}$ for determining the energy or temperature equation. In fvModels, an explicit contribution of the AHFM $\overline{{\theta u}_i}$ is added, which is modelled for incompressible, compressible and multiphase flow respectively. This implementation is only intended for use with the described non-linear turbulence model, while the Pr_t value in thermophysicalProperties is set to a very high value to override alphat contribution in α_eff (α_eff = alpha + alphat).

$$\begin{array}{cc}& S_{\mathtt{fvModels}}=-\frac{\partial }{\partial x_i}\left({\alpha }_f\rho C_p\overline{\theta u_i}\right)\hfill \end{array}$$(18)

$$\begin{array}{cc}& S_{\mathtt{fvModels}}=-\frac{\partial }{\partial x_i}\left(\overline{\theta u_i}\right).\hfill \end{array}$$(19)

As mentioned above, the temperature equation in OpenFOAM is no longer available for incompressible flows. Should incompressible flow be considered, then a passive scalar transport equation (T as temperature) is used in OpenFOAM. Using equation (19), the temperature equation is solved with user defined laminar and turbulent diffusivity using coefficients α_l and α_t.

$$\begin{array}{cc}& \frac{\partial T}{\partial t}+U_k\frac{\partial T}{\partial x_k}-\frac{\partial }{\partial x_i}\left[({\alpha }_l\nu +{\alpha }_t{\nu }_t)\frac{\partial T}{\partial x_i}\right]\hfill \\ \hfill & =-\mathrm{\nabla }·q+S_{\mathtt{fvModels}}\hfill \end{array}$$(20)

To close the incompressible flow momentum equation, a new fvModels class incompressibleBuoyancyForce is programmed, which calculates and applies the buoyancy force considering the Boussinesq hypothesis to the incompressible momentum equation corresponding to the specified incompressible flow velocity field:

$$\begin{array}{c}\hfill S_{\mathtt{fvModels}}=-\left((T-T_0)\beta \right)g_i\end{array}$$(21)

where, T₀ is the reference temperature.

3. Verification and validation of the new turbulence models with reference data

For the validation of the turbulence model, two test cases are considered, namely, the Turbulent Channel Flow and the Rayleigh–Bénard convection (RBC) cases. The Turbulent Channel Flow case provides a basis for comparing the flow field in a simple channel, while RBC features natural convection phenomena. In the present study, steady state, two-dimensional RANS (Reynolds-averaged Navier-Stokes) calculations are performed using a second order upwind scheme combined with the SIMPLE (Semi-Implicit Method for Pressure-Linked Equations) algorithm. It is worth mentioning that the methods applied here for flow regimes (i.e. liquid metal Pr ≪ 1) require accurate resolution of the boundary layer. Therefore, the first cell size is maintained at y⁺ <  1 for all simulations to ensure adequate near-wall resolution.

3.1. Turbulent channel flow case

In this verification case the flow is assumed to be fully developed and with fixed boundary conditions. For this, an infinite horizontal channel flow has been considered. Periodic boundary conditions are applied for the development of the momentum boundary layer. The computational domain is 2D (6.4 * δ × 2 * δ) with thickness of one cell in the third dimension, where δ = 1, as shown in Figure 1. The domain length is sufficient to avoid any influence of the inlet turbulence. Here Re_τ = 180 is considered, where the Reynolds number is based on friction velocity and a Prandtl number of 0.7 is selected. The mesh consists of 38400 hexahedral cells (160 × 240) in total.

Fig. 1. Sketch of turbulent channel flow case.

The coordinates and flow variables are normalized by the channel half-width δ, the friction velocity u_τ, and the kinematic viscosity ν. For this case, the extended models for compressible flows, basic LienCubicKEBase and SugaCubicKE turbulence models are compared with DNS data [36] as well as with the analytical calculation of [37] of the dimensionless velocity and the turbulent viscosity profile (see in Figs. 2 and 3).

Fig. 2. Evolution of the dimensionless velocity profile.

Fig. 3. Evolution of dimensionless turbulent viscosity profile.

In the region y⁺ ≈ 10–100 (log-layer), both models closely match the DNS data for the dimensionless velocity. Minor deviations appear at higher y⁺ values, where both developed models slightly overestimate the DNS values. The profile of the dimensionless turbulent kinematic viscosity also agrees well with DNS data in the viscous sublayer. In the outer region (y⁺ >  80), the LienCubicKEBase model tends to overestimate the turbulent viscosity compared to the DNS. In contrast, the SugaCubicKE model shows better agreement with the DNS trend, although it exhibits a small peak that slightly exceeds the DNS values.

3.2. Rayleigh–Bénard convection case

The RBC case is a simple configuration for studying heat transfer via natural convection. It is provided by a limitless fluid layer enclosed by two fixed, horizontal, isothermal walls. The upper one is cooled, while the lower one is heated. The physical problem is characterized by the Rayleigh (Ra) and Prandtl number (Pr):

$$\begin{array}{c}\hfill Ra=\frac{g_i\beta \mathrm{\Delta }T{D}^3}{\nu \alpha }=Gr·Pr\end{array}$$(22)

$$\begin{array}{c}\hfill Pr=\frac{\nu }{\alpha }\end{array}$$(23)

Fig. 4. Computational domain for the RBC case.

Fig. 5. Profiles of normalised temperature for RBC-I.

Fig. 6. Profiles of normalised temperature variance for RBC-I.

Fig. 7. Profile of normalised turbulent heat flux for RBC-I.

where, ΔT is wall temperature difference, D is the distance between two horizontal walls, α is the thermal diffusivity, and Gr is the Grashof number. The RBC calculations are performed on a slender geometry with a 1:8 aspect ratio, assuming constant fluid properties using Boussinesq model. Figure 4 illustrates the computational domain of the RBC case. The 2D mesh consists of 3240 hexahedral cells (82 × 41) in total. Mesh analysis was performed with different mesh sizes: 40 × 20, 82 × 41 and 164 × 82 elements. The variation in the temperature field was very low. Therefore, the presented results shall be considered as mesh independent. Two vertical walls are treated as periodic boundaries in the horizontal direction, while the lower and upper walls are defined with no-slip conditions and maintained at constant wall temperatures, where (T_h − T_c)≥1.

Fig. 8. Profiles of normalised temperature for RBC-II.

Fig. 9. Profiles of normalised temperature variance for RBC-II.

In this study, the LienCubicKEAHFM and SugaCubicKEAHFM turbulence models are compared with DNS data and the solution from [9], referred here as AHFM-2005. The selected test cases are summarised in Table 4. For all the cases, the plot of T+ profile represents a normalised mean temperature defined as T + =(T−T_min)/(T_max−T_min), the temperature variance is expressed as normalised temperature variance $\overline{{\theta }^2}/{\mathrm{\Delta }T}^2$, and the normalised turbulent heat flux is defined as $\overline{\theta u_i}/\mathrm{\Delta }T$. The channel height is represented in terms of the non-dimensional, wall-normal coordinate y/H. For the comparison of heat transfer, the Nusselt number (Nu) is derived using the following equation, where h_c is the heat transfer coefficient and L is the characteristic length:

$$\begin{array}{c}\hfill Nu=\frac{h_{c}L}{\kappa }·\end{array}$$(24)

Fig. 10. Profile of normalised turbulent heat flux for RBC-II.

Fig. 11. Profile of normalised temperature for RBC-III.

The simulation case RBC-I is characterised by Pr = 0.7 and Ra = 6.3 × 10⁵. The normalised mean temperature T+ profiles Figure 5 show good agreement between the turbulence models and DNS results. Both LienCubicKEAHFM and SugaCubicKEAHFM accurately capture the steep gradients near the thermal boundary layers and the linear behaviour in the bulk of the channel. Their performance compares well to DNS and is slightly better than AHFM-2005, although the DNS profile in the bulk is flatter than all numerical predictions. Figure 6 shows the thermal fluctuation behaviour represented by $\overline{{\theta }^2}/{\mathrm{\Delta }T}^2$. None of the models can quantitatively reproduce the DNS data, although AHFM-2005 and SugaCubicKEAHFM more or less follow the general profile trend. While SugaCubicKEAHFM is very close to the DNS data near the walls and captures the overall trend of the profile, it overestimates the temperature variance in the bulk by some 40%. The results with AHFM-2005 are better in the bulk region, while the peaks are underestimated significantly (approx. 35%). The performance of LienCubicKEAHFM is further away from the DNS data, in both, in the bulk and the near-wall regions and does not reproduce the peaks near the walls, which are present in the other three curves.

Fig. 12. Profiles of normalised temperature variance for RBC-III.

Fig. 13. Profile of normalised turbulent heat flux for RBC-III.

The DNS data in of the normalised turbulent heat flux $\overline{\theta u_i}/\mathrm{\Delta }T$ (Fig. 7) shows a flat profile in the bulk and strong decay near the walls. LienCubicKEAHFM achieves the closest match from all simulations across the domain with a maximum deviation of 12% in the middle of the channel and with accurate prediction in the central and near-wall regions. SugaCubicKEAHFM predicts peaks in the near-wall regions and in the bulk, which are not in-line with the relatively flat DNS profile.

Fig. 14. T+ scalar field over the CFD domain calculated with LienCubicKEAHFM for Pr = 0.7 (top), Pr = 0.025 (middle) and Pr = 0.006 (bottom).

The simulation case RBC-II is defined by Pr = 0.025 and Ra = 1 × 10⁵. In Figure 8, all models predict a somewhat flatter profile compared to the DNS one. Although the calculated trends are similar, AHFM-2005 and SugaCubicKEAHFM are closer to DNS. LienCubicKEAHFM predicts well in the central region, while slightly over- and underestimating the normalised temperature when moving towards bottom and top walls.

In the plot of $\overline{{\theta }^2}/{\mathrm{\Delta }T}^2$ (Fig. 9) and $\overline{\theta u_i}/\mathrm{\Delta }T$ (Fig. 10), LienCubicKEAHFM underestimates by 17% the variance and overestimates by 25% the turbulent heat flux, but still qualitatively accurately tracks the DNS profiles. SugaCubicKEAHFM and especially AHFM-2005 significantly overpredict the variance throughout the domain, with maximal deviations in the center of the channel. It was not easy to reach tight convergence in the RBC-II case with the developed models, specially LienCubicKEAHFM, which might explain the unexpected asymmetry in the profiles.

Figure 11 illustrates the profile of the normalised temperature for the RBC-III case with Pr = 0.006 and Ra = 2.4 × 10⁴. Here, all models including LienCubicKEAHFM and SugaCubicKEAHFM closely follow the DNS trend across the full height, accurately reproducing the nearly linear temperature gradient. In contrast, the normalized temperature variance throughout the domain is not reproduced by all models (Fig. 12). SugaCubicKEAHFM and AHFM-2005 models qualitatively represent the DNS profile shape, but significantly over predict the maximum value of the fluctuations. The LienCubicKEAHFM prediction of the temperature variance in the centre of the geometry is close to the DNS value, but unfortunately, the shape of the profile is flatter, and thus, differs from the reference one. Figure 13 presents the normalised turbulent heat flux, calculated with LienCubicKEAHFM and SugaCubicKEAHFM, DNS and AHFM-2005 were not available to compare with. The comparison of the two profiles is similar to the ones in RBC-I and RBC-II cases.

Fig. 15. Gradient of temperature field over the CFD domain calculated with LienCubicKEAHFM for Pr = 0.7 (top), Pr = 0.025 (middle) and Pr = 0.006 (bottom).

Fig. 16. Mean velocity distribution over the CFD domain calculated with LienCubicKEAHFM for Pr = 0.7 (top), Pr = 0.025 (middle) and Pr = 0.006 (bottom).

Figure 14 presents the T+ field distribution, calculated with the LienCubicKEAHFM turbulence model for the three RBC cases. This figure is nearly identical to the one with the SugaCubicKEAHFM results (not shown here). In the case of RBC-I (Pr = 0.7), thermal plumes clearly rise and fall, which is typical for convection-dominated flows, where thermal diffusivity and momentum diffusivity are comparable (indication of strong convective mixing). For this case, the calculated Nusselt number is 6.7, which deviates by approximately 12% from the empirical solution in the literature [15]: Niemela – 7.53 (Nu = 0.124Ra^0.309), Grossmann – 7.68 (Nu = 0.1Ra^1/4 + 0.05Ra^1/3). These correlations are valid for Prandtl number 0.7 and Rayleigh number between 10⁶ ≤ Ra ≤ 10¹⁷.

As the Prandtl number decreases, the thermal boundary layer gets thicker. The Rayleigh number is also insufficient to drive strong convection, which results in isothermal layers that are mostly horizontal, indicating the dominance of the conduction with stratified temperature gradient (also visible in Fig. 15). This behaviour closely follows the DNS profile seen in Figure 11. The calculated Nusselt number for RBC-II and RBC-III are 1.9 and 1.19, which are lower than the result of the empirical solution of Rossby [38] – 2.83 (Nu = 0.147Ra^0.257, Pr = 0.025, 10³ ≤ Ra ≤ 5 × 10⁵), Kek [38] – 1.45 (Nu = 0.117Ra^0.25, Pr = 0.006, 4 × 10⁴ ≤ Ra ≤ 2.5 × 10⁵). The possible reason for the lower Nu in our calculation might be the choice of AHFM model coefficients and parameters.

Figure 16 presents the velocity fields for RBC cases and illustrate how the flow structure changes with varying Ra and Pr numbers. In RBC-I (Ra = 6.3 × 105, Pr = 0.7), the velocity distribution shows strong, coherent convection rolls, with smoothly varying regions of higher velocity concentrated near the centers of the convection cells (i.e. Fig. 14, RBC-I case). Because momentum and thermal diffusion occur on similar time scales at moderate Pr, the flow organizes into stable large-scale circulation patterns with velocity magnitudes on the order of 10⁻¹. In RBC-II (Ra = 1 × 10⁵, Pr = 0.025), the dramatic drop in Prandtl number produces a very different behaviour. Thermal diffusion becomes dominant, weakening buoyancy-driven motion and generating thin, wavy, and irregular velocity structures (transitional or chaotic) rather than well-formed rolls. Velocities also decrease significantly, in the order of 10⁻⁵ to 10⁻⁴. RBC-III (Ra = 2.4 × 10⁴, Pr = 0.006) shows an even more diffusive (like laminar flow) and weakly convective regime, where both the lower Pr and lower Ra reduce flow intensity. The velocity field becomes smoother and more horizontally stretched, with minimal evidence of convection cells (smeared out), indicating a regime approaching conduction-dominated behaviour. Overall, the comparison highlights how decreasing Pr reduces the coherence of convection structures and how decreasing Ra further suppresses convective motion, resulting in progressively weaker and more diffused velocity fields across the three cases. The presented velocity distribution may vary depending on the selection of the governing parameters, see Table 3.

4. Conclusion and outlook

The commonly used eddy diffusivity approach has well-known limitations, which become particularly evident in advanced nuclear reactor applications that rely on liquid metal coolants such as sodium- and lead-cooled fast reactors. Deficiencies occur primarily in representation of natural and mixed convection regimes, being of particular importance in transient and accident conditions. The previously developed AHFM model at GRS was extended for compressible and incompressible flows within the framework of non-linear (cubic) turbulence models. The new developments were implemented in the HZDR’s NUSAR-RCS code, which is based on OpenFOAM.

It was shown that the predictions with the new models generally provide improved results compared to the ones, reported in the literature (AHFM-2005), but this is not always true for all investigated quantities. Deviations from DNS data in the shape and the magnitude of the evaluated profiles prove that further investigations are necessary. These shall aim at clarification of their origin, and further focus shall be put in the analysis of the different predictions of both implemented cubic k − ε turbulence models. Based on the discussed results, it is not easy to draw a final conclusion regarding the superiority of one of these models. Craft [40] proposed further correction in Suga’s model, which is currently under investigation at GRS. The analysis of the results may shed more light on the origin of deviations presented here.

The presented results are based on small scale academic test cases. In real life reactors, the geometries are very complex and the possible flow regimes, for example in normal and off-normal operational conditions, are quite different. Further, these regimes are not always steady and the transition between them further complicates the analysis of liquid metal flow and heat transfer phenomena. Therefore, these turbulence models must be further analysed and validated for reactor relevant conditions, before these can be applied for reactor safety analysis.

Funding

This work was funded by the German Federal Ministry for the Environment, Climate Action, Nature Conservation and Nuclear Safety within the UMRS 1607 project.

Conflicts of interest

The authors declare that they have no competing interests to report.

Data availability statement

The full description of performed work will be made publicly available within a report on the GRS website: https://www.grs.de/.

Author contribution statement

Conceptualization, Data curation, Formal analysis, Funding acquisition, Investigation, Methodology, Project administration, Resources, Supervision, Validation, Verification, Visualization, Writing – original draft, Writing – review & editing: – H. Mistry, J. Herb, A. Papukchiev.

Acknowledgments

This work was supported by the German Federal Ministry for the Environment, Climate Action, Nature Conservation and Nuclear Safety. The authors also acknowledge the kind support of HZDR in the implementation of the turbulence models in the NUSAR-RCS code.

Glossary

Glossary

AHFM: (Algebraic Heat Flux Model)

CFD: (Computational Fluid Dynamics)

DNS: (Direct Numerical Simulation)

GRS: (Gesellschaft f

HZDR: (Helmholtz-Zentrum Dresden-Rossendorf e.V.)

LFR: (Lead Fast Reactor)

LMFRs: (Liquid Metal Fast Reactors)

NLEVM: (Non-Linear Eddy Viscosity Model)

NUSAR-RCS: (Nuclear Safety Repository for OpenFOAM Foundation Software for Reactor Cooling System)

OECD/NEA: (Organisation for Economic Cooperation and Development (OECD)/Nuclear Energy Agency)

PWR: (Pressurized Water nuclear Reactor)

RANS: (Reynolds-Averaged Navier-Stokes equations)

RAS: (Reynolds Averaged Simulation)

RBC: (Rayleigh–B

SESAME: (Simulations and Experiments for the Safety Assessment of MEtal cooled rea-tor)

SFR: (Sodium Fast Reactor)

SGDH: (Simple Gradient Diffusion Hypothesis)

SIMPLE: (Semi-Implicit Method for Pressure-Linked Equations)

THINS: (Thermal-hydraulics of Innovative Nuclear Systems)

Roman symbols

a_ij: (Reynolds stress anisotropy tensor, −)

α_f: (Non-uniform volume fraction field, −)

α_t: (Turbulent thermal diffusivity, m²/s)

C_p: (Heat capacity of fluid, J/kg⋅K)

E: (Calculates near-wall turbulence energy dissipation, m²/s³)

e: (Internal energy, kg ⋅ m²/s²)

F_i: (Body force per unit mass, m/s²)

f₂: (Damping function, −)

G: (Production of turbulence, m²/s³)

Gr: (Grashof number, −)

g_i: (Gravitational acceleration, m/s²)

H: (Distance between two horizontal walls, m)

h: (Enthalpy, kg ⋅ m²/s²)

h_c: (Heat transfer coefficient, W/m²K)

K: (Specific kinetic energy, m²/s²)

k: (Turbulence kinetic energy, m²/s²)

L: (Characteristic length, m)

Nu: (Nusselt number, −)

P_t: (Production of temperature variance, K²/s)

p: (Pressure, Pa)

Pr: (Prandtl number, −)

Pr_t: (Turbulent Prandtl number, −)

q: (Heat flux vector, W/m²)

q_i: (Internal Energy Source, W/m³)

R: (Time-scale ratio, τ_θ/τ_m)

Ra: (Rayleigh number, −)

Re: (Reynolds number, −)

Re_τ: (Reynolds number based on friction velocity, −)

R_t: (Turbulent Reynolds number, −)

S_fvModels: (Source term using fvModels, −)

$\overline{S}$ : (Non-dimensional strain rate, −)

S_ij: (Mean strain rate tensor, 1/s)

T: (Temperature, K)

T+: (Normalized mean temperature, −)

T₀: (Reference temperature, K)

t: (Time, s)

U_i: (Velocity component, m/s)

$\overline{u_i{u}_j}$ : (Turbulent momentum flux, m²/s²)

x, y, z: (Direction, m)

y: (Wall normal distance, m)

y/H: (Non-dimensional wall-normal coordinate, −)

y^*: (Dimensionless distance from the wall, −)

Y_c: (Yap length scale correction m²/s⁴)

Greek symbols

α: (Thermal diffusivity, m²/s)

β: (Thermal expansion coefficient, 1/K)

δ_ij: (Kronecker delta, −)

ΔT: (Temperature difference, K)

ε: (Turbulent kinetic energy dissipation rate, m²/s³)

ε_θ: (Dissipation rate of temperature variance, K²/s)

η: (Suga turbulence model coefficient, −)

$\overline{{\theta }^2}$: (Temperature variance, K²)

$\overline{{\theta }^2}/{\mathrm{\Delta }T}^2$: (Normalized temperature variance, −)

$\overline{{\theta u}_i}$: (Turbulent heat flux, K⋅m/s)

$\overline{\theta u}/\mathrm{\Delta }T$: (Normalised turbulent heat flux, m/s)

κ: (Thermal conductivity of fluid, W/m⋅K)

κ_c: (Constant coefficient of wall function calculation, 0.41)

μ: (Dynamic viscosity, Pa⋅s)

μ_t: (Turbulent dynamic viscosity, Pa⋅s)

ν: (Kinematic turbulent viscosity, m²/s)

ν_t: (Eddy viscosity, m²/s)

ρ: (Fluid density, kg/m³)

τ_m: (Mechanical time-scale, s)

τ_θ: (Thermal time-scale, s)

$\overline{\mathrm{\Omega }}$: (Non-dimensional vorticity rate, −)

Ω_ij: (Mean vorticity rate tensor, 1/s)

References

Appendix

A. Turbulence model: LienCubicKEBase

Class:

Foam::RASModels::LienCubicKEBase

Header location:

$FOAM_SRC/MomentumTransportModels/momentumTransportModels/RAS/LienCubicKEBase/

Purpose:

The LienCubicKEBase model is a reformulated extension of the LienCubicKE turbulence model, constructed to ensure consistent behaviour in both compressible and incompressible flows and implemented in OpenFOAM through density- and volume-fraction-weighted transport equations for k and ε.

Main Components:

• Model coefficients: Defined in LienCubicKEBase.C and adjustable via RASProperties.

• Turbulent viscosity calculation: Involves cubic terms of the strain-rate and vorticity tensors.

• solve() routine: Handles sequential solution of k and ε equations each iteration.

• Correct() loop: Updates turbulent viscosity and non-linear stress contributions.

Users call:

 In constant/momentumTransport
  RAS
  {
    RASModel    LienCubicKEBase;
    turbulence  on;
  }

B. fvModels: energyAHFM

Class:

Foam::fv::energyAHFM

Header location:

$FOAM_SRC/fvModels/derived/energyAHFM/

Purpose:

The Algebraic Heat Flux Model (AHFM) provides a source term for the energy or temperature equation in compressible and incompressible flow simulations. It is specifically intended to operate consistently with the LienCubicKEAHFM and SugaCubicKEAHFM momentum transport models.

Users call:

 In constant/fvModels
  energyAHFM1
  {
    type    energyAHFM;
    libs    ("libaddonfvModels.so");
  }

C. fvModels: incompressible Buoyancy Force

Class:

Foam::fv::incompressibleBuoyancyForce

Header location:

$FOAM_SRC/fvModels/derived/incompressibleBuoyancyForce/

Purpose:

Evaluates the buoyancy contribution according to the Boussinesq approximation, expressed as (1 − β(T − T₀)) and incorporates it into the incompressible momentum equation for the specified velocity field. Intended to work with incompressible turbulence models.

Users call:

 In constant/fvModels
 incompressibleBuoyancyForce1
 {
  type    incompressibleBuoyancyForce;
  libs    ("libaddonfvModels.so");
  U       U; // Name of the velocity field
  T       T; // Name of the temperature field
  beta    beta; // Name of the thermal expansion field
 }

All Tables

All Figures

	Fig. 1. Sketch of turbulent channel flow case.
In the text

	Fig. 2. Evolution of the dimensionless velocity profile.
In the text

	Fig. 3. Evolution of dimensionless turbulent viscosity profile.
In the text

	Fig. 4. Computational domain for the RBC case.
In the text

	Fig. 5. Profiles of normalised temperature for RBC-I.
In the text

	Fig. 6. Profiles of normalised temperature variance for RBC-I.
In the text

	Fig. 7. Profile of normalised turbulent heat flux for RBC-I.
In the text

	Fig. 8. Profiles of normalised temperature for RBC-II.
In the text

	Fig. 9. Profiles of normalised temperature variance for RBC-II.
In the text

	Fig. 10. Profile of normalised turbulent heat flux for RBC-II.
In the text

	Fig. 11. Profile of normalised temperature for RBC-III.
In the text

	Fig. 12. Profiles of normalised temperature variance for RBC-III.
In the text

	Fig. 13. Profile of normalised turbulent heat flux for RBC-III.
In the text

	Fig. 14. T+ scalar field over the CFD domain calculated with LienCubicKEAHFM for Pr = 0.7 (top), Pr = 0.025 (middle) and Pr = 0.006 (bottom).
In the text

	Fig. 15. Gradient of temperature field over the CFD domain calculated with LienCubicKEAHFM for Pr = 0.7 (top), Pr = 0.025 (middle) and Pr = 0.006 (bottom).
In the text

	Fig. 16. Mean velocity distribution over the CFD domain calculated with LienCubicKEAHFM for Pr = 0.7 (top), Pr = 0.025 (middle) and Pr = 0.006 (bottom).
In the text