Numerical model to estimate subcooled flow boiling heat flux and to indicate vapor bubble interaction

There are numerous technical applications where hot components, with uneven temperature distribution, require cooling. In such applications, it is desired to provide efficient local cooling of the hot spots, while avoiding unnecessary over-cooling of the other regions. Such an approach, known as precision cooling, has several advantages. In addition to the fact that it reduces the effort for cooling, it limits the unintended heat lost to the cooling medium. In liquid cooled systems, such as Internal Combustion Engines (ICE), subcooled flow boiling offers immense potential for precision cooling. The primary challenges in extracting this potential are understanding the complexities in the subcooled flow boiling phenomenon and estimating the risk of encountering film boiling. The present study introduces a numerical model to estimate the wall heat flux in subcooled flow boiling and the model includes a mechanistic formulation to account for vapor bubble interaction. The formulation for vapor bubble interaction serves two purposes: (a) blends two well-established models in the literature, one in the isolated bubbles regime and other in the fully developed boiling regime, to estimate the wall heat flux; and (b) provides information to limit boiling in order to not encounter film boiling. The results from the new model are validated with two different experiments in the literature and the wall heat flux estimated by the model is in agreement with experimental results and responsive to different input parameters, such as bulk velocity, operating pressure and inlet subcooling. The new model requires only input of local flow quantities and hence implementation in Computational Fluid Dynamics (CFD) is straightforward. © 2021 Elsevier Ltd. All rights reserved.


Introduction
Technical applications involving hot components often require cooling for maintaining structural integrity and for optimum performance. Such components could experience uneven distribution of temperature. Therefore the requirement for cooling varies spatially along the surface of the component. Precision cooling, i.e. intense local cooling of hot spots, while avoiding unnecessary overcooling of other regions is an efficient way of saving energy and limiting the costly effort s f or cooling. There are a number of methods used to achieve precision cooling, such as impinging jet cooling, boosted forced convection cooling and surface structure modifications. If the cooling medium is a liquid coolant, the latent heat in the coolant provides an additional potential for local cooling. The subcooled flow of a liquid coolant could result in occurrence of local boiling on the hot spots, which in turn would remove a substantial amount of heat. However, such local boiling needs * Corresponding author.
to be controlled and limited in order to avoid film boiling which might cease the heat transfer process and result in undesired consequences. One of the technical applications, which might benefit from subcooled boiling for precision cooling purpose, can be found in the coolant jacket of the internal combustion engine (ICE).
The potential in forced convection boiling or flow boiling for precision cooling in ICE is discussed by Campbell in his PhD thesis [1] . Campbell distinguishes between the 'evaporative boiling engine', which involves coolant operating at saturation temperature resulting in significant vapor generation, and the 'nucleate boiling engine', in which the coolant operates at subcooled conditions. He further clarifies, the nucleate boiling engine promotes local boiling and avoids the problem of vapor generation. He also summarises the work of previous researchers who recommend the potential in nucleate boiling based cooling systems with warning against film boiling. More recently, the need for such a cooling system in high power density ICEs is discussed by Steiner [2] . Steiner emphasises that, significant improvement in heat transferred to the coolant is achieved for acceptable increase in wall temperature due to the https://doi.org/10.1016/j.ijheatmasstransfer.2021.121038 0017-9310/© 2021 Elsevier Ltd. All rights reserved. occurrence of nucleate boiling locally in vicinity of hotspots. In summary, flow boiling improves the performance of cooling systems compared to conventional forced convection cooling, but occurrence of film boiling impedes heat transfer and can be detrimental to the component. More over, boiling is a complicated phenomenon involving presence and interaction of multiple phases. The various mechanical and thermal interactions between the liquid coolant, vapor bubbles, and the solid heater are explained by Shoji [3] and Steiner [4] . An in-depth understanding of the underlying physics based on these interactions is essential to extract the potential in flow boiling for local or precision cooling.
The phenomenon of boiling has been analyzed both experimentally [1,[5][6][7][8][9] , to name a few, and numerically [10][11][12][13][14][15] , to name a few. These experiments made in simple channels, often with visual access, help understand the physics and provide useful data for validation of numerical models. Boiling is so complicated that a comprehensive numerical model, that accounts for all the interactions and physical mechanisms involved, does not currently exist and would be very challenging to develop. Often, numerical models are developed based on one or a few dominant physical mechanisms prevalent under certain flow and thermal conditions.
The Boiling Departure Lift-off (BDL) model by Steiner et al. [15] is a well established model for subcooled flow boiling conditions. The model is based on force balance on an individual vapor bubble. The heat flux estimated by the model is in good agreement with test data when the boiling involves individual vapor bubbles, which get convected by the flow after nucleation. However, with increased boiling, the bubble population increases and the vapor bubbles start interacting with each other. Under such conditions, the model predictions significantly deviate from the test results [15] . With further increase in boiling intensity, the bubble nucleation and interactions are the only dominant mechanisms of heat transfer. The heat transfer is then not influenced by forced convection. Under such conditions, the heat flux can be estimated using a pool boiling correlation, such as the one proposed by Rohsenow [11] . Clearly, the increase in bubble population and bubble interactions is the gap to be bridged between the conditions of applicability of the BDL model and Rohsenow's correlation.
In the present article, a numerical boiling model is proposed that includes a formulation for bubble interaction. A new blended model bridges over the predictions of the BDL model and that of Rohsenow's correlation based on this formulation. In addition, the formulation for bubble interaction provides an indication to limit boiling in practical applications, in order not to encounter film boiling. Furthermore, a modification is suggested to the BDL model that improves the accuracy of the predicted results. The results from the new blended model are validated with experiments by Steiner et al. [15] and Lee and O'Niell [8] . The results from the new model are in good agreement with the experiments.

Theory
Flow boiling can be divided into subcooled or saturated, depending on whether the bulk temperature of the liquid, T bulk , is lower than or equal to its saturation temperature, T sat . In applications where significant amount of vapor generation due to boiling is to be avoided, subcooled boiling is of more relevance. In such cases, the colder liquid bulk ensures condensation of vapour bubbles formed on the hot surface and thereby prevents significant increase in vapour fraction of the system. Different boiling regimes are identified in subcooled flow boiling based on the dominant physical mechanisms under given thermal and flow conditions. The regimes are conveniently represented in a boiling curve. A boiling curve was initially proposed by Nukiyama [16] for pool boiling of water at saturation temperature and at atmospheric pressure. Nukiyama's boiling curve depicts the relationship between the wall surface temperature in excess of the saturation temperature of the liquid, commonly known as wall superheat, and the wall heat flux. A representation of Nukiyama's boiling curve is shown in Fig. 1 a. First, the boiling regimes in pool boiling are discussed followed by flow boiling. When the wall surface temperature is below the saturation temperature of the liquid, heat is transferred purely by single phase convection. This is known as the free convection or natural convection regime in the case of pool boiling, A − B in Fig. 1 (a). The first bubble nucleates on the heated surface as the wall surface temperature rises above the liquid saturation temperature. This point is identified as onset of nucleate boiling, point B in Fig. 1 (a). The corresponding wall temperature is denoted by T ONB . In other words, the onset of boiling requires a few degrees of wall superheat. The fate of nucleating vapour bubbles depends on the bulk temperature of the coolant. If the bulk temperature is below the saturation temperature (subcooled boiling), the vapour bubbles condense into the bulk of the liquid coolant, thus promoting transport of latent heat from the heater surface to the liquid bulk. If the bulk temperature of the coolant is equal to its saturation temperature, i.e., saturated boiling, the vapour bubbles leaving the heated surface can no longer condense and thus contribute to a significant increase in vapour fraction of the system.
After the onset of boiling, the population of nucleating vapor bubbles increases steadily with rising wall superheat, in what is known as the nucleate boiling regime . Growth and transport of nucleating vapor bubbles enhance the heat transfer rate. With increasing wall superheat, the vapor bubbles begin to interact and coalesce into vapor patches and vapor columns. This eventually leads to formation of an unstable vapor blanket. The low thermal conductivity of vapor causes sharp decrease in heat transfer. The maximum heat flux, at point D in Fig. 1 (a), is known as critical heat flux ( q CHF ). Increase in wall superheat beyond q CHF results in the transition boiling regime. The path taken in the transition boiling regime (path D-E-D' or path D-D') depends on whether the wall surface temperature or the wall heat flux is controlled. When the wall surface temperature is controlled, increase in wall surface temperature beyond q CHF results in sharp decrease in wall heat flux along the path D − E in Fig. 1 (a), due to the low heat transferred by the unstable vapor blanket. This was initially guessed by Nukiyama and was later confirmed by other researchers [7] . Any futher increase in wall surface temperature leads to formation of a stable vapor blanket, denoted by regime E − F in Fig. 1 a, characterized with heat transfer by both natural convection and radiation through the vapour. This is known as the film boiling regime. When the wall heat flux is controlled, the increase in heat flux beyond q CHF results in sharp increase in wall temperature (along DF in Fig. 1 a) leading to film boiling. This behaviour was observed by Nukiyama [16] in his measurements where the heat flux was controlled.
In the case of flow boiling, the heat transfer in the single-phase regime before the onset of nucleate boiling is due to forced convection, which is higher than the heat transfer due to natural convection observed in pool boiling. The heat transfer by forced convection increases with increase in bulk velocity of the coolant, see Fig. 1 (b). Also, the onset of boiling, B 1 in Fig. 1 (b) is offset to a higher temperature for higher coolant bulk velocities. With increasing input heat flux or wall surface temperature, the different flow boiling curves eventually merge into the pool boiling curve at point C. The regime B − C is known as the partially developed boiling regime. In line with the numerical models dealt with in this work, a sub-regime within the partially developed boiling regime is highlighted. This sub-regime is known as the isolated bubbles regime, B − B in Fig. 1 (b). It is characterized by vapor bubbles that are isolated and not affected by neighbouring bubbles. Thereby, point B in the boiling curve denotes onset of bubble interactions. Increase in bubble interactions mitigate the effect of forced convection on heat transfer. Beyond point C, the effect of forced convection on heat transfer is marginal and the regime is known as the fully developed boiling , FDB, regime.
Understanding the physics governing each of the boiling regimes and the knowledge of an upper limit to avoid occurrence of film boiling is essential for extracting the potential in boiling and for development of relevant numerical models.

Heat flux partitioning and the BDL model
The presence of both forced convection and nucleate boiling in subcooled flow boiling enables the wall heat flux to be expressed using the heat flux partitioning approach. In this approach the total wall heat flux, q wall , is expressed as the sum of heat flux due to forced convection, q f c , and heat flux due to nucleate boiling, q nb . Chen's model [12] based on the heat flux partitioning approach was developed for analysing saturated flow boiling, where, the total wall heat flux is expressed as The forced convection heat flux, q f c , and nucleate boiling heat flux, q nb , are given by The forced convection heat transfer coefficient h f c in Eq. (2) is obtained from the Dittus-Boelter correlation [17] , Eq. (3) . The nucleate boiling heat transfer coefficient h nb in Eq. (4) is obtained from the Foster and Zuber correlation [18] , Eq. (5) . In Eqs. (2) -(7) , Re, P r l and Nu are the bulk Reynolds number, Prandtl number and Nusselt number, respectively, and L c is a characteristic length in the domain. The fluid properties k l , c p,l , ρ l , μ l denote the liquid thermal conductivity, specific heat, density and dynamic viscosity, respectively and P sat denotes saturation pressure of the liquid. The vapor density is denoted by ρ g , the surface tension by σ and latent heat of vaporization by h lg .
The terms F and S in Eq. (1) are two factors introduced by Chen, which denote the enhancement of forced convection due to presence of vapor bubbles and suppression of nucleate boiling due to the fluid motion, respectively. Chen obtained F and S based on bulk flow parameters. The BDL model by Steiner et al. [15] is an improvement of Chen's model that uses local parameters and a mechanistic approach to compute S.
While F = 1 can safely be considered for subcooled flow boiling, the BDL model estimates the flow induced suppression of nucleate boiling based on the following experimental observations: • In a horizontal channel flow with a heater at its bottom surface, isolated vapor bubbles initially attached to their respective nucleation sites are inclined along the direction of the flow; • They grow up to a certain size and then depart from their respective nucleation sites by sliding along the heater surface; and • The departed bubbles become upright, continue to grow until they depart from the heater surface and rise into the bulk liquid in the direction normal to the heater surface.
Klausner et al. [19] and Zeng et al. [20] observed this behaviour of vapor bubbles and numerically modelled this effect by performing a static force balance on an isolated bubble. Their models resulted in an estimation of the bubble radii at the instant when the bubble starts to slide -known as departure -and at the instant when it rises into the bulk liquid -known as lift-off. Similar behaviour of vapor bubbles was observed by Steiner et al. [15] in their own experiments. Based on these observations, they proposed the suppression of nucleate boiling due to the flow as the ratio of vapor bubble radius at the instance of departure ( r d ) to that at the instance of lift-off ( r l ).
In support of the above definition of the suppression factor, Steiner et al. [15] mentioned that, in case of pool boiling, the bubbles leave the nucleation site by lifting off the heater surface, i.e., without sliding, resulting in S f low = 1 . They used the static force balance equations proposed by Zeng et al. [20] . Fig. 2 shows the forces acting on an isolated bubble attached to its nucleation site. The static force balance equations read, at the instant of departure in the x and y directions, respectively. F d , F du , F sl and F bcy are the forces due to quasi-steady drag, unsteady drag due to asymmetric bubble growth, shear lift and buoyancy, respectively. The inclination of the vapor bubble with respect to the vertical axis is denoted by θ , see Fig. 2 . The force balance equation at the instant of lift-off is given by In line with the experimental observation discussed previously, the departed bubble becomes upright and hence, θ is not present in the force balance at lift-off, in Eq. (11) . Each of the aforementioned forces is a function of the vapor bubble radius ( r). The expressions for the forces involved in Eqs. (9) -(11) read, In these expressions, u denotes the velocity of the liquid at the centre of the bubble and du dy denotes the spacial derivative of u in the wall normal direction, y . These are obtained from the velocity profile of a single phase flow assuming there is no bubble. Properties of the liquid such as density, dynamic viscosity and thermal diffusivity are denoted as ρ l , μ l and α l , respectively. The density of the vapor is denoted by ρ g . The expressions include two model constants, b and C s . The non-dimensional numbers, the Jacob number, Ja, and the bubble Reynolds number, Re b , are written as where, c p,l and h lg are the specific heat of the liquid and latent heat of vaporization, respectively.
In addition to S f low , given in Eq. (8) , Steiner et al. [15] introduced another factor S subcool , which suppresses nucleate boiling due to subcooling of the liquid given by where, T wall , T sat and T bulk are the wall temperature, liquid saturation temperature and liquid bulk temperature, respectively. Finally, the suppression factor, S, used in Eq. (1) , is given by

Modified BDL model
In the original BDL model, the quasi-steady drag force ( Eq. (12) ) and shear lift force ( Eq. (13) ) are computed for a bubble in an unbounded flow field. The presence of the wall onto which the bubble is initially attached is not accounted for. Therefore, in the present work, the expressions proposed by Mazzocco et al. [21] , for drag and lift forces acting on a bubble in contact with the wall, are incorporated in the BDL model. The modified expressions for the drag and lift forces that replace Eqs. (12) and (13) , respectively, read,

Rohsenow's correlation for fully developed boiling
Fully developed boiling is encountered for higher wall superheats, see Fig. 1 (b). This boiling regime, as mentioned earlier, is characterized by that the heated surface is densely populated by vapor bubbles. As a result, the effect of forced convection on heat transfer becomes negligible. Therefore, a pool boiling correlation can be used to estimate the heat flux in this regime [11,22] . The pool boiling correlation by Rohsenow [11] is used in the present study to estimate fully developed boiling heat flux ( q F DB ). Rohsenow's correlation reads, where, T sat is the wall superheat, P r l and σ are the liquid Prandtl number and surface tension, respectively. The correlation also includes three model constants: c s f , n p and m . The constant C s f is

Blending of the models
Using different formulations for partially developed boiling and fully developed boiling regimes in modeling of subcooled flow boiling heat transfer is a common practice. The distinction between these two regimes is accounted for by identifying a transition point along the boiling curve. Kandlikars boiling model [14] , which uses the division description method, is a classic example for this approach. In Kandlikar's model, the intersection of the forced convection and fully developed boiling curves is identified and the heat flux at this point is multiplied by a factor of 1.4 and the resulting heat flux is defined as the transition point. Shah et al. [13] determined the transition point based on the boiling number. Prodanovic et al. [23] also used a boiling number based method to In the present study, a mathematical expression is sought for to blend the BDL model and Rohsenow's correlation. This expression handles the transition from the isolated bubbles regime to the fully developed boiling regime. The increase in population of vapor bub-bles on the heater surface and their resulting interaction is a dominant mechanism that distinguishes these two regimes. This should be accounted for in the mathematical expression that blends the two models.
Prior to arriving at this desired mathematical expression, the 'Dry spot model' developed by Ha and No [24] is discussed here in brief. Ha and No predict the occurrence of critical heat flux and transition boiling based on bubble interactions. They propose that, presence of a given number of vapor bubbles, or active nucleation sites, within a specified area on the heater can lead to formation of a dry spot. Subsequent merging of these dry spots eventually leads to critical heat flux and transition boiling. The bubble interactions are quantified based on probability density functions for spatial distribution of active nucleation sites. The nucleating vapor bubbles are assumed to follow a spatial Poisson distribution on the heated surface. Gaertner [25] initially represented the distribution of active vapor bubble nucleation sites on a heater surface using a spatial Poisson distribution. In the Dry spot model, the probability of presence of more than four active nucleation sites, within an area covered by two bubble diameters on the heated surface, is the criterion deduced for dry spot formation [24] . They support this criterion by stating that the surrounding bubbles cut-off the supply of liquid to the central vapor bubble through its microlayer and hence result in formation of a dry spot, as shown in Fig. 3 .
This idea of quantifying bubble interactions is used in the present study. It is proposed that the probability ( ) of occurrence of more than one active nucleation site within an area covered by two bubble diameters quantifies bubble interactions.
Given there exists an active nucleation site within the specified area A c on the heater surface, the expression for probability of having more than n c active sites [24] reads, where N denotes active bubble nucleation site density, the number of nucleation sites per unit area. Ha and No compute N using the Kocamustafaogullari and Ishii [26] correlation. More advanced models with improved prediction of active nucleation site density have been developed since then and are available in the literature [27][28][29][30][31] . In the present model, the correlation proposed by Li et al. [31] is used for computing N.
where, N 0 is an empirical constant, P is the pressure (in MPa ), φ is the contact angle of the vapor bubble, T c is the critical temperature at which the contact angle becomes zero, T 0 is the room tempera- In the present study, the probability of having more than one bubble within the area ( A c ) covered by two bubble diameters, indicating bubble interaction, reads, (31) where, N is obtained from Eq. (25) . d a v is a time averaged diameter. The factor S subcool accounts for the effect of liquid subcooling on growing vapour bubbles [15] , i.e., the growing vapor bubbles experience condensation at the bubble tip due to the subcooled bulk liquid. Hence, S subcool influences the bubble diameter. Note the expression for bubble growth rate, bubble radius as a function of time, is presented in Eq. (14) . Ha (32) where, S is provided in Eq. (19) . The probability of bubble interaction, , weighs the relative importance of the BDL model (in the isolated bubbles regime) and Rohsenow's pool boiling correlation (in the FDB regime), under given thermal and flow conditions. However, it is to be noted that the model does not result in any meaningful estimation of the wall heat flux once FDB is attained, i.e., = 1 . The model does not predict the CHF and boiling regimes beyond CHF. The procedure for calculating the wall heat flux using the new model is summarised in Fig. 4 .

Results and validation
The new blended model is validated with results from subcooled flow boiling experiments by Steiner et al. [15] , here after referred to as 'Experiment 1', and Lee and O'Niell [8] , here after referred to as 'Experiment 2'. Simplified academic geometries were used in these experimental investigations. The test sections are horizontal channels of rectangular cross section with a rectangular heater placed at the bottom face of the channel. While Steiner et al. [15] use an aluminium-alloy heater, Lee and O'Niell [8] use a copper heater. Both these experiments use water as coolant. The operating pressure, bulk coolant temperature and coolant bulk velocity are varied to study their effects on subcooled flow boiling heat flux. The wall heat flux versus the wall temperature data are available from the aforementioned experiments.
Empirical constants are present in the BDL model [15] , Rohsenow's correlation [11] and in Li's model [31] for active nucleation site density. Selected model constants are tuned to fit to available experimental data. In both Experiment 1 and Experiment 2, the constants m and C s f in the Rohsenow's model, and N 0 in In this study, an optimization framework, HAMON, based on evolutionary algorithms is used. HAMON is written in the programing language Python and can handle single-and multiobjective problems, as well as constrained and unconstrained ones [32,33] . It is freely available online via GitHub [34] . It uses evolutionary algorithms, either genetic algorithms or differential evolution, as the optimization method. These algorithms fall in the category of stochastic optimization methods, which in contrast to more conventional gradient descent or quasi-Newton methods, do not need the derivatives of the objective function(s) to be computed. They try to mimic the evolution process seen in nature where a population of individuals is advanced from generation to generation in order to try to improve the performance of the individuals.
For more details on the methods used, the interested reader is referred to Montero Villar et al. [32] .
In the present study, the objectives are to reduce the error and standard deviation with respect to the experimental data. For each of the experiments, the constants were optimized with a multiobjective differential evolution algorithm using 300 individuals and 500 generations. Even though this might be regarded as far beyond sufficient when it comes to the amount of individuals and generations needed, due to the low computational time required, this was not a problem.

Experiment 1
Different test conditions in experiments by Steiner et al. [15] , with an aluminium alloy heater and water as coolant, are summarized in Table 1 . In each of the plots in Figs. 5 and 6 , the experimental points bring out the initial linear forced convection regime followed by the non-linearity depicting the influence of nucleate boiling. This behaviour is well replicated by the numerical models too. In addition to this, the saturation temperature, T sat , and the temperature at onset of nucleate boiling, T ONB , are marked. The temperature at onset of nucleate boiling is computed using the criterion proposed by Hsu [35] . Delay in the onset of nucleate boiling with increase in bulk velocity is observed. The accuracy of the results has significantly improved after accounting for the presence of wall in the expression for the forces. This improvement is clearly visible in the regime after the onset of nucleate boiling for all cases except for the case of v bulk = 0 . 05 m/s in Fig. 5 . The departure radii predicted with the modified BDL model are within the limits of experimental measurements, except for the very low velocity case. This observation was in line with the agreement of departure radii predicted with the original BDL model by Steiner et al. [15] , with measurements.
In the case with v bulk = 0 . 05 m/s in Fig. 5 , the bulk velocity is very low and hence the boiling is more intense. This intense boiling characterized by multiple interacting vapor bubbles. Due to these bubble interaction, the suppression factor S f low in the BDL model, based on isolated bubble mechanics, has negligible influence on the wall heat flux. Thus, improvements made to the iso-   The improvement in estimation of wall heat flux for the v = 0.05m/s case is observed. The transition from the BDL model to the FDB model based on the probability of bubble interaction, also shown in plots, is evident. This probability of bubble interactions is high for the low bulk velocity cases and low for the high bulk velocity cases. The decrease in probability of bubble interactions in the high bulk velocity cases is attributed to enhanced suppression of nucleate boiling due to the flow. Moreover, owing to the higher saturation temperature at P operating = 2 . 0 bar, the onset of boiling is delayed compared to the low pressure case and the isolated bubbles regime is prevalent at higher wall temperatures. Therefore, for this operating pressure the results from the blended model are closer to the BDL model.

Experiment 2
Lee and O'Niel [8] studied experimentally the sensitivity of subcooled flow boiling heat flux to different parameters, such as liquid bulk velocity, system operating pressure and liquid subcooling. In this subsection the new boiling model is validated with the results from Lee and O'Niel's experiments, with a cooper heater and water as coolant, to showcase the sensitivity of the model to the aforementioned parameters. The different test conditions in their experimental study are summarized in Table 2 .
The Dittus-Boelter correlation with the default coefficients is shown in Eq. (3) . The agreement of experimental data points with the modified Dittus-Boelter model is shown in Fig. 10 b. In the BDL model, the default value of 20 / 3 is retained for the constant C s in Eq. (14) . There is no available data on bubble di-ameters for this experiment. Therefore, b = 1 . 0 , which is within the default range of values suggested by Zuber [37] in his diffusion controlled bubble growth model, is used in Eq. (14) . In the Rohsenow's correlation, n p = 1 is used since the coolant is water and the other two constants, C s f and m, are optimized using HA-MON to fit to the available experimental data. Furthermore, the constant N 0 in the model for active nucleation site density is also included in the optimization process. As a result of the optimization, C s f = 0 . 0145 , m = 2 . 9 and N 0 = 1120 are obtained. The constants remain the same for all the test cases in Experiment 2.
Good agreement of the results from the new blended model for one of the test conditions is shown in Fig. 11 . The new blended model transitions from the Dittus-Boelter model in the pure forced convection regime to the BDL model at wall temperatures just above the saturation temperature ( 100 • C at P = 1 . 0 bar) and fur-    Fig. 13 the liquid bulk temperature is kept constant at T bulk = 90 • C, in Fig. 14 the liquid subcooling is kept constant at T sub = T sat − T bulk = 10 • C. Lee and O'Niel visually observed decrease in bubble population with increasing operating pressure for constant liquid bulk temperature [8] . They established in their study that, this was due to increase in subcooling at the higher pressures, since the saturation temperature of the coolant increases with increasing operating pressure. The results from the new model are consistent with this observation; low values of probability of bubble interaction are estimated at higher pressures in Fig. 13 b. The effect of pressure alone is shown in Fig. 14 , where the subcooling is kept constant by varying the bulk liquid temperature. The model deviates significantly from a few data points. How-  Fig. 15 a. But, the heat flux in the FDB regime is over-predicted by the Rohsenow's correlation at P = 2 . 0 bar (see Fig. 16 a). Such an over-prediction of heat flux at higher pressure was not observed in Experiment 1. Furthermore, Lee and O'Niel visually observed fewer bubbles on the heated surface when the subcooling was higher. Figs. 15 b and 16 b confirm this observation, where the probability of bubble interaction is lower for higher subcooling.
Baring the over-prediction by Rohsenow's pool boiling correlation at higher values of operating pressure, the validation with results from Experiment 2 brings out the responsiveness of the model to the three different input parameters, i.e., velocity, subcooling and operating pressure.

Discussion
The new model estimates subcooled flow boiling heat flux across different boiling regimes and is responsive to different input parameters. The model includes a formulation that indicates vapor bubble interaction which forms the basis for the blending function. The blending function, in its current form, is a linear function that weighs the relative importance of the BDL model and Rohsenow's correlation, expressed by Eq. (32) . Improvement of the blending function with other mathematical formulations is for a future study.
An important concern in extracting the potential in boiling, in many practical applications, is the need for an upper limit in order not to encounter q CHF and film boiling. This concern is addressed in the present study, using a conservative approach, with the probability function. The value of the probability function represents proximity to different boiling regimes; a value close to zero represents proximity to the isolated bubbles regime and a value close to one represents proximity to FDB regime. We know from Section 2 that increasing the wall temperature beyond the FDB regime leads to encountering critical heat flux and eventually film boiling. Therefore, using the probability function to detect the intensity of boiling is proposed for practical applications. In this way, extracting the potential in boiling is limited to the partially developed boiling regime. It is important to note that the model does not provide any more information about the q CHF or proximity to q CHF once FDB ( = 1 ) is attained.
It is also observed that, the material of the heater surface impacts boiling heat flux and hence the correlations used to estimate it. For example, Hua et al. [9] and Abou-Ziyan [38] conducted subcooled flow boiling experiments with cast iron heating surfaces. The wall superheats involved in their experiments are as high as 100 • C, which is twice that encountered in the experiments with aluminium alloy by Steiner et al. [15] and copper by Lee and O'Niel [8] . In terms of modelling, the Foster and Zuber correlation [18] , with the default values of exponents, used in the BDL model, overpredicts the nucleate boiling heat transfer coefficient, h nb , for cast iron, shown in Fig. 17 . This over-prediction of h nb requires excessive flow induced suppression of boiling, i.e. a very small value for S f low , see Eqs. (8) and (19) . A small value of S f low results in un-physically small values of bubble departure diameter. Thus, the BDL model and the new blended model, in their current form, are not suitable to estimate the wall heat flux for cast iron, where high wall superheats are involved.
It is evident that the material of the heated surface influences boiling. Given below are the thermo-physical properties [36] (thermal conductivity k, density ρ, heat capacity c p and thermal diffusivity α) of the heated surface materials relevant to the current study, at 300 K. The BDL model is based on the force balance on a vapor bubble in horizontal channel flow with the heater at the bottom face of the channel. Therefore, we should have this in mind while using the model for other combinations of orientation of the heater surface with respect to the flow of coolant. This is critical, especially, when the heater surface is at the top face of the channel and buoyancy keeps the vapor bubble attached to the heater.
The new boiling model uses local input quantities except in the model for estimating forced convection heat flux, the Dittus-Boelter model. However, when implementing the boiling model within a CFD solver, the forced convection heat flux is obtained from the solver by solving the energy equation and the Dittus-Boelter model is not required. This also negates the need to use a different set of coefficients in the Dittus-Boelter model for different experiments (see Section 4 ). The requirement of local input quantities makes the integration of the new boiling model in a CFD solver straightforward.

Conclusions
A new numerical model to estimate subcooled flow boiling heat flux is introduced. The model is based on two well-established boiling models in the literature and additionally accounts for interaction of vapor bubbles. The concern for having an upper limit to avoid film boiling is addressed, in a conservative way, by using the probability of bubble interactions which indicates proximity to the FDB regime. Thereby, the boiling is limited to the partially developed boiling regime. Furthermore, there is scope for improving the mathematical formulation of the blending function in the model. The model employs local input quantities and therefore can be integrated with CFD simulations. The CFD simulations, as part of a complete thermal analysis framework, provide the platform to test and implement precision cooling strategies for high power density applications. The limitations of the model in estimating wall heat flux for different orientations of the heated surface and coolant flow is realized. Moreover, the Foster and Zuber correlation used in its current form over-estimates the nucleate boiling heat transfer coefficient when high wall super heats are encountered. Although the model can be implemented in CFD simulations, we should keep the limitations in mind while analysing the results.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.