Evaluation of Interfacial Heat Transfer Models for Flashing Flow with Two-Fluid CFD

Abstract

The complexity of flashing flows is increased vastly by the interphase heat transfer as well as its coupling with mass and momentum transfers. A reliable heat transfer coefficient is the key in the modelling of such kinds of flows with the two-fluid model. An extensive literature survey on computational modelling of flashing flows has been given in previous work. The present work is aimed at giving a brief review on available theories and correlations for the estimation of interphase heat transfer coefficient, and evaluating them quantitatively based on computational fluid dynamics simulations of bubble growth in superheated liquid. The comparison of predictions for bubble growth rate obtained by using different correlations with the experimental as well as direct numerical simulation data reveals that the performance of the correlations is dependent on the Jakob number and Reynolds number. No generally applicable correlations are available. Both conduction and convection are important in cases of bubble rising and translating in stagnant liquid at high Jakob numbers. The correlations combining the analytical solution for heat diffusion and the theoretical relation for potential flow give the best agreement.

1. Introduction

Flash boiling is a vaporization process triggered by depressurization instead of heating, which is relevant to a number of industrial economic and safety concerns. For example, in the automobile industry, the atomization of fuel spray in a combustion chamber is affected significantly by its flashing characteristics inside the injector nozzle [1]. In the nuclear industry, during the hypothetical loss of coolant accident of pressurized water nuclear reactors, the rate of coolant loss is determined by the critical flashing flow through the crack [2]. In the chemical industry, the severity of failure of pressurized vessels or pipes containing liquefied chemical hazardous gases is characterized by the external flashing flow [3]. An additional flashing evaporation phenomenon was reported in [4], which refers to an aerospace application and concerns the leading edge cooling of a space vehicle. Another similar phenomenon often encountered in case of pressure variation is cavitation. In general, cavitation occurs at relatively low temperature levels, where bubble growth is controlled mainly by the pressure difference across the interface. In contrast, flashing of hot fluids is more like a boiling process, which is driven principally by the thermal non-equilibrium. The complexity of flashing flows is represented by gas-liquid mixture with rapid phase change and bubble dynamics [5], and numerical studies are directed towards the determination of vapour generation rate. Good reviews have been given by Pinhasi [3] and Liao & Lucas [6]. In general, two methods have been used for the evaluation of the interfacial mass transfer rate in flashing flows. One is based on the observation of non-equilibrium mechanical and thermal processes. The other treats the transition of the thermodynamic system from non-equilibrium to equilibrium as a relaxation process. The two states are bridged by means of an empirical coefficient, i.e., the relaxation time [7,8,9,10]. The present paper will focus on the former one, which is consistent with the two-fluid framework. Under this category there are again two alternatives having been adopted for the estimation of vapour generation rate. One is based directly on the interfacial heat transfer process

$$\dot{m}=\frac{\dot{q}}{L},$$

(1)

where $\dot{m}$ is the mass flux, L the latent heat of vaporization, and $\dot{q}$ the heat flux transferring from the vapour and liquid to the phase interface. For vapour-liquid such as steam-water flows under most practical conditions, the interfacial heat transfer on the vapour phase is usually much smaller (less than $5\%$) than that on the liquid phase [11]. Therefore, it is usually neglected by assuming that the temperature is uniform inside the bubble and equal to that at the interface. This assumption is also made in the current work

$$\dot{q}=h^l{T}_{sup},$$

(2)

where $h^l$ is the overall heat transfer coefficient between the superheated liquid and the liquid-vapour interface, and $T_{sup}$ is the superheat degree of the liquid.

An alternative approach is formulated in terms of the resultant bubble growth rate

$$\dot{m}={\rho }_v\dot{R},$$

(3)

where $\dot{R}$ is the growth rate of bubble radius given by an analytical solution, and ${\rho }_v$ vapor density. Bubble growth in superheated liquid is known to be controlled successively by surface tension, liquid inertia and heat transfer [12]. The first stage is usually neglected in numerical analysis with the consideration of sufficient bubble size. The effect of liquid inertia is only important at the very early stage of depressurization [13] and for sufficiently small bubbles [14]. Therefore, the thermally controlled growth is of the greatest interest. In this domain the above two models are related to each other, and one gets

$$\dot{q}=h^l{T}_{sup}=L{\rho }_v\dot{R},$$

(4)

$$h^l=\frac{L{\rho }_v\dot{R}}{T_{sup}}.$$

(5)

As a result, a primary concern of the numerical study on flashing flows turns out to be a reliable prediction of the interphase heat transfer coefficient or bubble growth rate. It is generally recognized that high uncertainty is present in choosing an appropriate heat transfer coefficient correlation for two-fluid computational fluid dynamics (CFD) simulations [6]. One major reason is that the insufficiency and limitation of the correlations is not completely identified, and a quantitative evaluation is missing. This work aims to present a thorough literature survey on existing theories and correlations, and evaluate their performance by carrying out CFD simulations and making comparisons with available experimental and Direct Numerical Simulation (DNS) data. Note that although the background of the present work is flashing flow, the results and discussions are not restricted to it. Certain similarities have been identified in the transfer scenarios of a liquid sphere exposure to blowing hot gas and a bubble rising in superheated liquid or dissolving in liquid. The correlations for heat (or mass) transfer in condensation, evaporation and dissolution are often exchangeable. A variety of correlations are available in the literature. They take into account the effect of conduction, convection and turbulence partially or totally, but mostly in a pure empirical or semiempirical way. A good review was given by Mathpati & Joshi [15]. An overview of the classical theories, analytical solutions and empirical correlations is given below.

2. Theories and Analytical Solutions

2.1. Conduction

Numerous analytical solutions are available for the heat transfer between spherical vapour bubbles and the surrounding liquid. Most of them account for the simplified heat conduction problem and neglected the momentum effect in the liquid and gas phases on the bubble growth and its shape. It states that the bubble growth problem is analogous to a one-dimensional, unsteady state heat diffusion process with moving boundary, which is described by

$$\frac{\partial T_l}{\partial t}=a_l\frac{{\partial }^2{T}_l}{\partial x^2},$$

(6)

where $a_l$ is liquid thermal diffusivity, t the time coordinate, the direction x normal to the boundary surface, and $T_l$ is the temperature of the surrounding liquid. Under certain initial and boundary conditions, the temperature field of the liquid around the bubble surface may be solved from Equation (6) analytically. The assumption of a thin “thermal boundary layer”, i.e., the change of liquid temperature taking place only in a thin film adjacent to the interface, is often made in these solutions. The liquid temperature in the bulk, $T_0$, is uniform and constant. In addition, for constant pressure fields, the saturation temperature $T_{sat}$ at the interface remains unchanged, and the vapor inside the bubble is often assumed to have the saturation temperature. The theory of thermal boundary layer is shown schematically in Figure 1.

Knowing the temperature distribution, the heat flux transferring from the bulk to the bubble surface can be determined from the Fourier’s Law

$$\dot{q}={\lambda }_l{\left.\frac{\partial T_l}{\partial x}\right|}_{x=R},$$

(7)

where ${\lambda }_l$ is liquid thermal conductivity and R the bubble radius. Finally, the bubble growth rate $\dot{R}$ is obtained by substituting $\dot{q}$ into Equation (4). Fritz and Ende [16] solved the heat conduction across a semi-infinite plane slab under constant temperature boundary conditions. They derived the asymptotic bubble radius

$$R\left(t\right)=\frac{2}{\sqrt{\pi }}·{Ja}_T·{\left(a_{l}t\right)}^{1/2},$$

(8)

where the Jakob number, ${Ja}_T$, is defined as

$${Ja}_T=\frac{{\rho }_l{c}_{p,l}{T}_{sup}}{{\rho }_{v}L}.$$

(9)

Combining with Equation (5) the heat transfer coefficient for conduction is obtained

$$h_{cond}^l=\frac{{\lambda }_l}{\sqrt{\pi a_{l}t}},$$

(10)

wherein the time t is related to ${Ja}_T$ number and bubble radius via Equation (8). In terms of the dimensionless number, the Nusselt number Nu, above equation is expressed as

$${Nu}_{cond}=\frac{4}{\pi }{Ja}_T.$$

(11)

Unsteady heat conduction across spherical bubble surfaces were studied by Plesset & Zwick [17] and Forster & Zuber [13] independently. Both solutions are in the same form as Equation (8) with the exception of a so-called “spherical factor”, $K_s$, i.e.,

$$R\left(t\right)=K_s\frac{2}{\sqrt{\pi }}·{Ja}_T·{\left(a_{l}t\right)}^{1/2}.$$

(12)

The numerical constant $K_s$ is greater than 1, which means that under the same temperature difference, heat flux across a spherical bubble is larger than a planar surface because the temperature gradient in the “thermal boundary layer” is increased by the curvature. In [17]$K_s=\sqrt{3}$ while in [13] $K_s=\pi /2$. Olek et al. [18] derived an alternative expression for the heat flux at the boundary of a sphere by using the hyperbolic heat conduction equation. For long times, the asymptotic solution approaches those obtained by using the Fourier heat conduction, but with a correction factor like

$$K_s=\frac{1}{2}\left[1+{\left(1+\frac{2\pi }{{Ja}_T}\right)}^{1/2}\right].$$

(13)

With consideration of the correction factor $K_s$, the Nussellt number in Equation (11) turns into

$${Nu}_{cond}=K_s^2\frac{4}{\pi }{Ja}_T.$$

(14)

2.2. Convection

In principle, the above solutions without considering the effect of slip velocity are applicable for low void fraction and high superheat degrees. These conditions are expected to be satisfied only in a short time interval during a rapid depressurization, and at the transition of the flow changing from one-phase to two-phase [19]. For large bubbles one may expect a significant under-prediction by using these models. As observed and discussed in [20,21], the influence of slip velocity on the transfer rate is noticeable even in the case of bubbles rising in stagnant superheated liquid under normal gravity. The effect of slip velocity on interfacial heat transfer was firstly studied by Ruckenstein [22] and Sideman [23]. For the transfer between spherical independent vapour bubbles (influence of other bubbles and turbulence negligible) and the boiling liquid in motion, they suggested in a potential flow

$$h_{conv}^l=\frac{{\lambda }_l}{R}·\frac{1}{\sqrt{\pi }}·{Pe}^{1/2},$$

(15)

or in terms of Nusselt number

$${Nu}_{conv}=\frac{2}{\sqrt{\pi }}·{Pe}^{1/2},$$

(16)

where the Péclet number is defined by

$$Pe=\frac{d\left|{\overrightarrow{U}}_{rel}\right|}{a_l}={Re}_p{Pr}_l,$$

(17)

with the particle Reynolds number ${Re}_p=d\left|{\overrightarrow{U}}_{rel}\right|/{\nu }_l$ and the liquid Prandtl number ${Pr}_l={\nu }_l/a_l$. Equation (15) or (16) is often interpreted by the so-called “penetration theory” [24]. It states that the liquid molecules in contact with the bubble surface are replaced at a constant time interval, which can be expressed as a ratio of bubble diameter to relative velocity

$${\tau }_{conv}\approx \frac{d}{\left|{\overrightarrow{U}}_{rel}\right|}.$$

(18)

The schematic representation of the penetration theory is shown in Figure 2.

Similar expressions can be obtained from the “penetration” theory and the “thermal boundary layer” theory [23] discussed above. The heat transfer takes place in a thin laminar sublayer, where a constant velocity can be assumed regardless of the hydrodynamics in the bulk of liquid. By substituting Equations (17) and (18) to Equation (15), one gets

$$h_{conv}^l=2·\frac{{\lambda }_l}{\sqrt{\pi a_l{\tau }_{conv}}},$$

(19)

which has the same form as the conduction transfer coefficient given in Equation (10) except the factor 2 and the characteristic time scale.

2.3. Effect of Turbulence

The effect of wake and freestream turbulence on the transfer from spheres has received relatively less attention, and is still not well understood. Conflicting arguments and observations exist. It is commonly believed that the transfer is enhanced in the presence of wake interactions and freestream turbulence. The experimental investigation on heat transfer from solid spheres reported by Lavender & Pei [25] and Raithby & Eckert [26] showed that the Nusselt number increased with increasing turbulence intensity in the ambient flow. Yearling & Gould [27] measured the convective heat and mass transfer rates from liquid droplets in turbulent air flow and also observed that the Nusselt number increased with increasing turbulence intensity. However, the augmentation was not duplicated by the experiment on the evaporative heat and mass transfer of suspended heptane droplets performed by Buchanan [28]. Theoretical interpretation of the turbulence effect is mainly based on the so-called “surface renewal theory”, which is a modification of the “penetration theory” discussed above. As described in [22], if intense turbulent motions appear in the liquid, the contact surface between the liquid and the bubbles is continuously renewed by turbulence eddies, which brings the liquid from the bulk to the interface at average intervals of ${\tau }_{turb}$, see Figure 3.

The characteristic time is determined by the turbulent fluctuation velocity instead of the slip velocity used by the “penetration theory”. Dackwerts [29] suggested that turbulence renews the volume elements at the interface, but the turbulence dies out as it approaches the interface. As a result, the liquid element itself has a non-turbulent structure, and the heat or mass transfer at the interface still has a molecular character. The heat transfer coefficient may be expressed as

$$h_{turb}^l=2{\lambda }_l·\frac{1}{\sqrt{\pi a_l{\tau }_{turb}}}.$$

(20)

It has the same form as Equation (19) except the time scale ${\tau }_{turb}$, which is often estimated using the ratio between Kolmogorov length and velocity scales [30,31,32], i.e., ${\tau }_{turb}=\sqrt{\frac{\nu }{ϵ}}$. That means that the boundary layer adjacent to the interface is renewed by near-surface small scale eddies. In contrast, a large eddy model uses the scales of energy containing eddies [33]. Some researchers suggested that the small eddy model is more valid at higher Reynolds numbers while the large eddies dominate the surface renewal at lower Reynolds numbers, and thus proposed a two-regime model. However, there is no consistent definition of the Reynolds number and the criterion for transition. On the other hand, Sideman [34] approximated it as a ratio of the bubble diameter d to the fluctuation velocity at a distance of d

$$\overline{U^{\prime }}\sim \left\{\begin{array}{c}{\left(\epsilon d\right)}^{1/3}\mathrm{for}\mathrm{inertial}\mathrm{regime}(\eta <d\le l)\hfill \\ d{\left(\frac{\epsilon }{\nu }\right)}^{1/2}\mathrm{for}\mathrm{viscous}\mathrm{regime}(d\le \eta )\hfill \end{array}\right.,$$

(21)

where $\eta ,l$ is the Kolmogorov and integral length scale, respectively. Owing to advances in numerical algorithms and high performance computing, the transfer process occurring at the interface becomes amenable with the aid of DNS. Figueroa-Espinoza and Legendre [35] investigated the effect of bubble aspect ratio and bubble wake on the mass transfer from oblate spheroids by DNS solving the Navier-Stokes equations. They found that most of the transfer occurs on the front part of the bubble. The contribution in the wake region increases as the aspect ratio increases. The local transfer rate at the bubble surface as a function of the azimuthal angle deviates significantly from Equation (15) due to unsteady effects from vorticity production and wake destabilization. However, the total transfer rate expressed as the Nusselt number is shown to satisfy well the potential flow theory if the equivalent diameter is used as the characteristic length scale, see Figure 4.

Similar results about the wake effect on the transfer rate were obtained by Bagchi & Kottam [36] in their DNS simulation of heat transfer from a sphere in a turbulent flow. In addition, the freestream turbulence was shown to have a clear influence on the instantaneous and local Nusselt number. However, the time and surface averaged Nusselt number was found insensitive to the ambient turbulence and it can be predicted by correlations for steady and uniform flow.

3. Empirical Correlations

Besides analytical solutions, some researchers have proposed useful empirical correlations for calculating the heat transfer coefficient in flashing conditions. For conduction, a widely used expression in terms of the Jakob number was presented in [37]

$$Nu=\left[2+{\left(\frac{6{Ja}_T}{\pi }\right)}^{1/3}+\frac{12}{\pi }{Ja}_T\right]$$

(22)

The correlation was validated for bubble growth in uniformly heated liquid, and adopted in [38,39,40] for the modelling of various flashing flows. A slightly modified expression was presented later on in [41].

$$Nu=\left[2+{\left(2{Ja}_T\right)}^{1/3}+\frac{12}{\pi }{Ja}_T\right].$$

(23)

The Ranz-Marshall correlation [42], which was proposed based on experimental data for spherical water drops evaporating in blowing hot dry air

$$Nu=\left(2+C·{Re}_p^A{Pr}_l^B\right),$$

(24)

has been often used for estimating interphase heat transfer rates also in the case of flashing flows. Among many others examples are one-dimensional simulations presented by Richter [43], Bird et al. [44] and Dobran [45] and three-dimensional CFD simulations by Giese [46], Laurien [47] and Frank [48]. The empirical constants in Equation (24) are $A=1/2,B=1/3$ and $C=0.6$. Hughmark [49] suggested that these constants are valid for the range ${Re}_p<450$ and ${Pr}_l<250$. The exponents A and B may increase with the Reynolds number ${Re}_p$ and the Prandtl number ${Pr}_l$, respectively. Actually, correlations with slightly different constants have been widely used for convective transfer, e.g., in [50] $C=0.15$, while in [51,52] C was replaced by $0.46$ and $0.55$, respectively. At the same time, Lee and Ryley [53] found that the observations of water drops evaporating in superheated steam instead of air conform closely to Equation (24) but with $C=0.738$. For a vapour bubble freely oscillating in liquid, $A=B=0.5$, and $C=1.0$ according to Nigmatulin et al. [54] and Mahulkar et al. [55]. The expression of Aleksandrov et al. [56] was modified slightly by Saha et al. [57] for the calculation of heat transfer rate in flashing nozzle flows. They combine the conduction and convection transfer in the way

$$Nu={\left(\frac{{12}^2}{{\pi }^2}{Ja}_T^2+\frac{4}{\pi }Pe\right)}^{1/2}.$$

(25)

A similar expression was used by Wolfert [58] to simulate the rapid depressurization processes of high pressure pipes

$$Nu=\left(\frac{12}{\pi }{Ja}_T+\frac{2}{\sqrt{\pi }}{Pe}^{1/2}\right).$$

(26)

The above two correlations can be reformulated as

$$Nu=\sqrt{{\left({Nu}_{cond}\right)}^2+{\left({Nu}_{conv}\right)}^2},$$

(27)

and

$$Nu={Nu}_{cond}+{Nu}_{conv}.$$

(28)

The conduction and convection Nusselt number is evaluated by Equations (14) and (16), respectively. To account for the effect of turbulence, Wolfert et al. [59] introduced a so-called eddy conductivity, ${\lambda }_t$. The apparent thermal conductivity of liquid is given by

$${\lambda }_l^{\prime }={\lambda }_l+{\lambda }_t.$$

(29)

And the overall heat transfer coefficient is computed through

$$Nu=\left(\frac{12}{\pi }·{Ja}_T+\frac{2}{\sqrt{\pi }}\left(1+\frac{{\lambda }_t}{{\lambda }_l}\right)·{Pe}^{1/2}\right).$$

(30)

The eddy conductivity ${\lambda }_t$ was assumed to be dependent on the liquid velocity. In case of one velocity component $W_l$, ${\lambda }_t$ is expressed as

$${\lambda }_t={\lambda }_l·{\chi }_t·W_l,$$

(31)

where ${\chi }_t$ is an empirical constant. Based on a pressure release experiment on a test vessel and two blowdown experiments on the vessel and straight pipe, Wolfert et al. [59] found that ${\chi }_t=0.8$sm${}^{-1}$ gives the best agreement between calculated and experimental results. Whitaker [60] interpreted that the enhancement in transfer rates from a sphere due to the presence of turbulence comes purely from the wake contribution, while the transfer process at the front surface can be described by the law for potential flow, i.e., the Nusselt number $Nu\propto {Re}_p^{1/2}$, and the freestream turbulence has no effect. In the wake region, the functional dependence for the Reynolds number is ${Re}_p^{2/3}$, and it is cumulative to the laminar part

$$Nu=\left(2+\left(0.4·{Re}_p^{1/2}+0.06·{Re}_p^{2/3}\right){Pr}_l^{0.4}\right).$$

(32)

Another empirical method often used to account for the turbulence enhancement is to increase the constant A in Equation (24). Issa et al. [61] compared the Nusselt number correlations for the case of saturated steam bubbles condensing in subcooled water, and found that the dependency upon Reynolds number grows as it increases, e.g., for ${Re}_p$ up to 800, $Nu\propto {Re}_p^{0.5}$, and for ${Re}_p$ up to ${10}^4$, $Nu\propto {Re}_p^{0.7}$. They derived the following correlation for highly deformed large bubbles condensing in turbulent pipe flow,

$$Nu=0.0609{Re}_p^{0.89}{Pr}_l^{0.33}.$$

(33)

Based on a numerical study on the transient heat transfer from a sphere, Feng & Michaelides [62] obtained a simple correlation for the Nusselt number, which describes the dependence on Reynolds and Peclet numbers as follows

$$Nu=\left(0.922+{Pe}^{1/3}+0.1·{Re}_p^{1/3}{Pe}^{1/3}\right).$$

(34)

The DNS performed by Dani et al. [63] showed that the average mass transfer (heat transfer) from spherical bubble is affected not only by the Reynolds and Schmidt (Prandtl) numbers but also by the surface mobility and contamination. As shown in Figure 5a, for solid spheres or fully contaminated bubbles in a creeping flow (low ${Re}_p$ numbers) following expression proposed by Clift [64] reproduces the exact numerical solution with very high accuracy

$$Nu=1+{(1+Pe)}^{1/3}.$$

(35)

However, as the ${Re}_p$ number increases, the influence of ${Re}_p$ and ${Pr}_l$ are no more similar and have to be considered separately. In these cases, the Ranz & Marshall [42] correlation was shown to be able to give perfect agreements, see Figure 5b.

For clean spherical bubbles, Clift [64] suggested

$$Nu=1+{(1+0.564{Pe}^{2/3})}^{1/3},$$

(36)

which achieves excellent agreement with the DNS data at low Reynolds numbers, e.g., ${Re}_p\le 0.1$ (see Figure 6a). As shown in Figure 6b the aympotic solution for both ${Pr}_l\to \infty$ and ${Re}_p\to \infty$ agrees well with the potential theory, i.e., $Nu\sim {Pe}^{1/2}$ given by Equation (15). In this case, the Ranz & Marshall [42] correlation under-predicts the growth rate significantly.

4. CFD Simulation of Bubble Growth in Superheated Liquid

Bubble growing in stagnant liquid of uniform superheats has been adequately studied with high speed photography around the middle of the last century [20,65,66,67,68]. Dergarabedian [65] obtained bubble formation within the body of the liquid by heating a beaker of water slowly. In order to avoid wall nucleation, the Pyrex beaker surface was annealed carefully to be very smooth and free of pits. Nevertheless, thermal gradients were found to exist in a narrow boundary adjacent to the bottom of the beaker, and most of the bubbles were formed in this thermal layer. Hooper [66] and Florschuetz et al. [20] achieved uniform superheats by suddenly depressurizing heated pressurized water down to the atmospheric pressure within 5 ms. The bubble growth rates observed by Hooper [66] were found to be lower than the theoretical ones, and moreover, the deviation increased as the superheat was increased. Florschuetz et al. [20] investigated the effect of translational motion on the bubble growth rate at small superheat degrees under near zero- and normal-gravity conditions. The zero-gravity data was found to be suitable for the validation of theoretical solutions for heat-conduction controlled growth. In other words, the convective effects are negligible in these cases. On the other hand, the data taken at normal gravity indicate clearly that the enhancement of growth rates due to bubble translation becomes significant at later stages. Kosky [67] heated a tube of water uniformly in a silicone oil bath. The pressure in the test chamber was regulated with a vacuum pump, and recorded with a pressure transducer simultaneously.

In this section, the analytical and empirical correlations discussed above (see

Table 1. Correlations for estimating interphase heat transfer coefficient.

	Conduction
Reference	Correlation	Note
[16]	${\displaystyle Nu=\frac{4}{\pi }{Ja}_T}$	analytical
[17]	${\displaystyle Nu=\frac{12}{\pi }{Ja}_T}$	analytical
[13]	${\displaystyle Nu=\pi {Ja}_T}$	analytical
[18]	${\displaystyle Nu=\frac{{Ja}_T}{\pi }{[1+{(1+\frac{2\pi }{{Ja}_T})}^{1/2}]}^2}$	analytical
[37]	${\displaystyle Nu=2+{(\frac{6{Ja}_T}{\pi })}^{1/3}+\frac{12}{\pi }{Ja}_T}$	empirical
	Convection
[22]	${\displaystyle Nu=(\frac{2}{\sqrt{\pi }}){Pe}^{1/2}}$	potential theory
[56]	${\displaystyle Nu={\left(\frac{{12}^2}{{\pi }^2}{Ja}_T^2+\frac{4}{\pi }Pe\right)}^{1/2}}$	heuristic
[58]	${\displaystyle Nu=\left(\frac{12}{\pi }{Ja}_T+\frac{2}{\sqrt{\pi }}{Pe}^{1/2}\right)}$	heuristic
[42]	${\displaystyle Nu=2+0.6{Re}_p^{1/2}{Pr}_l^{1/3}}$	empirical
	Turbulence
[29]	${\displaystyle Nu=\frac{2}{\sqrt{\pi a_l{\tau }_{turb}}}}$	surface renew theory
[60]	${\displaystyle Nu=2+(0.4{Re}_p^{1/2}+0.06{Re}_p^{2/3}){Pr}_l^{0.4}}$	empirical
[59]	${\displaystyle Nu=\frac{12}{\pi }{Ja}_T+\frac{2}{\sqrt{\pi }}\left(1+\frac{{\lambda }_t}{{\lambda }_l}\right){Pe}^{1/2}}$	empirical
[62]	${\displaystyle Nu=0.922+{Pe}^{1/3}+0.1·{Re}_p^{1/3}{Pe}^{1/3}}$	empirical
[61]	${\displaystyle Nu=0.0609{Re}_p^{0.89}{Pr}_l^{0.33}}$	empirical

) are tested for the early stage of bubble growth in superheated liquid with the two-fluid CFD. Three configurations are considered, namely stationary bubble growth, translating bubble growth in stagnant and flowing liquid. These cases are chosen for validation considering the fact that the effect of heat conduction, convection and turbulence can be investigated separately to a certain extent, and many other uncertainties such as swarm effect and interphase momentum transfer can be excluded. Separate conservation equations are solved for the vapour and liquid phases. Further, the vapour is assumed to be saturated as done by many other researchers like Maksic [39] and Giese [46]. The particle model is applied for the mimic of interfacial morphology and computation of interfacial area density, where the vapour is modelled as spherical bubbles. A constant number concentration is assumed, which is consistent with the experimental observation by Florschuetz et al. [20] and valid for the early-stage of flashing. Microscopic or nearly microscopic bubbles are present naturally in the domain at the instant of pressure release, and provide the nuclei for subsequent phase growth. For further details about the numerical setup, the reader is referred to [69,70]. The presumed bubble concentration is found to have a negligible effect on the results, and a value of ${10}^4$m${}^{-3}$ is adopted in all the cases.

4.1. Stationary Bubble Growth

Two experimental test cases under atmospheric conditions at zero gravity are simulated with the CFD software ANSYS CFX (Version 18.0, ANSYS Inc., Canonsburg, PA, USA), where the buoyancy is deactivated and therefore there are no relative motion and convection effects. The experimental data are taken from the work of Florschuetz et al. [20] for steam-water systems. The water has an initial superheat degree of $2.9$K and $3.2$K, respectively. The simulation domain is a cube as shown in Figure 7. It is worth noting that the domain for the simulation of translating bubble growth in stagnant and flowing liquid is elongated in the stream direction.

No slip wall boundary conditions are applied to the bottom of the box. The top is set as opening, which allows both inflow and outflow of the gas phase. The other four sides (left, right, front, back) are treated as symmetrical planes. The results show that in the first ∼$500$ms the liquid temperature remains nearly constant. Meanwhile, the phase distribution and flow parameters are uniform in the domain with the exception of the region adjacent to the bottom wall. Therefore, the simulation condition conforms with the experiment that a bubble grows in stagnant uniform superheated liquid, for which an analytical solution of the growth rate is possible. As shown in Figure 8 the numerical results for the case $T_{sup}=2.9$K, which is averaged over the midplane of the domain(plane 1 in Figure 7), coincide with the analytical ones. This proves that the applied model is capable of simulating bubble growth in superheated liquid.

The comparison of the simulated and measured transient bubble size for the two cases is depicted in Figure 9a,b, respectively.

Note that in a log-log plot, the bubble growth exhibits a slope of 1/2 according to both theory and experiment. Nevertheless, linear-linear plots are used here in order to show the difference more clearly. For the case of $T_{sup}=2.9$ K the results obtained by using the correlations of Labuntzov et al. [37], Plesset & Zwick [17] and Forster & Zuber [13] agree well with the measurement, while the Olek et al. [18] and Fritz & Ende [16] correlations under-predict the bubble growth rate substantially. As the liquid superheat increases from $2.9$ K to $3.2$ K, at the early stage the experimental data are more consistent with the predictions of the former three correlations. However, asymptotically they approach to the latter two, see Figure 5b. It is interesting to note that the overall uncertainty of bulk liquid temperature values and equivalent bubble radii was estimated to be about $\pm 0.2$ K and $\pm 0.05$ mm [20]. Another DNS case of bubble growth in superheated under the zero-gravity condition presented in Ye [71] is simulated. It is based on the thermal properties of water under the atmospheric conditions. A liquid superheat of 1 K is considered, and the Jakob number is estimated as $3.0$. The predicted bubble growth rate is shown in Figure 10. The DNS results evidence that the asymptotic bubble growth follows the theoretical relation, $R\left(t\right)\propto t^{1/2}$. However, in the initial stage when the thermal boundary layer around the bubble is developing, the growth rate is evidently larger than the $1/2$ law (see Figure 6a). In contrast, the two-fluid simulation results obey the theoretical law well. The correlation of Fritz & Ende [16] is shown to be able to reproduce the DNS asymptotic results satisfactorily, while all the others are prone to over-predict the bubble size.

4.2. Translating Bubble Growth

Under the normal gravity condition, bubbles will be rising and translating simultaneously through the initially stagnant liquid, and the bubble growth rate is obviously larger than the diffusion limit $1/2$. In the simulation, the buoyancy model is activated, which results in a relative motion between the bubbles and the liquid. Both the conduction and convection play a role in the interphase heat transfer, but the turbulence effect is negligible. The momentum interaction is modelled by interphase drag force, and the drag coefficient is calculated according to the Ishii & Zuber correlation [72].

4.2.1. Florschuetz et al. Cases

Two cases from Florschuetz et al. [20] with $T=3.0K$ and $T=3.9K$ respectively are simulated by the two-fluid CFD with different correlations for the heat transfer coefficient. The results are presented in Figure 11. As expected, the relative motion accelerates bubble growing and results in a steeper slope of growth line than that under zero-gravity conditions. At the early stage ($t<5$ ms) the predictions given by Wolfert [58] and Aleksandrov et al. [56] agree well with the experimental data. Later on, the potential theory without consideration of heat conduction deliver the best agreement with the experimental data. It implies that heat conduction plays an important role at the initial stage, while convection becomes the dominant mechanism later. The results obtained by the correlation of Ranz & Marshall [42] under-predict the growth rate significantly in both cases. In addition, the asymptotic value of 2 is found to be insufficient to describe the heat transfer rate under the asymptotic condition of zero slip velocity, which is consistent with the findings of Walton [73]. In his experimental study on the evaporation of water droplets in hot air, Walton [73] measured the average Nusselt number of $3.8$ for natural convection and an asymptotic value of $6$ for forced convection.

4.2.2. Kosky Cases

Four cases from Kosky [67] with higher superheat degrees $(T=10.5\sim 23.2K)$ are also investigated. The results are shown in Figure 12. They have pressures other than $1atm$ and relatively high superheats in comparison with the last two cases. Here, a generally good agreement with the experimental data is demonstrated by the heuristic correlations of Aleksandrov et al. [56] and Wolfert [58]. On the other hand, the potential theory and the empirical correlation of Ranz & Marshall [42] under-predict obviously the transient bubble size.

The thermodynamic conditions of all above cases are summarized in

Table 2. Summary of test cases used for validation.

Case No.	$\mathit{p}\left[\mathit{atm}\right]$	${\mathit{T}}_{sup}\left[\mathit{K}\right]$	${Ja}_{\mathit{T}}$	error
1	$1.0$	$2.9$	$8.69$	$1.70$
2	$1.0$	$3.2$	$9.58$	$1.69$
3	$1.0$	$1.0$	$3.0$	$1.72$
4	$1.0$	$3.0$	$8.9$	$1.69$
5	$1.0$	$3.9$	$10.98$	$1.68$
6	$1.19$	$10.5$	$26.5$	$1.50$
7	$0.642$	$16.0$	$71.0$	$1.68$
8	$0.477$	$19.5$	$113.25$	$1.75$
9	$0.613$	$23.2$	$107.0$	$1.58$

.

4.3. Bubble Growth in a Flowing Liquid

Up to now the formulas obtained for stagnant liquid are commonly used for estimating the growth rate of vapor bubbles in high-velocity flowing liquid without checking the applicability. In his book, Avdeev [74] showed that these extrapolations are not justified and that they conflict with available measurements. Unfortunately, data on bubble growth rate in turbulent flows are quite limited because of the great difficulty of experiments. A set of experimental data regarding bubble growth in turbulent flow $(Re=1.8\times {10}^5\sim 1.9\times {10}^6)$ obtained by Kol’chugin et al. [74] and Lutovinov [75] was presented in [74]. The pressure of these experiments was in the range of $0.3\sim 4.0$ MPa, and the liquid superheat varied from $0.7$ to $2.5$ K. It was found that bubbles were always deformed and had an irregular shape. The bubble velocity and growth rate were measured by high-speed filming at different positions along the flow direction and the double-exposure method. The comparison between the bubble size predicted by the empirical correlations and the measured one is shown in Figure 13a. An under-prediction is given by the correlations of Whitaker [60], Feng & Michaelides [62] and Issa et al. [61]. In contrast, the Wolfert et al. [59] correlation given by Equation (30) accounts for the intensification due to flowing velocity and turbulence adequately by introducing an eddy conductivity. The prediction of bubble growth rate is in a good agreement with the measurement. Finally, the surface renewal theory is tested by combining the theoretical parts of conduction, convection and turbulence cumulatively, i.e.,

$$Nu=\left(\frac{12}{\pi }{Ja}_T+\frac{2}{\sqrt{\pi }}{Pe}^{1/2}+\frac{2d}{\sqrt{\pi a_l{\tau }_{turb}}}\right).$$

(37)

The results are shown in Figure 13b. The difference between the small and large eddy model lies in the time scale ${\tau }_{turb}$ for surface renewal as discussed above. It is evident that the prediction by the small eddy model is closer to the measurement. It may indicate that small eddies instead of large eddies are responsible for the surface renewal and interphase transfer. However, further data are required for a reliable evaluation of the correlations and theories for turbulent flows.

5. Conclusions

The interphase heat transfer in flashing flows is commonly believed to be a joint effect of heat conduction, convection and turbulence. The theory of bubble growth driven by each of the three mechanisms has been studied intensively. Whereas some general agreement exists on the contribution of conduction and convection, the effect of turbulence is still under debate. In spite of the large uncertainty, correlations obtained under stagnant or potential conditions are often extrapolated to turbulent cases. An evaluation of the applicability of these correlations in the two-fluid modelling of practical flashing flows is difficult, since there is a great deal of uncertainty related to other closures such as momentum interactions and turbulence modulation by bubbles. In the present work, heat transfer models are evaluated for the early stage of bubble growth in stagnant and high-velocity liquid with CFD simulations. The cases are ideal for the purpose of evaluation, since the effect of conduction, convection and turbulence can be investigated separately, and other uncertainties are reduced to a minimum. From the comparison with available DNS and experimental data, the following conclusions can be drawn:

• For creeping flows $({Re}_p\ll 1)$ the two correlations presented by Clift [64] provide excellent agreement for immobile (contaminated bubbles) and fully mobile (clean bubbles) interfaces, respectively.
• For high Reynolds number, the Ranz & Marshall [42] correlation reproduces well the transfer rate from (to) solid spheres, droplets and contaminated bubbles, while it gives under-predictions in the case of clean bubbles, for which the potential theory is more suitable at least for $Pe>10$.
• Stationary bubble growth follows the theoretical relation for heat conduction, i.e., $R\left(t\right)\propto t^{1/2}$. The numerical results are consistent with the analytical ones. The performance of the correlations is found to be dependent on the Jakob number. The correlation of Fritz & Ende [16] reproduces the bubble growth rate very well at low Jakob numbers, while those of Plesset & Zwick [17] and Forster & Zuber [13] give better predictions at moderate Jakob numbers. As the Jakob number increases further, the results of Olek et al. [18] get closer to the experimental data.
• For a reliable prediction of translating bubble growth, it is important to account for both heat conduction and convection. The conduction effect is evident in the initial stage even at moderate Jakob numbers. In cases with high Jakob numbers, the potential theory and the Ranz & Marshall [42] correlation under-predict the bubble size significantly, while the Wolfert [58] and Aleksandrov et al. [56] correlations, which account for both conduction and convection, deliver satisfying results.
• Wolfert et al. [59] is capable of reproducing the bubble growth rate in turbulent high-velocity flows by introducing an eddy conductivity, while significant under-prediction is given by other empirical correlations. The situation is improved by using the cumulative model proposed by Wolfert [58] supplemented with the surface renewal theory for turbulence. The time scale of small eddies is found to more suitable for the characterization of interfacial transfer than that of the large eddies. Nevertheless, acquisition of more detailed data is necessary for the quantitation of the turbulence effect.