CFD Two Fluid Model for Adiabatic and Boiling Bubbly Flows in Ducts

1. Introduction

The predictions of the CFD Two Fluid Model for bubbly flows in ducts depend on constitutive relations for various interfacial and wall transfer terms. The objective is to use state-of-the-art data covering a wide range of experimental conditions to validate the constitutive relations for two cases of practical interest: adiabatic flows and subcooled and saturated boiling flows. The key to the validation is that all the data was obtained with either double conductivity or double optical needle probes to measure the bubble sizes.

This chapter is divided into two parts. The first part considers adiabatic flows where the interfacial and wall momentum transfer is represented by various forces acting between the two phases and between the wall and each phase. There is a considerable amount of work on the development of constitutive relations for these forces in the open literature. However, these relations have been derived under special flow conditions which may not hold in general. Hence, every force term needs to be calibrated and this introduces the closure coefficients which are dependent on the flow conditions. The air-water experimental validation and calibration of lift, wall and turbulent diffusion forces, which determine the distribution of void fraction in the direction normal to the flow, is carried out.

The second part describes the implementation and validation of heat and mass transfer models required for boiling flows. In particular the wall heat transfer model is based on splitting the wall heat flux into three components: single phase convection, evaporation at the bubble surface and quenching after bubble departure. The evaporation and quenching components depend on three parameters: bubble nucleation site density, bubble departure diameter, and bubble departure frequency. These parameters depend on complex physical phenomena that are not completely understood. However, a considerable amount of experimental data, scaling and various models exist for each one. Nevertheless, it is impossible to test these three separate phenomena simultaneously and in particular the bubble departure diameter remains an open problem (Colin et al., 2009). Therefore only the other two parameters were validated and the validation was performed with two data sets where the bubble diameters were measured: one with high pressure R12 and another with water at atmospheric pressure. The high pressure Freon data corresponds to a fluid-vapor density ratio which is equivalent to that of water-steam at 150 bar. Thus the liquid to vapour density ratio between Freon and water data varies over two orders of magnitude. The surface tension also varies over a wide range from 0.0017 to 0.057 N/m. The ratio of the flow channel hydraulic diameter to the bubble diameter encountered during these simulations varied from 4 to 40. As anticipated, the two-fluid model performs very well when bubbles are small compared to the channel width. However, appropriate modifications have to be made to the two-fluid model to account for cases where the bubbles are too large.

The constitutive relations remain the main limitation of any Two Fluid Model. In particular the present model lacks constitutive relations for the bubble size at the wall and in the bulk. Hopefully, current research such as the development of interfacial area transport models (e.g., Ishii & Hibiki, 2006) will help to resolve this shortcoming in the foreseeable future.

The commercial CFD code ANSYS CFX (version 11 for adiabatic and 12-beta for boiling flows) was used for all the simulations reported in this chapter.

2. Nomenclature

3. Two fluid model mass and momentum equations

Mass Conservation

The ensemble averaged two-fluid model continuity equations (Ishii & Hibiki, 2006) governing the motion of n fields, k = 1 to n, has the following form:

$$\frac{\partial }{\partial t}{\alpha }_k{\rho }_k+\nabla \cdot {\alpha }_k{\rho }_k{\underset{\_}{\overline{\upsilon }}}_k=\Delta m_k$$ E1

For the bubbly flow analyzed during this study, the two-fluid model is comprised of two fields: liquid continuous (k = 1) and dispersed bubbles (k = 2) and the mass transfer across the interface is zero for adiabatic flows.

Momentum conservation

$$\frac{\partial }{\partial t}{\alpha }_k{\rho }_k{\underset{\_}{\overline{\upsilon }}}_k+\nabla \cdot {\alpha }_k{\rho }_k{\underset{\_}{\overline{\upsilon }}}_k{\underset{\_}{\overline{\upsilon }}}_k=-{\alpha }_k\nabla p_k+\nabla \cdot {\alpha }_k\left({\underset{\_}{\underset{\_}{\tau }}}_k+{\underset{\_}{\underset{\_}{\tau }}}_k^{\mathrm{Re}}\right)$$ E2

$$+{\alpha }_k{\rho }_k\underset{\_}{g}+{\underset{\_}{M}}_{ki}+{\displaystyle \sum _j\Delta m_{k,j}{\overline{\underset{\_}{\upsilon }}}_j}$$ E3

The interfacial momentum transfer force comprises of force terms due to drag, turbulent diffusion, lift and a wall force:

$${\underset{\_}{M}}_{ki}={\underset{\_}{M}}_{ki}^D+{\underset{\_}{M}}_{ki}^{TD}+{\underset{\_}{M}}_{ki}^L+{\underset{\_}{M}}_{ki}^W$$ E4

Other forces like the Virtual Mass force and Basset force are also present but are assumed to be negligible.

Drag force

The drag force of the bubbles is given by,

$$M_2^D=-\frac{3}{4}\alpha {\rho }_1\frac{C_D}{D_b}\left|{\overline{\underset{\_}{\upsilon }}}_R\right|{\overline{\underset{\_}{\upsilon }}}_R$$ E5

The drag coefficient is given by the Ishii-Zuber correlation (1979). The relative velocity is given by $${\overline{\underset{\_}{\upsilon }}}_R={\overline{\underset{\_}{\upsilon }}}_2-{\overline{\underset{\_}{\upsilon }}}_1$$

Lift force

The lift force is given by (Auton, 1987):

$${\underset{\_}{M}}_2^L=-C_L{\rho }_1\alpha \left({\upsilon }_2-{\upsilon }_1\right)\times \left(\nabla \times {\upsilon }_1\right)$$ E6

The lift coefficient used for the calculations is C _L = 0.1 (Lopez de Bertodano, 1992). This value is of the same order as the experimental value (Tomiyama et al., 2002) for a bubble in Couette flow (C _L = 0.288).

Wall force

The wall force model accounts for the effect that keeps the centers of the bubbles no closer than approximately one bubble radius from the wall. This force is important when the lift force is present. It is given by (Antal et al., 1991):

$${\underset{\_}{M}}_2^W=C_{wall}{\rho }_1\alpha {\left|{\upsilon }_2-{\upsilon }_1\right|}^2\underset{\_}{n},C_{wall}=\min \left\{\mathrm{0,}-\left(\frac{c_{w1}}{D_b}+\frac{c_{w2}}{y_{wall}}\right)\right\}$$ E7

where $$\underset{\_}{n}$$ and $$y_{wall}$$ are the normal inward vector and the distance from the wall respectively. The values for the coefficients that are used in the CFD code CFX are: $$c_{w1}$$ = -0.01 and $$c_{w2}$$ = 0.05. Note that the ratio $$\left|c_{w2}/c_{w1}\right|$$ determines the thickness of a layer next to the wall (in terms of the bubble diameter) in which the wall force is applicable.

Turbulent Diffusion force

The turbulent diffusion force is not yet well understood. The turbulent diffusion force model (Lopez de Bertodano, 1998) is based on a probability distribution function analysis of the two-fluid model (Reeks, 1991) which has the following form:

$${\underset{\_}{M}}_2^{TD}=-({\rho }_2+C_{VM}{\rho }_1)\frac{{\tau }_c}{{\tau }_b}\left(\frac{{\tau }_c}{{\tau }_c+{\tau }_b}\right)\overline{{\underset{\_}{\upsilon }}^{\prime }{}_1{\underset{\_}{\upsilon }}^{\prime }{}_1}\cdot \nabla \alpha$$ E8

The time constant of the bubbles is derived from the equation of motion of a single bubble as:

$${\tau }_b=\frac{4}{3}\frac{({\rho }_2+{\rho }_1{C}_{VM})}{{\rho }_1{C}_D}\frac{D_b}{{\upsilon }_R^{}}\cong \frac{4}{3}\frac{C_{VM}}{C_D}\frac{D_b}{{\upsilon }_R^{}}$$ E9

The time constant of the turbulent eddies is proportional to the turnover time of the eddies and it is calculated according to the model by Mostafa and Mongia (1988):

$${\tau }_{c1}=0.35\frac{k}{\epsilon }$$ E10

The effect of eddy cross over is accounted for by a simple model when the particle is moving very fast with respect to the eddies, the time constant is given by:

$${\tau }_{c2}=\frac{\lambda {}_e}{\underset{\_}{\upsilon }{}_R}$$ E11

where λ _e is the Eulerian length scale of the eddies. The effective time constant is obtained combining equations (9) and (10) as:

$$\frac{1}{\tau {\text{}}_c}=\frac{1}{\tau {\text{}}_{c1}}+\frac{1}{\tau {\text{}}_{c2}}$$ E12

This approximation has the right asymptotic behavior at both extremes when the bubbles are very small or very large. The model given by equations (7) to (11) may be simplified under several hypothetical asymptotic conditions. We consider two cases:

Case 1 - For $${\tau }_c{\tau }_d,{\tau }_{c1}{\tau }_{c2}$$ it can be shown that equation (7) is equivalent to the turbulent diffusion force model obtained by Lopez de Bertodano et. al. (2006):

$${\underset{\_}{M}}_2^{TD}=-\frac{1}{{\tau }_b}({\rho }_2+C_{VM}{\rho }_1){\nu }_2{}^{}\nabla \alpha$$ E13

where ${\nu }_2$ is the turbulent diffusivity of the bubbles. This equation was used by Lopez de Bertodano (2006) to model the diffusion of small bubbles by the very large vortices generated by cap bubbles in cap-bubbly flow. This model scales like $$\upsilon \text{'}(\propto \upsilon )$$ and is adequate when the turbulence is bubble induced. However for shear induced turbulence a diffusion model with stronger flow dependence is needed.

Case 2 - For $${\tau }_{c1}{\tau }_{c2}$$ it can be shown that equation (7) is equivalent to the model of Lopez de Bertodano et al. (1994a) for the diffusion of bubbles in turbulent duct flow, which is implemented in CFX as:

$${\underset{\_}{M}}_2^{TD}=-C_{TD}{\rho }_1{k}_1\nabla \alpha$$ E14

where the turbulent diffusion coefficient used for the present calculations is C _TD = 0.25. This model scales like $$k\propto \upsilon {\text{'}}^2(\propto {\upsilon }^2)$$ which provides more diffusion at higher flow velocities in the case of shear induced turbulence.

Turbulence Transport

The closure for the Reynolds stresses in equation (2) is based on the two-phase k-ε model developed by Lopez de Bertodano et al. (1994b). This model assumes that the shear induced (SI) and bubble induced (BI) turbulent stresses are added together:

$${\underset{\_}{\underset{\_}{\tau }}}_k^{\mathrm{Re}}={\left({\underset{\_}{\underset{\_}{\tau }}}_k^{\mathrm{Re}}\right)}_{SI}+{\left({\underset{\_}{\underset{\_}{\tau }}}_k^{\mathrm{Re}}\right)}_{BI}$$ E15

The equivalent expression for the turbulent diffusivity of momentum of the liquid phase is given by Sato et al. (1981):

$${\nu }_1=C_{\mu }\frac{k^2}{\epsilon }+C_{DB}\alpha D_b{\upsilon }_R$$ E16

where the first term on the RHS corresponds to the k-ε model for the shear induced diffusivity and the second term corresponds to Sato et al.’s (1981) model for the bubble induced diffusivity. The coefficient C _μ = 0.09 is the standard value for the k-ε model. However, the coefficient C _DB can be adjusted to fit the velocity profile data. A value of 0.6 is used for C _DB in the present simulations. Standard k-ε model transport equations are used for the continuous phase. For the dispersed phase, the kinematic eddy viscosity is assumed to be equal to that of continuous phase, and hence no turbulence transport equations are needed.

4. Adiabatic CFD model assessment with atmospheric air-water data

Experimental data and CFD model

Serizawa et al. (1986) performed measurements of void fraction and interfacial area concentration in upward air water bubbly flow in a 30 mm diameter tube. A double conductivity probe was used to measure void fraction and interfacial area concentration. The radial profiles of void fraction and interfacial area concentration were obtained at 2.5 m axial distance from the inlet. The flow rates covered the range j_L = 0.5 - 5 m/s and j_G = 0 - 0.54 m/s.

The data set was used to compare the spatial distribution of void fraction. In particular the two-fluid model assessment was done with Experiment III data for j_L ≤ 2 m/s, which was the range of interest for the present study.

The CFD calculations were performed with CFX-11. Geometric modeling and meshing was done using the meshing software GAMBIT. Following the CFX guidelines a one-node, 2-degree slice of the pipe was adopted as the solution domain as shown in Fig. 1. A hexahedral mesh with 30 uniform radial elements and 150 uniform axial elements was selected after performing a mesh sensitivity study. The mesh is shown in Fig. 1.

Figure 1.

The Computational Domain and Typical CFX Mesh

Turbulent Diffusion Force Model Comparison

The two turbulent diffusion force models described in the previous section were compared. The effect of the turbulent diffusion models is shown in Fig. 2. The two models available in CFX correspond to the models discussed in the section describing the turbulent diffusion force of the two-fluid model. The model of Lopez de Bertodano et. al. (2006) (i.e., Case 1) is similar to the Favre Averaged Drag Model of CFX and the Lopez de Bertodano Model of CFX is the model discussed in Case 2.

Figure 2.

Turbulent Diffusion Model Comparison for j_L = 1 m/s and 2 m/s

The results with these two models did not differ much at j_L = 1 m/s. A typical result at 1 m/s is shown in Fig. 2. However, Fig. 2 shows that at 2 m/s the k dependent model (i.e., Case 2) with C _TD = 0.25 is a significant improvement over the Favre type model (i.e., Case 1). The previously discussed inadequacy of the Case 1 model to produce sufficient diffusion for shear induced turbulence is evident in the figure and therefore the Lopez de Bertodano Model of CFX was selected for all further simulations.

Wall Force Model Comparison

Calculations were performed for two sets of wall force parameters, viz. Antal et al. (1991) (C _w1 = -0.01, C _w2 = 0.05) and Krepper and Prasser (2000) (C _w1 = -0.0064, C _w2 = 0.016). The uncertainties in the wall force are very significant and these two sets of parameters are used as the two bounds. The bubble size used for the calculations is D_b = 5 mm., which is bigger than the reported value of 3 mm. This is done so because the bubbles are oblate and the measurements correspond to the smaller diameter, whereas the distance from the wall is determined by the larger diameter.

Figures 3 to 5 show the prediction obtained using the two wall force models. A general trend noticed here is that the averaged experimental void fraction is higher than the computed results for all cases and the deviation is more for j_L = 1.5 m/s and 2 m/s cases. This has a significant impact and it is the main reason why data is consistently higher than the calculated results in the center of the pipe. One explanation for this systematic discrepancy could be that the volumetric gas flow rates were measured at the inlet pressure instead of atmospheric pressure, and in such case correcting the computations for the effect of pressure would give better results. However, this has not been done here.

It can be seen from Figures 3-5 that the effect of the wall force coefficients is significant both in the prediction of the wall peak and in the void distribution in the far wall region. In general, the wall peak data is bounded by Antal et al.’s coefficients from below and Krepper and Prasser’s coefficients from above.

Fig. 3 shows comparisons for j_L = 1.0 m/s. At this lower liquid flux, the void distribution is primarily controlled by the lift and the wall forces. The higher void fraction case (Fig. 3) shows that there is a wall peak in the near wall region and relatively uniform void distribution in the far wall region which implies that the turbulent diffusion force has a

Figure 3.

Wall Force Model Comparison (j_L = 1 m/s)

negligible effect. It is clear from the comparisons that Krepper and Prasser’s coefficients underestimate the void fraction in the region away from wall.

For j_L = 1.5 m/s, Krepper and Prasser’s coefficients predict the wall peaks accurately (Fig. 4) but the void fraction in the far wall region is under predicted. For the j_G = 0.083 m/s, the far wall region void fraction is well predicted with Antal et al’s coefficients.

Figure 4.

Wall Force Model Comparison (j_L = 1.5 m/s)

For j_L = 2 m/s, the data shows that the void peak moves further away from the wall than that in j_L = 1.5 m/s. This effect cannot be predicted by the present model without making some significant changes to the lift coefficient and it is an unresolved problem of two-phase flow CFD as of now. The void peak magnitude and bulk void fraction is better predicted by Antal et al.’s coefficients (Fig. 5).

Figure 5.

Wall Force Model Comparison (j_L = 2 m/s)

Thus, Antal et al.’s coefficients predict the wall peak magnitude and the trend of the void distribution in the far-wall region quite well. Furthermore, the far-wall results indicate that the turbulent diffusion force model is satisfactory. As a result of this study, following models were selected for the further analysis using One-Group IATE.

Turbulence Diffusion Force: Lopez de Bertodano Model of CFX (C _TD = 0.25)

Wall Force Model Coefficients: Antal et al’s (1991) coefficients. (c _w1 = -0.01, c _w2 = 0.05)

Lift Force: C _L = 0.1.

5. Two fluid model energy equation

Energy conservation

$$\frac{\partial }{\partial t}{\alpha }_k{\rho }_k{e}_k+\nabla \cdot {\alpha }_k{\rho }_k{\underset{\_}{\overline{\upsilon }}}_k{e}_k=-\nabla \cdot \left[{\alpha }_k\left(\underset{\_}{{q^{″}}_k}+{\underset{\_}{{q^{″}}_k}}^{\mathrm{Re}}\right)\right]$$ E17

$$+p\left[\frac{\partial {\alpha }_k}{\partial t}+\nabla \cdot \left({\alpha }_k{\underset{\_}{\overline{\upsilon }}}_k\right)\right]+\Delta m{}_k{H}_{ki}+q_{ki}$$ E18

$$\underset{\_}{{q^{″}}_k}=-{\lambda }_k\nabla T_k,{\underset{\_}{{q^{″}}_k}}^{\mathrm{Re}}=-\frac{{\rho }_k{\nu }_{tk}}{{\mathrm{Pr}}_t}\nabla e_k$$ E19

where, λ _k is the thermal conductivity and ν _tk is the turbulent viscosity for phase k. Pr _t is the turbulent Prandtl number, H _ki is the mean enthalpy of phase k near interfaces, and q _ki is the interfacial heat flux whose closure is described in detail later. In the present calculations of subcooled boiling, the vapor generated was assumed to be at saturated temperature and hence was treated as isothermal phase.

Interfacial heat transfer

The interface to liquid heat transfer (k=l) is expressed as,

$$q_{li}=h_{li}{a}_i\left(T_{sat}-T_l\right)$$ E20

where, h _li is the heat transfer coefficient between the liquid and the interface, closed as follows,

$$h_{li}=\frac{{\lambda }_l}{d_s}Nu$$ E21

In the above equation, λ _l is the liquid thermal conductivity and d _s is a length scale which is assumed equal to the bubble diameter. The Nusselt number (Nu) in the above closure equation is given by following expression (Ranz & Marshall, 1952):

$$Nu=2+0.6{\mathrm{Re}}_b^{0.5}{\mathrm{Pr}}_l^{0.33}$$ E22

Reynolds number Re _b used in the above closures is defined as,

$${\mathrm{Re}}_b=\frac{D_b\left|{\underset{\_}{\overline{\upsilon }}}_R\right|}{{\nu }_l}$$ E23

The interfacial mass transfer is given as,

$${\dot{m}}_{cond,L}=\max \left(\frac{h_{li}\left(T_{sat}-T_l\right)a_i}{H_{\lg }},0\right),{\dot{m}}_{evap,L}=\min \left(\frac{h_{li}\left(T_{sat}-T_l\right)a_i}{H_{\lg }},0\right)$$ E24

where, H _lg is the latent heat of phase change.

Wall heat transfer

Wall heat transfer is the most important aspect of boiling as it provides the sources for energy and phase mass balance equations. Most wall heat transfer models are based on partitioning the wall heat flux into three components. In particular, the multidimensional CFD model by Kurul and Podowski (1990) has been implemented in the CFX-12 beta version (Egorov & Menter, 2004) and some of the details of this implementation are reported in the next section.

The three components of the wall heat flux are:

Evaporation heat flux (q _e ) – The part of heat flux utilized in formation of vapor at the wall,

Forced convection to the liquid (q _c ),

Wall quenching by liquid phase (transient conduction) (q _q ) – This accounts for the heat transfer to the subcooled liquid that replaces the detached bubble at the wall.

The CFD Wall Boiling Model

The wall surface is assumed to be split into two parts (A ₁, A ₂) each under the influence of one phase. Fraction A ₂ is influenced by the vapour bubbles formed on the wall and participates in the evaporation and quenching heat transfer. Fraction A ₁ is the remaining part of the wall surface, (A ₁=1A ₂) and participates in the convective heat transfer to the liquid. A ₁ and A ₂ are related to the nucleation site density per unit wall area (n _site ) and to the influence area of a single bubble forming at the wall nucleation site. The CFX-12 model assumes, that the diameter of the bubble influence zone is twice as big as the bubble departure diameter (a = 2).

$$A_2=\frac{\pi }{4}{\left(a\cdot D_{dep}\right)}^2\cdot n_{site},A_1=1-A_2$$ E25

Evaporation Heat Flux:

The evaporation heat flux is obtained from the rate of vapor generation at the wall which is given as a product of mass of a bubble at detachment, the detachment (departure) frequency and the nucleation site density,

$${\dot{m}}_e=\frac{\pi D_{dep}^3}{6}{\rho }_G{f}_{dep}{n}_{site}$$ E26

and,

$$q_e={\dot{m}}_e{h}_{\lg }$$ E27

The specification of evaporation heat flux using above equation requires closure for bubble departure diameter, bubble departure frequency and the site density. The default model in CFX uses following closure relations:

The bubble departure diameter is given by the empirical correlation of Tolubinsky and Kostanchuk (1970), the Nucleation Site Density is obtained using the correlations of Lemmart and Chawla (1977) and the departure frequency is closed using the simplest available expression (Cole, 1960), which uses a characteristic bubble velocity (terminal velocity of bubble rise) and a characteristic bubble size (departure diameter).

The frequency of bubble departure obtained by Cole (1960) was derived under the assumption of bubble moving with terminal velocity once detached from the surface. This is a fair assumption for pool boiling scenario at low heat fluxes where bubble detachment is hydrodynamically governed, but may not be accurate for flow boiling at high heat fluxes where the thermodynamics governs the bubble growth and detachment. The closure relations are given below

$$D_{dep}=D_{ref}\text{exp}\left(-\left(T_{sat}-T_{liq}\right)/T_{ref}\right),D_{ref}=0.6mm\text{, }{T}_{ref}=\text{ 45 }K$$ E28

$$n_{site}=n_{ref}{((T_{wall}-T_{sat})/\Delta T_{ref})}^{1.805},n_{ref}=7.9384\times {10}^5{m}^{-2}\text{, }\Delta T_{ref}=\text{ 10 }K$$ E29

$$f_{dep}=\sqrt{\frac{4g\Delta \rho }{3D_{dep}{\rho }_{liq}}}$$ E30

An alternative to the Lemmart and Chawla (1977) correlation for the nucleation site density is the correlation of Kocamustaffaogullari and Ishii (1983) validated against water data for 0-50 bar. The Kocamustaffaogullari-Ishii (K-I) correlation relates the site density to the wall superheat as follows:

$$n_{site}^{*}=R_c^{*}{}^{-4.4}f\left({\rho }^{*}\right)$$ E31

where,

$$n_{site}^{*}=n_{site}{D}_{dep}^2$$ E32

$$R_c^{*}=R_C/(D_{dep}/2)$$ E33

R _c is the critical radius of the surface cavity which represents a minimum cavity size which can be activated at a given wall superheat. R _c is given as,

$$R_C=\frac{2\sigma \left(1+\left({\rho }_g/{\rho }_f\right)\right)/P}{\left(\exp \left(h_{fg}\Delta T_{sat,e}/\left(R{T}_g{T}_{sat}\right)\right)-1\right)}$$ E34

Under the following conditions:

$${\rho }_g{\rho }_l,h_{fg}(T_g-T_{sat})/\left(R{T}_g{T}_{sat}\right)1$$ E35

R _c can be simplified as,

$$R_C=2\sigma T_{sat}/\left({\rho }_g{h}_{fg}\Delta T_{sat,e}\right)$$ E36

The density dependent parameter in equation (29) is given as,

$$\begin{array}{l}f\left({\rho }^{*}\right)=2.157\times {10}^{-7}{\rho }^{*}{}^{-3.2}{\left(1+0.0049{\rho }^{*}\right)}^{4.13}\\ {\rho }^{*}=\left({\rho }_l-{\rho }_g\right)/{\rho }_g\end{array}$$ E37

The effective wall superheat in equation (32) is usually less than the actual wall superheat. The reason for this is as follows. A bubble nucleated at the wall grows through a liquid film adjacent to the wall where a considerably high temperature gradient exists. So, in reality, it experiences a lower mean superheat than the wall superheat. In case of pool boiling, the difference is not significant and the superheat based on the wall temperature can be taken as the effective superheat. In case of forced convective boiling, there exists a steeper temperature gradient in the near wall field due to thermal boundary layer. Hence, the effective superheat experienced by the bubble is smaller than the actual wall superheat and it is given as,

$$\Delta T_{sat,e}=S\Delta T_{sat}\text{where},\Delta T_{sat}=T_w-T_{sat}$$ E38

The multiplier S in the above equation is the superheat suppression factor and is given as (Chen, 1966),

$$S=1/\left(1+1.5\times {10}^{-5}{\mathrm{Re}}_{TP}\right)$$ E39

The two-phase Reynolds number is given as,

$${\mathrm{Re}}_{TP}=\left[G\left(1-x\right)d_h/{\mu }_l\right]F^{1.25}$$ E40

$$\begin{array}{l}F=1.0,for{X}_{tt}\ge 10\\ =2.35{\left(0.213+1/X_{tt}\right)}^{0.736},for{X}_{tt}10\end{array}$$ E41

The Martinelli parameter can be approximated as,

$$X_{tt}={\left(\frac{1-x}{x}\right)}^{0.9}{\left(\frac{{\rho }_g}{{\rho }_l}\right)}^{0.5}{\left(\frac{{\mu }_f}{{\mu }_g}\right)}^{0.1}$$ E42

Thus, equations (29)-(40) represent the complete model of Kocamustafaogullari-Ishii (1983), which has been implemented in CFX-12 for the present study.

Quenching Heat Flux:

As mentioned previously, the quenching heat flux is the component of wall heat flux utilized to heat the cold liquid replacing the detached bubble adjacent to the heated wall. In order to evaluate this component, Mikic and Rohsenow (1969) used an analytical approach starting with transient conduction in semi-infinite medium with heated wall being the only boundary where temperature is specified.

$$\frac{\partial T(x,t)}{\partial t}=\alpha \frac{{\partial }^{2}T(x,t)}{\partial x^2};T(x\mathrm{,0})=T_{\infty },T(\mathrm{0,}t)=T_w,T(\infty ,t)=T_{\infty }$$ E43

The above equation has a solution of the form,

$$\frac{T_w-T}{T_w-T_{\infty }}=erf\left(\frac{x}{2\sqrt{\alpha t}}\right)$$ E44

The heat flux at the wall boundary can be obtained as,

$$q^{″}\left(\mathrm{0,}t\right)=-{\lambda }_l{\frac{\partial T}{\partial x}|}_{(\mathrm{0,}t)}=\frac{{\lambda }_l\left(T_w-T_{\infty }\right)}{\sqrt{\pi \alpha t}}$$ E45

The above heat flux is averaged over a “characteristic” time period which is approximated by the inverse of bubble departure frequency. This leads to the original closure for quenching heat flux [Mikic & Rohsenow, 1969].

$$\overline{q^{″}}=f_{dep}{\displaystyle \underset{0}{\overset{1/f_{dep}}{\int }}{q}^{″}(\mathrm{0,}t)dt}=2{\lambda }_l\sqrt{\frac{f_{dep}}{\pi \alpha }}\left(T_w-T_{\infty }\right)=h_q\left(T_w-T_{\infty }\right)$$ E46

Kurul and Podowski (1990) made minor modification to the above expression with an assumption that quenching occurs between detachment of one bubble and appearance of next bubble (nucleation), and this time was assumed to be 80% of the detachment period. This is the final form of quenching heat transfer closure of the CFX-12 model.

$$h_q=2{\lambda }_{l}f\sqrt{\frac{t_w}{\pi \alpha }};t_w=\frac{0.8}{f_{dep}}$$ E47

$$q_q=A_2{h}_q\left(T_w-T_{\infty }\right)$$ E48

Single Phase Convection Heat Flux:

This is evaluated under the standard assumption of a logarithmic temperature profile across the turbulent boundary layer (Kader (1981), ANSYS CFX Solver Theory Guide (2006)).

$$q_{w\mathrm{,1}ph}=\frac{{\rho }_l{C}_{p,l}{u}^{*}}{T^{+}}\left(T_w-T_{liq}\right)$$ E49

$$\begin{array}{l}{T}^{+}={\mathrm{Pr}}_l{y}^{*}{e}^{-\Gamma }+\left(2.12\ln \left(y^{*}\right)+\beta \right)e^{-1/\Gamma };\beta =f({\mathrm{Pr}}_l),\Gamma =f(y^{*})\\ u^{*}=C_{\mu }^{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$4$}\right.}{k}_l^{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.},y^{*}={{\rho }_l{u}^{*}\Delta y/\mu }_l\end{array}$$ E50

where, Δy is the distance of the wall adjacent node from the wall, k _l is the liquid turbulence kinetic energy. In the above equation β is a function of liquid Prandtl number and Γ is a function of y ^* (ANSYS CFX Solver Theory Guide (2006)).

The convective heat flux component through the liquid area fraction is thus given as,

$$q_c=A_{1f}{h}_c\left(T_w-T_{liq}\right),h_c=\frac{{\rho }_l{C}_{p,l}{u}^{*}}{T^{+}}$$ E51

It should be noted here that all the wall heat flux model correlations above use a liquid temperature (T _liq ). The single phase heat transfer correlation above uses the wall adjacent node temperature. However, in the bubble departure diameter and the quenching heat transfer correlations which were based on one-dimensional models, the liquid temperature refers to the bulk mean temperature. As a good estimate, CFX approximates this temperature with the temperature at a fixed y^* (250). This is done in order to have a mesh size independent evaluation of wall heat flux partitions.

6. Boiling CFD Model assessment with high pressure R12 data

Experimental data and CFD model

Morel et al. (2003) performed boiling experiments in a vertical pipe of 19.2 mm internal diameter and a length of 5 m. The pipe had a heated section of 3.5 m preceded by and followed by unheated lengths of 1 m and 0.5 m respectively. The input parameters for the four tests reported were as summarized in Table 1. The liquid-vapor density ratio in these experiments corresponds to that for water-steam at 95-150 bar.

Radial Profiles of Void Fraction, Liquid Temperature and Interfacial Area Concentration were measured. Dual sensor fiber optics probes were used to measure void fraction and interfacial area concentration. Miniature thermocouples were used to measure the fluid temperature.

The CFD calculations were performed with CFX-12. Geometric modeling and meshing was done using the meshing software GAMBIT. A 2 wedge of the pipe is used as the domain. The mesh had 20 radial and 250 axial uniformly spaced elements. Fig. 6 shows the domain and cross-section of the mesh used.

Parameter	Test1 (Deb5)	Test2 (Deb6)	Test3 (Deb13)	Test4 (Deb10)
Mass Flux (kg/m ^2 )	1986	1984.9	2980.9	2027.0
Pressure (bar)	26.15	26.15	26.17	14.59
Inlet Subcooling ( ^o C)	18.12	16.11	18.12	23.24
Heat flux (W/m ^2 )	73.89	73.89	109.42	76.24

Table 1.

Simulation input parameters for R12 tests

Figure 6.

Domain and Mesh for R12 boiling simulation case

Comparisons

All the data sets include radial profiles of bubble diameter. Note that with Freon at these test pressures, the surface tension is 0.0017 N/m for 26 bar which results in very small bubble sizes (i.e., 0.5 mm). Instead of using bubble size correlations for the CFD calculations the mean of the measured radial profile value is used for the bubble diameter in the bulk, whereas, the value measured at the wall is used as the departure bubble size (Fig. 7).

It was found during all the cases that a better match with experimental data is obtained with the coefficient of the turbulent diffusion force C _TD = 0.5 instead of the value obtained for adiabatic air water flows, C _TD = 0.25. This change may be attributed to bubble diameters in these cases which are one order of magnitude smaller than that in atmospheric air-water systems (Prabhudharwadkar et al., 2009). Therefore there are more interactions of the bubbles with smaller turbulent eddies. Two comparisons are presented below for this particular benchmark problem: the nucleation site density model comparison and the bubble size model comparison.

Fig. 8 shows the prediction of radial variation of void fraction obtained using two different site density models. The Lemmart-Chawla (1977) model which is used in CFX was found to under-estimate the void fraction in all the cases. Kocamustafaogullari-Ishii model (termed as

Figure 7.

Bubble diameter approximation for the simulations

Figure 8.

Void Fraction prediction using the two Nucleation Site Density models

K-I Model in the figures below) predicted the data well in all the cases. The prediction however deteriorated at low pressure data (P = 14.59 bar) and this needs further investigation.

Fig. 9 shows the prediction of radial profiles of interfacial area concentration (IAC). The IAC is obtained using the expression for the Sauter mean diameter:

$$a_i=6\alpha /D_b$$ E52

As the bubble diameter is assumed to be uniform, the deviation of the CFD values from the experimental data near the wall, where the bubble diameter is smaller, is justified.

Fig. 10 shows prediction of liquid temperature for the first two cases (for which data is available). The simulation results show a good match with the experimental data in the bulk. However, towards the wall, the results with Lemmart-Chawla correlation showed steeper rise in temperature which caused over-prediction in this region.

Fig. 11 shows the bubble diameter prediction using the CFX default bubble departure diameter model with the Kurul and Podowski (1990) correlation for the bubble size in the bulk of the flow. The bubble departure diameter and the bulk diameter predicted by the models are approximately three times higher than the actual value, leading to a very high rate of vapor generation. This results in over-prediction of void fraction as shown in Fig. 11.

Figure 9.

Interfacial Area Concentration prediction using the two Nucleation Site Density models

Figure 10.

Liquid Temperature prediction using the two Nucleation Site Density models

Figure 11.

Bubble diameter prediction using Kurul-Podowski model and void fraction prediction.

7. Subcooled boiling CFD model assessment with atmospheric water data

Experimental data and CFD model

Lee et al. (2002) performed subcooled boiling experiments in a vertical annulus having a 19 mm diameter heating rod placed inside a tube with internal diameter of 37.5 mm. The total length of the test section was 2.37 m. The central heated length of 1.67 m was preceded and followed by unheated portions of 0.188 m and 0.518 m respectively as shown in Fig. 12.

Figure 12.

Domain and Mesh for water boiling simulation case

Yeoh and Tu (2005) have reported two cases for this experiment where the bubble diameter measurements have been done. These two cases were selected for the analysis here. The input parameters for the two tests reported are summarized in Table 2 below.

Parameter	Test1	Test2
Mass Flux (kg/m 2 )	474	1059.2
Pressure (bar)	1.42	1.43
Inlet Subcooling ( o C)	13.11	17.82
Heat flux (kW/m 2 )	152.9	251.5

Table 2.

Simulation input parameters for water boiling tests

Radial profiles of void fraction, vapor velocity and liquid velocity were measured at 1.61 m downstream of the start of the heated section. Two-point conductivity probe was used to measure void fraction and bubble velocity. Pitot tube was used for liquid velocity measurements with a correction factor for two phase flow.

The CFD calculations were performed with CFX-12. Geometric modeling and meshing was done using the meshing software GAMBIT. A 2 wedge of the annulus is used as the domain. The mesh had 20 radial and 142 axial uniformly spaced elements. Fig. 12 shows the domain and cross-section of the mesh used.

The bubble diameter measurements for the two cases are shown below in Fig. 13. It must be noted here that although a variation in the bubble size has been measured along the radial direction, however, the bubble layer is one bubble thick and the bubble size is of the order of 4-5 mm which occupies almost half of the annulus gap. This means the departure size can be assumed the same as the bubble size, which is equal to the maximum size of bubble measured across the cross-section. Hence, a bubble size (bulk and departure) of 4.7 mm is used in the low heat flux case and that of 4.1 mm is used in high heat flux case.

Figure 13.

Bubble Diameter data for the water boiling tests

The site density model comparison was done for these two cases also and the void fraction comparison is shown below in Fig. 14. It can be seen that the Lemmart-Chawla correlation over-predicted the void fraction at the wall by a significant margin. Although results using K-I model did not accurately match the void fraction distribution, the wall void fraction was consistently better predicted.

Figure 14.

Void Fraction prediction using different site density models (Water Data)

It is to be noted here that the K-I model did a better job even at the high pressure data and hence appears to be better than the Lemmart-Chawla model.

Modifications of the Two-Fluid model

As stated earlier, in this problem, the bubble size is comparable to the annulus gap width and hence further modifications were made to the two-fluid model to improve the results:

Addition of a wall drag force: The bubble moves touching the wall and hence experiences the shear due to the liquid velocity gradient near the wall. In such cases the bubble is further slowed by an additional wall drag which is essentially an additional shear stress on the bubble which is moving in the liquid shear layer. A model for this wall drag force has been derived by Lahey et al. (1993) and is given as follows:

$${\underset{\_}{M}}_2^{WD}=-{\underset{\_}{M}}_1^{WD}=-C_{DW}\frac{{\rho }_l}{100D_b}\alpha \left|{\underset{\_}{\overline{\upsilon }}}_l\right|{\overline{\underset{\_}{\upsilon }}}_g$$ E53

The value suggested for C _DW is 1.

Exclusion of Turbulent Dispersion: The turbulent dispersion force accounts for transverse dispersion of bubbles due to interaction with turbulent eddies. However, in the present case where there exists only one bubble across the cross-section, such a diffusion model is invalid and hence this force is turned off.

Fig. 15 shows the void fraction prediction for the two cases where the wall drag is included. It can be noticed that the wall drag effect is more in the high mass flux case. This is due to the fact that wall drag is comparable to the interfacial drag in this case.

Fig. 16 shows the void fraction predictions after turning off turbulent dispersion (M ^TD ) force. The peak void fraction increased substantially in both cases but it also resulted in slight narrowing of the vapor boundary layer near the boiling wall. Overall, the two modifications to the two-fluid model improved the results. However, these computations based on a model where the void distribution is the result of the balance between the lift and wall forces are not correct in a strict sense because the thickness of the bubble layer is determined instead by the diameter of one bubble. This fundamental problem has been addressed by Moraga et al. (2006).

Figure 15.

Void Fraction predictions after including the wall drag term

Figure 16.

Void Fraction predictions after excluding turbulent dispersion

Fig. 17 shows the vapor velocity comparison with two simulations, viz., the original two-fluid model and with wall drag and no turbulent dispersion. The low flux case showed no significant difference in the three cases. A deviation near the outer wall in the no TD case can be neglected since this region has no vapor. The vapor velocity did not change much even in the high flux case but showed a slight reduction when wall drag was added.

With this benchmark the two-fluid model with the wall-boiling module was validated over a wide range of pressures (~1-150 water-steam equivalent bar). It is worth noting again that the low pressure case presented above has a large bubble size (i.e., half the annulus size) which makes the two-fluid model constitutive relations invalid in a strict sense. Despite this shortcoming, the results are acceptable.

Figure 17.

Vapor velocity predictions without and with two-fluid model modifications

8. Conclusions

Several CFD Two Fluid Model constitutive relations were assessed with experimental data over a wide range of conditions for adiabatic and boiling flows to provide a guideline to choose the constitutive models.

For the air-water adiabatic flow data Auton’s (1987) formulation of the lift force was found appropriate for the prediction of the void fraction profiles. The turbulent diffusion force given by the Lopez de Bertodano et al. (1994a) produced satisfactory results, especially at higher liquid fluxes. Two sets of coefficients for Antal et. al.’s (1991) wall force were compared and it was found that Antal et al.’s coefficients produced more consistent predictions and with the combination of selected lift and turbulent diffusion models, it produced overall better predictions throughout the pipe cross-section. The interfacial momentum models selected using adiabatic data were then used for the boiling flow comparisons with satisfactory results.

For the high pressure R12 boiling flow data the nucleation site density model by Kocamustaffaogullari and Ishii (1983) and the bubble departure frequency model of Cole (1960) were found to predict the data well. For the atmospheric pressure subcooled boiling water data set, where the bubble size is about half the flow channel width and where the thickness of the bubble layer is one bubble diameter, use of a wall drag force and elimination of the turbulent diffusion force improved the predictions.

The present model does not include constitutive relations for the bubble size at the wall nor in the bulk. These remain as open problems to be resolved in the future.