Simulation of Single Vapor Bubble Condensation with Sharp Interface Mass Transfer Model

Abstract

Pure numerical simulation of phase-change phenomena such as boiling and condensation is challenging, as there is no universal model to calculate the transferred mass in all configurations. Among the existing models, the sharp interface model (Fourier model) seems to be a promising solution. In this study, we investigate the limitation of this model via a comparison of the numerical results with the analytical solution and experimental data. Our study confirms the great importance of the initial thermal boundary layer prescription for a simulation of single bubble condensation. Additionally, we derive a semi-analytical correlation based on energy conservation to estimate the condensing bubble lifetime. This correlation declares that the initial diameter, subcooled temperature, and vapor thermophysical properties determine how long a bubble lasts. The simulations are carried out within the OpenFOAM framework using the VoF method to capture the interface between phases. Our investigation demonstrates that calculation of the curvature of interface with the Contour-Based Reconstruction (CBR) method can suppress the parasitic current up to one order.

1. Introduction

Bubble condensation is an essential phenomenon for the description of heat and mass transfer in subcooled boiling. It is encountered in many industrial applications, such as microreactors or microchannels, where the bubble dynamics influence the cooling capacity and introduce instabilities to a robust operating condition [1,2]. The size and the shape of vapor bubbles change continuously during the condensation process, and this phenomenon significantly affects the flow structure. In order to better understand subcooled boiling, it is vital to obtain extensive knowledge of the condensing bubbles’ behavior.

Even though many experiments have been conducted on this topic [3], they are still limited to specific liquid properties or specific operating conditions. Moreover, experimental studies based on visualization with high-speed videography or PIV capture the bubble shape evolution but rarely provide detailed information on flow quantities such as temperature and pressure field. In the past decades, the numerical solvers evolved to provide more detailed information on the interface evolving phenomena such as vapor condensation. In reality, the interface between two phases is not sharp; it has a finite width, where the thermophysical properties change smoothly [4,5]. The interface thickness varies with temperature [6] or pressure [7] but remains in the range of a few nanometers [8]. It is not feasible to capture this thin transition region with conventional grid-based methods. Particle-based methods such as molecular dynamic (MD) can be employed to understand the process of formation of an equilibrium liquid-gas interface at the microscopic level [9,10]. The detailed physical insights obtained from MD simulations during phase-change [11,12,13] are utilized to apply appropriate boundary conditions or initial conditions around the liquid-gas interface in meso-scale or macro-scale simulations. Correspondingly, the interfacial region fits into one computational cell in these simulations, so the interface is assumed to be sharp. There are two general approaches for macro-scale simulations: two-fluid (Euler-Euler) methods and one-fluid methods [14]. The Euler-Euler method divides the gas-liquid into the continuous gas phase and the continuous liquid phase, and for each phase, the mass, momentum, and energy conservation equations are solved separately. It is the most practical method to simulate a two-phase system on a large macroscopic scale, as it requires lower grid resolution and is computationally less expensive compared to one-fluid methods [15]. However, this method is not in closed form and needs additional interaction terms depending on the two-phase flow regimes present. With different relationships under different flow regimes, it might be difficult to obtain an accurate solution for realistic scenarios [14,15]. The one-fluid method is sometimes termed direct numerical simulation for a two-phase system, as no assumptions or models for the interface shape are employed [16,17]. The one-fluid methods are classified into: sharp interface, such as Level-Set (LS) [18] and Volume-of-Fluid (VoF) [19,20], and diffuse interface, such as phase-field [21]. The one-fluid methods are applied for the simulation of a single vapor bubble condensation widely, but the main challenge in such numerical simulations is a reliable and physical computation of the mass transfer through the interface. Unfortunately, the mass transfer models are often based on semi-empirical correlations and hence are rarely generally applicable.

The Lee model [22] is one of the most popular mass transfer models. It assumes that mass is transferred at a constant pressure in a phase change flow system, and the model is derived for a quasi-thermo-equilibrium state:

$${\dot{m}}^{‴}=r_c{\alpha }_l{\rho }_l\frac{T-T_{sat}}{T_{sat}},\mathrm{for}\mathrm{condensation}\left(T<T_{sat}\right).$$

(1)

where $\alpha$(-) is the phase volume fraction and $\rho$ (kg/m³) is density. The subscripts l and g represent the liquid and vapor phases, respectively. The volumetric mass flux ${\dot{m}}^{‴}$ (kg/m³s) depends highly on the relaxation parameter $r_c$ (s⁻¹). A wide range between $0.1$ and ${10}^6$ s⁻¹ is proposed and successfully used for $r_c$ in previous studies [23]. Li et al. [24] derived a correlation for the relaxation parameter and showed the dependence of $r_c$ on the temperature, physical properties, and phase volume fraction of the grid element.

Another widespread model was derived by Tanasawa [25] based on the Schrage phase change model [26]. Schrage computed the interfacial mass flux ${\dot{m}}^{″}$ (kg/m²s) using Hertz–Knudsen equation assuming a jump in the temperature and pressure across the interface $T_{sat}\left(p_l\right)=T_l\ne T_{sat}\left(p_g\right)=T_g$:

$${\dot{m}}^{″}=\frac{2}{2-{\gamma }_c}\sqrt{\frac{M}{2\pi R}}\left[\frac{{\gamma }_c{p}_g}{\sqrt{T_g}}-\frac{{\gamma }_e{p}_l}{\sqrt{T_l}}\right],$$

(2)

where $R=8.314$ (J/mol K) is the universal gas constant, M (kg/mol) is the molar mass, and $\gamma$ is the fraction of molecules transferred from one phase to the other. The subscripts c and e refer to condensation and evaporation, respectively. ${\gamma }_c$ = 1 means all vapor molecules hitting the interface are converted to liquid. In the numerical simulation, usually, ${\gamma }_c={\gamma }_e$ is considered. Tanasawa assumed the interface is at saturation temperature and the mass flux varies linearly with temperature difference and the bulk temperature. He simplified Equation (2) to:

$${\dot{m}}^{″}=\frac{2\gamma }{2-\gamma }\sqrt{\frac{M}{2\pi R}}\frac{{\rho }_g{h}_{lg}\left(T-T_{sat}\right)}{T_{sat}^{3/2}},{\dot{m}}^{‴}={\dot{m}}^{″}\frac{A}{V},$$

(3)

where $h_{lg}$ (J/kg) is the latent heat, V (m³) is the cell volume, and A (m²) is the interfacial area in each cell obtained from the interface reconstruction technique (see Section 2.2). In Equations (2) and  (3), the computed mass flux depends on the empirical parameter $\gamma$. The $\gamma =0.1-1$ is suggested for dynamically renewing water surfaces such as jets and moving films and $\gamma <0.1$ for stagnant surfaces [27]. Samkhaniani and Ansari [28] report that the bubble lifetime is highly sensitive to the choice of $\gamma$ in vapor condensation simulations and suggest that an appropriate value must be selected for simulations in comparison with experiments.

A wider list of available mass transfer models is given in ref. [29]. Almost all available models suffer from a dependence on tuning parameters. In the present study, the sharp interface model is employed where the mass transfer is calculated based on a heat flux balance using the Fourier equation:

$${\dot{m}}^{″}=\frac{q^{″}}{h_{lg}}=\frac{k_l\nabla \mathrm{T}-k_g\nabla \mathrm{T}}{h_{lg}},{\dot{m}}^{‴}={\dot{m}}^{″}\frac{A}{V},$$

(4)

where k (w/m K) is the thermal conductivity. The main objective of the present study is to investigate the sharp interface model predictive capability for vapor bubble condensation in comparison with experiments and the Tanasawa model. Moreover, the model limitations are reported.

2. Methodology and Validation

Both liquid and vapor are treated as an incompressible and immiscible Newtonian fluid mixture. The interface between the two phases is resolved using the volume-of-fluid method (VoF) in OpenFOAM solver interFoam, which has been extended with the CBR method [30]. The reconstruction part is essential for an accurate mass flux rate and interface curvature calculation. This improves the surface tension estimation and reduces the parasitic current [31] by up to one order of magnitude (see Section 2.2). The two-phase Navier–Stokes equation in a single-fluid formulation is solved within the PIMPLE (merged PISO-SIMPLE) algorithm loop. The solver is extensively applied for simulations of boiling [30,32,33] and recently ported to OpenFOAM-6.

The following governing equations are utilized:

• Mass conservation

  $$\frac{\partial \rho }{\partial t}+\nabla ·\left(\rho \mathbf{u}\right)={\dot{m}}^{‴},$$

  (5)
• Momentum conservation

  $$\frac{\partial \rho \mathbf{u}}{\partial t}+\nabla ·\left(\rho \mathbf{uu}\right)-\nabla ·\left[\mu \left(\nabla \mathbf{u}+{\nabla \mathbf{u}}^T\right)\right]=-\nabla p+\rho \mathbf{g}+\kappa \sigma \nabla {\alpha }_l,$$

  (6)
• Energy conservation

  $$\frac{\partial \rho c_{p}T}{\partial t}+\nabla ·\left(\rho c_{p}T\mathbf{u}\right)-\nabla ·\left(k\nabla T\right)={\dot{m}}^{‴}{h}_{lg},$$

  (7)
• Phase-fraction transport equation

  $$\frac{\partial {\alpha }_l}{\partial t}+\nabla ·\left({\alpha }_l\mathbf{u}\right)={\dot{m}}^{‴}\frac{{\alpha }_l}{\rho }.$$

  (8)

Here ${\alpha }_l=\frac{V_l}{V_{cell}}$ is the liquid volume fraction, while the physical properties such as density, heat capacity, thermal conductivity, and viscosity $\theta \in \left[\rho ,c_p,k,\mu \right]$ are estimated with linear interpolation $\theta ={\theta }_l{\alpha }_l+{\theta }_g(1-{\alpha }_l)$ in the interfacial region. The last term in the momentum equation, Equation (6), represents the surface tension force between two phases using a continuous surface force (CSF) model [34], where the interface curvature $\kappa$ is obtained from the CBR method. In the present study, the volumetric mass transfer ${\dot{m}}^{‴}$ due to condensation is calculated with the Fourier model from Equation (4).

2.1. Stefan Problem

The classical one-dimensional Stefan problem is applied to validate the phase change models. The problem setup has been extensively used for the validation of mass transfer models [20,35,36]. The schematic of the problem and the comparison with the analytical solution are shown in Figure 1. The exact solution for interface position $x_i$ (m) is calculated as:

$$x_i\left(t\right)=2\eta \sqrt{d_{g}t},$$

(9)

where t (s) is time and $d_g=k/\rho c_p$ (m²/s) is the vapor thermal diffusivity, and $\eta$ is obtained from:

$$\eta exp\left(\eta \right)erf\left(\eta \right)=\frac{c_p\left(T_{wall}-T_{sat}\right)}{\sqrt{\pi }{h}_{lg}}.$$

(10)

The details of the test-case setup, including boundary conditions and the thermophysical properties, correspond to the ones reported in the previous study by Samkhaniani and Ansari [23]. The comparison with the analytical solution confirms that the sharp interface method is capable of accurately predicting the analytical solution.

2.2. Parasitic Current

The spurious (parasitic) current is described as nonphysical velocity created near the interface region as a result of the numerical imbalance between surface tension force and pressure gradient force. This current is intensified where the surface tension becomes dominant.

The strength of the parasitic current is often measured with the maximal magnitude of velocity $\left|u\right|$ for a single bubble placed in a stagnant liquid in zero gravity in the absence of phase change. The current may distort the interface [31] or impair the mass transfer estimation during phase change [28], rendering nonphysical results. Thus, it is crucial to improve surface tension modeling to avoid the parasitic current.

The conventional OpenFOAM VoF solver for simulations of the incompressible immiscible two-phase systems is interFoam [37]. It is an algebraic VoF method that does not reconstruct interface geometrically. The solver employs the Multidimensional Universal Limiter with Explicit Solution (MULES) technique with an additional interface compression term $\nabla ·\left({\alpha }_l(1-{\alpha }_l){\mathbf{u}}_c\right)$ to Equation (8) to keep the interface sharp within 2 or 3 computational cells [38]. The compression velocity ${\mathbf{u}}_c$ is calculated in the normal direction to the interface. This extra term acts only in the interfacial region to suppress numerical smearing in the $\alpha$-field. Recently, the geometric VoF methods, such as the iso-advector method [39,40], and some variants of Piecewise the Linear Interface Calculation (PLIC) method [41], like Multicut Piecewise-Linear Interface Calculation (MPLIC) [42], have been implemented in the OpenFOAM framework. Those methods geometrically reconstruct the interface on polyhedral mesh, which improves the interface sharpness and provides a more accurate curvature estimation. In the present study, the CBR method [30] reconstructs the interface at iso-surface ${\alpha }_l=0.5$, which is then used only for the calculation of curvature $\kappa$, not for the phase flux calculation correction.

In order to investigate the performance of various VoF methods regarding spurious current, a 2-dimensional gas bubble with the diameter $D_0=2$ mm is placed in the centre of the computational domain with the size of $2D_0\times 2D_0$ and filled with a quiescent liquid. Uniform hexahedral cells ($\Delta x=D_0/100$) are utilized for domain discretization. The liquid and gas physical properties are ${\rho }_l=1000$ kg/m³, ${\nu }_l={\mu }_l/{\rho }_l={10}^{-6}$ m²/s, ${\rho }_g=1$ kg/m³, ${\nu }_g={\mu }_g/{\rho }_g={10}^{-6}$ m²/s, $\sigma =0.1$ N/m. The pressure value at the boundaries is considered as uniformly constant and equal to zero and the Neumann boundary condition is assigned for velocity.

The spurious current contour is displayed in Figure 2 for four considered VoF methods. In geometric VoF methods (case B and D), the parasitic current influences a thicker region around the interface and the bubble deviates from its initial spherical shape. As the density of the gas is much smaller than the liquid, the imbalance force accelerates the gas, and a higher velocity is observed on the gas side. The magnitude of the parasitic current in the CBR method is at least one order of magnitude smaller compared to the other methods provided by the standard interFoam implementation, as shown in Figure 3, where the time history of the maximal magnitude of the parasitic current is plotted.

3. Results and Discussion

In this section, a single vapor bubble of condensation is simulated with the Fourier model. The bubble lifetime is compared with experiments and previous numerical simulation based on the Tanasawa mass transfer model. In conclusion, a semi-analytical correlation for the bubble’s life is proposed.

3.1. Problem Definition

The rising of a single vapor bubble is simulated similar to [28,43]. The bubble is introduced at the saturated temperature $T_{sat}=380.2$ K at $P_{sat}=0.13$ MPa and is surrounded by quiescent water at $T_{inf}=355.2$ K corresponding to a 25 K subcooling temperature. The thermophysical properties for vapor are $\rho =0.754$ kg/m³, $\nu =1.66\times {10}^{-5}$ m²/s, $k=0.0259$ W/m K, $c_p=2110.7$ J/kg K, and for liquid water $\rho =953.1$ kg/m³, $\nu =2.75\times {10}^{-7}$ m²/s, $k=0.68$ W/m K, $c_p=4224.4$ J/kg K. The surface tension $\sigma =0.057$ N/m, latent heat $h_{lg}=2237$ kJ/kg, and gravity $\mathbf{g}=9.81$ m/s² are used [43,44]. The computational domain size is $2D_0\times 4D_0$ and filled with $100\times 200$ uniform hexahedral cells, which corresponds to 50 cells per diameter. The initial diameter of the vapor bubble is $D_0=1.008$ mm, located at the middle line with distance $D_0$ from the bottom patch. There is a thin thermal region around the interface where the temperature smoothly changes from a saturated temperature inside the bubble to the subcooled temperature of the surrounding liquid. The temperature profile inside this region is initialized with [45]:

$$T\left(r\right)=a{r}^2-b\sqrt{\left(r\right)}+c,a=\frac{T_{sat}-T_{inf}}{{\delta }^2},b=2a\left(\frac{D_0}{2}+\delta \right),c=T_{inf}+a{\left(\frac{D_0}{2}+\delta \right)}^2,$$

(11)

where r (m) is the distance to the bubble centre and $\delta$ (m) is the thickness of the thermal region. At the boundaries, the uniform constant dynamic pressure ($p_d=0$ Pa) and temperature ($T=T_{inf}$) are introduced, while the velocity gradient and volume fraction gradient are set to be zero.

3.2. Validation

The bubble shape sequence is compared with an experiment from [43] and a previous numerical simulation [28] in Figure 4, and the result shows qualitatively good agreement. The vapor bubble condenses while moving upward and accelerates as it becomes smaller.

For quantitative comparison, the bubble’s lifetime is plotted against experimental data in Figure 5. Simulation of bubble evolution with the Fourier model for $\delta$∼$0.2D_0$ coincides well with the experimental data [43] and previous numerical simulation [28].

3.3. Thermal Boundary Layer

The bubble’s initial temperature profile is shown in Figure 5a. The bubble lifetime is influenced by this thin thermal region around the interface, as given in Figure 5b. It shows that when no thermal boundary region is defined, the bubble collapses too fast. In contrast, a thicker thermal boundary layer around the interface slows down the process and prolongs the bubble’s lifetime. The results indicate the importance of the sub-millimeter region around the interface for the bubble condensate rate. It is the most influencing parameter in the considered simulation with the Fourier model. Unfortunately, experimental studies often do not provide any information about this thermal boundary layer. A comparison with experimental data without the knowledge of the exact temperature field in the vicinity of the interface in the initial state of a simulation is a trial-and-error process since an accurate prescription of the temperature distribution in this tiny region is rather difficult.

In the present numerical simulation, a pure substance fluid is considered. Obviously, in real-world scenarios, there will practically always be some impurities involved, e.g., nitrogen, oxygen and carbon dioxide solved in the water. In such mixtures, a light boiling component (non-condensable gas) will usually accumulate at the interface [46,47]. This accumulation may act as a barrier to mass transfer [10,11,48]. It has been shown that even very small mole fractions of a non-condensable gas solved in the liquid might change the overall dynamics of the system and reduce the mass transfer rate significantly [49,50]. Thus, additionally to the thermal boundary layer, the non-condensable gases can hinder the mass transfer rate. As an alternative option, Tanasawa mass transfer model might be recommended since the mass coefficient $\lambda$ can indirectly account for the missing thermal boundary layer information and/or non-condensable gas effect. The effect of these gases might be considered in future work, but it remains out of the scope of the present study.

3.4. Bubble Lifetime

The lifetimes of bubbles at various subcooled temperatures $\Delta T_{sub}=\left[5-100\right]\mathrm{K}$ obtained from numerical simulations are plotted in Figure 6. The trend is similar for both diameters $D_0=4,16$ mm. The mass transfer rate is driven by the temperature difference between the vapor bubble at saturation temperature and the bulk liquid in the subcooled temperature. Therefore, at low subcooled temperatures, the bubble is maintained for a longer period of time. However, there is a limitation to the total heat capacity that can be transferred in a definite portion of time using the conduction and convection modes. Thus, at higher subcooled temperatures, the bubble’s lifetime becomes almost constant. The overall trend is also observed in Sideman et al.’s experiments [51]. A simplified correlation is derived for the bubble’s lifetime based on energy conservation around the bubble in Equation (12):

$$t_b=\frac{\rho \left(c_p{T}_{sat}+h_{lg}\right)D_0}{2f_1\Delta T_{sub}}+f_2,$$

(12)

where $f_1$ and $f_2$ are fitting coefficients obtained via comparing with numerical simulations. $f_1$ represents the effective heat transfer coefficient $h_{eff}$ at the beginning of bubble condensation. Assume $f_2=f_3\sqrt[3]{D_0}$ (s), then the coefficients are $f_1=32,440$ W/m²K and $f_3=0.12$ s/m^1/3. This correlation reveals how bubble lifetime varies with bubble diameter, subcooled temperature, and latent heat. The details of the derivation can be found in Appendix A.

4. Conclusions

In the present paper, single vapor bubble condensation is modeled with the sharp interface Fourier mass transfer model. The model is proven to be highly accurate via comparison with the analytical solution for the case of the Stefan problem. It does not require any fitting parameters. However, the numerical model is sensitive to the initial thermal field, which is usually unknown for a particular configuration case. Therefore, in future work, we aspire to identify a correlation for estimating the thermal boundary region’s thickness around the bubble based on further experimental or analytical studies.

Additionally, a correlation for the bubble’s lifetime is derived based on the energy balance. The correlation shows good agreement with numerical simulation and relates the dependence of the condensing bubble’s lifetime to the initial diameter, subcooled temperature, and physical properties of the considered fluids.