Prediction of Bubble Departure in Forced Convection Boiling with a Mechanistic Model That Considers Dynamic Contact Angle and Base Expansion

: A mechanistic model for bubble dynamics in ﬂow boiling that is based on a force balance approach for a growing bubble is introduced. It considers the evaporation of the microlayer underneath the bubble, thermal di ﬀ usion and condensation around the bubble cap as well as dynamic inclination and contact angles between the bubble and the heating wall. It requires no recalibration of parameters to predict the bubble growth. Validation against di ﬀ erent experimental ﬂow boiling data was carried out with no case-dependent recalibration and yielded good agreement. The simulations conﬁrmed the dependency of bubble departure and lift-o ﬀ diameters on di ﬀ erent parameters such as heat ﬂux, liquid properties, subcooling temperature, system pressure, inclination angle of channel, mass


Introduction
Boiling is an essential heat transfer mechanism that needs an appropriate modelling for scientific and engineering analyses. It is a complex multi-physics problem that mass, momentum and energy transfer occur at the liquid-vapour interface as well heat conduction in the solid heating wall. Computational Fluid Dynamics (CFD) for industrial-scale problems typically makes use of the Eulerian-Eulerian two-fluid framework of interpenetrating continua [1,2]. In this approach, the mass, momentum and energy equations are separately solved for each phase and weighted by volume fraction representing the local ensemble averaged probability of each phase occurrence. The accompanied loss of information on interface structure requires extra models for mass, momentum and energy at the interface in form of source terms [3]. Similarly, for the modelling of wall boiling a model is needed that describes the heat transfer and phase change near the heating wall. The most prominent one is the so called Rensselaer Polytechnic Institute (RPI) heat partitioning model developed by Kurul and Podowski [4]. There, the heat flux from the heating wall is divided into contributions of liquid convection, quenching and evaporation as shown in Figure 1. Heat transfer is determined by a number of parameters, such as active nucleation site density, bubble departure diameter, bubble waiting time and bubble departure frequency that they are typically modelled with empirical correlations. For instance, the correlation of Tolubinsky and Kostanchuk has been widely used to predict the bubble departure diameter [5]. It gives the bubble departure diameter as a function of subcooling, that is, d dp = min{d ref with f being the bubble departure frequency [5]. These correlations were developed from pool boiling experiments at ambient pressure. Recently, Krepper et al. coupled the RPI wall boiling model with a population balance model and applied it to flow boiling with consideration of coalescence and breakup and also evaporation and condensation [1,6]. In their critical review of the existing correlations they concluded that some of the parameters used in the RPI wall boiling model cannot be applied for a broad range of different fluids or pressure levels for flow boiling [1]. Instead, they require some careful recalibrations for the intended applications. As it limits the general prediction capability of such models, the derivation of improved models from first principles should be a mandatory research objective for the near future.

Bubble Dynamics
The bubble dynamics in boiling has been investigated since the 1900s, starting with pool boiling. In 1917, Rayleigh investigated the bubble growth rate in the region controlled by inertia forces at low pressure and high Jakob numbers [7]. In 1954, Plesset and Zwick modelled the bubble growth in a uniform superheated liquid for heat diffusion-controlled growth in pool boiling [8]. The heat diffusion controlled growth model was later extended to non-uniform temperature fields by others, such as Forster and Zuber [9], Grift [10], Savic [11] and Zuber [12]. Mikic et al. considered both inertiacontrolled and thermal diffusion-controlled bubble growths and introduced a very general model of bubble dynamics for these periods [13]. From these analyses the following model concept was developed.
When the bubble is still small, its growth is very fast and determined by the inertia of the liquid being displaced because of the bubble formation. Hence, this period is referred to as inertia-controlled growth and during this period, a microlayer postulated and proven by Cooper and Loyd in 1969 is formed underneath the bubble [14] (see Figure 2). They derived a model of the microlayer based on the experimental data for pool boiling at low pressure [14]. After a while, the growth of bubble becomes slower and is no longer limited by the liquid displacement but by the heat flux supporting Heat transfer is determined by a number of parameters, such as active nucleation site density, bubble departure diameter, bubble waiting time and bubble departure frequency that they are typically modelled with empirical correlations. For instance, the correlation of Tolubinsky and Kostanchuk has been widely used to predict the bubble departure diameter [5]. It gives the bubble departure diameter as a function of subcooling, that is, with f being the bubble departure frequency [5]. These correlations were developed from pool boiling experiments at ambient pressure. Recently, Krepper et al. coupled the RPI wall boiling model with a population balance model and applied it to flow boiling with consideration of coalescence and breakup and also evaporation and condensation [1,6]. In their critical review of the existing correlations they concluded that some of the parameters used in the RPI wall boiling model cannot be applied for a broad range of different fluids or pressure levels for flow boiling [1]. Instead, they require some careful recalibrations for the intended applications. As it limits the general prediction capability of such models, the derivation of improved models from first principles should be a mandatory research objective for the near future.

Bubble Dynamics
The bubble dynamics in boiling has been investigated since the 1900s, starting with pool boiling. In 1917, Rayleigh investigated the bubble growth rate in the region controlled by inertia forces at low pressure and high Jakob numbers [7]. In 1954, Plesset and Zwick modelled the bubble growth in a uniform superheated liquid for heat diffusion-controlled growth in pool boiling [8]. The heat diffusion controlled growth model was later extended to non-uniform temperature fields by others, such as Forster and Zuber [9], Grift [10], Savic [11] and Zuber [12]. Mikic et al. considered both inertia-controlled and thermal diffusion-controlled bubble growths and introduced a very general model of bubble dynamics for these periods [13]. From these analyses the following model concept was developed.
When the bubble is still small, its growth is very fast and determined by the inertia of the liquid being displaced because of the bubble formation. Hence, this period is referred to as inertia-controlled growth and during this period, a microlayer postulated and proven by Cooper and Loyd in 1969 is formed underneath the bubble [14] (see Figure 2). They derived a model of the microlayer based on the experimental data for pool boiling at low pressure [14]. After a while, the growth of bubble becomes slower and is no longer limited by the liquid displacement but by the heat flux supporting the evaporation through the vapour-liquid interface. This period is referred to as thermal diffusion-controlled growth. The microlayer formed in the inertia-controlled period also contributes to Energies 2019, 12,1950 3 of 29 the growth of the bubble in this period as the microlayer extends during the expansion of bubble base diameter (see Figure 3). the evaporation through the vapour-liquid interface. This period is referred to as thermal diffusioncontrolled growth. The microlayer formed in the inertia-controlled period also contributes to the growth of the bubble in this period as the microlayer extends during the expansion of bubble base diameter (see Figure 3).  In 1970, Mikic et al. derived the following equation for a growing bubble in a uniformly superheated liquid, which is applicable for the entire range of inertia-controlled and diffusioncontrolled growth [13]: with for bubble growth on a surface. T w , T sat , h fg , ρ g and ρ l are heating wall temperature, saturation temperature, latent heat, vapour density and liquid density respectively. α l is thermal diffusivity of liquid and Ja is the Jakob number defines as Ja = ρ l C pl (T w − T sat ) ρ g h fg (6) with C pl being the liquid specific heat capacity.  the evaporation through the vapour-liquid interface. This period is referred to as thermal diffusioncontrolled growth. The microlayer formed in the inertia-controlled period also contributes to the growth of the bubble in this period as the microlayer extends during the expansion of bubble base diameter (see Figure 3).  In 1970, Mikic et al. derived the following equation for a growing bubble in a uniformly superheated liquid, which is applicable for the entire range of inertia-controlled and diffusioncontrolled growth [13]: with and for bubble growth on a surface. T w , T sat , h fg , ρ g and ρ l are heating wall temperature, saturation temperature, latent heat, vapour density and liquid density respectively. α l is thermal diffusivity of liquid and Ja is the Jakob number defines as Ja = ρ l C pl (T w − T sat ) ρ g h fg (6) with C pl being the liquid specific heat capacity. In 1970, Mikic et al. derived the following equation for a growing bubble in a uniformly superheated liquid, which is applicable for the entire range of inertia-controlled and diffusion-controlled growth [13]: and with b = 2 3 for bubble growth in an infinite medium and b = π 7 for bubble growth on a surface. T w , T sat , h fg , ρ g and ρ l are heating wall temperature, saturation temperature, latent heat, vapour density and liquid density respectively. α l is thermal diffusivity of liquid and Ja is the Jakob number defines as with C pl being the liquid specific heat capacity. In 1975, Labuntsov developed a model of bubble growth in the thermal diffusion-controlled period [15]. After that, in 1975, Van Stralen et al. determined the initial thickness of the evaporating microlayer beneath a hemispherical vapor bubble on a horizontal heating wall [16]. They also derived Energies 2019, 12,1950 4 of 29 a heat and mass diffusion-type solution accounting for the contribution of relaxation microlayer and evaporation microlayer to the bubble growth.
More recently, Klausner et al. developed a model based on the balance of forces acting on the growing bubble to predict the departure and lift-off and obtained a satisfactory predictive accuracy against their own data for the flow boiling of refrigerant R113 [17]. In their model they considered a constant bubble base diameter (contact diameter), recommended as d w = 0.09 mm, and constant advancing and receding contact angles, recommended as β ad = π 4 and β re = π 5 . In reality, however, these parameters are not constant. Later on, the modified versions of Klausner's model have been used by others such as Situ et al. [18], Sugrue [19], Thorncroft et al. [20] and Chen et al. [21] to predict their own experimental results but with a calibration of d w , β ad and β re . Klausner et al. [17] used the model of Mikic et al. [13] to account for the bubble growth, while Situ et al. [18] and most of the latter researchers employed Zuber's formulation [12]. The latter defined a parameter b that considers the bubble sphericity during bubble growth [12]. Later, this parameter has been used by many other groups as a "calibration parameter" with values from 0.24 to 2.4 to fit the model to their experimental data [22]. Yun et al. improved the model of Klausner et al. by considering a bubble condensation and evaluating the model for a wider range of pressure, temperature, and flow rates for water as well [23]. Many of the existing bubble dynamics models require empirical constants for different conditions. In the following, we will briefly introduce exemplarily three of more recent models.
Colombo and Fairweather derived the following equation for the bubble growth for flow boiling [22]: In this equation, Pr is the Prandtl number of the liquid, C 2 is microlayer constant (determined as C 2 = 1.78 [22]) and f sub is the portion of the bubble surface in contact with the subcooled liquid. Re b is the bubble Reynolds number and T l is the liquid temperature distribution in the superheated layer. Following Klausner et al. they set β ad = π 4 and β re = π 5 [17]. In addition, for the bubble base diameter they applied the suggestion of Yun et al., that is, d w = d b 15 [23]. Raj et al. [24] used Equation (7) for their bubble growth model for flow boiling, but with β ad = 89 • , β re = 9 • , d w = d b 15 and a bubble inclination angle of θ b = 10 • . Most recently, Mazzocco et al. also used this model with β ad = π 4 , β re = π 5 and d w = 0.09 mm [25]. In addition, they introduced and optimized a correlation factor to account for the effects of saturated and subcooled flows [25].
In Section 2 of this paper, a submodel for the bubble departure during flow boiling is introduced, which is based on a force balance for a growing bubble. In Section 2.3, the role of bubble contact and inclination angles and bubble bottleneck formation is discussed. In Section 2.4, the calculation of heat flux for validation of results against some experimental data is shown. In Section 3, the performance of this new bubble dynamics model against different experimental cases is assessed and in Section 3.4, the prediction of the impact of different parameters such as pressure, mass flow rate and subcooling on departure and lift-off diameters by this model is discussed. The key novelties of this model are that (1) it considers the contributions of microlayer, superheated layer surrounding the bubble and condensation to the bubble growth process, (2) it considers the dynamic contact and inclination angles as well as dynamic deformation, (3) it is physics-based and requires no recalibration of parameters at different conditions while it (4) can still be efficiently implemented in the Eulerian CFD approach.

Bubble Growth and Detachment Models
We divide the bubble growth into two periods. In the inertia-controlled (initial) growth a microlayer is formed underneath of the bubble (see Figure 3a) and the bubble diameter evolution is given by a Rayleigh-type equation [14] (see Section 2.1). It is followed by thermal diffusion-controlled growth following a Labuntsov-type solution [15] (see Figure 3b and Section 2.1). The evaporation of the microlayer is considered as an additional contribution of heat transfer in addition to the thermal diffusion from the surrounding bulk liquid. The depletion of microlayer causes a dryout on the heating wall and the influence of the dryout on the surface tension force at the heating wall will also be considered (see Section 2.2). During thermal diffusion-controlled growth also the microlayer is expanded because of the expansion of the gas-liquid interface on the heating wall.
The force acting on a growing bubble is taken from the work of Chen et al. [21] who developed it from Klausner's model [17] (see Section 2.2). Other than in previous force-balance bubble departure models, the bubble in our model will not immediately depart from the heating wall when the force in perpendicular direction of the wall is balanced, but instead the formation of a bottleneck is considered. The bottleneck connects the bubble's main body with the heating wall and enlarges until the bubble base diameter, r w , becomes zero or the length of bottleneck becomes larger than what the internal pressure difference allows. With the occurrence of any of these conditions, the bubble lifts-off. When the total height of the bubble (diameter of main body plus the bottleneck) becomes larger than the thermal layer thickness (Figure 4), condensation at the surface of bubble occurs, which is outside the thermal layer occurs and is calculated by the model of Ranz and Marshall [26] (see Figure 4 and Section 2.1).
microlayer is formed underneath of the bubble (see Figure 3a) and the bubble diameter evolution is given by a Rayleigh-type equation [14] (see Section 2.1). It is followed by thermal diffusion-controlled growth following a Labuntsov-type solution [15] (see Figure 3b and Section 2.1). The evaporation of the microlayer is considered as an additional contribution of heat transfer in addition to the thermal diffusion from the surrounding bulk liquid. The depletion of microlayer causes a dryout on the heating wall and the influence of the dryout on the surface tension force at the heating wall will also be considered (see Section 2.2). During thermal diffusion-controlled growth also the microlayer is expanded because of the expansion of the gas-liquid interface on the heating wall.
The force acting on a growing bubble is taken from the work of Chen et al. [21] who developed it from Klausner's model [17] (see Section 2.2). Other than in previous force-balance bubble departure models, the bubble in our model will not immediately depart from the heating wall when the force in perpendicular direction of the wall is balanced, but instead the formation of a bottleneck is considered. The bottleneck connects the bubble's main body with the heating wall and enlarges until the bubble base diameter, r w , becomes zero or the length of bottleneck becomes larger than what the internal pressure difference allows. With the occurrence of any of these conditions, the bubble liftsoff. When the total height of the bubble (diameter of main body plus the bottleneck) becomes larger than the thermal layer thickness (Figure 4), condensation at the surface of bubble occurs,which is outside the thermal layer occurs and is calculated by the model of Ranz and Marshall [26] (see Figure  4 and Section 2.1).
When the force in the tangential direction of the heating wall is balanced and the bubble inclination angle reaches its maximum value, the bubble starts sliding. This is called departure (see Section 2.3). Similar to the bottleneck formation that delays the bubble lift-off after the forces are balanced in perpendicular direction, the bubble inclination that comes from balanced forces in tangential direction can also delay the bubble sliding. To calculate the expansion of the bubble base diameter, a hypothesis based on the work of Thorncroft et al. is adopted [20]. According to their hypothesis, the relation between bubble radius and bubble base radius is (see section 2.3.1) When the force in the tangential direction of the heating wall is balanced and the bubble inclination angle reaches its maximum value, the bubble starts sliding. This is called departure (see Section 2.3). Similar to the bottleneck formation that delays the bubble lift-off after the forces are balanced in perpendicular direction, the bubble inclination that comes from balanced forces in tangential direction can also delay the bubble sliding.

Bubble Growth Rate
To calculate the expansion of the bubble base diameter, a hypothesis based on the work of Thorncroft et al. is adopted [20]. According to their hypothesis, the relation between bubble radius and bubble base radius is (see Section 2.3.1)

Bubble Growth Rate
As explained above, the bubble growth is divided into two periods: the inertia-controlled and the thermal diffusion-controlled growth. For bubble growth during flow boiling the correlations of these two periods that have been proposed by Ding et al. for pool boiling are used [27].
Bubble growth in inertia-controlled growth period defines as [13], Energies 2019, 12, 1950 6 of 29 In the model calculation, we need to consider the maximum bubble radius in inertia-controlled period. It is given by the model of Mikic et al. as [13], with A and B that are defined in Equation (5). Bubble growth rate in the thermal diffusion controlled growth period is defined as [27] with . V mi,g being the total volume of generated gas underneath of the bubble that by considering the mass balance between the generated gas and the evaporated liquid, The microlayer evaporation rate, dδ mi dt , is given by considering the heat balance between the latent heat going into the bubble and conductive heat transport through the microlayer with the thickness of δ mi as [27], Bubble growth rate due to macrolayer can be calculated by the model of Labuntsov [15] as with B 1 that defines as Bubble growth due to condensation on the bubble cap is calculated by the model of Ranz and Marshal [26] as According to Cooper and Loyd [14], the initial microlayer thickness is given as with C mi = 0.8 and C = C mi 2 = 0.64 × Pr. Thereby, t is the time in the inertia-controlled growth period, ν l is the kinematic viscosity of liquid and τ g is the maximal inertia-controlled growth time according to Zhao and Tsuruta [28]. The thickness of microlayer is [28] δ 0 mi (r) = Cα l ρ g h fg 2k l ∆T sup r (18) with r being the distance from the cavity, k l being the liquid conductivity, ∆T sup being superheating and C = 0.64 × Pr (see Figure 2). Recently, from the measurement of Utaka et al. it was found that the initial microlayer thickness has a linear slope to the bubble base diameter as [29], that is, Energies 2019, 12, 1950 7 of 29 However, this is only valid for water. When a bubble is in the thermal diffusion-controlled growth period, the microlayer will be extended with the expansion of the bubble contact base. Therefore, it is better to adapt the formula of Cooper and Loyd [14] (Equation (17)) with the data of Utaka et al. [29] with a new constant as C = 0.075 × Pr.
During inertia-controlled growth, the bubble has a semi-spherical shape. The conductive heat flux from the superheated heating wall through the microlayer is given by Based on the Equations (9)-(20) the total bubble radius growth rate is There, A b is the bubble surface area and f sub is the portion of the bubble surface that is in contact with subcooled liquid.
When the bubble radius reaches the maximum inertia-controlled bubble radius, r m,g , the thermal diffusion-controlled growth period starts. The bubble growth in this period is mainly caused by the evaporation of the microlayer underneath and superheated liquid around the bubble. The growth rate of the bubble during evaporation of the microlayer is determined by the consumption of liquid in the microlayer. In this work, the effect of surface tension and inertia effect on the pressure of the bubble and further the saturation temperature of vapour is neglected.
In order to simplify the model, we introduce the concept of a fixed saturation thermal layer thickness, δ th,sat between T w and T sat ( Figure 5). Condensation starts when the height of bubble is larger than the saturation thermal layer thickness. The total thermal layer thickness is defined as αt √ π for pool boiling [30] and δ tl δ th = Pr 1 3 for flow boiling [31]. Assuming a linear slope of the temperature in the thermal layer ( Figure 5) one gets The turbulent layer thickness δ tl , which is defined when u = 0.99 u max , can be calculated from the one-seventh power law for pipes, that is, Here, u max is the centreline velocity and R is the radius of the pipe. The turbulent layer thickness is also dependent on the distance to the inlet, that is, [31] The height of the bubble without bottleneck ( Figure 5) is (25) and the bubble growth velocity in the heating wall perpendicular direction is

Forces Acting on a Growing Bubble
We follow the force balance analysis of Klausner et al. [17] and Chen et al. [21] for a bubble growing on a superheated heating wall as shown in Figure 6. Based on the conservation of momentum the forces acting on the bubble in parallel (tangential: t) and perpendicular (normal: n) directions of the heating wall are balanced according to F total,n = F growth,n + F drag,n + F surf,n + F b,n + F cp,n + F sl,n F growth is the bubble growth force, F drag is the quasi-steady drag force due to the viscous fluid flow around the bubble, F cp is the contact pressure force due to the effect of the heating wall, F b is the buoyancy force, F sl is the shear lift force resulting from the asymmetrical flow distribution in the tangential direction of the heating wall, F surf is the surface tension force due to the contact surface

Forces Acting on a Growing Bubble
We follow the force balance analysis of Klausner et al. [17] and Chen et al. [21] for a bubble growing on a superheated heating wall as shown in Figure 6.

Forces Acting on a Growing Bubble
We follow the force balance analysis of Klausner et al. [17] and Chen et al. [21] for a bubble growing on a superheated heating wall as shown in Figure 6. Based on the conservation of momentum the forces acting on the bubble in parallel (tangential: t) and perpendicular (normal: n) directions of the heating wall are balanced according to F total,n = F growth,n + F drag,n + F surf,n + F b,n + F cp,n + F sl,n F growth is the bubble growth force, F drag is the quasi-steady drag force due to the viscous fluid flow around the bubble, F cp is the contact pressure force due to the effect of the heating wall, F b is the buoyancy force, F sl is the shear lift force resulting from the asymmetrical flow distribution in the tangential direction of the heating wall, F surf is the surface tension force due to the contact surface Based on the conservation of momentum the forces acting on the bubble in parallel (tangential: t) and perpendicular (normal: n) directions of the heating wall are balanced according to F total,n =F growth,n + F drag,n + F surf,n + F b,n + F cp,n + F sl,n F growth is the bubble growth force, F drag is the quasi-steady drag force due to the viscous fluid flow around the bubble, F cp is the contact pressure force due to the effect of the heating wall, F b is the buoyancy force, F sl is the shear lift force resulting from the asymmetrical flow distribution in the tangential direction of the heating wall, F surf is the surface tension force due to the contact surface on the heating wall and F growth,b is added mass force due to bubble growth in the bulk liquid velocity field. The forces and their derivations are introduces here. The bubble growth forces in normal and tangential direction of the heating wall are given by [21] F growth,n = −ρ l πr w 2 (r b ..
Due to the relative motion between bubble and liquid phase, the drag force on the bubble in the perpendicular direction can be evaluated as, with v b being the velocity of the bubble in perpendicular direction of the heating wall and C D is the drag force coefficient, which is dependent on turbulence intensity, bubble Reynolds number and bubble shape. Because of the pre-assumption of spherical bubble shape, C D is simplified with the correlation of Moore [30] and Clift et al. [32] as In the tangential direction, the drag force is calculated by same equation as perpendicular direction. The different is only that v b is replaced by (v sl − v f ) where v sl is the sliding velocity of bubble and v f is the velocity of the supposed stream line through the mass center of the bubble. v sl is calculated by the force applied on the bubble and v f is calculated by the 1/7th power low of the turbulent flow. The tangential direction of drag force is calculated by With v sl being the sliding velocity of bubble and v f is the velocity of the supposed streamline through the mass center of the bubble. v sl is calculated from the force acting on the bubble and v f is calculated from the 1/7th power law of turbulent flow.
According to the investigations of Klausner [17] and Chen et al. [21] the surface tension force in the perpendicular and tangential direction of heating wall can be expressed as, Due to the different heating wall orientation angle, θ w , the buoyancy can be defined in the perpendicular and tangential direction of heating wall as [21], Equation (38) gives the contact pressure force acting only perpendicular to the heating wall.
It is from the correlation of Thorncroft et al. [20]. r c = 5r b is the radius of curvature at the reference point according to Klausner et al. [17]) and σ is surface tension, which is a function of temperature. Equation (39) represents the shear lift force F sl that is partially due to Bernoulli suction and partially due to vortices in the approaching flow. It acts only in the perpendicular direction of the heating wall and is according to Van Helden et al. [33] and Chen et al. [21].
In this equation, U c is the velocity of the supposed streamline through the mass centre of the bubble, ∆U c = U c − U b , U b is the bubble sliding velocity in the mass center, C L1 is a constant of Bernoulli suction, which is equal to 11/8 according to Van Helden et al. [33], A is the bubble cross-section area perpendicular to the flow, V b is the bubble volume, C L2 is an empirical constant which has been recommended as 0.53 by Auton [34], and ω = U(n) n is the vorticity which is also effective only in perpendicular direction.
Chen et al. [21] and Thorncroft et al. [20] added an additional growth force named as growth force bulk F growth,b in the wall tangential direction for flow boiling. This force acts on the bubble due to bubble growth in a bulk liquid velocity field which was evaluated as [21],

Bubble Contact Angle (θ b ) and Inclination Angle (β)
Though the contact angle plays an important role in the calculation of forces, experimental data and reliable models are rather scarce. Klausner et al. recommended β ad = π 4 and β re = π 5 from their measurements for R113 [17]. Sugrue and Buongiorno derived β ad = 90.63 • and β re = 8.03 • from water experiments [35]. Mukherjee and Kandlikar proved and described that the contact angle does vary during the ebullition cycle and it depends on the liquid and vapour properties and heating wall material as well [36]. The contact angle does not have a constant value even in equilibrium conditions. However, it is possible to calculate the value based on the analysis of forces. Ding et al. investigated the bubble growth for horizontal pool boiling and considered the dynamic contact angle and dynamic bubble base diameter [27]. They used a mechanistic model based on the force balance method to predict the bubble size and the result of validation against the experimental data was good [27].
In this model the dynamic deformation of the bubble is considered. At the transition from the inertia-controlled growth to thermal diffusion controlled growth the contact angle is 90 • and the surface tension force is zero because the dryout radius, r d , is zero (Equations (34) and (35)). If in this period, the mass center velocity of the bubble is larger than the expansion speed of bubble radius the bubble will depart. However, in most cases this will not happen due to fast expansion of bubble. The dryout radius, r d , increases in this period when the sum of the negative forces that point toward the heating wall (mainly surface tension force, sometimes also growth force) is much higher than the sum of positive forces. The resulting net negative force will lead to a deformation of the bubble and determine the curvature and contact angle. This in turn reduces the surface tension force in the negative direction until the forces on the bubble are balanced. From the force calculation, the expected contact angle (see Figure 6) is expressed as The effect of the force in the tangential direction (which is influenced by roughness and wettability) on the contact angle is still not clear and more investigations are needed. In this work, we consider the force in perpendicular direction as the dominant one that effectively determines the contact angle. The factor 2 in Equation (41) means that the contact angle, β, is twice the microlayer contact angle, θ, which is used to calculate the surface tension force. β s is changing with time due to the continues changing of forces during the bubble growth.
The deformation of a bubble causes an expansion of the base radius r w so that its growth becomes different from the growth of the bubble radius, r b . Klausner et al. considered d w = 0.09 mm as being constant [17]. Later, Thorncroft et al. adopted d w = 2r b sin(β re ) with the receding side contact angle, β r , to improve the modelling accuracy [20]. As we consider in this work a dynamic contact angle we prefer to use a relation between the expansion rate of base radius, . r w , and that of bubble, . r b , instead of the fixed values of r w and r b . During the bubble growth, when the contact angle, β, is 90 • , the difference between the expansion of r w and r b should be largest. Later, with decreasing β, this difference should become smaller and smaller as the driving force becomes smaller and smaller.
Our approach is based on the work of Thorncroft et al. [20] and so we express the expansion rate of r w as During the deformation of the bubble the contact angle, β, reduces until it becomes equal to the expected contact angle, β s . After that, it will increase with the increase of β s and when the net force in perpendicular direction becomes positive, the bubble starts departing and a bottleneck is formed.
The inclination angle is also here considered as being dynamic during the bubble growth, depending, however, not on force but on torque (see Figure 7). The fulcrum of the torque is the center of the contact surface on the heating wall (point F in Figure 7). The expected inclination angle, θ b,s , is calculated as Γ b is the torque of buoyancy force that acts at the bubble mass center. Γ drag is the torque of drag force that acts at the bubble surface which contacts with the fluid and the part in contact with the heating wall. Γ growth , Γ cp and Γ sl are the torques of growth, contact pressure and shear lift forces that act at the mass center of the bubble (point O in Figure 7b). The torque is computed from the acting forces as The deformation of a bubble causes an expansion of the base radius r w so that its growth becomes different from the growth of the bubble radius, r b . Klausner et al. considered d w = 0.09 mm as being constant [17]. Later, Thorncroft et al. adopted d w = 2r b sin (β re ) with the receding side contact angle, β r , to improve the modelling accuracy [20]. As we consider in this work a dynamic contact angle we prefer to use a relation between the expansion rate of base radius, ṙw , and that of bubble, ṙb , instead of the fixed values of r w and r b . During the bubble growth, when the contact angle, β, is 90°, the difference between the expansion of r w and r b should be largest. Later, with decreasing β, this difference should become smaller and smaller as the driving force becomes smaller and smaller.
Our approach is based on the work of Thorncroft et al. [20] and so we express the expansion rate of r w as During the deformation of the bubble the contact angle, β, reduces until it becomes equal to the expected contact angle, β s . After that, it will increase with the increase of β s and when the net force in perpendicular direction becomes positive, the bubble starts departing and a bottleneck is formed.
The inclination angle is also here considered as being dynamic during the bubble growth, depending, however, not on force but on torque (see Figure 7). The fulcrum of the torque is the center of the contact surface on the heating wall (point F in Figure 7). The expected inclination angle, θ b,s , is calculated as Γ b is the torque of buoyancy force that acts at the bubble mass center. Γ drag is the torque of drag force that acts at the bubble surface which contacts with the fluid and the part in contact with the heating wall. Γ growth , Γ cp and Γ sl are the torques of growth, contact pressure and shear lift forces that act at the mass center of the bubble (point O in Figure 7b). The torque is computed from the acting forces as Γ ⃗ x = F ⃗⃗ x × r ⃗ x where x stands for the respective type of forces/torque (see Figure 7b), and r ⃗ is the radius vector, that is, r ⃗ = OF ⃗⃗⃗⃗⃗⃗ except for the tangential drag force/torque, for which r ⃗ = FH ⃗⃗⃗⃗⃗⃗ . In this work, the influence of roughness is neglected and β w is the liquid contact angle at the equilibrium state [37] (Figure 8). Both the vertical and horizontal forces are balanced at the contact point, which is known as equilibrium. The horizontal component of the gas-liquid surface tension force is balanced by the adhesive force and the vertical one is balanced by the difference of the forces along the heating wall as In this work, the influence of roughness is neglected and β w is the liquid contact angle at the equilibrium state [37] (Figure 8). Both the vertical and horizontal forces are balanced at the contact point, which is known as equilibrium. The horizontal component of the gas-liquid surface tension force is balanced by the adhesive force and the vertical one is balanced by the difference of the forces along the heating wall as γ gl cos β w = γ sg − γ sl (44) with γ gl , γ sg and γ sl being the interphase gas-liquid, solid-gas and solid-liquid surface tension. with γ gl , γ sg and γ sl being the interphase gas-liquid, solid-gas and solid-liquid surface tension. The formation of bubble inclination angle, θ b , is due to the drag force acting on the top of the bubble. θ b is calculated using Newton's law with the movement of the top part of the bubble as with L top being the displacement of the bubble top side with respect to its horizontal position and h bubble is the height of the bubble (see Figure 7a).
The bubble will not leave the heating wall until θ b ≥ θ b,s , even when the force in the tangential direction becomes positive. The calculation of advancing and receding contact angle in this work is based on the contact angle, β , for horizontal pool boiling and an additional bubble inclination for flow boiling. As during the formation of the inclination angle the bubble shape in the model is hemispherical or perfectly spherical, the impact of buoyancy, F b,n , is neglected when F b,n < 0 , β ad = β and β re = β − θ b which is limited by β ad ≤ π 2 and β re ≥ 0.
The contact angle changes due to the force balance. When the force in perpendicular direction, F total,n , becomes positive, the contact angle becomes π 2 at both sides for horizontal pool boiling, while during flow boiling this is only at the advancing side. The reason is the dominance of the force in perpendicular direction at both sides of bubble in horizontal pool boiling while for flow boiling it is only at the advancing side. At the same time, the bottleneck is formed due to the evaporation of the microlayer and a bubble stretch by adhesion force on the heating wall. For flow boiling, the bubble contact angle at the advancing side, β ad , can also reach π 2 when the total force along the flow direction is large enough, that is, Because of the difference between β ad and β re during the bubble growth, the bubble has different base diameter expansion velocities at both sides of the bubble based on the Equation (42). The mass center of the bubble changes or the bubble can be considered as sliding already due to the fast expansion of the receding side and slow expansion of the advancing side, though the bubble still expands at both sides

Bottleneck
In our new model we try to account the dynamic bubble shape as close to the observed physics as possible in order to avoid the case-dependent parameter tuning. We start with a hemispherical bubble shape during inertia-controlled growth, then change from hemispherical to spherical shape during the thermal diffusion-controlled period, then continue with an inclined bubble to end in a The formation of bubble inclination angle, θ b , is due to the drag force acting on the top of the bubble. θ b is calculated using Newton's law with the movement of the top part of the bubble as with L top being the displacement of the bubble top side with respect to its horizontal position and h bubble is the height of the bubble (see Figure 7a). The bubble will not leave the heating wall until θ b ≥ θ b,s , even when the force in the tangential direction becomes positive. The calculation of advancing and receding contact angle in this work is based on the contact angle, β, for horizontal pool boiling and an additional bubble inclination for flow boiling. As during the formation of the inclination angle the bubble shape in the model is hemispherical or perfectly spherical, the impact of buoyancy, F b,n , is neglected when F b,n < 0, β ad = β and β re = β − θ b which is limited by β ad ≤ π 2 and β re ≥ 0. The contact angle changes due to the force balance. When the force in perpendicular direction, F total,n , becomes positive, the contact angle becomes π 2 at both sides for horizontal pool boiling, while during flow boiling this is only at the advancing side. The reason is the dominance of the force in perpendicular direction at both sides of bubble in horizontal pool boiling while for flow boiling it is only at the advancing side. At the same time, the bottleneck is formed due to the evaporation of the microlayer and a bubble stretch by adhesion force on the heating wall. For flow boiling, the bubble contact angle at the advancing side, β ad , can also reach π 2 when the total force along the flow direction is large enough, that is, Because of the difference between β ad and β re during the bubble growth, the bubble has different base diameter expansion velocities at both sides of the bubble based on the Equation (42). The mass center of the bubble changes or the bubble can be considered as sliding already due to the fast expansion of the receding side and slow expansion of the advancing side, though the bubble still expands at both sides.

Bottleneck
In our new model we try to account the dynamic bubble shape as close to the observed physics as possible in order to avoid the case-dependent parameter tuning. We start with a hemispherical bubble shape during inertia-controlled growth, then change from hemispherical to spherical shape during the thermal diffusion-controlled period, then continue with an inclined bubble to end in a sphere plus a bottleneck as illustrated in Figure 9. Eventually, after lift-off the bubble becomes spherical. The force balance in perpendicular direction directly leads to the bottleneck formation. At this moment, the bubble can still be in contact with the heating wall if the bubble main body moving velocity, v b , is less than the bubble expansion velocity due to the bubble volume increase coming from microlayer evaporation. Only when the bubble moves faster in perpendicular direction then it expands and a bottleneck is formed. The base diameter of bubble (now it is that of the bottleneck) starts to shrink when the evaporation of the microlayer becomes less and less. Unlike in the conventional force analysis model that the bubble lifts off when the forces in perpendicular direction are balanced, here, the criteria of lift-off is the bottleneck breakage or the base diameter shrinks to zero. sphere plus a bottleneck as illustrated in Figure 9. Eventually, after lift-off the bubble becomes spherical. The force balance in perpendicular direction directly leads to the bottleneck formation. At this moment, the bubble can still be in contact with the heating wall if the bubble main body moving velocity, v b , is less than the bubble expansion velocity due to the bubble volume increase coming from microlayer evaporation. Only when the bubble moves faster in perpendicular direction then it expands and a bottleneck is formed. The base diameter of bubble (now it is that of the bottleneck) starts to shrink when the evaporation of the microlayer becomes less and less. Unlike in the conventional force analysis model that the bubble lifts off when the forces in perpendicular direction are balanced, here, the criteria of lift-off is the bottleneck breakage or the base diameter shrinks to zero. The bubble base radius, r w , will shrink when v b πr w 2 > V̇m i,g . Due to volume conservation ( Figure   10) this gives. The bottleneck height,h bt , is calculated from the bubble mass center velocity, v b , and the time difference between the time when force in perpendicular direction becomes positive, t fp , and the time from this moment , t: The bubble base radius, r w , will shrink when v b πr w 2 > . V mi,g . Due to volume conservation ( Figure 10)  sphere plus a bottleneck as illustrated in Figure 9. Eventually, after lift-off the bubble becomes spherical. The force balance in perpendicular direction directly leads to the bottleneck formation. At this moment, the bubble can still be in contact with the heating wall if the bubble main body moving velocity, v b , is less than the bubble expansion velocity due to the bubble volume increase coming from microlayer evaporation. Only when the bubble moves faster in perpendicular direction then it expands and a bottleneck is formed. The base diameter of bubble (now it is that of the bottleneck) starts to shrink when the evaporation of the microlayer becomes less and less. Unlike in the conventional force analysis model that the bubble lifts off when the forces in perpendicular direction are balanced, here, the criteria of lift-off is the bottleneck breakage or the base diameter shrinks to zero. The bubble base radius, r w , will shrink when v b πr w 2 > V̇m i,g . Due to volume conservation ( Figure   10) this gives. The bottleneck height,h bt , is calculated from the bubble mass center velocity, v b , and the time difference between the time when force in perpendicular direction becomes positive, t fp , and the time from this moment , t: The bottleneck height, h bt , is calculated from the bubble mass center velocity, v b , and the time difference between the time when force in perpendicular direction becomes positive, t fp , and the time from this moment, t: Ding et al. claimed that in pool boiling when the microlayer underneath the bubble is completely consumed or the pressure difference along the bottleneck reaches its limit, the bottleneck will break [27].
Here also, we apply this methodology for bottleneck breakage. The pressure inside the bubble on the heating wall, P B , differs from the required pressure to hold up the bubble on the heating wall, P B , (Figure 9). P B is strongly dependent on the base radius, r w , and decreases with its shrinkage. The pressure difference has been estimated by Ding et al. [27] as When ∆P B B is larger than the total force in perpendicular direction acting on the bubble base area, that is, the bubble will return back to spherical shape and depart from the heating wall. Of course, if the base radius shrinks to zero before ∆P B B reaches the limit, the bubble will also depart. The complete bubble growth and departure model chart for subcooled flow boiling is summarized in Figure 11. Ding et al. claimed that in pool boiling when the microlayer underneath the bubble is completely consumed or the pressure difference along the bottleneck reaches its limit, the bottleneck will break [27]. Here also, we apply this methodology for bottleneck breakage. The pressure inside the bubble on the heating wall, P B , differs from the required pressure to hold up the bubble on the heating wall, P B ′ , (Figure 9). P B is strongly dependent on the base radius, r w , and decreases with its shrinkage. The pressure difference has been estimated by Ding et al. [27] as When ∆P B ′ B is larger than the total force in perpendicular direction acting on the bubble base area, that is, the bubble will return back to spherical shape and depart from the heating wall. Of course, if the base radius shrinks to zero before ∆P B ′ B reaches the limit, the bubble will also depart. The complete bubble growth and departure model chart for subcooled flow boiling is summarized in Figure 11.

Heat Transfer from Heating Wall to Liquid
In principle, it is not necessary to calculate the heat flux in this model as it will be calculated in the CFD simulation directly. However, for validation against experimental data, for which only heat fluxes but no superheated wall temperatures are available, it is required. As introduced, the bubble growth rate is strongly dependent on the wall superheat. Therefore, it is necessary to calculate the wall superheat under different heat fluxes. Similar to the heat partitioning model the heat flux is divided into different terms: evaporation of microlayer (Figure 12a), evaporation of macrolayer Figure 12b), heat transfer from the heating wall to gas in the dryout region (Figure 12c), quenching ( Figure 12d) and forced convection (Figure 12e). Quantitatively, this is given as Figure 11. Schematics of the model including submodels for bubble growth during flow boiling.

Heat Transfer from Heating Wall to Liquid
In principle, it is not necessary to calculate the heat flux in this model as it will be calculated in the CFD simulation directly. However, for validation against experimental data, for which only heat fluxes but no superheated wall temperatures are available, it is required. As introduced, the bubble growth rate is strongly dependent on the wall superheat. Therefore, it is necessary to calculate the wall superheat under different heat fluxes. Similar to the heat partitioning model the heat flux is divided into different terms: evaporation of microlayer (Figure 12a (Figure 12e). Quantitatively, this is given as Here, . m mi and . m ma are mass rates of evaporated liquid in the microlayer and macrolayer, k g and α g are conductivity and thermal diffusivity of vapor, τ d and τ q are the moments of the start of dryout and quenching, t d is time of departure and T b is the bulk liquid temperature. The quenching heat flux is from Zhao and Tsuruta [28] and h fc is from the Dittus-Boelter equation [31] h fc = 0.023 Re 0.8 l Pr 0.4 Re l is the Reynolds number of liquid flow that is calculated based on the bulk flow velocity and Pr l is the Prandtl number of the liquid. D h is the hydrodynamic diameter of channel which is simply the diameter of the tube for a circular tube.
Here, ṁm i and ṁm a are mass rates of evaporated liquid in the microlayer and macrolayer, k g and α g are conductivity and thermal diffusivity of vapor, τ d and τ q are the moments of the start of dryout and quenching, t d is time of departure and T b is the bulk liquid temperature. The quenching heat flux is from Zhao and Tsuruta [28] and h fc is from the Dittus-Boelter equation [31] h fc = 0.023 Re l 0.8 Pr l 0.4 ( K l D h ) (52) Re l is the Reynolds number of liquid flow that is calculated based on the bulk flow velocity and Pr l is the Prandtl number of the liquid. D h is the hydrodynamic diameter of channel which is simply the diameter of the tube for a circular tube.

Heat Transfer in the Heating Wall
The impact of heating wall thickness and material on the heat transfer in the boiling process has usually not been considered in the previous studies. However, the heating wall can act as a thermal buffer and affects the heat transfer near hot spots caused by the dryout under the bubble. The energy equation yields Q̇i n dt + Q̇o ut dt + Q̇n ,w dt = ρ w V w C p,w dT w . (53)

Heat Transfer in the Heating Wall
The impact of heating wall thickness and material on the heat transfer in the boiling process has usually not been considered in the previous studies. However, the heating wall can act as a thermal buffer and affects the heat transfer near hot spots caused by the dryout under the bubble. The energy equation yields . Q in dt + . Q out dt + . Q n,w dt = ρ w V w C p,w dT w .
(53) The temperature change in the heating wall is given as with the specific heat capacity, C pw , and the density, ρ w , of the heating wall. For the heating wall with the thickness of δ w and unit depth (1 m), the volume of one cell, V w , yields In Equation (54), Q. in , . Q out and . Q n,w are the heat rate into the heating wall, from the heating wall and conduction heat rate between two neighboring segments of wall as shown in Figure 13 and calculated as . .
. Q n,w = q n,w ·δ w (58) q in is the heat flux into the heating wall, q out comes from Equation (51) and the tangential heat flux in the heating wall, q n,w , is given by with k w being the thermal conductivity of heating wall, ∆T w is the temperature difference between two neighboring heating wall segments ( Figure 13). The temperature change in the heating wall is given as dT w dt = Q̇i n + Q̇o ut + Q̇n ,w ρ w V w C pw (54) with the specific heat capacity, C pw , and the density, ρ w , of the heating wall. For the heating wall with the thickness of δ w and unit depth (1 m), the volume of one cell, V w , yields In Equation (54), Q̇i n , Q̇o ut and Q̇n ,w are the heat rate into the heating wall, from the heating wall and conduction heat rate between two neighboring segments of wall as shown in Figure 13 and calculated as Q̇i n = q ′′ in . ∆L w (56) Q̇o ut = q ′′ out . ∆L w (57) q ′′ in is the heat flux into the heating wall, q ′′ out comes from Equation (51) and the tangential heat flux in the heating wall, q ′′ n,w , is given by with k w being the thermal conductivity of heating wall, ∆T w is the temperature difference between two neighboring heating wall segments ( Figure 13).

Discretization Dependency Study
The wall boiling sub-model requires a proper time discretization for bubble dynamics and a spatial discretization for the microlayer and the heat transfer inside the wall. Both time and space are discretized using a central differences scheme. In order to do the mesh independency, some cases with time steps from 1 µ s to 30 µ s and mesh sizes from 10 µ m to 50 µ m are tested. We found that the mesh size less than 50 µ m and a time step less than 30 µ s was found to yield a stable solution and the

Discretization Dependency Study
The wall boiling sub-model requires a proper time discretization for bubble dynamics and a spatial discretization for the microlayer and the heat transfer inside the wall. Both time and space are discretized using a central differences scheme. In order to do the mesh independency, some cases with time steps from 1 µs to 30 µs and mesh sizes from 10 µm to 50 µm are tested. We found that the mesh size less than 50 µm and a time step less than 30 µs was found to yield a stable solution and the model converged well. In addition, the deviation from the average value of bubble departure diameter predictions from these cases are less than 2%. The Courant-Friedrichs-Lewy (CFL) number is controlled to be less than 1. More information about the spatial and time discretization are available in [27].

Experimental Database
For validation, we used published data from different experiments. One is from Klausner et al. [17] who operated an experimental facility with a horizontal 25 × 25 mm square cross section for flow boiling with refrigerant R113. They measured the mean bubble departure diameter at different mass flux, heat fluxes and wall superheat. In their experiments, there was a stratified two-phase flow in the boiling section.
Prodanovic et al. [38] studied boiling in a vertical annular test section (quartz glass, 22-mm inner diameter) with upward water flow. Boiling was on a heated hollow stainless-steel tube with 12.7-mm outer diameter welded to solid copper rods. The heated rod was 0.22-m downstream of the inlet to get a fully developed turbulent flow. The authors monitored the bubble behaviour with a high-speed camera at 6000 to 8000 frames per second located 0.44-m downstream of the start of the heated section. They captured the bubble behaviour from inception to collapse and obtained bubble shapes during lifetime, times of detachment from the heating wall and typical bubble size.
Situ et al. [18] conducted different studies for forced convective subcooled flow boiling in vertical upward water flow in an annular channel of 38.1-mm inner diameter with a 19.1-mm diameter heater rod. They also used a high-speed camera to assess bubble lift-off diameters, growth rate, and velocity after lift-off at different mass fluxes between 466 to 900 kg m −2 s −1 and heat fluxes from 54 to 206 kWm −2 at ambient pressure.
Sugrue investigated the effect of the heating wall orientation angle on bubble departure diameter for subcooled boiling at different mass flux, pressures and heat fluxes with a high-speed camera [19]. The section was a rectangular stainless-steel channel with 1.43-cm width and 1.99-cm height and a heater on one side. The author carefully controlled the Froude number (based on the channel hydraulic diameter) between 0.42 to 1.06 and Reynolds number between 11,800 and 34,500 to eliminate the impact of buoyancy and wall orientation angle. Moreover, the author controlled the subcooled flow and heat flux to ensure an "isolated bubbles" regime for which interactions between bubbles at adjacent nucleation sites can be neglected.
All the experiments introduced above cover nearly all physical parameters, which can impact the bubble behaviour. The present model has been validated for the departure or lift-off diameter with these data as shown in Table 1.

Validation of Bubble Growth and Departure (Lift-Off) Models
For flow boiling, the departure diameter is not equal to the lift-off diameter when the bubble slides on the heating wall. Departure means that the bubble leaves the originating cavity in tangential direction while lift-off means that the bubble leaves the wall. We compare the calculated departure diameter against experimental data. There are two experiments dealing with departure diameter: Klausner et al. [17] and Sugrue [19]. Figure 14 shows the comparison with the data from Klausner et al. [17] for R113 as the working fluid. The average deviation is 21%, which can be considered as a good agreement. Because of the open channel in the experimental setup, fully developed turbulent flow is considered in the model calculation. Moreover, the wettability of R113 has also been taken into account by considering the contact angle β w =0.174 rad [39].

Validation of Bubble Growth and Departure (lift-off) Models
For flow boiling, the departure diameter is not equal to the lift-off diameter when the bubble slides on the heating wall. Departure means that the bubble leaves the originating cavity in tangential direction while lift-off means that the bubble leaves the wall. We compare the calculated departure diameter against experimental data. There are two experiments dealing with departure diameter: Klausner et al. [17] and Sugrue [19]. Figure 14 shows the comparison with the data from Klausner et al. [17] for R113 as the working fluid. The average deviation is 21%, which can be considered as a good agreement. Because of the open channel in the experimental setup, fully developed turbulent flow is considered in the model calculation. Moreover, the wettability of R113 has also been taken into account by considering the contact angle β w =0.174 rad [39]. An even better agreement for bubble departure diameter we get for the data of Sugure [19] ( Figure 15). It is about 12.44%. These experiments [19] cover a wide range of parameters: wall orientation angle, mass flux, heat flux and subcooling temperature. The model of Klausner et al. [17] has a reasonable agreement with their own experimental data for R113. However, it predicts the departure diameter for Sugure's experiments [19] with water only with an average deviation of 65.5 %. The modified model version that was introduced for water by Yun et al. [23] shows a better agreement but still has an average deviation of 37.9 %. Both models can predict the experimentally observed trends but they overpredict the bubble departure diameter (see Figure 16). Only at the pressure of 5 bar the model shows a small overprediction. An even better agreement for bubble departure diameter we get for the data of Sugure [19] ( Figure 15). It is about 12.44%. These experiments [19] cover a wide range of parameters: wall orientation angle, mass flux, heat flux and subcooling temperature. The model of Klausner et al. [17] has a reasonable agreement with their own experimental data for R113. However, it predicts the departure diameter for Sugure's experiments [19] with water only with an average deviation of 65.5%. The modified model version that was introduced for water by Yun et al. [23] shows a better agreement but still has an average deviation of 37.9%. Both models can predict the experimentally observed trends but they overpredict the bubble departure diameter (see Figure 16). Only at the pressure of 5 bar the model shows a small overprediction.   The dynamic contact and inclination angles are shown in Figure 17. The advancing contact angle β ad strongly decreases from π 2 in the hemisphere to ~0.65 rad due to bubble growth and then increases back to π 2 due to the force acting on the bubble in the tangential direction. The receding contact angle, β re , firstly decreases in a similar way as at the advancing side due to bubble growth. Then it increases due to the force acting in the perpendicular direction and also the bubble inclination. Further, it decreases again due to the force acting on the bubble in the tangential direction resulting in the bubble inclination. Finally, it increases again to π 2 due to the force applied on the bubble in the The dynamic contact and inclination angles are shown in Figure 17. The advancing contact angle β ad strongly decreases from π 2 in the hemisphere to~0.65 rad due to bubble growth and then increases back to π 2 due to the force acting on the bubble in the tangential direction. The receding contact angle, β re , firstly decreases in a similar way as at the advancing side due to bubble growth. Then it increases due to the force acting in the perpendicular direction and also the bubble inclination. Further, it decreases again due to the force acting on the bubble in the tangential direction resulting in the bubble inclination. Finally, it increases again to π 2 due to the force applied on the bubble in the perpendicular direction. The bubble inclination angle, θ b , increases from zero when the bubble is generated as a symmetric hemisphere until 1.43 rad when the forces acting on the bubble are still not balanced in the tangential direction. Then it decreases slightly due to the base diameter expansion. Finally, it decreases back to zero due to bubble sliding, speed up and the decreased drag force resulting from the relative velocity between bubble and liquid. Since the dynamic contact and inclination angles are obtained from the model there is no need to predetermine them like the works done by Klausner et al. [17], Hong et al. [40], Colombo and Fairweather [22] and Raj et al. [24]. generated as a symmetric hemisphere until 1.43 rad when the forces acting on the bubble are still not balanced in the tangential direction. Then it decreases slightly due to the base diameter expansion. Finally, it decreases back to zero due to bubble sliding, speed up and the decreased drag force resulting from the relative velocity between bubble and liquid. Since the dynamic contact and inclination angles are obtained from the model there is no need to predetermine them like the works done by Klausner et al. [17], Hong et al. [40], Colombo and Fairweather [22] and Raj et al. [24]. In 2013, Krepper et al. [1] assessed the application of Lemmert and Chawla correlations [41] (see Equation (60) to model the nucleation site density in a CFD simulation. The nucleation site density is too complex to predict as it is strongly dependent on the heating wall material properties. However, with some recalibration it can be described by a simple correlation based on wall superheat: Here N ref is a referenced nucleation site density (to be recalibrated). ∆T refN = 10 K and P = 1.805 were found to yield satisfactory results in the overall model [1]. This correlation has been applied in the present model to calculate the heat flux based on the known wall superheat for the data of Klausner et al. [17]. The recalibration gives N ref = 7.0E+4 m −2 and comparison of calculated heat flux and measured heat flux is shown in Figure 18. The average deviation is 30%, which proves the applicability of the Lemmert and Chawla correlation [41] and the assessment of Krepper et al. [1]. In 2013, Krepper et al. [1] assessed the application of Lemmert and Chawla correlations [41] (see Equation (60) to model the nucleation site density in a CFD simulation. The nucleation site density is too complex to predict as it is strongly dependent on the heating wall material properties. However, with some recalibration it can be described by a simple correlation based on wall superheat: Here N ref is a referenced nucleation site density (to be recalibrated). ∆T refN = 10 K and P = 1.805 were found to yield satisfactory results in the overall model [1]. This correlation has been applied in the present model to calculate the heat flux based on the known wall superheat for the data of Klausner et al. [17]. The recalibration gives N ref = 7.0E + 4 m −2 and comparison of calculated heat flux and measured heat flux is shown in Figure 18. The average deviation is 30%, which proves the applicability of the Lemmert and Chawla correlation [41] and the assessment of Krepper et al. [1]. generated as a symmetric hemisphere until 1.43 rad when the forces acting on the bubble are still not balanced in the tangential direction. Then it decreases slightly due to the base diameter expansion. Finally, it decreases back to zero due to bubble sliding, speed up and the decreased drag force resulting from the relative velocity between bubble and liquid. Since the dynamic contact and inclination angles are obtained from the model there is no need to predetermine them like the works done by Klausner et al. [17], Hong et al. [40], Colombo and Fairweather [22] and Raj et al. [24]. In 2013, Krepper et al. [1] assessed the application of Lemmert and Chawla correlations [41] (see Equation (60) to model the nucleation site density in a CFD simulation. The nucleation site density is too complex to predict as it is strongly dependent on the heating wall material properties. However, with some recalibration it can be described by a simple correlation based on wall superheat: Here N ref is a referenced nucleation site density (to be recalibrated). ∆T refN = 10 K and P = 1.805 were found to yield satisfactory results in the overall model [1]. This correlation has been applied in the present model to calculate the heat flux based on the known wall superheat for the data of Klausner et al. [17]. The recalibration gives N ref = 7.0E+4 m −2 and comparison of calculated heat flux and measured heat flux is shown in Figure 18. The average deviation is 30%, which proves the applicability of the Lemmert and Chawla correlation [41] and the assessment of Krepper et al. [1]. and departure diameters is shown in Figures 19 and 20. The average deviation from the case of Situ et al. [18] for more than 90 data points is 35%. The location of the observation point (where the flow was observed a high-speed camera) affects the data via the thickness of turbulent layer on the heating wall [18]. For the experiments of Situ et al. [18] the location of observation points were specified which we found helpful for analysis and interpretation. Both Situ et al. [18] and Prodanovic et al. [38] did not measure the local superheat in their experiments. However, for the present model the superheating is a necessary parameter for calculation of nucleation site density. The recalibration of N ref is 9.0E+4 m −2 for the case of Situ et al.'s [18] and 4.8E+4 m −2 for the case of Prodanovic et al. [38] is considered. The comparison for lift off and departure diameters is shown in Figure 19 and Figure 20. The average deviation from the case of Situ et al. [18] for more than 90 data points is 35%. The location of the observation point (where the flow was observed a high-speed camera) affects the data via the thickness of turbulent layer on the heating wall [18]. For the experiments of Situ et al. [18] the location of observation points were specified which we found helpful for analysis and interpretation.  The highest deviation of 37% on average is found for the experiment of Prodanovic et al. [38]. In particular at low velocity (u in = 0.08 m/s) the model strongly overpredicts the lift-off diameter which is up to 3 times of the experimental value. Similar findings were reported by Colombo and Fairweather [22] and Sugrue and Buongiorno [35]. Obviously, the data of Prodanovic et al. [38] is the most difficult data to predict. A possible reason may be the high subcooling that can lead to cases of bubble collapse or merging [22]. Another reason maybe is the impact of the location of observation point which is not considered in this work. Both Situ et al. [18] and Prodanovic et al. [38] did not measure the local superheat in their experiments. However, for the present model the superheating is a necessary parameter for calculation of nucleation site density. The recalibration of N ref is 9.0E+4 m −2 for the case of Situ et al.'s [18] and 4.8E+4 m −2 for the case of Prodanovic et al. [38] is considered. The comparison for lift off and departure diameters is shown in Figure 19 and Figure 20. The average deviation from the case of Situ et al. [18] for more than 90 data points is 35%. The location of the observation point (where the flow was observed a high-speed camera) affects the data via the thickness of turbulent layer on the heating wall [18]. For the experiments of Situ et al. [18] the location of observation points were specified which we found helpful for analysis and interpretation.  The highest deviation of 37% on average is found for the experiment of Prodanovic et al. [38]. In particular at low velocity (u in = 0.08 m/s) the model strongly overpredicts the lift-off diameter which is up to 3 times of the experimental value. Similar findings were reported by Colombo and Fairweather [22] and Sugrue and Buongiorno [35]. Obviously, the data of Prodanovic et al. [38] is the most difficult data to predict. A possible reason may be the high subcooling that can lead to cases of bubble collapse or merging [22]. Another reason maybe is the impact of the location of observation point which is not considered in this work. The highest deviation of 37% on average is found for the experiment of Prodanovic et al. [38]. In particular at low velocity (u in = 0.08 m/s) the model strongly overpredicts the lift-off diameter which is up to 3 times of the experimental value. Similar findings were reported by Colombo and Fairweather [22] and Sugrue and Buongiorno [35]. Obviously, the data of Prodanovic et al. [38] is the most difficult data to predict. A possible reason may be the high subcooling that can lead to cases of bubble collapse or merging [22]. Another reason maybe is the impact of the location of observation point which is not considered in this work.

Sensitivities Analysis
In this section, the effect of different parameters such as mass flux, heat flux, subcooling temperature, geometry, heating wall orientation angle and pressure on departure and lift-off diameters for flow boiling are being analysed with our model. An increase of the mass flux (inlet velocity) usually causes a decrease of the bubble departure diameter (see Figure 21a) which is also agreed by Sugrue [19], Klausner et al. [17] and Yun et al. [23] as shown in Figure 16b. The reason is that drag and lift forces increase with mass flux, which decreases lift-off diameter. This trend is confirmed by the model and shown in Figure 22a.

Sensitivities Analysis
In this section, the effect of different parameters such as mass flux, heat flux, subcooling temperature, geometry, heating wall orientation angle and pressure on departure and lift-off diameters for flow boiling are being analysed with our model. An increase of the mass flux (inlet velocity) usually causes a decrease of the bubble departure diameter (see Figure 21a) which is also agreed by Sugrue [19], Klausner et al. [17] and Yun et al. [23] as shown in Figure 16b. The reason is that drag and lift forces increase with mass flux, which decreases lift-off diameter. This trend is confirmed by the model and shown in Figure 22a. Next, the heat flux from the heating wall determines the growth speed of the bubble. A higher heat flux leads to a faster growth of the bubble, which impacts the growth force that opposes the bubble departure. That is, a high heat flux leads to larger bubble departure diameter. In addition, for flow boiling the pressure on the distorted bubble is reduced in the direction facing the flow. This leads to a delay of the bubble lift-off, which has been confirmed by the model calculations. The increase is even relatively small (see Figure 21b). However, Sugrue's experiment only confirmed this effect at a wall orientation angle of 90° [19]. At an orientation angle of 45° it is even the other way round when heat flux increases from 50-100 kWm -2 (see Figure 21b). When the heat flux increases to 500 kWm -2 the departure diameter increases sharply in the model calculation for the two orientation angles. The same effect was found in the modelling results based on the setup of Situ et al. [18] shown in Figure 22a. That is, also for a wider range of heat flux. The calculation is for an axial location of 0.57 m from the inlet. At low heat flux, i.e. q ′′ ≤ 60 kWm -2 , the bubble even does not leave the heating wall. When the heat flux is between 60 kWm -2 and 120 kWm -2 , the bubble lift-off diameter decreases slightly from 0.2 mm to 0.17 mm. Afterwards, with an increase of the heat flux the bubble lift-off diameter rapidly increases to 0.54 mm at 180 kWm -2 . A credible explanation for these results comes from the dynamic contact angle and dynamic inclination angle. The dynamic inclination angle is not only determined by the bubble size but also by the flow velocity. A high heat flux can lead to a fast growth of the bubble, but the inclination angle forming time can be about the same due to the same Next, the heat flux from the heating wall determines the growth speed of the bubble. A higher heat flux leads to a faster growth of the bubble, which impacts the growth force that opposes the bubble departure. That is, a high heat flux leads to larger bubble departure diameter. In addition, for flow boiling the pressure on the distorted bubble is reduced in the direction facing the flow. This leads to a delay of the bubble lift-off, which has been confirmed by the model calculations. The increase is even relatively small (see Figure 21b). However, Sugrue's experiment only confirmed this effect at a wall orientation angle of 90 • [19]. At an orientation angle of 45 • it is even the other way round when heat flux increases from 50-100 kWm −2 (see Figure 21b). When the heat flux increases to 500 kWm −2 the departure diameter increases sharply in the model calculation for the two orientation angles. The same effect was found in the modelling results based on the setup of Situ et al. [18] shown in Figure 22a. That is, also for a wider range of heat flux. The calculation is for an axial location of 0.57 m from the inlet. At low heat flux, i.e., q ≤ 60 kWm −2 , the bubble even does not leave the heating wall. When the heat flux is between 60 kWm −2 and 120 kWm −2 , the bubble lift-off diameter decreases slightly from 0.2 mm to 0.17 mm. Afterwards, with an increase of the heat flux the bubble lift-off diameter rapidly increases to 0.54 mm at 180 kWm −2 . A credible explanation for these results comes from the dynamic contact angle and dynamic inclination angle. The dynamic inclination angle is not only determined by the bubble size but also by the flow velocity. A high heat flux can lead to a fast growth of the bubble, but the inclination angle forming time can be about the same due to the same flow conditions. In other words, the bubble lift-off period is more or less the same while the growth speed is much higher at high heat flux which means a larger lift-off diameter. However, from Equation (42) the base diameter expansion is dependent on the bubble diameter expansion rate, . r b , and the contact angle, β. A smaller contact angle leads to a faster expansion of the base diameter. A higher heat flux leads to r b , a larger lift off diameter but also a larger contact angle, β, and a smaller lift-off diameter. A high heat flux results in two contrary effects on the base diameter, which is determining the bubble lift-off and departure. These different effects of heat flux on inclination and contact angles are the reasons that the lift off diameter slightly decreases and further increases with increase of heat flux as shown in Figure 22a. flow conditions. In other words, the bubble lift-off period is more or less the same while the growth speed is much higher at high heat flux which means a larger lift-off diameter. However, from Equation (42) the base diameter expansion is dependent on the bubble diameter expansion rate, ṙb, and the contact angle, β. A smaller contact angle leads to a faster expansion of the base diameter. A higher heat flux leads to a higher bubble expansion rate, ṙb, a larger lift off diameter but also a larger contact angle, β, and a smaller lift-off diameter. A high heat flux results in two contrary effects on the base diameter, which is determining the bubble lift-off and departure. These different effects of heat flux on inclination and contact angles are the reasons that the lift off diameter slightly decreases and further increases with increase of heat flux as shown in Figure 22a. Sugrue [19] also investigated the effect of the heating wall orientation angle on the bubble departure diameter, which is accurately reproduced by the present model (see Figure 21c). The maximal bubble departure diameter occurs in the downward facing horizontal configuration, while the smallest occurs in the vertical case. For example, for G= 250 kg m -2 s -1 and 10 K subcooling, the departure diameter is maximal (0.668 mm) when the heating wall is downward facing (wall orientation angle = 0°). It decreases to 0.33 mm at 90° wall orientation angle and then it increases again to 0.6 mm at 180° orientation angle. The force difference between the downward and the upward facing horizontal heating wall is just the direction of buoyancy force relative to the heating wall. In the downward facing case, the total force in the positive perpendicular direction is smaller than that in the upward facing case. It means that the expected inclination angle, θ b,s , (Figure 7 and Equation (43)) is larger in the downward facing case than in the upward facing case when all other conditions are same. Therefore, sliding comes later and bubble size is larger in the downward facing case. In the vertical case, buoyancy enforces bubble sliding which reduces departure diameter. Sugrue [19] also investigated the effect of the heating wall orientation angle on the bubble departure diameter, which is accurately reproduced by the present model (see Figure 21c). The maximal bubble departure diameter occurs in the downward facing horizontal configuration, while the smallest occurs in the vertical case. For example, for G = 250 kg m −2 s −1 and 10 K subcooling, the departure diameter is maximal (0.668 mm) when the heating wall is downward facing (wall orientation angle = 0 • ). It decreases to 0.33 mm at 90 • wall orientation angle and then it increases again to 0.6 mm at 180 • orientation angle. The force difference between the downward and the upward facing horizontal heating wall is just the direction of buoyancy force relative to the heating wall. In the downward facing case, the total force in the positive perpendicular direction is smaller than that in the upward facing case. It means that the expected inclination angle, θ b,s , (Figure 7 and Equation (43)) is larger in the downward facing case than in the upward facing case when all other conditions are same. Therefore, sliding comes later and bubble size is larger in the downward facing case. In the vertical case, buoyancy enforces bubble sliding which reduces departure diameter. However, in the model calculation, the departure diameter decreases from 0.56 mm to 0.351 mm and further increases to 0.57 mm. The difference between model and experiment may be caused by the impact of the bubble geometry. In the experiment, the bubble is an ellipsoid that is pushed by buoyancy against the heating wall while in the model it is still hemispherical or spherical.
The subcooling temperature affects the departure diameter due to condensation. Higher liquid subcooling can reduce bubble departure and lift-off diameter. In the Sugrue's case [19] this impact is small (see Figure 21c). At Sugrue's conditions the bubble remains in the thermal layer after departure and so no condensation occurs. The thickness of the thermal layer is~6 mm for Sugrue's case [19] while it is~1 mm for the case of Situ et al. [18]. This effect becomes obvious in the simulation results for the case of Situ et al. [18] where the effective parameters span a wider range. The lift-off diameter decreases for larger subcooling (see Figure 22c). The lift-off diameter reduces from 0.52 mm at 4 K subcooling to 0 at 20 K subcooling. The latter means that the bubble shrinks and collapses due to condensation even before it leaves the heating wall.
In addition, pressure significantly affects the bubble departure and lift-off diameter. With the pressure predominantly the vapour to liquid density ratio increases. That is, at the same heat flux the volume of bubble is smaller at higher pressure resulting in a smaller departure diameter as confirmed by Sugrue [19]   However, in the model calculation, the departure diameter decreases from 0.56 mm to 0.351 mm and further increases to 0.57 mm. The difference between model and experiment may be caused by the impact of the bubble geometry. In the experiment, the bubble is an ellipsoid that is pushed by buoyancy against the heating wall while in the model it is still hemispherical or spherical. The subcooling temperature affects the departure diameter due to condensation. Higher liquid subcooling can reduce bubble departure and lift-off diameter. In the Sugrue's case [19] this impact is small (see Figure 21c). At Sugrue's conditions the bubble remains in the thermal layer after departure and so no condensation occurs. The thickness of the thermal layer is ~6 mm for Sugrue's case [19] while it is ~1 mm for the case of Situ et al. [18]. This effect becomes obvious in the simulation results for the case of Situ et al. [18] where the effective parameters span a wider range. The lift-off diameter decreases for larger subcooling (see Figure 22c). The lift-off diameter reduces from 0.52 mm at 4 K subcooling to 0 at 20 K subcooling. The latter means that the bubble shrinks and collapses due to condensation even before it leaves the heating wall.
In addition, pressure significantly affects the bubble departure and lift-off diameter. With the pressure predominantly the vapour to liquid density ratio increases. That is, at the same heat flux the volume of bubble is smaller at higher pressure resulting in a smaller departure diameter as confirmed by Sugrue [19]   For the whole database of four different cases, the average relative deviation for prediction of the bubble lift off diameter is 23%, which proves the good predictability of the present model (see Table 2). Most important issue is that no recalibration is necessary for bubble growth except the nucleation site density that is due to the lack of superheat data. For the whole database of four different cases, the average relative deviation for prediction of the bubble lift off diameter is 23%, which proves the good predictability of the present model (see Table 2). Most important issue is that no recalibration is necessary for bubble growth except the nucleation site density that is due to the lack of superheat data.

Conclusions and Outlook
A mechanistic bubble dynamics model for subcooled flow boiling has been developed. The model integrates several well-developed conventional theories and models with or without modification and some new conceptual ideas. It covers the whole bubble life cycle including inertia-controlled and thermal diffusion-controlled growths and departure as well. It includes the microlayer and base diameter expansion contributions to the bubble growth and also condensation at the bubble cap. The force balance equations based on the works of Klausner et al. [17] and Chen et al. [21] are applied to determine the departure and lift-off of the bubble. As a novelty, the bottleneck in the late bubble formation period is modelled. Further, the modelling of dynamic contact angle and inclination angle makes the model applicable under different thermal-hydraulic and geometric conditions without recalibration. The model has been successfully validated against different experimental data from the literature. As such, the model can be well implemented into an Euler-Euler multiphase CFD frame to simulate large-scale heat transfer problems.
A sensitivity analysis proved the model can accurately predict the dependency of bubble departure and lift-off diameter on different parameters such as mass flux, heat flux, subcooling, pressure and heating wall orientation. The results also confirm previous findings such as a reduction of departure and lift-off diameter with increasing inflow velocity and increasing system pressure. The model can also interpret the contradictory findings, such as an increase of departure and lift-off diameter with higher heat flux as found in the experiments of Situ et al. [18] and Sugrue [19] while an opposite behaviour was found in the case of Prodanovic et al. [38]. Quantitatively, a total average deviation of 23% for bubble lift off diameter prediction against data from four different experiments was found which proves a good predictability of the new model. In particular, the comparison covers different mass fluxes, heat fluxes, subcooling, pressure, heating wall orientation angle, pipe designs and fluids.
Beside deviations with respect to sliding length and time, which will be a subject of future analyses, the model still has some shortcomings with respect to an implementation in the Eulerian-Eulerian multiphase CFD approach. The first is the base diameter expansion hypothesis. Although it is a commonly accepted approach, it lacks physical support. In the future, this approach will be critically assessed and further developed to improve the accuracy of the model. The second shortcoming is the different time scale of this submodel and that of the global CFD model. The third is bubble merging during sliding, which is not considered at present. A solution to that may lead to a better prediction of the Prodanovic et al. case [38].
Author Contributions: All the authors have contributed their efforts to complete the paper. H.S. and W.D. were involved in doing the simulation, data curation, model validation, results analysis and writing the paper. D.L. and U.H. supervised the work, reviewed and edited the paper. Prandtl number (-) q heat flux (W/m 2 ) q in heat flux into the heating wall (W/m 2 ) q out total heat flux from heating wall to fluid (W/m 2 ) q e,mi heat flux due to evaporation of microlayer (W/m 2 ) q e,ma heat flux due to evaporation of macrolayer (W/m 2 ) q total total heat flux of a heating wall segment (W/m 2 ) q n,w conduction heat flux between two neighboring segments of heating wall (W/m 2 ) . Q heat transfer rate (W) r coordinate/position r b bubble radius (m) r c radius of curvature (m) r d bubble dryout radius (m) r m,g maximum radius in the initial growth period (m) r w bubble contact radius (base radius) (m) R pipe diameter (m) Re b bubble Reynolds number (-) T b bulk temperature (K) T l liquid temperature (K) T w heating wall temperature (K) T ∞ temperature in the bubble in the inertia-controlled growth period (K) T sat saturation temperature (K) T sub subcooling temperature (K)