Numerical simulation of annular flow boiling in microchannels

Abstract

Flow boiling in microchannels has been researched for decades due to its potential for heat transfer applications in micro structured devices. The large number of existing experimental datasets and fitted correlations show considerable scatter and lack a mechanistic basis. The computational fluid dynamics (CFD) model developed here applies to boiling in the annular flow regime and provides insights into the underlying behaviour. Analysis of the heat transfer augmentation by interfacial waves is studied. Parametric investigations are performed using the model to understand the impact of important flow and system parameters. The calculated heat transfer coefficients are consistent with those from existing correlations, showing no more scatter than the data used to create them.

Introduction

Development and application of micro structured devices in modern science and engineering technologies have brought major challenges for heat transfer in microchannels. ¹ Flow boiling has received increasing interest due to its potential for addressing the thermal management challenges arising from the increasing heat flux density in micro structured devices. ² However, the vast number of experimental and analytical studies of microchannel flow boiling conducted by many researchers in the past few decades show contradictory results. Computational fluid dynamics (CFD) approach has been used for microchannel flow boiling research in recent years, but little CFD work focuses on annular flow, which is the prevalent flow regime. This paper proposes a CFD model for the detailed study of annular flow boiling in microchannels. Verification of the CFD heat transfer model is presented and detailed analysis of the annular flow boiling mechanism is performed. Studies of the effects of some important flow and geometric parameters are also performed and discussed.

Literature review

Comprehensive reviews of microchannel flow boiling are available in Thome, ³ Garimella and Sobhan, ⁴ Bertsch et al., ⁵ Kandlikar ² and Baldassari and Marengo, ⁶ where many datasets from experimental studies are compared with the available correlations and theoretical models. However, not only do the experimental results show contradictory outcomes, but also the correlations do not provide satisfactory predictions. Experiments showed that the two-phase flow regime in microchannels transitions from bubbly to slug, slug-annular and annular flow at very low vapour qualities (x < 0.1). ⁷ Data from Lee and Mudawar ⁸ show that heat transfer is still strong in the annular regime (x > 0.1), where very few bubbles are observed and the heat transfer coefficient does not decrease much compared with that in the bubbly region. Particularly, data from Bertsch et al., ⁹ Lin et al. ¹⁰ and Yun et al. ¹¹ yield the peak heat transfer coefficients at vapour qualities in the range 0.1–0.6, in which the two-phase flow regime is annular.

Apart from the experimental and analytical studies, flow boiling has also been studied numerically using different approaches. Guo et al. ¹² reviewed the CFD work on microchannel flow boiling and presented the main numerical approaches used for the analyses. However, most of the flow boiling simulations focus on bubble growth,^13–15 bubble formation and the subsequent bubble dynamics,^16–21 bubble merging ²² and liquid evaporation in Taylor bubbles.^23,24 Only a few numerical studies of phase change in microchannel annular flow^25–27 have been reported and these are mostly for condensation.

Da Riva and Del Col ²⁵ performed a three-dimensional (3D), steady-state simulation of R134a condensation in a 75 mm long, 1 mm diameter circular microchannel using the commercial CFD solver Fluent 6.3. The grid spacing used in the core region was approximately 15 µm and was refined to 0.8 µm in the near wall region. An implicit volume of fluid (VOF) method ²⁸ and a modified HRIC scheme ²⁹ were used to model the multiphase interface. A modified k–ω turbulence model ³⁰ was applied in the gas phase. The phase change mass source (S) was computed based on the normalised temperature difference between the interface temperature and a constant saturation temperature as $S=-f_0{\alpha }_i{\rho }_i(T_{\text{sat}}-T)/T_{\text{sat}}$ . Here, ${\alpha }_i$ and ${\rho }_i$ are the volume fraction and the density of the ith phase, with $i=g,\text{or}l$ for condensation and evaporation processes, respectively, and f₀ = 750,000 s⁻¹ is an arbitrary coefficient chosen based on numerical testing. The vapour quality varied from 1 to 0.45 in this work so the flow is in the annular flow region. The simulation results showed the effect of gravity on the film thickness and the heat transfer coefficient. The effect of vapour quality on the heat transfer rate during the condensation process was also investigated. The computed heat transfer results were shown to be in good agreement with the experimental data of Matkovic et al. ³¹ The use of an artificial coefficient to calculate the phase change process does not reflect the real kinetics of interfacial mass transfer and may result in an incorrect temperature distribution in the liquid film.

Chen et al. ²⁶ used Fluent 6.1 to perform a three-dimensional transient condensation simulation of FC-72 in a 300 mm long microchannel with 1 mm² cross section. The computational domain contained around 200,000 cells, but the details of the grids in each spatial direction are not provided. However, it is estimated from the geometry and the cell number that the number of grids across the side length is less than 10. According to the analysis of Gupta et al. ³² which showed that at least five cells are needed to resolve a thin liquid film in microchannel slug flow, the mesh used in Chen et al. ²⁶ is clearly too coarse to perform accurate microchannel slug, transition and annular flow simulations. The VOF method was used to capture the interface and a second-order implicit method is used for temporal discretisation. A k–ɛ model ³³ was used to model turbulence. The interphase mass transfer was calculated using $S=-f_0{\alpha }_i{\rho }_i(T_{\text{sat}}-T)$ with f₀ = 100 s⁻¹K⁻¹ determined based on numerical testing. Although the researchers have shown the modelling results in slug formulation and transition to annular flow, the use of an arbitrary coefficient for phase change modelling is not a physical representation of the phase change process. Moreover, the k–ɛ model has been applied to the entire flow region in these simulations, whereas the formed liquid film and slugs during condensation processes are in the laminar region due to their low velocity.

A quasi-3D numerical approach for the simulation of annular flow evaporation in rectangular microchannels is proposed by Jesseela and Sobhan. ²⁷ The numerical approach used solves the equations on a cross section and then marches in the flow direction to obtain a 3D solution. In their work, annular flow phase change was modelled in a 15 mm long channel with a rectangular cross section. The edge lengths of the rectangular cross section vary from 30 µm to 300 µm. The mesh size in the cross-sectional domain was 41 × 81. It was assumed that the flow was laminar, incompressible and steady, for both the liquid and the vapour phases. Subcooled water enters at the inlet and then the phase change is calculated based on the heat input assuming that no superheating occurs in the gas and liquid phases. The numerical results were shown to be in agreement with the experimental data of Qu and Mudawar ³⁴ for the pressure drop along the flow direction. The results for the effect of vapour quality on heat transfer qualitatively match the observation by Ong and Thome ³⁵ for refrigerant flow. Although this approach provides some qualitative predictions of annular flow boiling, the assumption of steady annular flow is not universally valid. Furthermore, the assumption that all the input heat is used to evaporate liquid neglects the dynamics of the mass transfer at the phase interface.

The above review shows that despite the importance of annular flow boiling in microchannels, there is a complete lack of detailed numerical studies in this area. In our previous work, ³⁶ a numerical model for isothermal annular flow was proposed and validated. In this paper, a phase change model is added to form a complete model for annular flow boiling.

Mathematical modelling

The Eulerian-based multiphase flow model used previously ³⁶ is extended to include thermal and phase change models.

Governing equations

The governing equations (1) to (3) are based on the following assumptions:

All the fluids are Newtonian and incompressible.

Homogeneous velocities and temperatures are assumed, i.e. when the two fluids coexist they have the same velocity and temperature.

The gravitational/buoyancy force is neglected.

Viscous dissipation is neglected.

The flow is laminar, with velocities low enough to permit the use of a thermal energy, as opposed to a total energy, equation.

Continuity

$\frac{\partial \tilde{\rho }}{\partial t}+\nabla .(\tilde{\rho }v)=S_{\tilde{\rho }}$

(1)

Momentum

$\frac{\partial (\tilde{\rho }v)}{\partial t}+\nabla .(\tilde{\rho }v\otimes v)=-\nabla P+\nabla .(\tilde{\mu }(\nabla v+\nabla v^T))+F_{\sigma }+S_v$

(2)

Energy

$\frac{\partial \tilde{E}}{\partial t}+\nabla .(v\tilde{H})=\nabla .(\tilde{k}\nabla T)+S_{\tilde{H}}$

(3)

where $\tilde{\rho }$ , $\tilde{\mu }$ and $\tilde{k}$ are the bulk density, dynamic viscosity and thermal conductivity; v , P and T are the velocity, pressure and temperature; $\tilde{E}$ and $\tilde{H}$ are the bulk internal energy and enthalpy; t is time and $F_{\sigma }$ represents the interfacial force due to surface tension. $S_{\tilde{\rho }}$ is the phase change mass source, $S_v$ is the momentum transfer due to the mass source, and $S_{\tilde{H}}$ is the phase change energy source.

Additionally, a transport equation for a colour fraction is solved to identify the interface location and to calculate the volume fractions

$\frac{\partial {\rho }_{c}C}{\partial t}+v.\nabla {\rho }_{c}C=S_{{\rho }_c}$

(4)

where C, ${\rho }_c$ and $S_{{\rho }_c}$ are the colour fraction, the density and the phase change mass source for the phase used as the colour function. In the VOF model, the volume fraction (α) of a fluid phase is used as the colour function (C).

The volume fraction for one fluid is obtained using equation (4) and the other is obtained from ${\alpha }_j$  = 1 −  ${\alpha }_i$ . The bulk properties of the fluids in the governing equations are then calculated based on the phase volume fraction weighted values, as shown in equation (5) below

$\tilde{\phi }={\alpha }_l{\phi }_l+{\alpha }_g{\phi }_g$

(5)

where $\phi =\rho ,\mu ,k$ denote the density, viscosity and thermal conductivity, respectively. l and g denote the liquid and gas phases, respectively. The internal energy and enthalpy are defined using equation (6)

$\tilde{E}={\alpha }_l{\rho }_l{e}_l+{\alpha }_g{\rho }_g{e}_g$

(6a)

$\tilde{H}={\alpha }_l{\rho }_l{h}_l+{\alpha }_g{\rho }_g{h}_g$

(6b)

where e and h denote the specific internal energy and the specific enthalpy and are calculated based on the temperature.

The continuum surface force (CSF) model of Brackbill et al. ³⁷ is used for surface tension modelling. In this model, the surface tension force is calculated via

$F_{\sigma }=\sigma \kappa \delta (r-r_{\text{int}})n$

(7)

$n=\frac{\nabla C}{|\nabla C|};\kappa =-\nabla .n$

(8)

where $\delta (r-r_{\text{int}})$ is a numerical approximation of the Dirac delta function. A height function method is implemented in the conventional CSF model to alleviate the problem of spurious currents associated with simpler formulations, including the in-built model in ANSYS Fluent. ³⁶

The objective of this implementation is to replace the curvature calculation method of the original CSF method with a height function algorithm which yields more accurate curvature values than equation (8). The detailed development and implementation of the height function method into ANSYS Fluent is presented in Guo et al. ³⁸

Note that the bulk mass, momentum and energy sources in equations (1) to (3) are zero as conservation of mass, momentum and energy needs to be maintained during the phase change processes. However, local sources may become non-zero due to the numerical smoothing procedure used in their solution, described in the source terms used in the modelling of the phase change process section.

Phase change model

The phase change model implemented in the simulations is used to calculate the local interphase mass flux, depending on the local conditions at the liquid–gas interface. In the current simulations, the Tanasawa model³⁹ is modified for the flow boiling process, as shown in equation (9)

$j_m=\frac{2\gamma }{2-\gamma }(\frac{M}{2\pi R_{\text{u}}{T}_v}{)}^{1/2}(P_{\text{sat}}(T_{\text{int}})-P_{\infty })$

(9)

where j_m is the boiling mass flux, M is the molecular weight of the fluid, $R_{\text{u}}$ is the gas constant, and γ is the phase change or accommodation coefficient that represents the proportion of the liquid molecules absorbed by the gas phase of those leaving the liquid phase. T_v is the vapour temperature, $P_{\infty }$ is the system pressure and $P_{\text{sat}}$ is the saturation pressure evaluated at the temperature of the interface. The accommodation coefficient is set to 1 for water.^40,41 The saturation pressure is calculated using an Antoine equation

$P_{\text{sat}}=P_{\text{scale}}{10}^{(a-\frac{b}{T+c})}$

(10)

where the constants $P_{\text{scale}}$ , a, b and c are component-specific constants. For water over the temperature range 273–473 K, the values for $P_{\text{scale}}$ , a, b and c are 10⁵Pa, 5.12, 1687.54 K and −42.98 K, respectively. ⁴²

For boiling at an interface, the mass flux depends on the difference between the saturation pressure evaluated using the interface temperature and the local gas phase pressure. However, in any numerical simulation using a VOF method, the interface is smeared, so that the boiling source must be applied across a number of computational cells. It is straightforward to evaluate the saturation pressure in each cell but it is not so clear how to choose a local pressure, as the gas phase pressure is modified by surface tension effects inside the liquid region. There is only a small pressure drop along the flow direction (<1.3 kPa in the current simulations) which is small compared with the system pressure. Therefore, it was decided to use the constant system pressure $P_{\infty }$ to provide the evaporation driving force. There are clearly some numerical challenges with identifying the local gas phase pressure on a finite volume grid and applying it in such a model.

Boundary conditions

The inlet velocity profiles are obtained from the analytical solution of Fouilland et al. ⁴³

$v_{z_l}=-\frac{{\mathrm{BR}}^2}{4{\mu }_l}(1-{\xi }^2)$

(11a)

$v_{z_g}=-\frac{{\mathrm{BR}}^2}{4{\mu }_g}[({\xi }_i^2-{\xi }^2)+\frac{{\mu }_g}{{\mu }_l}(1-{\xi }_i^2)]$

(11b)

where ${\mu }_g$ and ${\mu }_l$ are the viscosity of gas and liquid, respectively, B is the axial pressure gradient, and R is the radius of the channel. $\xi =r/R$ is the non-dimensional radial coordinate with the subscript i denoting the location of the interface determined for a given set of gas and liquid superficial velocities using equation (12) below

$U_l=-\frac{{\mathrm{BR}}^2}{4{\mu }_l}[\frac{1}{2}+\frac{{\xi }_i^4}{2}-{\xi }_i^2]$

(12a)

$\mathrm{\hspace{1em}\hspace{1em}\hspace{1em}}{U}_g=-\frac{{\mathrm{BR}}^2}{4{\mu }_l}[(\frac{1}{2}\frac{{\mu }_l}{{\mu }_g}-1){\xi }_i^4+{\xi }_i^2]$

(12b)

In order to analyse how interfacial instabilities affect flow boiling, various interface perturbation approaches can be used. For the current work, the time-dependent interface location at the inlet is perturbed as in Guo et al. ³⁶ via

${\xi }_i(0,t)={\xi }_i+{\xi }_{\mathrm{iA}}\sin (2\pi f_{p}t)$

(13)

where f_p and ${\xi }_{\mathrm{iA}}$ are the frequency and dimensionless amplitude of the inlet wave, respectively.

Once the perturbed interface location is calculated using equation (13), the velocity profiles at the inlet boundary depend on time and cannot be specified using equations (12a) and (12b) which provides two constraints on the three remaining variables (U_l, U_g, and B). To achieve a closed solution, either the liquid superficial velocity (“fixed U_l perturbation”) or the pressure gradient (“fixed B perturbation”) was fixed. The inlet fluid velocity profiles were then determined according to equations (11) and/or (12). Note that the inlet mean flow rates of the smooth and two perturbed cases presented later are not identical but the differences between them are typically less than 1%.

A zero average pressure boundary condition is used for the outlet, together with a damping treatment to maintain the interface at the outlet parallel to the flow axis, in order to avoid reversed flow. The detailed implementation of the outlet boundary condition and the modified damping model of Park et al. ⁴⁴ are described in Guo et al. ³⁶

For the inlet thermal boundary condition, the gas and liquid temperatures are set to the saturation temperature for the system pressure under consideration. A constant heat flux is applied at the channel wall.

Initial conditions

The velocity and volume fraction fields are initialised using equations (11) and (12). The temperature and pressure fields are initialised with the constant system pressure and the corresponding saturation temperature. As the heat applied at the wall conducts into the liquid, it becomes superheated and the evaporation mass flux in equation (9) becomes non-zero. However, evaporation only occurs in those cells for which there is a non-zero interfacial area, with the evaporation source being the product of the flux and interfacial area.

CFD methodology

The commercial software ANSYS Fluent (14.5 and 15.0.7) is used to model 2D, axisymmetric gas–liquid annular flow in a microchannel. The double precision transient solver is used.

The volume of fluid method

As noted previously, the VOF method is used to capture the interface between the liquid film and gas core. The explicit Piecewise-Linear Interface Calculation (PLIC) scheme developed by Youngs ⁴⁵ is used to track the interface with sub-stepping to ensure the Courant number remains below 0.25. However, the actual time-step used for the simulations is smaller as it is limited by stability of the boiling source terms.

Source terms used in the modelling of the phase change process

To implement the phase change mass flux numerically, it must be transformed into a volumetric source. As noted above, this is done as follows:

$S_{{\rho }_g}=-S_{{\rho }_l}=j_m|\nabla {\alpha }_i|$

(14)

where $S_{{\rho }_g}$ and $S_{{\rho }_l}$ are source terms for the gas and liquid phases, respectively, and $|\nabla {\alpha }_i|$ is the simplest approximation to the interfacial area density. The integration of the interfacial area density over the computational domain gives the interfacial area between the two phases, as shown in equation (15)

${\int }_V|\nabla {\alpha }_i|\text{d}V=A_{\text{int}}$

(15)

where $A_{\text{int}}$ is the interfacial area between the two phases. ⁴⁶

As fluid is moved from the liquid to gas phase, it takes its momentum with it so that we have

$S_{v_l}=S_{{\rho }_l}v\mathrm{\hspace{1em}\hspace{1em}\hspace{1em}}{S}_{{\rho }_l}<0$

(16a)

$S_{v_g}=-S_{v_l}=S_{{\rho }_g}v$

(16b)

The energy sources also need to be calculated depending on the direction of phase change. For evaporation, the energy sources are

$S_{h_l}=S_{{\rho }_l}(h_l+h_{\mathrm{lv}})\mathrm{\hspace{1em}\hspace{1em}\hspace{1em}}{S}_{{\rho }_l}<0$

(17a)

$S_{h_g}=S_{{\rho }_g}(h_l+h_{\mathrm{lv}})$

(17b)

where h_l and h_g are the liquid and gas specific enthalpy, respectively, and $h_{\mathrm{lv}}$ is the latent heat of vaporisation.

A highly superheated interface will generate a high mass transfer rate which may cause numerical instabilities to occur. This is particularly so if a high evaporation flux is calculated in a cell that contains relatively little liquid. A numerical smoothing procedure proposed by Hardt and Wondra ⁴⁶ used in their phase change simulations was implemented to improve the numerical stability by means of shifting such a source to the liquid side of the interface and smearing the mass source terms ( $S_{{\rho }_i}$ ) over a few cells near the interface.

First, the boiling mass source term is shifted to the liquid side and the latent heat source is then calculated using the shifted mass source. This step is shown in equation (18) below

$S_{{\rho }_g}\text{'}=N_{{\rho }_g}{\alpha }_l{S}_{{\rho }_g}$

(18a)

$N_{{\rho }_g}=\frac{{\int }_V{S}_{{\rho }_g}\text{d}V}{{\int }_V{\alpha }_l{S}_{{\rho }_g}\text{d}V}$

(18b)

$S_{\mathrm{lv}}\text{'}=-S_{{\rho }_g}\text{'}{h}_{\mathrm{lv}}$

(18c)

where $S_{{\rho }_g}\text{'}$ is the shifted gas mass source and $N_{{\rho }_g}$ is a normalisation factor. Note that the choice of gas or liquid sources to be shifted does not affect the absolute value of the shifted source. Here the gas mass source, which is positive for boiling, is shown. The latent heat source $S_{\mathrm{lv}}\text{'}$ calculated using this shifting procedure will only be non-zero in cells containing liquid and the use of a normalisation factor conserves the mass source and therefore the latent heat source in the computational domain.

The second step of the procedure is smoothing of the mass source terms to avoid the generation of unphysical local velocities due to large source term, and the mass source term is smeared over a few cells near the interface, using the diffusion approach shown in equation (19) to generate the diffused mass source ${\varphi }_g$

$\nabla ·(\nabla {\varphi }_g)=\frac{1}{Z}({\varphi }_g-S_{{\rho }_g}\text{'})\mathrm{\hspace{1em}\hspace{1em}\hspace{1em}}\text{with}\mathrm{\hspace{1em}\hspace{1em}\hspace{1em}}\nabla {\varphi }_g=0\text{atallboundaries}$

(19)

where Z is a diffusion parameter (equivalent to the parameter combination $\Delta {\mathrm{tD}}_i$ in a discretised diffusion equation). Here, the value of Z is adjusted to give a diffusion length of ∼0.02D. This diffusion length is short and therefore does not have a significant effect on the calculated solution. However, this smoothing can compromise the shifting performed in step one, so the shifting procedure described earlier is repeated on the diffused sources to prevent this.

An example of the original and the modified source terms is shown in Figure 1. The legends in the figure show that smoothed gas and liquid mass sources ( $S_{{\varphi }_i}$ ) have smaller magnitudes than the original mass source ( $S_{{\rho }_g}$ ).

Figure 1.

Typical original and modified phase change source terms in a boiling simulation: (a) original gas mass source; (b) shifted latent heat source; (c) diffused gas source; (d) final gas mass source; and (e) final liquid mass source. The solid black lines denote the gas–liquid interface.

After the mass source terms are smoothed, the corresponding bulk mass, momentum and energy sources need to be calculated based on the smoothed mass source ( $S_{{\varphi }_i}$ ). The bulk mass and momentum source become

$S_{\tilde{\rho }}=S_{{\varphi }_g}+S_{{\varphi }_l}$

(20a)

$S_v=S_{{\varphi }_g}v+S_{{\varphi }_l}v$

(20b)

As the latent heat sources are calculated individually, the energy source now consists of the latent heat source and the enthalpies of the removed liquid and added gas

$S_{\tilde{H}}=S_{{\varphi }_g}{h}_g+S_{{\varphi }_l}{h}_l-S_{{\rho }_g}\text{'}{h}_{\mathrm{lv}}$

(20c)

As shown in Figure 1, the final mass sources of the gas and liquid phases are diffused over several cells into each phase separately, so the bulk sources calculated using equation (20) may become locally non-zero. Therefore, the use of equations (20b) and (20c) is critical to prevent artificial acceleration/deceleration and heating/cooling. Global conservation is guaranteed as the smoothing procedure is performed under this constraint. However, as shown in Figures 1(d) and (e), the final gas and liquid sources are located in different cells which may have different local velocity and temperature values. Therefore, the summed momentum and energy sources calculated as volume integrals of the bulk source terms may become non-zero. However, the mass source is controlled to be smeared over only a few cells near the interface to minimise this effect in order to conserve momentum and energy during the phase change processes.

Coding and parallel considerations

The boundary and initial conditions, the outlet damping model, the height function method and the phase change model including the numerical smoothing procedure are coded in User-Defined Functions (UDFs) of ANSYS Fluent. The simulations are run in parallel on multiple computational nodes and parallelisation was ensured as described in Guo et al. ³⁸ and the ANSYS Fluent 15.0 UDF Manual. ⁴⁷

Numerical solver options

The SIMPLE algorithm of Patankar ⁴⁸ is used for the pressure–velocity coupling. The implicit body force option is enabled to take into account the partial equilibrium of the pressure gradient and the body forces because the surface tension force is large compared with the viscous force ( ${\mathrm{Ca}}_g$  < 0.005, ${\mathrm{Ca}}_l$  < 0.0008) and is significant compared with the inertial force ( ${\mathrm{We}}_g$  < 3.2, ${\mathrm{We}}_l$  < 0.25) in the current work. The Green-Gauss node-based method is used to calculate gradients. The PRESTO! scheme is used for pressure interpolation. The QUICK scheme is used for the discretisation of the momentum and energy equations. A first-order implicit scheme is used for the temporal discretisation of the transient terms because the use of second-order implicit scheme is not compatible with the explicit PLIC scheme for the volume fraction equation. Globally scaled residuals are used and the normalised residual targets for all equations are set to 1 × 10⁻⁶ except for the continuity equation in some high heat flux cases where 1 × 10⁻⁵ is used to accelerate convergence.

Simulation results

The annular flow simulations are performed using circular channels of diameter D = 0.25 to 2 mm and length L = 10D to 50D. The 2D, axisymmetric computational domain is shown in Figure 2. The initial section of length L_A is the region of interest where the phase change model is enabled. The final region of length L_D serves as a numerical beach for wave damping to minimise any effect of the exit boundary conditions on the boiling zone. The two are separated by a short (∼D) adiabatic zone. The physical properties of water and steam are extracted from the NIST chemistry webbook. ⁴⁹

Figure 2.

Schematic representation of the computational domain.

Grid-independence study

In our previous work, the dynamics of the interfacial waves in microchannel annular flow have been shown to be independent of grid spacing at Δ = D/200. ³⁶ To perform the simulations with heat transfer and phase change, the mesh density needs to be checked with the energy equation and the phase change model enabled. An annular flow boiling simulation in a 1 mm diameter channel of L = 13D (L_A = 10D and L_D = 2D) at 101.3 kPa is used for this purpose. Liquid and gas superficial velocities of U_l = 0.11 ms⁻¹ and U_g =  6 ms⁻¹ are used to calculate the velocity profile for the inlet boundary and the initialisation. A fixed U_l perturbation of amplitude ${\xi }_{\mathrm{iA}}$  = 0.002 and frequency f_p = 100 Hz is also applied at the inlet gas–liquid interface. A heat flux of 25 kW m⁻² is applied at the channel wall. A diffusion parameter Z  = 4 × 10⁻¹⁰ m² is used for the numerical smoothing procedure and a time-step of 5 × 10⁻⁷ s is used.

The boiling simulation was firstly performed using four square meshes with grid spacing Δ = D/100 to Δ = D/300. To characterise the boiling process, the overall boiling mass source over the computational domain is calculated as

$S_{{\rho }_g\text{overall}}={\int }_V{S}_{{\rho }_g}\text{d}V$

(21)

Figure 3 illustrates the overall vaporisation mass source as a function of the simulation time. It is clear from Figure 3 that when the mesh is refined, the calculated evaporation rate decreases until the mesh is refined to Δ = D/200 beyond which further refinement has no impact. Moreover, the simulations for the two low mesh resolution cases diverged shortly after t = 20 ms.

Figure 3.

Overall vapour mass source obtained using various mesh densities.

As shown in Figure 3, a rectangular mesh of aspect ratio $\Delta z/\Delta r$  = 2, with $\Delta r$  = D/200, yields the same results as the more refined square meshes. Further increases in the aspect ratio to $\Delta z/\Delta r$  =  4 and 8 led to the generation of unphysical velocities in the damping region that eventually destroyed the simulations. Therefore, the mesh resolution used in the present work was chosen as $\Delta z/\Delta r$  = 2 with $\Delta r$  = D/200.

In the above simulations, a time-step of 5 × 10⁻⁷ s corresponds to a maximum Courant number of 0.10 for the grid spacing Δ = D/200. The time-step size for the Δ = D/200 case was increased from 5 × 10 ^− 7 s to 1 × 10⁻⁶ s and the results showed no difference. The time-step size was not increased beyond 1 × 10⁻⁶ s as the Courant number reached 0.20 at $\Delta t$  = 1 × 10⁻⁶ s. This test shows that time-steps corresponding to maximum Courant numbers smaller than 0.25 give identical results. The time-step size used in the current work varies with different cases and was determined for each case based upon extensive numerical testing.

Model testing

Steady, laminar, fully developed annular flow heat transfer without phase change

The model is tested initially for sensible heat transfer to an annular flow water/steam flow at 101.3 kPa in 1 mm diameter channel (L = 13D, L_A = 10D and L_D = 2D). The inlet liquid and gas superficial velocities are U_l = 0.11 ms⁻¹ and U_g = 6 ms⁻¹, with the corresponding mass flows and quality remaining unchanged through the section in the absence of evaporation. A wall heat flux of 25 kW m⁻² is applied. A time-step of 1 × 10⁻⁶ s is used and the simulation is performed for a physical time of 100 ms. The axial wall temperature gradient becomes constant after t = 50 ms, indicating thermally fully developed flow.

Figure 4 shows the axial distribution of the local heat transfer coefficient calculated using equation (22)

$h_{\text{tp},z}=\frac{q}{T_{\text{w},z}-T_{\text{b},z}}$

(22)

where $h_{\text{tp},z}$ is the local two-phase heat transfer coefficient and $T_{\text{w},z}$ is the local wall temperature. $T_{\text{b},z}$ is the local bulk temperature at an axial location and is calculated as

$T_{\text{b},z}=\frac{{\int }_0^R|v_z|\widetilde{\rho C_p}{T}_{z,r}2\pi r\text{d}r}{{\int }_0^R|v_z|\widetilde{\rho C_p}2\pi r\text{d}r}$

(23)

with $T_{z,r}$ being the local temperature. $\widetilde{\rho C_p}$ is the bulk volumetric heat capacity calculated as given by equation (5). In keeping with the expected variation of $h_{\text{tp},z}$ in the direction of flow in a laminar, thermally developing flow, the results in Figure 4 are presented as a function of the scaled axial length, or inverse Graetz number calculated as ${\mathrm{Gz}}^{-1}=z/{\mathrm{DRe}}_l{\mathrm{Pr}}_l$ . As shown in Figure 4, the heat transfer coefficient is initially high and decreases to a constant value, and is in precise agreement with the analytical solution for fully developed annular flow (neglecting axial conduction) provided by Gupta. ⁵⁰

Figure 4.

Validation of the CFD heat transfer model of steady, laminar, fully developed annular flow using the analytical model of Gupta. ⁵⁰ The CFD heat transfer coefficient profile is obtained at t = 100 ms.

Assessment of the impact of interface perturbations

Boiling simulations are performed using the dimensions and inlet conditions for the case described in the ‘Steady, laminar, fully developed annular flow heat transfer without phase change' section with the phase change model enabled. In addition, inlet perturbations (amplitude ${\xi }_{\mathrm{iA}}$  = 0.002 and frequency f_p = 100 Hz) are examined, with both the fixed liquid flow U_l and fixed pressure gradient B perturbation methods described earlier employed. A diffusion parameter Z  = 4 × 10⁻¹⁰ m² is used.

Figure 5(a) shows the evolution of the total vapour mass source with time for the steady and perturbed inlets. After t = 40 ms, the vapour source for the constant inlet case reaches a plateau with only small fluctuations. With the imposition of inlet fluctuations, the overall vapour mass sources become periodic after t = ∼50 ms, indicating that the flows are quasi-steady. The oscillation frequencies of the source for both perturbation cases are similar, but the magnitudes of the source of the fixed liquid flow rate case are slightly higher than those of the fixed pressure gradient case. Fast Fourier Transform (FFT) analysis shows that the source frequencies for both perturbation cases are ∼10 ms which is equivalent to the perturbation frequency. The average magnitudes of the sources in the quasi-steady phase of the three cases are not exactly the same because the mean mass flow in the three cases is not identical.

Figure 5.

Exploration of the annular flow boiling model for water at 101.3 kPa with U_l = 0.11 ms⁻¹ and U_g = 6 ms⁻¹; D = 1 mm and wall heat flux of 25 kW m⁻². (a) The evolution of vapour mass sources with time for different inlet perturbation conditions. The axial wall temperature at t = 100 ms profiles for all four cases is shown in (b).

Figure 5(b) shows the axial wall temperature distribution for all three boiling cases at t = 100 ms, with the non-boiling results also shown for comparison. At low values of ${\mathrm{Gz}}^{-1}$ , all four cases show the same profile, suggesting that the establishment of the film temperature profile by conduction is the same in all cases. However, towards the end of the conductive thermal development region, the evaporative cooling associated with the boiling cases sees their rates of rise of the wall temperature drop below that of the non-boiling case. The residual slow rise in the wall temperature at the end of the computational domain reflects continuing thermal development. We note that the wall temperatures for the three boiling cases are essentially indistinguishable in Figure 5(b) and that the imposed inlet perturbations have not had a significant effect on heat transfer through the annular liquid layer. In the following boiling simulations, only the flat interface initialisation and non-perturbed inlet are used.

Boiling heat transfer

In annular flow boiling, the vapour quality increases along the flow direction due to the continuous heat input. For a channel with diameter D and length L_A, assuming all heat input is used to vaporise the liquid, the boiling mass transfer rate ( ${\stackrel{·}{m}}_{\text{bo},q}$ ) is

${\stackrel{·}{m}}_{\text{bo},q}=\frac{q\pi {\mathrm{DL}}_A}{h_{\mathrm{lv}}}$

(24)

The change in vapour quality ( $\Delta x$ ) from the inlet to the outlet is therefore

$\Delta x=\frac{{\stackrel{·}{m}}_{\text{bo},q}}{{\stackrel{·}{m}}_{\text{in}}}=\frac{q\pi {\mathrm{DL}}_A/h_{\mathrm{lv}}}{G\pi (\frac{D}{2}{)}^2}=\mathrm{Bo}·\frac{4L_A}{D}$

(25)

where Bo is the boiling number. The calculated $\Delta x$ in the current work ranges from 0.003 to 0.024. As the vapour quality used in the initialisation of the current simulations ranges from 0.1 to 0.3, this change may have an effect on the heat transfer rate.

In the following sections, we consider boiling of water at 160 kPa ( $T_{\text{sat}}$  = 386.45 K) in a 0.5 mm diameter channel (L = 20D, L_A = 15D and L_D = 4D) with constant fluid mass flux G = 60 kgm⁻²s⁻¹ and inlet vapour quality x = 0.1. The inlet superficial velocities calculated based on the above parameters are U_l = 0.057 ms ^− 1 and U_g = 6.55 ms⁻¹. A constant heat flux of 20 kW m⁻² is applied at the channel wall, giving a maximum vapour quality change $\Delta x$  ∼ 0.009 which is ∼9% of the inlet vapour quality. The time-step for these cases is fixed at 5 × 10⁻⁷ s.

Transient behaviour of annular flow boiling in microchannels

Figure 6 shows the evolution of the global mass source calculated using equation (21). Beyond the initial transient discussed previously, the global source reaches a peak value at t = ∼12 ms after which it reaches an approximate (±10%) plateau. The dashed line, denoting the theoretical boiling mass source ${\stackrel{·}{m}}_{\text{bo},q\text{w}}$ calculated using equation (24), is slightly higher than this plateau value because, in the simulation, some heat input has gone into the liquid phase to create a superheated liquid film, as evident in Figure 7. The interfacial waves seen in Figure 7 are a result of the phase change process at the interface giving rise to the periodic fluctuations in the global vapour source seen in Figure 6. The origin of these waves is discussed later.

Figure 6.

Evolution of the amount of vapour mass source in the computational domain (D = 0.5 mm, U_g = 6.55 ms⁻¹, U_l = 0.057 ms⁻¹, q = 20 kWm⁻², x = 0.1). The dashed line denotes the calculated boiling source ${\stackrel{·}{m}}_{\text{bo},q}$ assuming all heat input is used to vaporise the liquid.

Figure 7.

Characterisation of local boiling heat transfer at t = 30 ms for a microchannel boiling annular flow (D = 0.5 mm, U_g = 6.55 ms⁻¹, U_l = 0.057 ms⁻¹, q = 20 kWm⁻², x = 0.1). (a) The temperature field and the local boiling heat transfer coefficient distribution along the axial direction. (b) Local heat transfer coefficient distribution and the temperature and velocity fields near the wave crest in the dashed rectangle in (a). The black solid lines in the temperature contours denote the gas–liquid interface. The axial length of the channel was scaled by a factor of 0.25 and 0.5 for (a) and (b). The display density of the vector field in (b) was reduced for better presentation.

Figure 7(a) also shows the distribution of local boiling heat transfer coefficient calculated using equation (26)

$h_{\text{bo},z}=\frac{q}{T_{\text{w},z}-T_{\text{sat}}}$

(26)

where $T_{\text{w},z}$ is the local wall temperature. The peaks in $h_{\text{bo},z}$ are closely associated with the location of the interfacial wave troughs because the wall temperature is lowered by enhanced evaporative cooling due to better heat conduction through a thinner film. This effect is seen in more detail in Figure 7(b) which is an enlargement of the dotted region around $z/D$  = 9. However, Figure 7(b) also reveals that recirculation in the troughs provides an additional convective contribution to the film heat transfer in this region.

Determination of developed flow and average boiling heat transfer coefficient

As the interfacial waves shown in Figure 7(a) are advected by the flow, the local heat transfer coefficient becomes time-dependent due to the oscillating interface radius at the same axial location. To characterise the time-dependent interfacial radius and the corresponding local heat transfer, monitor point pairs are placed at $z/D$  = 3.35, 6.70, 10.05, 13.40 (corresponding to ${\mathrm{Gz}}^{-1}$  = 0.02, 0.04, 0.06 and 0.08, respectively). The results for ${\mathrm{Gz}}^{-1}$  = 0.06 are shown in Figure 8 where it is seen that, after t = ∼ 15 ms, both the local wall temperature and the interface radius show regular periodic fluctuations. Qualitatively, the same situation arises at the other monitor points.

Figure 8.

The monitored wall temperature and interface radius at the same axial location ${\mathrm{Gz}}^{-1}$  = 0.06.

Fast Fourier transform (FFT) analyses of the respective periodic behaviours are shown in Figures 9(a) and 9(b) for interface location and wall temperature, respectively. The same dominant frequency of ∼250 Hz is seen in every case, indicating the generation and advection of regular waves at the gas–liquid interface. The amplitudes of the component frequencies are shown in Figures 9(c) and 9(d). These are smallest at the first axial location ( $T_{\text{w}1}$ and $r_{\text{i}1}$ ), in keeping with the observation in Figure 7 that the film thickness remains at its inlet value some distance into the analysis section. The amplitudes for the remaining points are similar, with those for the third and fourth points being essentially indistinguishable, suggesting that the flow, while unsteady, is quasi-steady after the third axial location.

Figure 9.

FFT analysis for the interface radius and the wall temperature variations at four axial locations. The power spectral densities (PSD) for the interface radius and the wall temperature are shown in (a) and (b), respectively. The amplitudes of the component frequencies are shown in (c) and (d).

The average boiling local heat transfer coefficient for a monitored axial location can be calculated using equation (27)

${\overline{h}}_{\text{bo},z}=\frac{q}{{\overline{T}}_{\text{w},z}-T_{\text{sat}}}$

(27)

where ${\overline{T}}_{\text{w},z}$ is the time-averaged local wall temperature calculated as shown in equation (28)

${\overline{T}}_{\text{w},z}=\frac{1}{(t_{\text{e}}-t_{\text{s}})}{\int }_{t_{\text{s}}}^{t_{\text{e}}}{T}_{\text{w},z}\text{d}t$

(28)

where $t_{\text{s}}$ and $t_{\text{e}}$ denote the start and end times separated by an integral number of wave periods, as shown in Figure 8. The calculated value of ${\overline{h}}_{\text{bo},z}$ at the four monitor point locations, shown in Figure 10, clearly shows the approach to quasi-steady behaviour. Here it must be recalled, however, that the vapour quality increases in the flow direction in boiling flows and that the flow and heat transfer behaviour must change along the duct. However, in the present case, there is very little change in vapour quality between the third and fourth monitor points (x = 0.106 and 0.108, respectively) and it is therefore not surprising that ${\overline{h}}_{\text{bo},z}$ appears to approach a quasi-steady value.

Figure 10.

Average heat transfer coefficients at different monitored axial locations.

Heat transfer enhancement caused by the interfacial waves

For steady-state conduction across a thin laminar film of thickness ${\delta }_{\text{F},z}$ , the heat transfer coefficient can be estimated as

$h_{\delta ,z}=\frac{k_l}{{\delta }_{\text{F},z}}$

(29)

and it is instructive to compare this simple result with the detailed simulations since, as described in the ‘Transient behaviour of annular flow boiling in microchannels’ section, the boiling heat transfer coefficient in the simulations is enhanced by the transient effects.

Figure 11 shows this comparison for the axial location of ${\mathrm{Gz}}^{-1}$  = 0.06. At simulation time t₁, the film is at its thinnest and the film-conduction heat transfer coefficient takes its highest value. However, the simulated boiling heat transfer coefficient does not reach its peak until shortly after t₁, due to the time delay required for the effect of the passage of the interfacial wave to be observed on the wall temperature. From t₁ to t₂, the film conduction heat transfer coefficient falls sharply as the wave trough passes through but this effect is muted in the simulations, presumably because of flow recirculation – such recirculation provides additional radial convection immediately after the passage of the wave trough, as shown in Figure 7(b). After t = t₂, the film conduction heat transfer increases as the next wave trough approaches, but the simulated boiling heat transfer coefficient continues to decrease until it reaches its lowest value at t = t₃ where it is lower than the film conduction value. From t₃ to t₄, both boiling and film conduction heat transfer coefficients increase by similar amounts due to the approach of the next wave trough which arrives at t₄, thus completing a cycle.

Figure 11.

Comparison of local heat transfer coefficient between annular flow boiling and film conduction model at the axial location of ${\mathrm{Gz}}^{-1}$  = 0.06. The interface coordinate is also plotted to track the passage of the wave.

The average local film thickness calculated using an approach similar to the average wall temperature (equation (28)) yields ${\overline{\delta }}_{\text{F},z}$  = 30.1 µm and $h_{\delta }$  = 22.7 kWm⁻² K⁻¹. This film thickness is 25% thinner than the inlet (non-boiling smooth) film thickness (40.2 µm). Compared with the simulation average $h_{\text{bo}}$  = 32.4 kWm⁻² K⁻¹, it is apparent that the flow patterns associated with the interfacial waves have enhanced the heat transfer process by a factor of 1.43.

Origin of the waves

The mechanism by which waves form in the absence of any introduced perturbation requires investigation. The question is whether this is a numerical artefact or something that represents a physical phenomenon. In order to investigate this, a number of test cases were performed as follows:

A simulation was performed in which the boiling source was calculated based on the boiling number and applied at the cell in which the interfacial area was highest at each radial location. Source term smoothing and shifting were disabled. No waves were generated.

The above case was modified to apply the source at all cells with a non-zero interfacial area and the source was weighted by the cell value of the interfacial area. This saw the source being spread over several cells but no waves were generated.

Visualisation of the source terms in a standard boiling simulation showed that prior to shifting and smoothing locally, high evaporation sources could occur in cells where the interfacial area was non-zero in a cell closer to the wall than its neighbours. Examination of the components of the source showed that the evaporation flux j_m increases rapidly towards the wall so that any small interfacial area term in a cell nearer to the wall can cause a large source term locally. Even though this effect is reduced in the smoothing step, there are still locally high sources. This potentially explains a mechanism in which locally high source terms can generate instability at the interface.

Quantitative examination of the source components of five adjacent interfacial cells at the liquid side shows an average temperature jump of ∼0.015 K between these neighbouring cells, corresponding to a mass flux difference of ∼0.1 kgm⁻²s⁻¹. Compared with the average evaporation mass flux of ∼0.012 kg m⁻² s⁻¹ over the interface, this difference shows the high sensitivity of the evaporation mass flux to the interface temperature variation and could cause local perturbations to the interface by changing the volume fraction values of the interfacial cells non-uniformly.

A standard boiling simulation was run in which the momentum source terms were suppressed. These can potentially play a role in generating instabilities as explained in Granovskii et al. ⁵¹ However, even in this case, wave formation was still observed.

The interfacial area is represented using the gradient of the volume fraction and it is known that this will diffuse the interface across several numerical cells. It is possible to use other functions that provide a much sharper interface. To test the impact of this, the interfacial area density was calculated using

${\int }_{V}6{\alpha }_i(1-{\alpha }_i)|\nabla {\alpha }_i|\text{d}V=A_{\text{int}}$

(30)

which gave an interface that was only one cell thick. When coupled with the boiling source term (with smoothing and shifting), this also gave rise to waves at the interface but the boiling source was observed to oscillate wildly with time.

A calculation of local mass flux using the evaporation model used by Hardt and Wondra ⁴⁶ yields identical results to that obtained using the current phase change model, confirming that no additional instabilities are introduced by using the pressure difference form of the Tanasawa ³⁹ model instead of using the temperature difference form.

A boiling case performed with the phase change model disabled in the near-inlet region (0 <  $z/D$  < 4) give similar waves to the standard boiling case, showing that the near-inlet interface region plays no role in the initial perturbation of the interface.

The boiling simulation was run with the phase change kinetic coefficient ( $(2\gamma /2-\gamma )(M/(2\pi R_{\text{u}}{T}_v){)}^{1/2}$ ) scaled by a factor of 0.1. The interface temperature values became higher than those of the standard boiling case, keeping an almost identical interphase mass transfer rate and still giving waves very similar to those of the standard boiling case.

In summary, these tests have shown that neither the phase change model nor the numerical smoothing procedure have introduced numerical artefacts that generate waves. However, the observed high sensitivity of evaporation mass flux to the interface temperature has shown large potential for the generation of local perturbation when the interfacial area moves into a cell near the wall. If this happens, the rapid local increase of the evaporation mass flux causes a very high local evaporation flux and thus acts as a perturbation to the interface. In any real system, there are a variety of sources of possible perturbations, such as non-uniformity of the wall heat flux, wall surface finish variations, 3D perturbations to the flow, etc. that could provide the slight asymmetries needed to induce wave formation. We therefore proceed on the basis that the waves generated are a consequence of the physical modelling assumption.

Parametric analysis of annular flow boiling heat transfer

In this section, simulations for various conditions are performed to evaluate the effects of the mass flux, channel diameter, system pressure, vapour quality and wall heat flux on the boiling heat transfer coefficient. The conditions for the different cases are listed in Table 1, where Case 2 is the base case used in the ‘Boiling heat transfer’ section. The length of channel is chosen to have an analysis section long enough to obtain quasi-steady flow ( ${\mathrm{Gz}}^{-1}$  > 0.08). The number of cases studied is limited by the computational cost, with average wall clock times of around 160 h when using 8 to 16 cores (depending on the length of the tube) on a modern Linux cluster.

Table 1.

List of cases used for the parametric analysis.

Case no.	G(kgm−2s−1)	D (mm)	$P_{\infty }$ (kPa)	x	q (kWm−2)	L (D)	LA (D)	Bo	Z (m2)	$\Delta t$ (s)
1	21	0.5	160	0.1	20	10	5	4.2 × 10−4	1 × 10−10	1 × 10−6
2	60	0.5	160	0.1	20	20	15	1.5 × 10−4	1 × 10−10	5 × 10−7
3	120	0.5	160	0.1	20	32	26	7.5 × 10−5	1 × 10−10	5 × 10−7
4	179	0.5	160	0.1	20	50	40	5.0 × 10−5	1 × 10−10	3 × 10−7
5	60	0.25	160	0.1	20	10	7	1.5 × 10−4	2.5 × 10−11	2 × 10−7
6	60	0.9	160	0.1	20	20	15	1.5 × 10−4	3.2 × 10−10	5 × 10−7
7	60	1.5	160	0.1	20	26	21	1.5 × 10−4	9 × 10−10	1 × 10−6
8	60	2	160	0.1	20	32	26	1.5 × 10−4	1.6 × 10−9	1 × 10−6
9	60	0.5	210	0.1	20	20	15	1.5 × 10−4	1 × 10−10	5 × 10−7
10	60	0.5	310	0.1	20	20	15	1.5 × 10−4	1 × 10−10	1 × 10−6
11	60	0.5	160	0.15	20	20	15	1.5 × 10−4	1 × 10−10	5 × 10−7
12	60	0.5	160	0.2	20	18	13	1.5 × 10−4	1 × 10−10	5 × 10−7
13	60	0.5	160	0.25	20	18	13	1.5 × 10−4	1 × 10−10	5 × 10−7
14	60	0.5	160	0.3	20	18	13	1.5 × 10−4	1 × 10−10	5 × 10−7
15	60	0.5	160	0.1	5	20	15	3.8 × 10−5	1 × 10−10	5 × 10−7
16	60	0.5	160	0.1	10	20	15	7.5 × 10−5	1 × 10−10	5 × 10−7
17	60	0.5	160	0.1	15	20	15	1.1 × 10−4	1 × 10−10	5 × 10−7
18	60	0.5	160	0.1	30	20	15	2.3 × 10−4	1 × 10−10	5 × 10−7
19	60	0.5	160	0.1	40	20	15	3.0 × 10−4	1 × 10−10	5 × 10−7

Effect of mass flux

Figure 12 shows the effects of varying the mass flux from 60 to 179 kgm⁻²s⁻¹. The data from Case 1 are not plotted because the low superficial velocity led to a slug-annular flow in the simulation. Overall, it is found that both boiling and film conduction heat transfer coefficients increase slightly with increasing mass flux. This effect is correlated with higher interfacial wave amplitudes relative to the average film thickness, In addition, the frequency of the waves increases with increasing mass flux, as shown in Figure 12(b). The enhancement factor ( $h_{\text{bo}}/h_{\delta }$ ) remains constant at a value of ∼1.43 over the mass flux range considered.

Figure 12.

Effect of fluid mass flux (cases 2–4) on annular flow boiling (D = 0.5 mm, $P_{\infty }$  = 160 kPa, x = 0.1, q = 20 kWm⁻²): (a) boiling and film conduction heat transfer coefficients; and (b) main frequency and the corresponding amplitude of the interfacial waves.

Effect of tube diameter

Figure 13(a) shows that the boiling heat transfer coefficient decreases significantly as the channel diameter is increased from 0.25 mm to 2 mm. The increase in film thickness with increasing tube diameter is the main reason for this effect but the smaller relative amplitude of the waves, Figure 13(b), also contributes. Figure 13(b) shows that the frequency of the interfacial waves varies inversely with channel diameter, suggesting an approximately constant value of Strouhal number, $\mathrm{St}=\mathrm{fd}/U_{\mathrm{TP}}.$ Indeed, for all of the runs we have carried out, 0.011 $<\mathrm{St}<$ 0.025, but we have not investigated this further here. The enhancement factor lies between 1.36 and 1.43 and decreases slightly with increasing diameter.

Figure 13.

Effect of channel diameter (cases 2 and 5–8) on annular flow boiling (G = 60 kgm⁻²s⁻¹, $P_{\infty }$  = 160 kPa, x = 0.1, q = 20 kWm⁻²): (a) boiling and film conduction heat transfer coefficients; and (b) main frequency and the corresponding amplitude of the interfacial waves.

Effect of system pressure

Figure 14(a) shows that the simulated boiling heat transfer coefficient increases with pressure in the range 160–310 kPa, even as the film thickens and the film conduction heat transfer coefficient decreases slightly. These opposing trends arise because the interfacial wave amplitude increases as the pressure increases, as shown in Figure 14(b). At the same time, the inlet superficial gas velocity decreases from 6.55 ms⁻¹ to 3.52 ms⁻¹ and the wave frequency shown in Figure 14(b) also decreases. The enhancement factor between the boiling and film conduction heat transfer coefficients increases from 1.43 to 1.91 with increasing pressure.

Figure 14.

Effect of system pressure (cases 2 and 9–10) on annular flow boiling (G = 60 kgm⁻²s⁻¹, D = 0.5 mm, x = 0.1, q = 20 kWm⁻²): (a) boiling and film conduction heat transfer coefficients; and (b) main frequency and the corresponding amplitude of the interfacial waves.

Effect of vapour quality

In the simulations, the mean liquid film thickness decreases significantly with increasing vapour quality due to higher gas superficial velocity and lower liquid superficial velocity. Therefore, as shown in Figure 15(a), the film conduction heat transfer coefficient increases with increasing vapour quality, as does the simulated boiling heat transfer coefficient. The non-dimensional wave amplitude is shown to decrease with increasing vapour quality in Figure 15(b), indicating that the film is relatively smoother at higher quality. Figure 15(b) also shows that higher wave frequencies are predicted with higher vapour quality. The enhancement factor ranges from 1.43 to 1.59.

Figure 15.

Effect of vapour quality (cases 2 and 11–14) on annular flow boiling (G = 60 kgm⁻²s⁻¹, D = 0.5 mm, $P_{\infty }$  = 160 kPa, q = 20 kWm⁻²): (a) boiling and film conduction heat transfer coefficients; and (b) main frequency and the corresponding amplitude of the interfacial waves.

Effect of wall heat flux

The results from the variation of the wall heat flux are shown in Figure 16. Both boiling and film conduction heat transfer coefficients increase with increasing heat flux as the non-dimensional interfacial wave amplitude increases slightly with increasing heat flux (Figure 16(b)). The frequencies shown in Figure 16(b) have little dependence on the heat flux – since the gas velocity is also only slightly affected by the magnitude of the heat flux in these simulations. The enhancement factor stays in a narrow region between 1.40 and 1.44.

Figure 16.

Effect of wall heat flux (cases 2 and 15–19) on annular flow boiling (G = 60 kgm⁻²s⁻¹, D = 0.5 mm, $P_{\infty }$  = 160 kPa, x = 0.1): (a) boiling and film conduction heat transfer coefficients; and (b) main frequency and the corresponding amplitude of the interfacial waves.

Comparison with correlations

In this section, the CFD results are compared with some of the existing correlations for the prediction of film thickness and the boiling heat transfer coefficient.

Liquid film thickness

Bang et al. ⁵² calculated the film thickness using the gas void fraction determined by the Lockhart–Martinelli approach. In their model, instead of using the conventional turbulent–turbulent form for the Martinelli parameter, an empirical expression proposed by Yu et al. ⁵³ was used. In addition, Bao et al. ⁵⁴ showed the CISE approach, recommended by Hewitt, ⁵⁵ provided good estimates of void fractions in two-phase microchannel flows. In this section, the film thickness is calculated following the procedure of Bang et al., ⁵² with the gas void fraction determined by the Lockhart–Martinelli and CISE approaches above.

Figure 17 shows that the Lockhart–Martinelli approach predicts the CFD film thickness within ±20% error, whereas the CISE approach underestimates some data points by more than 20%. The increasing/decreasing trends in the parametric groups are captured correctly except for the wall heat flux cases where both methods predict a constant film thickness value because the heat flux effect is not included in either of them. The ratio of film thickness calculated in the CFD simulations to the predicted value ( ${\delta }_{\text{F},\text{CFD}}/{\delta }_{\text{F},\text{pred}}$ ) yields an average value of 1.07 (standard deviation 0.12, 18 points) for the Lockhart–Martinelli approach and an average value of 1.18 (standard deviation 0.11, 18 points) for the CISE approach. The better performance of Lockhart–Martinelli approach is in accord with the comparison between experimental data and various correlations performed by Bao et al. ⁵⁴

Figure 17.

Comparison of the CFD film thickness with those predicted by: (a) Lockhart–Martinelli, and (b) CISE approaches. Various symbol types are used to denote the values from different parameter groups: base case (▪), mass flux (□), tube diameter (^), system pressure (Δ), vapour quality (▽) and wall heat flux (⋄).

Heat transfer coefficient

A large number of correlations have been developed for the prediction of heat transfer coefficient in boiling flows. Reviews and discussion of the correlations are available in many works, such as Bertsch et al. ⁵ and Kim and Mudawar. ⁵⁶ In the current work, the Kim and Mudawar, ⁵⁶ Kandlikar and Balasubramanian, ⁵⁷ and the Bang et al. ⁵² models are chosen for comparison with the current CFD results, because of their wide usage and validity ranges consistent with the current geometric and flow parameters. The Kim and Mudawar ⁵⁶ model was constructed through reviewing a large number of existing correlations and covers various fluids and parameter ranges. The Kandlikar and Balasubramanian ⁵⁷ model calculates the boiling heat transfer coefficient as an enhanced single-phase laminar heat transfer coefficient. The Bang et al. ⁵² model is a film conduction model with the film thickness being estimated using a modified Lockhart–Martinelli approach. Despite the consideration of both nucleate and convective boiling contributions in the first two correlations, only the convective parts are used here because the current CFD work assumes a completely convective boiling process. The film conduction heat transfer coefficient calculated using the average CFD film thickness is also used in the comparison.

Figure 18 shows the comparison between the CFD boiling heat transfer coefficients with the predicted values from the correlations. The Kim and Mudawar ⁵⁶ correlation (Figure 18(a)) underestimates the CFD boiling heat transfer coefficient for most data points except the low values (<20 kW m⁻²K⁻¹), whereas the Kandlikar and Balasubramanian ⁵⁷ correlation (Figure 18(b)) overestimates the CFD data except for those from two pressure cases. The Bang et al. ⁵² film model (Figure 18(c)) underestimates all the CFD data by around 30%. It is noted that the vapour quality effect is weakly predicted by the Kim and Mudawar ⁵⁶ correlation. More strikingly, the pressure effect predicted by the Kandlikar and Balasubramanian ⁵⁷ correlation and the Bang et al. ⁵² film model is contrary to that for the CFD boiling heat transfer coefficients. As pointed out by Liu et al. ⁵⁸ the correlations are developed based on a large number of experimental datasets obtained with various levels of upstream compressibility (“soft inlet” vs “hard inlet”) and have often been shown to provide poor predictions of any particular dataset. It is therefore not surprising that there are significant discrepancies between the predicted values and the current CFD results (which are obtained with a “hard inlet”).

Figure 18.

Comparison of the CFD boiling heat transfer coefficients with those predicted from correlations: (a) Kim and Mudawar correlation ⁵⁶ ; (b) Kandlikar and Balasubramanian ⁵⁷ correlation; (c) Bang et al. ⁵² film model; and (d) CFD film thickness model. Various symbol types are used to denote the values from different parameter groups: base case (▪), mass flux (□), tube diameter (^), system pressure (Δ), vapour quality (▽), and wall heat flux (⋄).

Notably, the predictions made by the CFD film thickness model (Figure 18(d)) are worse than those made by the Bang et al. ⁵² film model. In fact, it is shown in the ‘Heat transfer enhancement caused by the interfacial waves' section that the boiling heat transfer coefficient is higher than the film conduction value due to the heat transfer enhancement made by the interfacial waves. The average enhancement factor for the film conduction coefficient found in the simulations is 1.48 (standard deviation 0.14, 18 points).

It is evident that none of the correlations or the simple film thickness model can capture the behaviour observed in the CFD simulations. This is hardly surprising given the massive scatter seen when the correlations are compared with experimental data and the simple film model ignores the enhancement arising from the interfacial waves. Work is underway using the apparatus described in Liu et al. ⁵⁸ to obtain experimental data for comparison using a hard inlet.

Conclusions

The CFD model developed in the current work has been applied to the analysis of convective flow boiling heat transfer in laminar, quasi-steady annular flow with interfacial instabilities caused naturally by either phase change or imposed numerical perturbation. Small inlet perturbations have a negligible effect on boiling heat transfer using the current model. Natural instabilities, which are caused by variations in the phase change mass flux at the interface, have shown significant enhancements to the overall heat transfer coefficient through the advection of interfacial waves that continuously reduce the film thickness and generate flow recirculation around them. The time-averaged boiling heat transfer coefficients obtained in the quasi-steady region have shown significant enhancement compared with the values obtained from a steady film conduction model.

Parametric analysis using the numerical model have shown positive effects of increasing mass flux, system pressure, vapour quality and wall heat flux, but negative effects of increasing channel diameter on the heat transfer coefficient. Comparisons have shown reasonable agreements between the CFD results and some existing correlations for the film thickness and the boiling heat transfer coefficient, even though some of the correlations were developed based on large number of experimental data obtained with unknown upstream compressibility. However, it failed in capturing effects of some parameters, such as vapour quality and pressure. A relatively constant enhancement factor (∼1.5) is obtained between the boiling model and the film conduction heat transfer coefficient.

Although the current model is only 2D and applies only to laminar flow, it has significant potential for the analysis of detailed mechanisms of microchannel flow boiling and the effects of various flow parameters. This model is also an excellent basis which could be extended to 3D to study asymmetric behaviour of boiling annular flows, such as wave crest stripping of droplets, and their effects on microchannel heat transfer.