Numerical Simulation Research of Bubble Characteristics and Bubble Departure Diameter in Subcooled Flow Boiling

Abstract

Three-dimensional subcooled flow boiling of R134a in a horizontal tube was simulated by a VOF (volume of fluid) model combined with the level set method. Bubble characteristics were explored at heat flux of 0.3 MW/m², inlet subcooling of 3 K, and inlet velocity of 0.4 m/s. It was observed that five representative bubbles occurred in subcooled flow boiling, including sliding bubble, coalescing bubble, non-departed bubble, bouncing bubble, and continuous-boiling bubble. The results showed that the bubble radial velocity was an important factor of bubble departure after a sliding process. Moreover, the effect of heat flux, inlet velocity, and inlet subcooling on bubble departure diameter were investigated. The departure diameter increased with increasing inlet velocity from 0.2 to 0.4 m/s and heat flux from 0.2 to 0.4 MW/m², while diameter decreased with inlet subcooling from 3 to 10 K. Finally, based on the influence of heat flux, inlet velocity, and inlet subcooling on average departure diameter of the bubble except the coalescing bubble, a model was proposed to predict the average departure diameter. The deviation of the model was within 5%.

1. Introduction

Compared with saturated boiling, subcooled flow boiling shows better heat transfer ability. Bubble generation requires lots of latent heat and transports energy from then hot region to cold bulk liquid in subcooled flow boiling [1,2]. Subcooled flow boiling is applied to many fields, for instance, refrigeration, electronic cooling, aerospace, nuclear, and biological engineering [3,4]. Researchers have paid more attention to study bubbles in subcooled flow boiling, because better understanding and prediction of bubbles are critical for heat transfer mechanism and enhancement in subcooled flow boiling [5].

The visualization is absolutely necessary for studying bubbles in subcooled flow boiling. In recent decades, the process of subcooled flow boiling was observed and investigated through high-speed camera visualization, which made progress in understanding onset of boiling (ONB), bubble growth, bubble sliding, and bubble departure. Yang et al. [6] studied upward subcooled flow boiling in a narrow channel through experiments at mass flux of 122 to 657 kg/(m²·s), subcooling of 4.7 to 33.3 K, and heat flux of 0.17 to 2.89 kW/m². It was concluded that the ONB relied strongly on wall superheat and heat flux. The mass flux, wall superheat, and subcooling had considerable impact on bubble departure diameter. Liu et al. [7] carried out a visual experiment to study bubble size in vertical subcooled flow boiling. It was observed that bubbles slid long distances along the heating wall without increasing size at heat flux of 0.15 MW/m². However, the bubble size increased with sliding at heat flux of 0.3 MW/m². It was found that bubble size increased with increasing heat flux of 0.15 to 0.3 MW/m², while it decreased with increasing system pressure of 0.1 to 1.0 MPa, subcooling of 10 to 30 K, and mass flux of 800 to 900 kg/(m²·s). Yuan et al. [8] researched the bubble sliding process in a vertical rectangular narrow channel in subcooled flow boiling. The working conditions were at mass flux of 117.49 to 300.84 kg/(m²·s), subcooling of 20 to 36 K, and heat flux of 64.57 to 147.82 kW/m². The results showed that both bubble speed and sliding distance increased with increasing mainstream velocity. Cao et al. [1] gave a very detailed description of bubble deformation at a velocity of 0.1 to 0.8 m/s and subcooling of 5 to 50 K. It was found that the very weak deforming phenomena was observed at low subcooling. Zhou et al. [9] established a bubble dynamics model on the basis of visual experimental results, which could predict the bubble departure diameter well. Omar et al. [10] investigated bubble departure in subcooled flow boiling, and a new empirical correlation was proposed to predict bubble departure. There are more researchers exploring bubble through experimental visualization, such as Lie and Lin [11], Fu et al. [12], and so on. As reviewed above, exploring subcooled flow boiling through high-speed camera visualization has been reported widely. However, there are many limiting and uncertain factors in the high-speed visualization experiment; the cost of the experiment is high, the experimental period is long, and experimental results are not controllable.

In recent years, it has been popular to develop numerical solutions to achieve the visualization of subcooled flow boiling [13]. Compared with experiments, the numerical simulation has great advantages, including low cost, good controllability, and short time frame. The software CFD (Computational Fluid Dynamics) is used in the research of subcooled flow boiling because the bubble is tracked by the volume of fluid (VOF) model [14,15,16]. In addition, the level set method also plays a role in bubble tracking [17,18]. Zhuan and Wang [19] simulated subcooled flow boiling in a microchannel by a VOF model; the process of bubble growth, condensation, and collapse were observed. It was concluded that bubble growth, condensation, and collapse were affected by subcooling, departure size, heat flux, and mass flux. At high subcooling, the surface tension acting on a bubble is very great. The flow pattern observed in simulation corresponded to the experiment in a certain range. Guo et al. [20] researched bubble dynamics in a one-side heating wall and built a bubble growth model with a VOF model. The results showed that the fluid field was divided into three zones, including bubble generating zone, growth zone, and annihilation zone, and the three zones changed with bubble motion. Wei et al. [21] used a VOF model to investigate bubbles in subcooled flow boiling under swing motion. Bubble coalescence, sliding, and departure were observed; meanwhile, bubble departure and reattachment were occurring occasionally. Lee et al. [22] simulated bubble growth on a microfinned face by using the level set method. The bubble shape was tracked well. The effect of a fin on bubble growth was studied to find the optimal conditions for heat transfer enhancement. The VOF model is widely used compared with the level set method, which is especially applied to shape changes. In addition, more studies tend to use water as the working fluid. Though researchers have made some progress in understanding subcooled flow boiling through the VOF model or the level set method, simulations are still being conducted.

In this paper, both the VOF model combined and the level set method were adopted in a numerical simulation. In addition, the refrigerant R134a was selected as the working fluid for its wide heat transfer applications. R134a has been widely used in commercial refrigeration, but few studies have examined subcooled boiling of R134a. The study of R134a subcooled boiling is of great significance for improving system performance in the field of low carbon and energy saving for our study, the boiling process of subcooled flow boiling of R134a in a three-dimensional horizontal tube was simulated. The bubble characteristics in subcooled flow boiling were discussed, and, the influence of inlet velocity, heat flux, and subcooling on bubble departure diameter were investigated.

2. Numerical Solution

2.1. Geometrical Configuration

A three-dimensional horizontal tube model was established to simulate the subcooled flow boiling as shown in Figure 1. The diameter of the inner heating wall and the outer adiabatic wall were 10 mm and 16 mm, respectively, and the length was 100 mm. The mesh used in the simulation was a structured O-grid, which ensured the precision and stability of the numerical solution. Considering the minimal value of bubble departure diameter, the height of the first layer of mesh on the heating wall was determined by experimental data.

2.2. The VOF Model

In this paper, the visual simulation of flow boiling was performed by using a VOF model combined with the level set method. The VOF model has the advantage of mass conservation in interface capture aspect, but it is not ideal in a complex interface simulation. The level set method has a strong ability to handle complex interfaces, but the disadvantage is that mass is not conserved. The combination of a VOF model and level set can simulate the bubble characteristics well.

In a VOF model, the interface is tracked by obtaining the volume fraction α in each computational cell [14], and the volume fraction α is,

$$\alpha =\left\{\begin{array}{c}0 \\ 0<\alpha <1 \\ 1 \end{array}\right.\begin{array}{c}\hspace{1em}\hspace{1em}\mathrm{liquid} \mathrm{phase}\\ \mathrm{interface}\\ \hspace{1em}\hspace{1em}\mathrm{vapor} \mathrm{phase}\end{array}$$

(1)

where α = 1 means the computational cell is occupied by liquid. α = 0 means the computational cell is occupied by vapor; α between 0 and 1 represents the vapor–liquid interface.

The continuity equation for the vapor phase is,

$$\frac{\partial {\alpha }_v}{\partial t}+\nabla ·(\overrightarrow{u}{\alpha }_v)=\frac{S_{M-lv}}{{\rho }_v}$$

(2)

The continuity equation for the liquid phase is,

$$\frac{\partial {\alpha }_l}{\partial t}+\nabla ·(\overrightarrow{u}{\alpha }_l)=\frac{S_{M-vl}}{{\rho }_l}$$

(3)

$${\alpha }_l+{\alpha }_v=1$$

(4)

where α_l is the liquid fraction, α_v is the vapor fraction, $\overrightarrow{u}$ is the fluid velocity vector [m/s], ρ_v is the vapor density [kg/m³], ρ_l is the liquid density [kg/m³], S_M-lv is the mass source from liquid phase to vapor phase [kg/(m³s)], t is time [s], and S_M-vl is the mass source from vapor phase to liquid phase [kg/(m³s)].

The velocity field obtained by solving the single momentum equation in the whole domain is shared among the phases. The momentum equation shown below depends on the volume fraction of all phases of the physical properties ρ and μ. The liquid phase and vapor phase share the same set of momentum and energy equations. The momentum equation is,

$$\frac{\partial \rho \overrightarrow{u}}{\partial t}+\nabla ·(\rho \overrightarrow{u}\overrightarrow{u})=-\nabla p+\nabla ·[{\mu }_{eff}(\nabla \overrightarrow{u}+\nabla {\overrightarrow{u}}^T)]+\rho \overrightarrow{g}+{\overrightarrow{F}}_{vol}$$

(5)

The energy equation is,

$$\frac{\partial \rho c_{p}T}{\partial t}+\nabla (\overrightarrow{u}(\rho c_{p}T+p))=\nabla ({\lambda }_{eff}\nabla T)+S_E$$

(6)

where p is the pressure of interface [N/m²], ${\overrightarrow{u}}^T$ is the transpose part of the strain rate tensor, g is the acceleration of gravity [m/s²], μ_eff is the effective dynamic viscosity [Ns/m²], λ_eff is the effective thermal conductivity [W/mK], ${\overrightarrow{F}}_{vol}$ is the surface tension per unit volume of fluid [N/m³], and S_E is the energy source term that is loaded by the UDF (user defined function) [W/m³].

The effective viscosity μ_eff is,

$${\mu }_{eff}=\mu +\rho {\nu }_t$$

(7)

The effective thermal conductivity λ_eff is,

$${\lambda }_{eff}=\lambda +\rho c_p\frac{{\nu }_t}{p{r}_t}$$

(8)

where v_t is the turbulent kinematic viscosity [m²/s], and Pr_t is the turbulent Prandtl number. The Pr_t is obtained by referring to Weigand et al. [23].

In Equations (5) and (6), where ρ is the mixture density [kg/m³], λ is the mixture conductivity [W/mK], $\overrightarrow{u}$ is the mixture velocity [m/s], μ is the mixture viscosity [Ns/m²], T is the mixture temperature [K], and c_p is the mixture specific heat capacity [J/kgK]. They are acquired by the volume average method, that is,

$$\rho ={\alpha }_l{\rho }_l+{\alpha }_v{\rho }_v$$

(9)

$$\lambda ={\alpha }_l{\lambda }_l+{\alpha }_v{\lambda }_v$$

(10)

$$\mu ={\alpha }_l{\mu }_l+{\alpha }_v{\mu }_v$$

(11)

$$T=\frac{{\alpha }_l{\rho }_l{T}_l+{\alpha }_v{\rho }_v{T}_v}{{\alpha }_l{\rho }_l+{\alpha }_v{\rho }_v}$$

(12)

$$c_p={\alpha }_l{c}_{pl}+{\alpha }_v{c}_{pv}$$

(13)

2.3. The Level Set Method

The level set method defines a symbolic distance function ϕ to track the phase interface [24], and the bubble surface or liquid–vapor interface is tracked by the level set function ϕ, which is defined as a signed distance from the interface. The negative sign is chosen for the vapor phase and the positive sign for the liquid phase. In each simulation, the level set looks for ϕ = 0 to track the gas–liquid interface, and the expression equation is,

$$\varphi (x,t)=\left\{\begin{array}{c}d(x,\Gamma (t))\\ 0\\ -d(x,\Gamma (t))\end{array}\right.\begin{array}{c}x\in R\\ x=\Gamma (t)\\ x\in R\end{array}$$

(14)

where d(x, Γ(t)) is the distance from x to Γ(t), and Γ(t) is the zero iso-surface of ϕ(x,t), which is defined as,

$$\Gamma (t)=\left\{x\subset \left.\mathsf{\Omega }\right|\varphi (x,t)=0\right\}$$

(15)

The control equation for vapor based on the level set method is,

$$\frac{\partial \varphi }{\partial t}+\nabla ·(\overrightarrow{u_v}\varphi )=0$$

(16)

The curvature ${\kappa }_{level}$ and unit normal vector ${\overrightarrow{n}}_{level}$ of the interface are,

$${\kappa }_{level}=\nabla \cdot {\overrightarrow{n}}_{level}$$

(17)

$${\overrightarrow{n}}_{level}=\frac{\nabla \varphi }{\left|\nabla \varphi \right|}$$

(18)

Considering numerical difficulty caused by a sudden change of a physical property on the interface, the smooth Heaviside function is,

$$H(\varphi )=\left\{\begin{array}{c}0 \\ \frac{1}{2}\\ 1 \end{array}\right.\left[1+\frac{\varphi }{\epsilon }+\frac{1}{\pi }\sin (\frac{\pi \varphi }{\epsilon })\right]\begin{array}{c}\hspace{1em} \mathrm{if} \phi <-\epsilon \\ \hspace{1em} \mathrm{if} \phi \le \left|\epsilon \right|\\ \mathrm{if} \phi >\epsilon \end{array}$$

(19)

where ε is the thickness of the interface [m].

The density and viscosity smoothed are,

$$\rho \left(\varphi \right)={\rho }_v+({\rho }_l-{\rho }_v)H(\varphi )$$

(20)

$$\mu \left(\varphi \right)={\mu }_v+({\mu }_l-{\mu }_v)H(\varphi )$$

(21)

The fluid momentum equation combined with the level set method is [25],

$$\frac{\rho (\varphi )\partial \overrightarrow{u}}{\partial t}+\rho (\varphi )\nabla ·(\overrightarrow{u}\overrightarrow{u})=-\nabla p+\mu (\varphi )\nabla ·[(\nabla \overrightarrow{u}+\nabla {\overrightarrow{u}}^T)]+\rho (\varphi )\overrightarrow{g}+{\overrightarrow{F}}_{level}$$

(22)

2.4. Surface Tension

The surface tension is the most important force to impact bubble behaviors, including nucleation of small-sized, activation of nucleation cavities, bubble growth, and bubble departure. The continuous surface force (CSF) model takes into account the impact of surface tension, which transfers the surface tension into body force [26]. The CSF model is added to the momentum equation as a source term. The equation of the source term in Equation (5) is defined as:

$${\overrightarrow{F}}_{vol}=\sigma \frac{\rho {\kappa }_l\nabla ·{\alpha }_l}{0.5({\rho }_l+{\rho }_v)}$$

(23)

where σ is the coefficient of surface tension [N/m], k_l is the surface curvature [1/m], which is defined in terms of the divergence of the unit normal vector $\overrightarrow{n_l}$, k_l and $\overrightarrow{n_l}$ defined as,

$${\kappa }_l=\nabla \cdot \overrightarrow{n_l}$$

(24)

$$\overrightarrow{n_l}=\frac{\nabla \alpha }{\left|\nabla \alpha \right|}$$

(25)

the equation of the south term in Equation (23) is,

$${\overrightarrow{F}}_{level}=\sigma {\kappa }_{level}\delta (\varphi )\nabla \varphi$$

(26)

where δ is defined by averages of the Heaviside function, and the Heaviside function is,

$$H(\varphi )=\left\{\begin{array}{l}\begin{array}{cc}\frac{1}{2}(1+\frac{\cos (\frac{\pi \varphi }{\epsilon })}{\epsilon })& \left|\varphi \right|\ge \epsilon \end{array}\\ \begin{array}{cc}0& \left|\varphi \right|\ge \epsilon \end{array}\end{array}\right.$$

(27)

2.5. Phase Changing Model

In the study of pure working fluid boiling simulation, researchers have proposed several phase changing models from different perspectives, such as a dynamic model, a two-fluid model, an energy model, etc. Considering the motion of bubbles, the Lee [27] model was adopted as the phase change model. The Lee model is a dynamic model which is driven by temperature difference and has obvious physical sense. Whether the boiling appears or not depends entirely on the difference between the fluid and the saturated temperature. The specific form of the source item is as follows [27]:

The liquid phase source term is,

$$S_{M-lv}=-\eta {\alpha }_l{\rho }_l\frac{T-T_{sat}}{T_{sat}} (if T>T_{sat})$$

(28)

The vapor phase source term is,

$$S_{M-vl}=\eta {\alpha }_v{\rho }_v\frac{T_{sat}-T}{T_{sat}} (if T<T_{sat})$$

(29)

The energy source term is,

$$S_E=h{S}_M$$

(30)

where η is the phase change mass transfer coefficient [s⁻¹], and T_sat is the saturation temperature [K], h is the latent heat [J/kg].

The phase change mass transfer coefficient is important for convergence of the simulation and the temperature of the interface, which is small for simulation of boiling while large for condensation [21,28], De Schepper et al.’s [29] value $\eta$ = 0.1 was used in this work. The results showed that the temperature of the interface kept at the saturated temperature, and the simulation could converge well.

2.6. Turbulent Flow Model

The transition Reynolds number in the small channel is smaller than that in the conventional channel. Meanwhile, the turbulence model is also suitable for the small Reynolds number due to the disturbance of bubbles [30]. We define the Reynolds number as follows:

$$Re=\frac{\rho Vd}{\mu }$$

(31)

where V is the inlet velocity [m/s], and d is the hydraulic diameter [m].

Considering the shear effect caused by turbulence near the wall during subcooled boiling, it is necessary to choose a turbulent model for the prediction of boiling. The SST k − ω model considers the effect of shear stress and combines the k − ω model in the near-wall low Reynolds number region, which greatly improves the accuracy and stability of turbulence simulation [31]. Therefore, the SST k − ω model was adopted in the simulation, and the transport equation regarding k and ω are,

$$\frac{\partial \rho k}{\partial t}+\nabla \cdot (\rho k\overrightarrow{u})=\nabla \cdot [(\nu +\frac{{\nu }_t}{{\sigma }_k})\nabla \cdot k]+G_k-Y_k$$

(32)

$$\frac{\partial \rho \omega }{\partial t}+\nabla \cdot (\rho \omega \overrightarrow{u})=\nabla \cdot [(\nu +\frac{{\nu }_t}{{\sigma }_{\omega }})\nabla \cdot \omega ]+G_{\omega }-Y_{\omega }+D_{\omega }$$

(33)

where k is the turbulent kinetic energy [m²/s²], ω is the turbulent vortex frequency, ${\sigma }_{\omega }$ is the reciprocal of effective turbulent Prandtl number of the turbulent vortex frequency, ${\sigma }_k$ is the reciprocal of the effective turbulent Prandtl number of turbulent kinetic energy, G_k is the generation term of turbulent kinetic energy caused by time-average velocity gradient, G_ω is the generation term of turbulent vortex frequency, Y_k is the dissipation term of turbulent kinetic energy, Y_ω is the dissipation term of turbulent vortex frequency, and D_ω is the cross diffusion term. For the specific definition of the above content, see reference [32].

The turbulent kinematic viscosity ${\nu }_t$ is,

$${\nu }_t=\rho \frac{k^2}{\omega }\frac{1}{\max [\frac{1}{a^{*}},\frac{S{F}_2}{a_1\omega }]}$$

(34)

where $a_1$ is 0.31, S is the strain rate, a* = 1, and F₂ refers to [32].

The $y^{+}$ of the first cell was controlled to $4 < y^{+} < 20$ by adjusting the cell distance to adjacent walls. According to the characteristics of flow boiling, the bubble departs from the boundary layer that belongs to the buffer layer. The buffer layer region is the region where viscosity and Reynold stress are obvious. Therefore, the first layer cell uses $4 < y^{+} < 20$ to ensure that the near-wall region is within the buffer region [33]. The equation of $y^{+}$ is,

$$y^{+}=\frac{y{u}_{\tau }}{\nu }$$

(35)

where u_τ is the shear velocity [m/s], y is the boundary layer length [m], and v is the kinematic viscosity [m²/s].

2.7. Operation and Boundary Conditions

The R134a was chosen as the working fluid in the simulation, and the saturated temperature was 288.15 K. The thermophysical properties of the working fluid was acquired by consulting REFPROP.

Table 1. The properties of the working fluid.

Fluid	R134a Liquid	R134a Vapor
Reference temperature (K)	288.15	288.15
Density (kg/m3)	1243.4	23.758
Thermal conductivity (W/(m.K))	0.085	0.0128
Viscosity (lbm/(ft.s))	0.0014828	7.5871 × 10−6

shows the properties of the working fluid. The specific heat capacity c_p varied with fluid temperature, so it was set by means of piecewise-linear in the simulation.

In the simulation, the second boundary condition was selected for the inner wall and the outer wall; that is, the given heat flux q was 200, 300, and 400 KW/m², respectively, the outer wall was adiabatic, and the given heat flux was 0 KW/m². In addition, the wall contact angle was 17°, and no slip was along the wall. The inlet boundary was set as mass-flow-inlet. The mass flux was calculated by combining the inlet velocity, fluid density, and cross-sectional area. The inlet velocities were 0.2, 0.3, and 0.4 m/s, respectively. The inlet temperature was set to 3 K, 6 K, and 10 K, respectively. The outlet boundary was the pressure outlet and set to 1 atmospheric pressure.

2.8. Solution Methods

Solution methods including PISO [34] were chosen for calculation, PRESTO! [35] was chosen for pressure interpolation, and geo-reconstruct was chosen for interface reconstruction. A better simulation scheme was obtained after debugging, and more details are shown in

Table 2. More details of setting.

Item	Content
Discrete scheme of momentum and energy equation	Second-order upwind
Turbulent kinetic energy	First-order upwind
Specific dissipation rate	First-order upwind
Level set function	Second-order upwind
Time step	1 × 10−5
Gradient scheme	Green–Gauss node-based method

.

3. Results and Discussion

3.1. Validation of the Numerical Model

It was discovered that bubbles were not regular spheres when departed from theheating wall; therefore, the equivalent departure diameter was used to calculate departure diameter. The bubble departure diameter is defined as the diameter of the bubble measured by a digital image obtained by the CFD software at the departure time. The departure time is when the bubble stays on the heating wall in the previous frame image and leaves the heating wall in the latter frame. The bubble departs between two frames, and the previous frame is used as the bubble departure time. The time interval between frames was 0.5 ms in the simulation. The calculation of departure diameter followed Han and Griffitth [36], and the equation is,

$$D_d=\sqrt{ab}$$

(36)

where a is the bubble width [m], and b is the bubble height [m].

Average departure diameter is,

$$D_{d,\mathrm{ave}}=\frac{1}{n}{\displaystyle \sum D_d}$$

(37)

The validity of the numerical model was tested by experimental data from Chen et al.’s [37] experiment and was carried out at an inlet velocity of 0.4 m/s, inlet subcooling of 3 K, saturated temperature of 288.15 K, and narrow gap of 1 mm. The R134a was used as the working fluid. The visualization of the boiling process was observed by a high-speed camera, and the experiment showed the change of average departure diameter under different heat fluxes. The simulation of subcooled flow boiling was performed in accordance with the above experimental condition. Firstly, the height of the first layer mesh on the heating wall was set to 0.1 mm by referring to minimal departure diameter. The turbulent model met the simulation requirement due to the value of $y^{+}$ which varied from 4 to 16, which belongs to the buffer region. Then, the average bubble departure diameter under different heat fluxes was compared between the simulation results and experimental data as shown in Figure 2. It was found that the average departure diameter from the simulation had a good agreement with the experimental data. It was proved that the numerical solution was effective.

3.2. Characteristics of Boiling Process

A three-dimensional horizontal tube model was established to simulate the subcooled flow boiling as shown in Figure 1. A simulation case was finished to explore bubble characteristics in subcooled flow boiling under heat flux of 0.3 MW/m², inlet subcooling of 3 K, and inlet velocity of 0.4 m/s.

The developing process of 13 bubbles was observed in the case. Furthermore, 5 representative bubbles were investigated, including a sliding bubble, a coalescing bubble, a non-departed bubble, a bouncing bubble, and continuous-boiling bubbles. Figure 3 shows the developing process of the sliding bubble (in the red rectangular frame). The image of the bubble shows that the bubble experienced nucleation, sliding, departure, and collapse. It can be seen from Figure 3 that the bubble nucleated at 6.5 ms. After that, the bubble slid until at 14.5 ms it began to depart and collapse. Obviously, the sliding dominated the growth process of the sliding bubble. Figure 4 shows the variation of size of the sliding bubble; it was found that the sliding bubble size increased with sliding. There is an isothermal line of 288.818 K near the heating surface from Figure 5 indicating that an overheating layer provided the bubble sliding process with an overheating environment. At the same time, the bubble collapsed after departing from heating, which reflects the characteristics of subcooled flow boiling. Figure 6 presents the developing process of a coalescing bubble (in the red rectangular frame). Compared with the sliding bubble, the bubble coalesces with other bubbles during the sliding process. From Figure 6, two adjacent small bubbles were independent at 7 ms, which then coalesced in one bubble at 8 ms. It is inferred that adequate bubble nucleation causes bubble coalescence.

Coalescing bubbles and sliding bubbles are common in subcooled boiling. An occasional bubble behavior also found is the non-departed bubble. Figure 7 exhibits the non-departed bubble (in the red rectangular frame) that does not depart from the heating wall but collapses directly on the wall. This may be related to the local heat transfer difference and flow field difference. The trend of axial and radial velocities of the above three bubbles during sliding is shown in Figure 8. The variation trends of the axial velocity of the abovementioned three bubbles are consistent, while the radial velocity variation is different. The radial velocity of sliding and coalescing bubbles increase during the sliding process, but the radial velocity of the non-departed bubble decreases during sliding process. It is concluded that increasing radial velocity is the direct reason of bubble departure. In addition, the change of surface tension caused by the change of temperature field affects bubble departure. Obviously, it is shown in Figure 5 that the temperature field of the bubble decreases as a whole when the bubble is departing, which weakens the hindrance of surface tension. Surface tension is the main reason for non-departed bubbles. Analysis of Figure 8 shows that the axial bubble velocity decreases during sliding; this deceleration behavior may provide an environment for bubble growth during sliding.

With bubble departure as the time node, the developing process of two kinds of remaining bubbles occurs after bubble departure, including bouncing and continuous-boiling bubbles. The developing process these bubbles is shown in Figure 9 and Figure 10. The bouncing bubble experiences complex changes after departing as shown in Figure 9. The bubble begins to depart at 14.5 ms, but does not collapse immediately after departing. It keeps close to the heating wall and the bubble size increases from 15.5 to 30 ms. However, the bouncing bubble is condensed and its size greatly reduces from 30 to 41 ms. It can be observed that the bouncing bubble reattached to the heating face at 2 ms, and the bubble continued to slide and grow after reattaching. The bouncing bubble was also observed in early research of Okawa et al.’s [38] experiment. Experimental and simulation results show that the bubble size is condensed before reattaching with the heating wall. The surface tension acting on the bubble after condensation is obvious, which made the bubble reattach to the heating wall. Finally, the continuous-boiling bubble is described in Figure 10. At 17.5 ms, the bubble begins to depart from the heating wall, leaving a small amount of vapor attached to the wall, which will grow up to be a bubble. A portion of vapor was also left behind when the bouncing bubble departed in Figure 9, but the continuous boiling did not appear, which may be related to wall temperature and flow field conditions. In the existing theoretical research, a small amount of bubble vapor remaining in the heating wall is also one of the causes of bubble nucleation.

Except for the five bubbles described above, the characteristic of all bubble departure was investigated. Figure 11 shows the bubble departure diameter of 12 bubbles at a heat flux of 0.3 MW/m², inlet velocity of 0.4 m/s, and inlet subcooling of 3 K. As can be seen from Figure 11, the bubble departure diameter is widely in the range of 0.25 to 0.35 mm. However, there are bubbles with diameters above 0.4 mm and below 0.25 mm. This is because different bubble nucleation points have different positions, and the difference between the local heat transfer and the flow field velocity in different positions makes the characteristics of bubble departure quite different. The heat flux, inlet velocity, and inlet subcooling are important factors affecting bubble departure, which are discussed in the next section.

3.3. Bubble Departure Diameter

The complete boiling process of bubbles was observed to explore the influence of different factors on bubble departure diameter. Firstly, the effect of heat flux on bubble departure diameter was investigated. Figure 11 and Figure 12 show bubble departure diameter at different heat flux, inlet subcooling of 3 K, and inlet velocity of 0.4 m/s. The bubble departure diameter is between 0.15 and 0.475 mm when the heat flux is 0.2 MW/m² as shown in Figure 12a, and most bubbles are between 0.2 and 0.35 mm. The bubble departure diameter is in the range of 0.2 and 0.45 mm when the heat flux increases to 0.3 MW as shown in Figure 11, and most bubbles are in range of 0.25 and 0.35 mm. The bubble departure diameter is between 0.25 and 0.5 mm when the heat flux is up to 0.4 MW as shown in Figure 12b, and most bubbles are between 0.3 and 0.35 mm. By comparison, the bubble departure diameter increases as the heat flux increases. Apparently, the heat absorbed by bubbles increases as heat flux increases, which causes the increased departure diameter. In addition, the maximum departure diameter of some bubbles can exceed 0.4 mm, which are attributed to coalescence.

Figure 13a–c illustrate bubble departure diameter at different inlet velocities, simultaneously, keepings the heat flux of 0.2 MW/m² and subcooling of 3 K constant. At the lower inlet velocity of 0.2 m/s, bubbles have departure diameters more than 0.15 mm and less than 0.35 mm. With the velocity raised to 0.3 m/s, bubbles have departure diameters more than 0.2 mm and less than 0.35 mm. Most of the bubbles are in the range of 0.25 to 0.35 mm when inlet velocity is up to 0.4 m/s, and the maximum diameter of the bubble is over 0.4 mm. By comparing departure diameter of the three inlet velocities, it is inferred that departure diameter is increased by larger inlet velocity. It is explained that increasing inlet velocity gives bubbles the energy to obtain larger sliding velocity. Therefore, the bubble with larger velocity slides a longer distance and further increases the bubble departure diameter.

The effect of subcooling on bubble departure diameter is investigated at three inlet subcooling of 3 K, 6 K, and 10 K with constant heat flux of 0.3 MW/m² and velocity of 0.4 m/s, as shown in Figure 14. In low subcooling of 3 K, most bubbles range from 0.25 to 0.35 mm as shown in Figure 14a. In high subcooling of 10 K, most bubbles range from 0 and 0.3 mm as shown in Figure 14c. Under the condition of 6 K, the bubble departure diameter ranges from 0.2 to 0.4 mm. According to the comparison of three inlet subcooling of 3 K, 6 K, and 10 K, the bubble departure diameter decreases with increasing inlet subcooling. The departure diameter decreases with the increase of inlet subcooling due to the condensation of subcooled mainstream. However, several bubbles have large departure diameter at subcooling of 6 K and 10 K. On the one hand, the coalescence of bubbles increases the bubble size during sliding; on the other hand, the influence of inlet velocity and heat flux are prominent for departure diameter.

3.4. The Prediction of Departure Diameter

Considering the difference of departure diameter among bubbles, as well as the complexity of departure diameter caused by the coalescence of bubbles, the average departure diameter of bubbles is used to study predicted correlation. Figure 15 displays the effect of different factors on average departure diameter, heat flux, inlet velocity, and inlet subcooling included. The results show that departure diameter of sliding bubbles increases with increasing heat flux, inlet velocity, and decreasing inlet subcooling. Through Figure 15, on the one hand, it can be seen from a single operating condition (red dotted line or blue dotted line) that the x-axis variable has an effect on the average bubble detachment diameter. On the other hand, it illustrates that heat flux, inlet velocity, and subcooling have a common effect on bubble departure diameter. When the velocity, heat flux, and subcooling increase to a certain value, there will be a certain interaction between them, resulting in either an increase or a slowdown of the trend. Therefore, the dimensionless number Bo, Re, Jα_sub, and D* regarding heat flux, velocity, subcooling, and bubble departure diameter are introduced to investigate the combined effects.

The model of Zeitoun and Shoukri [39] was used to predict bubble departure diameter regarding sliding. The equation of Zeitoun and Shoukri [39] includes the above dimensionless numbers as given in Equation (35). Figure 16 presents the comparison of departure diameter of the simulation with the predicted departure diameter obtained by the Equation (38). It is inferred that the equation cannot predict the bubble departure diameter well due to overpredicting departure size and showing great deviation compared with 5% deviation. The main possible reason is that the equation was developed for water, and the working condition was at mass flux of 156.15 to 403 kg/(m²·s), subcooling of 11.6 to 31.1 K, and heat flux of 286.68 to 603.24 kW/m².

$$D_d^{*}=\frac{0.0683{\rho }^{*1.326}R{e}^{-0.324}}{J{\alpha }_{sub}+\frac{149.2{\rho }^{*1.326}}{B{o}^{0.487}R{e}^{1.6}}}$$

(38)

A model is proposed for predicting the bubble departure diameter. The dimensionless bubble diameter in horizontal subcooled flow boiling of R134a can be correlated as Equation (36). The model proposed predicts the bubble departure diameter of the simulation results and the experiment data very well as shown in Figure 17, and nearly all values are within the deviation of ±5%. The straight line indicates that the predicted value is the same as the experimental or simulated value. Furthermore, the equation is applied to the condition of heat flux ranging from 0.2 to 0.4 MW, inlet velocity ranging from 0.2 to 0.4 m/s, and inlet subcooling varying from 3 to 10 K.

$$D_d^{*}=10.3{\rho }^{*1.658}R{e}^{0.443}B{o}^{0.195}J{\alpha }_{sub}^{0.033}$$

(39)

4. Conclusions

The boiling process of subcooled flow boiling of R134a in a three-dimensional horizontal tube was investigated by a VOF model combined with the level set method at conditions of heat flux from 200 to 300 kW/m², velocity from 0.2 to 0.4 m/s, inlet subcooling from 3 to 10 K, and saturated temperature of 288.15 K. The major conclusions are given as follows.

(1)

Five major bubbles, including sliding, coalescing, non-departed, bouncing, and continuous-boiling bubbles, were observed in a subcooled boiling process at heat flux 0.3 MW/m², inlet subcooling of 3 K, and inlet velocity of 0.4 m/s. The sliding bubble underwent nucleation, sliding, departure, and collapse, and sliding dominated the growth process. The heat flux input provides an overheated environment for bubble growth in the sliding process. The increasing of bubble radial velocity accelerates bubble departure.

(2)

The bubble departure diameter of R134a varied from 0.2 to 0.4 mm at different working conditions. The larger departure diameter was related to coalescence during sliding. The departure diameter increased with increasing heat flux because the bubble absorbs heat from the heating wall during the sliding process. The bubble velocity is subjected to the inlet velocity, and the departure diameter increases with increasing inlet velocity. Because the mainstream temperature affects bubble growth during the sliding process, the bubble departure diameter decreases with increasing inlet subcooling.

(3)

Based on the influence of heat flux, inlet velocity, and inlet subcooling on average departure diameter of bubble except coalescing bubbles, a model is proposed to predict the average bubble departure diameter. The deviation of the model is within 5%, which can be used to predict bubble departure diameter.

(4)

In the future, the developing process and mechanism of coalescing bubbles will be further studied to better understand and control the subcooled boiling process.