The critical heat flux (CHF) corresponding to the departure from nucleate boiling (DNB) crisis is a regulatory limit for the licensing of pressurized water reactors (PWRs) worldwide. Despite the abundance of predictive tools available to the reactor thermal-hydraulics community, the path for an accurate CHF model remains elusive. This work approaches the prediction of DNB through a physics-informed machine learning-aided framework (PIMLAF) with the objective of achieving superior predictive capabilities for a rod bundle. In view of the limitations in existing macro-scale physics-driven tools, an improved mechanistic model is first proposed, leveraging key concepts in the liquid sublayer dryout and bubble crowding mechanisms. The proposed mechanistic model is able to predict DNB in different heater geometries for a broad range of flow conditions without the need for recalibration. This model is then incorporated as the physics-informed component of the hybrid framework PIMLAF, which takes advantage of established understanding in the field (i.e., domain knowledge [DK]) and uses machine learning (ML) to capture undiscovered information from the mismatch between the actual and DK-predicted output. Two bundle-related case studies using the PWR subchannel and bundle tests (PSBT) database are carried out to illustrate the PIMLAF’s improved performance over traditional approaches for both interpolation and extrapolation purposes. In light of the PIMLAF’s promising potential to reduce prediction error, reactor vendors are encouraged to leverage their in-house experimental efforts and apply the hybrid framework to potentially achieve margin reductions in the minimum DNB ratio (MDNBR) for the designs of interest. Keywords: CHF; DNB; Physics-informed machine learning-aided framework; Mechanistic model; Rod bundle; MDNBR. 3 Nomenclature 𝐴 channel cross-sectional area [m2] 𝐶! drag coefficient [-] 𝑐" specific heat [J/kg-K] 𝐷 channel diameter (tube) [m] 𝐷# vapor blanket thickness (or equivalent diameter) [m] 𝐷$ bubble departure diameter [m] 𝐷% channel hydraulic (or equivalent) diameter [m] 𝐷& channel heated diameter [m] 𝐹 force [N] 𝑓 friction factor [-]; function 𝐺 mass flux [kg/m2-s] 𝑔 gravitational acceleration [m/s2] ℎ'( latent heat of vaporization [J/kg] 𝐿 length [m] 𝑚̇ mass flow rate [kg/s] 𝑃 pressure [Pa] 𝑞)) heat flux [W/m2] 𝑅𝑒 Reynolds number [-] 𝑇 temperature [K] 𝑈 velocity [m/s] 𝑈* friction velocity [m/s] 𝒙 input feature vector 𝑥 flow quality [-] 𝑥% equilibrium quality [-] 𝑦 distance from heated wall (to blanket centerline) [m]; target/actual output 𝑦+ non-dimensional distance from heated wall [-] 𝑦6& predicted output with hybrid framework 𝑦6" predicted output with prior model 4 Greek symbols 𝛼 void fraction [-] 𝛿 liquid sublayer thickness [m] 𝛿, bubble region thickness [m] 𝜀 surface roughness [m]; residual 𝜀̂- predicted residual with ML 𝜂 fraction of cross-sectional area [-] 𝜆 critical wavelength [m] 𝜇 dynamic viscosity [Pa-s] 𝜌 density [kg/m3] 𝜎 surface tension [N/m] 𝜏. wall shear stress [Pa] Subscripts 𝐵 vapor blanket 𝑏 bubble region 𝑐 core region 𝑐ℎ channel-averaged 𝑒𝑥𝑝 experimental 𝑓 liquid at saturation 𝑔 vapor at saturation 𝑖𝑛 channel inlet 𝑙 liquid 𝑜𝑢𝑡 channel outlet 𝑟 relative 𝑠𝑎𝑡 saturation 𝑠𝑢𝑏 subcooling 𝑣 vapor 5 Acronyms ANN artificial neural network BWR boiling water reactor CFD computational fluid dynamics CHF critical heat flux DK domain knowledge DNB departure from nucleate boiling EPRI Electric Power Research Institute IQR interquartile range LUT look-up table MDNBR minimum departure from nucleate boiling ratio ML machine learning NN (feed-forward) neural network NRC Nuclear Regulatory Commission (U.S.) PCC Pearson correlation coefficient PIMLAF physics-informed machine learning-aided framework PSBT PWR subchannel and bundle tests PWR pressurized water reactor ReLU rectified linear unit RF random forest RIA reactivity-initiated accident rRMSE relative root-mean-square error SA sensitivity analysis SRCC Spearman rank correlation coefficient T/H thermal-hydraulics 6 1. Introduction The departure from nucleate boiling (DNB) crisis in a two-phase flow boiling system is widely encountered in thermal management of high power-density electronic devices, in the refrigeration industry, in medical technology fields, and more typically in nuclear power plants. In a pressurized water reactor (PWR), DNB is perceived as a grand challenge to the integrity of fuel components but remains one of the least understood thermal-hydrodynamic instabilities [1]. Under DNB conditions, the heated surface is blanketed by a stable vapor film, transforming high-efficiency nucleate boiling into low-efficiency film boiling and leading to a potentially catastrophic escalation of the heater temperature. The corresponding heat flux, or critical heat flux (CHF), is a regulatory limit for the licensing of commercial PWRs worldwide [2]. The minimum departure from nucleate boiling ratio (MDNBR)—the minimum value of the ratio between the CHF and the actual operating local heat flux—defines various PWR thermal limits and margins accounting for correlation uncertainties, engineering and process parameter uncertainties, core anomalies, and transients. While the numerical values for MDNBR depend on the individual reactor design and on the specific CHF model utilized, the U.S. Nuclear Regulatory Commission (NRC) requires that the MDNBR be at least 1.3 at 112% power transient [2]. The importance of CHF has led to intensive experimental and theoretical investigations over several decades. Unfortunately, the path for its accurate prediction has been elusive due to lack of consensus on the mechanism that triggers DNB [3]. This lack of established agreement arises from an incomplete understanding of the very complex nature of two-phase flow and heat transfer, which was caused primarily by the absence of experimental capabilities to perform measurements at the space and time scales characteristic of the DNB phenomenon [4]. Nevertheless, the reactor thermal-hydraulics (T/H) community has proposed a large number of predictive tools for CHF, ranging from data-driven correlations and look-up tables (LUTs) to physics-driven mechanistic models. On one hand, the best-fit correlations and LUTs are easy to implement (thus having been used extensively in system and subchannel codes) but are entirely empirical. They may result in relatively good agreement with specific experimental datasets, but they often fail to extend beyond their ranges of validity in terms of geometry and operating conditions. Among the commonly used empirical tools in flow boiling are: the Tong-68 correlation (developed from tube data) [2], the 7 Biasi correlation (tube) [5], the Groeneveld 2006 LUT (tube) [6], the W-3 correlation (rod bundle) [7], and the Electric Power Research Institute (EPRI) correlation (rod bundle) [8]. On the other hand, the mechanistic models rely on assumptions made based on reasonable yet limited understanding of the underlying physics and are supplemented with mostly empirical constitutive relations to close the conservation equations. Since the mid-1960s, numerous macro-scale mechanistic models in flow boiling have been developed and can be generally grouped into six categories [9] based on their respective DNB triggering mechanisms, among which the liquid sublayer dryout and the near-wall bubble crowding have received considerable attention [9–11]. Unfortunately, all the established macro-scale models were originally developed for circular flow channels in round tubes, and very few efforts have been made recently to extend them to rod bundle subchannels. The rod bundle experiments, especially those in large-scale simulation of PWR fuel assemblies, have been relatively scarce and are mostly proprietary. In terms of measuring techniques, the vast majority of published CHF data were based on the heater temperature measurements obtained using thermocouples [9]. Recent advancement of infrared thermometry diagnostics has allowed for high-resolution measurements of temperature and heat flux distributions on the boiling surface [4,12]. However, the applications have so far been limited to pool boiling or flat plate heaters at low pressure conditions [12,13]. Such limitations also hinder adoption of higher resolution micro-scale models using computational fluid dynamics (CFD) because their validation requires high-fidelity experiments at the prototype scale and operating parameters. With advances in computational capabilities and optimization techniques, predictive methods based on machine learning (ML) have become a promising alternative to existing data-driven tools. The ML approach can be particularly useful for engineering fields in which complex physical phenomena are challenging to model. Within this category, an artificial neural network (ANN) is one of the most popular choices due to its superior performance when dealing with strongly nonlinear relations [14,15]. Depending on the expected outcomes, different types of ANNs can be applied to a variety of disciplines. To predict CHF, the use of a feed-forward neural network (or multilayer perceptron, or simply NN) is described in several relevant publications [16–19], but 8 none of these efforts have attempted to differentiate DNB from dryouta when selecting appropriate data sources; nor have they employed cross-validation techniques to assess their neural network architectures. Another popular ML branch for regression—tree-based ensemble learning such as random forest (RF)—is highly efficient and robust [20], although its applications in nuclear engineering have been scarce to date. While standalone ML-based predictive tools require minimal prior knowledge and almost no explicit mathematical modeling, they can be prone to largely scattered, physically undesired solutions (especially for sparsely distributed datasets) due to their purely data-driven “black-box” nature and lack of interpretability. Regardless of how advanced any ML algorithm has become, prior knowledge in the field of interest—or domain knowledge (DK)—is still deemed important and useful by many researchers [21–25]. Accordingly, the idea of combining ML and DK arose through the concept of a hybrid/integrated framework; this initially was seen in chemical process industries [21–24], and later it was considered for applications in electrical and aerospace engineering [25,26]. The hybrid framework can be further classified into two main types: the series approach, in which ML is applied to estimate intermediate variables that are key closures in a DK-based mechanistic model; and the parallel approach—also known as the physics-informed ML-aided framework (PIMLAF)—which takes advantage of the credible baseline solutions provided by established DK and uses ML to compensate for the bias (i.e., to capture undiscovered information from the mismatch between the actual and DK-predicted target). According to the open literature, the hybrid framework—especially the parallel approach—has seen little use for nuclear engineering applications. Zhao et al. [27] recently proposed to use the PIMLAF for predicting DNB-type CHF in tubes, annuli and one-side heated plates. The objective of this paper is to extend prior work on mechanistic modeling of DNB and to achieve superior predictive capabilities for rod bundle applications through the use of the hybrid framework PIMLAF. Section 2 describes the PIMLAF methodology for the prediction of CHF and highlights the role of domain knowledge in the framework. Section 3 surveys existing macro-scale a The term CHF is sometimes confusingly used for both DNB and dryout. While DNB occurs at low-void fractions typical of PWR-operating conditions with subcooled or slightly saturated bulk flow, the dryout crisis is triggered in annular flow at high-void fractions and lower-than-DNB heat fluxes that are typical of operating conditions in a boiling water reactor (BWR). 9 physics-driven approaches and proposes an improved model that leverages a relatively well-established understanding of the DNB triggering mechanism. Model optimization and sensitivity analysis are conducted. Section 4 incorporates the proposed mechanistic model as the physics-informed component of the hybrid framework and uses the PWR subchannel and bundle tests (PSBT) database [28] to illustrate the PIMLAF’s improved performance over its standalone DK or ML counterparts for both interpolation and extrapolation purposes. The implication of a more accurate CHF prediction for the MDNBR margin reduction is briefly discussed. Finally, Section 5 summarizes the key findings and presents recommended future research tasks. 2. PIMLAF methodology As discussed in Section 1, the conventional DK-based approach has its limits. The data-driven empirical correlations and LUTs often fail to extend beyond their ranges of validity due to limited calibration data, and correlation/table updating can be cumbersome and error-prone when new data become available, especially for LUTs. The physics-driven models rely on assumptions deduced from a limited understanding of the underlying physics, and their mechanistic foundation is inevitably diluted by making use of best-fit correlations for the closure parameters. Due to lack of modeling and restriction from laws of physics, the standalone ML-based approach is prone to largely scattered, physically undesired solutions. For regression problems, some researchers claim that ML is simply an advanced data fitting tool that offers trivial novelty, while others consider it as an omnipotent panacea [22]. The reality no doubt lies somewhere in between. One point that most people would agree upon is that domain knowledge should not be thrown away when ML is employed, no matter how advanced any ML algorithm has become. In recent work by Zhao et al. [27], a fresh approach to the prediction of DNB has been presented through the concept of PIMLAF—physics-informed ML-aided framework—in which ML is used to assist (and not replace) DK. In this hybrid parallel framework,b an established DK-based data- or physics-driven model is selected as the fixed-structure prior model. While DK b As discussed in Zhao et al. [27], regarding the prediction of DNB, the hybrid series framework could be applied to a liquid sublayer dryout mechanistic model, where ML would be used to forecast such closure terms as bubble velocity, drag coefficient and bubble departure diameter. At this stage, a much more comprehensive database is required for proper ML training on any of these intermediate variables, so the series approach is not explored in this study. 10 serves to lay the groundwork and provide a baseline solution, ML is used to learn from the residual between actual and DK-predicted output (and not from the final output, as in the case of a standalone ML method). Figure 1 shows the simplified workflow of the PIMLAF during training and validation/testing. Input information is available to both components of the prior and ML model, although their features do not need to be identical. The prior model is essentially a nonlinear function (𝑓) of input features (𝒙) in the form of an empirical correlation/table or a system of first-principles equations. In the training stage (Fig. 1a), the residual (𝜀) is obtained by subtracting the prior model’s predicted output (𝑦6"=𝑓(𝒙)) from the actual/measured output (𝑦), the experimental CHF, and then ML is trained on the residual. The ML-predicted residual (𝜀̂-) is compared with 𝜀 through a loss function. The objective of the training process is to minimize the loss function, which is presented as the mean squared error in this work. The final predicted output (𝑦6&) by the PIMLAF is the sum of 𝑦6" and 𝜀̂-. Similarly, in the validation or testing stage (Fig. 1b), the prior model and the trained ML model are arithmetically combined to make predictions on unseen data. (a) training (b) validation/testing Fig. 1. Simplified workflow of the hybrid framework PIMLAF [27]. 11 The feed-forward neural network (NN) and the random forest (RF) are adopted for use in the ML-aided component of the hybrid framework. The NN is essentially a collection of multilayer,c fully-connected units (see Fig. 2a for a simplified architecture) capable of nonlinear mapping via activation functions between two layers. Weights and biases are randomly initiated from a uniform distribution and are iteratively updated during training by the back-propagation algorithm following a gradient descent approach that exploits the chain rule. To ensure model generalization and to prevent overfitting, various regularization techniques can be applied to help reduce the test error. The RF is a highly efficient, flexible tree-based ensemble learning method that produces robust results with minimum hyperparameter tuning. By randomly selecting subsets of observations and features using the bootstrap aggregation technique (also known as bagging, a model averaging approach designed to improve accuracy and reduce variance), multiple decision trees are combined and their predictions are then averaged, as depicted in Fig. 2b. (a) NN (b) RF Fig. 2. Simplified ML structures [27]. More details on the specifics of both ML methods—including feature engineering, hyperparameter tuning, cross-validation, and model evaluation—are provided in Zhao et al. [27], which presents an assessment of the PIMLAF performance with a CHF test matrix covering a broad range of flow conditions for DNB in tubes, annuli, and one-side heated plates. This work aims to focus on DNB prediction in rod bundles. As prior knowledge has been proved to be of c An NN consists of an input layer, an output layer, and intermediate layers which are also called hidden layers. The shallow NNs have only one hidden layer as opposed to the deep NNs which have multiple hidden layers. In general, the power of the NN increases with the number of hidden layers at the cost of computational complexity. input layerhidden layer 1hidden layer 2output layer 12 essential importance in the hybrid framework—especially for extrapolation and when the training dataset is relatively small—a robust DK-based model is deemed necessary. An improved mechanistic DNB model will be proposed in Section 3 and then integrated into the PIMLAF in Section 4. The resulting hybrid models will be evaluated against their standalone DK and ML counterparts for both interpolation and extrapolation purposes through two bundle-related case studies. 3. Improved mechanistic DNB model 3.1. Overview of previous research Over the past several decades, numerous macro-scale mechanistic CHF models in flow boiling have been proposed. Based on their respective DNB triggering mechanisms, the major approaches can be grouped into six categories [9], as briefly summarized below in chronological order: 1) The liquid layer superheat limit mechanism (1965) assumes that CHF occurs when the liquid layer in the immediate vicinity of the heated surface reaches its critical superheat. 2) The boundary layer separation mechanism (1968–1975) assumes that CHF takes place when the rate of vapor steaming increases beyond a critical value and induces stagnation of the liquid flow close to the heater. 3) The liquid flow blockage mechanism (1980–1981) postulates that CHF is triggered when the supply liquid flow normal to the wall is blocked by the vapor flow at a critical velocity. 4) The near-wall bubble crowding mechanism (1981–1985) considers the limited turbulent interchange between the bubble region and the core region. The CHF occurs when the bubble region void fraction exceeds the critical value at which an array of ellipsoidal bubbles can be maintained without significant contact between the bubbles, thus preventing access of the supply liquid. 5) The liquid sublayer dryout mechanism (1988–2000) hypothesizes that the onset of CHF is caused by the depletion of a thin superheated liquid sublayer underneath a vapor blanket flowing over the heated surface. The liquid sublayer dryout occurs when the liquid entering the sublayer falls short of balancing the sublayer evaporation. 6) The interfacial lift-off mechanism (1993) assumes that bubbles coalesce into a continuous wavy vapor layer. Wetting fronts occur periodically due to instabilities, supplying fresh 13 liquid to the heater. The CHF is triggered when vapor effusion lifts the first wetting front from the heated surface. Among the six model categories, the liquid sublayer dryout and near-wall bubble crowding have received considerable attention. The liquid sublayer dryout mechanism is supported by experimental evidence at pressures up to 2.4 MPa [10,11], including the existence of a very thin liquid layer trapped beneath the vapor blankets and observed fluctuant instabilities prior to DNB. Celata and Mariani [29] and Kandlikar [30] discussed the inadequacies in the other model categories. Some of their arguments—such as the fact that there would be no abrupt visible change in the macro-layer bulk flow pattern during the transition from nucleate boiling to film boiling—have been confirmed by recent experimental studies [9]. Furthermore, Katto [10] claimed that the DNB in subcooled or low-quality flow boiling has behavioral similarities to the pool boiling DNB, as both phenomena seem to be governed by the hydrodynamic instability of the liquid-vapor interface. It is therefore worth reviewing the state of art of the sublayer dryout-based models proposed by various authors. The original model of the liquid sublayer dryout mechanism was proposed by Lee and Mudawar (1988) [11] for subcooled flow boiling in a round tube, assuming that vertically distorted cylindrical vapor blankets are formed at near-CHF conditions as a consequence of small bubbles accumulating along the near-wall region. As illustrated in Fig. 3, a thin superheated liquid sublayer is postulated to be enclosed underneath the vapor blankets. The equivalent diameter (or thickness) of the vapor blanket is approximately equal to the departure diameter of individual bubbles under the respective T/H conditions and is expected to remain constant. The blanket length is assumed to be the critical wavelength of Helmholtz instability at the sublayer-blanket interface. The DNB occurs when the rate of liquid sublayer evaporation exceeds that of the liquid entering the sublayer, forming persistent dry patches on the heated surface. A force balance on the vapor blanket in the radial direction was applied to calculate the sublayer thickness, and three empirical constants were introduced in the process. The Lee model was later found to behave poorly at low pressures [9,31], as no relevant tests had been included in its calibration database. 14 Fig. 3. Generic sketch of the liquid sublayer dryout mechanism. Katto (1990) proposed a vapor stem approach [10] based on an earlier model by Haramura and Katto (1983) [32] for saturated pool boiling on an infinite horizontal flat plate. In this approach, the vapor blanket is continuously replenished by a large number of small vapor stems normal to the heated surface. In flow boiling, the same vapor stem coverage area used in pool boiling was assumed, and the effects of forced convection and subcooling were considered. No supply liquid was included during the time that the blanket passed over the sublayer. An empirical coefficient that related blanket velocity to the local velocity of a homogeneous two-phase flow was optimized with the model’s calibration dataset. The Katto model covered a pressure range of 0.1–20.0 MPa. However, it was found to be incapable of predicting accurate CHF at relatively high void fraction [31], and the existence of the presumed vapor stems was questionable due to lack of experimental support. Celata et al. (1994) [3] followed the same calculation sequence as that used in the Lee model for the vapor blanket thickness, length, and velocity. The Celata study included an effort to develop a rationalized model to reduce empiricism, in which the vapor blanket was postulated to lie on the superheated layer boundary: that is, the blanket existed only in the near-wall region where the local liquid was superheated. This strong assumption would eliminate all empiricism, but is deemed inapplicable to (near-)saturation conditions under which DNB may still occur, especially in PWR-like scenarios [33]. Furthermore, the model remained deficient in predictions at low channel length-to-diameter ratio or high-pressure conditions [31]. Celata et al. (1999) [34] later derived a so-called superheated layer vapor replenishment approach which used the same assumptions and closures as in their previous model. According to the new model, the CHF occurs when the vapor blanket replenishes the entire superheated layer, coming into contact with the heated surface. This 15 model conceptually resembles the bubble crowding-based models first proposed by Weisman and Pei [35]. The validation performance of the new Celata model was similar to that with the previous version, and the superheated layer boundary assumption was also defective at low subcooling. Liu et al. (2000) [31] further considered instabilities at the sublayer-blanket and the blanket-bulk interfaces. The Liu study proved that both interface wavelengths were equal to the vapor blanket length as long as the blanket remained stable, and their model marginally outperformed the 1994 Celata model on a tube validation dataset of over 2,400 test cases. However, in the Liu model, the void fraction in the postulated core region was simplified as the channel-averaged void fraction (i.e., all bubbles were assumed to appear in the bulk) since the liquid sublayer and the vapor blanket were expected to be very thin. This assumption contradicts experimental observations [10] and the bubble crowding mechanism in which most bubbles accumulate in the vicinity of the heated wall. Like all the aforementioned models, the Liu model only focused on subcooled boiling. In view of the limitations in existing studies, an improved macro-scale model is presented here. The proposed model is intended to leverage the key concepts in both mechanisms of sublayer dryout and bubble crowding while minimizing additional empiricism. 3.2. Model description 3.2.1. Governing equation According to Lee and Mudawar [11], the DNB is triggered at the depletion of the liquid sublayer. The near-CHF cylindrical-like vapor blanket touches the heated wall and forms a dry patch as a result of the Helmholtz instability. The dry patch persists and spreads over the wall if the enthalpy of the liquid entering the sublayer falls short of balancing the evaporation of the sublayer. As illustrated in Fig. 4, using a control volume that surrounds the sublayer and moves at the same velocity as the vapor blanket, the heat balance can be written as 𝑞/01))𝐿#𝐷#=𝐺23𝛿𝐷#[ℎ'(+𝑐"2(𝑇456−𝑇2)], (1) where 𝑞/01)) is the critical heat flux, 𝐿# is the length of the cylindrical-like vapor blanket, 𝐷# is the equivalent diameter of the blanket, 𝐺23 is the relative mass flux of the liquid entering the sublayer, 16 𝛿 is the sublayer thickness, ℎ'( is the latent heat of vaporization, 𝑐"2 is the liquid specific heat, 𝑇456 is the saturation temperature, and 𝑇2 is the temperature of the liquid entering the sublayer, for which the subcooling (in subcooled flow) can be approximated as 𝑇456−𝑇2=𝑐7(𝑇456−𝑇,829), (2) where 𝑇,829 denotes the channel local bulk temperature, and the empirical constant 𝑐7 = 0.35 was optimized in the Lee model, suggesting that the supply liquid subcooling would be equal to 35% of the bulk subcooling. The effect of 𝑐7 on the predicted CHF was evaluated through sensitivity analysis and is described in Section 3.3. Fig. 4. Schematic view of a near-CHF vapor blanket in the present model. Since the liquid velocity in the thin sublayer is small compared to the blanket velocity, the relative mass flux of the liquid entering the sublayer (downward) can be simplified as 𝐺23≈𝜌2𝑈#, (3) where 𝜌2 is the liquid density and 𝑈# is the blanket velocity. Combining Eqs. (1)–(3) and calculating both 𝜌2 and 𝑐"2 at saturation conditions, one can obtain: 𝑞/01))≈:!;<&!"+=#=$!(?%&'@?()*+)BC,𝑈#. (4) The key for the sublayer dryout mechanism turns to the determination of 𝐿#, 𝑈#, and 𝛿. 3.2.2. Calculation of vapor blanket length 17 As mathematically derived in Liu et al. [31], at the DNB location, the blanket length should be equal to the critical wavelength at interfaces I and II (see Fig. 4). That is, 𝐿#=𝜆7=𝜆D, (5) where 𝜆7 can be expressed as 𝜆7=DEF(:!+:"):!:"(G,@G%)(*&-./)0≈DEF(:!+:"):!:"G,0, (6) with 𝑈48,25H%3≪𝑈#. Similarly, the wavelength at interface II is written as 𝜆D=DEF(:0+:"):!:"(G0@G,)0=DEF(:(+:"):!:"(G(@G,)0, (7) where the subscript 2 represents the local region on the right side of interface II and the subscript b denotes bubble region, as sketched in Fig. 5. As proposed in the original bubble crowding model [35], just before DNB, the flow channel can be radially divided into two regions: the near-wall bubble region and the core region (Fig. 5). Based on the modeling of the bubble region thickness described in Section 3.2.3, the local region 2 shall fall within the bubble region since the sublayer thickness can be neglected as compared to the blanket thickness. From Eqs. (5)–(6), the calculation of 𝐿# requires knowledge of the vapor blanket velocity 𝑈#. Fig. 5. Conceptual view of the present model in a tube channel (adapted from Liu et al. [31]). 3.2.3. Calculation of vapor blanket velocity Combining Eqs. (5)–(7) produces 18 𝑈#=G(7+I1!1!21" 1(21"1(, (8) where the mixture density in the bubble region 𝜌, can be calculated from Eq. (9), assuming two-phase flow at thermal equilibrium (i.e., saturation): 𝜌,=𝜌(𝛼,+𝜌'(1−𝛼,). (9) The void fraction in the bubble region 𝛼,, also known as the critical wall-void fraction, is a key parameter inherent to the bubble crowding mechanism. While various constants or correlations have been proposed [35–40], there are not sufficient data [37] to accurately determine 𝛼, at this point. Kodama and Kataoka [38] suggested that the predicted CHF should not be strongly dependent on the critical wall-void fraction value as long as it stays within a reasonable range set in the literature. To avoid introducing extra empiricism and to remain consistent with prior work, its value is set to be 𝛼, = 0.8 in the present model. Another key parameter in the bubble crowding mechanism is the bubble region thickness 𝛿,. It has widely been modeled as a constant multiplier of the bubble departure diameter 𝐷$; i.e., 𝛿,=𝑐D∙𝐷$, where the multiplier 𝑐D varies from 1.0 up to 10.0 when being fitted with datasets at different operating conditions and with different bubble departure diameter models [35–37,41]. It has also been suggested that in the fully-developed flow of a tube channel, the bubble region thickness would increase with the channel heated length until a plateau is reached. A correlation for 𝑐D as a function of the channel length-to-diameter ratio was developed in this work.d That is, 𝛿,=𝑐D∙𝐷$, where 𝑐D=𝑚𝑖𝑛 (0.1𝐿/𝐷,10) and 𝐿/𝐷>10. (10) Bubble departure diameter 𝐷$ The determination of the bubble departure diameter 𝐷$ is essential not only for 𝛿, (Eq. 10), but also for the vapor blanket thickness (or equivalent diameter) 𝐷#, because 𝐷#≈𝐷$ is one of the basic assumptions in the sublayer dryout mechanism [11]. Its modeling has been the focus of numerous efforts, and various (semi-)empirical or mechanistic models have been proposed, among which the Staub model [42] has received special attention in the scope of mechanistic DNB d The values of 𝑐3 were evaluated as part of the sensitivity analysis as presented in Section 3.3. 19 modeling [3,34,38,40]. Based on the force balance on a hemispherical-shape bubble attached to the wall, the Staub model can be expressed as 𝐷$=KD'F=4:!L0, (11) where the friction factor 𝑓 is calculated with the Colebrook–White equation (in tube channels): 7M'=1.14−2.0𝑙𝑜𝑔7NaO!+P.KRS%M'b=1.14−2.0𝑙𝑜𝑔7N(=5!6!+P.KRS%M'). (12) The surface roughness 𝜀 is on the order of 0.75 times the bubble departure diameter [3], or 𝑐T = 0.75. The channel Reynolds number 𝑅𝑒 is calculated with saturated liquid properties. The empirical coefficient 𝑐K in Eq. (11), defined as 𝑐K=𝐹483'5=% 6%U4VWU/(𝜋𝐷$𝜎), was found to lie within the interval [0.015, 0.17] based on pool boiling experiments. Staub [42] also calibrated 𝑐K with limited measurements in flow boiling and suggested a value of 0.02–0.03. Following the recommendation, Celata et al. [3,34] adopted 𝑐K = 0.03, and Pan et al. [40] used a much smaller value of 𝑐K = 0.0075. Kodama and Kataoka [38] suggested that 𝑐K would vary widely with the local T/H conditions and introduced a linear function of local subcooling based on qualitative parametric trends. Recent efforts have focused on reassessing the acting forces on the bubble and proposed more generic formulae to calculate the bubble diameter at departure to bulk or by sliding, such as that presented in Mazzocco et al. [43]. However, such mechanistic models still require calibration to determine multiple closure terms, and have only been validated against low mass flux, high subcooling data. At this stage, very limited bubble diameter measurements are available under high pressure, high flow rate or near-CHF conditions. Moreover, in the proximity of CHF, the size of individual bubbles may become practically unmeasurable since they tend to quickly coalesce and form agglomerated vapor clusters. Hence, in the scope of DNB modeling for this work, the vapor blanket thickness might not correspond to the actual bubble departure diameter, or bubble diameter when leaving the wall; instead, it should rather be modeled as an “equivalent” diameter. Therefore, a modified Staub model has been adopted in this work where the coefficient 𝑐K is optimized with a tube dataset (see Section 3.3). Note that 𝑐K is the only tunable term in the CHF prediction procedure, and its calibrated values are well bounded within the interval suggested by previous studies. The resulting bubble diameter is subject to sensitivity analysis in Section 3.3. Bubble region velocity 𝑈, 20 Using Eqs. (10)–(12), the bubble region thickness 𝛿, can be calculated, and then the fraction of the cross-sectional area occupied by the core region 𝜂==𝐴=/𝐴 is obtained, where 𝐴 denotes the channel total cross-sectional area. The void fraction and mixture density in the core region, 𝛼= and 𝜌=, are calculated from 𝛼==X78Y7−X((7@Y7)Y7, and (13) 𝜌==𝜌(𝛼=+𝜌2,,829(1−𝛼=), (14) where the channel-averaged void fraction at the DNB location 𝛼=& can be obtained from a commonly used drift-flux approach such as the Chexal–Lellouche model, which requires the flow quality (𝑥=&) information. 𝑥=& is calculated either from Levy [44] or from the Saha and Zuber model [45]. Assuming homogeneous flow in the core region [37], its quality 𝑥= turns to 𝑥==:"[7":"[7"+:*,()*+[7*=:"X7:"X7+:*,()*+(7@X7)=:"X7:7. (15) The total liquid mass flow rate 𝑚̇2 is written as 𝑚̇2=𝑚̇(1−𝑥=&)=𝐺𝐴(1−𝑥=&). (16) To calculate the liquid mass flow rate in each region, one may start with the average liquid velocity in the bubble region, 𝑈d,2. This can be derived by assuming that liquid within the bubbly layer follows the Karman velocity distribution [37]:e 𝑈,2+=e𝑦+, 0≤𝑦+<55.0𝑙𝑛𝑦+−3.05, 5≤𝑦+<302.5𝑙𝑛𝑦++5.5, 𝑦+≥30, (17) where 𝑈,2+=G(*G:, 𝑈*=k*;:!, 𝑦+=HG::!\!, and 𝜏.='L0]:!. 𝑈* is the friction velocity, 𝜏. is the wall shear stress, 𝑦 here is a variable representing the distance from the heated wall, and 𝜇' is the dynamic viscosity of saturated liquid. The average velocity then becomes 𝑈d,2=𝑈d,2+∙𝑈*, with 𝑈d,2+ obtained from the integral of 𝑈,2+ with respect to 𝑦+ (from 0 to 𝑦;(+=;(G::!\!). Once 𝑈d,2 is known, the liquid mass flow rate in the bubble region 𝑚̇,2 can be written as e This assumption was used by all the previous sublayer dryout-based models. The Karman velocity profile was originally developed for single-phase turbulent flow; whether it is applicable to regions where significant bubbles are present remains debatable. Its effect on the predicted CHF was evaluated through sensitivity analysis (Section 3.3). 21 𝑚̇,2=𝜌'𝑈d,2𝐴,2=𝜌'𝑈d,2𝐴(1−𝜂=)(1−𝛼,), (18) and the liquid mass flow rate in the core region can be written as 𝑚̇=2=𝑚̇2−𝑚̇,2. Then, the total mass flow rate in the core region is 𝑚̇==𝜌=𝑈=𝐴𝜂=, (19) where 𝑈=, by assuming homogeneous flow, is 𝑈==𝑈=2=-̇7*:*,()*+[7*=-̇7*:*,()*+[Y7(7@X7), (20) and the total mass flow rate in the bubble region becomes 𝑚̇,=𝑚̇−𝑚̇==𝐺𝐴−𝑚̇=. The bubble region quality 𝑥, turns to 𝑥,=-̇(<-̇(=1−-̇(*-̇(. (21) Finally, the bubble region (average) velocity, 𝑈,, is expressed by 𝑈,=-̇(:([(7@Y7)=L(_78@_7):((_(@_7)(7@Y7). (22) The vapor blanket velocity (𝑈#) can then be obtained from Eq. (8). 3.2.4. Calculation of liquid sublayer thickness As illustrated in Fig. 4, the sublayer thickness 𝛿 is approximately the difference between the distance from the wall to the blanket centerline 𝑦 and one half of the blanket thickness 𝐷#; i.e., 𝛿≈𝑦−0.5𝐷#≈𝑦−0.5𝐷$. (23) Making use of Eq. (17), one can obtain 𝑦=\!G::!∙e𝑈#2+, 0≤𝑈#2+<5exp [(𝑈#2++3.05)/5], 5≤𝑈#2+<13.96exp [(𝑈#2+−5.5)/2.5], 𝑈#2+≥13.96, (24) where 𝑈#2+=𝑈#2/𝑈* and 𝑈#2 is the liquid velocity at the blanket centerline. 𝑈#2 is calculated via a force balance on the vapor blanket between buoyancy and drag force, and can be arranged as 𝑈#2=𝑈2q𝑦=𝛿+!,Dr=𝑈#−kDC,((:!@:"):!/=. (25) 22 The drag coefficient 𝐶! is obtained by either the Harmathy model (if pressure is smaller than 1 MPa) or the Chan and Prince expression (otherwise) [31]. Through the equations described above, the CHF can be predicted by an iterative procedure for a given set of geometric and flow conditions. All fluid properties are obtained at saturation conditions unless otherwise specified. An initially assumed CHF value is explicitly involved in the calculation of the channel flow quality and void fraction, propagating all the way to the determination of the three key parameters: vapor blanket length, vapor blanket velocity, and liquid sublayer thickness. Convergence is reached when both sides of Eq. (4) are numerically equal. 3.3. Model optimization and sensitivity analysis To optimize the proposed mechanistic model, an extensive DNB-specific tube dataset was collected. Tube heaters are used in the stage of model optimization because of their relatively high data availability and ease of modeling. As shown in Table 1, the tube dataset consists of 1,439 uniformly heated test cases from five data sources, and it covers a broad range of flow conditions. The test matrix is DNB-specific such that for each data point, the channel void fraction is less than 70%, or the local equilibrium quality does not exceed 0.2 [2]. Table 1. Parameter ranges of tube CHF data. Reference No. of data P [MPa] 𝐺 [kg/m2-s] 𝐷 [mm] 𝐿/𝐷 [-] 𝑇!"#,%& [°C] 𝑥',("),'*+ [-] Inasaka [46] 7 0.31 to 0.91 4,300 to 6,700 3.0 33 74 to 148 -0.15 to -0.04 Peskov [47] 17 10.0 to 20.0 750 to 5,361 10.0 40 to 165 5 to 205 -0.23 to 0.13 Thompson [48] 1,202 0.1 to 20.7 542 to 7,975 1.1 to 37.5 12 to 366 0.2 to 338 -0.45 to 0.21 Weatherhead [49] 162 13.8 332 to 2,712 7.7 to 11.1 41 to 59 2 to 299 -0.49 to 0.19 Williams [50] 51 5.5 to 15.2 670 to 4,684 9.5 193 21 to 251 -0.03 to 0.17 Total 1,439 0.1 to 20.7 332 to 7,975 1.1 to 37.5 12 to 366 0.2 to 338 -0.49 to 0.21 To calibrate the tunable coefficient 𝑐K for the bubble departure (or rather equivalent) diameter 𝐷$, as in Eq. (11), polynomial regression was applied. For each data point, given the value of the experimental CHF, the bisection method was used to find the target value of 𝑐K. The 10-fold cross-validation technique was employed to optimize the polynomial structure and to assess how the regression results would generalize to out-of-sample data. The tube dataset was randomly 23 partitioned into ten subsets of equal size. Of the ten subsets, a single subset was retained as the validation set, and the remaining nine subsets were used as training data. This process was repeated ten times, with each of the ten subsets used exactly once as the validation set. In the optimized regression model, 𝑐K was expressed as a function of channel length-to-diameter ratio and local equilibrium quality, as the other measurable variables were already captured in the Staub equation. Through the cross-validation process, the set of polynomial coefficients that yielded the lowest validation error could be found. The resulting polynomial regression model with the corresponding (optimized) coefficients was used to represent 𝑐K in its calibrated function form. The calibrated function was then implemented in the iterative procedure for making new CHF predictions. In Fig. 6, the present model is compared with the mechanistic Liu model [31] and the widely-used Groeneveld 2006 LUT [6], a normalized data bank of 8-mm-diameter tubes containing over 30,000 CHF test cases and covering the full range of conditions of practical interest. The LUT values are adjusted by multiplicative correction factors to account for a number of specific geometric conditions and different heating profiles [2,6,51]. One can see from Fig. 6 that the present model clearly outperforms the other two methods: over 90% of predicted tube CHF results fall within ± 10% uncertainty (vs. 78% of predicted data with LUT and 35% with Liu for all cases; 71% with LUT and 62% with Liu for subcooled boiling cases). The cross-validated results (and not the regression-fitted training solutions) are shown for the present model. Table 2 further confirms the superior predictive capability of the present model with respect to at least three out of the five data sources, and the relative root-mean-square error (rRMSE) for all tube data is 0.07, vs. 0.11 with LUT and 0.27 with Liu. Unsurprisingly, the Liu model becomes much more accurate in subcooled boiling, yet the present model still outperforms. 24 (a) (b) Fig. 6. Cumulative data fraction of absolute CHF relative errorf with LUT, Liu, and the present model on (a) all tube data and (b) subcooled boiling tube data. Table 2. DNB model comparison on tube data (all / subcooled cases). Data source No. of data LUT rRMSE Liu rRMSE Present rRMSE Inasaka 7 / 7 0.31 / 0.31 0.18 / 0.18 0.23 / 0.23 Peskov 17 / 11 0.10 / 0.09 0.14 / 0.05 0.09 / 0.06 Thompson 1,202 / 544 0.12 / 0.14 0.28 / 0.16 0.07 / 0.09 Weatherhead 162 / 89 0.05 / 0.05 0.18 / 0.13 0.06 / 0.05 Williams 51 / 5 0.04 / 0.05 0.21 / 0.03 0.03 / 0.04 Total 1,439 / 656 0.11 / 0.13 0.27 / 0.16 0.07 / 0.08 Sensitivity analysis (SA) is the study of how different input variations of a model influence the variability of its output. There are two main classes of SA methods based on local or global definitions. The local SA analyzes the effect of one parameter on the model output at a time while keeping the other parameters fixed, whereas the global SA aims at quantifying the relative significance of uncertain input parameters over the entire variation range. The latter approach is better suited for nonlinear input–output relationships, and it is more realistic since it allows all f Absolute CHF relative error = absolute value of CHF relative error, defined as (𝐶𝐻𝐹>?@ABCD@A−𝐶𝐻𝐹@E>)/𝐶𝐻𝐹@E>. 25 targeted input factors to be varied simultaneously. In this work, the regression-based global SA of the present DNB model was performed on the key closure terms to assess their associations with the model output (i.e., CHF value) and to evaluate how variation in the output can be apportioned to different sources of input uncertainty. To measure the effect of each input parameter on the output, the Pearson correlation coefficient (PCC) and the Spearman rank correlation coefficient (SRCC) have been used as sensitivity measures in nuclear engineering applications [52,53]. The PCC between two variables is their covariance divided by the product of their standard deviations; the SRCC is defined as the PCC between two ranked variables. Both coefficients can take a range of values from +1.0 to -1.0. As shown in Table 3, seven closure terms (or their leading coefficients) were selected as the input parameters. For each case in the tube dataset, all seven parameters were simultaneously sampled 1,000 times by assuming that their perturbations followed normal distributions. The Latin hypercube sampling method was used for its ability to generate controlled random samples and to make sampling point distribution close to the true distribution. The input uncertainty bounds were determined using expert judgement with appropriate conservatism—based on discussions in the literature and as described in Section 3.2—and they corresponded to three standard deviations of the postulated normal distributions, with the nominal values as the distribution means. Table 3. Selected closure terms for model SA. Parameter Description Uncertainty bound Eq(s). in Section 3.2 𝐷! bubble departure diameter ± 100 % 10, 11, 12, 23 𝛼" bubble region void fraction ± 20 % 9, 13, 18 𝑐# supply liquid subcooling constant ± 100 % 4 𝑐$ bubble region thickness coefficient ± 100 % 10 𝑐% surface roughness fraction ± 100 % 12 𝑈"& bubble region liquid velocity ± 30 % 17, 18 𝑦 distance from wall to blanket centerline ± 30 % 23, 24 To illustrate the correlation coefficients, two example cases were selected from the tube dataset, representing a PWR-like scenario (Case No.1) and a low pressure, high mass flux (Case No.2) scenario. Table 4 summarizes their nominal input/output values. The correlation coefficients are compared in Fig. 7 with respect to each of the seven closures. The length of the horizontal bars 26 indicates the sensitivity of the respective input parameter on the predicted CHF. A larger absolute value (closer to 1.0) indicates that the uncertain parameter is more significant to the output. The positive sign suggests that the model output tends to move in the same direction as the corresponding input, and vice versa. For these two cases, it can be observed from Fig. 7 that (1) the PCC and the SRCC values agree closely, and (2) the model output tends to move in the same direction as four of the input parameters (𝐷$, 𝑦, 𝑐D, 𝑈,2) and in the opposite direction as two others (𝛼,, 𝑐T). The suggested correlation directions are physically reasonable, and the predicted CHF is insensitive to the supply liquid subcooling constant 𝑐7. Table 4. Key closure sensitivity: nominal values for the selected tube test cases. Case No. 𝐷A [mm] 𝛼F [-] 𝑐G [-] 𝑐3 [-] 𝑐H [-] 𝑈FI [m/s] 𝑦 [mm] Predicted CHF [MW/m2] 1 2.18e-2 0.8 0.35 10.0 0.75 4.30 1.10e-2 2.53 2 7.56e-3 0.8 0.35 3.3 0.75 3.96 5.50e-3 15.36 (a) (b) Fig. 7. Correlation coefficients with respect to seven key closure terms in the present DNB model for (a) Case No.1 and (b) Case No.2. In Fig. 8, the output results of all 1,000 samples are plotted as a normalized histogram—such that the sum of all bins is one—for both test cases. The histogram profile can be closely fitted by 27 a normal probability density function with almost the same mean and standard deviation as the sample statistics. In Case No.1, large input uncertainties which were hypothetically bounded at up to ± 100% led to less than 6% deviation from the nominal CHF for 997 out of 1,000 samples. In Case No.2, the simulated outputs were more spread out, yet all samples fell within 13% of the nominal value. As for the entire tube dataset, such large closure perturbations have led to on average less than 9% deviation from the nominal CHF values. The output deviations can be mainly apportioned to variations in the distance from the heated wall to the blanket centerline (𝑦) and the bubble departure diameter (𝐷$). Given the conservatism in the input perturbations, the relatively small spread indicates a convincing stability and robustness of the proposed mechanistic model. (a) (b) Fig. 8. Key closure sensitivity: normalized histogram for (a) Case No.1 and (b) Case No.2. 4. PIMLAF performance in rod bundles In this section, the freshly proposed mechanistic DNB model is integrated into the hybrid framework PIMLAF as one candidate for the physics-informed component. In the scope of hybrid modeling, the proposed mechanistic model is referred to as the Zhao model. The resulting hybrid models are assessed against their standalone DK and ML counterparts for both interpolation and extrapolation purposes through two bundle-related cases studies. As summarized in Table 5, the bundle dataset consists of 218 five-by-five rod bundle DNB test cases from the PSBT benchmark [28], in which the heating profiles and the spacer grid configurations differed among five test series. 28 All the assemblies contained 25 heater rods except for assembly type A8, which had a central thimble rod (see Fig. 9). Three types of spacers were instrumented along the axial length: simple spacers, spacers with no mixing vanes, and spacers with mixing vanes. More details can be found in the benchmark specifications report [28]. Within the selected dataset, the outlet pressure ranged from 7.5 to 16.8 MPa, the inlet mass flux from 1,367 to 4,817 kg/m2-s, and the inlet temperature subcooling from 17 to 165 K. While experimental uncertainties were not reported, other studies [54][55] suggest that the CHF measurement error would be on the order of 5–20% on average for a test case with uniform axial power profile.g Table 5. PSBT DNB benchmark bundle test matrix. Test series Assembly Axial power shape No. of spacers Grid spacing [mm] No. of data TS0 A0 uniform 13 305 61 TS2 A2 uniform 17 229 33 TS4 A4 cosine 16 227; 254 36 TS8 A8 cosine 16 227; 254 43 TS13 A4 cosine 16 254 45 Fig. 9. PSBT DNB benchmark bundle cross-sectional view. 4.1. Interpolation: the PSBT bundle example In order to determine local T/H information in a rod bundle where the flow area is modeled as a multitude of interconnected subchannels, the well-recognized and validated subchannel code CTF [56–58] was used in this work. Relevant subchannel-level boundary conditions were provided by CTF as inputs for ML and prior models. To save computational time, only the peak (i.e., non-peripheral) subchannel candidates in the selected PSBT dataset were considered. All the peak g Local CHF in a flow channel with non-uniform axial power profile may have a higher measurement error. 29 subchannels shared the same geometry: a typical channel-centered cross-section with a 9.5 mm rod outer diameter and a 12.6 mm pitch. Therefore, the channel hydraulic and heated diameters were not included in the input feature vector. The ML input vectorh included the following features: - local pressure (7.6–16.8 MPa) - local mass flux (1,409–4,958 kg/m2-s) - local equilibrium quality (-0.12–0.24) - heated length (2.69–3.58 m) - grid spacing (227–305 mm) - distance between DNB location and closest upstream spacer (151–227 mm) - number of upstream mixing vane spacers (5–7) As summarized in Table 6, the DK-based prior models are the LUT, the W-3 correlation, and the mechanistic Zhao model. Both LUT and W-3 are available in CTF, so they can be called online within the subchannel code. For each case, the predicted CHF value (using LUT or W-3) was found via subchannel analysis by iterations in CTF, and several correction multipliers were included to account for effects such as heat flux nonuniformity and spacer grids. The Zhao model was validated with internally heated annulus data and the PSBT bundle data to demonstrate its ability to make reliable predictions without recalibration or correction factors for non-tube heater geometries under a variety of operating conditions over the database considered in this study. Its validation results are summarized in Appendix A. The ML methods are NN and RF, both trained and tested using the high-level Keras application programming interface with TensorFlow backend and the scikit-learn library in Python 3.6. While the RF generally performs well with raw inputs, feature scaling is highly recommended for NNs. Therefore, standardization (which normalizes an input to have zero mean and unity variance) was applied to all features except the local equilibrium quality. Based on hyperparameter tuning with 10-fold cross-validation, the best-estimate ML structure could be simplified within the hybrid framework: the NN had a single hidden layer of eight units (i.e., a 7/8/1 network architecture vs. 7/20/1 if standalone), and only ten estimators (instead of 30) were needed in the RF. h ML input features are usually selected using expert judgement based on prior knowledge, although more rigorous feature selection techniques may be needed when dealing with strongly correlated high-dimensional inputs. 30 Table 6. Hybrid framework configuration for PSBT bundle data. Prior model LUT; W-3; Zhao ML-aided component NN; RF Best-estimate NN model • network architecture 7/8/1 (7/20/1 if standalone NN) • optimization algorithm Adami • activation function (hidden layer) rectified linear unit (ReLU)j Best-estimate RF model • no. of trees/estimators 10 (30 if standalone RF) • percentage of features at each split 40–60% Figure 10 compares the performance of the standalone predictive tools (ML models: NN and RF; DK models: LUT [Fig. 10a], W-3 [Fig. 10b] and Zhao [Fig. 10c]) with that of the hybrid models (NN+LUT and RF+LUT [Fig. 10a], NN+W-3 and RF+W-3 [Fig. 10b], NN+Zhao and RF+Zhao [Fig. 10c]). All ML (standalone and hybrid) results were taken from cross-validation and averaged across 100 repeated training–validation processes to account for the inherent randomness in ML. The hybrid approach is perfectly in line with measurements and clearly outperforms its standalone counterparts: the rRMSE values are smaller at less than 4%, regardless of the prior/ML model, whereas the best standalone model for this dataset (i.e., Zhao) has an rRMSE of 6.5%. The corresponding cumulative data fraction curves rise faster than those with standalone models. Within the PIMLAF, the prior model is fixed throughout the training and validation stage, and the ML-aided component can guarantee (for interpolation) that the final overall prediction performance is at least as good as the one obtained with the prior model only,k regardless of the ML model’s complexity. Additionally, such a framework can leverage admitted laws of physics and well-established empirical relations in the prior model to help inform ML. It is expected to outperform standalone ML when the input/output causality is not fully established by the attempted ML structure alone or by the data sources at hand, or when the prior model is at least as accurate as standalone ML. One may argue that with an extensive high-quality dataset and i Adam [70] is one of the most popular gradient-based optimization algorithms and is based on adaptive estimation of the first- and second-order moments. j ReLU is one of the most commonly used activation functions in ANNs. It is a universal approximator, as any function can be approximated with combinations of ReLU. k This conclusion does not necessarily hold for individual data points because the ML-aided component is expected to minimize the loss function, an optimization objective with respect to the entire training set. 31 an ideally configured structure, a standalone ML model, especially an NN, would theoretically be able to emulate similar performance. However, such efforts are deemed cumbersome and would introduce an extra layer of complexity. (a) (b) (c) Fig. 10. Cumulative data fraction of absolute CHF relative error with standalone and hybrid models on PSBT bundle data: (a) LUT, (b) W-3, and (c) Zhao as DK model. 32 Another observation based on the results presented in Fig. 10 is that while the standalone LUT and W-3 perform worse than the Zhao model, their hybrid versions—combined with either NN or RF—actually outperform marginally. This paradoxical finding is consistent with the previous assessment on tube data [27] and implies that the ML-aided component in the PIMLAF is likely to play a greater role when dealing with relatively more scattered or biased residuals, when the prior model is less accurate. A reasonable interpretation is that randomness may prevail over undiscovered regular trends in the residual matrix if the prior model is very accurate itself [22], making it challenging for ML to capture the “true” information. This finding further suggests, in the scope of interpolation, that the quest for a perfectly developed DK-based model would be deemed unnecessary from the practical engineering viewpoint. 4.2. Implication for MDNBR margin reduction The enhanced predictive capabilities of the hybrid framework are further confirmed by the results presented in Fig. 11, which compares the hybrid RF+Zhao model against the standalone W-3/Zhao/RF in the form of boxplots. In a boxplot, the box displays the lower (or first) and upper (or third) quartile; the horizontal band inside the box is the median (or second quartile). The end of the lower whisker represents the lowest datum within 1.5 times the interquartile range (IQR) of the lower quartile, and that of the upper whisker represents the highest datum within 1.5 times IQR of the upper quartile. The outliers are plotted as individual points beyond the two whiskers. It can be observed that the hybrid model yields the smallest scatter and the closest-to-zero medians in terms of relative error. All the predicted CHF values with RF+Zhao fall within ± 15% error, vs. the biased W-3 or the largely scattered RF solutions.l With the hybrid model, no performance difference is observed between uniform (TS0, TS2) and cosine (TS4, TS8, TS13) axial power shape cases. l The large scatter in standalone ML models may be due to the presence of outliers in the training/validation set or the measurement uncertainty for each individual data point. In addition, the training and validation/test data should follow the same distribution, a prerequisite to statistically meaningful out-of-sample solutions in ML. 33 Fig. 11. Boxplots of PSBT CHF relative error with W-3, Zhao, RF, and hybrid RF+Zhao. While the present dataset is limited, such low scatter in predicted values with the PIMLAF could potentially translate into a less conservative MDNBR margin (or limit) for model uncertainties and therefore offer options for reactor power upgrade. The design criterion established for DNB—the so-called 95/95 tolerance rule—is that there is at least 95 percent probability that the DNB will not occur on the limiting fuel rod during normal operation and anticipated transients with 95 percent confidence. To meet this criterion, for a typical Westinghouse PWR, the MDNBR limit was set at 1.17 using the proprietary WRB-1 correlation or at 1.3 using W-3 [59].m These values correspond to the one-over-MDNBR(-limit) slope line on a measured vs. predicted CHF plot, and at least 95% of the plotted points shall fall above this line. As for the PSBT bundle with the present dataset, when using a standalone model, the MDNBR limit is equal to 1.19 with W-3, 1.16 with LUT, 1.14 with Zhao, 1.21 with NN, and 1.23 with RF; when using the hybrid RF+Zhao, the MDNBR margin is reduced to 1.07. The other combinations of prior model and ML yield very similar MDNBR margin values. To potentially achieve an m According to email correspondence with a senior project manager at Framatome, the U.S. vendors currently adopt a limit between 1.12 and 1.2 using different correlations for DNB, which is designed to cover the full range of PWR operating conditions. 34 MDNBR margin reduction for industrial applications, reactor vendors are encouraged to leverage their in-house experimental efforts and apply the PIMLAF on the designs of interest. 4.3. Extrapolation: from tube to bundle In this section, the extrapolation performance of the PIMLAF is explored by assuming that the PSBT bundle experimental CHF values were not available during training (i.e., they were treated as test data only). Only the tube data from Table 1 were accessible for training. In the hybrid framework, the prior model is LUT or Zhao, both of which can predict tube and bundle CHF for this application relatively well. The ML-aided component is NN or RF. In both standalone and hybrid versions, for a proper learning procedure, the ML training and test data must have the same number of input features, which was six in this case: local pressure, local mass flux, local equilibrium quality, channel hydraulic diameter, channel heated diameter, and heated length. Note that the spacer related information was not retained because those bundle-specific features did not appear in tubes. For this example, since the tube dataset covers a wider range of operating conditions than the PSBT benchmark (except for the heated length, which is a relatively weak feature), extrapolation is assessed in terms of channel geometry: from tube channel to bundle subchannel. The standalone DK models LUT and Zhao serve as reference solutions since they do not extrapolate here. The best-estimate ML models are secured through hyperparameter tuning with tube data. The train–test process is first run with standalone NN and RF. As can be seen from Table 7, the standalone ML extrapolates poorly, as less than 65% of predicted bundle data fall within ± 20% error. As mentioned in Zhao et al. [27], classical ML is only designed for and capable of interpolating data, but it is incapable of extrapolating. It learns to fit the known situations (no matter if data are used for training, validation, or testing) locally as closely as possible, regardless of how it performs outside these situations [60]. The hybrid approach results in a much more accurate prediction for test data (i.e., PSBT bundle). It can be concluded from Table 7 that the PIMLAF is guaranteed to be more extrapolative than standalone ML, and the margin of performance improvement is directly related to the prior model robustness and accuracy. In other words, the choice of prior model is much more important 35 for extrapolation than for interpolation, especially when the training dataset is within a limited range. A robust prior model provides additional reliable input–output relations that help inform the ML-aided component. Table 7. Performance of different CHF models on PSBT bundle data: extrapolation. Model Test data rRMSE Test data within 10% error Test data within 20% error LUT (reference) 0.09 67% 100% Zhao (reference) 0.07 86% 100% Standalone NN 0.22 30% 64% Standalone RF 0.33 23% 47% Hybrid NN+LUT 0.08 74% 100% Hybrid RF+LUT 0.08 73% 100% Hybrid NN+Zhao 0.07 87% 100% Hybrid RF+Zhao 0.06 89% 100% Finally, it is worthwhile to mention that the PIMLAF slightly outperforms standalone DK even though their extrapolation capabilities cannot be directly compared here. The role of the ML-aided component seems to be significantly more crucial in interpolation (see Fig. 10) than in extrapolation, since ML is unlikely to be able to capture the correct and complete bundle input/output causality simply from limited tube data, not to mention that the unique features for bundles (e.g., spacer grids) were excluded from the process in this case study. 5. Conclusions and recommendations This paper extends prior work on the prediction of DNB and achieves superior predictive capabilities of CHF models for the PSBT database through a physics-informed ML-aided framework (PIMLAF). This first-of-a-kind hybrid framework relies on established DK to provide credible baseline solutions and uses ML to capture undiscovered information from the mismatch between the actual and DK-predicted output. To enhance the DK robustness, an improved physics-driven model is proposed in the first place, incorporating key concepts and assumptions in the liquid sublayer dryout and near-wall bubble crowding mechanisms. The proposed mechanistic model is then integrated into the hybrid framework as one candidate for the physics-informed component. Two bundle-related case studies have been carried out to illustrate the PIMLAF’s 36 improved performance over its standalone DK or ML counterparts for both interpolation and extrapolation purposes. In light of PIMLAF’s promising potential to reduce prediction error and scatter, reactor vendors are encouraged to leverage their in-house experimental efforts and apply the hybrid framework to potentially achieve MDNBR margin reductions for the designs of interest. Despite its demonstrated accuracy and robustness, the proposed mechanistic model was calibrated with tube CHF data, which inevitably diluted its mechanistic foundation. In the future, instead of tuning the model directly with specific datasets, insights may be sought from the higher resolution micro-scale DNB models based on wall heat partitioning [1,12] or from the freshly presented percolative scale-free model [61] which stochastically described near-wall bubble interactions. Those models rely on different triggering mechanisms from those assumed for macro-scale models, and their robustness may still be challenged since they rely strongly on high-fidelity data and require closure tuning. However, they provide additional theoretical considerations that could be potentially leveraged, such as heater-specific closures to include surface characteristics and bubble interaction dynamics for a better description of the coalesced bubble scenario at near-CHF conditions. Moreover, the predicting accuracy of CHF for rod bundles should also account for location accuracy (and not just the magnitude of CHF). This is not an issue for uniform axial heating profile cases without mixing vanes at the top of the bundle where DNB always occurs at the channel outlet, but in other cases—such as cosine power shape—the DNB will occur at different axial locations. While assessing CHF location accuracy is beyond the scope of this workn (mainly due to lack of experimental results), it will be important to validate the proposed model against location data once they become available and to quantify the associated uncertainties (which can be very large for measurements using thermocouples). Future work will also couple the proposed mechanistic model with recent progress in transient CHF modeling [62] to help improve T/H and fuel performance capabilities under fast transient scenarios such as during a reactivity-initiated accident (RIA). A desired next step for the PIMLAF is to demonstrate its industrial usefulness in a production environment with a code like CTF. The implementation may be achieved through a look-up table n In this work, the CHF axial locations for the PSBT bundle dataset were obtained from predicted location results using W-3 (in CTF), the most widely used correlation for evaluating DNB in PWRs. 37 format or an external wrapper. The data complexity of the hybrid framework can be further extended on the fly if more data sources and CHF-related features become available, including specific spacer grid geometries and heater surface characteristics. Similarly, the model complexity of the ML-aided component will be extended if there is interest for uncertainty quantifications using Bayesian inference or more advanced extrapolation using transfer learning. Due to the purely data-driven nature of ML, the PIMLAF still depends on data quality and does not guarantee physically realistic parametric trends for all individual test points when the data at hand are noisy or fluctuant. Adding a pre-processing stage of DK-informed data smoothing or boosting with synthetic data could potentially fix this issue, and the PIMLAF’s full predictability will be further demonstrated. Finally, the generic framework proposed in this work is suitable for a broad spectrum of applications for nuclear engineering problems in which domain knowledge and experimental (or high-fidelity numerical) data are available. Acknowledgments The first author would like to acknowledge the insightful discussions with Xu Wu, Assistant Professor at North Carolina State University, on the model sensitivity analysis work. This research was supported by the Consortium for Advanced Simulation of Light Water Reactors (www.casl.gov), an Energy Innovation Hub for modeling and simulation of nuclear reactors under U.S. Department of Energy Contract No. DE-AC05-00OR22725. Appendix A. Validation of the improved mechanistic DNB model A.1. Annulus data The mechanistic DNB model proposed in Section 3 is herein validated against a dataset containing 378 internally heated annulus test cases from three data sources (see Table A.1). Like the tubes (Table 1), the annulus dataset also covers a broad range of flow conditions, and the axial power shape is also uniform. Note that the tube diameter (𝐷) is substituted by the hydraulic diameter (𝐷%) and heated diameter (𝐷&) of the annulus. Table A.1. Parameter ranges of annulus CHF data. 38 Reference No. of data P [MPa] 𝐺 [kg/m2-s] 𝐷' [mm] 𝐷, [mm] 𝐿/𝐷, [-] 𝑇!"#,%& [°C] 𝑥',("),'*+ [-] Beus [63] 77 5.5 to 15.5 671 to 3,721 5.6 15.2 140 20 to 243 -0.32 to 0.20 Janssen [64] 282 4.1 to 9.7 381 to 5,913 4.6 to 22.2 11.3 to 96.3 18 to 159 3 to 212 -0.12 to 0.21 Mortimore [65] 19 8.3 to 13.8 677 to 3,637 5.0 13.3 161 20 to 242 -0.13 to 0.19 Total 378 4.1 to 15.5 381 to 5,913 4.6 to 22.2 11.3 to 96.3 18 to 161 3 to 243 -0.32 to 0.21 While the governing equation (Section 3.2.1) of the present model is insensitive to flow channel geometry, several closure terms require modifications when dealing with non-tubes, such as those for calculating the channel-averaged flow quality and void fraction. At this stage, there is no well-used annulus-specific closure model available in the open literature for calculating flow quality or void fraction. A best-practice recommendation is to adapt tube correlations for annulus geometry. The core region cross-sectional area also needs to be adjusted. No recalibration is conducted. Figure A.1 compares the present model with LUT in the form of boxplots for all three annulus data sources. On average, the present model yields closer agreement with experiments and smaller scatter as compared to LUT: closer-to-zero medians (except for the Beus data) and fewer outliers. As summarized in Table A.2, the rRMSE values are also marginally in favor of the present model. Fig. A.1. Boxplots of annulus CHF relative error with LUT and present model. Table A.2. DNB model comparison on annulus data. 39 Data source No. of data LUT rRMSE Present rRMSE Beus 77 0.15 0.13 Janssen 282 0.17 0.14 Mortimore 19 0.19 0.13 Total 378 0.16 0.14 A.2. PSBT bundle data For the selected PSBT bundle test series (see Table 5), the present model is compared with LUT and the W-3 correlation. The subchannel-level information is provided by CTF. Like most other mechanistic models, the present model requires a number of iterations to achieve solution convergence. A fully-coupled procedure with CTF would easily become computationally too intensive, because each iteration (i.e., each attempted heat flux value) requires a new run of the subchannel code. Instead, a semi-coupled approach is adapted. Such a simplified iterative procedure is based on a commonly accepted assumption: the local flow distribution computed by CTF remains unchanged at different heat flux levels close to the experimental CHF [66,67]. The core region cross-sectional area is adjusted by assuming that the heater rods surrounding the peak subchannels are covered with bubbles of the same bubble region thickness at near-CHF conditions. The subchannel hydraulic diameter is used in the corresponding closure relations. No recalibration is performed on bundle data. As can be observed in Fig. A.2 and Table A.3, the LUT and W-3 results are closely in line with each other. While all cases fall within ± 20% prediction error using either of the three approaches, the present model yields the smallest scatter without introducing any correction factor. It consistently provides unbiased predictions, although being slightly higher than measurements, for test series with different heating profiles and spacer configurations, as further demonstrated in Fig. A.3 by the CHF relative error vs. (a) outlet pressure, (b) inlet mass flux, (c) local equilibrium quality, and (d) inlet temperature subcooling. For this dataset, both LUT and W-3 tend to overpredict at relatively low quality (𝑥% < 0) and underpredict at relatively high quality (𝑥% > 0.1); the reverse tendency holds for the inlet subcooling dependence. 40 Fig. A.2. Boxplots of PSBT CHF relative error with LUT, W-3, and present model. Table A.3. DNB model comparison on bundle data. PSBT test series No. of data LUT rRMSE W-3 rRMSE Present rRMSE TS0 61 0.11 0.12 0.08 TS2 33 0.05 0.06 0.07 TS4 36 0.07 0.08 0.08 TS8 43 0.09 0.10 0.07 TS13 45 0.09 0.08 0.05 Total 218 0.09 0.09 0.07 41 (a) (b) (c) (d) Fig. A.3. PSBT CHF relative error with LUT, W-3, and present model vs. (a) exit pressure, (b) inlet mass flux, (c) local equilibrium quality, and (d) inlet subcooling. Regarding the impact of mixing vane spacer grids on the CHF, Weisman et al. [66] suggested that the prediction procedure of a mechanistic DNB model would remain unchanged except for the enhanced flow mixing. Lin et al. [67] argued that the effect of mixing vanes might not be covered by the turbulent mixing coefficient alone: in their mechanistic DNB model, doubling the coefficient increased the predicted CHF by less than 2%. According to Arment [68], another effect 42 of mixing vanes is the local suppression of void fraction due to the swirl flow-induced bubble size change. Physical modeling of these local parameters is specific to spacer design and may require assistance from high-fidelity CFD. A preliminary sensitivity study indicates that the predicted CHF with the present model is increased by ~2.1% on average for the selected PSBT dataset following a void fraction decrease of 5% at the predicted DNB location (upstream of a mixing vane, where the enhancement effect is expected to be marginal [69]).