Numerical Investigation of Rising Vapour Bubble in Convective Boiling Using an Advanced 3D Hybrid Numerical Method

Abstract

This chapter introduces an advanced and new type of Three-Dimensional (3D) numerical method called the InterSection Marker (ISM) method. The ISM method - a hybrid Lagrangian–Eulerian 3D front-tracking algorithm specifically crafted for multi-phase flow simulation. The method was used to simulate rising vapour bubble behaviour in Convective boiling conditions. Two applications: bubble growth and bubble condensation due to the convective action, were investigated. Numerically obtained bubble properties, such as size, shape and velocity, are compared well against the past works, and the ISM method proved to be an efficient numerical tool for the interface tracking of multi-phase flow CFD simulations involving heat and mass transfer.

1. Introduction

Multi-phase flow simulation, such as vapour bubbles is complicated. Various numerical methods, such as Front-tracking, Marker and Cell, Volume of Fluid, Level Set, Lattice-Boltzmann, were invented to mimic such phenomena – see Figure 1. These methods, however, have their own characteristics strengths and weaknesses (e.g. see summary in [1, 2]). Hybrid methods, as a result, have emerged to harness the merits of their parent methods to improve accuracy and to tackle complex multi-phase flow applications. Combining the strength of the VOF method and the Front-Tracking method, Aulisa et al. [3] Three-Dimensional (3D) method tracks the interface as a Lagrangian but finds the intersection of the surface mesh with control volume faces and locally remeshes the surface contour while preserving the tracked volume. Aulisa’s method, however, requires permanent markers which cannot be seeded or removed after the simulation is executed, and leads to resolution issues for spherical bubble expansion problem. To improve the method prescribed in Aulisa et al. [3], InterSection Marker (ISM) method [4] was devised. The ISM method eliminated the need for permanent markers and addressed the local surface resolution issue in volume inflationary type problems. The ISM method was previously successfully applied to air bubble rise simulations which were adiabatic in nature [5]. However, ISM method’s ability to calculate interfacial area more accurately (uncertainty in the order of 1–2%) than conventional VOF methods proved it an ideal candidate for multi-phase simulations involving heat and mass transfers across the interface, such as rising vapour bubbles in superheated or sub-cooled water. During the simulation, the predicted vapour bubble properties such as size, shape and velocity were compared against the past works and found to be in good agreement.

Figure 1.

Taxonomy of selective numerical methods for multi-phase flows, and the introduction to the ISM method.

2. A brief introduction to the InterSection marker (ISM) method

In the search for higher surface tracking fidelity, the InterSection Marker (ISM) method [4, 5] was developed where the proper determination of the interfacial area is critical, such as for the heat and mass transfers process across the interface separating two-phase fluids. An in-depth description of the ISM method is out of the scope of this chapter, and the reader should consult [4, 5] for details. Below, however, highlight the key features of the ISM method to provide the reader with a basic understanding.

The ISM method uses a Lagrangian surface mesh co-located within a uniform Eulerian mesh where upon flow-field qualities such as pressure, velocity and temperature are calculated. The total surface is modelled as a connected series of discrete interfaces (planar polygons), each located within their own cubic control volume (CCV). Each planar polygon intersects the edges of the control volume, and the combination of cell-edge intersections uniquely identifies the type of polygon a control volume holds.

The ISM method identifies the type of interface residing in a cell by the combination of cell-edge intersections that interface makes. Total of 51 combinations of basic set of planar-type interfaces had been identified: 8 intersection marker combinations for 3 sided interfaces, 15 for 4 sided, 24 for 5 sided, and 4 for 6 sided – see Figure 2. Given the combination of cell edge intersections is unique, a look-up table can be used to identify the type of interface located within each cell [4, 5] in a manner similar to that used in the marching cube method [6]. Further subdivisions of these polygons are carried out to maintain planar surface during translation/deformation – see Figure 3. A triangular tessellation pattern is the preferred option because three points randomly translated will always form a plane. Additional intersection-marker combinations of non-planar-type interfaces were also identified (details in [4, 5]), which are necessary to prevent the modelled interface from collapsing and folding onto itself.

Figure 2.

Planar surfaces co-located within cubic cells can be of 3 to 6 sides.

Figure 3.

Subdivision of planar polygons.

After identification of the planar polygons and their sub-divisions, the next step in the ISM method is to identify the component points of the interface – as shown in Figure 4: (i) the intersection markers where the interface crosses the control volume cell edges, (ii) the cell face conservation points which allow composite curves to be modelled, and (iii) the raised centroid whose position is calculated to satisfy volumetric conservation. The Volume-of-fluid (VOF) is then calculated by summing the volumes of individual triangular columns – refer to Figure 5. For the surface translation and the remeshing process, see Ho et al. [4].

Figure 4.

ISM Interface points.

Figure 5.

VOF calculation.

3. Simulation of rising vapour bubble in convective boiling conditions using the ISM method

Here, two examples of rising vapour bubble simulations due to the convective boiling conditions are illustrated. Firstly, vapour bubble growth in superheated water; secondly, bubble condensation in sub-cooled boiling conditions.

3.1 Vapour bubble growth during rising due to the natural-convective action

In nucleate pool boiling, Vapour bubble typically attains its maximum size at the moment of departure. After lift-off, vapour bubble, however, could continue to grow during its ascent in the presence of favourable condition – a layer of superheated liquid which exists close to the heated surface and is metastable in nature [7]. Although the bubble growth rate in this region is not that significant with compared to the initial pre-departure growth, the bubble will grow for the convective action (convective boiling) with the presence of superheated liquid. Vapour bubble size, shape, and rise velocity for this superheated liquid region can significantly affect the heat and mass transfer mechanism involve. It is thus critical to understand and predict the behaviour and properties of rising and simultaneously growing vapour bubbles. In this section, the growth of a rising vapour bubble is numerically investigated in quiescent superheated water under the influence of buoyancy and surface tension forces with special emphasis given to heat and mass transfers due to the convective action.

3.1.1 Numerical features

Both the water and the vapour phases can be assumed to experience the same ‘mixture velocity’ at any local point within the computational domain, and the two-fluid system can be approximated as a one-fluid mixture. The mixture density and viscosity of each control volume can be calculated based on the volume fraction (α), which has the following values:

$$\alpha =\left\{\begin{array}{c}1\\ 0<\alpha <1\\ 0\end{array}\right.\begin{array}{c}\text{liquid phase}\\ \text{interface}\\ \mathrm{gas}\text{phase}\end{array}$$E1

The variable density and viscosity are then estimated using the α value:

$$\rho =\left(1-\alpha \right){\rho }_g+\alpha {\rho }_l$$E2

$$\mu =\left(1-\alpha \right){\mu }_g+\alpha {\mu }_l$$E3

Where: Subscripts l and g indicate liquid (water) and gas (vapour) phases.

When the mass transfer is considered, a source term needs to be added to the α-transport equation:

$$\frac{\partial \alpha }{\partial t}+\nabla .\left(\mathit{\alpha V}\right)=\left(\frac{1}{\rho }\right)S_{\mathit{\text{mass}}}$$E4

Where: S_mass is the interfacial mass transfer source term.

Similarly, the continuity equation becomes:

$$\nabla .V=\left(\frac{1}{{\rho }_g}-\frac{1}{{\rho }_l}\right)S_{\mathit{\text{mass}}}$$E5

The Momentum equation is:

$$\frac{\partial \mathit{\rho V}}{\partial t}+\nabla \left(\mathit{\rho V}.V\right)=-\nabla p+\mathit{\rho g}+\nabla .\mu \left(\nabla V+\nabla V^T\right)+F_{\sigma }$$E6

Where: p, g and F_σ are the pressure, gravity and surface tension force respectively.

The Energy equation is:

$$\frac{\partial }{\partial t}\left(\rho \mathrm{E}\right)+\nabla \left(V\left(\rho \mathrm{E}+p\right)\right)=\nabla \left(k\nabla T\right)+S_{\mathit{\text{heat}}}$$E7

Where: S_heat, E, T and k are the interfacial heat transfer source term, energy, temperature and thermal conductivity respectively.

S_heat can simply be obtained by multiplying S_mass by the change of enthalpy (h_fg) for phase-change as [8]:

$$S_{\mathit{\text{heat}}}=S_{\mathit{\text{mass}}}\times h_{\mathit{fg}}$$E8

Vapour bubble growth was simulated using the convective heat transfer mechanism and change in enthalpy (phase-change), and the source term can be presented as [9]:

$$S_{\mathit{\text{mass}}}=a_{\mathit{\text{if}}}.\frac{h_{\mathit{\text{if}}}{\mathit{\Delta T}}_{\mathit{\text{super}}}}{h_{\mathit{fg}}}$$E9

Where: $a_{\mathit{\text{if}}}$ is the interfacial bubble surface area per unit volume; hif is the convective heat transfer coefficient; $∆T_{\mathit{\text{super}}}$ is the liquid superheat. Using the ISM method’s capability of calculating interfacial surface area more precisely (than the conventional VOF models), local bubble interfacial surface areas (cell-wise) were used to evaluate the total mass transfer onto the growing bubble.

Thus, Eq. (9) can be written as:

$$S_{\mathit{\text{mass}}}=\sum _n{S}_{\mathit{\text{mass}}\_\mathit{\text{cell}}\left(n\right)}={\sum }_n\frac{h_{\mathit{\text{if}}}\times a_{b_{\mathit{\text{cell}}\left(n\right)}}\times ∆T_{\mathit{\text{super}}}}{h_{\mathit{fg}}}$$E10

Where: $S_{\mathit{\text{mass}}\_\mathit{\text{cell}}\left(n\right)}$ is the 3D spatial interfacial cell-by-cell mass transfer rate onto the growing bubble. $a_{b_{\mathit{\text{cell}}\left(n\right)}}$ is the local bubble surface (interfacial) area at the interface cell and can be obtained from the ISM simulation (see Figures 4 and 5). $∆T_{\mathit{\text{super}}}$ is the liquid superheat. h_if can be calculated by the following evaporation correlation:

$${\mathit{Nu}}_{\mathit{\text{evap}}}=\frac{h_{\mathit{\text{if}}}{D}_b}{k_l}\to h_{\mathit{\text{if}}}=\frac{k_l}{D_b}\times {\mathit{Nu}}_{\mathit{\text{evap}}}$$E11

Evaporative Nusselt number (Nu_evap) can be calculated from the correlations. Nu_evap depends on the mechanism of fluid flow, the properties of the fluid, and the geometry. Numerous heat transfer correlations have been proposed for convective heat and mass transfer from the sphere for specific applications and conditions. To identify an appropriate Nu_evap, selective and widely acceptable correlations were considered (see Table 1), and were plotted against the Bubble Reynolds number (Re_b) [10]. It was found, for small Re_b, there is not much difference among the correlations. However, for higher Re_b, a significant discrepancy exists among the correlations. Hughmark [13] is not only chosen for its broad-application and popularity, but also for providing median range values (not too high or low) for higher Reynolds number. Bubble Reynolds number (Re_b) in the correlation can be expressed as:

Reference	Correlation Proposed (Nuevap = Nu)	Valid For
Ranz and Marshall [11]	$\mathit{Nu}=2+0.6R{e}^{1/2}P{r}^{1/3}$	0 ≤ Re < 200
Whitakar [12]	$\mathit{Nu}=2+\left(0.4R{e}^{1/2}+0.06R{e}^{2/3}\right)P{r}^{0.4}{\left(\frac{\mu }{{\mu }_s}\right)}^{1/4}$	3.5 ≤ Re ≤ 7.6 x 10^4; 0.71 ≤ Pr ≤ 380; 1.0 ≤ (μ/μs) ≤ 3.2
Hughmark [13]	$\mathit{Nu}=2+0.6R{e}^{0.5}P{r}^{0.33}$	0 ≤ Re < 776.06; 0 ≤ Pr < 250
	$\mathit{Nu}=2+0.27R{e}^{0.62}P{r}^{0.33}$	776.06 ≤ Re; 0 ≤ Pr < 200
Akiyama [14]	$\mathit{Nu}=0.37R{e}^{0.6}P{r}^{1/3}$	Laminar Flow
McAdams [15]	$\mathit{Nu}=0.37R{e}^{0.6}$	17 < Re < 70,000

Table 1.

Evaporative Nusselt number correlations [10].

$${\mathit{Re}}_b=\frac{{\rho }_l{U}_b{D}_b}{{\mu }_l}$$E12

Where Bubble rise velocity (U_b) can be calculated as:

$$U_b=\frac{C_z^{n+1}-C_z^n}{∆t}$$E13

Where C_z is the location of bubble centre in upward, z-direction.

Sphere-equivalent Bubble Diameter (D_b) is calculated as:

$$D_b=\sqrt[3]{\frac{6V_b}{\pi }}$$E14

In-a-nut-shell h_if depends on the variables below:

$$h_{\mathit{\text{if}}}=f\left({\rho }_l,{\mu }_l,k_l,\mathit{Pr},U_b,D_b\right)$$E15

First four variables are for water properties at saturation temperature and are constant (for isothermal condition). However, the last two variables are for the vapour bubble and will change continuously for added mass onto the bubble and corresponding varying rise velocity. As such h_if needs to be calculated at each time-step for varying bubble diameter and velocity. As a result, values for the interfacial mass transfer source term (S_mass) will also change in each time step. This demonstrates, even for the isothermal condition, the complex physics behind a growing vapour bubble.

Coupled with an in-house variable-density and variable-viscosity single-fluid flow solver, the ISM interface tracking method was employed to simulate single vapour bubble growth (test sizes 2.5 mm, 3 mm, 4 mm) in quiescent water under the influence of gravity and surface tension forces. Detail descriptions of all the numerical features are not the scope of this chapter. Table 2, however, shows the salient features used during the numerical simulation. Interested readers could get further information from the relevant references.

Description	Method/Mechanism/Platform	Reference
Surface construction, Interface advection and remeshing	ISM method (Hybrid Eulerian–Lagrangian)	[4]
Coupling between ISM interface tracking method and in-house variable-density and variable-viscosity single-fluid flow solver	Immersed Boundary Method (IBM)	[16]
Surface (interfacial) tension	Continuum Surface Force (CSF)	[17]
3D surface curvature	Paraboloid Least Square fitting method	—
Pressure–velocity coupling	Semi-Implicit Method for Pressure-Linked Equations (SIMPLE) algorithm	[18]
Discretisation schemes	Finite volume formulation - hybrid and central schemes	—
Complier for ISM interface tracking algorithm and the flow solver program	Intel Visual ForTran Composer XE 2011	—
CFD result Visualisation	Techplot	—

Table 2.

Numerical features.

Simulations were carried out in a computational domain of 31 × 51 × 31 Cubic Control Volumes (CCV) with an initial spherical bubble of radius 5 h (where h is the width of the non-dimensional cubic control volume) – see Figure 6. Other mesh sizes, such as 21 x 31 x 21 and 41 x 61 x 41, were also investigated, but the mesh size of 31 x 51 x 31 is maintained the same as the previous successful ISM application of Ho et al. [5] to minimise numerical error and optimise computational time. See Ho et al. [4, 5] for ISM fidelity and sensitivity testing. The centre of the bubble was located in line with the centre of the cavity, at a distance of 15.5 h from each side wall and at a distance of 15.5 h from the bottom boundary. All thermos-physical properties were taken at the saturation temperature of 100 °C. To check the effect of liquid superheat on the bubble growth, a wide variety of liquid superheat temperatures, ΔT_super (1 °C, 15 °C, and 35 °C) were considered during the testing. Variable time steps (in the range of $1\times {10}^{-4}$ to $1\times {10}^{-5}$) were also used and found to have insignificant effects on the numerical results. In terms of the computational efforts by using the ISM method, it took 1–3 days to simulate the transient bubble growth process on a personal computer with 2.2 GHz quad-core processor and 16 GB RAM.

Figure 6.

Schematic diagram of test setup (not-to-scale).

3.1.2 Numerical results

3.1.2.1 Size

Normalised bubble volumes of growing bubble over time are plotted in Figure 7. As the interfacial mass transfer source term (S_mass) is directly proportional to the liquid superheat (See Eq. (9)), the bubble was growing at higher rates for larger liquid superheats. For all cases, bubble growth started to vary from 1–2 ms, as in addition to liquid superheat bubble velocity also begun to play a critical role in the bubble growth. Bubble volume growth ratios obtained during the numerical simulation are compared with the theory in Figure 8 and found to be in good agreement. The deviation is less than 1% and is obvious, as fixed h_if values are used in the entire theoretical calculations; on the other hand, h_if values in the numerical simulations are always changing for varying bubble velocities. It is to be noted, normalised bubble volume of a growing bubble due to the convective action can be evaluated analytically as [10]:

Figure 7.

Normalised bubble volume over time (D_b0 = 3 mm) [10].

Figure 8.

Comparison with theory (D_b0 = 3 mm, ΔT_super = 35 °C) (adapted from [10]).

$$\frac{V_b}{V_{b0}}={\left(1+\frac{2h_{\mathit{\text{if}}}∆T_{\mathit{\text{super}}}t}{{\rho }_l{D}_{b0}{h}_{\mathit{fg}}}\right)}^3$$E16

3.1.2.2 Shape

Aspect Ratio (AR) is used to quantify the bubble shape, and is defined by the bubble height by width. AR value of 1.0 indicates the bubble is in a perfect spherical shape. Values less than 1.0 designate the bubble is an oblate spheroid. AR for various bubbles are plotted in Figure 9. With the increase of initial bubble size, bubble deformed at a higher rate and became flattened/oblate spheroid. The liquid superheat also has direct effects on the bubble shape. For the same size bubble, with the increase of liquid superheat, the bubble deformed at a higher rate (i.e. became oblate) due to the increased rate of mass transfer.

Figure 9.

Comparison of aspect ratio, AR (ΔT_super = 1 °C).

Numerically obtained bubble shapes are also compared with the shape regimes of Clift et al. [19] and found to be in good agreement (see Table 3). Here, Eötvös number (Eo) is defined as

Table 3.

Bubble shape validation (cases with selective simulation time are chosen for demonstration). All bubbles are on the same scale for comparison.

$$\mathit{Eo}=\frac{\left({\rho }_l-{\rho }_g\right)g{D}_b^2}{\sigma }$$E17

As Eo is fixed for the same-sized individual bubble, during its ascend the bubble will go through different shape remiges based on its velocity or Re (see [19] for details). For, 2.5 mm, 3 mm and 4 mm bubbles, this transformation is from spherical-to-ellipsoidal-to-wobbling. It is to be noted that the boundary of these shape regimes is not strictly defined, and for the small bubble, such as 2.5 mm bubble could transform from spherical shape to the wobbling direct based on Re. At the later stage, numerical bubble shapes also reveal unstable features and disturbance on the interface. Larger bubble, 4 mm, in this case, became reverse-heart like shape with disturbance on the bottom.

3.1.2.3 Velocity

When an adiabatic bubble is released from a stationary position, its velocity continues to rise until it reaches to its terminal velocity regime where the bubble continues to rise at a constant velocity, i.e. no acceleration [19, 20, 21]. Growing vapour bubble also showed similar velocity profile during its ascend during the numerical simulation [10]. Initially, the velocity rose exponentially after the release of the bubble from its stationary position. The bubble then entered into a relatively stable velocity regime. Generally, with the increase of the bubble size, the bubble reaches to its stable velocity regime more quickly for the larger buoyancy force resulted from larger bubble volume.

The rise velocities (U_stb) of growing bubbles in the stable regime are benchmarked against similar terminal velocity (U_ter) regime of adiabatic bubbles in Table 4. Overall, with the increase of the bubble size, U_stb decreases, and the simulation results compared well with the high-low ranges reported by Clift et al. [19]. ISM’s results follow the lower limit, as for the growing bubbles, U_stb is lower than a comparable adiabatic bubble’s U_ter. This is because the drag force increases at a larger rate due to the deformed frontal area (i.e. bubble shape) of a growing bubble for continually added mass which reduces the bubble rise velocity [22, 23]. Sideman and Taitel’s [24] experiments also exhibited similar trends for evaporated growing bubbles.

Initial Bubble Diameter, D0 [mm]	Liquid Superheat, ΔTSuper [$°C$]	Numerically obtained average rise velocity of a growing bubble in stable regime, UStb [cm/s]	Terminal Velocity, Uter [cm/s] of Adiabatic Bubble
			Experimental (High - Low) [10]
2.5	1	20.28	28–17
	15	19.98
	35	19.60
3	1	20.98	26–17.5
	15	20.12
	35	19.70
4	1	17.12	25–18
	15	16.98
	35	16.94

Table 4.

Bubble velocity comparison (adapted from [10]).

Although small in magnitude, it is observed, with the increase of the liquid superheat, U_stb decreases for a same-sized bubble. This is due to the same reason: with the increase of added mass at a higher rate for higher liquid superheat (see Eq. (9)), the bubble deformed at a faster rate causing more drag force and reduced velocity. Sideman and Taitel’s [24] also overserved similar phenomena of reduced bubble velocity with the increase of temperature.

3.1.3 Conclusion

Below summarised the major findings based on the numerical results obtained using the ISM method:

• As the interfacial mass transfer source term (S_mass) is directly proportional to the liquid superheat, the bubble was growing at higher rates for larger liquid superheats.

• The bubble deforms more with the increase of the size. The liquid superheat also has a direct effect on the bubble shape. For a same-sized bubble, with the increase of liquid superheat, the bubble deforms at a higher rate due to the increased rate of mass transfer.

• When an adiabatic bubble is released from a stationary position, its velocity continues to rise until it reaches to its terminal velocity regime where the bubble continues to rise at a constant velocity (i.e. no acceleration). A growing bubble also shows similar trends. In the stable velocity regime, however, the growing bubble velocity (U_stb) is generally lower than the terminal velocity (U_ter) of an adiabatic bubble. For the continuously added mass, the growing bubble shape deforms more rapidly than an adiabatic bubble causing higher drag force and reduced velocity.

• With the increase of the liquid superheat, U_stb decreases for a same-sized bubble. This is due to the same reason: with the increase of added mass at a higher rate for higher liquid superheat, the bubble deformed at a faster rate causing more drag force and reduced velocity.

3.2 Vapour bubble condensation in sub-cooled boiling condition

Vapour bubble condensation in subcooled liquid where the water temperature is below saturation is an important physical phenomenon [25]. Subcooled boiling flow is one of such examples and can be found in boilers, steam generators, nuclear reactors, and other engineering systems. To optimise the design and to make these systems safe, it is critical to understand and predict the behaviour of bubbles behaviour in subcooled boiling flows. As the presence of vapour bubbles has a significant effect on the heat transfer characteristics of a system as well as pressure drops and flow instability [8, 26]. Vapour bubble life and collapse during condensation can be either inertia (for high liquid subcooling and high Jacob number) or heat transfer controlled (for low liquid subcooling and low Jacob number). Using the ISM method, this section discusses the heat transfer controlled vapour bubble condensation in quiescent water where the bubble reduction rate is longer, and the process is controlled by the heat transfer at the interface. In order to simulate the condensing bubble, the source terms were modelled in the CFD governing equations to account for heat and mass transfers from the bubble. During the simulation, the predicted condensing vapour bubble properties such as size reduction rate, shape and velocity were compared against the past works and found to be in good agreement.

3.2.1 Numerical features

Likewise, bubble growth due to the convective action in the previous section, vapour bubble size reduction (i.e. condensation) was simulated using the same convective heat transfer mechanism due to the temperature gradient between vapour and water phase. Eq. (1)–(8) are also applicable here. S_mass, however, needs to be treated in the opposite way for mass loss from the vapour bubble.

Thus,

$$S_{\mathit{\text{mass}}}=\sum _n{S}_{\mathit{mas}{s}_{\mathit{\text{cell}}\left(n\right)}}=\sum _n\frac{h_{\mathit{\text{if}}}\times a_{b_{\mathit{\text{cell}}\left(n\right)}}\times ∆T_{\mathit{sub}}}{h_{\mathit{fg}}}$$E18

Where: $S_{\mathit{\text{mass}}\_\mathit{\text{cell}}\left(n\right)}$ is the 3D spatial interfacial cell-by-cell mass transfer rate from the condensing bubble; $a_{b_{\mathit{\text{cell}}\left(n\right)}}$ is the local bubble surface area at the interface cell and can be obtained from the ISM simulation (see Figure 4). $∆T_{\mathit{sub}}$ is the liquid subcooling. Interfacial (convective) heat transfer coefficient (h_if) can be calculated by the following condensation correlation:

$${\mathit{Nu}}_{\mathit{\text{cond}}}=\frac{h_{\mathit{\text{if}}}{D}_b}{k_l}\to h_{\mathit{\text{if}}}=\frac{k_l}{D_b}\times {\mathit{Nu}}_{\mathit{\text{cond}}}$$E19

Eq. (12)–(14) can also be used here to determine Reynolds number, Bubble Velocity and Bubble Diameter. Next using this Re_b, Condensate Nusselt number (Nu_cond) can be calculated from the correlations. Table 5 shows some of the notable relations found in the literature. A preference was given to the relations having Jacob number (Ja) – a dimensionless number which is usually used for boiling, evaporation and condensation applications. To check the model’s sensitiveness, Condensing bubble volume reduction with time for various correlations was investigated [31]. The correlations exhibit a wide discrepancy in volume reduction rates, especially with the progress of simulation time. The reasons being these co-relations were developed for specific media and test setups and are valid for a range of parameters only. Warrier et al. [30], for instance, reported an uncertainty of ±25% when they benchmarked their correlation with others.

Reference	Correlation Proposed	Valid For
Zeitoun et al. [25]	${\mathit{Nu}}_{\mathit{\text{cond}}}=2.04{\mathit{Re}}_b^{0.61}{\alpha }^{0.328}{\mathit{Ja}}_l^{-0.308}$	For Steam-Water flow at near atmospheric pressure and for void fraction up to 30 percent.
Kim & Park [27]	${\mathit{Nu}}_{\mathit{\text{cond}}}=0.2575{\mathit{Re}}_b^{0.7}{\mathit{Pr}}_l^{-0.4564}{\mathit{Ja}}_l^{-0.2043}$	Reb = 1000–6000; 18 < Ja < 36; 1.87 < Pr < 2.03; 0.8 mm < D < 6 mm
Lucic & Mayinger [28]	${\mathit{Nu}}_{\mathit{\text{cond}}}=1.46{\mathit{Re}}_b^{0.61}{\mathit{Pr}}_l^{0.33}{\mathit{Ja}}_l^{-0.31}$	${\mathit{Re}}_b\approx 1000-3400$; $\mathit{Ja}\approx 10-30$
Yuan et al. [29]	${\mathit{Nu}}_{\mathit{\text{cond}}}=0.6{\mathit{Re}}_b^{1/2}{\mathit{Pr}}_l^{1/3}\left(1-{\mathit{Ja}}_l^{0.1}{\mathit{Fo}}_{b0}\right)$	For Narrow channel; Reb = 335–1770; Pr$\approx$ 1.7; Ja = 20–60
Warrier et al. [30]	${\mathit{Nu}}_{\mathit{\text{cond}}}=0.6{\mathit{Re}}_b^{1/2}{\mathit{Pr}}_l^{1/3}\left(1-1.2{\mathit{Ja}}_l^{9/10}{\mathit{Fo}}_{b0}^{2/3}\right)$	20 < Reb < 700; 1.8 < Pr < 2.9; 12 < Ja < 100

Table 5.

Condensate Nusselt number correlations.

Considering the scope of current work and relevance, and for validation purposes, Kim and Park [27] is applied during the numerical simulations. Jacob number (Ja) is defined as:

$${\mathit{Ja}}_l=\frac{{\rho }_l{C}_{\mathit{pl}}∆T_{\mathit{sub}}}{{\rho }_g{h}_{\mathit{fg}}}$$E20

Where: $C_{\mathit{pl}}$ is the liquid (water) specific heat.

For condensation, vapour bubble size and corresponding bubble rise velocity will change continuously. As such, h_if needs to be calculated at each time-step for varying bubble size and velocity. As a result, values for the interfacial mass transfer source term (S_mass) will also change in each time step.

Numerical test setups are the same as the previous application. For numerical test cases, 2 mm, 3 mm and 4 mm initial sized bubbles with various liquid subcooling were considered to check the effect of liquid subcooling on the bubble condensation rates. To validate ISM simulation results further, other odd bubble sizes (e.g. 1.008 mm and 4.9 mm) were also considered.

3.2.2 Numerical results

3.2.2.1 Size

The ISM simulation results were validated against past-correlations. In literature, Bubble history ($\beta$) has been expressed as the transient form consisting of the Fourier number (Fo), and is formulated as, for instance [25, 27]:

$$\beta ={\left(1-5.67{\mathit{Re}}_b^{0.61}{\alpha }^{0.328}{\mathit{Ja}}_l^{0.692}{\mathit{Fo}}_0\right)}^{0.72}$$E21

$$\beta ={\left(1-0.6695{\mathit{Re}}_b^{0.7}{\mathit{Ja}}_l^{0.7957}{\mathit{Pr}}^{0.4564}{\mathit{Fo}}_0\right)}^{0.7692}$$E22

Where: Bubble history ($\beta$) which is defined as:

$$\beta =\frac{D_b}{D_{b0}}$$E23

Where: the subscript 0 indicates the initial state. Therefore, D_b0 is the initial bubble diameter, and D_b is the instantaneous bubble diameter.

The Fourier number (Fo₀) is based on the initial bubble size and is written as:

$${\mathit{Fo}}_0=\frac{\mathit{at}}{D_0^2}$$E24

Where: a is the thermal diffusivity.

Bubble history ($\beta$) obtained during the numerical simulation is plotted against past correlations in Figure 10. Depending on the bubble size, liquid subcooling and Reynolds number (Re), the bubble size reduced at varying rates. The bubble was condensing at a higher rate for larger liquid subcooling and higher Reynolds number. The ISM simulation, however, shows discrepancy at the beginning of the condensation stage. The reason being: fixed Reynolds numbers were used to calculate the correlated $\beta$values. The interfacial (convective) heat transfer coefficient (h_if) is proportional to the bubble velocity, so higher the bubble velocity higher the bubble condensation rate. As the bubble was released from a stationary position (with zero velocity and Reynolds number), it took some time for the bubble to reach a reasonable rise velocity (similar to the terminal velocity regime of an adiabatic bubble) and corresponding Reynolds number. Numerical results show close agreement with past correlations at the later stage of condensation after the bubble achieved relatively higher velocity values.

Figure 10.

Bubble history comparison between correlations and ISM simulation [31]. [D_b0 = 4 mm, ΔT_sub = 10 °C, Re = 500, Ja = 30].

The ISM simulation results were also compared against past experimental results – see Figure 11. The ISM simulation results show good agreement, and the overall bubble condensation trends closely follow other benchmarked case of Kim and Park [27] (deviation in the range of 5–15%). The differences between the experimental and ISM numerical simulation are noticeable for different initial test conditions, bubble shapes, and so on.

Figure 11.

Comparison of normalised bubble volume [31]. [D_b0 = 4.9 mm, ΔT_sub = 12 °C].

3.2.2.2 Shape

Generally, small bubble keeps its spherical shape during its rise for high surface tension forces. Bubble shape deformation accelerates with the increase of its size and becomes ellipsoidal or spherical cap or wobbling. Since 2 mm and 3 mm condensing bubbles quickly loss their mass and become small, their shape is generally limited to spherical or ellipsoidal. A larger bubble, e.g. 4 mm condensing bubble, on the other hand, is going through a series of interesting shape evolution, as such, it is considered to compare against Clift et al. [19] shape regimes. A series of instantaneous bubble status points were considered, and their corresponding Re were calculated using bubble velocity obtained from the ISM simulation for comparison. With the increase of bubble velocity and it corresponding higher Re, the bubble was deforming from spherical to ellipsoidal to wobbling shape regimes. Figure 12 demonstrates the bubble shapes obtained from the ISM simulation have excellent agreement with Clift et al. [19] shape regimes.

Figure 12.

Bubble shape validation with Clift et al. [19]. Test case: D_b0 = 4 mm, ΔT_sub = 25 °C, Eötvös number, Eo = 2.55 (adapted from [31]). All bubbles are on the same scale for comparison.

Bubble shapes obtained in the ISM simulation were also compared with Kamei and Hirata’s [32] experimentation and found to be a good agreement – see Table 6. The condensing bubble was keeping its spherical shape because of small size and high surface tension forces during its rise. Table 6 also shows the comparison of shapes with past numerical works (Zeng et al. [33] and Samkhaniani and Ansari [34]).

Table 6.

Bubble shape comparison between past experimental/numerical results and ISM simulation (D_b0 = 1.008 mm, $∆T_{\mathit{sub}}$ = $25°\mathrm{C}$) (adapted from [31]).

3.2.2.3 Velocity

Likewise growing bubble, the rise velocity of condensing bubbles is different from adiabatic bubbles [8]. For continuous reduction in bubble size (i.e. volume) in subcooled boiling flow condition, bubble rise velocity and shape are also always changing. From Figure 13, it is evident that with the increase of liquid Subcooling bubble rise velocity continues to increase. The findings are consistent with [8] numerical results. The deviation is the result of different test setups; however, Figure 13 overall indicates the trends of higher the liquid Subcooling higher the bubble rise velocity. Bubble buoyancy force decreases for continuous reduction of bubble volume. The drag force is also reduced for smaller bubble frontal area; however, the resulted net effect is positive buoyancy force acting upwards, and the bubble rise velocity increases continuously.

Figure 13.

3 mm bubble rise velocity comparison with Jeon et al. [8]. (as demonstrated in [31]).

3.2.2.4 Effects of fluid flow and varying temperature fluid field

Effects of fluid flow and varying temperature fluid field on bubble condensation were also investigated. With the increase of fluid flow velocity, the bubble was condensing at a higher rate (see Figure 14). This is because the relative bubble velocity increases for same-directional (positive) fluid flow velocity resulting in the higher mass transfer or mass loss from the bubble (see Eq. (18)–(19)). Table 7 depicts the effects of fluid flow velocity on the bubble shape. With a rapid mass loss for a higher fluid flow field, the bubble was deforming at a higher rate and becoming unstable with relatively shorter life span.

Figure 14.

Bubble condensation rates for various fluid flow velocity [D_bo = 2 mm, ΔT_sub = 5 °C].

Table 7.

Condensing bubble shape comparison for various fluid flow velocity [D_bo = 2 mm, ΔT_sub = 5 °C]. All bubbles are on the same scale for comparison.

For varying temperature fluid field, two models were considered: (i) linear and (ii) exponential – see Figure 15. The Linear model can be applied to thermal stratification of hot water tanks and the natural systems, such as lakes and ponds. The exponential model is more suitable for the convective boiling application and hence was applied to current numerical study. Because of the computational domain and the relative higher heat transfer coefficient (h) value for convective boiling condition (h = 8,000 W/m²$°C$ was used in Figure 15), the temperature of fluid field is rapidly achieving the liquid bulk-temperature. As a result, the effects of varying temperature field, in this instance, were minimal on both condensing bubble size and shape. See Figure 16 and Table 8.

Figure 15.

Temperature models.

Figure 16.

2 mm bubble condensation rates [D_bo = 2 mm, ΔT_sub = 5 $°C$].

Table 8.

Condensing bubble shape comparison for varying temperature field [D_bo = 2 mm, ΔT_sub = 5 $°C$]. All bubbles are on the same scale for comparison.

3.2.3 Conclusion

Following critical observations were observed during Condensing bubble simulation:

• Bubble reduction rate depends on the liquid subcooling and bubble velocity. Higher the liquid subcooling higher the bubble reduction rate.

• Same is true for the bubble velocity. The interfacial (convective) heat transfer coefficient (h_if) is proportional to the bubble velocity. As such, higher the bubble velocity higher the bubble condensation rate.

• Precaution should be taken during the selection of an appropriate Condensate Nusselt number (Nu_cond) correlation. As the numerical results can vary significantly and should be chosen based on the applications.

• Bubble deforms at a higher rate for larger bubble sizes during ascend. For higher surface tension force, smaller bubbles tend to keep their original spherical shape.

• Liquid subcooling also contributes to the bubble shape determination as heat and mass transfers across the interface cause distortion. With higher liquid subcooling, bubbles get smaller quickly and maintain spherical shape.

• When an adiabatic bubble is released from a stationary position, the rise velocity of the bubble continues to rise until it reaches a stable velocity regime. Whereas for condensing bubble, the rise velocity is always changing for continual heat and mass transfers.

• With the increase of liquid subcooling, rise velocity of the condensing bubble increases.

• With the increase of fluid flow velocity, both bubble condensation rates and shape deformation increase.

• Varying temperature field in this particular instance had minimal effects on condensing bubble size and shape.

4. Summary

The InterSection Marker (ISM) – a new type of 3D interface tracking method was used to simulate the evaporative growth and condensation of a rising vapour bubble due to convective boiling conditions. The ISM method’s ability to calculate interfacial area more accurately than conventional VOF methods proved it an ideal candidate for multi-phase flow simulations involving heat and mass transfers across the interface. During the simulation, the predicted vapour bubble properties such as size, shape and velocity were compared against the past works and found to be in good agreement.

Conflict of interest

The authors declare no conflict of interest.

Nomenclature

Romana

Thermal diffusivity, m²/s

a_if

Interfacial area between phases per unit volume, m²/m³

c_p

Specific heat, kJ/kg$°\mathrm{C}$

C_z

Location of bubble centre in the z-direction

D_b

Sphere-equivalent Bubble Diameter, m

Eo

Eötvös number, −

Fo

Fourier number, −

g

Gravity, m/s²

h_if

Interfacial (convective) heat transfer coefficient, W/m²$°\mathrm{C}$

h_fg

Enthalpy for vaporisation, kJ/kg

Ja

Jacob number, −

k_l

Thermal conductivity, W/m$°\mathrm{C}$

Nu

Nusselt number, −

p

Pressure, Pa

Re

Reynolds number, −

S_heat

Interfacial heat transfer source term, W/m³

S_mass

Interfacial mass transfer source term, kg/m³s

t

Time, s

ΔT_sub

Liquid subcooling, $°\mathrm{C}$

ΔT_super

Liquid superheat, $°\mathrm{C}$

U_b

Bubble rise velocity, m/s

U_ter

Bubble terminal velocity, m/s

U_stb

Evaporating bubble velocity in the stable regime, m/s

V_b

Bubble volume, m³

Greek symbolsα

Volume fraction

β

Bubble History

μ

Dynamic viscosity, kg/ms

ρ

Density, kg/m³

σ

Surface tension force, N/m

Subscripts0

Initial/Original

b

Bubble

cell

Cell

cond

Condensation

evap

Evaporation

if

Interfacial

g

Gas (vapour) phase

heat

Heat

l

Liquid (water) phase

mass

mass

stb

Stable

ter

Terminal

Acronym3D

Three-Dimensional

AR

Aspect Ratio

CCV

Cubic Control Volumes

CFD

Computational Fluid Dynamics

IBM

Immersed Boundary Method

ISM

InterSection Marker

SIMPLE

Semi-Implicit Method for Pressure-Linked Equations

VOF

Volume-of-fluid