Development of Explainable Data-Driven Turbulence Models with Application to Liquid Fuel Nuclear Reactors

: Liquid fuel nuclear reactors offer innovative possibilities in terms of nuclear reactor designs and passive safety systems. Molten Salts Reactors (MSRs) with a fast spectrum are a particular type of these reactors using liquid fuel. MSFRs often involve large open cavities in their core in which the liquid fuel circulates at a high speed to transport the heat generated by the nuclear reactions into the heat exchangers. This high-speed ﬂow yields a turbulent ﬁeld with large Reynolds numbers in the reactor core. Since the nuclear power, the neutron precursor’s transport and the thermal exchanges are strongly coupled in the MSFR’s core cavity, having accurate turbulent models for the liquid fuel ﬂow is necessary to avoid introducing signiﬁcant errors in the numerical simulations of these reactors. Nonetheless, high-accuracy simulations of the turbulent ﬂow ﬁeld in the reactor cavity of these reactors are usually prohibitively expensive in terms of computational resources, especially when performing multiphysics numerical calculations. Therefore, in this work, we propose a novel method using a modiﬁed genetic algorithm to optimize the calculation of the Reynolds Shear Stress Tensor (RST) used for turbulence modeling. The proposed optimization methodology is particularly suitable for advanced liquid fuel reactors such as the MSFRs since it allows the development of high-accuracy but still low-computational-cost turbulence models for the liquid fuel. We demonstrate the applicability of this approach by developing high accuracy Reynolds-Averaged Navier–Stokes (RANS) models (averaged ﬂow error less than 5%) for a low and a large aspect ratio in a Backward-Facing Step (BFS) section particularly challenging for RANS models. The newly developed turbulence models better capture the ﬂow ﬁeld after the boundary layer tipping, over the extent of the recirculation bubble, and near the boundary layer reattachment region in both BFS conﬁgurations. The main reason for these improvements is that the developed models better capture the ﬂow ﬁeld turbulent anisotropy in the bulk region of the BFS. Then, we illustrate the interest in using this turbulence modeling approach for the case of an MSFR by quantifying the impact of the turbulence modeling on the reactor key parameters.


Introduction
The technology of liquid nuclear fuels was originally developed in the 1960s with molten salt, aqueous solutions, fluoride, chloride, and liquid metal technologies [1]. Very important progress on the molten salt liquid fuel technology was achieved during the Molten Salt Reactor Experiment (MSRE) [2] developed at ORNL. Molten Salt Reactors (MSRs) experimental studies were practically halted during the following decades until a resurgence in the late 2000s [3]. Liquid fuel Molten Salt Reactors (MSRs) offer the possibility of maintaining an actinide burning and fuel breeding cycle. Moreover, they can operate with a small reactivity excess when using continuous reprocessing [4]. These are attractive qualities from safety and economics points of view. Moreover, innovative passive safety systems can be designed for these reactors thanks to the solidification/melting possibility of molten salt fuel [5] and its fission gas retention properties [6].
On the downside, having a liquid nuclear fuel greatly complexifies the MSRs modeling efforts and the safety studies due to the strong coupling between neutronics, thermalhydraulics and thermomechanics aspects and also the novel phenomena introduced by the circulating liquid nuclear fuel. To cite two examples: the neutron precursors advective transport and internal heat source in the flowing molten fuel salt. In standard solid fuel reactors, such as the Light Water Reactors (LWRs) or the Sodium Fast reactors (SFRs), the power production in the solid fuel is only loosely coupled to the external cooling of the fuel rods. Hence, errors in modeling the core coolant thermal hydraulics generally mildly affect the neutron flux and thus the fission power calculations [7]. This is not the case for liquid fuel reactors, in which errors arising from the thermal-hydraulics modeling will affect the predictions for the neutron precursors transport and the fuel density/Doppler effects and thus the nuclear cross-sections. Moreover, nuclear power and neutron precursors concentration will, in turn, affect the liquid fuel temperature field. When developing a multiphysics model for liquid nuclear fuels one should therefore target a similar level of accuracy for the determination of the molten fuel salt flow distribution and heat transfer on one side and for the neutron transport resolution on the other side. However, when modeling these reactors, significantly larger errors are usually introduced from the thermalhydraulics field rather than from the neutronics and thermomechanics ones [8,9]. Much of this thermal-hydraulics modeling error arises from the use of inaccurate models of the turbulent fields in the MSRs core cavity. As an example, the turbulent flow field generated at the entrance to the Molten Salt Fast Reactor (MSFR) [10] core cavity is depicted in Figure 1. The flow Reynolds number in the cavity of this reactor, computed using the mean vertical velocity and the diameter at the mid-height section of the core cavity, is about 4 × 10 6 and, hence, Direct Numerical Simulations (DNS) and even Large Eddy Simulations (LES) are usually not practical with multiphysics tools when performing transient calculations and considering the full geometry of the fuel circuit. Therefore, in this work, we propose a novel method for automatically developing high-fidelity and low-cost Reynolds-Averaged Navier-Stokes (RANS) turbulence models from flow data obtained either from experimental measurements or from higher fidelity thermal-hydraulics numerical simulations yielded by a super-computer. Originally RANS turbulence research focused on developing adapted transport models for modeling the effective diffusivity due to turbulence usually assuming the turbulent Boussinesq hypothesis [11]. Different turbulence surrogated fields were thus proposed, e.g., turbulent intensity, and then closure coefficients describing the transport of these fields were fitted from the experimental database of flows. Nonetheless, it was soon realized that the Boussinesq hypothesis was not accurate in many flow problems, e.g., mixing, boundary later detachment, etc., and research proceeded towards the development of algebraic Reynolds stress models [12,13]. In these models, an algebraic equation was determined for the components of the turbulent Reynolds Stress Tensor (RST), which provides the model the ability to capture the anisotropy and second-order diffusivity effects due to turbulence [14]. Nonetheless, these models are still very dependent on the hypothesis made by the algebraic closure equation and thus on the flow field for which the algebraic stress balance was calibrated. Therefore, algebraic stress models have issues in yielding high accuracy in complex, realistic geometries such as those encountered in the liquid-fueled reactors. Owing to the limitations of algebraic Reynolds stress models, it did not come as a surprise that the rise of machine learning applications brought attached the proposal of new techniques for developing data-driven turbulence models [15].
Most works in the machine-learned turbulence field have used different versions of artificial neural networks (subrogate models) to either fit the turbulent diffusivity or to model the components of the Reynolds Stress Tensor as a function of the resolved velocity and pressure fields [15][16][17][18][19]. Even though high accuracy has been achieved with this approach, neural networks present two key limitations. First, the turbulent model in the neural network cannot be directly interpreted or analyzed (i.e., they are non-physics-based models). Second, owing to the large numbers of parameters involved in the neural network, the fitted models yield large errors when extrapolated to different flow regimes than the ones used for the original fitting. An approach that solves these limitations has been recently proposed separately by Akolekar et al. [20], Gimenez and Bre [21], and the authors in this work [22], among others. In this approach, genetic algorithms are rather employed to optimize algebraic closure in the turbulence model or to build functions for the Reynolds Stress Tensor based on a closed set of tensorial bases. This method solves the limitations of machine-learned turbulence models developed with neural networks because it yields an expression for the Reynolds Stress Tensor that can be interpreted. Moreover, with adequate extrapolations, it can be analyzed using standard methods. Nonetheless, it introduces two new drawbacks. Firstly, contrary to neutral networks and due to the randomness of genetic algorithms, a different "optimal" turbulence model is generally obtained for the same dataset at each new evaluation of the genetic optimization. This generates uncertainty on the appropriateness of the models learned. Secondly, standard genetic algorithms frequently converge to local minimums. Hence, different tailored variations need to be introduced to the standard process for trying to enforce its convergence to a global minimum [23] and not just a local one. This modification is frequently not well adapted for turbulent model optimization.
In the present work, we improve upon our previous attempts using genetic algorithms by proposing a new method that solves these standard issues of genetic algorithms. This method is inspired by the multi-level grammar genetic programming technique recently developed for heterogenous problem optimization [24]. First, we solve the randomness issue in genetic algorithms by building a Pareto priority stack that orders models from best to worse fit and then we apply the genetic operations deterministically on this stack. A set of grammar rules is then defined to ensure that each genetic operation builds a mathematically correct turbulence model. The order in the operations to be applied in the stack is optimized at each iteration with a surrogate Bayesian optimization model. This later converges to global minima provided that a large enough exploration rate is included in the algorithms [25]. Thus, the problem of convergence to local minima is also eliminated, or at least, significantly reduced.
The overall approach of this novel method and the results are discussed in detail in the next sections. In Section 2, the numerical model for the incompressible turbulent flows used in this work is introduced. In Section 3, the new deterministic method for data-driven modeling of the turbulence stress tensor is described. In Section 4, the application of this model for two backward-facing steps with a small and a large aspect ratio is presented. The results from the current model are compared with other standard RANS and LES models.
Finally, in Section 5, the impact of turbulence modeling in other multiphysics coupled fields is evaluated for the case of the MSFR.

Numerical Model
The fluid flows considered in this work can be described via the incompressible Navier-Stokes equations, which read as follows [26]: where u = u(x, t) is the Eulerian flow velocity, ρ the density, p = p(x, t) the mechanical pressure, µ the molecular viscosity, x is the position and t the time. To avoid resolving all scales in the turbulent flow, it is customary to decompose the velocity and pressure into a resolved space-time averaged field and a non-resolved sub-grid scale field. For this purpose, a firstorder space-time average filter defined as φ( are the boundaries of the computational cell. In Large Eddy Simulations (LES), one generally relies on the implicit filter provided by the mesh, whereas in transient Reynolds-Averaged Navier-Stokes (RANS) simulations, the space dependency in the filter is dropped and one defines a box filter in time being a Heaviside function and ∆t significantly larger than the largest turbulent timescales. Applying the filtering operation to the Navier-Stokes equations we obtain [26]: where τ turb = ρ( uu − u u ) is the RST. This tensor must be modeled in both RANS and LES approaches. Hence, our proposed optimization method for τ turb is applicable to both modeling methods. However, LES modeling is frequently less complicated as more flow scales are resolved than in RANS modeling. For the latter, it is common practice to separate the isotropic from anisotropic part of this tensor. The former results in an equivalent turbulent pressure, whereas the latter in an equivalent turbulent shear stress rate tensor. Defining the turbulent kinetic energy as k = 0.5 u ·u , the decomposition of the Reynolds Stress Tensor reads as [26]: where I is the identity tensor, ν t is the turbulent viscosity, and T is the structural form for the anisotropic part of the Reynolds Stress Tensor (or the Reynolds Shear Stress tensor) that for simplicity we call here the turbulent anisotropy tensor. The tensorial base for modeling the anisotropic part (both ν t and T) can be expressed as a function of the averaged strain-rate tensor S = 0.5 ∇ u + ∇ u T and the averaged vorticity (or rotation rate) tensor Ω = 0.5 ∇ u − ∇ u T , i.e., T = T(S, Ω). Different turbulent models seek an approximation for ν t and T. Namely, some of the models that will be compared in this work are the improved k − (implemented as in [26]), k − ω (implemented as in [27]), and k−ω SST models (implemented as in [28] with curvature corrections proposed in [29]). These models determine the turbulent viscosity ν t as function of the turbulent kinetic energy and its non-specific and specific dissipation rates, respectively, named and ω, and model the structural form of the Reynolds Stress Tensor via the Boussinesq turbulent hypothesis, i.e., T = S. Additionally, we will examine the Reynolds Stress Transport model (implemented with the analytical closures in [30]), where a separate model is solved for each of the independent terms in T along with an equation for the turbulent dissipation. Regarding LES, we will also include in the comparison the multiple gradient LES model [31], which seeks a model ν t = ν t ( S ) and T = S, but requires a significantly finer mesh.
Each of these turbulent models generally involves a calibration process where several coefficients are fit to yield simple turbulent problems. However, even fitted models can yield inaccurate results if the structural approximation used to model the turbulent anisotropy tensor T is not physically correct. For instance, for recirculation problems, one would expect to have the rotational component Ω in the turbulent anisotropy tensor T if the full scale of the recirculation is not resolved in the numerical mesh. For this purpose, we developed a deterministic method to automatize the fitting of the structural form of the turbulent anisotropy tensor T. It is important to remark that this optimization model deals exclusively with the structural form of the anisotropy and thus allows us to retain standard models (i.e., simple models) for the turbulent viscosity. Since the coefficients of ν t for standard turbulence models have already been calibrated on a large database of problems, we choose not to change them since this would change the applicability of the calibrated model. For example, modifying the fitting coefficients in the standard k − model for calculating k, and ν t will lead to suboptimal predictions of the return to isotropy in turbulent shear flows and in the decay of turbulent kinetic energy in isotropic turbulence. These issues will ultimately limit the applicability of the fitted model in more complex flow configurations in which several turbulence production and destruction mechanisms are present but could be modified if needed.
In this work, we will use as a baseline model (i.e., reference for the comparisons) for calculating ν t the standard k − model proposed by Mohammadi et al. [26]. In this one, k and are solved via the following transport equations: and the closure model used for the turbulent viscosity is simply ν t = C µ k 2 with the coefficients presented in Table 1. However, contrary to the standard k − model, we will not assume T = S (Boussinesq turbulent hypothesis) but optimize this tensor to improve the accuracy. The proposed method for structural optimization of the turbulent anisotropy tensor T is described in the next section.

Structural Optimization Tool for the Turbulent Anisotropy
The following section presents the method proposed in this work for optimizing a mathematical expression for the turbulent anisotropy tensor T. First, the set of independent bases and invariants for the turbulent anisotropy originally derived by Pope [32] is introduced. Then, the mechanics of the optimization process is described. In this later section, we introduce the procedure for the optimization of the free-fitting constants in the turbulence models and the genetic modification of the turbulent anisotropy terms.

Independent Bases and Invariants for the Description of the Turbulent Anisotropy
A functional representation using tensor sets can be proposed for the turbulent anisotropy tensor as follows [32]: where T ik (S, Ω) are sets of tensors that all together form a complete base, where i is the order of the tensor set, k is the index of the tensor in the set of order i and c ik (I S , I Ω ) the projection functions. These projection functions are sought to be a function of the flow invariants I obtained from S and Ω to provide scaling, translation, and rotational invariance to the fitted model. To reduce the number of tensor sets and invariants, we can use the Cayley-Hamilton theorem [33], which states that a tensor must satisfy its own characteristic equation: where I a , I I a , and I I I a are the coefficients of the characteristic equations of T ik and c ik = c ik (I S , I Ω ) (see [32] for the details of these coefficients). Developing this equation, we can find a set of constraints for the turbulent anisotropy tensor T [32], which yields the unique set of independent tensor sets and invariants presented in Table 2. In this table, the characteristic time scales of the eddies τ S have been defined as τ S = k/ . Note that the tensor sets and invariants presented in Table 2 form a complete, closed base, which means that any physically realizable turbulence anisotropy tensor can be projected on these tensor sets and invariants without any residual error. Once the independent tensor sets and invariant sets have been defined, we use a method inspired by prioritized grammar expression programming [24] to seek an expression for the turbulent anisotropy tensor. Table 2. Tensor base defined by the tensor T p and invariants I q for the turbulent anisotropy tensor. For simplicity, we dropped the two index notations and define them in the table: 43 and T 10 = T 51 .

Order (Base i) Member T k and the Associated Invariants I
First Order (i = 1)

Procedure for the Turbulent Anisotropy Optimization
The procedure for finding an optimized expression for the turbulent anisotropy tensor is presented in Figure 2. The goal is to use the set of tensors T p and the fitting functions c p (I 1 , . . . , I 5 , x) to yield an optimal turbulent anisotropy tensor. For consciences, the set of fitting functions and tensors describing each turbulence model T (q,j) are encapsulated in a vector product as follows: where q is the numbering of the turbulence model within a population of candidate models, j refers to the iteration of the optimization process for this model, and T (q,j) r refers to a particular tensor set formed by a functional combination of the original tensor base T p outlined in Table 2, i.e., T (q,j) r = f (T 1 , . . . , T 10 ).

Figure 2.
Schematic diagram of the optimization process to find an expression for the turbulent anisotropy tensor. This process is called structurally optimization of the turbulent anisotropy tensor.
As depicted in Figure 2, the optimization procedure goes as follows. First, an initial random population of N models T (q,1) for the turbulence anisotropy tensor is proposed (Step A). Using a standard turbulence model in the literature, e.g., the k − turbulence model, and the description of the turbulence anisotropy tensor provided by each candidate, we compute the flow field predicted by each candidate. The predicted flow field for each candidate is then compared against a reference flow field (e.g., from an experimental database) and the constants in each c (q,1) r function are fitted to minimize the difference between the predicted and reference flow field (Step B). This leads to a series of fitted candidates, which are ordered in an array from lower to higher fitting error. The best fitted model is then compared against an error bound criterion. If the error of the best fitted candidate is smaller than this criterion, we select this candidate as the turbulence model (Step C). Otherwise, we move into the genetic optimization process. During the genetic optimization process, a series of deterministic, prioritized genetic operations are applied to ordered candidates in the array. This leads to the creation of 2N new candidate models (Step D). These models are then ordered in an array with increasing error. The error is obtained by comparing the models against the reference flow. Following, the N best models of this 2N-models population are selected (Step E). In the next iteration (j = j + 1), their constant coefficients are now fitted with the reference flow solution (Step B), and a new series of genetic operations is applied if the error criterion is not satisfied (Step D). The loop continues until one model satisfies the error-bound criterion. Details about each step in this optimization process are provided in the following paragraphs.

Model Initialization
An initial population on N model candidates is proposed at the beginning of the optimization cycle. The number of candidates to be used in this population can be defined by the user. In general, we observed that better results are obtained for large values of N, i.e., N ≈ 10. This is because this allows us to obtain higher order nonlinear tensor models by increasing the number of possible combinations in the genetic optimization process, which yield a more succinct description of the turbulent anisotropy with a better extrapolation ability. Therefore, a set of initial models (i.e., combinations of tensors of the base) are constructed for describing the anisotropy in the flow, for example, at iteration 1, the first model proposed could be s 1 I 1 T 1 , where s 1 is a fitting constant, I 1 = I 1 (x) is the first invariant of the turbulence anisotropy, and c 1 = s 1 I 1 (x); we will use the symbol s for fitting constants in this discussion. The set of candidate models at one iteration is named the test population. The fitting functions in c (q,j) can have a spatial dependency via an independent invariant set (I n = [I 1 , I 2 , . . . , I 5 ]), which depends on space (i.e., the position vector x). This provides the fitting functions c (q,j) with spatial description power, while avoiding overfitting. In summary, for each model (q) at iteration (j), selecting for example the k − for the RANS turbulent viscosity, we find a Reynolds Stress Tensor of the form: As previously discussed, the optimization process then starts with the random generation of a set of N turbulent model candidates, which entails defining a set of candidates from the bases { T p , invariants {I i }, and constants {s m } sets. Each candidate model T (q,j) is constructed by randomly combining these sets by a group of allowed mathematical operations defined in the expressions operator set: In the present case, the operator set ϑ contains scalar-matrix multiplication ( * R 2 ,R 0 ), scalar multiplication ( * R 0 ,R 0 ) and division (/ R 0 ,R 0 ), matrix addition (+ R 2 ,R 2 ) and subtraction (− R 2 ,R 2 ), and scalars addition (+ R 0 ,R 0 ) and subtraction (− R 0 ,R 0 ). This is the only random step in the process, and it can be turned into computationally deterministic by fixing the seed in the random number generator.

Definition of the Error Metric
Each of the generated turbulent candidates is then implemented in the RANS averaged Navier-Stokes Equations (3), (4) and (6), and a simulation is performed to yield an averaged velocity field u (q,j) (x, t), for model candidate q at iteration optimization iteration j. The predicted velocity field is then used for the definition of an error metric against available flow data for optimizing the candidates. This error metric is defined as follows: where L (q,j) w is the error function for the current model candidate, u data (x, t) is the velocity of the dataset used for calibrating the turbulent model, · is some scalar metric on the difference between the velocity fields, and f is a function of this metric. Moreover, the second term on the right-hand side is a Levesque regularization correction [34] that penalizes turbulent models as they increase their length (l), where the length l is the total number of monomials in model c (q,j) T (q,j) , i.e., number of T p composing each T (q,j) tensor function plus the number coefficients and invariants used in the c (q,j) fitting functions. For example, the model s 1 I 1 T 1 will have a length l = 3. This penalty has been introduced because the extrapolation behavior of the turbulence model has been observed to deteriorate with the number of T p composing each T (q,j) tensor function. The error metric function f and the penalty parameters µ and l MAX are user defined parameters. In this work, we perform grid-search studies to find the optimal values of these parameters using the envelopes µ ∈ [0.001, 1.0] and l MAX ∈ [1200]. Note that we focus the optimization metric on the velocity since this is the field that introduces advective transport in the temperature and precursors concentration fields in liquid fuel reactors. Nonetheless, the method accepts different kinds of metric functions. Moreover, in practice, one should define the error function according to the desired output from the turbulence model. For example, if a turbulence model is used to calculate the pressure drop in a section, then a suitable error function would weight predominantly the velocities next to walls, whereas if a turbulent model is used to describe the bulk mixing effect, then the error function weight should be given to the bulk values of the velocity. In a similar form, the penalty parameters can also be defined according to the need. Moreover, if a very precise model is necessary for a specific flow condition, then µ → 0 , whereas if a robust extrapolable model wants to be developed over a large dataset the value of µ will be increased.

Optimization of the Constants in Each Candidate Model
Once the error function and penalty parameters have been defined, then they can be used to optimize the set of constant coefficients ({s m }) for each of the fitting functions c (q,j) of each candidate model T (q,j) . In mathematical terms, this is an unconstrained optimization problem for a fitness function L Hessian for an optimization step n, the iterative quasi-newton optimization problem for all constant coefficients in a candidate is defined as follows: Note that the constants included in the c (q,j) m can define a non-convex optimization problem or, equivalently, the Hessian matrix is not necessarily positively defined. Therefore, the Hessian is factorized by a Cholesky algorithm with row pivoting before the solution. During factorization, the Hessian is decomposed as H = LDL T , replacing the negative values in the matrix diagonal with small positive values during the factorization process. This process allows us to find a descendent direction towards a minimum independently of the convexity of the problem. This technique is also known as the active set method in the literature [35].
Once the constant coefficients have been fitted, the calculated error of each turbulent model is checked against an error acceptance criterion. If none of the models is satisfactory (i.e., the error for all the models is too large), a genetic optimization process is performed using all candidate models.

Genetic Optimization of Candidate Models
To reduce the number of optimization cycles, we apply a deterministic tensorial optimization method inspired by the symbolic functional optimization procedure developed by Worm and Chiu [36], which we will name Prioritized Grammar Tensorial Regression (PGTR) in this paper. This optimization method is described in the following paragraphs.
The key idea of the PGTR is to replace the randomness in the selection, survival, and replication operations of genetic algorithms with an optimized set of deterministic operations. The first step for the modification of the turbulent models consists of building a priority stack with the candidate models whose constants were fitted in the previous step. We order the stack from largest (top) to smallest (bottom) error regarding error metric L (q,j) w . The goal of this stack is to prioritize combinations with the best-fitted models to progressively obtain better-fitted models during iterations. Once ordered in the stack, standard genetic operations are applied to the mathematical expression of the models in the stack to generate new model candidates. The operations used in this work are onepoint mutation, term mutation, point transposition, term transposition and recombination. Examples of these operations are illustrated in Figure 3. However, contrary to standard genetic algorithms, these operations are applied in a deterministic way in the PGTR. The key to the method is expressing the genetic operations in a way that can be optimized by standard Bayesian optimization algorithms. To generate the ordering of genetic operations in the stack optimizable, we label each operation with a TOTC (Term Operation Type Complement) label. The first letter T ∈ N identifies the term where the modification operation will be applied. This is a signalization to the term in the expression representing the turbulent model that will be subjected to modification. The second letter O ∈ N refers to the number of the type of operation that will be applied. The third letter T ∈ N refers to the type of term in which the operation is applied. This term can be either a constant, a tensor, or an operator (sum, subtraction, multiplication, or division). Finally, the letter C ∈ N refers to the flag indicating the complement of the operation. This flag will depend on the type of operation to be applied. For instance, in a transposition operation, the complement flag will indicate the length of the term being transposed and the new position where the transposed term should be inserted. In mutation operations, the complement flag is composed of a signaling number indicating the position of mutation, the length of the term mutated and a set of numbers indicating the mutations that should be introduced in the bases or in the constants. Finally, the complement flag in the recombination operation is composed of two numbers that signalize the terms combined and a third signalization number that indicates the position in the first term where the second term needs to be introduced.
All models are represented through Open Reading Frames (ORFs) for simplicity when performing Bayesian optimization. We then use the Bayesian optimization implemented in the scikit-optimize Python library [37] to find optimal values for the TOTC flags during the optimization process. During the Bayesian optimization process, for each new model created in each Bayesian iteration, we fit the constants and insert it in the priority stack according to its fitting error. Then, a new Bayesian optimization step is performed. This method allows us to prioritize genetic operations in the priority stack that yield high accuracy improvements in the turbulent models. We observed that starting with large exploration rates in the Bayesian method always guaranteed the convergence of our models to a global minimum.

Final Remarks
Normally, we could continue this optimization process until finding a model with the required accuracy according to L (q,j) w . However, in our experience, this generally yields a lengthy process since the error metrics of the best-fitted models improve slowly with subsequent iterations after a few Bayesian steps have been performed. Therefore, once N new models have been generated, yielding a prioritized stack of 2N models, we select the N best-fitted models as new candidates and perform a new global optimization step as shown in Figure 2.
The developed optimization procedure was implemented in the OpenFOAM ® openssource software v5.0. It is interesting to note that in all the test cases performed, the maximum number of global iterations required to find a turbulent model that satisfied the error criterion was 8 when using n = 10 models as candidates. This means that, at most, 8 × 10 = 80 RANS steady-state simulations were performed to obtain an accurate turbulent model. In most test cases, the global optimization process converged after three iterations.
In this section, we have presented a methodology for obtaining data-driven turbulent models to minimize an error metric L w . This methodology is implemented for two test problems in the next section.

Results
In the following sub-sections, the previously described methodology to optimize turbulence models is implemented for two cases. First, the case presented in Section 3.1, a Backward-Facing Step (BFS) with a small aspect ratio, in which shear-driven boundary layer detachment and reattachment dominates over the vorticity-induced mixing after the tripping of the boundary layer. The effect of turbulence anisotropy in this type of flow is close to what can be found in some simple components of a liquid fuel reactor (such as a pipe with a section change). This first part of the study seeks to answer the two following questions: (i) Can an optimized RANS model yield high accuracy for a simple boundary layer detachment problem? and (ii) How important is the structural optimization of the turbulent anisotropy in yielding accurate results? A second case is presented in Section 3.2 and consists of a large aspect ratio Backward-Facing Step (BFS). This case is more representative of the recirculating flow field that will exist at the entrance of the core cavity in an MSFR. We then evaluate if we can develop and optimize RANS to yield high accuracy in a complex flow field after a boundary layer detachment into an opened cavity.

Optimizing Turbulence Model for Backward-Facing Step (BFS) with a Low Aspect Ratio
The flow geometry for this problem consists of the small aspect ratio Backward-Facing Step (BFS), first proposed by Kopera et al. [38], that is presented in Figure 4. After this sudden expansion, the boundary layer generated in the inlet throat detaches, producing a shear layer that yields high turbulence production downstream. After a short distance, this layer reattaches to the down part of the channel. In the inferior part of the channel after the expansion, flow recirculation is produced as a product of the interaction of the detached shear layer with the fluid in this region. The flow data used to optimize the RANS models with the method introduced in the previous section was generated using Direct Numerical Simulation (DNS) [39] for a Reynolds number in the bulk of the step Re ∈ {8000; 9000; 10, 000}. In order to check the accuracy of our DNS simulations, the results of the BFS obtained for Re = 9000 were compared against the DNS results published by Kopera et al. [40]. It was found a maximum difference below 1% in the root-mean-squared velocity field. The turbulence RANS models were trained according to the proposed method for Re = 9000 and then they were tested for two different Reynolds numbers. For training the turbulence models, we define the merit function as the quadratic weighted average of the difference between the stream-wise predicted velocity (u MODELil ) and the time-averaged DNS reference velocity (u DNSil ), for a set equidistant points i = 1, . . . , 100 over the four vertical lines l = 1, . . . , 4 (shown in the zoom of the BFS section displayed on Figure 4) as follows: where y il is the vertical position i in line l and H l is the total height of this line. For the RANS and LES simulations, we have assumed the BFS geometry to be 2dimensional (assuming plane symmetry) and the domain to consist of a cut across the x − y plane. For the discretization process, we have adopted a regular structured quadrilateral mesh, which has a maximum aspect ratio of 5 and maximum skewness ratio of 10 −6 . We refined the mesh until the accuracy of the predictions in the bulk velocities measured by the error L w were unchanged by further refinements for each turbulence model studied in the section. Then, we use in the simulations the finer grid of all RANS models for performing all RANS simulations. As expected, the mesh in the region next to the sudden expansion was key for the accurate description of the turbulence phenomena in the bulk region of the BFS. In addition, the predictions in the tripping of the boundary layer were very sensitive to the wall functions adopted in the simulation. Therefore, since our goal is to make the results only dependent of the turbulence model used, we implemented an enhanced wall treatment procedure [40]. In this regard, we adjusted the centers of the mesh cells along the walls of the BFS, using a uniform inflation ratio boundary layer with between five layers, to obtain a dimensionless wall distance (y + ) varying between 0.5 and 1.5 for these cells. The meshes were generated using the snappyHexMesh utility in OpenFOAM. The RANS mesh had approximately 1.3 × 10 5 computational cells, whereas LES one had approximately 1.7 × 10 6 .
The following turbulence models were used as baseline models to predict the averaged flow velocity in the Backward-Facing Step (BFS):

•
RANS k − standard model [27] • RANS Wilcox k − ω model [28] • RANS Reynolds Shear-Stress Transport (RSS) model [29] • LES with Multiple-Gradient model [30] These models were evaluated according to the cited references. The first three RANS models were already implemented in OpenFOAM ® v5.0. We implemented the LES model using OpenFOAM libraries.
We now proceed to the analyses of the results in light of the questions discussed at the beginning of this section. As a baseline nonlinear turbulence model, we used the turbulent anisotropy tensors of the base presented in Table 2 up to the third order (i.e., from T 1 to T 6 ) and we fit the projection constants for the Re = 9000 data using the metric function defined in Equation (13). The parameters for the genetic optimization process were N = 10, µ = 0.2, and l MAX = 10. Note that this model, called non-linear cubic model (NLCM), is a reference of the performance that can be expected from a nonlinear turbulence model without structural optimization. This means that this is the best nonlinear cubic model that can be obtained without optimizing the tensor bases describing the turbulence anisotropy tensor. The fitted model reads as follows: The results of all models were then compared against the DNS simulation for Re = 9000 (named DNS data) in Figure 5. The general performance of the k − ε model for the BFS is poor. This is because of the overestimation of ε and k after the tripping of the boundary layer (too much diffusion). Thereupon, the k − ε model predicts an early development of the velocity profile past the BFS expansion. The results predicted by the k − ω model are, on the contrary, in much better agreement with the DNS data. This is because when calculating ω = ε/k both over-predictions of ε and k compensate remarkably well. The non-linear viscosity model has a good agreement with the DNS data close to the middle stream, but it performs badly next to the walls. This is because the closure coefficients of this model have been fitted using a merit function Equation (13) in which points next to the walls have low weight. Hence, the flow behavior near the wall is not solved with high accuracy. The RSS transport model shows good results in the region immediately after the detachment of the boundary layer. However, it underestimates the dissipation of the turbulent kinetic energy in the streamwise direction, causing an over diffusion of the shear layer and significant errors in the predicted velocities downstream. This is mainly because the model assumed for the turbulent dissipation in the RSS simulations does not perform accurately in regions with a large amount of swirl. Finally, the multi-gradient LES model shows an almost perfect agreement with the DNS data. We then apply the PGTR method to fit a turbulent model in the BFS geometry that is able to perform with the desired accuracy at a reduced computational cost. The parameters for the genetic optimization process were N = 10, µ = 0.2, and l MAX = 10. We set the target accuracy for the model at 2.5% after experimentation. We take the k − model as a reference for the turbulent viscosity field ν t and optimize the anisotropy tensor by taking the error function L w defined by Equation (13) as a base function with penalty parameters µ = 0.05 and l MAX = 50. The resulting merit function reads as follows: The penalty parameters µ and l MAX were defined via a grid-search method on the testing set. The final anisotropy model constructed via the optimization process for this case is the following: Note that the coefficients of the optimized model are close to the ones of optimized for the cubic model (NLCM) presented in Equation (14) with the exception of the new factor (0.98 + I 2 /I 1 ) ≈ (1 + I 2 /I 1 ) that pre-multiply the anisotropy second order terms. The term I 2 /I 1 is the norm of the rate of rotation tensor normalized with the norm of the strain rate of the velocity. Hence, this term can be thought of as a projection factor modulating the projection of the second-order terms into the anisotropy according to the rotation rate. This term is close to one in most of the bulk of the flow but larger than one in the flow regions close to the walls. This dumps the effect of the second-order description of the anisotropy near the walls, which produces inaccurate results in the NLCM. Furthermore, note that third-order tensors are not present in the optimized model. This is logical since they do not contribute much to the description of the flow field (3D effects in this flow are neglected).
The optimized model in Equation (16), referred to as the non-linear k − model, is presented in Figure 6 for the streamwise velocity over the measurement lines along with the standard k − model and the DNS data. A very good agreement is observed between the optimized turbulent model and the DNS simulations, which is significantly better than the NLCM model without structural optimization fitted in Equation (14). However, a small over-prediction in the streamwise velocity is found due to the action of the damping term in the bulk of the fluid that slightly decreases the values of the quadratic terms in this region. Nevertheless, the global error for the optimized model L w = 2.47% is much lower than the other models (except for the LES) and the computational time is 1187 s, which is reasonable and satisfies the required error criteria. This demonstrates that a high-accuracy RANS model can be developed for the Backward-Facing Step (BFS) with a low aspect ratio.
The comparison of the average merit function defined in Equation (13) for the training and validation of Reynolds numbers is presented in Table 3. It can be noted that the performance of the nonlinear k − model fitted without structural optimization (NLCM), Equation (14), and the one with structural optimization (from GETFOAM), Equation (16), outperform the standard k − model. Considering that the former two models use standard k − as the underlying model for computing the turbulent viscosity, it is interesting that the fitting of the anisotropy tensor has a considerable impact on the performance of the models. Moreover, the structurally optimized turbulence models outperform more sophisticated RANS models such as the k − ω and the RSS transport ones. Regarding extrapolation of the fitted models to Re = 8000 and Re = 10,000 , we observe that the extrapolation performance of the model with structural optimization is significantly better than the one without structural optimization. This shows that the larger order terms included in the non-structurally optimized model cause overfitting issues. This observation also shows the need of including structural optimization in higher-order structural models. Finally, we observe that the structurally optimized model yields the best performance among the training and validation sets for all RANS models. Nonetheless, its performance is still worse than what can be expected with the LES model. However, the computational runtime of the LES model is approximately 77 times larger than the structurally optimized model. Therefore, the structurally optimized model can be a good fit for applications in which the computational cost of an LES model is prohibitive and only the averaged flow values are needed. In the next section, these observations are re-analyzed for a backward-facing step with a large aspect ratio.

Optimizing Turbulence Model for Backward-Facing Step of Large Aspect Ratio
In the previous section, we studied the advantages of structural optimization for a small aspect ratio Backward-Facing Step (BFS). The flow field in this problem is regulated by a tripping of the boundary layer in the sudden expansion and a shear-induced reattachment of the boundary layer relatively close to the detachment point. This type of BFS can be used to improve turbulence models for some of the components of liquid fuel reactors but not for the core cavity existing in the MSFR for example. Indeed, in the MSFR core cavity, once the boundary layer detaches at the entry of the cavity, the production of coherent turbulence in the shear layer, the vortex merging, and the anisotropic turbulence stresses govern the flow field before the re-attachment further downstream. Therefore, as our second test case, we have decided to study the Backward-Facing Step (BFS) section with a large aspect ratio.
The design of the large aspect ratio backward-facing step is presented in Figure 7a. It consists of a backward-facing step with a 4:1 aspect ratio. The width of the section was designed to maintain a quasi-2D flow profile while minimally affecting the turbulent dissipation perpendicular to the 2D plane. Given the importance of this reactor design, an experimental section was built in this case to test the flow in this BFS geometry. The experimental setup is shown in Figure 7b. The flow in the section was established by setting a pressure difference between upstream and downstream tanks rather than using a pump. This ensures that the pressure-driven flow perturbations normally generated by mechanical pumps are absent in the flow field. A proportional-integral control system was designed in the LabVIEW software [41] to fix the flow velocity which was measured by an ultrasonic flowmeter right after the exit of the upstream tank. The action of the control system was an automatic pressure release valve in the downstream tank. An outward pipe of a length of 2.0 m was used for the inlet to the test section to ensure that the flow field was fully developed before entering the Backward-Facing Step (BFS). The section was built in transparent, low-rugosity plexiglass to allow for optical measurements. Thus, a developed Poiseuille flow boundary condition can be assumed in the inlet of our simulations, no rugosity zero velocity conditions have been applied to the enhanced wall treatment, and standard a pressure-outlet condition has been assumed for the section's outlet. A Particle Image Velocimetry (PIV) technique was used in the BFS test section to measure the flow field. In our PIV setup, the laser light sheet entered the channel by its top part and the high-speed, high-resolution camera was placed in front of the test section. A squared grid-pattern metal sheet with a spacing of 1cm was placed behind the test section for calibration and correction of the light diffraction errors during the data postprocessing. The measurement procedure was carried out using an overlapping sliding window technique. Each of the overlapping measurement windows used is depicted in Figure 7c. Time measurements were set between 2-5 min per window, which allowed us to obtain steady time-averaged profiles when postprocessing the measurements. An example of the postprocessed PIV measurements for Re = 1111 is presented in Figure 7d. All Reynolds numbers reported here are computed in the inlet section, before the sudden expansion. Implementation errors in PIV measurements can be estimated by comparing measurements in the overlapping region between windows. An average mean-squared error of~0.5% was obtained with maximum errors of~2% for the time-averaged velocity fields.
We measured three different Reynolds numbers in the test section, Re = 222, Re = 1111 and Re = 3889. The flow field produced for Re = 222 after the sudden expansion presented turbulent characteristics. We use the measurements obtained at these Reynolds numbers to train the structurally optimized RANS model with GETFOAM. Furthermore, we performed nine DNS simulations for building the test set for Re = i * 1000, ∀i ∈ [1, 2, . . . , 9]. Moreover, we performed three extra DNS simulations at the experimentally measured Reynolds numbers to validate these simulations against the experiment. We obtained a mean squared average difference of 1.2% in the 2D measuring plane with a maximum error of 1.8%, 3.4%, and 2.5% for the Reynolds numbers 222, 1111, and 3889, respectively. This is very consistent with the previously estimated experimental errors.
As in the previous case, for the training set, we have extracted the streamwise velocity over vertical lines in the mid-plane of the test section, using 100 equidistant points over each line, to build the error function. The seven vertical lines used for this case are presented in Figure 7c. These flow lines capture the profile just before the tripping of the inlet boundary layer, the turbulent recirculation phenomena after detachment, and the development of the flow post-attachment of the boundary layer. The error function for this case is then defined as follows: where the reference averaged velocity can either be the experimentally measured or DNScomputed one depending on whether the training or testing values are reported, respectively. In the analyses, we have performed 3D simulations using RANS and LES models (contrary to the previous section where 2D cases were studied). We implemented an equivalent mesh convergence procedure for these simulations than the one in the previous case and the L w defined in Equation (17) as an error function. The procedure resulted in 1.5 × 10 5 computational cells for the RANS simulations and 2.7 × 10 6 for the LES one. Next, we proceeded to build a structurally optimized turbulence model. The parameters for the genetic optimization process were N = 10, µ = 0.1, and l MAX = 10. Like in the study reported in the previous section, we use the standard k − model to compute the turbulent viscosity field. The metric function used for the optimization of this model is the penalized average L w error over the three training Reynolds number, which is defined as follows: As in the previous case, the penalty parameters were obtained by a grid search algorithm over the testing set. The structurally optimized turbulence model obtained reads as follows: Contrary to the previous case, higher-order anisotropy terms are selected in the closure model. We can further analyze the fitted model to find physical meaning in its terms. We note that the third-order term 0.511 I 1 I 3 s Ω·S 2 − S 2 ·Ω , captures the shear strain rate unbalance produced by fluid rotation Ω·S 2 − S 2 ·Ω, modulated with the amount of shear strain rate over the amount of rotation I 1 I 3 s Ω 2 ·S 2 + S 2 ·Ω 2 − 2 3 S 2 : Ω 2 I − (Ω : Ω)· S 2 − 1 3 (S : S)I corrects the total shear strains rates produced by flow rotation Ω 2 ·S 2 + S 2 ·Ω 2 with the total shear and rotation rates 2 3 S 2 : Ω 2 I and the non-isotropic compounded rotation-shear rate (Ω : Ω)· S 2 − 1 3 (S : S)I . Since the aspect ratio of the backward-facing step is now larger and the recirculation after the boundary layer tripping now plays a key role in determining the flow field it is logical that the basis describing the interaction between rotation and shear rate, T 6 and T 7 , appear in the optimized turbulence model.
As an example of the key role of these terms in describing the correct behavior of the flow field in the recirculation region, the contour plots of the streamwise velocity fields predicted by a laminar solver, a standard and an optimized RANS model, the LES model, and the DNS model for Re = 1111 are presented in Figure 8. When comparing with the DNS model, it is observed that standard RANS models overpredict the development of the recirculation region. Further analysis showed us that this is because the added turbulent diffusivity in the shear strain layer produced after the boundary layer detachment was overpredicted. This can be explained by the fact that the RANS models were not capturing the modulation of shear strain produced by the rotation rates expected in the recirculation region. One can see that, when compared with the DNS solution, standard RANS models over-predict the extent of the recirculation. However, the structurally optimized model correctly captures the sharp descent in the shear stress layer towards reattachment.
Contrary to the previous example, instead of analyzing the model performance per sampled line, average errors are presented when analyzing the performance of the structurally optimized models. This is because the detailed analysis of the predicted solutions over each line could lead to erroneous interpretations as there are many physical effects playing a role over each line. The average errors over the training and validation set are presented in Table 4 along with the computational runtime for each of the models. It is observed that the structurally optimized turbulence model presents an error that is significantly below the other standard RANS models over the training set. As expected, the performance of the optimized model deteriorates over the testing set. Nonetheless, the observed deterioration in performance is not large and the optimized turbulence model still significantly outperforms the standard RANS models over the testing set. Like the previous case, the LES model presents a smaller error than the optimized RANS model. Nonetheless, the computational time for the LES model, in this case, is about 382 times larger than the optimized RANS model. The added computational cost in the LES simulation is due to a finer mesh and larger simulation time required in the transient 3D simulation to yield steady time-average results. As in the previous case, the structurally optimized RANS model is a good compromise between accuracy and computational speed.  To understand more in detail the applicability of the RANS model beyond the fitting domain, the error measured by L w defined in Equation (17) for each Reynolds number in the testing domain is presented in Figure 9. We observe that the error is smaller for Reynolds numbers in the [1000-5000] domain. This is expectable since the optimized turbulence model has been developed for Re = {222, 1111, 3889}. Then, the errors suddenly increase for Re = 6000. This is also expectable since as the Reynolds number increases, the ratio between the strain rate tensor norm and the vorticity rate tensor norm also changes in the recirculation. DNS studies showed that for Re ≥ 6000 a counterclockwise recirculation bubble appears in the top part of the BFS. This bubble can already be seen for the k − ω and RSS transport model in Figure 9. As the models have not been trained to reproduce the recirculation phenomena in this bubble and the associated flow modification in the boundary layer reattachment region, it is logical that the optimized model cannot capture the flow field phenomena with a similar accuracy to that in the fitting range. Nonetheless, as we keep increasing in Reynolds number, we observe that the error function starts to stagnate. This is mainly because the turbulent dissipation in the shear layer starts playing a key role in turbulence dissipating as we increase the Reynolds number and, hence, capturing the accurate ratio of shear strain rate over rotation rate becomes less important on the model performance. In brief, the error yield by the optimized RANS turbulence model is acceptable across the whole range of Reynolds numbers of interest.

Assessment of the Impact of the Turbulence Model in the Molten Salt Fast Reactor
For justifying the importance of improving the accuracy of the RANS turbulence models, we assess the impact that the new improved turbulence models can have in liquid fuel reactors. For this purpose, we use as an example the Molten Salt Fast Reactor (MSFR) and we simulate the steady state operating conditions of this reactor [42]. Note that although the optimized turbulence model should perform better than standard turbulence models when capturing the recirculation in the MSFR cavity, we are by no means claiming that the models developed from the Backward-Facing Step (BFS) in Section 3.2 would provide better results in the MSFR. This is because of the large difference in Reynolds numbers and geometries between the BFS and the MSFR cavity. In this exercise, we just intend to quantify the impact that an optimized turbulence model has on the coupled velocity, pressure, neutronics, neutron precursors, and energy fields in the nuclear liquid fuel to justify future developments.
To solve the coupled multiphysics problem in the liquid fuel in this reactor, one must solve the equations of conservation of mass and linear momentum, Equations (1) and (2), coupled with an equation for the balance of energy, neutrons, and neutron precursors. For the neutron balance equation, we have used a six-group diffusion model, while a linear-log interpolation is applied in the cross sections as a function of temperature for modeling the thermal expansion or contraction of the fuel and the Doppler effects. For the neutron precursor concentration equations, we used an eight (8) families model while we considered precursor advection in the flow. The complete multiphysics model is not discussed in this paper for conciseness but we refer the reader to our previous work in [43,44] for more details. Of particular interest for this work is that the optimized turbulence model will predict a different turbulence stress tensor than the standard model. The modified viscosity field will affect the velocity field and the associated energy and neutron precursor's turbulent diffusivities. Thus, temperature and delayed neutron precursor concentration fields will be directly affected, while neutronics will be indirectly affected via the modifications in the delayed neutron precursor source and the density-Doppler effects via the temperature field.
The simulation domain and computational mesh used in the simulations are depicted in Figure 10. We impose standard wall boundary conditions everywhere in the reactor for the mass, linear momentum, energy, and delayed neutron precursors conservation equations. We implemented a diffusive albedo boundary condition for neutronics with an albedo coefficient of 0.3. We model the pump and the heat exchanger as uniform porous media, with a uniform volumetric heat extraction rate and momentum source for the pump and the heat exchanger, respectively. We calibrated the heat extraction rate to maintain the design mean operation temperature of the reactor (~1000 K), while the momentum source was set to yield the required loop time (~4 s). We set the thermal power of the reactor to 3 GW. We call this state the nominal operating conditions. The two optimized turbulence models developed in the previous section were tested. These are the one developed for the Backward-Facing Step (BFS) with a small aspect ratio in Equation (16) and the one developed for the Backward-Facing Step (BFS) with a large aspect ratio in Equation (19). However, instead of using the standard k − model for computing the turbulent viscosity field, we use the realizable k − model proposed by Shih et al. [45]. The reason is that this latter model performs better for flow fields regulated by boundary layer detachment under adverse pressure gradients and where buoyancy forces exist (i.e., natural convection contribution), such as the ones that may exist in the MSFR during some operating conditions. Since the flow phenomena in the large aspect ratio Backward-Facing Step (BFS) should be closer to the one observed in the large expansion of the flow into the MSFR core cavity, we analyze the results of this model in more detail. Nonetheless, note that due to the larger flow velocity in the MSFR, the actual balance between shear strain and rotation rates in the heat exchanger region follows more closely the ones obtained in the low aspect ratio Backward-Facing Step (BFS).
The contour plots for axial-cut thermal-hydraulics flow fields obtained in the simulation of the MSFR are depicted in Figure 11. The temperature rises as the flow circulates from the bottom to the top of the reactor cavity. Nonetheless, the temperature field in the rapid expansion into the entry of the cavity is approximately uniform. Hence, buoyancy should have little effect on the turbulent flow field at this expansion and the flow conditions approximate the ones obtained in the high aspect ratio backward-facing step. By analyzing the velocity field, one observes a boundary layer detachment in the core cavity with a recirculation towards the cavity side walls. The predicted relatively rapid re-attachment of the boundary layer to the sidewall after this recirculation is due to the ability of the optimized turbulence model to capture shear strain dissipation due to rotation rate. Finally, we observe an excess production of turbulence after the flow injection in the core cavity. The turbulent kinetic energy then extends into the center of the core cavity and causes a large turbulence viscosity in the bottom part of the reactor. This behavior is produced by capturing the imbalance between shear strain and rotation rates in the optimized turbulence model, which does not overpredict turbulence diffusivity and, thus, extends the effects of turbulence beyond the entrance to the core cavity. The domain-averaged mean-squared deviations in the molten fuel salt velocity, pressure, and temperature fields and the effective neutronic parameters obtained using the two structurally optimized turbulence models from the large and small aspect ratio BFS are compared in Table 5. The deviations of both models are calculated with respect to the standard realizable k − model. It is observed that significant changes are produced in the thermal-hydraulics field. In fact, the maximum temperature in the reactor increases by 27 K and 22 K when using the turbulence models optimized in Equations (16) and (19), respectively. This temperature rise is not negligible and should be accounted for in the thermal design constraints. However, effective nuclear parameters present small changes regarding the turbulent model used. The uncertainty generated by the turbulence model on these parameters is below the one arising from the uncertainties in the nuclear crosssections. This can be partially explained by the fact that this reactor has a fast spectrum and thus the neutron mean free paths are long enough to average the flow field heterogeneities. Finally, Table 5 should not be interpreted as a well-defined measure of the uncertainty that will be introduced by turbulence modeling, but rather as an indication of the impact that turbulence models may have on the coupled multiphysics fields. Table 5. Mean-squared relative difference over the computational domain for the coupled fields in the Molten Salt Fast Reactor and total change in the effective neutronics parameters using the structurally optimized turbulence models proposed in Equations (16) and (19) vs. a standard k − realizable model [45].

Optimized Model Given by
Equation (

Conclusions
In the present work, we have proposed a method based on genetic algorithms for the deterministic development of explainable, data-driven turbulence models. This method performs the structural optimization of the Reynolds anisotropy tensor with respect to a high-fidelity dataset, while it preserves the baseline models for the turbulent viscosity. Moreover, the integration of a Bayesian optimization loop for determining the order of the genetic operations in our optimization method guides the convergence of the method towards the global minima. The convergence to a global minimum has been observed for high and low aspect-ration backward-facing step flows. However, the convergence to a global minimum in more complex flow configurations, such as the one in the MSFR, must still be demonstrated.
We have applied this method to build optimized turbulence models for Backward-Facing Steps (BFS) sections with low and high aspect ratios. We observed that the developed turbulence models could not only yield high accuracy results for the streamwise velocity fields but also presented a good extrapolation ability beyond the training set. For instance, for the low aspect ratio backward-facing step, the structurally optimized RANS model yielded average mean-squared errors below 3%, whereas its baseline k − model presented errors of~14%. Moreover, we observed the need for structural optimization of the turbulent models in this geometry, as the direct fitting of nonlinear turbulence models is prompt to overfitting. A similar conclusion was obtained for the large aspect ratio BFS section for which a RANS turbulence model could be structurally optimized to produce errors of 3.2% and 4.1% over the training and validation sets, respectively. In all cases and for the training and validation sets, the accuracy of the structurally optimized RANS models was the best among all RANS models and was comparable to the one of LES.
Finally, the Molten Salt Fast Reactor has been used as a multiphysics test case for the impact of turbulence models on the multiphysics coupled fields in liquid fuel reactors. For this purpose, we applied the two structurally optimized turbulence models in a multiphysics simulation of the MSFR to solve for the steady-state operation of the reactor. It was observed that turbulence is unlikely to play a determinant role in affecting the averaged effective neutronics parameters. However, it has been shown that it can significantly change the temperature and precursors concentrations fields in the reactor, which have considerable implications in the thermomechanical design and species tracking in this one.
Future work will aim at extending the current tool for optimizing turbulent models for energy and species diffusivities and analyzing the effects that optimized turbulence models can have on species tracking and corrosion rates. The present approach could be used to design relatively simpler experiments to improve turbulence modeling for these advanced reactors.