## Abstract

Subcooled flow boiling is widely used as a mode of heat transfer in many industries, especially in nuclear reactors. Despite its advantages, the heat transfer is hampered beyond a certain flux due to a phenomenon known as departure from nucleate boiling (DNB). It is important to determine the void fraction profiles, especially the near-wall void fractions, to evaluate the limiting heat flux conditions. The two-fluid Eulerian model, coupled with the heat flux partitioning model, is widely used to predict subcooled flow boiling characteristics. Over the years, many researchers have not considered lift and wall lubrication forces in their modeling of subcooled flow boiling. Few researchers have considered the Tomiyama model for lift force; however, their results were not encouraging. Moreover, there is no systematic study in evaluating the impact of lift and wall lubrication forces on subcooled flow boiling. In this paper, various lift and wall lubrication models are compared to understand the implications of these forces on void distribution. The advantages and limitations of the models are discussed in detail.

## 1 Introduction

Boiling has been used as a mode of heat transfer in several industries, especially in nuclear reactor core cooling and steam generator cooling. Heat transfer is significantly enhanced in boiling systems compared to single-phase flow. Some degree of subcooled boiling is experienced during the normal operation of pressurized water reactors. However, there exists a limitation on the heat flux beyond which heat transfer deteriorates, leading to failure of the fuel. It is of great importance to determine such limiting heat flux values, which is also known as critical heat flux. At low quality flows, during the subcooled flow boiling, bubbles crowd near the surface, forming a vapor film which in turn hampers the heat transfer, leading to the failure of the heater. This phenomenon is also referred to as departure from nucleate boiling (DNB). Accurate prediction of the near-wall void fraction is essential when it comes to the prediction of DNB. Over the years, computational fluid dynamics (CFD) modeling using the two-fluid Eulerian model, coupled with the heat flux portioning model, has been used to predict subcooled flow boiling characteristics [1–4]. However, the predictive capability of the two-fluid model depends on boiling and momentum closures.

It is essential to understand the implications of different forces like drag force, lift force, wall lubrication force, and turbulent dispersion force, on the subcooled flow boiling, as some might lead to numerical divergence problems and others might give unrealistic solutions due to inapplicable coefficients.

Most of the research is focused on the development of boiling closures [1,5]. Some studies have been performed on the population balance models and their implications on flow boiling predictions [6,7].

Previous studies by Krepper and his colleagues [1,7] indicate that whenever the lift model is considered, the near-wall void fraction is overestimated. According to them, the deviation from the experimental data could be due to a lack of consideration of the swarm effect in models of drag and lift forces. They used the Tomiyama model for lift force and did not consider the wall lubrication force.

Most of the researchers ignored wall lubrication force; only some researchers considered wall lubrication force [8,9]. Colombo and Fairweather [3] did not consider both lift and wall lubrication forces in their models.

Despite many research works related to subcooled flow boiling prediction using the two-fluid model, there is no common consensus among the researchers regarding the momentum closures, especially on the lift and wall lubrication models. Moreover, there is no comprehensive study to understand the influence of lift models on radial void distribution. It is of great importance to have a clear understanding of the applicability and limitations of the available closure models. In this study, various lift and wall lubrication forces models are compared. Limitations of the models are discussed in detail.

## 2 Computational Fluid Dynamics Modeling

The two-fluid Eulerian model, coupled with the wall boiling model, is used for modeling subcooled flow boiling.

### 2.1 Governing Equations.

In the two-fluid Eulerian approach, the continuity, momentum, and energy equations are solved for each phase separately [10].

where $\alpha k$, $\rho k$, $uk$, $qk\u2033$, $Hk$, $\Gamma k$, **M**_{k}, and $Qk$ represent volume fraction, density, velocity, heat flux, specific enthalpy, mass transfer, momentum transfer, and energy transfer, respectively. *P* indicates pressure, $\tau k$ and $\tau kt$indicates viscous and turbulent stresses, respectively.

### 2.2 Wall Partitioning Model.

In subcooled flow boiling, vapor generation occurs at the heater surface, and condensation occurs in the subcooled bulk liquid away from the surface. Kurul and Podowski [11] developed a wall boiling model, where they partitioned the heat flux into three components such as

Evaporative heat flux

Liquid phase convective flux

Quenching heat flux

where $q\u2033total$ is total heat flux, $q\u2033C$ is convective heat flux, $q\u2033Q$ is quenching heat flux, and $q\u2033E$ is evaporative heat flux.

where *h*_{C} is the liquid heat transfer coefficient, $Ab$ is the fraction of the heater surface area influenced by the evaporation process, and $Tw$ and $Tl$ are wall and liquid temperatures, respectively.

*t*

_{w}); according to Kurul and Podowski [11], the waiting period is 80% of the bubble detachment period. Quenching heat flux is given by the following expression [11]:

where $kl$ and $\lambda l$ are the thermal conductivity and diffusivity in the liquid phase.

where $Vb$ is the volume of the bubble based on the bubble departure diameter, $Nw$ is active nucleate site density, $hfg$ is the latent heat of evaporation, and *f* is bubble departure frequency.

Area of influence:

where *K* is the area influence factor.

### 2.3 Boiling Closures

#### Bubble Departure Diameter.

Bubble departure diameter plays an essential role as it affects all three terms (convective, evaporative, and quenching heat flux values). It depends on various parameters like the contact angle, mass flowrate, pressure, and subcooling.

where $dref$ = 0.0006 m, $\Delta Tsub=Tsat\u2212Tl$, and $\Delta Tref$ = 45 K.

In this study, the Tolubinsky and Kostanchuk model is used for the bubble departure diameter. The reference diameter and temperatures are chosen according to the case, based on the experimental data.

#### Bubble Departure Frequency.

#### Nucleation Site Density.

where $Tw$ is the wall temperature and $Tsat$ is the saturation temperature of the liquid.

### 2.4 Momentum Closures

#### Drag Force.

where $\sigma $ is surface tension.

#### Turbulent Dispersion Force.

In this study, the Burns model is used to take into account the turbulent dispersion.

#### Lift Force.

In upward flow, bubbles are pushed toward the wall, and in downward flow, bubbles are pushed away from the wall.

where $Re\omega =\rho l|\u2207\xd7ul|Db2/\mu l.$

**Sr**). The following expression gives the lift coefficient:

The models described above were developed for solid bodies in a shear field, but they have not considered the wobbly nature of the bubble. Tomiyama et al. [23] derived a model for lift coefficient based on their experimental results and found that bubble shape and size have a significant impact on the lift force. They performed experiments on single bubbles under various flow configurations and correlated the lift coefficient based on the Eötvös number. They had observed experimentally that the lift coefficient's sign changed from positive to negative when the bubble size increased, as the bubble deforms.

In Moraga and Tomiyama models, the sign of the lift coefficient changes depending on the flow parameters and bubble size; however, the sign of lift coefficient obtained from Klausner and Legendre models' is always positive irrespective of the bubble size. As a consequence, bubbles always move toward the wall in the upward flow when Klausner and Legendre models are considered.

#### Wall Lubrication Force.

Wall lubrication force is based upon the idea that the liquid velocity between the bubble and the wall is lower than the liquid velocity between the bubble and the core, which results in hydrodynamic pressure difference and drives the bubble away from the wall.

In this study, the coefficients provided by Frank et al. [25] are used (*C*_{w1} = −0.01 and *C*_{w2} = 0.05).

*C*

_{w1}with the second-order term to take into account the influence of bubble motion. Moreover, a functional dependence on the diameter of the pipe to dampen the impact of the force at the center was also introduced. After calibrating the constants with experimental data, they found

*C*

_{w2}to be zero. The final form of the equation is as follows:

Experimental results of Tomiyama [26], Hosokawa et al. [27], and Takemura et al. [28] suggest the presence of wall lubrication force. At low Reynolds numbers, they observed that the bubble drifted away from the wall. Moreover, it is important to note that the force decreased with an increase in Reynolds number. However, if one goes by the recent experiments of Hassan et al. [29], it is evident that there is no liquid film (even if it exists, it is negligible compared to bubble size) between the bubble and the wall, which indicates that the there is no macroscopic force acting on the bubble to push it away from the wall. Moreover, Vaidheeswaran et al. [30] and Lubchenko et al. [31] proposed that the shape of the void fraction profile has to do with the shape of the bubble, not the hydrodynamic pressure gradient.

### 2.5 Turbulence Closures

#### Turbulence Model.

In this paper, the shear stress transport (SST) $k\u2212\omega $ model with enhanced wall function is used, which can deal with *y ^{+}* from 1 to 100 order of magnitude [32]. Zhang et al. [32] have analyzed the performance SST $k\u2212\omega $ model, where they noted that with finer meshes (

*y*close to 1), the predictions are not accurate for two-phase systems involving boiling. Furthermore, Rzehak and Krepper [33], Mali et al. [34] have shown the applicability of the SST $k\u2212\omega $ model to simulate turbulence in boiling two-phase flows.

^{+}#### Bubble Induced Turbulence.

Bubble agitation (movement) induces turbulence in the liquid; this force was first considered by Sato et al. [35]. The most common practice is to add a term to the turbulent viscosity to account for the bubble induced turbulence (BIT).

Another way to consider the bubble effects is by adding source terms to *k* and $\omega $ equation, respectively. In this study, the Sato model is used to take into account the bubble induced turbulence.

## 3 Benchmark Data

To understand the implications of lift and wall lubrication forces, proper validation with experimental data is essential. For validation, three tests are selected from Garnier et al. [36], one each from Roy et al. [37] and Bartolomei and Chanturiya [38]. The experimental details are tabulated in Table 1.

Flow boiling experiments were performed by Garnier et al. [36] over a wide range of pressure, mass flux, and heat flux. The test loop consisted of a vertical heated pipe of 19.2 mm inner diameter and length of 3.5 m. Freon-12 is used working fluid. Void fractions were measured using optical probes at the end of the test section. The motivation for selecting three cases from their work is to understand the influence of lift and wall lubrication forces on wall peaking scenarios.

Roy et al. [37] conducted experiments in a vertical annulus with heated inner surface of length 3.66 m and inner and outer diameters are 15.8 mm and 38.1 mm, respectively. The two-conductivity probe method is used to obtain radial void fraction and vapor velocity profiles. All the measurements were done at an elevation of 1.99 m.

Bartolomei and Chanturiya [38] case is selected to understand the influence of lift and wall lubrication forces on average void fraction along the length of the channel. They performed experiments in vertical tubes of inner diameter 15.4 mm and length of 2 m. At different axial locations, the average void fraction was measured. This is the most commonly used test data for CFD model validation [4,39]

## 4 Base Case Model

In the base case model, Tolubinsky and Kostanchuk's model is used for bubble departure diameter, where reference diameters are taken as 0.35 mm for case-1 and case-2 and 0.5 mm for case-7 (based on the experimental bubble diameter near the wall) and $\Delta Tref$ is taken as 45 °C. However, for cases 4 and 5, the default model is used due to the lack of experimental data. In all the cases, Cole model is used for bubble departure frequency and Lemmert and Chawla model for and nucleation site density. Ishii model is used for the drag force, and the Burns model is used for turbulent dispersion force. Bulk bubble diameter is estimated using Anglart et al. [40] model for all the cases except for case 4. Based on the experimental evidence, Roy et al. [37] reported that the most probable bubble diameter for case 4 was in the range of 0.7 to 1.2 mm. Anglart et al. [40] model underpredicted the bulk bubble diameter, so in this study 0.7 mm is used as bulk bubble diameter. However, it has to be noted that the effect of bulk bubble diameter on the void distribution is not significant. Lift and wall lubrication models are not considered in the base case. SST *k*–*ω* model is used for turbulence modeling along with bubble induced turbulence model by Sato. Ranz and Marshall [41] model is used for interfacial heat transfer. The details of the base case CFD model are tabulated in Table 2. Both liquid and vapor properties are taken from the National Institute of Standards and Technology standard reference database. Two-dimensional axis-symmetric geometries are used for simulations.

## 5 Grid Sensitivity Study

In this work, equidistant grids are used for the simulations. Grid sensitivity analysis is carried out for all five cases using the base model. The influence of axial grid size is studied by using 200, 400, and 600 grid cells while maintaining 20 radial grid cells. Figure 2 shows that the variation in axial grid cells did not impact the radial or axial void distribution. Moreover, to understand the influence of radial grid size, while maintaining the axial grid cells at 400, radial grid cells are varied. Three radial grid cells, 15, 20, and 25, are studied. For all the cases, the *y ^{+}* is in the range of 30–150. Further refinement of mesh led to numerical stability. Figure 2 shows that the variation in grid size did not any significant impact on the radial and axial void distribution.

## 6 Results and Discussions

The results are discussed in two sections, assessing both the impact of lift and wall lubrication forces. Experimental data of the radial void fraction of cases 1 to 4 are compared against the CFD predictions to understand the influence of lift and wall lubrication models on radial void distribution. Moreover, to understand the influence of the models on axial void distribution, the CFD predictions are compared against the axial cross-sectional averaged void fractions. In Fig. 2, along with grid sensitivity, experimental data is also reported. Experimental uncertainties are also plotted in Fig. 2; however, for case-5, the uncertainty data is not available. In cases 1, 2, and 3, the experimental uncertainties involved in void fraction measurement are ±0.02 [36], and in case-4 the experimental uncertainty is ±0.02 when void fraction is greater than 0.1 and ±0.01 when void fraction is less than 0.1 [37].

Moreover, it can also be observed from Fig. 2 that far away from the wall, the value of void fraction in cases 1, 2, and 4 are comparable with uncertainties involved in the experiments.

### 6.1 Lift Force.

The impact of the lift model is studied by comparing various lift models, as shown in Fig. 3. All five cases are simulated with different lift models and are compared with experimental data and the base case model. It can be observed that using Moraga's model for lift force, no significant difference is obtained in the radial void and axial distributions compared to the base case where no lift force is considered, implying that the lift coefficient value obtained using the Moraga model is negligible for the test cases investigated in this study. Notably, in cases 1, 2, and 4, the void fraction is overpredicted away from the wall. The maximum deviation in the near-wall void fraction is about 40% (Fig. 4), when the Moraga model is used, for the investigated cases.

However, when the Tomiyama model is used, it even overpredicted the near-wall void fraction in cases 1, 2, and 4. The deviation from the experimental data in case 1, in the near-wall region is about 250% in case 1 (Fig. 4(a)), even in cases 2 and 4, the deviation in the near-wall void fraction is higher compared to the base case (Figs. 4(b) and 4(d)). Moreover, when the Tomiyama model is considered for the lift force, it was found difficult to achieve convergence, especially in case 5, where the near-wall void fraction in the last computational cell exceeded 0.8.

When the Saffman model is used, the near-wall void fraction deviation increased in cases 1, 2, and 4 (Fig. 4) compared to the base case. Notably, in case 4, the deviation from experimental data is less than 5%. Furthermore, deviation in the core void fraction prediction decreased from 225% (base case) to 150% in case-1 and from 170% (base case) to 30% in case 2 when the Saffman model is considered for the lift force.

Legendre and Magnaudet model followed the same trend as the Tomiyama model; it overpredicted the near-wall and core void fraction (cases 1, 2, and 4).

However, in case-3, where the near-wall void fraction is high, the influence of the lift model is not significant.

It is also clearly evident from the above discussion that, in most cases, the void fraction profile deviated from the experimental data, especially near-wall void fraction when the lift force is considered. It is important to note that the lift coefficients were developed based on experiments with single bubbles. Moreover, when the lift model is not considered, the void fraction predictions are better near the wall, this implies that the present lift models are not fine-tuned for applicability to subcooled flow boiling cases expect for the Moraga model which almost gave the same result as the base case. It is imperative to develop lift models taking into account the swarm effect to have a realistic model for the lift force.

The deviation in the void fraction in the core away from the wall (Fig. 4) is significantly high irrespective of the choice of lift model. It could also be due to a simplistic model of bubble diameter based on liquid subcooling, and Rans and Marshall heat transfer correlation, which may not be appropriate. However, discussion on that is not in the scope of this paper.

### 6.2 Wall Lubrication Force.

To understand the impact of wall lubrication models on void distribution, while considering the Moraga model for the lift force, various wall lubrication models are compared. Three wall lubrication models, i.e., Antal, Tomiyama, and Frank models, are compared (Fig. 5). Figure 5 shows that when the Antal model is used for wall lubrication force, the near-wall void fraction is predicted better than the other cases. Especially in case 2, it captures the shape of the void profile, but underpredicts the void fraction near the wall (Fig. 5(b)). However, a slight positive lift coefficient might improve the near-wall prediction.

Figure 5(c) shows that in the case 3, i.e., when the near-wall void fraction is high and close to 0.8, the wall lubrication model did not have any impact on the radial void distribution. Figure 6(c), which shows the deviation in predicted radial void fractions from the experimental data also indicates that there is no change in deviation when compared with the base case.

Even in case 4, when the wall lubrication force is considered, the void near the wall is underpredicted, which is not conservative (Fig. 6(d)). Tomiyama and Frank wall lubrication models underpredicted the near-wall void fraction in all the cases (Fig. 6). However, the deviation in axial cross-sectional average void fraction improved from around 48% (base case) to 24% (Fig. 6(e)) when the Antal model is considered for wall lubrication force in combination with the Moraga model for the lift force. Furthermore, when the Moraga model is used for lift force, the lift coefficient is slightly negative. When a model like Tomiyama is used for the lift force, for the bubbles of small size, a significant positive lift coefficient is generated, this implies the predictions deviate even further from the experimental data. Improvement in results can be expected when a small value of a positive lift coefficient is used along with the Antal model for wall lubrication.

It is important to note that when wall lubrication models are considered, the near-wall void fraction is underpredicted in all the cases except case 1. However, even in case 1, the void profile is significantly different from that of the experimental data (Fig. 5(a)).

Even though the void fraction prediction near the wall is improved in cases 1 and 2 when the Antal model is considered for wall lubrication force, the existence of such a force in boiling flows, especially at high-pressure conditions, is still debatable.

In summary, the existing models of lift and wall lubrication models are not fine-tuned for subcooled boiling flows. And the variation among the existing models is evident from the results. As noted earlier, most of the models are developed based on experiments on single bubbles/solid particles in a simple shear flow at specific flow conditions. So, the models are not generic and lack the complexities involved in boiling flows.

## 7 Conclusion

A comprehensive analysis is performed to analyze the impact of lift and wall lubrication forces on subcooled flow boiling predictions. For this, some of the most widely used models for lift and wall lubrication forces have been considered in this work, and the implication of these forces on the near-wall void fraction has been studied using five different test data.

The following conclusions are obtained from this study:

Consideration of lift force did not improve the void fraction prediction. Using lift models of Tomiyama et al. [23] and Legendre and Magnaudet [22], the maximum error in near-wall void fraction prediction is found to be more than 200%.

The performance of the lift model by Saffman [21] is better in comparison with Legendre and Magnaudet [22] and Tomiyama et al. [23] models. However, even it overpredicted the near-wall void fraction.

The lift model by Moraga et al. [19] did not affect the radial and axial void distributions for the investigated cases.

Using wall lubrication force model by Antal et al. [24] improved the predictions slightly; however, the existence of such force in the flow blowing situations is still debatable.

Wall lubrication force did not have any considered impact when the near-wall void fraction is high.

Fine-tuning of lift and wall lubrication coefficients is essential to predict the radial void profiles with high accuracy.

## Nomenclature

*A*=_{i}interfacial area (1/m)

*A*=_{b}fraction of the heater surface area influenced by the evaporation process (dimensionless)

*C*=_{D}drag coefficient (dimensionless)

*C*=_{L}lift coefficient (dimensionless)

*C*_{L,}_{lowRe}=lift coefficient at low Reynolds numbers (dimensionless)

*C*_{L,}_{highRe}=lift coefficient at high Reynolds numbers (dimensionless)

*C*_{TD}=turbulent dispersion coefficient (dimensionless)

*C*_{WL}=wall lubrication coefficient (dimensionless)

*C*_{w1}=coefficient (dimensionless)

*C*_{w2}=coefficient (dimensionless)

*C*=_{μ}shear-induced turbulence coefficient (dimensionless)

*D*=_{b}bubble diameter (m)

*D*=_{d}maximum horizontal bubble diameter (m)

*D*=_{bw}bubble departure diameter (m)

*d*_{ref}=reference diameter of bubble (m)

- F
=_{D}drag force (N/m

^{3}) **F**=_{L}lift force (N/m

^{3})**F**_{TD}=turbulent dispersion force (N/m

^{3})**F**_{WL}=wall lubrication force (N/m

^{3})*f*=bubble departure frequency (1/s)

*g*=acceleration due to gravity (m/s

^{2})*H*=specific enthalpy (J/kg)

*h*=_{C}liquid heat transfer coefficient (W/m

^{2}K)*h*=_{fg}latent heat of vaporization (J/kg)

*K*=area influence factor (dimensionless)

*k*=turbulent kinetic energy (m

^{2}/s^{2})*k*=_{l}thermal conductivity of liquid (W/m K)

**M**=interfacial momentum transfer (N/m

^{3})**n**=unit normal vector perpendicular to the wall (dimensionless)

*N*=_{w}nucleation site density (1/m

^{2})*P*=pressure (Pa)

*Q*=energy transfer (W/m

^{3})*q"*=heat flux (W/m

^{2})*r*=radial position (m)

*r*=_{i}inner radius of the annular channel (m)

*r*=_{o}outer radius of the annular channel (m)

*t*=time (s)

*t*=_{w}waiting time (s)

*T*=temperature (K)

**u**=velocity vector (m/s)

*V*=_{b}volume of the bubble (m

^{3})*y*=distance to the wall (m)

*y*=^{+}non-dimensional wall distance (dimensionless)

*z*=axial position (m)

### Greek Symbols

*α*=void fraction (dimensionless)

*Γ*=mass transfer (kg/m

^{3}s)*ΔT*_{ref}=reference subcooling (K)

*ΔT*_{sub}=subcooling (K)

*λ*=_{l}thermal diffusivity of the liquid (m

^{2}/s)*μ*=_{l}viscosity of the liquid (Pa s)

*μ*=_{t}turbulent viscosity (Pa s)

*ρ*=density (kg/m

^{3})*σ*=surface tension (N/m)

*σ*=_{t}turbulent Schmidt number (dimensionless)

**τ**=viscous stress tensor (N/m

^{2})**τ**=^{t}turbulent stress tensor (N/m

^{2})*ω*=specific rate of turbulence dissipation (1/s)