Data-driven turbulent heat flux modeling with inputs of multiple fidelity
Matilde Fiore Enrico Saccaggi Lilla Koloszar Yann Bartosiewicz Miguel Alfonso Mendez


APS/123-QED




matilde.fiore@vki.ac.be
Environmental and Applied Fluid Dynamics Department, von Karman Institute for Fluid Dynamics, Belgium


Department of Mechanical and Aerospace Engineering, Politecnico di Torino, Italy
Environmental and Applied Fluid Dynamics Department, von Karman Institute for Fluid Dynamics, Belgium


Environmental and Applied Fluid Dynamics Department, von Karman Institute for Fluid Dynamics, Belgium


Institute of Mechanics, Materials and Civil Engineering (IMMC), Universite catholique
de Louvain (UCLouvain), Place du Levant 2, 1348 Louvain-la-Neuve, Belgium


Environmental and Applied Fluid Dynamics Department, von Karman Institute for Fluid Dynamics, Belgium


             

§ ABSTRACT 
Data-driven RANS modeling is emerging as a promising methodology to exploit the information provided by high-fidelity data. However, its widespread application is limited by challenges in generalization and robustness to inconsistencies between input data of varying fidelity levels. This is especially true for thermal turbulent closures, which inherently depend on momentum statistics provided by low or high fidelity turbulence momentum models. This work investigates the impact of momentum modeling inconsistencies on a data-driven thermal closure trained with a dataset with multiple fidelity (DNS and RANS).
The analysis of the model inputs shows that the two fidelity levels correspond to separate regions in the input space. It is here shown that such separation can be exploited by a training with heterogeneous data, allowing the model to detect the level of fidelity in its inputs and adjust its prediction accordingly. In particular, a sensitivity analysis and verification shows that such a model can leverage the data inconsistencies to increase its robustness. Finally, the verification with a CFD simulation shows the potential of this multi-fidelity training approach for flows in which momentum statistics provided by traditional models are affected by model uncertainties.  



                              
Data-driven turbulent heat flux modeling with inputs of multiple fidelity
    Miguel A. Mendez
    September 9, 2024
=========================================================================







§ INTRODUCTION


Data-driven methods are progressively entering the field of turbulence modeling as one of the most promising avenues to overcome the modeling barrier reached by traditional approaches in the last decade. A variety of data-driven approaches for the modeling of Reynolds stresses have been recently proposed <cit.>, among which artificial neural networks <cit.> are widely adopted to handle large databases and generate complex, non-linear input/output mappings. Despite the merits of these approaches, significant limitations still restrict their use for a wide range of Computational Fluid Dynamics (CFD) simulations.

<cit.> highlights the most important of these limitations, including the lack of generality in data-driven closures, i.e., their limited applicability to flows that deviate from the conditions covered by the training. The derived closures usually target specific classes of flows <cit.> and existence of general or universal data-driven turbulence models is currently being questioned by the community <cit.> due to the complex dependence of coherent structures on geometry, Reynolds number, and many other local/global conditions. Moreover, the literature points out limitations when data-driven closures are trained offline <cit.>, i.e., with high-fidelity data only in a frozen mode. This training strategy neglects the influence of the numerical methods and additional sub-models involved in the CFD setup. The frozen training mode causes model-data inconsistencies <cit.>, where the quantities taken as input by the model differ during training and testing phases. This limits the accuracy when the model is applied a posteriori. A typical example of the consistency problem is the turbulent time scale k/ϵ <cit.>, which is systematically different in high-fidelity and modeled (RANS) data, leading to inconsistencies in the data-driven closure if this quantity is one of the model inputs.

Current research proposes mitigating this issue by partially replacing the critical high-fidelity turbulent statistics with their RANS counterparts. This approach is followed, for example, by <cit.>, who developed a k-corrective frozen training strategy where the turbulent dissipation rate ε is computed by solving its RANS transport equation, taking all the other statistics from the high-fidelity database and adding a corrective term accounting for errors in the k production term. <cit.> proposed a similar approach in which ω=ε/k is computed from RANS equations evaluated with high-fidelity mean flow and Reynolds stresses. <cit.> included the inconsistency between RANS and DNS estimates of k in the outputs of the data-driven model, and used this to inform the data-driven model when it operates in conjuction with a RANS solver.   

In general, addressing model inconsistencies in turbulence modeling requires balancing the physical consistency of the data-driven model against its applicability in standard RANS solvers. On one hand, using modeled statistics enhances the applicability of the closures but inevitably introduces structural errors in the mapping to compensate for input inaccuracies. On the other hand, utilizing high-fidelity data promotes physical consistency and leverages the real correlation among the statistics, but this comes at the expense of applicability when implemented in RANS solvers.

Moving to non-isothermal problems, the modeling of turbulent heat fluxes depends on the momentum turbulence modeling, and the Reynolds stress tensor is inevitably one of the critical inputs for turbulent heat flux models. Consequently, a thermal turbulence model trained with Reynolds stresses derived from high-fidelity data makes significant errors when deployed in a RANS solver, which uses low-fidelity modeling of Reynolds stresses along with the Boussinesq approximation. 




This problem was highlighted in <cit.>, focusing on the data-driven modeling of the turbulent heat flux at near unity and low Prandtl numbers. Model inconsistencies in the Reynolds stresses were identified to be the major challenge in data-driven thermal closures, partially addressable only by the use of second order momentum turbulence closures. Modeling the Reynolds stresses with the eddy viscosity concept poses significant limitations to thermal turbulence modeling, because the turbulence anisotropy is essential to represent the heat flux vector, and is indeed the foundation of most algebraic thermal closures <cit.>. Such modelling choices are typical of both data-driven and theory-driven RANS modeling. In fact, traditional thermal models are often developed to be combined with specific momentum closures <cit.>, or variants are proposed to account for different momentum models <cit.>. 

This work explores the hypothesis that machine learning tools could identify the source of the input (e.g. high fidelity vs low fidelity momentum modeling) and adapt accordingly, or find compromises by mitigating its sensitivities to critical inputs. Both approaches could ensure an optimal compromise between physical consistency of the model when this is deployed on high-fidelity data, and robustness when this is deployed with traditional RANS solvers. 

To explore this hypothesis, this paper proposes a training method based on inputs of multiple fidelity multi-objective optimization.
More specifically, the training of a data-driven model for the turbulent heat flux, introduced in section <ref>, was carried out with a hybrid dataset consisting of high-fidelity (DNS) and low fidelity (RANS) data. The peculiarities of such heterogeneous dataset are discussed in section <ref>. The learning performances are analyzed in section <ref>: from the model choices within the input space, to the verification of the learned closure and its sensitivity analysis.
The validation of the model is primarily aimed at estimating the robustness gain obtained with this strategy. Based on the results presented in section <ref>, conclusions and future development are presented in section <ref>. 






§ DATA-DRIVEN AHFM


The data-driven Algebraic Heat Flux Model (AHFM) considered in this work is a physics-constrained neural network that predicts the turbulent heat flux at near unity and low Prandtl numbers. Details of its architecture are provided in <cit.> and briefly recalled in this section for completeness.

The network embeds rotational invariance properties by construction, as proposed by <cit.>, and its predictions are realizable, i.e. they satisfy the second law of thermodynamics. Figure <ref> provides a schematic of the neural network structure and its layers. The network predicts variable closure coefficients a_i and w_i based on a set of dimensionless invariants, denoted as π_i and defined on the left column of Table <ref>. Among the quantities listed in the table, k represents the turbulent kinetic energy, ϵ the turbulent dissipation rate, k_θ the variance of the thermal fluctuations, ϵ_θ the thermal turbulent dissipation rate. The quantities ν and α_l denote the molecular viscosity and diffusivity, respectively, while the tensors 𝐒, Ω and 𝐛 are defined as:
	
    𝐒 =1/2(∂𝐔/∂𝐱+∂𝐔/∂𝐱^T),



	
    Ω=1/2(∂𝐔/∂𝐱-∂𝐔/∂𝐱^T) ,


    𝐛=𝐮𝐮/k-2/3𝐈,

in which 𝐔 is the mean velocity field, 𝐮 the velocity fluctuation and the operator · denotes the Reynolds averaging. 

The coefficients a_i and w_i are used to compute the following expansions of tensors:

    𝐀   =∑_i=1^n a_i 𝐓^i,    𝐖   = ∑_i=1^n w_i 𝐓^i,

in which 𝐓^i are the element of the tensor basis indicated in Table <ref>. The dispersion tensor is then computed as a sum of a symmetric, positive-definite tensor and a skew-symmetric tensor:

    𝐃=[(𝐀+𝐀^T) (𝐀^T+𝐀) + k/ϵ^0.5 (𝐖-𝐖^T)],




In Ref. <cit.>, the network was trained with high-fidelity DNS data for Prandtl numbers ranging from 0.71 to 0.01. The loss function applied for the training promotes the smoothness of the predicted heat flux fields:


    ℒ(𝐪)= 1/N(∑_i=1^N ∑_j=1^3 (q̂_i,j -q_i,j)^2 )   
    
        + λ/N( ∑_i=1^N ∑_j,k=1^3  | ∂q̂_i,j/∂ x_k - ∂ q_i,j/∂ x_k| Δ x_k ),


where i∈[1,… N] is the index spanning across the N data points contained in each mini-batch, q_i,j=u_j θ is the prediction of the ANN and  q̂_i,j is the corresponding flux provided by DNS data. Both values q_i,j and  q̂_i,j are normalized with respect to the maximum DNS value achieved in each flow configuration.

Once trained, the data-driven AHFM was verified with some non-isothermal simulations run in OpenFoam, for which the model was compared with other theory-driven closures <cit.>. The verification highlighted the sensitivity of the data-driven AHFM to the momentum turbulence model applied in combination with it. As an example, Figure <ref>
presents the validation of the data-driven AHFM (denoted as ANN in the Figure) for turbulent channel flow, and the comparison with other theory-driven models such as the Manservisi model (MM) <cit.>, the Kays correlation <cit.> and the AHFM developed by Shams (AHFM) <cit.>. In the figures on the left, the AHFM is coupled with the Elliptic Blending Reynolds Stress Model (EBRSM) <cit.> while in the figure on the right the AHFM is coupled with the Launder-Sharma k-ϵ. Note that when the Reynolds stresses are accurately modelled with second order closures, the data-driven AHFM outperforms the traditional AHFM considered for comparison <cit.>. However, the accuracy deteriorates when the network receives Reynolds stresses modelled with the Linear Eddy Visocisty Model (LEVM), because of its dependency on the Reynolds stress anisotropy 𝐛, defined by eq. (<ref>).

Hence, the comparison shown in Figure <ref> highlights that the data-driven AHFM cannot handle deficiencies of the momentum modeling. This is a significant limitation of the machine learning closure, which restricts its applicability in CFD solvers.























§ MULTI-FIDELITY DATABASE AND PRELIMINARY ANALYSES

The database employed to train the data-driven AHFM consists of high-fidelity (DNS) data for non-isothermal turbulent channel flow <cit.> and non-isothermal backward-facing step <cit.> at various Reynolds and Prandtl numbers. An additional dataset of a non-isothermal planar impinging jet <cit.> was employed for testing. The details about the reference DNS data are provided in Table <ref>. In the present work, the training database is extended with RANS counterparts of the same flows simulated in the OpenFoam environment <cit.>. For the RANS simulations, the Launder-Sharma k-ϵ model was selected as the turbulence model. The computational domains were discretized to achieve y^+ ranging from 0.1 to 1.0. The RANS simulations were limited to the momentum field, i.e. the temperature field was not resolved. The 
reader is referred to the Appendix <ref> for further details about the computational setup for the training data.




The collected DNS and RANS data were manipulated to compute the tensors 𝐓^i and the invariants π_i indicated in Table <ref>. For turbulent channel flow, the discrepancy between the π_i computed with the two datasets is shown in Figure <ref>, representing the invariants in their 10-dimensional space with red (DNS) and green (RANS) lines. Specifically, each line in the plot indicates the values of the input features at each point of the computational mesh. Only 500 out of 425400 points constituting the training database are depicted for visualization purposes. 


As expected, the distribution of the invariants π_1 and π_2, π_7 and π_8 and Re_t is similar between RANS and DNS since these only involve isotropic turbulent quantities and the derivatives of the mean fields (𝐔 and T). On the other hand, significant differences appear on the invariants related to the anisotropic part of the Reynolds stress tensor 𝐛, i.e. π_3, π_4, π_5 and π_6. In fact, this behavior of the invariants reflects the differences in the tensor 𝐛 that can be appreciated in Figure <ref>, depicting the Reynolds stress anisotropy in the barycentric map <cit.> for the two datasets in case of turbulent channel flow.

This representation shows that the maximum distance between the turbulent states predicted by the two approaches occurs at y^+<30, at which the anisotropy is severely misrepresented. The gap decreases at higher y^+, for which the points move towards the state of pure isotropy. 

 





§ METHODOLOGY


The proposed methodology for a model handling heterogeneous input data combines dimensionality reduction, multiobjective optimization, sensitivity analysis, and uncertainty quantification. We describe the various steps in the following. Concerning the dimensionality reduction, Principal Component Analysis (PCA) was used to identify the optimal combination of features describing the different levels of fidelity as described in Section <ref>. Section <ref> reports on the multiobjective optimization at the foundation of the model training. The performances of the model obtained with this training are analyzed and compared with the original data-driven AHFM by inspecting its output layer and estimating its sensitivities in terms of Shapley values, introduced in section <ref>. Finally, the data-driven model is implemented in a CFD solver (OpenFoam), and there tested for a non-isothermal planar impinging jet at Pr=0.01. The details of this validation test case and the simulation setup are provided in section <ref>. 



 §.§ Principal Component Analysis (PCA)


The PCA was applied on the feature matrix 𝐗∈ℝ^N × F collecting the value of all the F features on the N grid points of both the training datasets. The reader is referred to <cit.> for an introduction to the PCA.

The goal was to identify P<F linear combinations of the F input features that explain most of the input variance and hence restrict the dimensionality of the model input's space along the associated principal directions, here denoted as 𝐯_i, with i=[1,..P]. These directions are eigenvectors of the covariance matrix 𝐒


    𝐒 = 1/N-1 (𝐗-𝐱_μ)^T (𝐗-𝐱_μ),
 where 𝐱_μ∈ℝ^N is the vector containing the average of each column. Denoting as 𝐕=[𝐯_1 …𝐯_P]∈ℝ^N× P the matrix collecting the leading P eigenvectors and Λ∈ℝ^P× P the diagonal matrix with the associated eigenvalues sorted in decreasing order, the reduced feature space can be written as 
𝐗̃=𝐗𝐕∈ℝ^N× P.

The correlations between the original variable and the principal components can be evaluated from the load matrix 𝐋=Cov(𝐗, 𝐗̃) ∈ℝ^F × P. Introducing the Singular Value Decomposition of the feature matrix 𝐗=𝐔√(Λ)𝐕^T, the load matrix can be computed as defined as 

 
    𝐋=Cov(𝐗, 𝐗̃) = 𝐗^T 𝐗̃/N-1 = 𝐕√(Λ)𝐔^T 𝐔√(Λ)/N-1 =  𝐕Λ/N-1.






 §.§ Multi-fidelity training


We denote as 𝐗_lf the feature matrix collecting data from low-fidelity simulations and as 𝐗_hf the one collecting data from high-fidelity simulations. To ensure reasonable predictions with both sets of data, the artificial neural network must be trained on the inconsistencies of the input data between the two fidelity levels. Denoting with 𝐲_lf and 𝐲_hf the predictions of the model with the inputs 𝐗_lf and 𝐗_hf, and with 𝐲̂_hf the reference values provided by the high-fidelity database, the cost function driving the training in case of multi-fidelity inputs is taken as:

    ℒ=ℒ(𝐲_hf,𝐲̂_hf)+αℒ(𝐲_lf,𝐲̂_hf),

in which ℒ(·) is defined by (<ref>). This choice of the loss translates the previous regression problem introduced in section <ref> into a multi-objective optimization problem in which the errors computed with both families of inputs need to be minimized. The scalar α in (<ref>) is a training hyperparameter that varies between 0.001 and 1000 to build the Pareto front of the the multi-objective problem. 




 §.§ Shapley values


We analyze the model sensitivity to its feature using Shapley values <cit.>. Shapley values, originating from cooperative game theory, offer a way to distribute the total gain (or cost) among the players (or features) based on their marginal contributions. For a model f with input features 𝐗=(𝐱_1,𝐱_2,…,𝐱_F), the Shapley value ϕ_i for the feature 𝐱_i is defined as:

    ϕ_i = ∑_S ⊆ K ∖{i}|S|! (|K| - |S| - 1)!/|K|![ f(S ∪{i}) - f(S) ],

in which K is the set including all F features,  S is a subset of N features not containing X_i, |S| is the cardinality of subset S, and f(S) denotes the model prediction using the subset of features S. The term f(S ∪{i}) - f(S) represents the marginal contribution of feature 𝐱_i when added to the subset S. The Shapley value calculation involves averaging these marginal contributions over all possible subsets S, weighted by the combinatorial factors |S|! (|K| - |S| - 1)!/|K|!. Note that, compared to other approaches for sensitivity analysis (e.g. the Integrated Gradient method <cit.> applied in <cit.>) this method provides more complete information about the sensitivity of the model, since the discrepancies in eq. (<ref>) is evaluated over the whole input space instead of targeting specific trajectories. 
Specifically, the predictions f(S) for each subset are computed by averaging the model evaluations obtained by perturbing the features excluded from the subset S around their predefined baseline values. In the present case, the baseline values of the ablated features is considered their average over the entire range of flows considered in the database indicated in Table <ref>. 

Note that based on eq. (<ref>) the computational cost of the Shapley method is given by the product of the number of features with their number of permutations, i.e., F · F!, which would lead to 18· 18! model evaluations based on the list of features reported in Table <ref>. To reduce the computational cost, the features were assigned to four main groups that will also ease the interpretation of the results. These groups, reported in Table <ref>, are:  



  * Group 1: Momentum isotropic features (MI) depending on the velocity gradient and isotropic momentum statistics (e.g. k, ϵ); 

  * Group 2: Momentum anisotropic features (MA) depending on the anisotropic part of the Reynolds stress tensor 𝐛; 

  * Group 3: Thermal-based features (TH) depending on the molecular Prandtl number and thermal related statistics;

  * Group 4: The basis tensors (TE) indicated in Table <ref>. 







 §.§ Uncertainty propagation

The uncertainty analysis of the trained data-driven AHFMs was carried out to evaluate their robustness with respect to the momentum treatment. The method is based on the perturbation of the turbulent state in the barycentric map introduced by Figure <ref> and, specifically, between the true state (DNS) and the one predicted by Linear Eddy Viscosity Models (LEVMs). The coordinates of the points in the map are computed from the eigenvalues ψ_i (i=1,2,3) of the anisotropic part of the Reynolds stress tensor 𝐛. The vertices of the triangle (𝐳_1C,𝐳_2C,𝐳_3C) represent the limiting states of turbulence:


  * one component (𝐳_1C=[1,0]) for which ψ_i = 2/3, -1/3, -1/3.

  * two components (𝐳_2C=[0,0]) for which ψ_i = 1/6, -1/6, -1/3.

  * isotropic (𝐳_3C=[1/2,√(3)/2]), for which ψ_i are all zero. 


Each point 𝐳∈ℝ^2 in this plane is associated to a specific set of eigenvalues ψ∈ℝ^3 by a linear mapping <cit.>:


    𝐳^*= 𝐳_1C(ψ_1-ψ_2) + 𝐳_2C(2ψ_2-2ψ_3) + 𝐳_3C (3ψ_3 + 1 ),
 complemented with the conditions: 

    ψ_1 + ψ_2 + ψ_3 = 0,


Writing the mapping ψ→𝐳 as 𝐳=𝐁ψ, its inverse becomes ψ=𝐁^-1𝐳.

This mapping was here used to propagate uncertainties in the baricentric map to uncertainties in the associated anisotropic stress and then the turbulent heat flux through the model. To this end, the high-fidelity (DNS) state 𝐳_hf in the barycentric map is perturbed towards the state of turbulence given by the low fidelity (RANS) momentum treatment 𝐳_lf, leading to a modified location 𝐳^*:

    𝐳^* = 𝐳_hf + Δ (𝐳_lf-𝐳_hf),

in which Δ is a uniform random variable with range [0,1]. An example of a perturbed state for turbulent channel flow (Re_τ=640) is given in Figure <ref>.


Based on the perturbed location 𝐳^* in the barycentric map, the new vector of eigenvalues ψ^* is computed as

    ψ^* = 𝐁^-1𝐳^* .
 

The associated perturbed anisotropy tensor is then computed as:

    𝐛^* = 𝐄diag(ψ^*) 𝐄^-1,

in which 𝐄 is the matrix of eigenvectors of the unperturbed tensor 𝐛. The tensors 𝐓^i and invariants π_i are computed based on the perturbed 𝐛^* and propagated through the model.






For each point of the original DNS dataset, 50 perturbed states are sampled and submitted to the data-driven models trained in single fidelity mode <cit.> and in the multi-fidelity mode described in section <ref>. The predictions from both models are averaged, and confidence intervals are constructed to assess how uncertainties in the Reynolds stress components propagate through the data-driven closure.




 §.§ Verification test case

To test the performance of the data-driven AHFMs generated with single fidelity <cit.> and multi-fidelity training modes, the networks were implemented in a CFD solver (OpenFoam) following the procedure proposed by Maulik et al. <cit.>. The test case proposed is a non-isothermal planar impinging jet simulated by Duponcheel et al. <cit.> at Pr=0.01 for the geometry depicted in Figure <ref>. 



This test case is ideal for testing the data-driven closures, as this flow configuration is excluded from the training database. Additionally, for this flow configuration, the momentum field obtained with low-order turbulence models (e.g., k-epsilon, k-omega) is affected by large uncertainties due to the extent of anisotropy, secondary recirculation, and transport effects. 

The computational domain for this test case consists of two infinite parallel plates with a slit in the middle injecting the flow conditions reached in a fully developed turbulent channel flow at Re_τ=181 (see Figure <ref>). The heat transfer is triggered by the temperature gap imposed on the walls. The structured computational grid leads to y^+ values ranging from 0.32 to 1.6 in the first cell near the wall, allowing a wall-resolved treatment for both momentum and thermal fields. The boundary conditions imposed are summarized in Table <ref>. 






 







The simulations employed the Launder-Sharma k-ϵ and k-ω SST as turbulent momentum closures. The energy equations was closed with the analytical model developed by Manservisi <cit.>, the original data-driven model introduced in section <ref> and the model generated with the multi-fidelity training mode. 




§ RESULTS

This section presents the results of the analysis of the data-driven AHFM discussed in section <ref>, comparing the model trained with multi-fidelity inputs to the one trained exclusively with DNS data. For brevity, we refer to the former as the Hybrid ANN and the latter as the High-Fidelity ANN.






 §.§ Results of the PCA

The PCA was carried out on a feature matrix combining both the high-fidelity and the low-fidelity data. 






The variance explained by each principal component is illustrated in Figure <ref>, which also reports the cumulative sum of the variance explained. These results show that the first three principal components are able to explain 77.5% of the total variance of the input data. 

The loading matrix (cf. (<ref>)) was computed to analyze the principal components' contribution to the original features' variance. The loadings of the first three principal components are shown in Figure <ref>. The loading values show that the first component (PC1) is mainly composed of momentum features depending on the mean velocity gradients and isotropic turbulence statistics (e.g., k, ϵ), which are similar in both RANS and DNS datasets. The second principal component (PC2) is primarily constituted by thermal-based features depending on the Prandtl number of the fluid. The third principal component (PC3) depends on momentum-based features that describe the turbulence anisotropy through the deviatoric part of the Reynolds stress tensor. These latter features are misrepresented in the RANS input dataset due to the Boussinesq approximation. Interestingly, the PCA algorithm splits the entire input dataset in three clusters, including isotropic momentum-based features (PC1), thermal-based features (PC2), and anisotropic momentum-based features (PC3). 





The two datasets (RANS and DNS) are projected onto the PCs axes and compared in Figures <ref> and <ref>. In particular, Figure <ref> compares the two datasets in terms of PC1 and PC2. The DNS inputs achieve higher peaks of PC1 for the same PC2 values due to the higher turbulent kinetic energy retrieved in DNS simulations compared to the RANS counterpart. However, the distributions of the data on this plane are similar. A huge separation between the two datasets is observed when projected on the PC1 and PC3 plane, as done in Figure <ref>. This data separation reflects the effect of the Reynolds stress modeling and suggests that the current parameter space allows the detection of the type of momentum modeling based on the distance between critical features (mainly π_5 and π_3). 
This cluster separation is crucial to the present work, as it demonstrates that specific features constituting PC3 can be utilized by a data-driven algorithm to detect the fidelity level of the momentum turbulence model and to adapt predictions accordingly.











 §.§ Performance of the training with multi-fidelity inputs

The neural network architecture introduced in section <ref> was trained with the multi-fidelity approach introduced in section <ref>. Specifically, for each value of α, 30 statistical trainings were carried out. The Pareto front of the multi-objective problem is reported in Figure <ref>, which is built from the final values achieved by the two losses at the end of each training. It is worth stressing that the resulting Pareto front is sharp, meaning that the training finds solutions that do not compromise the accuracy of the heat flux predictions with both types of inputs. This is an interesting result because it proves that increasing the robustness of the data-driven AHFM to the momentum modeling is possible. 




Regarding accuracy, a comparison between the hybrid ANN and the high-fidelity ANN is given in Figure <ref>, which relates the network output obtained with DNS and RANS predictions. Note that hybrid ANN significantly reduces the dispersion of the heat flux predictions with the change of the momentum treatment. 





 §.§ Results of the sensitivity analysis

The Shapley values defined in (<ref>) were computed considering the streamwise and wall-normal components of the heat flux as the outputs of the model. 
Figures <ref> and <ref> presents the obtained Shaply values for Pr=0.71 and Pr=0.025, respectively. The values obtained in the presence of DNS inputs are indicated with solid lines, of RANS inputs with dashed lines. At near unity Prandtl numbers, the data-driven model is shown to be highly sensitive to the thermal features. This sensitivity reduces at Pr=0.025, meaning that the dependence of the model on the thermal gradients and features increases with the Prandtl number. Clearly, the ANN trained with the multi-fidelity approach is insensitive to the anisotropy-based invariants in the presence of RANS inputs, though it preserves a small sensitivity to this group when the model is fed with high-fidelity momentum data. 

The comparison of the Shapley values in the presence of DNS and RANS data (Figure <ref> and <ref>) shows that the sensitivity of the model to all the groups of features decreases when low fidelity data are detected. Hence, the hybrid ANN relies more on high-fidelity data than low-fidelity ones, i.e., the robustness in case of inaccurate Reynolds stress modeling is higher. The profiles are in fact smoother, characterized by lower peaks and fewer changes of sign. This behavior means that the model derived with multi-fidelity training simplifies when low-fidelity input data are detected. The network adapts to the uncertain inputs by reducing its sensitivity and non-linearities.







Figure <ref> depicts the distribution of the Shapley values of the original, high-fidelity ANN for the same flow configuration at Pr=0.71. The comparison between Figure <ref> and <ref> highlights that the hybrid training mode reduces the sensitivity of the ANN to the features belonging to the last two groups indicated in Table <ref>. In particular, the sensitivity drop is evident for y^+ ≈ 10 where the momentum turbulent production reaches its peak. The barycentric triangle in Figure <ref> shows that, at this distance from the wall, the gap between RANS and DNS turbulent states is maximum. Hence, the hybrid training mode moderates the sensitivity to the anisotropic part of the Reynolds stresses where the separation between the two categories of inputs is significant. This confirms that this training strategy effectively leads to a more robust thermal model with respect to the one trained with high-fidelity data only, which would better interface with standard momentum closures based on the eddy viscosity concept.   



 §.§ Layer output analysis

Additional insights into the differences between the models obtained from single-fidelity and multi-fidelity training can be gained by analyzing the output layer of the network schematized in Figure <ref>. Figure <ref> depicts the coefficients a_i and w_i at the output of the merge layer (see Figure <ref>) of the ANN trained with the multi-fidelity approach, in case of non-isothermal turbulent channel flow at Re_τ=640 and Pr=0.025 and with high-fidelity and low-fidelity momentum input data.[Note that the sign of the coefficients a_i does not alter the final output, since 𝐀=∑ a_i 𝐓^i is the Cholesky factorization of the symmetric part of 𝐃, as indicated by eq.(<ref>).]. For both kinds of input data, the dominant coefficients are a_1, a_2, a_6, w_4, and w_6. For y^+>100, the values of the predicted coefficients are similar for both DNS and RANS inputs. At y^+≃ 10, where the turbulence anisotropy reaches its peak, the predicted coefficients are much higher in the presence of DNS input data than of RANS inputs. This explains the lower sensitivity to the RANS inputs detected by the Shapley value analysis in section <ref>. The reduction is significant, especially for a_2, the coefficient multiplying the tensor 𝐓_2 that is proportional to the anisotropic part of the Reynolds stress tensor 𝐛, as reported in Table <ref>. 







The output of the last layer of the high-fidelity network is significantly different, as shown in Figure <ref>. The coefficients a_1 and a_2 dominate above all the other terms and achieve a single peak in the region of maximum momentum turbulence production (y^+ ≈ 10). Hence, compared to the high-fidelity network, the new model increases the dependence on 𝐓_4,  𝐓_6, and  𝐓_8 which depend on the strain and rotation tensors. This implies that training the network with data of multiple fidelity replaces the dependence on the true Reynolds stress tensor with linear and quadratic functions of the velocity gradient. In other words, the network tries to reconstruct the true Reynolds stress anisotropy as a function of the mean velocity gradient and its uncertain estimate given by the combined momentum turbulence model. This explains the low dispersion of the predictions shown in Figure <ref> with both kinds of momentum treatments. 



 §.§ Propagated uncertainties

The uncertainties of the Reynolds stress tensor caused by the inaccuracies of the combined momentum treatment were propagated to the output of the data-driven AHFM with the method explained in section <ref>. The statistics obtained for Re_τ=640 and Pr=0.71 and 0.025 are shown in figures <ref> and <ref>, respectively. The confidence intervals, highlighted with shaded areas, quantify the uncertainty of the predictions due to the inconsistency between reference and modeled input data. The comparison between the hybrid ANN and the high-fidelity ANN reveals a significant reduction of the model uncertainty when the training is enriched with RANS data. This is true at both near unity and low Prandtl numbers, albeit more evident for Pr=0.71, at which a stronger dependence of the heat flux field on the Reynolds stresses is expected.

The standard deviation computed for the two models at both values of the Prandtl numbers is shown in Figure <ref>. The standard deviation of the high-fidelity model achieves its peak in the near wall region, where the dependence of the model on the anisotropy state is maximum. In this region, the standard deviation of the hybrid model predictions is much lower since the training with both data mitigates the dependency on the Reynolds stress anisotropy. This evidence agrees with the considerations drawn for the Shapley values for y^+ close to 10.0. This region of peak turbulence production, where temperature gradients are most significant, and the influence of the temperature field on heat flux is at its maximum, is critical. Therefore, reducing the uncertainty in model predictions in this area is of utmost importance. These considerations lead to the conclusion that the hybrid ANN is significantly more robust than the high-fidelity ANN when it comes to handling inaccuracies in the combined momentum turbulence model.

 

   


 




 §.§ Results of the simulation test case

The numerical setup described in section <ref> was employed to verify the behavior of the hybrid model when combined with LEVMs, and to compare the new formulation with the original data-driven AHFM. As anticipated in section <ref>, for this test case, the computed velocity field highly depends on the type of momentum treatment and the choice of the specific LEVM. This is confirmed by Figure  <ref>, which presents the velocity fields computed with the k-ϵ and k-ω models. The recirculation bubble's size and shape generated after the impinging point significantly differs from the two closure models. The difference can also be appreciated from the velocity and k profiles reported in Figure <ref> and <ref>, where they are compared with the reference DNS. In particular, the k-ϵ and k-ω models underestimate and overestimate the recirculation bubble's extension.  








The data-driven thermal model and the Manservisi model <cit.> were applied over the two underlying momentum fields to compute the heat flux and temperature distributions. The high-fidelity ANN, whose predictions are indicated with dashed blue lines, completely mismatches the reference heat flux profiles due to the inconsistency between training (DNS) and simulation (RANS) momentum data. The hybrid ANN gives instead accurate heat flux predictions, regardless of the combined momentum closure (k-ϵ or k-ω). These results show that the applicability of the hybrid ANN is not restricted to the momentum closure that generated the training data (Launder-Sharma k-ϵ), but could potentially extend to the entire family of LEVMs. The heat flux fields given by the hybrid ANN are close to the one computed with the Manservisi model, although a slight improvement in the accuracy can be appreciated close to the slit. 

The comparison among the thermal models in terms of temperature distribution is presented in Figure <ref>, in which the hybrid ANN combined with the k-ω SST shows the best accuracy, while the hybrid ANN and the Manservisi model significantly overestimate the temperature, especially far from the slit. However, the better thermal field computed with the k-ω and hybrid ANN does not originate from a significantly more accurate heat flux modeling. Indeed, the comparison of the profiles in Figure <ref> does not indicate the superiority of the heat flux predictions compared to the other two setups (k-ϵ-hybrid ANN and k-ϵ-Manservisi model). The better temperature agreement obtained with the high-fidelity ANN than with the hybrid ANN with the same momentum closure suggests that this temperature field is a result of compensation errors of both momentum and thermal closures. This consideration further underlines the importance of combining accurate thermal and momentum closures for heat transfer RANS simulations, as both the turbulent treatments contribute substantially to the resulting thermal field. This is especially true for low Prandtl numbers, at which the turbulent heat transport can be comparable, or even lower, than molecular transport and convection contributions in the overall energy balance. 









§ CONCLUSIONS AND OUTLOOK


This article investigated the problem of model-data inconsistency for a data-driven thermal turbulence model trained with high-fidelity data. Specifically, we explored ways to improve the robustness of the data-driven model to extend its coupling with low-fidelity models for Reynolds stresses (LEVMs), generally preferred in industrial contexts for their stability and computational cost. 

The crucial point was understanding whether the modeled (RANS) input statistics could be employed in training to inform the data-driven model about the inconsistencies between low and high-fidelity data. Detecting the quality of the momentum treatment could allow to adapt to different inputs or mitigate the sensitivities to the most critical ones. For such purpose, the work proposes the analysis of the multi-fidelity input space and a training strategy based on a hybrid dataset consisting of DNS and RANS input data. 

The work demonstrates that increasing the robustness of the thermal model with respect to momentum modeling is possible. The systematic deviation between high-fidelity and low-fidelity statistics can be utilized to identify the type of momentum modeling and adapt the output accordingly or to find alternative relationships of the input statistics. The multi-fidelity training mode generates a model less sensitive to Reynolds stress anisotropy than the model trained with only high-fidelity data, especially in regions characterized by high anisotropy. The network naturally adapts its sensitivity to low-fidelity inputs (RANS) and is thus more robust than the original. From a structural standpoint, the hybrid model appears to incorporate Reynolds stress anisotropy modeling through algebraic expansions of tensors dependent on velocity gradients.

The validation of the new data-driven model for a non-isothermal planar impinging jet shows that it is significantly more robust than the original high-fidelity network and provides accurate predictions with different eddy viscosity models. The uncertainty propagation analysis demonstrates, more generally, the greater robustness of the new model against perturbations of the Reynolds stresses in the barycentric plane.

While more test cases are needed to develop a truly general-purpose model, the promising results of this work open the door to new training strategies for data-driven turbulence closures that can leverage information from databases of varying sizes and fidelity levels. As demonstrated here, multi-fidelity inputs enhance model robustness and resilience to perturbations. Additionally, incorporating multi-fidelity targets can expand the training database to include conditions where high-fidelity data are challenging to obtain, such as high Reynolds or Grashof numbers or flows involving multiple physics or complex geometries. This approach can improve the generality and applicability of the learned closures. Future work will focus on extending the current framework to a broader range of data, both in inputs and outputs, to further advance data-driven turbulence modeling and its application in industrial contexts.







§ NOMENCLATURE




[t]0.5



𝐔     [m/s]     Mean velocity 

T	    [K]     Mean Temperature

𝐮     [m/s]     Velocity fluctuation 

θ	    [K]     Thermal fluctuation

k     [m^2/s^2]	    Turbulent kinetic energy

ϵ     [m^2/s^3]     Turbulent dissipation rate

k_θ     [K^2]     Thermal variance

ϵ_θ	    [K^2/s]     Thermal dissipation rate

𝐮𝐮	    [m^2/s^2]     Reynolds stress 

 

 




[t]0.5



𝐮θ	    [mK/s]    	Turbulent heat flux


ν	    [m^2/s]     Molecular viscosity

α_l	    [m^2/s]     Molecular diffusivity


δ	    [m]     Half channel width


𝐈	    [-]     Identity tensor

𝐒	    [1/s]     Strain rate tensor

Ω	    [1/s]     Rotation tensor

𝐛	    [-]     Reynolds stress anisotropy tensor

 










This work was supported by an F.R.S.-FNRS FRIA grant, and the authors gratefully acknowledge Prof. Iztok Tiselj and Dr. Mathieu Duponcheel for providing their datasets.







§ RANS SETUP FOR THE TRAINING FLOWS
 
The simulations of the RANS counterpart of the flows in the training database were conducted in the openfoam environment with the Launder-Sharma k-ϵ model. The choice of the computational domains follow the setups described in Ref. <cit.> and <cit.>. For both setups, the mesh is structured and consistent with the wall-resolved approach, i.e., the values of y^+ range between 0.1 and 1.0. Table <ref> and <ref> indicate the boundary conditions applied for the two flows. 

For the channel, the fully developed flow is obtained by imposing cyclic boundary conditions at the inlet and outlet and a pressure gradient source term that adapts to the bulk velocity prescribed. For the backward-facing step flow, a fully developed flow in the channel preceding the step is obtained using recycling conditions, i.e., the iterative remapping of the fields at a certain distance from the inlet to the inlet boundary.  















