J. Cent. South Univ. (2016) 23: 1515-1524

DOI: 10.1007/s11771-016-3203-0

Numerical simulation on boiling heat transfer of evaporation cooling in a billet reheating furnace

FENG Ming-jie(冯明杰), WANG En-gang(王恩刚), WANG Hai(王海), LI Yan-dong(李艳东), LIU Bing(刘兵)

Key Laboratory of Electromagnetic Processing of Materials, Ministry of Education

(Northeastern University), Shenyang 110819, China

Central South University Press and Springer-Verlag Berlin Heidelberg 2016

Abstract: The boiling heat transfer of evaporation cooling in a billet reheating furnace was simulated. The results indicate that the bubbles easily aggregate inside of the elbow and upper side of the horizontal regions in the π shaped support tubes. The circulation velocity increasing helps to improve the uniformity of vapor distribution and decrease the difference of vapor volume fraction between upper and down at end of the horizontal sections. With the increase of circulation velocity, the resistance loss and the circulation ratio both increase, but the former will decrease with the increase of work pressure.

Key words: furnace; evaporation cooling; subcooled flow boiling; support tube; two-fluid model

1 Introduction

Evaporation cooling technique is widely used in the pusher type continuous reheating billet furnace due to saving water, avoiding scaling on tube wall, prolonging service life of furnaces and improving the quality of billet reheated. In an evaporative cooling system, the saturation state water will become un-saturation when it flows from the steam drum to the front of the furnace along the downcomer because of heat dissipation of the outer wall. The water nearby the wall will firstly reach saturation state and then generate bubbles because of continuous absorbing heat, but that in mainstream zone is still in un-saturation state after flowing into the furnace. The bubbles slid forward along the wall due to drag force, then detach into water when their diameters gradually become large enough, and condense into water again in the mainstream zone. Thereby, this is a typical 3D sub-cooled flow boiling heat transfer process. About boiling heat transfer phenomenon, as early as in 1934, NUKIYAMA [1] traced out the first boiling heat transfer curve by experiment. After that, boiling heat transfer technique was further developed and more extensively applied, thanks to more people’s research finding [2-4]. In recent decades, lots of numerical simulation on sub-cooling flow boiling heat transfer phenomenon was developed with the development of the computer technique. In the 1990s, KURUL and PODOWSKI [5] established the wall heat flux partition model (RPI model) and two fluid models for numerically analyzing boiling heat transfer. These models were then modified by ISHIMOTO and KAMIJO [6] and BASU et al [7]. In recent years, CHEN [8] and LI et al [9] corrected the bubble departure diameter model on the basis of the UNAL’s mechanism [10] and they confirmed it by their experimental results. DOU [11] and FAN et al [12] simulated the heat transfer characteristics of water sub-cooled flow boiling under different conditions. However, the conditions concerned by above mentioned research are far away from those of the evaporative cooling in furnace. The numerical simulation on the evaporation cooling in furnace has not attracted much attention from researchers. Some present studies only choose one-dimensional two-liquid models and several empirical formulas which are difficult to reflect the actual work conditions.

In this work, an Euler vapor-liquid two-fluid model was used to solve the boiling heat transfer in the π shaped support tubes of an evaporation cooling system in furnace under different conditions. The model consists of a series of sub-models, such as inter-phase heat transfer model, inter-phase mass transfer model, inter-phase momentum transfer model, mean bubble diameter model, bubble departure diameter model, bubble depar ture frequency model and wall heat flux partition model. Furthermore, the influence of pressure loss on saturation temperature and latent heat was also taken into account. Our research objective is to provide reference for the design and operation of evaporative cooling system.

2 Mathematics and physical models

2.1 Geometry model and mesh

Because the temperature and heat transfer conditions are various along the furnace length direction and a part of vertical tubes locate at the out of the furnace, the π shaped support tubes were chosen to be the research objects because of their complex structure, harsh flow and poly-tropic heat transfer conditioning. As shown in Fig. 1, a short horizontal tube at the lower part is connected with the downcomer and a longer horizontal tube at the upper part is connected with the riser tube. The un-saturation water coming from the downcomer firstly flows into the furnace along the lower short horizontal tube, then becomes a mixture of vapor and water because of heat absorption and partial evaporation, and finally flows into the riser tube along the upper long level tube. After separation of steam and water in the drum, vapor is fed into the vapor network and non vaporization water continues the next cycle.

2.2 Mathematics models

In building the models, the basic assumptions were as follows: 1) The influence of external insulation layer on heat transfer can be negligible. 2) The shielding effects of billets and longitudinal tubes can be negligible. 3) The temperature difference between top and bottom in the furnace can be negligible. 4) Flow pattern of water or vapor is incompressible.

Taking the inlet section center of the support tube as origin point and the direction of flow as positive in the x axis direction, a Cartesian coordinate system was established in order to build the mathematics models. According to symmetry, take the 1/2 as the calculation region. Euler vapor–liquid two-fluid models are adopted to describe the two-phase flow, and RNG k-ε turbulence model and steady-state heat transfer model are adopted to describe turbulent flow and wall heat transfer,separately. The details of these models have been described in many previous studies [11-15]. The other models used in this work are described as follows.

2.2.1 Wall heat flux partition model

In sub-cooled flow boiling, heat is transported from the hot wall to the fluid by three different mechanisms. On parts of the wall where no bubbles reside, the mechanism of single phase convective heat transfer is dominant. On parts of the wall where bubbles grow, two mechanisms of evaporation and quenching are dominant, on which grown bubbles flow out and cold liquid flow in and instantly reheated. Thus, the total heat flux of the inner wall of the π shaped support tubes (q_w) will include three parts: liquid phase evaporation heat flux (q_e), quenching heat flux (q_q) and single phase convective heat transfer flux (q_c) between liquid and the wall.

                             (1)

where

                         (2)

                (3)

                      (4)

The nucleation site density n and the bubble detachment diameter d_bw can be given by Eqs. (5) and (6) proposed by LEMMERT and CHAWLA [16] and FRITZ [17], separately.

                               (5)

                   (6)

The bubble detachment frequency f in Eqs. (2) and (3) can be calculated by Eq. (7) given by COLER [18].

                           (7)

Fig. 1 Geometry model of π shaped support tube for calculation of boiling heat transfer and meshes (Unit: mm)

2.2.2 Inter-phase momentum transfer model

The vapor and liquid momentum transfer equation can be described by inter-phase force which mainly involves drag force F_D, lift force FL, turbulence dispersion force F_T, and wall slide force F_W.

                      (8)

The drag force F_D can be calculated by Eq. (9) given by ISHII and ZUBER [19].

              (9)

The lift force FL and the turbulence dispersion force F_T can be calculated by Eqs. (10) and (11), respectively, given by DREW and LAHEY [20].

                (10)

                            (11)

The wall slide force F_W can be calculated by Eq. (12) given by TOMIYAMA [21].

 (12)

The bubble average diameter d_b of Eq. (9) can be calculated by Eq. (13) given by ANGLART and NYLUND [22].

             (13)

2.2.3 Inter-phase mass transfer model

Inter-phase mass transfer model can be given by Eq. (14):

          (14)

The first term on the right side of Eq. (14) is evaporation mass transfer rate of liquid around the wall, and the second one is mass transfer rate in the main flow area. The inter-phase transfer rate in the main flow area depends on the temperature of liquid. Under under- cooling condition, the mass transfer is from vapor to liquid owing to vapor condensation, and contrarily, the mass transfer is from liquid to vapor due to liquid evaporation.

2.2.4 Inter-phase heat transfer model

Inter-phase heat transfer model can be given by Eq. (15):

                         (15)

where inter-phase heat transfer coefficient h_fg can be calculated by Eq. (16) given by RANZ and MARSHALL [23].

                   (16)

2.3 Boundary conditions

1) Inlet boundary conditions

The inlet boundary conditions of liquid phase are determined as follows: taking the average velocity of water and the temperature lower 3 K than the saturation one under work pressure at the entrance as the inlet velocity and inlet temperature, respectively. The liquid volume rate is 1, and the inlet turbulent kinetic energy k and turbulent dissipation rate ε can be calculated by Eqs. (17) and (18), separately.

                                 (17)

                                   (18)

The inlet boundary conditions of vapor phase are as follows the inlet velocity and vapor volume rate are both 0, and the inlet temperature for vapor is the same as that of water.

2) Outlet boundary condition

The outlet boundary condition is set as outflow boundary.

3) Wall boundary conditions

The no-slip boundary condition and standard wall functions are used for the wall boundary conditions of inter wall of the support tubes. For the outer wall, take the mixed heat transfer (Eq. (19)) boundary conditions including radicalization and convection as wall heat boundary conditions in furnace zone, but the convection boundary condition (Eq. (20)) as wall heat boundary conditions in the outer zone.

               (19)

                            (20)

4) Symmetry boundary condition

The boundary condition for symmetry face is determined as symmetry boundary condition, so the normal components of all variables in the face are 0.

3 Research and calculation method

3.1 Research plan

Because the π shaped support tubes studied in this work are made from 20A steel, the data of their thermal properties can be calculated by the thermodynamic calculation module of Procast 2008 according to the composition of 20A steel, then the quadratic function relation between the thermal properties of the π shaped support tubes and temperature can be fitted. The thermodynamic parameters of water can be calculated according to the work pressure, inlet temperature and pressure drop by IAPWS-IF97 which is a thermodynamic industrial calculation module for steam and water. Firstly, the π shaped support tubes were divided into 4 groups according to their installation positions, i.e. the first heating zone, soaking zone, second heating and preheating zone, and separately denoted by positions1-4. The temperatures of the furnace corresponding to these four positions are 1673, 1473, 1273 and 1073 K, respectively. The work pressures of the furnace corresponding to these four positions are 0.4, 0.5, 0.6 and 0.7 MPa, respectively. The circulation velocities of the furnace corresponding to these four positions are 0.5, 1.0, 1.5 and 2 m/s, respectively. Therefore, 64 cases were calculated for each temperature, work pressure, and circulation velocity values.

3.2 Calculation methods

According to Fig. 1, the calculation geometry was created and meshed with hexahedral grid using Gambit 2.3 software. To guarantee the computational efficiency and simulation accuracy, a dense mesh was used in the interface zone between solid and liquid. The mesh file was exported after setting continuum types and boundary types. The sub-models were inputted using user-defined functions based on fluent software. During calculation, the constant temperature flow field had been firstly obtained by closing the energy and boiling model, but opening RNG k-ε turbulent model. Then, take this field as the initialization for the subsequent calculation, and then open the energy and boil model to continue the calculation to convergence. The residuals are 10^-6 for energy and mass conservation functions and 10^-3 for others functions. When the residual of all functions were less the specified value, the flow and temperature field was accepted.

4 Models verification

The applicability of the models and calculation method were verified by experimental results reported in Ref. [24]. In that work, the experimental section was a vertical annular channel having an adiabatic outer tube with an inner diameter R_i=0.0375 m and a heated inner tube with an outer diameter R_o=0.019 m. The heated tube had a length of 1.67 m. The measuring position was located in 1.89 m distance from the inlet. Figure 2 shows the comparison between the simulation results of this work and experimental results of the reference 16 cases under the conditions with operation pressure of 0.1339 MPa, heating flux of 133.6 kW/m², mass flow rate of 478.9 kg/(m²·s) and inlet temperature of 368.85 K. In the figure, the variation trends of simulated values are in agreement with the measured ones. Therefore, the models and calculation method are credible.

Fig. 2 Comparison between simulation results and reference results

5 Results and discussion

5.1 Vapor volume fraction

In a furnace, because almost all of the billet weight is supported by the support tubes, as well as the furnace temperature is quite high, the temperature of the tubes will increase continuously and finally cause the tube wall burning through and the furnace collapsing if the tube wall is cooled inadequately. For this reason, the most important thing for evaporation cooling system is to keep the safe running of evaporation cooling system then to think about producing and offering vapor. By studying the distribution behavior of the gas phase in the tubes, not only the cooling condition of the tubes can be predicted, but the boiling heat transfer process in the tubes can be profoundly understood as well. Figure 3 shows the vapor distribution in the different positions of the support tubes at a work pressure p=0.4 MPa but different circulation velocities. Generally, bubbles easily aggregate in the inside of the bending tubes and upper side of the horizontal tubes, particularly in the horizontal interval of the lower middle of the support tubes and the upper side of the end of the upper horizontal tubes. In addition, the volume fraction of the vapor is also higher in the inside of the vertical tubes located in the outside of the furnace. Figure 3 indicates that the position of the π shaped support tubes and circulation velocity both have great influence upon the vapor distribution. If the circulation velocity keeps constant, with the increase of the furnace temperature around the support tube, the volume fraction of the vapor in the same region increases; however, the whole distribution tendency is nearly unchanged. When the circulation velocity is 1.0 m/s, for the case of the second heating zone, the maximum volume fractions of the vapor in the horizontal interval of the lower middle and in the upper side of the end of the upper horizontal tube are 0.913 and 0.742, respectively (Fig. 3(a)); for the case of the first heating zone, they are 0.931 and 0.947, respectively (Fig. 3(b)). At the same position of the support tube, with increasing circulation velocity, both the volume fraction and distribution of the vapor show an obvious change. When the circulation velocity is 0.5 m/s, at the soaking zone of the support tube, the maximum volume fractions of the vapor in the horizontal interval of the lower middle and in the upper side of the end of the upper horizontal tube are 0.963 and 0.982, respectively (Fig. 3(c)). However, when the circulation velocity increases to 1.5 m/s, they decrease to 0.681 and 0.763 (Fig. 3(d)), respectively. In addition, with increasing circulation velocity, the uniformity of the vapor distribution increases and the difference between the volume fractions of the vapor in the upper part and lower part of the end of the horizontal intervals decreases. This is because with the increase of the circulation velocity, the turbulent kinetic energy of cooling medium increases, as well as the gas and liquid two-phase turbulent mixing intensifies. This in turn weakens the influence of inertia force and gravity. Therefore, the cooling of the tubes can be suitably enhanced by increasing the circulation velocity.

Fig. 3 Vapor volume fraction in π shaped support tubes under different conditions:

5.2 Cross-section average vapor fraction

In the design of evaporation cooling system, the cross-section average vapor fraction is an important design parameter. Figure 4(a) shows the change of cross-section average vapor fraction at Position 1 along flow path under different circulation velocities when the working pressure is 0.7 MPa. Figure 4(b) shows the change of cross-section average vapor fraction along flow path in different support tubes at the workingpressure p=0.7 MPa and the circulation velocity u_L=1.0 m/s. It can be seen that corresponding to different heat transfer processes, each curve can roughly be divided into five sections, i.e. the initial horizontal section, the transition section, the quickly rising section, the slowly rising section, and the rollback section. In the initial horizontal section, the sub-cooling water exhibits a super-cooled state, so the cross-section average vapor fraction is zero and the shape of the corresponding curve is a horizontal line. In this section, the heat transfer shows a liquid single phase convective behavior and the heat transfer capacity is low. Thus, the water temperature will increase quickly for the near wall region but increase very slowly for the main flow region due to the lower heat transfer capacity in the radial direction. In the transition section, the water in the near wall region will quickly reach saturation and form a small number of bubbles, and the water in the main flow region is still at the super-cooled state. However, the degree of the super-cooling will gradually decrease because of the detachment of bubbles which will release a great deal of condensation heat. With the water flowing forward, the generation frequency of the bubbles in the near wall region increases sharply to form nucleate boiling state. At this state, the heat transfer capacity is strongly enhanced because of the transformation of the heat transfer mode from the fluid single phase convective heat transfer to the phase transition heat transfer. This then causes the increase of the cross-section average vapor fraction. When the water in the main flow region is also under saturation states, the bubbles will be mixed in the water and flow forward without condensation. With the increase of the heat transfer, the amount and velocity of the vapor increase gradually, leading to the change in flow pattern. The volume fraction distribution of the vapor along radial direction is different due to the inertia force and gravity force acting on them, causing a change in the heat transfer rate along radial direction of the tube. With increasing boiling time, the amount of the vapor increases, whereas the amount of the water decreases. Thus, the most region of wall is filled with vapor, which causes the decrease of the heat transfer intensity because the most heat needs pass through the vapor film to transfer to interior area. Therefore, in this stage, the cross-section vapor fraction increases slowly, as shown in the figure. Furthermore, a part of the π shaped support tube is out of the furnace, whose wall temperature is higher than the ambient one and causes the heat transfer to the environment. This further causes the small decrease of the cross-section average vapor fraction. In a word, with the increase of the circulation velocity and furnace temperature, the decrease of the temperature gradually decreases. It can also be seen from Fig.4 that with the increase of the circulation velocity and decrease of the furnace temperature, the length of the initial horizontal section gradually increases. This indicates that the water needs to flow a longer distance from the inlet before the bubble generating. Furthermore, the cross- section vapor fraction of the outlet also changes with the change of the position where the support tube locates and the change of the operation parameters. In the practice, if the maximum cross-section vapor fraction ω<0.7, the tubes can be cooled effectively; however, if the maximum cross-section vapor fraction ω>0.82, the water hammer is easy to occur. Therefore, it is of great importance to ensure the ω to be less than 0.8 by adjusting operating parameters. Moreover, if the circulation velocity is small, it is easy to occurrence of the gas-liquid two-phase stratified flow in longer horizontal tubes with the vapor distributing in the upper part of the tube and water distributing in the lower part of the tube. This may causes the damage of the water membrane and lead to overheat of the wall, which may usually cause the crack in the tube, vapor-water corrosion, etc. Therefore, it is also of great importance to control the circulation velocity to avoid the separation of the gas and liquid.

Fig. 4 Change of crosses-section average vapor fraction along flow path under different operation conditions:

5.3 Circulation ratio

The circulation ratio has an important influence on the stable operation of the evaporation cooling system. So, it is an important parameter to be considered in the design and operation of the furnace. According to a previous study, the circulation ratio can be calculated by Eq. (21), compared with the general single support tube, the π shaped support tube requires a higher circulation ratio because the flow distance of water in the π shaped support tube is longer than that in the general single support tube, as well as that it has two dropping vertical tubes. In general, it can assure the demand for the safe operation of the system only when the circulation ratio is more than 10 in a circulation loop.

                             (21)

Figure 5 shows the circulation ratio of the support tubes under different operation conditions. The influence of the circulation velocity on the circulation ratios for four different positions in the support tubes at a work pressure (P) of 0.4 MPa is shown in Fig. 5(a). It can be seen from Fig. 5(a) that the support tubes at high temperature zone in the furnace have a higher circulation ratio when the circulation velocity keeps constant. For a fixed position of the support tube, the circulation ratio increases with increasing circulation velocity. In addition, because the heat transfer is different in the different positions of the support tubes, the amount of the increaseof the circulation ratios with increasing circulation velocity is thus different in the different positions of the support tubes. Under the same condition, the amount of the increase of the circulation ratios at high temperature zone is much larger than that at low temperature zone. For instance, for the preheating zone, when the circulation velocity increases from 0.5 to 2.0 m/s, the circulation ratio increases only from 16.08 to 18.41, but increases sharply from 21.09 to 75.51 for the second heating zone.

Fig. 5 Circulation ratio of support tubes under different operation conditions:

Figure 5(b) shows the influence of the work pressure on the circulation ratios for four different positions when the circulation velocity (u_L) is 1.0 m/s. It can be seen from Fig. 5(b) that the circulation ratio shows a slight increase with increasing work pressure. For the second heating zone, when the work pressure increases from 0.4 to 0.7 MPa, the circulation ratio only increases from 24.84 to 26.11. This is because the primary parameters which can be influenced by the work pressure are saturation temperature and evaporation heat of water. In general, with the increase of the work pressure, the saturation temperature increases, but the evaporation heat decreases. The former means the decrease of the heat transfer flow owing to the decrease in different temperatures between the water and the walls; the later means that a less heat will be needed for evaporating the same amount of water. The changes in saturation temperature and evaporation heat show a contradictory effect on the circulation ratio, leading to a small change of the circulation ratio with work pressure.

5.4 Resistant loss

Figure 6 shows the resistant loss for four different positions under different circulation velocities when the work pressure (p) is 0.4 MPa. It can be seen from Fig. 6 that the resistant loss increases with increasing circulation velocity, but they exhibit a non-quadratic function relation. This is much different from the evolution of the resistant loss for the single phase flow or two-phase flow with a stable gas phase volume fraction. It can also be seen from Fig. 6 that if the circulation velocity keeps constant, the resistant loss at the high temperature zone is larger than that at the low temperature zone. Furthermore, their difference will become larger with the increase of the circulation velocity. Figure 7 shows the change in resistant loss with the work pressure at the same position (Position 3) but different circulation velocities. It can be seen that for allcirculation velocities, the resistant loss decreases with increasing work pressure. At the circulation velocity of 1.5 m/s, the resistant loss decreases from 21.86 to 17.91 kPa when the work pressure increases from 0.4 to 0.7 MPa. In fact, in a natural circulation evaporation cooling system, the driving force for the water flow mainly depends on the fluid density difference between the upper tube and the downcomer and the height difference between the steam drum and the bottom tube. Because the influence of the work pressure on the vapor volume fraction at the outlet of support tubes is not obvious, the work pressure can be enhanced and the installation height of the steam drum can be shortened to reduce the investment and production costs.

Fig. 6 Resistant loss of four different positions support tubes under different circulation velocities as p=0.4 MPa

Fig. 7 Resistant loss of Position 3 support tube under different circulation velocities and work pressures

5.5 Heat load

Some circulation loop in the evaporation cooling systems of the furnace sometimes may emerge flow pulsation, which commonly manifests as periodic fluctuations of circulation flow rate, vapor volume fraction and evaporation site. Some severe flow pulsation even manifests as fluctuations of steam blanketing, of which the whole loop is filled with water; this water may suddenly boil and cause the spraying of water or steam to the drum. The flow pulsation can produce the alternative variation of wall temperature and thermal stress and lead to plastic deformation or fatigue crack. At present, it is difficult to simulate the formation process of the flow pulsation. In general, it can be avoided if the following condition is satisfied.

           (22)

It is necessary to study the influence of some parameters on the heat load in order to accurately predict the stability of circulation loop. Figure 8(a) shows the change in heat load with circulation velocity at four different positions and the work pressure (p) of 0.4 MPa. Figure 8(b) shows the change in heat load with work pressure at the second reheating zone and different circulation velocities. It can be seen from two figures that at different positions, there is a linear relationship between the heat load and the circulation velocity when the work pressure keeps constant. The heat load increases with the increases of the circulation and the slope of each curve is almost same. Furthermore, the heat load increases with the increase of the temperature, at which the support tube locates. For the same position, when the circulation velocity keeps constant, with the increase of the work pressure, the heat load decreases at first and tends to keep unchanged when the work pressure is above 0.6 MPa.In general, it is helpful for the stable operation of circulation loop by fittingly enhancing work pressure.

Fig. 8 Heat load of support tubes under different conditions:

6 Conclusions

1) Bubbles easily aggregate in the inside of the bending tubes and upper side of the horizontal tubes in the π shaped support tubes. With the increase of the circulation velocity, the uniformity of vapor distribution increases. There is a linear relationship between heat load and circulation velocity. Furthermore, the heat load increases with the increase of the temperature, at which the support tube locates.

2) The circulation ratios are affected by support tube position and circulation velocity. When the circulation velocity keeps constant, the support tubes at high temperature zone have a larger circulation ratio. For the same support tube, the circulation ratio increases with the increase of work pressure and circulation velocity.

3) The resistant loss increases and decreases with the increase of circulation velocity and work pressure, separately. If the circulation velocity keeps constant, the resistant loss at high temperature zone is larger than that at low temperature zone. Furthermore, their difference will become larger with the increase of the circulation velocity.

Nomenclature

q_w          Inner wall heat flux, W/m²

q_e          Heat flux due to evaporation, W/m²

q_q          Heat flux due to quenching, W/m²

q_c          Heat flux due to single phase convection, W/m²

n         Nucleation site density, m^-2

d_bw         Bubble detachment diameter, m

h_fg         Heat transfer coefficient for bulk evaporation/ condensation, W/(m²·K)

τ_w         Waiting time before bubble detachment, s

λ_f         Heat transfer coefficient for liquid, W/(m²·K)

A₁        Liquid influence wall area factor

T_w          Wall temperature, K

T_f          Liquid temperature, K

St        Stanton number

c         Liquid specific heat capacity, J/(kg·K)

μ_f         Liquid dynamic viscosity, Pa·s

A₂          Bubble influence wall area factor

f         Bubble detachment frequency, s^-1

θ        Sub-cooling temperature, K

φ_c         Experiment constant

σ        Surface tension, N/m

ρ_f         Liquid density, kg/m³

ρ_g         Steam density, kg/m³

g         Acceleration of gravity, m/s²

F_fg         Gaseous dispersed phase force, N/m³

F_D          Drag force, N/m³

FL        Lift force, N/m³

F_T          Turbulence dispersion force, N/m³

F_W          Wall slide force, N/m³

C_D          Drag coefficient

C_L          Lift coefficient

C_T          Dispersion coefficient

β_g         Gas phase rate

d_b          Bubble average diameter, m

V_f          Liquid phase velocity vector, m/s

V_g          Gas phase velocity vector, m/s

          Laplacian

k_f          Liquid phase turbulent kinetic energy, m²/s²

C_W          Wall slide force coefficient

y_w          Distance to the wall, m

D         Pipe diameter, m

n_w          Wall unit normal vector

d₀, d₁      Reference bubble diameter as sub-cooling temperatures are θ₀ and θ₁, m

A_fg          Interface between two phases, m²

Q_k           Energy exchange between two phases, W

Re_b          Reynolds number for bubble

Pr         Prandtl number

k          Inlet turbulent kinetic energy, m²/s²

ε         Turbulent dissipation rate, m²/s³

u_in          Inlet velocity, m/s

r          Hydraulic radius, m

q₁           Wall absorption heat flux in furnace, W/m²

q₂           Wall dissipate heat flux out of furnace, W/m²

α₁          Heat transfer coefficient between furnace gas and pipe wall, W/(m²·K)

α₂          Heat transfer coefficient between pipe wall and environment, W/(m²·K)

T_w           Wall temperature, K

T₀           Furnace temperature, K

T_f           Environment temperature, K

u_L           Circulation velocity, m/s

p          Work pressure, MPa

ω         Cross-section average vapor volume fraction

L          Flow path length, m

x          Circulation ratio

△p        Resistant loss, Pa

Q          Heat load, W

G          Mass flux, kg/s

△i        Subcooling enthalpy, J/kg

△i′      Enthalpy drop due to pressure decreasing, J/kg

H          Altitude scale, m

ρ′       Apparent density, kg/m³

dF         Area element scale, m

References

[1] NUKIYAMA S. Maximum and minimum values of heat transmitted from a metal to boiling water under atmospheric pressure [J]. Journal of Soc Mechanics Eng, 1934, 37(206): 367-374.

[2] BARTOLOMEJ G G, CHANTURIYA V M. Experimental study of true void fraction when boiling subcooled water in vertical tubes [J]. Thermal Engineering, 1967, 14: 123-128.

[3] MIKIC B, ROHSENOW W M. A new correlation of pool-boiling data including the fact of heating surface characteristics [J]. ASME Journal of Heat Tansfer, 1969, 91(2): 245-250.

[4] KLIMENKO V V, SUDARCHIKOV A M. Investigation of forced flow boiling of nitrogen in a long vertical tube [J]. Cryogenics, 1982, 23(7): 379-385.

[5] KURUL N, PODOWSKI M Z. On the modeling of multidimensional effects in boiling channels [C]// 27th National Heat Transfer Conference. Minneapolis, Minnesota, USA: ANS, 1991: 28-31.

[6] ISHIMOTO J, KAMIJO K. Numerical simulation of cavitating of liquid heliun in vertical channel [J]. Cryogenics, 2003, 43(1): 53-58.

[7] BASU N, WARRIER G R, DHIR V K. Onset of nuclrate boiling and active nucleation site density during subcooled boiling [J]. Journal of Heat Transfer, 2002, 124(4): 717-728.

[8] CHEN R F, LI Y Z, WANG S M. Numerical simulation of subcooled boiling water in vertical concentric annulus under low pressure [J]. Journal of Xi’an Jiaotong Unoversity, 2008, 42(7): 855-859.

[9] LI S W, ZHANG H, JIANG S Y, YU J Y. Numerical simulation research of subcooled boiling water in vertical concentric annulus pipe under low pressure [J]. Nuclear Power Engineering, 2014, 35(1): 46-50.

[10] UNAL H C. Maximum bubble diameter, maximum bubble- growth time and bubble-growth rate during the subcooled nucleate flow boiling of water up to 17.7 MN/m² [J]. International Journal of Heat and Mass Transfer, 1976, 19(5): 643-696.

[11] DOU C C, MAO Y, WANG J, WANG J Y. Numerical simulation of subcooled boiling flow under high pressure and high subcooling condition [J]. Journal of Chemical Industry and Engineering, 2010, 51(3): 580-586.

[12] FAN P, JIA D N, QIU S Z. CFD investigation of suncooled flow boiling model under low pressure [J]. Atomic Energy Science and Technology, 2011, 45(4): 412-420.

[13] BYONG Y, ANDREW S, SIMON L, SONG C H. Prediction of a subcooled boiling flow with advanced two-phase models [J]. Nuclear Engineering and Design 253, 2012: 351-359.

[14] ECKHARD K, ROLAND R. CFD for subcooled flow boiling: Simulation of DEBORA experiments [J]. Nuclear Engineering and Design, 2011, 241: 3851-3866.

[15] FENG M J, WANG E G, HE J C. Numerical simulation on unsteady state heat transfer of work rolls during continuous hot strip rolling [J]. Acta Metallurgica Sinica, 2010, 46(8): 1009-1017.

[16] LEMMERT M, CHAWLA J M. Influence of flow velocity on surface boiling heat transfer coefficient [M]. New York and Washington D. C: Academic Press and Hemisphere, 1977.

[17] FRITZ H C. Maximum bubble diameter, maximum bubble growth time and bubble-growth rate [J]. International Journal of Heat and Mass Transfer, 1935, 36(5): 379-384.

[18] COLER R. A photographic study of pool boiling in the region of the critical heat flux [J]. AIChE Journal, 1960, 6(4): 533-538.

[19] ISHII M, ZUBER N. Drag coefficient and relative velocity in bubbly, droplet or particulate flows [J]. AIChE Journal, 1979, 25(5): 843-855.

[20] DREW D A, LAHEY R T. The virtual mass and lift force on a sphere in rotating and straining inviscid flow [J]. International Journal of Multiphase Flow, 1987, 13(1): 113-121.

[21] TOMIYAMA A. Struggle with computational bubble dynamics [C]// Third International Multiphase Conference. Lyon, France, 1998: 1-18.

[22] ANGLART H, NYLUND O. CFD application to predicition of void distribution in two-phase bubbly flows in rod boundles [J]. Nuclear Engineering and Design, 1996, 163(1/2): 81-98.

[23] RANZ W E, MARSHALL W R. Evaporation from drops [J]. Chemical Engineering Program, 1952, 48(1): 141-146.

[24] LEE T H, PARK G C, LEE D J. Local flow characteristics of subcooled boiling flow of water in a vertical annulus [J]. International Journal of Multiphase Flow, 2002, 28(8): 1351-1368.

(Edited by YANG Hua)

Foundation item: Project(51171041) supported by the National Natural Science Foundation of China

Received date: 2015-04-20; Accepted date: 2015-06-25

Corresponding author: FENG Ming-jie, Associate Professor; Tel: +86-24-83690770; E-mail: fengmj@epm.neu.edu.cn