An Overall Bubble Diameter Model for the Flow Boiling and Numerical Analysis through Global Information Searching

: Bubble diameter is important for the three-dimensional nucleate boiling two-phase ﬂow simulation. However, the bubble diameter is rarely considered as a variable because it is difﬁcult to estimate. In this paper, a novel bubble diameter model is proposed for the boiling ﬂow. The heat transfer growth, the breakup or coalescence of collision and the liquid impacting separation are considered as the factors affecting the bubble diameter, and three bubble diameter sub-models are developed to calculate the overall bubble diameter based on a departure diameter. In the model, the heat transfer growth is calculated after estimating the bubble location. The collision variation is deduced from the interfacial area concentration equation. The liquid impacting variation is calculated through phase interaction force and gas viscous force. The proposed model is conducted through global information searching in the ﬂow boiling simulation and validated using the void fraction and temperature data from the literature. The effects of the heat transfer growth, the collision, and the separating volume on the bubble diameter are discussed.


Introduction
Nucleate boiling is a highly efficient heat transfer process due to the phase change latent heat transfer and the two-phase forced convection heat transfer in the near-wall region. It has been applied in many industrial facilities, such as nuclear reactor, chemical reactor, pipeline system for oil and gas transport and automotive industry [1].
An appropriate operating condition is the key to ensure the safety and the efficiency of the nucleate boiling technique. The gas-liquid flow with nucleate boiling is usually simulated through the Eulerian-Eulerian approach, in which the two phases are considered continuous, even though it should be a multiphase flow. However, the Eulerian-Eulerian approach is widely accepted and many multiphase interaction models are proposed based on it, such as the drag model, the wall boiling model and the inner-phase heat transfer model. High mass flux subcooled flow boiling was investigated by Kenarsari [2] through two-phase Eulerian numerical simulation and experimental study in horizontal micro-tubes. The effects of heat flux, mass flux, vapor quality, void fraction, and inlet subcooling on heat transfer characteristics are obtained and in a good agreement with the experimental results. Gao [3] introduced a wall heat flux partition algorithm in the water subcooled flow boiling simulation to decide the onset of the significant void, which is important for predicting the void fraction profile. Comparing with Bartolomei and Chanturiya's experiment data [4] and Saha and Zuber's [5] correlation, Gao's work obtains a comparatively accurate liquid subcooling at the onset of the significant void.
Tu [6] pointed out that the bubble size correlation is essential for calculating the total rate of the interfacial momentum, energy and mass transport between the two phases in a low-pressure subcooled flow boiling. Due to the lack of experimental bubble size, the constant bubble diameter or the linear bubble diameter correlations [7,8] are usually used in the flow boiling simulations, which affects the simulation accuracy greatly. Anglart and Nylund [7] modeled the bubble diameter as a linear function of the local supercooling and the maximum bubble diameter in the near-wall region, which is validated only for relatively high-pressure flows. In Končar's [8] research, a simple model of radial distribution of the bubble diameter is proposed, which introduces the exponent function of the local subcooling to describe the bubble diameter decrease due to the condensation in the subcooled flow field. Similarly, the model has a good agreement with the low pressure experimental data. Byong [9] also indicated that the bubble size is an importance factor for accurately predicting the subcooled flow boiling, adopted a volumetric conserved bubble size model and introduced new wall boiling and two-phase logarithmic wall function models. The volumetric conserved bubble size model proposed by the population balance equation [10] considers the breakup and the coalescence based on pre-assumed lognormal bubble size distribution. The benchmark calculation shows that the model can predict the high pressure experimental data fairly well. Colombo [11] also adopted the above model to predict the bubble diameter distribution, and employed the heat flux partitioning approach to establish a wall boiling model, which can provide relatively accurate Eulerian-Eulerian CFD (Computational Fluid Dynamics) predicted results except for the average bubble diameter and the mean velocity. Although the bubble size models can provide acceptable simulations, they are only verified partially in the subcooled flow boiling either in high pressure case or in low pressure case.
Focusing on the variation mechanism of the bubble diameter and its effect on the flow boiling simulation, a novel mathematical model of the bubble diameter is developed in this paper. The overall bubble diameter consists of the departure diameter (the initial diameter) and three terms contributed from the heat transfer growth, the collision and the fluid impact. On this basis, a two-phase Eulerian flow boiling simulation framework is proposed. Instead of using constant or temperature-dependent linear bubble diameter, multi-phase flow subcooled flow boiling simulation is carried out with consideration of the nonlinearity of bubble variation. The simulation approach is verified by experimental results from the literature. Then, a case study is carried out under different operating conditions.

Bubble Behavior in Flow Boiling
As shown in Figure 1, the flow domain is usually divided into two regions for the wall flow boiling: (Region 1) the near-wall region where the gas bubble generates, grows, slides and lifts up due to superheating temperature, and (Region 2) the main-flow region where the gas bubble flows following the liquid and condensing to transfer heat to the supercooled liquid.
Energies 2018, 11, x FOR PEER REVIEW 2 of 19 Saha and Zuber's [5] correlation, Gao's work obtains a comparatively accurate liquid subcooling at the onset of the significant void. Tu [6] pointed out that the bubble size correlation is essential for calculating the total rate of the interfacial momentum, energy and mass transport between the two phases in a low-pressure subcooled flow boiling. Due to the lack of experimental bubble size, the constant bubble diameter or the linear bubble diameter correlations [7,8] are usually used in the flow boiling simulations, which affects the simulation accuracy greatly. Anglart and Nylund [7] modeled the bubble diameter as a linear function of the local supercooling and the maximum bubble diameter in the near-wall region, which is validated only for relatively high-pressure flows. In Končar's [8] research, a simple model of radial distribution of the bubble diameter is proposed, which introduces the exponent function of the local subcooling to describe the bubble diameter decrease due to the condensation in the subcooled flow field. Similarly, the model has a good agreement with the low pressure experimental data. Byong [9] also indicated that the bubble size is an importance factor for accurately predicting the subcooled flow boiling, adopted a volumetric conserved bubble size model and introduced new wall boiling and two-phase logarithmic wall function models. The volumetric conserved bubble size model proposed by the population balance equation [10] considers the breakup and the coalescence based on pre-assumed lognormal bubble size distribution. The benchmark calculation shows that the model can predict the high pressure experimental data fairly well. Colombo [11] also adopted the above model to predict the bubble diameter distribution, and employed the heat flux partitioning approach to establish a wall boiling model, which can provide relatively accurate Eulerian-Eulerian CFD (Computational Fluid Dynamics) predicted results except for the average bubble diameter and the mean velocity. Although the bubble size models can provide acceptable simulations, they are only verified partially in the subcooled flow boiling either in high pressure case or in low pressure case.
Focusing on the variation mechanism of the bubble diameter and its effect on the flow boiling simulation, a novel mathematical model of the bubble diameter is developed in this paper. The overall bubble diameter consists of the departure diameter (the initial diameter) and three terms contributed from the heat transfer growth, the collision and the fluid impact. On this basis, a two-phase Eulerian flow boiling simulation framework is proposed. Instead of using constant or temperature-dependent linear bubble diameter, multi-phase flow subcooled flow boiling simulation is carried out with consideration of the nonlinearity of bubble variation. The simulation approach is verified by experimental results from the literature. Then, a case study is carried out under different operating conditions.

Bubble Behavior in Flow Boiling
As shown in Figure 1, the flow domain is usually divided into two regions for the wall flow boiling: (Region 1) the near-wall region where the gas bubble generates, grows, slides and lifts up due to superheating temperature, and (Region 2) the main-flow region where the gas bubble flows following the liquid and condensing to transfer heat to the supercooled liquid. In the flow boiling simulation, a constant bubble diameter or a linear correlation is easy to use, but it would be difficult to adequately describe the bubble size variation flow. After the bubble generation at the wall, the bubble size can be divided into three parts according to the bubble  In the flow boiling simulation, a constant bubble diameter or a linear correlation is easy to use, but it would be difficult to adequately describe the bubble size variation flow. After the bubble generation at the wall, the bubble size can be divided into three parts according to the bubble transformation manner. Under this consideration, the bubble departure diameter can be regarded as the initial bubble size. Phase change is the major manner changing the bubble size. Evaporation occurs in the near wall region, and condensation occurs in the subcooled liquid with mass transfer. The other two forms changing the bubble size can be seen in Figure 2. One is the collision that leads to coalescence or breakup, and the other is the crush because the liquid impact makes the phase interaction force exceed the maximum gas viscosity force. In multiphase simulation, the bubbles are usually equivalent to spheres and the equivalent diameter is to describe the bubble size.
Energies 2018, 11, x FOR PEER REVIEW  3 of 19 transformation manner. Under this consideration, the bubble departure diameter can be regarded as the initial bubble size. Phase change is the major manner changing the bubble size. Evaporation occurs in the near wall region, and condensation occurs in the subcooled liquid with mass transfer. The other two forms changing the bubble size can be seen in Figure 2. One is the collision that leads to coalescence or breakup, and the other is the crush because the liquid impact makes the phase interaction force exceed the maximum gas viscosity force. In multiphase simulation, the bubbles are usually equivalent to spheres and the equivalent diameter is to describe the bubble size.

Bubble Diameter Model
The bubble diameter, db, is defined as follows: where d d is the bubble departure diameter taken as the initial size of the bubble, heat d is the bubble size variation contribution due to the heat transfer, collision d is the bubble size variation contribution due to the collision, and V is the separating volume function representing the bubble size variation due to the phase interaction force.

Bubble Departure Diameter
For the critical situation of bubble departure, many correlations [12][13][14][15] have been presented by deriving the force equilibrium equation including the surface tension force, the drag force, the buoyancy force, the lift force, etc. Some modifications [16][17][18] have also been done considering the heat transfer effect through the Jacob number and the contact angle. According to the lumped energy balance during the bubble growth and lift up process, correlations [19][20][21] of the bubble departure diameter are proposed as functions of the pressure, the saturation temperature and latent heat in a superheated liquid. Although some of them achieve good predictions in the bubble departure diameter, they are more proper for the pool boiling case rather than the flow boiling case. Comparing with the models mentioned above, the Tolubinsky model [22] is more appropriate for the flow boiling [1,3,23], which is given as follows 45 min 0.0014, 0.0006 where w T  is the temperature difference between the wall and the liquid.

Bubble Diameter Model
The bubble diameter, d b , is defined as follows: where d d is the bubble departure diameter taken as the initial size of the bubble, d heat is the bubble size variation contribution due to the heat transfer, d collision is the bubble size variation contribution due to the collision, and V is the separating volume function representing the bubble size variation due to the phase interaction force.

Bubble Departure Diameter
For the critical situation of bubble departure, many correlations [12][13][14][15] have been presented by deriving the force equilibrium equation including the surface tension force, the drag force, the buoyancy force, the lift force, etc. Some modifications [16][17][18] have also been done considering the heat transfer effect through the Jacob number and the contact angle. According to the lumped energy balance during the bubble growth and lift up process, correlations [19][20][21] of the bubble departure diameter are proposed as functions of the pressure, the saturation temperature and latent heat in a superheated liquid. Although some of them achieve good predictions in the bubble departure diameter, they are more proper for the pool boiling case rather than the flow boiling case. Comparing with the models mentioned above, the Tolubinsky model [22] is more appropriate for the flow boiling [1,3,23], which is given as follows where ∆T w is the temperature difference between the wall and the liquid.

Bubble Size Due to Heat Transfer
In the subcooled flow boiling, the temperature of gas or liquid varies with time and depends on the bubble flow route. The changing temperature leads to a varying heat transfer growth rate of the bubble. To calculate the bubble heat transfer growth rate, two key variables should be given: (1) the temperature variation; and (2) the capacity of bubble growth rate. The relation between the heat transfer growth rate and the temperature difference is established based on Zuber's model [19]. The bubble is assumed to be a sphere considering the two-way heat transfer from the liquid thermal boundary layer to the bubble and the cooler liquid in the stabilized supersaturated temperature environment, where t is the time of heat transfer, k l is the heat conductivity, q b is the heat transfer from the liquid thermal boundary layer to the liquid. Especially, q b equals zero in the constant liquid temperature field. Based on the constant fluid property assumption, the bubble size increases if the temperature difference is positive and vice versa. In the model, the environment temperature is uniform and the heat transfers are the same in different directions. To take into account the nonuniformity of the environment temperature, an approximate average heat loss of 50% from the total is used. The bubble diameter growth rate due to heat transfer can be calculated as The Eulerian-Eulerian multiphase simulation cannot identify the bubble route, so the temperature difference and the square root of time are not available. An approximate method for the temperature difference and the time is applied in this paper to solve the problem. A variable n is introduced as the ratio between the flow distance S and the unit distance dS.
Then, the average temperature difference can be obtained where T wl and T wg are the liquid and gas temperatures at the wall, and T l and T g are the liquid and gas temperatures in the main flow. The temperature difference of the near-wall region is the sum of the initial temperature difference and the average temperature variation. In the main-flow region, the bubble enters the single phase flow as an external object, so the initial temperature difference is zero and only the average temperature variation is considered. The liquid or gas temperature will be replaced by the saturation temperature T sat considering the unsaturation heat transfer. The time term can be calculated by where S is the flow distance, and u g is the velocity. The flow distance S is a function of location consisted of two parts, the departure and departure-lifting distance S d and the flow distance S f , as shown in Figure 1. Because of the complexity of the flow boiling, the criterion of the bubble merging with the main flow is difficult to determine exactly. In this work, the liquid velocity is used as a simplified criterion to distinguish whether the bubble has left the near-wall region. In the near-wall region, the bubble path is the syntheses of the radial and axial distances. In the main-flow region, the bubble path is simplified as the mean axial distance without disturbances. Thus, the flow distance S can be calculated as u l u l,radial 2 , u l < C × u l,max for the near − wall region (region 1) (9) where R is the pivot distance, R tube is the tube radius, S 0 is the axial distance from inlet, u l,max is the liquid phase main stream velocity, and C is the bubble location factor.

Bubble Size Due to Collision
The contributions of bubble collision coalescence and breakup are extracted from the interfacial area concentration and Hibiki [24] research. The bubble coalescence rate and the bubble breakup rate can be described in an interfacial area concentration transport equation [24] as follows, where φ bc is the changing rate of the interfacial area concentration due to bubble breakup or coalescence, φ phase is that due to phase change, and φ pressure is that due to pressure change.
In the case of only considering the breakup or coalescence due to collision, the last two terms on the right hand side of Equation (11) can be neglected. Then, it becomes where φ RC is the bubble coalescence rate and φ TI is the bubble breakup rate, which can be given as [24]: where f C is the bubble random collision frequency, f B is the bubble eddy random collision frequency, λ C is the coalescence efficiency due to bubble random collision induced by turbulence in liquid phase, and λ B is the breakup efficiency due to collision of the turbulent eddy with the bubble. The collision and break up are ignored due to low rigidity.
where Kc is 1.29 and K B is 1.37 [24]. The bubble diameter is considered as the variable in the transport equation, and the influence of bubble collision is represented by the interface concentration variation. It accurately represents the change of the interfacial area concentration but ignores the change of the bubble diameter itself. Interfacial area concentration is the interfacial area between two phases of per unit mixture volume. To calculated bubble collision, the interfacial area concentration is calculated in algebra by the Particle model which is derived from the surface area to volume ratio [3,25], Through Equation (19), the collision effect is considered in the bubble diameter calculation directly. When a i is replaced by d b and only the time gradient is considered, the contributions of coalescence and breakup mechanisms by collision to the averaged bubble size variations can be obtained: where n/2 is introduced to consider the average distance from 0 to n times of unit distance, dt = dS/u g is the unit time.

Bubble Size Due to Phase Interaction Force
If the external force difference ∆F exceeds the viscosity force and the deformation reaches the critical size, the bubble will break up because the intermolecular forces cannot maintain the bubble shape anymore, and the external force difference has a similar tearing effect on the bubble. In the flow boiling, the external force is mainly caused by the liquid impact and can be calculated by the drag force difference. Assuming the separating gas velocity is equal to the liquid velocity, the separating volume function can be derived by the momentum theorem (the force on a particle is equal to the rate of change of momentum of that particle) as follows where ∆u is the relative velocity between liquid and gas, F vis is the bubble viscous force, grad(F) is the gradient of phase interaction force, and ∆F is calculated by the phase interaction force difference of two adjacent units. Equation (20) considers the eddies influence to supplement Equation (21) which cannot consider the liquid cutting effect in different directions.

Overall Bubble Diameter and Information Searching
The overall bubble diameter can be described as follows, where d b is deduced by the bubble volume. The bubble diameter model is compiled in Fluent software through C language. The calculation frame consists of four parts, as shown in Figure 3: the blue route is for the bubble departure diameter, the red route is for the heat transfer bubble size contribution, the yellow route is for the collision bubble size contribution and the green route is for the separating volume. Information of general multi-phase flow is not sufficient for the proposed new bubble diameter model. As a result, adding some data from the whole flow field into the computation is essential in the calculation, which tightly connects the computation variables and the flow condition together. The information searching scheme can be thought as a reform of the computation element regarding the data collection. Firstly, the main flow velocity is needed to determine the element location in the flow, for which only nearby elements are enough. Given as Loop 1, the calculated range is obtained according to the distance in the flow direction under the first global search, and the element velocities in the range are compared to fix the location. The second information searched is the temperature difference at the wall for determining the initial temperature. Seen as Loop 2, boundary elements nearby the presented element determined by the second global search are used to average the wall temperature. The phase interaction force difference ∆F between elements is the third information needed to represent the tear force of liquid impact. In Loop 3, a series of subtractions are done using the phase interaction forces of neighbor elements determined by the third global search to obtain the maximum force difference. Another global searching task over logic level of the three loops searching is conducted before every iteration of the computation to allocate information to the computation elements. In other words, information searching and allocating are the cores of the present computation apart from solving the multi-phase equations and the turbulence equations.
The bubble diameter model is compiled in Fluent software through C language. The calculation frame consists of four parts, as shown in Figure 3: the blue route is for the bubble departure diameter, the red route is for the heat transfer bubble size contribution, the yellow route is for the collision bubble size contribution and the green route is for the separating volume. Information of general multi-phase flow is not sufficient for the proposed new bubble diameter model. As a result, adding some data from the whole flow field into the computation is essential in the calculation, which tightly connects the computation variables and the flow condition together. The information searching scheme can be thought as a reform of the computation element regarding the data collection. Firstly, the main flow velocity is needed to determine the element location in the flow, for which only nearby elements are enough. Given as Loop 1, the calculated range is obtained according to the distance in the flow direction under the first global search, and the element velocities in the range are compared to fix the location. The second information searched is the temperature difference at the wall for determining the initial temperature. Seen as Loop 2, boundary elements nearby the presented element determined by the second global search are used to average the wall temperature. The phase interaction force difference F  between elements is the third information needed to represent the tear force of liquid impact. In Loop 3, a series of subtractions are done using the phase interaction forces of neighbor elements determined by the third global search to obtain the maximum force difference. Another global searching task over logic level of the three loops searching is conducted before every iteration of the computation to allocate information to the computation elements. In other words, information searching and allocating are the cores of the present computation apart from solving the multi-phase equations and the turbulence equations.

Gas-Liquid Two Phase Flow Conservation Equation
The continuity equation for phase q: The momentum equation for phase q: The energy equation for phase q:

Inner-Phase Mass Transfer and Heat Transfer
In the interior of the flow, the mass transfer rate between the two phases depends upon the temperature of the liquid and the gas, the interphase heat transfer coefficient and the latent heat. The temperature affects the mass transfer direction. The mass transfer rate where a i is the interfacial area concentration and can be calculated by Equation (19). The interfacial heat transfer coefficient h inter is modeled by Ranz-Marshall [26] correlation, where Re p is the relative Reynolds number regarding the diameter of the p phase and the relative velocity and Pr is the Prandtl number of the q phase.

Wall Boiling Model
The manners of the wall boiling heat transfer include the liquid convection heat transfer, the gas convection heat transfer, the latent heat transfer by liquid evaporation, and the quench heat transfer, as shown in Figure 4. The gas convection is not obvious when the gas bubble is little. The quench heat transfer is resulted by the liquid and wall contact and the gas and wall contact successively. The wall boiling heat equation is where Q w is the total heat flux, Q con is the convective heat flux, and Q e is the gas evaporation latent heat flux. In fact, the bubble nucleation starts at the wall due to the local overheating caused by the wall roughness or other factors. The bubble keeps growing for phase exchange and heat expansion until it departs from the wall for force equilibrium. Usually, the bubble departure diameter, the bubble departure frequency and the nucleation site density are used to describe the wall heat transfer. Thus, the heat flux can be given as where h con is the convection heat transfer coefficient between the liquid and the wall, A b is the wall surface covered by gas, (1 − A b ) is the wall surface covered by liquid, m b is the bubble mass regarding the bubble departure diameter, N w is the nucleation site density, and f is the bubble departure frequency. The wall surface covered by gas, A b , is a function of the bubble departure diameter, d d , and the nucleation sites density, N w , as follows where K is the empirical constant calculated by Del Valle and Kenning correlation [27] and Ja is the Jakob number.
Energies 2018, 11, x FOR PEER REVIEW 9 of 19 Usually, the bubble departure diameter, the bubble departure frequency and the nucleation site density are used to describe the wall heat transfer. Thus, the heat flux can be given as where con h is the convection heat transfer coefficient between the liquid and the wall, b A is the wall surface covered by gas,   where K is the empirical constant calculated by Del Valle and Kenning correlation [27] and Ja is the

Interphase Momentum Transfer
The velocity of each phase changes with the phase interaction force through the momentum exchange. The phase interaction force where → F li f t is the lift force acted on a particle, droplet or bubble mainly due to velocity gradients in the primary-phase flow field, → F vm is the virtual mass force that occurs when a secondary phase accelerates relative to the primary phase, → F TD is the turbulent dispersion force that acts as a turbulent diffusion in the dispersed flow accounting for the interphase turbulent momentum transfer, → F W L is the wall lubrication force that acts in the lateral direction and tends to push the secondary phases away from the wall, and → F drag is the drag force aroused by the interfacial friction force and pressure differences that acts on the vapor bubble. Among these interphase forces, the drag force, → F drag , is the dominant component. The drag force per unit volume can be expressed as [25,30] where ∆u is the relative velocity between the liquid and the gas, and C D is the drag force coefficient that can be given as [31] C D = 24 1 + 0.15Re 0.687 /Re Re ≤ 1000 The lift force can be calculated as follows [25,30] → where C lift is lift force coefficient that can be calculated [32].
The virtual mass force can be obtained [25,30] where C vm is the virtual mass force coefficient and given as 0.5. The wall lubrication force is [25,30] where C W L is the wall lubrication force coefficient [33].
where y w is the distance to the nearest wall.

Model Verification
The first step in the model verification is to determine the bubble location. The experimental data from Reference [35] is used to calculate the value of the coefficient, C, which reflects the bubble location information. The subcooled fluid without the second phase gas flows into the heating pipeline after a free flow, and the experiment is observing the heating pipeline fluid flow. The simulation ignores the free flow, and the uniform subcooled liquid is set as the initial flow. The sketch of the experiment apparatus is shown in Figure 5. The flow rate G is 283.1 kg/m 2 s, the subcooled temperature is less than the saturation temperature of Pressure 1.23 bar 19.7 K.

Model Verification
The first step in the model verification is to determine the bubble location. The experimental data from Reference [35] is used to calculate the value of the coefficient, C , which reflects the bubble location information. The subcooled fluid without the second phase gas flows into the heating pipeline after a free flow, and the experiment is observing the heating pipeline fluid flow. The simulation ignores the free flow, and the uniform subcooled liquid is set as the initial flow. The sketch of the experiment apparatus is shown in Figure 5. The flow rate G is 283.1 kg/m 2 s, the subcooled temperature is less than the saturation temperature of Pressure 1.23 bar 19.7 K. Four levels of C are chosen to compare the corresponding void fraction with the experimental results, which are 0.6, 0.7, 0.8 and 0.9, as shown in Figure 6. Along the flow direction, the experimental void fraction increases from the inlet to the outlet. For the C = 0.6 and 0.8 cases, there are steep decreases of the void fraction around 75 mm away from the inlet, which are opposite to the experimental trend. For the C = 0.7 and 0.9 cases, the distributions of the void fraction along the tube length are similar to the experiment. Comparing with the C = 0.7 case, the C = 0.9 case provides a closer prediction to the experimental data, in which the relative error is smaller than 15%. The biggest gap between the C = 0.9 simulation and the experiment is in the near-inlet section. Since the simulation ignores the free flow section in front of the heating section in the experiment, the error can be accepted. As a result, C = 0.9 is adopted in the bubble diameter model. The experiment conducted in Reference [4] is used for the model verification, as shown in Figure 7. The void fraction and the temperature along the tube length are calculated under the same operating condition using the proposed model and compared with the experiment. Figure 8 shows the calculated void fraction and temperature overlaid by the experimental results. It can be seen in Figure 8a that the calculated void fraction agrees with the experiment data well with relative errors no more than 20% in most cases. The biggest difference between the simulation and the experiment is at the position 1.5 m away from the inlet, where the liquid temperature has a higher value than the Four levels of C are chosen to compare the corresponding void fraction with the experimental results, which are 0.6, 0.7, 0.8 and 0.9, as shown in Figure 6. Along the flow direction, the experimental void fraction increases from the inlet to the outlet. For the C = 0.6 and 0.8 cases, there are steep decreases of the void fraction around 75 mm away from the inlet, which are opposite to the experimental trend. For the C = 0.7 and 0.9 cases, the distributions of the void fraction along the tube length are similar to the experiment. Comparing with the C = 0.7 case, the C = 0.9 case provides a closer prediction to the experimental data, in which the relative error is smaller than 15%. The biggest gap between the C = 0.9 simulation and the experiment is in the near-inlet section. Since the simulation ignores the free flow section in front of the heating section in the experiment, the error can be accepted. As a result, C = 0.9 is adopted in the bubble diameter model.

Model Verification
The first step in the model verification is to determine the bubble location. The experimental data from Reference [35] is used to calculate the value of the coefficient, C , which reflects the bubble location information. The subcooled fluid without the second phase gas flows into the heating pipeline after a free flow, and the experiment is observing the heating pipeline fluid flow. The simulation ignores the free flow, and the uniform subcooled liquid is set as the initial flow. The sketch of the experiment apparatus is shown in Figure 5. The flow rate G is 283.1 kg/m 2 s, the subcooled temperature is less than the saturation temperature of Pressure 1.23 bar 19.7 K. Four levels of C are chosen to compare the corresponding void fraction with the experimental results, which are 0.6, 0.7, 0.8 and 0.9, as shown in Figure 6. Along the flow direction, the experimental void fraction increases from the inlet to the outlet. For the C = 0.6 and 0.8 cases, there are steep decreases of the void fraction around 75 mm away from the inlet, which are opposite to the experimental trend. For the C = 0.7 and 0.9 cases, the distributions of the void fraction along the tube length are similar to the experiment. Comparing with the C = 0.7 case, the C = 0.9 case provides a closer prediction to the experimental data, in which the relative error is smaller than 15%. The biggest gap between the C = 0.9 simulation and the experiment is in the near-inlet section. Since the simulation ignores the free flow section in front of the heating section in the experiment, the error can be accepted. As a result, C = 0.9 is adopted in the bubble diameter model. The experiment conducted in Reference [4] is used for the model verification, as shown in Figure 7. The void fraction and the temperature along the tube length are calculated under the same operating condition using the proposed model and compared with the experiment. Figure 8 shows the calculated void fraction and temperature overlaid by the experimental results. It can be seen in Figure 8a that the calculated void fraction agrees with the experiment data well with relative errors no more than 20% in most cases. The biggest difference between the simulation and the experiment is at the position 1.5 m away from the inlet, where the liquid temperature has a higher value than the The experiment conducted in Reference [4] is used for the model verification, as shown in Figure 7. The void fraction and the temperature along the tube length are calculated under the same operating condition using the proposed model and compared with the experiment. Figure 8 shows the calculated void fraction and temperature overlaid by the experimental results. It can be seen in Figure 8a that the calculated void fraction agrees with the experiment data well with relative errors no more than 20% in most cases. The biggest difference between the simulation and the experiment is at the position 1.5 m away from the inlet, where the liquid temperature has a higher value than the experiment, as shown in Figure 8b. The gas phase generates from axial distance about 0.6 m, and has a sharp increase at 1.1 m. The void fraction arrives at 40.74% at the outlet. Similar errors occur in the liquid temperature, as shown in Figure 8b. The biggest error is at the void fraction of 0.2 because the wall boiling heat transfer becomes non-equilibrium subcooled boiling when the void fraction is more than 0.2. It means the gas convection heat transfer cannot be ignored in the non-equilibrium subcooled boiling. The relative error of the wall temperature is smaller. The wall temperature reaches a peak of 543 K at 1.54 m, then falls to 541 K and keeps stable until the outlet. Comparing with the experimental data, the precision of the proposed model can be accepted in the flow boiling evaluation. Moreover, the bubble diameter model developed in this paper is appropriate for different operation conditions, in particular for different pressures. experiment, as shown in Figure 8b. The gas phase generates from axial distance about 0.6 m, and has a sharp increase at 1.1 m. The void fraction arrives at 40.74% at the outlet. Similar errors occur in the liquid temperature, as shown in Figure 8b. The biggest error is at the void fraction of 0.2 because the wall boiling heat transfer becomes non-equilibrium subcooled boiling when the void fraction is more than 0.2. It means the gas convection heat transfer cannot be ignored in the non-equilibrium subcooled boiling. The relative error of the wall temperature is smaller. The wall temperature reaches a peak of 543 K at 1.54 m, then falls to 541 K and keeps stable until the outlet. Comparing with the experimental data, the precision of the proposed model can be accepted in the flow boiling evaluation. Moreover, the bubble diameter model developed in this paper is appropriate for different operation conditions, in particular for different pressures. experiment, as shown in Figure 8b. The gas phase generates from axial distance about 0.6 m, and has a sharp increase at 1.1 m. The void fraction arrives at 40.74% at the outlet. Similar errors occur in the liquid temperature, as shown in Figure 8b. The biggest error is at the void fraction of 0.2 because the wall boiling heat transfer becomes non-equilibrium subcooled boiling when the void fraction is more than 0.2. It means the gas convection heat transfer cannot be ignored in the non-equilibrium subcooled boiling. The relative error of the wall temperature is smaller. The wall temperature reaches a peak of 543 K at 1.54 m, then falls to 541 K and keeps stable until the outlet. Comparing with the experimental data, the precision of the proposed model can be accepted in the flow boiling evaluation. Moreover, the bubble diameter model developed in this paper is appropriate for different operation conditions, in particular for different pressures.

Case Study and Discussion
With the validated model, a case study is carried out to discuss the availability of the proposed model in different conditions. The simulations are conducted with a 300 mm long horizontal tube with inner diameter of 7 mm, which is long enough to produce bubbles and the gravity effect can be considered, as shown in Figure 9. The operation conditions are listed in Table 1.  The calculated bubble size distribution along the tube is shown in Figure 10. For the horizontal tube studied, the flow in this specific case is affected by the gravity, which makes the gas phase concentrate in the upper part of the tube when the velocity is relatively low. Overall, the distribution of the bubble size is different from case to case. It is because different boundaries lead to different superheat, and the evaporation/condensation rate varies accordingly. The amount of the gas phase changes the flow type and the phase interaction force decides the velocity. Consequently, the bubble diameter components vary along with the flow development and result in totally different bubble size distributions. The maximum of bubble is approximately 1.4 mm in all three cases, and the differentiation is the distance that produces the large bubbles. In other words, the maximum size of bubbles that can be formed is 1.4 mm under such a pipe diameter and thermal conditions. From the comparison between Cases (a) and (b), it can be seen that a higher wall temperature tends to increase the bubble size except for the inlet area. The comparison between Cases (b) and (c) indicates the bubble size decreases with increasing operating pressure, especially in the middle section of the

Case Study and Discussion
With the validated model, a case study is carried out to discuss the availability of the proposed model in different conditions. The simulations are conducted with a 300 mm long horizontal tube with inner diameter of 7 mm, which is long enough to produce bubbles and the gravity effect can be considered, as shown in Figure 9. The operation conditions are listed in Table 1.

Case Study and Discussion
With the validated model, a case study is carried out to discuss the availability of the proposed model in different conditions. The simulations are conducted with a 300 mm long horizontal tube with inner diameter of 7 mm, which is long enough to produce bubbles and the gravity effect can be considered, as shown in Figure 9. The operation conditions are listed in Table 1.  The calculated bubble size distribution along the tube is shown in Figure 10. For the horizontal tube studied, the flow in this specific case is affected by the gravity, which makes the gas phase concentrate in the upper part of the tube when the velocity is relatively low. Overall, the distribution of the bubble size is different from case to case. It is because different boundaries lead to different superheat, and the evaporation/condensation rate varies accordingly. The amount of the gas phase changes the flow type and the phase interaction force decides the velocity. Consequently, the bubble diameter components vary along with the flow development and result in totally different bubble size distributions. The maximum of bubble is approximately 1.4 mm in all three cases, and the differentiation is the distance that produces the large bubbles. In other words, the maximum size of bubbles that can be formed is 1.4 mm under such a pipe diameter and thermal conditions. From the comparison between Cases (a) and (b), it can be seen that a higher wall temperature tends to increase the bubble size except for the inlet area. The comparison between Cases (b) and (c) indicates the bubble size decreases with increasing operating pressure, especially in the middle section of the tube. The effect of pressure can also be regarded as the effect of the subcooled temperature because a higher operation pressure can bring a higher saturation temperature. Due to the constant physical parameters in these cases with, the Eotvos number related to the diameter of the bubble is a diverse distribution with a maximum of 3. 19 and minimum of about 0, and the Morton  The calculated bubble size distribution along the tube is shown in Figure 10. For the horizontal tube studied, the flow in this specific case is affected by the gravity, which makes the gas phase concentrate in the upper part of the tube when the velocity is relatively low. Overall, the distribution of the bubble size is different from case to case. It is because different boundaries lead to different superheat, and the evaporation/condensation rate varies accordingly. The amount of the gas phase changes the flow type and the phase interaction force decides the velocity. Consequently, the bubble diameter components vary along with the flow development and result in totally different bubble size distributions. The maximum of bubble is approximately 1.4 mm in all three cases, and the differentiation is the distance that produces the large bubbles. In other words, the maximum size of bubbles that can be formed is 1.4 mm under such a pipe diameter and thermal conditions. From the comparison between Cases (a) and (b), it can be seen that a higher wall temperature tends to increase the bubble size except for the inlet area. The comparison between Cases (b) and (c) indicates the bubble size decreases with increasing operating pressure, especially in the middle section of the tube.
The effect of pressure can also be regarded as the effect of the subcooled temperature because a higher operation pressure can bring a higher saturation temperature. Due to the constant physical parameters in these cases with, the Eotvos number Eo = g∆ρd 2 b σ related to the diameter of the bubble is a diverse distribution with a maximum of 3. 19 and minimum of about 0, and the Morton number Mo = µ 4 g∆ρ ρ 2 σ 3 is a constant, equal to 2.885 × 10 −6 . Thus, only some smaller shape deformations of the bubble might happen in these cases. The external forces mainly change the size of the bubbles.  Six sections along the tube length are chosen to show the bubble diameter distribution in the axial view in Case (b), as shown in Figure 11. In Sections 1 and 6, the large bubble size concentrates at the center of the tube. In the four other sections, the big bubbles tend to generate in the upper half of the tube. A pattern of the bubble size distribution can be concluded that bubbles of large diameter move to the center of the tube when the flow is close to the outlet, which is affected by the exit ejection of the tube. If the constant value bubble diameter or approximately linear and temperature-dependent bubble diameter model applied in these cases, there will be a fixed value in entire flow field regardless of temperature and speed change or a series of rings bubble diameter distribution, which are all inconsistent with actual fact. Bubble diameter distributions due to heat transfer are shown in Figure 12, those due to collision are shown in Figure 13, and those due to fluid impact are shown in Figure 14. In Figure 12, the maximum bubble diameter appears in the near-wall region where the liquid temperature is the highest. The flow boiling occurs in the subcooled condition, so the growth rate is less than zero according to the bubble diameter decrease. The upper part of the front half of the tube owns a greater heat transfer growth rate, as shown in Figure 12 (Sections 2 and 3). The high heat transfer growth rate region varies along the tube and the area reduces as the flow approaches to the outlet because the temperature difference gets smaller. Ignoring the wall (the blue boundary), the bubble diameter distribution along the radial direction has local changes due to turbulence. Apart from that, the heat transfer results in an annular distribution of the bubble diameter in the subcooled boiling. The components due to heat transfer change from 0 mm to 1.4 mm, which is close to the maximum bubble diameter. In other words, the heat transfer contributes most to the bubble size in this case. Six sections along the tube length are chosen to show the bubble diameter distribution in the axial view in Case (b), as shown in Figure 11. In Sections 1 and 6, the large bubble size concentrates at the center of the tube. In the four other sections, the big bubbles tend to generate in the upper half of the tube. A pattern of the bubble size distribution can be concluded that bubbles of large diameter move to the center of the tube when the flow is close to the outlet, which is affected by the exit ejection of the tube. If the constant value bubble diameter or approximately linear and temperature-dependent bubble diameter model applied in these cases, there will be a fixed value in entire flow field regardless of temperature and speed change or a series of rings bubble diameter distribution, which are all inconsistent with actual fact.   Figure 11. In Sections 1 and 6, the large bubble size concentrates at the center of the tube. In the four other sections, the big bubbles tend to generate in the upper half of the tube. A pattern of the bubble size distribution can be concluded that bubbles of large diameter move to the center of the tube when the flow is close to the outlet, which is affected by the exit ejection of the tube. If the constant value bubble diameter or approximately linear and temperature-dependent bubble diameter model applied in these cases, there will be a fixed value in entire flow field regardless of temperature and speed change or a series of rings bubble diameter distribution, which are all inconsistent with actual fact. Bubble diameter distributions due to heat transfer are shown in Figure 12, those due to collision are shown in Figure 13, and those due to fluid impact are shown in Figure 14. In Figure 12, the maximum bubble diameter appears in the near-wall region where the liquid temperature is the highest. The flow boiling occurs in the subcooled condition, so the growth rate is less than zero according to the bubble diameter decrease. The upper part of the front half of the tube owns a greater heat transfer growth rate, as shown in Figure 12 (Sections 2 and 3). The high heat transfer growth rate region varies along the tube and the area reduces as the flow approaches to the outlet because the temperature difference gets smaller. Ignoring the wall (the blue boundary), the bubble diameter distribution along the radial direction has local changes due to turbulence. Apart from that, the heat transfer results in an annular distribution of the bubble diameter in the subcooled boiling. The components due to heat transfer change from 0 mm to 1.4 mm, which is close to the maximum bubble diameter. In other words, the heat transfer contributes most to the bubble size in this case. Bubble diameter distributions due to heat transfer are shown in Figure 12, those due to collision are shown in Figure 13, and those due to fluid impact are shown in Figure 14. In Figure 12, the maximum bubble diameter appears in the near-wall region where the liquid temperature is the highest. The flow boiling occurs in the subcooled condition, so the growth rate is less than zero according to the bubble diameter decrease. The upper part of the front half of the tube owns a greater heat transfer growth rate, as shown in Figure 12 (Sections 2 and 3). The high heat transfer growth rate region varies along the tube and the area reduces as the flow approaches to the outlet because the temperature difference gets smaller. Ignoring the wall (the blue boundary), the bubble diameter distribution along the radial direction has local changes due to turbulence. Apart from that, the heat transfer results in an annular distribution of the bubble diameter in the subcooled boiling. The components due to heat transfer change from 0 mm to 1.4 mm, which is close to the maximum bubble diameter. In other words, the heat transfer contributes most to the bubble size in this case.  Figure 13 shows the distribution of the bubble diameter variation due to collision. collision d is positive in the main-flow region while it is negative in the near-wall region. That is to say, when the collision occurs in the main flow region, the possibility of coalescence is higher than that of breakup. In contrast, the collision in the near-wall region would result in more breakups than coalescences. The breakup contribution is in the range of 10 −8 to 10 −7 m while the coalescence contribution is 10 −9 to 10 −12 m. Even if the integral is considered, those two terms are still not significant (10 −5 to 10 −4 m for breakup and 10 −6 to 10 −9 m for coalescence). In the front half tube, the influence of collision on the bubble size distribution is very slight even if the integral is considered. In the rear half, the proportion of the gas phase increases and the break up efficiency during the bubble collisions increases as well in the tube top. It is because the collision becomes more intensive as the gas phase increases to a certain extent, which has a higher efficiency of breaking up. The coalescence efficiency reduces due to the collisions frequency decrease of the full developed flow in the rear half. In short, the collision prompts the bubble diameter distribution to form two regions in the nuclear flow boiling, a ring high incidence area of bubble break up to decrease bubble diameter and a round high incidence area of bubble coalescence to increase bubble diameter.  Figure 14 is the distributions of the separating volume at different sections. A negative V region cannot separate the gas so that the separating volume region in Figure 14 is outside the isoline of V =   Figure 13 shows the distribution of the bubble diameter variation due to collision. collision d is positive in the main-flow region while it is negative in the near-wall region. That is to say, when the collision occurs in the main flow region, the possibility of coalescence is higher than that of breakup.
In contrast, the collision in the near-wall region would result in more breakups than coalescences. The breakup contribution is in the range of 10 −8 to 10 −7 m while the coalescence contribution is 10 −9 to 10 −12 m. Even if the integral is considered, those two terms are still not significant (10 −5 to 10 −4 m for breakup and 10 −6 to 10 −9 m for coalescence). In the front half tube, the influence of collision on the bubble size distribution is very slight even if the integral is considered. In the rear half, the proportion of the gas phase increases and the break up efficiency during the bubble collisions increases as well in the tube top. It is because the collision becomes more intensive as the gas phase increases to a certain extent, which has a higher efficiency of breaking up. The coalescence efficiency reduces due to the collisions frequency decrease of the full developed flow in the rear half. In short, the collision prompts the bubble diameter distribution to form two regions in the nuclear flow boiling, a ring high incidence area of bubble break up to decrease bubble diameter and a round high incidence area of bubble coalescence to increase bubble diameter.  Figure 14 is the distributions of the separating volume at different sections. A negative V region cannot separate the gas so that the separating volume region in Figure 14 is outside the isoline of V = 0 where the phase interaction force exceeds gas viscous force. At the inlet of the tube (Section 1), the separating volume region widely distributes around the wall. In this area, the separating volume  Figure 13 shows the distribution of the bubble diameter variation due to collision. d collision is positive in the main-flow region while it is negative in the near-wall region. That is to say, when the collision occurs in the main flow region, the possibility of coalescence is higher than that of breakup. In contrast, the collision in the near-wall region would result in more breakups than coalescences. The breakup contribution is in the range of 10 −8 to 10 −7 m while the coalescence contribution is 10 −9 to 10 −12 m. Even if the integral is considered, those two terms are still not significant (10 −5 to 10 −4 m for breakup and 10 −6 to 10 −9 m for coalescence). In the front half tube, the influence of collision on the bubble size distribution is very slight even if the integral is considered. In the rear half, the proportion of the gas phase increases and the break up efficiency during the bubble collisions increases as well in the tube top. It is because the collision becomes more intensive as the gas phase increases to a certain extent, which has a higher efficiency of breaking up. The coalescence efficiency reduces due to the collisions frequency decrease of the full developed flow in the rear half. In short, the collision prompts the bubble diameter distribution to form two regions in the nuclear flow boiling, a ring high incidence area of bubble break up to decrease bubble diameter and a round high incidence area of bubble coalescence to increase bubble diameter.
transfer component of 1.4 mm, the collision component of −0.1 mm and the separating volume component of about −10 mm totally break the bubble in the near-wall region. The separating volume makes the overall bubble diameter decrease more obviously compared with the other two diameter components. In the flow, the heat transfer is the main contribution of the increment of the bubble diameter, but the overall diameter cannot increase infinitely because of the departure criterion. It can be described as a process that a growing bubble of a certain size of diameter bears the external force which reduces the bubble diameter and the process consists of growing and shrinking at the same time. The specific values of the increment and the decrement depend on the interaction between factors. The non-equilibrium force produces an annular break-up region surrounding that without break-up to decrease the bubble diameter. The break-up region slightly moves up due to the horizon flow and the gravity.

Conclusions
A new bubble diameter model is developed for the flow boiling simulation. The bubble diameter is divided into four components: the departure diameter, the diameter due to heat transfer, the diameter due to collision, and the diameter due to phase interaction force. Based on the proposed bubble diameter model, a flow boiling simulation scheme is established and carried out through C language and Fluent software. Main conclusions can be drawn: (1) The proposed overall diameter model can be regarded as a combination of the initial diameter and contributions from heat transfer, collision and phase interaction force. It is suitable for more operation conditions and can describe the void fraction and the temperature of the flow boiling accurately in different operating pressures and temperatures.  Figure 14 is the distributions of the separating volume at different sections. A negative V region cannot separate the gas so that the separating volume region in Figure 14 is outside the isoline of V = 0 where the phase interaction force exceeds gas viscous force. At the inlet of the tube (Section 1), the separating volume region widely distributes around the wall. In this area, the separating volume contribution to the overall bubble diameter is around 1 mm, which is big enough to cover the heat transfer term. The orange area with around 1.1 mm in bubble size in Section 1 is the combination of the heat transfer term (1.07 mm) and the collision component (0.01 mm). As the flow develops along the tube length, the gas separating phenomenon becomes increasingly obvious in the first half of the tube. The maximum separating volume occurs at the bottom of the tube and then expands as a wide ring around the wall, as shown in Sections 2 and 3. After that, the gas separating area continues to increase and the maximum region falls back to the top of the near-wall region. The maximum of separating volume component increases from 10 −8 m to 10 −6 m and maintain the size to the outlet, which means the speed difference gradually increases until stability. The separating volume changes also indicate that the flow is fully developed in the rear half of the tube. In Section 5, the heat transfer component of 1.4 mm, the collision component of −0.1 mm and the separating volume component of about −10 mm totally break the bubble in the near-wall region. The separating volume makes the overall bubble diameter decrease more obviously compared with the other two diameter components. In the flow, the heat transfer is the main contribution of the increment of the bubble diameter, but the overall diameter cannot increase infinitely because of the departure criterion. It can be described as a process that a growing bubble of a certain size of diameter bears the external force which reduces the bubble diameter and the process consists of growing and shrinking at the same time. The specific values of the increment and the decrement depend on the interaction between factors. The non-equilibrium force produces an annular break-up region surrounding that without break-up to decrease the bubble diameter. The break-up region slightly moves up due to the horizon flow and the gravity.

Conclusions
A new bubble diameter model is developed for the flow boiling simulation. The bubble diameter is divided into four components: the departure diameter, the diameter due to heat transfer, the diameter due to collision, and the diameter due to phase interaction force. Based on the proposed bubble diameter model, a flow boiling simulation scheme is established and carried out through C language and Fluent software. Main conclusions can be drawn: (1) The proposed overall diameter model can be regarded as a combination of the initial diameter and contributions from heat transfer, collision and phase interaction force. It is suitable for more operation conditions and can describe the void fraction and the temperature of the flow boiling accurately in different operating pressures and temperatures. (2) Comparing to linear correlations, the proposed model shows a more dispersed distribution of the bubble diameter inside the channel. Seen from the section distribution, the bubbles of large size spread over a wide asymmetric area instead of the evenly annular distribution of linear correlations. (3) Among the components of the overall bubble diameter, the variation due to heat transfer is considered as the dominant part controlling the bubble size. Meanwhile, the bubble diameter is also limited by the phase interaction force. The contribution from the collision is the smallest in the overall bubble diameter.
Author Contributions: Z.X. designed the model and implemented programming; J.Z. and J.L. conducted the simulation; J.Z., T.X., J.W. and Z.L. analyzed the data; and J.L. wrote the paper.