## 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 [14]. 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].

Continuity equation
$∂∂t(αkρk)+∇·(αkρkuk)=Γk$
(1)
Momentum equation
$∂∂t(αkukρk)+∇.(αkρkukuk)=αk∇P+αkρkg+αk∇·(τk+τkt)+Γkuk+Mk$
(2)
Energy equation
$∂∂t(αkρkHk)+∇.(αkρkukHk)=αk∂P∂t−∇.qk″+Qk+ΓkHk$
(3)

where $αk$, $ρk$, $uk$, $qk″$, $Hk$, $Γk$, Mk, and $Qk$ represent volume fraction, density, velocity, heat flux, specific enthalpy, mass transfer, momentum transfer, and energy transfer, respectively. P indicates pressure, $τk$ and $τ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

1. Evaporative heat flux

2. Liquid phase convective flux

3. Quenching heat flux

An illustrative drawing of the wall boiling model is shown in Fig. 1. The following expression gives the total heat flux:
$q″total=(q″C+q″Q+q″E)$
(4)

where $q″total$ is total heat flux, $q″C$ is convective heat flux, $q″Q$ is quenching heat flux, and $q″E$ is evaporative heat flux.

Heat transferred to the liquid is given by convective heat flux as
$q″C=hC(Tw−Tl)(1−Ab)$
(5)

where hC 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.

Liquid replaces the detached bubble, and heat is transferred to the liquid via transient conduction up until the formation of the new bubble; this process is referred to as quenching. Quenching depends on the waiting period (tw); 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]:
$q″Q=2klπλltw(Tw−Tl)Ab$
(6)

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

The fraction of heat is used to create bubbles over the heater surface; this is taken into account in evaporative flux and is given by the following expression [11]:
$q″E=VbNwρvhfgf$
(7)

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:

The fraction of the heater surface affected by the evaporation process depends on nucleation site density, bubble departure diameter, and area influence factor as given by [11]
$Ab=min(1, KNwπDbw24)$
(8)

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.

Tolubinsky and Kostanchuk [12] gave a correlation for bubble departure diameter based on liquid subcooling, and it goes as follows:
$Dbw=min(0.0014,drefe−ΔTsubΔTref)$
(9)

where $dref$ = 0.0006 m, $ΔTsub=Tsat−Tl$, and $Δ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.

Bubble departure frequency calculation mainly depends upon the bubble departure diameter and interfacial forces. Many models are available for determining bubble frequency. Most widely used model for bubble departure frequency is given by Cole [13]. The expression goes as follows:
$f= [4g(ρl−ρν)3ρlDbw]1/2$
(10)

#### Nucleation Site Density.

Nucleation site activation depends on wall surface roughness, local fluid parameters, wall superheat, and wettability of the wall surface with the liquid. It is widely accepted that the nucleation site density increases with an increase in wall superheat. Lemmert and Chawla [14] proposed a correlation for nucleation site density based on wall superheat as following:
$Nw=(210·(Tw−Tsat))1.805$
(11)

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

### 2.4 Momentum Closures

#### Drag Force.

Drag force takes into the resistance on the bubble when it flows in the fluid. It is a function of relative velocity (drift) between vapor and liquid. It also depends on the size and shape of the bubble and fluid through which it is moving. The following equation gives drag force:
$FD=−CDμlAiReb8Db(uv−ul)$
(12)
where $Reb=ρl|uv−ul|Db/μl$, $CD$ is the drag coefficient, $μl$ is the viscosity of the liquid, $Db$ is the bubble diameter, $Ai$ is the interfacial area concentration, $uv$ and $ul$ represents the velocity of bubble and liquid, respectively. Ishii and Zuber [15] developed the widely used model for drag coefficient for two-phase systems, considering three different shape regimes (sphere, ellipse, and cap), the expression goes as follows [2]:
$CD=max(24(1+0.1Reb0.75)Reb, min(23Db(g|ρl−ρv|σ),83))$
(13)

where $σ$ is surface tension.

#### Turbulent Dispersion Force.

Turbulent fluctuations in the continuous phase have a profound impact on the redistribution of the dispersed phase. This effect is considered in the form of a turbulent dispersion force. This force is proportional to the gradient of the void fraction. The two models, Lopez de-Bertodano et al. [16] and Burns et al. [17], are widely used. Lopez de-Bertodano et al. [16] specifies the turbulent dispersion force as follows:
$FTD=−CTDρlk∇α$
(14)
Burns et al. [17] derived the expression by time-averaging the interfacial drag term, the expression goes as follows:
$FTD=−3CDμt4Dbσt(uv−ul)∇ α1−α$
(15)

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

#### Lift Force.

Lift force accounts for the lateral motion of bubbles, which is mainly due to the interaction between the bubble and the shear stress in the liquid. Drew and Lahey [18] derived an expression for lift force on a sphere as
$FL=−CLρlα(uv−ul)×rot(ul)$
(16)

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

Moraga et al. [19] found that the lateral lift force correlates with Reynolds number based on relative velocity, $Reb$ and average shear, $Reω$. They noted that the lift force is a balance between two competing components, inviscid lift forces, and the vortex shedding-induced lateral forces. The coefficient is given by the following expression:
$CL=0.0767 RebReω≤6000−(0.12−0.2e−RebReω3.6×10−5)eRebReω3×10−7 6000
(17)

where $Reω=ρl|∇×ul|Db2/μl.$

Klausner et al. [20] extended the Saffman [21] model, which was developed for spherical solid particles in a simple shear field. Primarily, this model was developed for low bubble Reynolds number. According to the Saffman force, the particles drift toward the side where the relative velocity magnitude is high. They correlated the lift coefficient in terms of $Reb$ and $Reω$. The expression goes as following:
$CL=6.46×f(Reb,Reω) Reb≤40 6.46×0.0524(Reω2)12 40
(18)
where
$f(Reb,Reω)=(1−0.3314(Reω2Reb)0.5)e−0.1Reb+0.3314(Reω2Reb)0.5$
Legendre and Magnaudet [22] derived numerically by solving the Navier–Stokes equations. They found the lift coefficient to depend strongly on both Reynolds number and shear rate at low Reynolds numbers. However, they found that for moderate to high Reynolds numbers, these dependencies are not considerable. They proposed a lift model based on the bubble Reynolds number and nondimensional shear rate (Sr). The following expression gives the lift coefficient:
$CL=(CL,lowRe)2+(CL,highRe)2$
(19)
where
$CL,lowRe=6π2(RebSr)−0.52.55(1+0.2RebSr)1.5$

$CL,highRe=121+16Reb−11+29Reb−1$

$Sr=ReωReb$

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.

The following expression gives the lift coefficient:
$CL=min[0.288 tanh(0.121 Reb),f( Eod)] for Eod<4f( Eod) for 4≤Eod≤10.7−0.27 for Eod>10.7$
(20)
where
$Eod=(ρl−ρv)gDd2σ$

$Dd=Db(1+0.16Eo0.757)1/3$

$Eo=(ρl−ρv)gDb2σ$

$f(Eod)=0.00105 Eod3−0.0159 Eod2−0.0204 Eod+0.474$

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.

Antal et al. [24] proposed a model by assuming asymmetric drainage of liquid around the spherical bubble in the vicinity of the wall; it goes as follows:
$FWL=CWLρLα(uv−ul)2n$
(21)

$CWL=max(Cw1+Cw2Db2y,0)$
(22)

In this study, the coefficients provided by Frank et al. [25] are used (Cw1 = −0.01 and Cw2 = 0.05).

Tomiyama [26] proposed a new model by replacing Cw1 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 Cw2 to be zero. The final form of the equation is as follows:
$CWL=Db2(1y2−1(D−y)2).e−0.933Eo+0.179, 1
(23)
Frank et al. [25] further extended Tomiyama [26] model by removing the dependence on the diameter of the pipe. The equation goes as follows:
$CWL=max(0,16.81−y10Dby(y10Db)0.7).e−0.933Eo+0.179 , 1
(24)

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−ω$ 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−ω$ 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−ω$ 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).

According to Sato et al. [35], the additional term is given by the following expression:
$μBIT=CμρlDbα|uv−ul|$
(25)

Another way to consider the bubble effects is by adding source terms to k and $ω$ 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 $Δ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].

As mentioned in Sec. 5, the base case model is used where lift and wall lubrication models are not considered. In Secs. 6.1 and 6.2, CFD predictions with different models are also compared with the base model to understand the significance of the models. Moreover, for better understanding, the deviation of CFD predictions from the experimental data is also plotted. The following expression gives deviation (in %):
$Deviation (%)=αpre−αexpαexp×100$

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:

1. 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%.

2. 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.

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

4. 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.

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

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

## Nomenclature

• Ai =

interfacial area (1/m)

•
• Ab =

fraction of the heater surface area influenced by the evaporation process (dimensionless)

•
• CD =

drag coefficient (dimensionless)

•
• CL =

lift coefficient (dimensionless)

•
• CL,lowRe =

lift coefficient at low Reynolds numbers (dimensionless)

•
• CL,highRe =

lift coefficient at high Reynolds numbers (dimensionless)

•
• CTD =

turbulent dispersion coefficient (dimensionless)

•
• CWL =

wall lubrication coefficient (dimensionless)

•
• Cw1 =

coefficient (dimensionless)

•
• Cw2 =

coefficient (dimensionless)

•
• Cμ =

shear-induced turbulence coefficient (dimensionless)

•
• Db =

bubble diameter (m)

•
• Dd =

maximum horizontal bubble diameter (m)

•
• Dbw =

bubble departure diameter (m)

•
• dref =

reference diameter of bubble (m)

•
• FD=

drag force (N/m3)

•
• FL =

lift force (N/m3)

•
• FTD=

turbulent dispersion force (N/m3)

•
• FWL=

wall lubrication force (N/m3)

•
• f =

bubble departure frequency (1/s)

•
• g =

acceleration due to gravity (m/s2)

•
• H =

specific enthalpy (J/kg)

•
• hC =

liquid heat transfer coefficient (W/m2 K)

•
• hfg =

latent heat of vaporization (J/kg)

•
• K =

area influence factor (dimensionless)

•
• k =

turbulent kinetic energy (m2/s2)

•
• kl =

thermal conductivity of liquid (W/m K)

•
• M =

interfacial momentum transfer (N/m3)

•
• n =

unit normal vector perpendicular to the wall (dimensionless)

•
• Nw =

nucleation site density (1/m2)

•
• P =

pressure (Pa)

•
• Q =

energy transfer (W/m3)

•
• q" =

heat flux (W/m2)

•
• r =

•
• ri =

inner radius of the annular channel (m)

•
• ro =

outer radius of the annular channel (m)

•
• t =

time (s)

•
• tw =

waiting time (s)

•
• T =

temperature (K)

•
• u =

velocity vector (m/s)

•
• Vb =

volume of the bubble (m3)

•
• y =

distance to the wall (m)

•
• y+ =

non-dimensional wall distance (dimensionless)

•
• z =

axial position (m)

### Greek Symbols

Greek Symbols

• α =

void fraction (dimensionless)

•
• Γ =

mass transfer (kg/m3 s)

•
• ΔTref =

reference subcooling (K)

•
• ΔTsub =

subcooling (K)

•
• λl =

thermal diffusivity of the liquid (m2/s)

•
• μl =

viscosity of the liquid (Pa s)

•
• μt =

turbulent viscosity (Pa s)

•
• ρ =

density (kg/m3)

•
• σ =

surface tension (N/m)

•
• σt =

turbulent Schmidt number (dimensionless)

•
• τ =

viscous stress tensor (N/m2)

•
• τt =

turbulent stress tensor (N/m2)

•
• ω =

specific rate of turbulence dissipation (1/s)

### Subscripts

Subscripts

• b =

bubble

•
• BIT =

bubble induced turbulence

•
• C =

convective

•
• E =

evaporative

•
• exp =

experimental value

•
• k =

kth phase

•
• l =

liquid

•
• pre =

predicted value

•
• Q =

quenching

•
• sat =

saturation

•
• tot =

total

•
• v =

vapor

•
• w =

wall

### Non-Dimensional Numbers

Non-Dimensional Numbers

• Eo =

Eötvös number

•
• Eod =

modified Eötvös number

•
• Reb =

bubble Reynolds number

•
• Reω =

vorticity Reynolds number

•
• Sr =

non-dimensional shear rate

### Acronyms and Abbreviations

Acronyms and Abbreviations

• CFD =

computational fluid dynamics

•
• DNB =

departure from nucleate boiling

•
• SST =

shear stress transport

## References

References
1.
Krepper
,
E.
, and
Rzehak
,
R.
,
2011
, “
CFD for Subcooled Flow Boiling: Simulation of DEBORA Experiments
,”
Nucl. Eng. Des.
,
241
(
9
), pp.
3851
3866
.10.1016/j.nucengdes.2011.07.003
2.
Rzehak
,
R.
, and
Krepper
,
E.
,
2013
, “
CFD for Subcooled Flow Boiling: Parametric Variations
,”
Sci. Technol. Nucl. Install.
,
2013
, pp.
1
22
.10.1155/2013/687494
3.
Colombo
,
M.
, and
Fairweather
,
M.
,
2016
, “
Accuracy of Eulerian–Eulerian, Two-Fluid CFD Boiling Models of Subcooled Boiling Flows
,”
Int. J. Heat Mass Transfer
,
103
, pp.
28
44
.10.1016/j.ijheatmasstransfer.2016.06.098
4.
Murallidharan
,
J. S.
,
,
B. V. S. S. S.
,
Patnaik
,
B. S. V.
,
Hewitt
,
G. F.
, and
,
V.
,
2016
, “
CFD Investigation and Assessment of Wall Heat Flux Partitioning Model for the Prediction of High Pressure Subcooled Flow Boiling
,”
Int. J. Heat Mass Transfer
,
103
, pp.
211
230
.10.1016/j.ijheatmasstransfer.2016.06.050
5.
Colombo
,
M.
, and
Fairweather
,
M.
,
2015
, “
Prediction of Bubble Departure in Forced Convection Boiling: A Mechanistic Model
,”
Int. J. Heat Mass Transfer
,
85
, pp.
135
146
.10.1016/j.ijheatmasstransfer.2015.01.103
6.
Yeoh, G. H., and Tu, J. Y., 2006, “
Two-Fluid and Population Balance Models for Subcooled Boiling Flow
,”
Appl. Math. Model.
,
30
(
11
), pp.
1370
1391
.10.1016/j.apm.2006.03.010
7.
Krepper
,
E.
,
Rzehak
,
R.
,
Lifante
,
C.
, and
Frank
,
T.
,
2013
, “
CFD for Subcooled Flow Boiling: Coupling Wall Boiling and Population Balance Models
,”
Nucl. Eng. Des.
,
255
, pp.
330
346
.10.1016/j.nucengdes.2012.11.010
8.
Li
,
H.
,
Vasquez
,
S. A.
,
Punekar
,
H.
, and
Muralikrishnan
,
R.
,
2011
, “
Prediction of Boiling and Critical Heat Flux Using an Eulerian Multiphase Boiling Model
,”
ASME Paper No. IMECE2011-65539
.10.1115/IMECE2011-65539
9.
Yun
,
B. J.
,
Splawski
,
A.
,
Lo
,
S.
, and
Song
,
C. H.
,
2012
, “
Prediction of a Subcooled Boiling Flow With Advanced Two-Phase Flow Models
,”
Nucl. Eng. Des.
,
253
, pp.
351
359
.10.1016/j.nucengdes.2011.08.067
10.
Drew
,
D. A.
,
1992
, “
Analytical Modeling of Multiphase Flows
,”
Boiling Heat Transfer: Modern Developments and Advances
,
R.T.
Lahey
, ed.,
Elsevier
,
Amsterdam, The Netherlands
, pp.
31
84
.
11.
Kurul
,
N.
, and
Podowski
,
M. Z.
,
1990
, “
Multidimensional Effects in Forced Convection Subcooled Boiling
,”
Proceedings of the Ninth International Heat Transfer Conference
, Vol.
2
,
Jerusalem, Israel
, Aug. 19–24, pp.
21
26
.
12.
Tolubinsky
,
V. I.
, and
Kostanchuk
,
D. M.
,
1970
, “
Vapour Bubbles Growth Rate and Heat Transfer Intensity at Subcooled Water Boiling
,”
Proceedings of Fourth International Heat Transfer Conference
, Vol.
5
,
Paris, France
, Aug. 31–Sept. 5, Paper No. B-2.8, pp.
1
11
.
13.
Cole
,
R.
,
1960
, “
A Photographic Study of Pool Boiling in the Region of the Critical Heat Flux
,”
AIChE J.
,
6
(
4
), pp.
533
538
.10.1002/aic.690060405
14.
Lemmert
,
M.
, and
Chawla
,
J.
,
1977
, “
Influence of Flow Velocity on Surface Boiling Heat Transfer Coefficient
,”
Heat Transfer and Boiling
,
E.
Hahne
and
U.
Grigull
, eds.,
,
New York
, pp.
237
247
.
15.
Ishii
,
M.
, and
Zuber
,
N.
,
1979
, “
Drag Coefficient and Relative Velocity in Bubbly, Droplet or Particulate Flows
,”
AIChE J.
,
25
(
5
), pp.
843
855
.10.1002/aic.690250513
16.
Lopez de Bertodano
,
M.
,
Lahey
,
R. T.
, and
Jones
,
O. C.
,
1994
, “
Turbulent Bubbly Two-Phase Flow Data in a Triangular Duct
,”
Nucl. Eng. Des.
,
146
(
1–3
), pp.
43
52
.10.1016/0029-5493(94)90319-0
17.
Burns
,
A. D.
,
Frank
,
T.
,
Hamill
,
I.
, and
Shi
,
J.-M. M.
,
2004
, “
The Favre Averaged Drag Model for Turbulent Dispersion in Eulerian Multi-Phase Flows
,”
Proceedings of the Fifth International Conference on Multiphase Flow (ICMF-2004)
,
Yokohama, Japan
, May 30–June 4, Paper No. 392.https://www.researchgate.net/publication/261761053_The_Favre_Averaged_Drag_Model_for_Turbulent_Dispersion_in_Eulerian_Multi-Phase_Flows
18.
Drew
,
D. A.
, and
Lahey
,
R. T.
,
1987
, “
The Virtual Mass and Lift Force on a Sphere in Rotating and Straining Inviscid Flow
,”
Int. J. Multiphase Flow
,
13
(
1
), pp.
113
121
.10.1016/0301-9322(87)90011-5
19.
Moraga
,
F. J.
,
Bonetto
,
F. J.
, and
Lahey
,
R. T.
,
1999
, “
Lateral Forces on Spheres in Turbulent Uniform Shear Flow
,”
Int. J. Multiphase Flow
,
25
(
6–7
), pp.
1321
1372
.10.1016/S0301-9322(99)00045-2
20.
Klausner
,
J. F.
,
Mei
,
R.
,
Bernhard
,
D. M.
, and
Zeng
,
L. Z.
,
1993
, “
Vapor Bubble Departure in Forced Convection Boiling
,”
Int. J. Heat Mass Transfer
,
36
(
3
), pp.
651
662
.10.1016/0017-9310(93)80041-R
21.
Saffman
,
P. G.
,
1965
, “
The Lift on a Small Sphere in a Slow Shear Flow
,”
J. Fluid Mech.
,
22
(
2
), pp.
385
400
.10.1017/S0022112065000824
22.
Legendre
,
D.
, and
Magnaudet
,
J.
,
1998
, “
The Lift Force on a Spherical Bubble in a Viscous Linear Shear Flow
,”
J. Fluid Mech.
,
368
, pp.
81
126
.10.1017/S0022112098001621
23.
Tomiyama
,
A.
,
Tamai
,
H.
,
Zun
,
I.
, and
Hosokawa
,
S.
,
2002
, “
Transverse Migration of Single Bubbles in Simple Shear Flows
,”
Chem. Eng. Sci.
,
57
(
11
), pp.
1849
1858
.10.1016/S0009-2509(02)00085-4
24.
Antal
,
S. P.
,
Lahey
,
R. T.
, and
Flaherty
,
J. E.
,
1991
, “
Analysis of Phase Distribution in Fully Developed Laminar Bubbly Two-Phase Flow
,”
Int. J. Multiphase Flow
,
17
(
5
), pp.
635
652
.
25.
Frank
,
T.
,
Shi
,
J. M.
, and
Burns
,
A. D.
,
2004
, “
Validation of Eulerian Multiphase Flow Models for Nuclear Safety Application
,”
Proceedings of the Third International Symposium on Two-Phase Flow Modelling and Experimentation
,
Pisa, Italy
, Sept. 22–24, p.
9
.
26.
Tomiyama
,
A.
,
1998
, “
Struggle With Computational Bubble Dynamics
,”
Multiphase Sci. Technol.
,
10
(
4
), pp.
369
405
.
27.
Hosokawa
,
S.
,
Tomiyama
,
A.
,
Misaki
,
S.
, and
,
T.
,
2002
, “
Lateral Migration of Single Bubbles Due to the Presence of Wall
,”
ASME Paper No. FEDSM2002-31148
.10.1115/FEDSM2002-31148
28.
Takemura
,
F.
,
Takagi
,
S.
,
Magnaudet
,
J.
, and
Matsumoto
,
Y.
,
2002
, “
Drag and Lift Forces on a Bubble Rising Near a Vertical Wall in a Viscous Liquid
,”
J. Fluid Mech.
,
461
, pp.
277
300
.10.1017/S0022112002008388
29.
Hassan
,
Y. A.
,
,
C. E.
, and
Yoo
,
J. S.
,
2014
, “
Measurement of Subcooled Flow Boiling Using Particle Tracking Velocimetry and Infrared Thermographic Technique
,”
Nucl. Eng. Des.
,
268
, pp.
185
190
.10.1016/j.nucengdes.2013.04.044
30.
Vaidheeswaran
,
A.
,
,
D.
,
Guilbert
,
P.
,
Buchanan
,
J. R.
, and
de Bertodano
,
M. L.
,
2016
, “
New Two-Fluid Model Near-Wall Averaging and Consistent Matching for Turbulent Bubbly Flows
,”
ASME J. Fluids Eng.
,
139
(
1
), p.
011302
.10.1115/1.4034327
31.
Lubchenko
,
N.
,
Magolan
,
B.
,
Sugrue
,
R.
, and
Baglietto
,
E.
,
2018
, “
A More Fundamental Wall Lubrication Force From Turbulent Dispersion Regularization for Multiphase CFD Applications
,”
Int. J. Multiphase Flow
,
98
, pp.
36
44
.10.1016/j.ijmultiphaseflow.2017.09.003
32.
Zhang
,
R.
,
Cong
,
T.
,
Tian
,
W.
,
Qiu
,
S.
, and
Su
,
G.
,
2015
, “
Effects of Turbulence Models on Forced Convection Subcooled Boiling in Vertical Pipe
,”
Ann. Nucl. Energy
,
80
, pp.
293
302
.10.1016/j.anucene.2015.01.039
33.
Rzehak
,
R.
, and
Krepper
,
E.
,
2013
, “
CFD Modeling of Bubble-Induced Turbulence
,”
Int. J. Multiphase Flow
,
55
, pp.
138
155
.10.1016/j.ijmultiphaseflow.2013.04.007
34.
Mali
,
C. R.
,
Vinod
,
V.
, and
Patwardhan
,
A. W.
,
2017
, “
Comparison of Phase Interaction Models for High Pressure Subcooled Boiling Flow in Long Vertical Tubes
,”
Nucl. Eng. Des.
,
324
, pp.
337
359
.10.1016/j.nucengdes.2017.09.010
35.
Sato
,
Y.
,
,
M.
, and
Sekoguchi
,
K.
,
1981
, “
Momentum and Heat Transfer in Two-Phase Bubble Flow—I: Theory
,”
Int. J. Multiphase Flow
,
7
(
2
), pp.
167
177
.10.1016/0301-9322(81)90003-3
36.
Garnier
,
J.
,
Manon
,
E.
, and
Cubizolles
,
G.
,
2001
, “
Local Measurements on Flow Boiling of Refrigerant 12 in a Vertical Tube
,”
Multiphase Sci. Technol.
,
13
(
1–2
), p.
111
.10.1615/MultScienTechn.v13.i1-2.10
37.
Roy
,
R. P.
,
Kang
,
S.
,
Zarate
,
J. A.
, and
Laporta
,
A.
,
2002
, “
Turbulent Subcooled Boiling Flow—Experiments and Simulations
,”
ASME J. Heat Transfer
,
124
(
1
), pp.
73
93
.10.1115/1.1418698
38.
Bartolomei
,
G. G.
, and
Chanturiya
,
V. M.
,
1967
, “
Experimental Study of True Void Fraction When Boiling Subcooled Water in Vertical Tubes
,”
Therm. Eng.
,
14
(
2
), pp.
123
128
.
39.
Kumar
,
M.
,
Moharana
,
A.
,
Nayak
,
A. K.
, and
Joshi
,
J. B.
,
2018
, “
CFD Simulation of Boiling Flows Inside Fuel Rod Bundle of a Natural Circulation BWR During SBO
,”
Nucl. Eng. Des.
,
338
, pp.
300
329
.10.1016/j.nucengdes.2018.08.011
40.
Anglart
,
H.
,
Nylund
,
O.
,
Kurul
,
N.
, and
Podowski
,
M. Z.
,
1997
, “
CFD Prediction of Flow and Phase Distribution in Fuel Assemblies With Spacers
,”
Nucl. Eng. Des.
,
177
(
1–3
), pp.
215
228
.
41.
Ranz
,
W. E.
, and
Marshall
,
W. R.
,
1952
, “
Evaporation From Drops 1
,”
Chem. Eng. Prog.
,
48
(
3
), pp.
141
146
.