Numerical Investigation of the Effect of Grids and Turbulence Models on Critical Heat Flux in a Vertical Pipe

Numerical simulation has been widely used in nuclear reactor safety analyses to gain insight into key phenomena. The Critical Heat Flux (CHF) is one of the limiting criteria in the design and operation of nuclear reactors. It is a two-phase flow phenomenon, which rapidly decreases the heat transfer performance at the rod surface. This paper presents a numerical simulation of a steady state flow in a vertical pipe to predict the CHF phenomena. The detailed Computational Fluid Dynamic (CFD) modeling methodology was developed using FLUENT. Eulerian two-phase flow model is used to model the flow and heat transfer phenomena. In order to gain the peak wall temperature accurately and stably, the effect of different turbulence models and wall functions are investigated based on different grids. Results show that O type grid should be used for the simulation of CHF phenomenon. Grids with Y+ larger than 70 are recommended for the CHF simulation because of the acceptable results of all the turbulence models while Grids with Y+ lower than 50 should be avoided. To predict the dry-out position accurately in a fine grid, Realizable k-ε model with standard wall function is recommended. These conclusions have some reference significance to better predict the CHF phenomena of vertical pipe. It can also be expanded to rod bundle of Boiling Water Reactor (BWR) by using same pressure condition.

Introduction

Critical Heat Flux (CHF) is a two-phase flow phenomenon that is characterized by a heat transfer mechanism change which rapidly decreases the efficiency of the heat transfer performance at the heater surface. When it occurs, heated surfaces are no longer wetted by boiling liquid and the vapor phase start to occupy the heat surface. As a result, the energy is directly transferred from the heat surface to vapor. It results in rapid reduction of the heat removal ability and sharp rise of the vapor temperature, as well as the rod surface temperatures which are important for the nuclear safety. During the design and operation of the nuclear reactors, CHF value should be calculated in advance. But experience and thousands of data points have shown the complexity of CHF phenomenon for different conditions. In the past decades, both experiment and numerical simulation are widely used to predict the CHF value during the design and operation process. During the experimental process, the wall temperature is monitored by use of thermocouples. Along with the increasing of heat power, CHF phenomenon is detected when the temperature of one thermocouple has a sharp rise. Similar to the experiment, this sharp rise of the wall temperature is also a signal of CHF in the numerical simulation. So it is of great significance to predict the wall temperature accurately and stably.

Nowadays, experimental method is widely used in the CHF prediction. But most of them are focused on the vertical pipe. Large scale experiments of rod bundle are few reported or open-access around the world. In the reactor design area, CHF look-up table is the main method to determine the CHF value. But both these two methods will give a large safety margin, and then reduce the power level. Compared to experimental measurements and CHF look-up table, numerical simulation has its own advantages. Thermal-hydraulics system or subchannel codes can efficiently characterize bulk flow behavior, while computational fluid dynamics (CFD) analysis may provide relatively accurate results for the velocity and temperature profiles around the rods. Although the subchannel analysis is more widely used for CHF prediction, CFD method can work as a supplement to provide more accurate results. Furthermore, numerical simulation can reduce the cost significantly when we use different geometries, materials, and test conditions. More accurate simulation techniques using these combined analysis tools can be used to optimize experimental designs and improve the test-analysis design cycle for advanced fuel rods.

Grids and turbulence models are found to be very important in the Computational Fluid Dynamic (CFD) simulation. At the moment, turbulence models used in two phase flow are still same with that of single phase flow. Especially in Eulerian two phase model, the liquid and vapor phases are applied with same turbulence model separately. In single phase flow, the results of Chen et al. (2017a) show different turbulence models will deliver to different secondary flow status. This conclusion is more obvious for isotropic and anisotropic turbulence models. However, the conclusion about applicability of these turbulence models is not consistent for the simulation of two phase flow. Standard k-ε model has been chosen by some researchers for its stability on the two phase calculation. Shin and Chang (2009) has used this model to study the effect of mixing vane on CHF enhancement by use of two-fluid model. It is also used by Shirvan (2016) to study the boiling crisis at high pressures. Shirvan's Results show same trend with the Russian CHF data at the condition close to Departure of Nucleate Boiling (DNB). In addition, a new mechanistic model was developed by Mimouni et al. (2016) to predict CHF in boiling flow of water and R12 refrigerant by use of k-ε models in Neptune-CFD code. Mimouni's results are compared with 150 tests and the mean relative error is equal to 8.3%. Apart from standard k-ε model, Shear Stress Transport (SST) k-ω model is also used by some researchers. Krepper et al. (2007, 2013) and Krepper and Rzehak (2011) predicted the subcooled boiling process in vertical heated tubes with water and R12 based on the SST k-ω turbulence model. Subcooled boiling flow in Internal Combustion (IC) engine cooling passages was studied based on SST k-ω turbulence model (Hua et al., 2015). The simulation results showed that present two-fluid model could get an accurate temperature field for cylinder head.

Although different turbulence models are used for the calculation of two phase flow, few researchers has ever compared the effect of different grids and turbulence models. In general, a grid dependency study should always be performed in high quality CFD simulation. Even though it is hard to obtain fully grid independent results in some cases, it is necessary to discuss the effect of key grid parameters. Thakrar et al. (2016) has investigated the nucleate boiling in rectangular geometries at high pressure by use of CFD code STAR-CCM+. Standard k-ε model and Reynolds Stress Transport models were compared to investigate the influence of turbulence. Results show the most mechanistic configurations produced remarkably good agreement with measurements of the area-averaged void. But the grid-dependent wall treatments are still needed to develop. Zhang et al. (2015) has analyzed the effect of grids and turbulence models for subcooled boiling flow in a vertical pipe. Her results show that the turbulence models could be chosen appropriately according to the Y+. This paper shows a relatively overall investigation for the options of grids and turbulence models in the subcooled boiling flow.

However, as we can see above, few researchers have ever shown the sensitivity analysis on turbulence models and wall functions for CHF phenomenon. The grid dependency study is also ignored in the CFD simulation of boiling two-phase flow, especially on the CHF phenomenon. It is still not clear that if the turbulence models have the same performance on CHF phenomena as that of subcooled boiling flow, especially on the value and position of highest wall temperature. So based on Eulerian two-phase model, this paper presents numerical simulations of CHF phenomenon by the use of CFD method to investigate the effect of different grids, turbulence models and wall functions. Lack of open CHF experimental data of rod bundle, especially the wall temperature distribution, we use vertical pipe instead. Considering both the flows in rod bundle and vertical pipe are vertical upward flow, the flow and heat characteristics are similar to each other as long as we use same pressure condition. The type of boiling crisis is dryout which can be identified by the distribution of wall temperature and void fraction. Results are focused on the wall temperature distribution to give a reference for the detector of CHF phenomena.

Mathematical Models

Eulerian two-phase models are used as the basic models to simulate the two phase flow in a vertical pipe. All the interfacial mass, momentum and energy transfer models are based on the interfacial area density model. Furthermore, CHF model will also be shown below.

Conservation Equation

The conservation equation of Eulerian two phase model includes six equations which are mass, momentum and energy equations for two phases separately. These six equations will be solved for six parameters which include pressure, velocities of liquid and vapor, temperatures of liquid and vapor, volume fraction.

(1) Mass equations

$\begin{array}{l}\{\begin{array}{c}\frac{\partial }{\partial t}\left({\alpha }_l{\rho }_l\right)+\nabla ·\left({\alpha }_l{\rho }_l{v}_l\right)=S_l+m_{vl}-m_{lv}\\ \frac{\partial }{\partial t}\left({\alpha }_v{\rho }_v\right)+\nabla ·\left({\alpha }_v{\rho }_v{v}_v\right)=S_v+m_{lv}-m_{vl}& (1)\end{array}\end{array}$

(2) Momentum equations

$\{\begin{array}{c}\frac{\partial }{\partial t}\left({\alpha }_l{\rho }_l{v}_l\right)+\nabla ·\left({\alpha }_l{\rho }_l{v}_l{v}_l\right)=-{\alpha }_l\nabla p+\nabla ·{\tau }_l\\ +{\alpha }_l{\rho }_{l}g+m_{vl}{v}_{vl}-m_{lv}{v}_{lv}+F_l\\ \frac{\partial }{\partial t}\left({\alpha }_v{\rho }_v{v}_v\right)+\nabla ·\left({\alpha }_v{\rho }_v{v}_v{v}_v\right)=-{\alpha }_v\nabla p+\nabla ·{\tau }_v\\ +{\alpha }_v{\rho }_{v}g+m_{lv}{v}_{lv}-m_{vl}{v}_{vl}+F_v& (2)\end{array}$

(3) Energy equations

$\{\begin{array}{c}\frac{\partial }{\partial t}\left({\alpha }_l{\rho }_l{h}_l\right)+\nabla ·\left({\alpha }_l{\rho }_l{v}_l{h}_l\right)={\alpha }_l\frac{\partial p}{\partial t}-\nabla ·q_l\\ +S_l+Q_{vl}+m_{vl}{h}_{vl}-m_{lv}{h}_{lv}\\ \frac{\partial }{\partial t}\left({\alpha }_v{\rho }_v{h}_v\right)+\nabla ·\left({\alpha }_v{\rho }_v{v}_v{h}_v\right)={\alpha }_v\frac{\partial p}{\partial t}-\nabla ·q_v\\ +S_v+Q_{lv}+m_{lv}{h}_{lv}-m_{vl}{h}_{vl}& (3)\end{array}$

Where α,ρ,v,S, p,τ, g, and h refer to volume fraction, density, velocity, source term, pressure, stress-strain tensor, gravity and specify enthalpy. m and Q are the interfacial mass and energy transfers from phase v to phase l. The force F includes five parts, F_D, F_lift, F_wl, F_vm, F_td. They are the drag force, lift force, wall lubrication force, virtual mass force and turbulence dispersion force which will be introduced in detail respectively.

Interfacial Momentum Transfers

The interfacial momentum transfers between liquid and vapor phases are decided by the five forces which are the drag force, lift force, wall lubrication force, virtual mass force and turbulence dispersion force.

For fluid-fluid flow, each secondary phase is assumed to form droplets or bubbles. This assumption gives a method to derivate the drag force which is the most important force of these five parameters. Drag force refers to the force that droplets or bubbles suffer because of the different velocities. As we know, because of viscosity, the general form of the drag force is:

$\begin{array}{ll}F=C_D·\frac{1}{2}\text{ρ}A{V}^2& (4)\end{array}$

When it is used for a particle in the fluid, this form will be changed to

$\begin{array}{ll}F=C_D·\frac{1}{2}A{\rho }_l\left|V_l-V_v\right|\left(V_l-V_v\right)& (5)\end{array}$

Where A is the projected area of the particle, ρ_l is the density of the continuous phase, V_l and V_v are the velocities of different phase. C_D is a coefficient. For the droplets or bubbles, A can be expressed as $A=\frac{1}{4}\pi d^2$, where d is the diameter.

Associated with Re number whose form is

$\begin{array}{ll}Re=\frac{{\rho }_{l}d}{{\mu }_l}\left|V_l-V_v\right|& (6)\end{array}$

We can gain a new form of drag force of the droplets or bubbles which is

$\begin{array}{ll}F=C_D·Re·\frac{1}{8}\pi d{\mu }_l(V_l-V_v)& (7)\end{array}$

So in per unit mixture volume, the drag force can be expressed as

$\begin{array}{l}{F}_D=\frac{F}{\frac{1}{6}\pi d^3}=\frac{C_{D}Re}{24}·\frac{{\rho }_{v}d}{6·\frac{{\rho }_v{d}^2}{18{\mu }_l}}·\frac{\pi d^2}{\frac{1}{6}\pi d^3}(V_l-V_v)\\ =\frac{{\rho }_{v}f}{6{\tau }_l}d{A}_i(V_l-V_v)& (8)\end{array}$

Where $f=\frac{C_{D}Re}{24}$ is always called drag force coefficient which is estimated by Ishii model (Ishii, 1979). ${\tau }_v=\frac{{\rho }_v{d}^2}{18{\mu }_l}$ is called particulate relaxation time. $A_i=\frac{6}{d}$ is the interfacial area concentration which is an important value.

Lift force act on a particle mainly due to velocity gradients in the continuous phase flow field. It is calculated by the function from Drew and Lahey (1993).

$\begin{array}{ll}{F}_{lift}=-C_{lift}{\rho }_l{\alpha }_v(v_l-v_v)\times (\nabla \times v_l)& (9)\end{array}$

In this equation, C_lift is the lift force coefficient provided by Moraga model (Moraga et al., 1999). This model is applicable mainly to the lift force on spherical solid particles, though it can be applied to liquid drops and bubbles.

Wall lubrication force which tends to push the bubbles away from the wall can be expressed as

$\begin{array}{ll}{F}_{wl}=C_{wl}{\rho }_l{\alpha }_v{\displaystyle |(v_l-v_v){|}_z^2}{n}_w& (10)\end{array}$

Where $|(v_l-v_v){|}_z$ is the phase relative velocity component tangential to the wall surface, and n_w is the unit normal pointing away from the wall. The coefficient C_wl is calculated by Antal Model (Antal et al., 1991).

Virtual mass effect should be included when the secondary phase accelerates relative to the primary phase. It is defined as

$\begin{array}{ll}{F}_{vm}=0.5{\alpha }_v{\rho }_l\left(\frac{\partial v_l}{\partial t}+(v_l·\nabla )v_l-\left(\frac{\partial v_v}{\partial t}+(v_v·\nabla )v_v\right)\right)& (11)\end{array}$

The turbulent dispersion force acts as a turbulent diffusion in dispersed flows. It plays an important role in driving the vapor away from the vicinity of the wall to the center of the channel. Bertodano (1991) proposed the following formulation:

$\begin{array}{ll}{F}_{td,l}=-F_{td,v}=C_{TD}{\rho }_l{k}_l\nabla {\alpha }_v& (12)\end{array}$

Where C_TD is a constant, it is equal to 1 by default in the following calculation.

Interfacial Energy Transfers

The interfacial energy transfers includes two parts which are the heat from liquid to vapor phase at the near wall region and the heat transfer between vapor and liquid phases in the subcooled bulk. As the bubbles depart from the wall and move to the subcooled region, there is heat transfer from the bubble to liquid that is defined as

$\begin{array}{ll}{q}_{lt}=h_{sl}(T_{sat}-T_l)& (13)\end{array}$

The heat transfer coefficient can be computed using Ranz and Marshall (1952) model.

The interface to vapor heat transfer is calculated using the constant time scale return to saturation method which was provided by Lavieville et al. (2005). It is assumed that the vapor retains the saturation temperature by rapid evaporation and condensation. The formulation is

$\begin{array}{ll}{q}_{vt}=\frac{{\alpha }_v{\rho }_v{C}_{p,v}}{\delta t}(T_{sat}-T_v)& (14)\end{array}$

In this equation, δt is the time scale.

Interfacial Mass Transfers

Depends on the theory, interfacial mass transfer can also be divided in to two parts: liquid evaporation near the wall and liquid evaporation or vapor condensation in the main stream. The former is calculated on the basis of evaporation heat flux which will be introduced below.

$\begin{array}{ll}{m}_E=\frac{q_E}{h_{fv}+C_{p,l}{T}_{sub}}& (15)\end{array}$

The latter can be calculated directly based on the interfacial energy transfer and latent heat h_fv.

$\begin{array}{ll}m=m_{lt}+m_{vt}=\frac{q_{lt}+q_{vt}}{h_{fv}}& (16)\end{array}$

Wall Boiling Model and CHF Model

Different from the heat transfer of single phase flow on the near wall region, the heat transfer of two-phase flow includes three different types heat transfer. As we can see from the subcooled boiling flow, the energy is transferred directly from the wall to the liquid. Part of this energy named convective heat flux will cause the temperature of the liquid to increase and part which is called evaporative heat flux will generate vapor. Interphase heat transfer will also cause the average liquid temperature to increase; however, the saturated vapor will condense. In addition, there is a quenching heat flux which model the cyclic averaged transient energy transfer related to liquid filling the wall vicinity after bubble detachment. These basic mechanisms are the foundations of the Rensselaer Polytechnic Institute (RPI) models developed by Kurul and Podowski (1991).

To model boiling departing from the nucleate boiling regime, or to model it up to the critical heat flux and post dry-out condition, it is necessary to include the vapor temperature in the solution process. In such cases, some of the energy may be transferred directly from the wall to the vapor. So in total, the wall heat partition will be expressed as

$\begin{array}{ll}{q}_W=f({\alpha }_l)·(q_C+q_Q+q_E)+(1-f({\alpha }_l))·q_V& (17)\end{array}$

The four heat fluxes q_C, q_Q, q_E, q_V is defined as

$\begin{array}{ll}{q}_C=h_C(T_W-T_l)(1-A_b)& (18)\end{array}$

$\begin{array}{ll}{q}_Q=\frac{2k_l}{\sqrt{\pi {\lambda }_{l}T}}(T_W-T_l)& (19)\end{array}$

$\begin{array}{ll}{q}_E=V_d{N}_w{\rho }_v{h}_{fv}f& (20)\end{array}$

$\begin{array}{ll}{q}_V=h_V(T_W-T_V)& (21)\end{array}$

Where T_W, T_l, T_V are the temperature of the wall, liquid and vapor. h_C and h_V are the heat transfer coefficient. k_l is the thermal conductivity. ${\lambda }_l=\frac{k_l}{{\rho }_l{C}_{pl}}$ is the diffusivity. h_fv is the latent heat of evaporation. ρ_v is the vapor density. V_d is the volume of the bubble based on the bubble departure diameter D_w. T is the periodic time. A_b is the area of influence which is covered by nucleating bubbles. N_w is the nucleate site density.

The nucleate site density N_w is usually represented by a correlation based on the wall superheat.

$\begin{array}{ll}{N}_w=C^n{(T_w-T_{sat})}^n& (22)\end{array}$

The empirical parameters C = 210 and n = 1.805 are from Lemmert and Chawla (1977) model.

The bubble departure diameter D_w means the maximum diameters of bubbles when bubbles leave the wall. The empirical correlation is given by Tolubinski and Kostanchuk (1970).

$\begin{array}{ll}{D}_w=min\left(0.0014,0.0006e^{-\frac{T_{sub}}{45}}\right)& (23)\end{array}$

The area of influence A_b is based on the departure diameter and the nucleate site density.

$\begin{array}{ll}{A}_b=min\left(1,K\frac{N_w\pi {D_w}^2}{4}\right)& (24)\end{array}$

The coefficient K is given by Del Valle and Kenning (1985).

$\begin{array}{ll}K=4.8e^{\left(-\frac{J_{a_{sub}}}{80}\right)}& (25)\end{array}$

$\begin{array}{ll}{J}_{a_{sub}}=\frac{{\rho }_l{C}_{pl}{T}_{sub}}{{\rho }_v{h}_{fv}}& (26)\end{array}$

The last parameter is the frequency of bubble departure $f=\frac{1}{T}$. In general, the bubble in the wall suffers three forces, which are buoyancy force, surface tension and drag force. When the heat flux is high and bubbles are large, the surface tension could be ignored. So according to the equilibrium of buoyancy force and drag force which was described in paper (Cole, 1960).

$\begin{array}{ll}g({\rho }_l-{\rho }_v)\frac{\pi D_w^3}{6}=F_D\frac{\pi D_w^2}{4}\frac{{\rho }_l{u}^2}{2}& (27)\end{array}$

Associated with u = fD_w and F_D = 1, f can be expressed as.

$\begin{array}{ll}f=\sqrt{\frac{4g({\rho }_l-{\rho }_v)}{3{\rho }_l{D}_w}}& (28)\end{array}$

In the model, the critical heat flux is defined to occur when the volume fraction of vapor reach given values which is provided by Tentner et al. (2006).

$f\left({\alpha }_v\right)=1-f\left({\alpha }_l\right)=\{\begin{array}{c}\begin{array}{cc}0& {\alpha }_v<{\alpha }_{v,1}\end{array}\\ \begin{array}{c}\frac{1}{2}\left(1-\cos \left(\pi \frac{{\alpha }_v-{\alpha }_{v,1}}{{\alpha }_{v,2}-{\alpha }_{v,1}}\right)\right){\alpha }_{v,1}\le {\alpha }_v\le {\alpha }_{v,2}\end{array}\\ \begin{array}{cc}1& {\alpha }_v>{\alpha }_{v,2}\end{array}& (29)\end{array}$

Where α_{v, 1} = 0.9 and α_{v, 2} = 0.95.

Numerical Scheme

The numerical analysis is performed using the commercial CFD code FLUENT 16.0 (Ansys Fluent, 2015) which is part of ANSYS WORKBENCH 16.0. The modules we use are based on the finite volume method, and uses RANS (Reynolds averaged Navier-Stokes) equations for the solution of mass, momentum, energy and turbulence model.

Geometry and Boundary Condition

CHF experiment data in an upward flow vertical pipe from Becker (1983) are used to validate the simulation method. Although there are much data in the Becker's experiment, one case with 7Mpa pressure which is condition of boiling water reactor (BWR) is investigated about the effect of grids and turbulence models. The facility has a 7m long heated pipe and the inner diameter is 10mm. Detail information about boundary condition are shown in Figure 1.

FIGURE 1

Figure 1. Geometry and boundary condition.

Grids, Turbulence Models and Wall Functions

It has been tested that the grid type has limit influence on the temperature distribution in the single-phase flow (Chen et al., 2017b). But in two-phase flow, this conclusion is doubtful for the complexity of wall boiling status. So the effect of different grid type should be tested first before the investigation of turbulence models. As we can see from Figure 2, O type grid, arbitrary Hexaprism type grid and Triprism grid are calculated under the same boundary condition. Results are shown in section Effect of Grid Type.

FIGURE 2

Figure 2. O type grid, arbitrary Hexaprism grid and Triprism grid.

All the turbulence models can be divided into high Reynolds and low Reynolds turbulence models according to its application. The standard k-ε model (Launder and Spalding, 1972) is typical high Reynolds turbulence model. It should be used for fully developed turbulence and isotropic flow. Because of its reduced computational complexity, it has been widely used for decades of years. An important weak point is that the standard k-ε model cannot deal with rotational flow or flow in curved pipe. To solve these problems, RNG (Re-Normalization Group) k-ε model (Yakhot and Orszag, 1986) and Realizable k-ε model (Shih et al., 1995) are developed in 1986 and 1995 separately. These two models expanded the application of k-ε model to a larger field. The secondary flow status of shear flow and rotational flow could be predicted to a certain extent. The k-ω models are low Reynolds turbulence model which have a larger application. The standard k-ω model (Wilcox, 1998) is more suitable for wall limited flow and free shear flow. In SST (Shear Stress Transport) k-ω model (Menter, 1994), the effect of turbulent shear force is included when solve the turbulent viscosity. It is also worth to note that SST k-ω model is same as standard k-ω model in near wall region and similar to standard k-ε model in the fully turbulent region.

The turbulence model should be used with suitable wall function or resolve the boundary layer with a fine-enough mesh without any wall functions. As we know, the near-wall region can be subdivided into three sublayers, i.e., viscous sublayer, buffer layer and fully turbulent region (Pope, 2000). For k-ε model with standard wall function, it uses an empirical correlation to solve the equations in viscous sublayer and buffer layer. It also requires the first node should be located in fully turbulent region where viscosity is of less effect. So according to reference (Pope, 2000), Y+ larger than 30 is recommended for high Reynolds turbulence models to ensure accuracy and stability. However, the enhanced wall function maybe helpful to use high Reynolds turbulence model in the region where Y+ is smaller than 30. The results sometimes show good agreement with experiment data. For k-ω model, it is essential to ensure the viscous sublayer to be covered with several cells to maintain the near wall Y+ to 1. However, k-ω model could be Y+-independent if works with enhanced wall function. This is also how the k-ω model is used in FLUENT.

In order to validate the correlation of Y+, grids and turbulence models, five sets of grids are used under different turbulence models and wall functions. The detailed information about grid is listed in Table 1. The computational matrix of grid, turbulence model and wall function is shown in Table 2. StdWF, EnhWF, NeqWF, and ScaWF represent Standard Wall Function (Launder and Spalding, 1974), enhanced wall function, Non-Equilibrium Wall Function (Kim et al., 1999) and Scalable Wall Function, respectively. The detail information about range of application is shown in Table 3.

TABLE 1

Table 1. Grids information.

TABLE 2

Table 2. Computational matrixes of grids, turbulence model and wall function.

TABLE 3

Table 3. Wall functions.

Not all the possible combinations are calculated in this paper for the numerous computational time and resource. Because of the same basic mechanism for the three k-ε models, the results of each can be used to the others. So the simulation of ICEM1 to ICEM5 on Realizable k-ε models and enhanced wall function are calculated to investigate the effect of different grids. This work is also done for the standard k-ω models on ICEM3~ICEM5. The effect of turbulence models are compared based on ICEM3 and ICEM4. We believe the conclusion can be extended to the other grids. For the wall function part, different wall functions can only be used on different k-ε models. So in the k-ω models, no wall functions are calculated except the enhanced wall function.

Results and Discussion

The effect of grid type, grids number, turbulence and wall functions are discussed below. Results are compared among the highest temperature values and dry-out location of the inner surface of the pipe. The temperature distribution is the key safe-relative parameter in the CHF experiment and should be paid more attention to discuss.

Effect of Grid Type

As has been validated, the grid type has limit influence on the temperature distribution in single phase flow. But for two phase flow, a little difference is found among three grid types. The grid types are shown in Figure 2, and the numbers of nodes are 1344000, 1298500, and 1228500 separately. Figure 3 shows the temperature distribution of the inner surface of the pipe. Results show all the three grids have the same dry-out position while the highest temperature values are different. More detail information can be found in Figure 4 which is the radial temperature distribution at 4.5 m. The distribution of Hexaprism grid shows a little decentered phenomenon which is not in accordance with theory. Furthermore, the distribution in Triprism grid shows a zigzag in the center of the channel. Though this phenomenon can be avoided through refining the grid, it will cost more computing resource and time. So to predict the dry-out position and highest temperature value accurately, O type or axis-symmetric grid should be chosen for the calculation.

FIGURE 3

Figure 3. Temperature distribution of the inner surface of the pipe.

FIGURE 4

Figure 4. Radius distribution of temperature at the 4.5 m height. (A) O type grid. (B) Hexaprism grid. (C) Triprism grid.

Effect of Grids Number on Turbulence Models

Grid independent calculation is always an important requirement in the numerical simulation. It is necessary to estimate the quality of gird through the calculations of different grids' level. However, it is not easy to gain fully grid-independent results for all the simulation, especially for two phase flow. But It has been tested that all the k-ε models and k-ω models with enhanced wall function can be used for a large scale of Y+ values from one to hundreds in single phase flow. Zhang et al. (2015) had also investigated the effect of grids and Y+ values in the subcooled boiling flow. It is concluded that results of k-ε models match well with experiment data when near-wall grid Y+ is larger than 11.25 and enhanced wall function cannot cope with the grids with wall adjacent grid Y+ near to 1. However, results of CHF calculation show a different conclusion which has a more strict limitation for Y+ value.

Figure 5 show the comparison of the surface temperature under different grids on standard k-ε model and enhanced wall function. Results show little difference is detected among different grids except the highest temperature. The simulation results show a good match with experimental data. The dry-out positions of simulation results and experimental data are also in accord with each other. What's more, the difference of highest temperature decreased along with the refinement of the grids. So it may have the trend to gain grid-independent results.

FIGURE 5

Figure 5. Surface temperatures under different grids on standard k-ε model and enhanced wall function.

The surface temperature results including ICEM4 are shown in Figure 6, and it is calculated based on RNG k-ε model and enhanced wall function. As we can see from Figure 6, the dry-out position of ICEM4 falls behind that of the other three grids. The highest value of surface temperature of ICEM4 is also much lower than that of the other three grids and the experimental data. So the results start to be worse along with the refinement of grid and the decrease of the Y+ value.

FIGURE 6

Figure 6. Surface temperatures under different grids on RNG k-ε model and enhanced wall function.

For a further validation, the comparison results including the ICEM5 are shown in Figure 7 which is calculated based on Realizable k-ε model and enhanced wall function. The temperature distribution of ICEM 4 shows same trend as that in Figure 6. The dry-out position is found behind the others and the highest temperature value is lowest. We can also see from Figure 6 that when compared with good performance of the other grids, the temperature distribution of ICEM5 shows a large fluctuation in 0.7-2.5m. What's more, the highest temperature is lower than the other four. After applying ICEM3, ICEM4 and ICEM5 with Standard k-ω model, we also found same phenomena in the same position. The results are shown in Figure 8. The Y+ value of ICEM 5 is shown in Figure 9. After comparing this two figure, we found that this unreasonable phenomena mainly appear near the wall when the Y+ value is lower than 50. The reason is that the transient behavior of bubbles has a great influence on the space averaged property when the near wall cells are too small.

FIGURE 7

Figure 7. Surface temperatures under different grids on Realizable k-ε model and enhanced wall function.

FIGURE 8

Figure 8. Surface temperatures under different grids on Standard k-ω model.

FIGURE 9

Figure 9. Y+ value of ICEM5 on Realizable k-ε model and enhanced wall function.

Effect of Turbulence Models

Different turbulence models are applied to the CHF calculation with enhanced wall function based on ICEM3 and ICEM4. The results are shown in Figures 10, 11. Some of the results show good agreement with experimental data while the others not.

FIGURE 10

Figure 10. Surface temperatures under different turbulence models on ICEM3.

FIGURE 11

Figure 11. Surface temperatures under different turbulence models on ICEM4.

Figure 10 shows the inner surface temperature distribution of the pipe which is calculated under different turbulence models on ICEM3. Little difference is detected between Standard k-ω model and SST k-ω model. Only the highest value of temperature shows a difference about 2-5K which can be ignored. The results of Standard k-ε model and RNG k-ε model are nearly the same while both of these two are slightly lower than the results of Realizable k-ε model. The highest value of temperature of k-ω models is 30K higher than that of k-ε models which is a considerable difference. However, when compared with experimental data, all the simulation results are acceptable. The predicted dry-out position and temperature distribution are in accord with experimental data.

When ICEM4 is employed to the CHF calculation, large differences are found among the different turbulence models. The results are shown in Figure 11. The two k-ω models still show similar results in dry-out position and temperature distribution. Furthermore, these two results are also in agreement with experimental data except the post-dryout position. For k-ε models, the dry-out positions are predicted behind the experimental data and the highest temperatures are also lower than the measurement. It is worth to say that all the above calculations are under enhanced wall function. The effect of wall functions will be discussed in next section.

Effect of Wall Functions on k-ε Models

Heat transfer and flow characters in near wall regions are important to the calculation of temperature and velocity in two-phase flow. The wall functions and wall boiling model are also influenced by each other. Because the wall functions use different empirical correlations to solve the equations in viscous sublayer and buffer layer, different wall functions will have different performance. In FLUENT, k-ω models are low Reynolds models and set to be used with wall functions unless the grid Y+ doesn't meet the requirements of the low Reynolds number modeling. If the requirements are met, then FLUENT doesn't apply wall function but resolves the near wall region. So there is no need to discuss the effect of different wall functions when use k-ω models.

Based on ICEM3, the CHF calculation results of k-ε models under different wall functions are shown in Figure 12 separately. All these three pictures show same trend for different wall functions. The temperature of non-equilibrium wall functions is 20K higher than that of the others while the enhanced wall function has the lowest temperature. The temperature distributions of standard wall function and scalable wall function are the same. This is determined by the correlation of scalable wall function.

FIGURE 12

Figure 12. Surface temperatures under different wall functions on (A) Standard k-ε model, (B) RNG k-ε model, (C) Realizable k-ε model.

The purpose of scalable wall functions is to force the usage of the log law in conjunction with the standard wall functions approach. This is achieved by adding a limiter in the Y+ calculations Y + is calculated by

$\begin{array}{ll}Y+=\mathrm{MAX}\left(Y+,11.25\right)& (30)\end{array}$

So when the Y+ is larger than 11.25, Y+ is the same with original value, and then the scalable wall function is the same as standard wall function. Otherwise, Y+ is equal to 11.25 which bring different results. For the cases in this paper, all the Y+ values are larger than 11.25, so there is no difference for the results between these two wall functions. When compared with experimental data, the standard wall function and scalable wall function show the best performance while results of the other two are also acceptable.

Based on ICEM4, the comparison is only investigated between standard wall function and enhanced wall function on two k-ε models. The results are shown in Figure 13. As mentioned above, the predicted dry-out positions of enhanced wall function are behind the measurement when ICEM4 is employed to the simulation. Among these four cases, only Realizable k-ε model with standard wall function can give an accurate dry-out position and acceptable temperature distribution.

FIGURE 13

Figure 13. Surface temperatures under different wall functions on ICEM4.

Conclusion

CHF phenomenon in a vertical pipe is investigated by use of CFD simulation. Simulations are calculated among different grids type, grids number, different turbulence models and wall functions to detect the effect of these parameters. Several conclusions are obtained after comparing the simulation results and experimental data. The conclusion about the effect of grid and turbulence models can help someone reduce their workload on mesh validation and set up solve algorithm quickly.

(1) type or axis-symmetric grid is recommended for the CFD simulation of vertical pipe or the other geometries.

(2) Grids with Y+ larger than 70 are recommended for the CHF simulation while Grids with Y+ lower than 50 should be avoided for the simulation of CHF phenomenon in the vertical pipe.

(3) The temperature distributions of k-ω models are 20K higher than that of k-ε models in all the cases. But they are all acceptable when comparing with experimental data.

(4) The predicted dry-out position of enhanced wall function is always behind that of the others if used with a fine grid. To predict the dry-out position accurately in a refine grid, Realizable k-ε model with standard wall function is recommended.

Author Contributions

XD is the main author of the paper, he finished the research and write. ZZ is the supervisor and give the guide for the other authors. DL is the co-author who works for the comparison of the experiment and simulation. ZT is the associated supervisor and give the reference too. GC works with the other authors for the discussion part.

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

This work is supported by the Research on Key Technology of Numerical Reactor Engineering (J121217001), and the authors also appreciate deeply the suggestions from Dr. Tenglong Cong.

References

Ansys Fluent (2015). Ansys Fluent 16.0 Documentation.

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, 635–652. doi: 10.1016/0301-9322(91)90029-3

CrossRef Full Text | Google Scholar

Becker, K. M. (1983). An Experimental Investigation of Post Dryout Heat Transfer. Department of Nuclear Reactor Engineering, Royal Institute of Technology, KTH-NEL-33.

Google Scholar

Bertodano, L. D. (1991). Turbulent Bubbly Flow in a Triangular Duct. Ph.D. Thesis. Rensselaer Polytechnic Institute, Troy, New York, NY.

Chen, G. L., Zhang, Z. J., Tian, Z. F., Dong, X. M., and Ju, H. R. (2017a). Design of a CFD scheme using multiple RANS models for PWR. Ann. Nucl. Energy. 102, 349–358. doi: 10.1016/j.anucene.2016.12.030

CrossRef Full Text | Google Scholar

Chen, G. L., Zhang, Z. J., Tian, Z. F., Li, L., and Dong, X. M. (2017b). Challenge analysis and schemes design for the CFD simulation of PWR. Sci. Technol. Nucl. Installations. 2017:5695809. doi: 10.1155/2017/5695809

CrossRef Full Text | Google Scholar

Cole, R. (1960). A photographic study of pool boiling in the region of the critical heat flux. AIChE J. 6, 533–542. doi: 10.1002/aic.690060405

CrossRef Full Text | Google Scholar

Del Valle, V. H., and Kenning, D. B. R. (1985). Subcooled flow boiling at high heat flux. Int. J. Heat Mass Transf. 28, 1907–1920. doi: 10.1016/0017-9310(85)90213-3

CrossRef Full Text | Google Scholar

Drew, D. A., and Lahey, R. T. (1993). In Particulate Two-Phase Flow. Boston, MA: Butterworth-Heinemann.

Google Scholar

Hua, S. Y., Huang, R. H., and Zhou, P. (2015). Numerical investigation of two-phase flow characteristics of subcooled boiling in IC engine cooling passages using a new 3d two-fluid model. Appl. Therm. Eng. 90, 648–663. doi: 10.1016/j.applthermaleng.2015.07.037

CrossRef Full Text | Google Scholar

Ishii, M. (1979). “Two-fluid model for two-phase flow,” in Proceedings of the 2nd International Workshop on Two-phase Flow Fundamentals. Troy, NY: RPI.

Google Scholar

Kim, S. E., Choudhury, D., and Patel, B. (1999). “Computations of complex turbulent flows using the commercial code fluent,” in Modeling Complex Turbulent Flows, eds M. Salas, J. Hefner, and I. Sakell (Dordrecht: Springer), 259–276.

Google Scholar

Krepper, E., Koncar, B., and Egorov, Y. (2007). CFD modelling of subcooled boiling—concept, validation and application to fuel assembly design. Nucl. Eng. Des. 237, 716–731. doi: 10.1016/j.nucengdes.2006.10.023

CrossRef Full Text | Google Scholar

Krepper, E., and Rzehak, R. (2011). CFD for subcooled flow boiling: simulation of debora experiments. Nucl. Eng. Des. 241, 3851–3866. doi: 10.1016/j.nucengdes.2011.07.003

CrossRef Full Text | Google Scholar

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, 330–346. doi: 10.1016/j.nucengdes.2012.11.010

CrossRef Full Text | Google Scholar

Kurul, N., and Podowski, M. Z. (1991). “On the modeling of multidimensional effects in boiling channels,” in Proceedings of the 27th National Heat Transfer Conference (Minneapolis, MN).

Google Scholar

Launder, B. E., and Spalding, D. B. (1972). Lectures in Mathematical Models of Turbulence. London: Academic Press.

Launder, B. E., and Spalding, D. B. (1974). The numerical computation of turbulent flows. Comput. Methods Appl. Mech. Eng. 3, 269–289. doi: 10.1016/0045-7825(74)90029-2

CrossRef Full Text

Lavieville, J., Quemerais, E., Mimouni, S., Boucker, M., and Mechitoua, N. (2005). NEPTUNE CFD V1.0 Theory Manual. EDF.

Lemmert, M., and Chawla, L. M. (1977). “Influence of flow velocity on surface boiling heat transfer coefficient,” in Heat Transfer Boiling, eds E. Hahne and U. Grigull (New York, NY: Academic Press and Hemisphere).

Menter, F. R. (1994). Two-equation eddy-viscosity turbulence models for engineering applications. AIAA J. 32, 1598–1605. doi: 10.2514/3.12149

CrossRef Full Text | Google Scholar

Mimouni, S., Baudry, C., Guingo, M., Lavieville, J., Merigoux, N., and Mechitoua, N. (2016). Computational multi-fluid dynamics predictions of critical heat flux in boiling flow. Nucl. Eng. Des. 299, 28–36. doi: 10.1016/j.nucengdes.2015.07.017

CrossRef Full Text | Google Scholar

Moraga, F. J., Bonetto, R. T., and Lahey, R. T. (1999). Lateral forces on spheres in turbulent uniform shear flow. Int. J. Multiphase Flow 25, 1321–1372. doi: 10.1016/S0301-9322(99)00045-2

CrossRef Full Text | Google Scholar

Pope, S. B. (2000). Turbulent Flows. Cambridge, UK: Cambridge University Press.

Google Scholar

Ranz, W. E., and Marshall, W. R. (1952). Evaporation from drops, chemical engineering. Progress 48, 141–146.

Google Scholar

Shih, T. H., Liou, W. W., Shabbir, A., Yang, Z., and Zhu, J. (1995). A new k-ε Eddy-viscosity model for high reynolds number turbulent flows-model development and validation. Comput. Fluids 24, 227–238. doi: 10.1016/0045-7930(94)00032-T

CrossRef Full Text | Google Scholar

Shin, B. S., and Chang, S. H. (2009). CHF experiment and CFD analysis in a 2 × 3 rod bundle with mixing vane. Nucl. Eng. Des. 239, 899–912. doi: 10.1016/j.nucengdes.2009.01.011

CrossRef Full Text | Google Scholar

Shirvan, K. (2016). Numerical investigation of the boiling crisis for helical cruciform-shaped rods at high pressures. Int. J. Multiphase Flow 83, 51–61. doi: 10.1016/j.ijmultiphaseflow.2016.03.014

CrossRef Full Text | Google Scholar

Tentner, A., Lo, S., and Kozlov, V. (2006). “Advances in computational fluid dynamics modeling of two-phase flow in a boiling water reactor fuel assembly,” in International Conference on Nuclear Engineering (Miami, FL).

Google Scholar

Thakrar, R., Murallidharan, J., and Walker, S. P. (2016). CFD Investigation of nucleate boiling in non-circular geometries at high pressure. Nucl. Eng. Des. 312, 410–421. doi: 10.1016/j.nucengdes.2016.08.020

CrossRef Full Text | Google Scholar

Tolubinski, V. I., and Kostanchuk, D. M. (1970). “Vapor bubbles growth rate and heat transfer intensity at subcooled water boiling,” in Proceedings of the 4th International Heat Transfer Conference (Paris).

Google Scholar

Wilcox, D. C. (1998). Turbulence Modeling for CFD. La Canada: DCW Industries, Inc.

Google Scholar

Yakhot, V., and Orszag, S. A. (1986). Renormalization group analysis of turbulence i basic theory. J. Sci. Comput. 1, 1–51. doi: 10.1007/BF01061452

CrossRef Full Text | Google Scholar

Zhang, R., Cong, T. L., Tian, W. X., Qiu, S. Z., and Su, G. H. (2015). CFD analysis on subcooled boiling phenomena in PWR coolant channel. Prog. Nucl. Energy 81, 254–263. doi: 10.1016/j.pnucene.2015.02.005

CrossRef Full Text | Google Scholar