- Email: [email protected]
Numerical simulation of heat transfer deterioration phenomenon in supercritical water through vertical tube
Annals of Nuclear Energy 37 (2010) 1272–1280
Contents lists available at ScienceDirect
Annals of Nuclear Energy journal homepage: www.elsevier.com/locate/anucene
Numerical simulation of heat transfer deterioration phenomenon in supercritical water through vertical tube Q.L. Wen a,b,*, H.Y. Gu a a b
School of Nuclear Science and Engineering, Shanghai Jiao Tong University, Shanghai 200240, China Reactor Engineering Research Division, Nuclear Power Institute of China, Chengdu 610041, China
a r t i c l e
i n f o
Article history: Received 23 March 2010 Received in revised form 15 May 2010 Accepted 31 May 2010 Available online 1 July 2010 Keywords: Supercritical water Heat transfer deterioration Turbulence model Numerical analysis
a b s t r a c t In this study, a numerical investigation of heat transfer deterioration (HTD) in supercritical water flowing through vertical tube is performed by using six low-Reynolds number turbulence models. All low-Reynolds models can be extended to reproduce the effect of buoyancy force on heat transfer and show the occurrence of localized HTD. However, most k–e models seriously over-predict the deterioration and do not reproduce the subsequent recovery of heat transfer. The V2F and SST models perform better than other models in predicting the onset of deterioration due to strong buoyancy force. The SST model is able to quantitatively reproduce the two heat transfer deterioration phenomena with low mass flux which have been found in the present study. Ó 2010 Elsevier Ltd. All rights reserved.
1. Introduction Supercritical water-cooled reactor (SCWR) has been regarded as an innovative reactor and selected as one of candidates of Generation IV reactor systems based on the criteria of economics, safety, and sustainable production of efficient nuclear power. During the last fifteen years, research activities are ongoing worldwide for developing of the advanced nuclear power plants as SCWR (Oka et al., 2007; Schulenberg et al., 2008). One of the main features of supercritical water is the strong variation of its thermal–physical properties in vicinity of the pseudo-critical temperature. Fig. 1 shows the thermal–physical properties of water at 23.5 MPa. The large variation of thermal–physical properties with temperature leads to strong coupling between the flow and the temperature fields, and creates significant effect by buoyancy force or thermal-induced flow acceleration which affect both the mean flow fields and the turbulence structures. These conditions lead to unusual heat transfer characteristics such as HTD phenomenon (Jackson, 2001). Such phenomenon is generating gradually, but can result in a substantial increase of the fuel cladding temperature. Therefore, an in-depth understanding of the phenomenon is essential to specify cladding temperature limits for the design of the SCWR core. Studies on the thermal–hydraulic behavior of supercritical fluids have been performed since the 1950s. Existing experimental * Corresponding author at: School of Nuclear Science and Engineering, Shanghai Jiao Tong University, Shanghai 200240, China. Tel.: +86 21 34204917. E-mail address: [email protected] (Q.L. Wen). 0306-4549/$ - see front matter Ó 2010 Elsevier Ltd. All rights reserved. doi:10.1016/j.anucene.2010.05.022
and theoretical studies on HTD at supercritical pressure conditions were reviewed and published by several authors (Cheng and Schulenburg, 2001; Pioro and Duffey, 2005). As mentioned in the previous paragraph, the increase in the wall temperature under heat transfer deterioration condition is much slower and smoother than that in the case of the boiling crisis in classical Light Water Reactors (LWRs). Therefore, it is difficult to define the onset of HTD. In the literature, different definitions were used, most of which are based on the ratio of the heat transfer coefficient to a reference value (Yamagata et al., 1972; Styrikovich et al., 1967; Jackson and Hall, 1979; Ogata and Sato, 1972; Petuhkov and Kurganov, 1983; McEligot and Jackson, 2004). A detailed derivation in predictions between different correlations is obtained as indicated by Cheng et al. (2007). Due to the measurement technique limitation and costs associated with experimental study, detailed information on experimental flow field and turbulence structure is rare for water at supercritical pressure (Licht, 2008). With the recent developed technologies in computational methods, numerical investigations using Computational Fluid Dynamics (CFD) codes have achieved a significant attention with the purpose to obtain prediction and to provide a better understanding of the heat transfer mechanism. In recent years, some successful cases have been achieved. The main difficulties in the numerical analysis are related to the turbulence modeling at supercritical pressures. During the last few years, the low-Reynolds number type eddy viscosity turbulence models are used to investigate the HTD. Koshizuka et al. (1995) used JL turbulence model (Jones and Launder, 1972) to study the HTD phenomenon in supercritical
Q.L. Wen, H.Y. Gu / Annals of Nuclear Energy 37 (2010) 1272–1280
1273
Nomenclature Bo* buoyancy parameter Cl, Ce1, Ce2 turbulence models’ constants specific heat transfer (J kg1 K1) Cp D pipe diameter (m)/damping function defined in Table 4 E damping function defined in Table 4 f1, f2, fl damping functions Gr Grashof number g gravity accelerations (m s2) G mass flux (kg m2 s1) h heat transfer coefficient (W m2 K1)/enthalpy (J kg1) m mass flow rate (kg s1) k turbulent kinetic energy (m2 s2) Nu Nusselt number P pressure (MPa) Pk turbulent shear production (W m3) Pr Prandtl number r, x radial and axial coordinates (m) q heat flux (W m2) Re Reynolds number
T U, V y y+, y*
temperature (K) velocity components in the x-, r-directions (m s1) distances in coordinates (m) non-dimensional distance from pipe wall
Greek letters rate of dissipation of k (m2 s3) k thermal conductivity (W m1 K1) l viscosity (kg m1 s1) s shear stress (N m2) q density (kg m3)
e
Subscripts b bulk t turbulent w wall pc pseudo-critical
Fig. 1. Thermal–physical properties for water at 23.5 MPa.
water. The buoyancy effect at low flow rate and flow acceleration effect at high flow rate on heat transfer were analyzed. Calculation results qualitatively agreed with the experimental data of Yamagata et al. (1972). He et al. (2004) applied six low-Reynolds number turbulence models to simulate the turbulent mixed convection heat transfer experiments with carbon dioxide at supercritical pressure. The results showed that all models were able to simulate the influence of buoyancy to heat transfer, but those models have significantly different performances to predict onset and magnitude of this effect. Sharabi et al. (2008) used five low-Reynolds turbulence models, i.e. ABID model (Abid, 1993), LB (Lam and Bremhorst, 1981), AKN model (Abe et al., 1994), CHC model (Chang et al., 1995), YS (Yang and Shih, 1993) and RNG model with the two-layer approach (Yakhot and Orszag, 1986), implementing by FLUENT, for prediction of three-dimensional turbulent heat transfer to carbon dioxide at supercritical pressure flowing upward through heated square and triangular channels as presented by Kim et al. (2005). The results showed that the low-Reynolds turbulence models were able to reproduce the trends of heat transfer deterioration due to buoyancy influence which may lead to relatively large variations of wall temperatures. Among the low-Reynolds number models, the YS model (Yang and Shih, 1993) was able to qualitatively reproduce
the second localized deterioration phenomenon observed in the experimental data downstream of the location of first one. On the other hand, the RNG model with the two-layer approach was unable to capture even the expected qualitative trends. It was found that the flow acceleration due to density non-uniformity has insignificantly affected on heat transfer deterioration for the cases with strong buoyancy influence as considered in the study. He et al. (2008a) made simulations for vertical upward and downward flow of carbon dioxide in a heated tube at a pressure just above the critical value using the V2F turbulence model (Lien et al., 1998) and AKN model (Abe et al., 1994). The results showed that both models can qualitatively show the general trends of HTD due to strong buoyancy influence in upward flow, but significant discrepancies in detailed wall temperature distribution are found. The AKN model is the best of the two for the conditions considered. In the study, they found that buoyancy influence on heat transfer was almost entirely due to the shear production effect caused by the distortion of the mean flow as a result of the influence of buoyancy. Palko and Anglart (2008) used the low-Reynolds k–x turbulence model to make simulations of heat transfer deterioration. The results showed for the low mass flux of coolant, the influence of buoyancy forces on the heat transfer in heated pipes is signifi-
1274
Q.L. Wen, H.Y. Gu / Annals of Nuclear Energy 37 (2010) 1272–1280
cant. For the high mass flux, buoyancy influence could be neglected and there are clearly other mechanisms causing the decrease in heat transfer at high coolant mass flux. Recently, He et al. (2008b) assessed the performance of seven low-Reynolds number turbulence models in responding to the modifications of the turbulence field due to influences of flow acceleration and buoyancy to carbon dioxide at supercritical pressure by comparing model predictions with DNS results by Bae et al. (2005). The results showed that the low-Reynolds number turbulence models whose damping functions are based on variables readily responding to buoyancy and flow acceleration (LS, YS, AKN), significantly over-predict flow laminarization and therefore, the HTD. However, turbulence models whose damping functions are based on variables not responding to buoyancy/flow acceleration (CH (Chien, 1982) and MK (Myoung and Kasagi, 1990) models) reproduce closely the variations of wall temperature exhibited in the DNS in flow laminarization cases due to some canceling effects. The V2F model produces the best predictions among all the turbulence models tested. This paper is devoted to computational simulations of experiments on HTD phenomena with a well-validated commercial CFD software FLUENT, using k–e type and k–x type low-Reynolds turbulence models. The purpose of the work presented here is to numerically study the interesting behavior exhibited in experimental results using detailed information on the flow and thermal fields. 2. Experimental data used The experimental results considered in this study (Pis’menny et al., 2006; Shitsman, 1963) were obtained with water at supercritical pressure flowing in a uniformly heated vertical tube. The experimental conditions which have been simulated in this study are summarized in Table 1. Run1 and run3 are in the HTD region due to strong buoyancy influence. The experimental condition of run2 is the same as that of run1, except for the flow direction. More experimental details can be found in Pis’menny et al. (2006) and Shitsman (1963).
V-momentum:
1 @ @ ðr quv Þ þ ðr qv 2 Þ r @x @r @p 1 @ @ v @u @ @v lv r le þ2 r le þ 2 e2 ¼ þ @x r @x @r @x @r @r r ð3Þ where le is the effective viscosity defined by the turbulent viscosity defined as
lt ¼ qC l fl
k
le ¼ l þ lt and lt is
2
ð4Þ
e
in which fl is a damping function to account for near-wall effects and C l is a constant. Energy:
1 @ @ ðr quhÞ þ ðrqv hÞ r @x @r
¼
1 @ l lT @h r þ r @x Pr rT @x @ l lT @h þ r þ @r Pr rT @r
ð5Þ
where Pr is molecular Prandtl number and rT is the turbulent Prandtl number, the constant 0.9 is used in this study. Six turbulence models have been selected, namely, five k–e type low-Reynolds turbulence models due to Abe et al. (1994) (AKN), Yang and Shih (1993), Chang et al. (1995) (CHC) and Abid (1993) (AB), the k–e–v2–f model of Behnia et al. (1998) (V2F) and one k– x type low-Reynolds turbulence model, i.e., shear stress transport model (SST) developed by Menter (1993), the SST model combines the robustness of the k–e model for the bulk flow with the low-Re treatment of the boundary layer using the modified Wilcox k–x model. The general form of k–e models can be expressed as follows: Turbulent kinetic energy:
3. Governing equations and turbulence models
@ 1 @ ðqukÞ þ ðr qv kÞ @x r @r @ l @k 1 @ l @k lþ t þ r lþ t þ Pk þ Gk qðe DÞ ¼ @x rk @x r @r rk @r
In the FLUENT code, the flow is assumed to be 2D axisymmetry and the governing equations of the continuity, momentum and energy for such a flow written in cylindrical coordinates are as follows: Continuity:
ð6Þ where the production term due to the mean velocity gradient is:
1 @ @ ðrquÞ þ ðr qv Þ r @x @r
¼0
ð1Þ
" (
P k ¼ lt 2
2 2 2 ) 2 # @u @v v @u @ v þ þ þ þ @x @r @x @r r
ð7Þ
The gravitational production term is computed by:
U-momentum:
1 @ @ ðr qu2 Þ þ ðr qv uÞ r @x @r @p 1 @ @u @ @u @ v 2 r le þ r le þ ¼ þ qg þ @x r @x @x @r @r @x
Gk ¼ q0 u0 g x ¼ ch bqgx ð2Þ
k
e
lT
@u @ v @T @u 2 @T þ þ 2lT qk @r @x @r @x 3 @x
ð8Þ
where g x ¼ g for upward flow, g x ¼ g for downward flow and ch ¼ 0:3.
Table 1 Experimental conditions.
a
Cases
P (MPa)
G (kg/m2 s)
qw (kW/m2)
D (mm)
Dir.
Boa 107
Authors
Run1 Run2 Run3
23.5 23.5 23.5
509 509 430
390 390 319.87
6.28 6.28 8
Up Down Up
5.39–5.75 – 9.14–4.57
Pis’menny et al. (2006) Pis’menny et al. (2006) Shitsman (1963)
2
Bo ¼ Gr=ðRe3:425 Pr0:8 Þ, Gr ¼ gbD4 qw =kv .
Q.L. Wen, H.Y. Gu / Annals of Nuclear Energy 37 (2010) 1272–1280
Turbulence dissipation rate:
@ 1 @ ðqueÞ þ ðr qv eÞ @x r @r
physical properties are generated using IAPWS-95 and are incorporated in the CFD code.
@ l @e ¼ lþ t @x re @x 1 @ l @e þ r lþ t r @r re @r 1 qe þ qE þ C e1 f1 ðPk þ Gk Þ C e2 f2 Tt Tt
4. Results and discussion
ð9Þ where T t ¼ k=e. The V2F model has two additional equations for v 2 and f, which are shown below: scale v 2: Turbulent velocity @ 1 @ 2 r qvv 2 quv þ @x r @r " # " # @ l @v 2 1 @ l @v 2 e þ þ kf 6v 2 r lþ t lþ t ¼ @x r @r re @x re @r k
ð10Þ Production f:
0¼
@ @f 1 @ @f 1 ð2=3 v 2 =kÞ þ r 2 f ðC 1 1Þ @x @x r @r @r Tt L þ
C2 1 1 5v 2 =k P þ 2 2 qk k Tt L L
1275
ð11Þ 2
where the turbulent viscosity is defined as: lt ¼ qC l v 2 T t ¼ qC l fl ke , h h 3=2 i qffiffiffiffii 3 where T t ¼ max ke ; 6 qle , L ¼ C L max k e ; C l ððl=eqÞ Þ1=4 , C 1 ¼ 1:4, 2
C 2 ¼ 0:3, C l ¼ 0:22, C L ¼ 0:23, fl ¼ v 2 T t =ðk =eÞ. Other model constants (C l , C e1 , C e2 , rk , re ), damping function (f2 , f2 , fl ) and boundary conditions at wall are specified in Tables 2–4. The complete computational domain, which covered the whole unheated and heated lengths of the test section and ranged from the centre of the tube to the inner wall, was discretized into a mesh of grids, typically, 120 90 (axial radial). The mesh was refined in the radial direction towards the tube wall and in the axial direction in the region where the heating commenced. The grid size is adjusted in each individual run to satisfy the requirement y+ < 0.2 for the first node close to the wall. The sensitivity of calculation results on mesh structure is carefully checked and the nearly grid independent solutions are obtained. The QUICK scheme is used for discretization of momentum and energy equations and the first order upward scheme is used to discretize the turbulence equations to achieve numerical stability. The SIMPLE scheme is used for coupling the pressure and the velocity fields. The convergence criterion for normalized residual of individual equation is set to be less than 106. The effect of temperature is considered, while the effect of pressure on physical properties is neglected since the pressure drop along the ducts is very small. The tables of thermal–
Table 2 Constants in the turbulent models. Models
Cl
C e1
C e2
rk
re
Abe–Kondoh–Nagano (1994) (AKN) Yang–Shih (1993) (YS) Chang–Hsieh–Chen (1995) (CHC) Abid (1993) (AB) Behnia–Parneix–Durbin (1998) (V2F) Menter (1993) (SST)
0.09 0.09 0.09 0.09 0.22 0.09
1.50 1.44 1.5 1.45 1.4 1.55
1.90 1.92 1.90 1.83 1.90 1.83
1.4 1.0 1.4 1.0 1.0 2.0
1.4 1.3 1.4 1.4 1.3 2.0
Figs. 2 and 3 show the comparison between the predictions and the measurements of wall temperature for run1 with upward flow and run2 with downward flow, respectively, where the wall temperature is plotted against the bulk enthalpy. It can be seen that, for run1, the experimental wall temperature increases smoothly at first and then exhibits an abrupt increase when Tw reaches Tpc, which is typical HTD phenomenon. After the peak temperature is achieved, the Tw reduces rapidly along the tube. All turbulent models, except CHC model, are able to predict the wall temperature well in the normal heat transfer region and to capture the general trend of HTD. All k–e type models seriously over-predict the deterioration while the SST model considerably underestimates the peak in wall temperature and delays its appearance, which is in a good agreement with the condition presented by He et al. (2008a). The V2F model performs the best in terms of the prediction of deterioration onset. The deterioration also occurs for run2 with downward flow after Tw > Tpc, but the magnitude of wall temperature is considerably lower than that for run1 as shown in Fig. 3. All models can qualitatively reproduce this heat transfer behavior in downward flow. But CHC, ABID, YS, and AKN models significantly over-estimate the wall temperature in downstream region. The V2F and SST model can predict the wall temperature quite well. Run1 and run2 have same flow rate and wall heat flux. The only difference between the two runs is the buoyancy influence due to the opposite flow direction. Consequently the buoyancy influence must be responsible for the large localized increase of wall temperature observed in the upward flow case. According to Mikielewicz et al. (2002), the effect of buoyancy is negligible for normal fluids such as air or water at atmospheric pressure if the value of buoyancy parameter Bo* based on bulk fluid properties is less than 6.0E7. However, the value of Bo* for run1 is in the range of 5.39E7 to 5.75E7. Clearly, the onset of significant buoyancy influence for supercritical water is very different from that for normal fluid. This is due to the special fluid properties variation and the flow fields in supercritical water. Fig. 4 gives the predictions by V2F model for run1 and run2. The buoyancy effect due to the flow direction on wall temperature is significant. Five axial positions are selected for both runs, as shown in Fig. 4, to present details of the flow and thermal fields based on the modeling results. For run1 with upward flow, P1 is the point selected from the normal heat transfer region and P2 is chosen from the onset of HTD. P3 and P4 are in the deterioration region, and P5 is in the heat transfer recovery region. Fig. 5 presents the predicted radial distributions of fluid temperature, fluid properties variation and flow field at five axial locations obtained using the V2F model for run1 with upward flow. Although the fluid temperature varies sharply near the wall, the sharp variation of fluid properties occurs in the region near the pseudo-critical value, referred to as the large-property-variation (LPV) region (He et al., 2008a). This is a special feature for fluids at supercritical pressure. To specify the region of LPV at different positions, a band of temperature specified by Tpc ± 5 °C is marked in Fig. 5a. The LPV is located within a narrow layer very close to the wall at upstream positions P1 and P2. The LPV moves away from the wall and spreads more widely at downstream locations as P3, P4 and P5. The specific heat distributions at different locations are shown in Fig. 5b. The range of large specific heat at downstream positions is much wider than that at upstream positions due to the spread of LPV region.
1276
Q.L. Wen, H.Y. Gu / Annals of Nuclear Energy 37 (2010) 1272–1280
Table 3 Functions in the turbulent models. Models AKN
f1
f2
fl
1.0
n 2 oh i2 y a t 1 0:3 exp Re 1 exp 3:1 6:5 pffiffiffiffiffiffi t pReffiffiffiffiffi ffi
n
pffiffiffiffiffiffi t pReffiffiffiffiffi ffi
YS
1þ
a
CHC
1.0
ABID
1.0
Ret
1þ
V2F
1 þ 0:045
SST
1.0 2
1=2
Ret ¼ qlke , Rey ¼ yqlk , yþ ¼ ylq
Ret
h i2 1 0:0631 exp 0:0631Rey 1 0:01 exp Re2t h 2 i2 h i Re t 1 29 exp Re 1 exp 12y 36
qffiffiffiffi k
v2
2 oh i2 y Ret 1 þ Re10:75 exp 200 1 exp 1:4 t 2 0 130:5 1:5 104 Rey pffiffiffiffiffiffiffiffi 6 B 5:0 107 Re3 C7 1 þ 1= Ret 41 exp @ y A5 1:0 1010 Re5y 2 1 0:01 exp 0:0215Rey 1 þ 31:66Re1:125 t 3=4
tanhð0:008Rey Þð1 þ 4Ret
1.0
v
1.0
1.0
Þ
2
k
pffiffiffiffiffiffiffiffiffiffiffiffi yq u ffi , sw =q, y ¼ l ue , uþ ¼ pffiffiffiffiffiffiffiffi sw = q
sw ¼ l du , ue ¼ ðle=qÞ0:25 . dy
Table 4 Damping terms D and E and wall boundary conditions for k and e. Models
D
E
AKN
0.0
0.0
Wall BC
YS
0.0
CHC
0.0
ABID
kw ¼ 0,
ew ¼ 2 lq yk2
kw ¼ 0,
ew ¼ 2 lq
0.0
kw ¼ 0,
ew ¼ lq
0.0
0.0
kw ¼ 0,
ew ¼ lq
V2F
0.0
ew ¼ lq @@y2k
0.0
0.0 h i l @k @ e @ e 2 lq þ 2qt @k @y þ @x @y k þ @x k
kw ¼ 0,
SST
ll
2 q2t
2
@ v @x2
2
þ
2
@ v @y2
2
pffiffi 2 @ k @y
@2 k @y2 @2 k @y2
2
2 kw ¼ 0, ew ¼ lq @@y2k
Fig. 4. Comparison of predictions by V2F model with experiment for upward and downward flow (run1 and run2).
Fig. 2. Comparison of predictions with experiment for upward flow (run1).
Fig. 3. Comparison of predictions with experiment for downward flow (run2).
The buoyancy influence in a heated tube is mainly determined by the non-uniformity distribution of density along the radial direction. Fig. 5c shows the density distributions at different axial
locations. At location P1, the density varies sharply in vicinity very close to the wall, which leads to a strong local buoyancy influence. However, the influence is limited to the viscous sub-layer (y+ < 5) and turbulence is not significantly modified. So the buoyancy effect on the heat transfer at P1 is not significant, as supported by wall temperature results shown in Fig. 4. Proceeding to downstream, at location P2, the LPV region moves away from the wall to the buffer layer region. Further downstream, at locations P3, P4 and P5, the LPV region moves further away from the wall and the spread is much wider. The movement of LPV region leads to such a distribution that the lighter water occupy the near-wall region while heavier water stay in the core flow region as shown in Fig. 5c. The density difference between the two regions generates strong buoyancy force to modify the flow field. Fig. 5d gives the velocity profiles at different locations. At location P1, the buoyancy influence is negligible, so the velocity profile is a typical forced convection type. Proceeding to downstream, the velocity profile becomes flattened at locations P2 and P3 where the effect of buoyancy begins to be evident. As the flow proceeds to downstream at location P4, the distortion of the velocity profile due to buoyancy effect is more evident, with the profiles becoming inverted in the core flow region and presenting an M-shape. At P5 of downstream, the velocity profile maintains the M-shape but the extent to which the profile is inverted reduces. The velocity distribution significantly affects the turbulent characteristics and subsequently changes the heat transfer behavior. Fig. 5e–g shows the radial profiles of turbulent shear stress, turbulent shear production and turbulent kinetic energy, respectively, at different axial locations. At location P1, the buoyancy force is limited in the viscous sub-layer, the buoyancy effect is negligible. Significant buoyancy influence can be observed at location P2, where the turbulent shear tress, turbulent shear production and
Q.L. Wen, H.Y. Gu / Annals of Nuclear Energy 37 (2010) 1272–1280
(a) Fluid temperature
1277
(b) Specific heat capability
(c) Density
(d) Velocity
(e) Turbulent shear stress
(f) Turbulent shear production
(g) Turbulent kinetic energy
(h) Turbulent conductive coefficient
Fig. 5. Fluid properties and flow field near-wall region for upward flow (run1).
turbulent kinetic energy are obviously reduced. Subsequently, the onset of deterioration occurs. Stronger buoyancy influence can be observed in following positions P3 and P4, where the shear stress
and turbulent kinetic energy are significantly suppressed. The turbulent conductive coefficient decreases sharply at the locations of P3 and P4 due to the buoyancy effect as shown in Fig. 5h, resulting
1278
Q.L. Wen, H.Y. Gu / Annals of Nuclear Energy 37 (2010) 1272–1280
(a) Fluid temperature
(b) Specific heat capability
(c) Density
(d) Velocity
(e) Turbulent shear stress
(f) Turbulent shear production
(g) Turbulent kinetic energy
(h) Turbulent conductive coefficient
Fig. 6. Fluid properties and flow field near-wall region for downward flow (run2).
in strong HTD. Proceeding further to downstream, at the location P5, turbulent shear stress is negative with the help of even stronger buoyancy effect. The turbulent shear stress, turbulent production
and turbulent kinetic energy start to recover near the wall as seen in Fig. 5. As a result, the turbulent recovery increases the heat transfer and decreases the wall temperature as shown in Fig. 4.
Q.L. Wen, H.Y. Gu / Annals of Nuclear Energy 37 (2010) 1272–1280
Fig. 6 shows the fluid properties and flow field near-wall region for run2 with downward flow. The fluid properties vary sharply near the wall. For the upstream positions of P1, P2 and P3, the LPV regions are limited to a region very close to the wall and subsequently the buoyancy effect on turbulence is negligible. Further downstream, at the positions P4 and P5, the LVP spreads to the buffer layer and the buoyancy effect on turbulence becomes evident as shown in Fig. 6e–g. For downward flow, the direction of the buoyancy force is opposite to flow direction, so the buoyancy force decreases the velocity of lighter water near the wall. On the other hand, the volume expansion of heated water leads to flow accelerations in the near-wall region. The turbulent shear and kinetic energy are reduced, caused by the above two canceling effects. However, the reduction is much smaller than that for upward flow. So the HTD is also observed at downstream position for downward
Fig. 7. Comparison of predictions with experiment for upward flow (run3).
(a) Specific heat capability
(c) Turbulent shear stress
1279
flow, although the deterioration is much weaker than that for upward flow, as shown in Fig. 4. The experimental data of Shitsman (1963) clearly shows two HTD phenomena with two peaks in the wall temperature: one close to the entrance and the other broader. The wall temperature decrease downstream the location of each peak indicates recovery of heat transfer as shown in Fig. 7. It can be seen that each deterioration phenomenon occurs under the condition that the wall temperature is higher than pseudo-critical point temperature Tpc while the bulk fluid temperature is lower than Tpc. All k–e type models used in this study can predict the onset of the first deterioration phenomenon. However, those models under-predict the heat transfer recovery significantly and subsequently seriously overpredict the wall temperature at the downstream of the first peak location. However, only the V2F and YS models show a little deeper valley between the two peaks Remarkably, the SST model can predict interesting heat transfer behavior quantitatively and qualitatively as shown in Fig. 7. The properties and flow field at four axial positions, as marked in Fig. 7, are presented based on the numerical results by SST model. Fig. 8 gives the properties variation and flow fields at different axial locations. It can be seen that the flow behavior due to the buoyancy influence in both HTD phenomena are similar to that in run1 as described above. With the help of the strong buoyancy, the velocity profile with M-shape promotes very effective renewed turbulent production in the near-wall region and the heat transfer recovery. Consequently, the wall temperature decreases significantly from the wall to the first peak and the LVP moves to the near-wall region again as can be seen from the property profile distributions shown in Fig. 8a. At the same time, the velocity in core increases due to the increase in bulk temperature, which makes the velocity gradient in the near-wall region decreases again as
(b) Velocity profile
(d) Turbulent kinetic energy
Fig. 8. Fluid properties and flow field near-wall region for upward flow (run3).
1280
Q.L. Wen, H.Y. Gu / Annals of Nuclear Energy 37 (2010) 1272–1280
shown in Fig. 8b. The turbulent production near the wall decreases, causing the second HTD phenomenon. After the wall temperature reaches the second peak point, the strong buoyancy force recovers the heat transfer again. As the bulk fluid temperature increases, the non-uniformity of density distribution along the radial direction decreases, so the effect of buoyancy on the second heat transfer is much weaker than that on the first one. As a result, the deterioration of heat transfer with low fluid bulk temperature is stronger than that with relatively high fluid bulk temperature. This is also proved by the experimental results by Kim et al. (2005).
5. Conclusions In this study, HTD phenomenon is numerically studied using six low-Reynolds models including AKN, CHC, ABID, YS, V2F and SST to provide detailed information on flow and turbulence. The main conclusions obtained are as follows: (1) All low-Reynolds models considered are to some extent able to reproduce the effect of buoyancy on heat transfer. But all k–e type models seriously over-predict the heat transfer deterioration and do not reproduce the subsequent heat transfer recovery. V2F and SST models perform better than other models in predicting the onset of deterioration due to strong buoyancy force. (2) The second heat transfer deterioration phenomenon may occur if the full recovery of the heat transfer after the first heat transfer deterioration is achieved with the help of very strong buoyancy force. Among the six models used in present study, only the V2F and YS models show a little deeper valley between the two peaks. In particular, the SST model can qualitatively predict this interesting phenomenon. (3) Large property variation in the vicinity of the pseudo-critical temperature causes significant axial density non-uniformity between the near-wall region and the core region in tube, resulting in deformation of velocity profile and impairment of turbulent production. Strong buoyancy effect may occur when the buoyancy parameter Bo* based on bulk flow properties is less than the limiting value of Bo* for onset of buoyancy effects for a normal fluid.
Acknowledgements The authors are grateful to National Natural Science Foundation (No. 50806044) and National Basic Research Program of China (No. 2007CB209804) for providing the financial support for this study. References Abe, K., Kondoh, T., Nagano, Y., 1994. A new turbulence model for predicting fluid flow and heat transfer in separating and reattaching flows – I. Flow field calculations. Int. J. Heat Mass Transfer 37, 139–151. Abid, R., 1993. Evaluation of two-equation turbulence models for predicting transitional flows. Int. J. Eng. Sci. 31, 831–840. Bae, J.H., Yoo, J.Y., Choi, H., 2005. Direct numerical simulation of turbulent supercritical flows with heat transfer. Phys. Fluids 17, 1–24. Behnia, M., Parneix, S., Durbin, P.A., 1998. Prediction of heat transfer in an axisymmetric turbulent jet impinging on a flat plate. Int. J. Heat Mass Transfer 41, 1845–1855. Chang, K.C., Hsieh, W.D., Chen, C.S., 1995. A modified low-Reynolds number turbulence model applicable to recirculating flow in pipe expansion. Trans. ASME J. Fluid Eng. 117, 417–423.
Cheng, X., Schulenburg, T., 2001. Heat transfer at supercritical pressure. Literature review and application to a HPLWR. Wissenschaftliche Berichte FZKA 6609, Forschungszentrum, Karlsruhe, May, 2001. Cheng, X., Kuang, B., Yang, Y.H., 2007. Numerical analysis of heat transfer in supercritical water-cooled flow channels. Nucl. Eng. Des. 237, 240–252. Chien, K.Y., 1982. Predictions of channel and boundary-layer flows with a lowReynolds number turbulence model. AIAA J. 20, 33–38. He, S., Kim, W.S., Jiang, P.X., Jackson, J.D., 2004. of mixed convection heat transfer to carbon dioxide at supercritical pressure. J. Mech. Eng. Sci. 218, 1281–1296. He, S., Kim, W.S., Jackson, J.D., 2008a. A computational study of convective heat transfer to carbon dioxide at a pressure just above the critical value. Appl. Therm. Eng. 28, 1662–1675. He, S., Kim, W.S., Bae, J.H., 2008b. Assessment of performance of turbulence models in predicting supercritical pressure heat transfer in a vertical tube. Int. J. Heat Mass Transfer 51, 4659–4675. Jackson, J.D., 2001. Some striking features of heat transfer with fluids at pressures and temperatures near the critical point. In: International Conference on Energy Conversion and Application (ICECA-2001), Wuhan, China, June 6–12. Jackson, J.D., Hall, W.B., 1979. Forced convection heat transfer to fluids at supercritical pressure. Turbulent Forced Convection in Channels and Bundles, vol. 2. Hemisphere, New York, pp. 563–611. Jones, W.P., Launder, B.E., 1972. The prediction of laminarization with two-equation models of turbulence. Int. J. Heat Mass Transfer 15, 301–314. Kim, J.K., Jeon, H.K., Yoo, J.Y., Lee, J.S., 2005. Experimental study on heat transfer characteristics of turbulent supercritical flow in vertical circular/non-circular tubes. In: Proceedings of the 11th NURETH-11, Avignon, France, October 2–6. Koshizuka, S., Takano, N., Oka, Y., 1995. Numerical analysis of deterioration phenomena in heat transfer to supercritical water. Int. J. Heat Mass Transfer 38, 3077–3084. Lam, C.K.G., Bremhorst, K.A., 1981. Modified form of k–e model for predicting wall turbulence. ASME J. Fluids Eng. 103, 456–460. Licht, J., 2008. Heat Transfer and Fluid Flow Characteristics in Supercritical Water. Doctoral dissertation. University of Wisconsin-Madison, USA. Lien, F.S., Kalitzin, G., Durbin, P.A., 1998. RANS modeling for compressible and transitional flow. In: Proceedings of the Summer Program. Stanford University, USA, p. 267. McEligot, D.M., Jackson, J.D., 2004. ‘‘Deterioration” criteria for convective heat transfer in gas flow through non-circular ducts. Nucl. Eng. Des. 232, 327–333. Menter, F.R., 1993. Zonal two equation k–x turbulence models for aerodynamic flows. In: AIAA 24th Fluid Dynamics Conference, Orlando, Florida USA, July 6–9. Mikielewicz, D.P., Shehata, A.M., Jackson, J.D., McEligot, D.M., 2002. Temperature, velocity and mean turbulence structure in strongly heated internal gas flows: comparison of numerical predictions with data. Int. J. Heat Mass Transfer 45, 4333–4352. Myoung, H.K., Kasagi, N., 1990. A new approach to the improvement of k–e turbulence model for wall bounded shear flows. JSME Int. J. 33, 63–72. Ogata, K., Sato, S., 1972. Measurement of forced convection heat transfer to supercritical helium. In: Proc. of the 4th. Int. Cryogenic Engineering Conference, Eindhoven, Netherlands, May 24–26. Oka, Y., Ishiwatari, Y., Koshizuka, S., 2007. Research and development of super LWR and super fast reactor. In: 3rd Int. Symp. Supercritical Water-cooled Reactors – Design and Technology, Shanghai, China, March 12–15. Palko, D., Anglart, H., 2008. Theoretical and Numerical Study of Heat Transfer Deterioration in High Performance Light Water Reactor, Science and Technology of Nuclear Installations, vol. 2008. Hindawi Publishing Corporation, 5 p. doi:10.1155/2008/405072 (Article ID 405072). Petuhkov, B.S., Kurganov, V.A., 1983. Heat transfer and flow resistance in the turbulent pipe flow of a fluid with near-critical state parameters. Teplofiz. Vys. Temp. 21 (January (1)), 92–100. Pioro, I.L., Duffey, R.B., 2005. Experimental heat transfer in supercritical water flowing inside channels (survey). Nucl. Eng. Des. 235, 2407–2430. Pis’menny, E.N., Razumovskiy, V.G., Maevskiy, A.E., Pioro, I.L., 2006. Heat transfer to supercritical water in gaseous state or affected by mixed convection in vertical tubes. In: Proceedings of the ICONE14, Miami, USA, July 17–20. Schulenberg, T., Starflinger, J., Heinecke, J., 2008. Three pass core design proposal for a high performance light water reactor. Prog. Nucl. Energy 50, 526–531. Sharabi, M., Ambrosini, W., He, S., Jackson, J.D., 2008. Prediction of turbulent convective heat transfer to a fluid at supercritical pressure in square and triangular channels. Ann. Nucl. Energy 35, 993–1005. Shitsman, M.E., 1963. Impairment of the heat transmission at supercritical pressures. High Temp. 1, 237–244. Styrikovich, M.A., Margulova, T.K., Miropolskiy, Z.L., 1967. Problem in the development of designs of supercritical boilers. Teploenergetika 14, 4–7. Yakhot, V., Orszag, S.A., 1986. Renormalization group analysis of turbulence. I. Basic theory. J. Sci. Comput. 1, 1–51. Yamagata, K., Nishikawa, K., Hasegawa, S., Fuji, T., Yoshida, S., 1972. Forced convection heat transfer to supercritical water flowing in tubes. Int. J. Heat Mass Transfer 15, 2575–2593. Yang, Z., Shih, T.H., 1993. New time scale based k–e model for near-wall turbulence. AIAA J. 31, 1191–1198.
Contents lists available at ScienceDirect
Annals of Nuclear Energy journal homepage: www.elsevier.com/locate/anucene
Numerical simulation of heat transfer deterioration phenomenon in supercritical water through vertical tube Q.L. Wen a,b,*, H.Y. Gu a a b
School of Nuclear Science and Engineering, Shanghai Jiao Tong University, Shanghai 200240, China Reactor Engineering Research Division, Nuclear Power Institute of China, Chengdu 610041, China
a r t i c l e
i n f o
Article history: Received 23 March 2010 Received in revised form 15 May 2010 Accepted 31 May 2010 Available online 1 July 2010 Keywords: Supercritical water Heat transfer deterioration Turbulence model Numerical analysis
a b s t r a c t In this study, a numerical investigation of heat transfer deterioration (HTD) in supercritical water flowing through vertical tube is performed by using six low-Reynolds number turbulence models. All low-Reynolds models can be extended to reproduce the effect of buoyancy force on heat transfer and show the occurrence of localized HTD. However, most k–e models seriously over-predict the deterioration and do not reproduce the subsequent recovery of heat transfer. The V2F and SST models perform better than other models in predicting the onset of deterioration due to strong buoyancy force. The SST model is able to quantitatively reproduce the two heat transfer deterioration phenomena with low mass flux which have been found in the present study. Ó 2010 Elsevier Ltd. All rights reserved.
1. Introduction Supercritical water-cooled reactor (SCWR) has been regarded as an innovative reactor and selected as one of candidates of Generation IV reactor systems based on the criteria of economics, safety, and sustainable production of efficient nuclear power. During the last fifteen years, research activities are ongoing worldwide for developing of the advanced nuclear power plants as SCWR (Oka et al., 2007; Schulenberg et al., 2008). One of the main features of supercritical water is the strong variation of its thermal–physical properties in vicinity of the pseudo-critical temperature. Fig. 1 shows the thermal–physical properties of water at 23.5 MPa. The large variation of thermal–physical properties with temperature leads to strong coupling between the flow and the temperature fields, and creates significant effect by buoyancy force or thermal-induced flow acceleration which affect both the mean flow fields and the turbulence structures. These conditions lead to unusual heat transfer characteristics such as HTD phenomenon (Jackson, 2001). Such phenomenon is generating gradually, but can result in a substantial increase of the fuel cladding temperature. Therefore, an in-depth understanding of the phenomenon is essential to specify cladding temperature limits for the design of the SCWR core. Studies on the thermal–hydraulic behavior of supercritical fluids have been performed since the 1950s. Existing experimental * Corresponding author at: School of Nuclear Science and Engineering, Shanghai Jiao Tong University, Shanghai 200240, China. Tel.: +86 21 34204917. E-mail address: [email protected] (Q.L. Wen). 0306-4549/$ - see front matter Ó 2010 Elsevier Ltd. All rights reserved. doi:10.1016/j.anucene.2010.05.022
and theoretical studies on HTD at supercritical pressure conditions were reviewed and published by several authors (Cheng and Schulenburg, 2001; Pioro and Duffey, 2005). As mentioned in the previous paragraph, the increase in the wall temperature under heat transfer deterioration condition is much slower and smoother than that in the case of the boiling crisis in classical Light Water Reactors (LWRs). Therefore, it is difficult to define the onset of HTD. In the literature, different definitions were used, most of which are based on the ratio of the heat transfer coefficient to a reference value (Yamagata et al., 1972; Styrikovich et al., 1967; Jackson and Hall, 1979; Ogata and Sato, 1972; Petuhkov and Kurganov, 1983; McEligot and Jackson, 2004). A detailed derivation in predictions between different correlations is obtained as indicated by Cheng et al. (2007). Due to the measurement technique limitation and costs associated with experimental study, detailed information on experimental flow field and turbulence structure is rare for water at supercritical pressure (Licht, 2008). With the recent developed technologies in computational methods, numerical investigations using Computational Fluid Dynamics (CFD) codes have achieved a significant attention with the purpose to obtain prediction and to provide a better understanding of the heat transfer mechanism. In recent years, some successful cases have been achieved. The main difficulties in the numerical analysis are related to the turbulence modeling at supercritical pressures. During the last few years, the low-Reynolds number type eddy viscosity turbulence models are used to investigate the HTD. Koshizuka et al. (1995) used JL turbulence model (Jones and Launder, 1972) to study the HTD phenomenon in supercritical
Q.L. Wen, H.Y. Gu / Annals of Nuclear Energy 37 (2010) 1272–1280
1273
Nomenclature Bo* buoyancy parameter Cl, Ce1, Ce2 turbulence models’ constants specific heat transfer (J kg1 K1) Cp D pipe diameter (m)/damping function defined in Table 4 E damping function defined in Table 4 f1, f2, fl damping functions Gr Grashof number g gravity accelerations (m s2) G mass flux (kg m2 s1) h heat transfer coefficient (W m2 K1)/enthalpy (J kg1) m mass flow rate (kg s1) k turbulent kinetic energy (m2 s2) Nu Nusselt number P pressure (MPa) Pk turbulent shear production (W m3) Pr Prandtl number r, x radial and axial coordinates (m) q heat flux (W m2) Re Reynolds number
T U, V y y+, y*
temperature (K) velocity components in the x-, r-directions (m s1) distances in coordinates (m) non-dimensional distance from pipe wall
Greek letters rate of dissipation of k (m2 s3) k thermal conductivity (W m1 K1) l viscosity (kg m1 s1) s shear stress (N m2) q density (kg m3)
e
Subscripts b bulk t turbulent w wall pc pseudo-critical
Fig. 1. Thermal–physical properties for water at 23.5 MPa.
water. The buoyancy effect at low flow rate and flow acceleration effect at high flow rate on heat transfer were analyzed. Calculation results qualitatively agreed with the experimental data of Yamagata et al. (1972). He et al. (2004) applied six low-Reynolds number turbulence models to simulate the turbulent mixed convection heat transfer experiments with carbon dioxide at supercritical pressure. The results showed that all models were able to simulate the influence of buoyancy to heat transfer, but those models have significantly different performances to predict onset and magnitude of this effect. Sharabi et al. (2008) used five low-Reynolds turbulence models, i.e. ABID model (Abid, 1993), LB (Lam and Bremhorst, 1981), AKN model (Abe et al., 1994), CHC model (Chang et al., 1995), YS (Yang and Shih, 1993) and RNG model with the two-layer approach (Yakhot and Orszag, 1986), implementing by FLUENT, for prediction of three-dimensional turbulent heat transfer to carbon dioxide at supercritical pressure flowing upward through heated square and triangular channels as presented by Kim et al. (2005). The results showed that the low-Reynolds turbulence models were able to reproduce the trends of heat transfer deterioration due to buoyancy influence which may lead to relatively large variations of wall temperatures. Among the low-Reynolds number models, the YS model (Yang and Shih, 1993) was able to qualitatively reproduce
the second localized deterioration phenomenon observed in the experimental data downstream of the location of first one. On the other hand, the RNG model with the two-layer approach was unable to capture even the expected qualitative trends. It was found that the flow acceleration due to density non-uniformity has insignificantly affected on heat transfer deterioration for the cases with strong buoyancy influence as considered in the study. He et al. (2008a) made simulations for vertical upward and downward flow of carbon dioxide in a heated tube at a pressure just above the critical value using the V2F turbulence model (Lien et al., 1998) and AKN model (Abe et al., 1994). The results showed that both models can qualitatively show the general trends of HTD due to strong buoyancy influence in upward flow, but significant discrepancies in detailed wall temperature distribution are found. The AKN model is the best of the two for the conditions considered. In the study, they found that buoyancy influence on heat transfer was almost entirely due to the shear production effect caused by the distortion of the mean flow as a result of the influence of buoyancy. Palko and Anglart (2008) used the low-Reynolds k–x turbulence model to make simulations of heat transfer deterioration. The results showed for the low mass flux of coolant, the influence of buoyancy forces on the heat transfer in heated pipes is signifi-
1274
Q.L. Wen, H.Y. Gu / Annals of Nuclear Energy 37 (2010) 1272–1280
cant. For the high mass flux, buoyancy influence could be neglected and there are clearly other mechanisms causing the decrease in heat transfer at high coolant mass flux. Recently, He et al. (2008b) assessed the performance of seven low-Reynolds number turbulence models in responding to the modifications of the turbulence field due to influences of flow acceleration and buoyancy to carbon dioxide at supercritical pressure by comparing model predictions with DNS results by Bae et al. (2005). The results showed that the low-Reynolds number turbulence models whose damping functions are based on variables readily responding to buoyancy and flow acceleration (LS, YS, AKN), significantly over-predict flow laminarization and therefore, the HTD. However, turbulence models whose damping functions are based on variables not responding to buoyancy/flow acceleration (CH (Chien, 1982) and MK (Myoung and Kasagi, 1990) models) reproduce closely the variations of wall temperature exhibited in the DNS in flow laminarization cases due to some canceling effects. The V2F model produces the best predictions among all the turbulence models tested. This paper is devoted to computational simulations of experiments on HTD phenomena with a well-validated commercial CFD software FLUENT, using k–e type and k–x type low-Reynolds turbulence models. The purpose of the work presented here is to numerically study the interesting behavior exhibited in experimental results using detailed information on the flow and thermal fields. 2. Experimental data used The experimental results considered in this study (Pis’menny et al., 2006; Shitsman, 1963) were obtained with water at supercritical pressure flowing in a uniformly heated vertical tube. The experimental conditions which have been simulated in this study are summarized in Table 1. Run1 and run3 are in the HTD region due to strong buoyancy influence. The experimental condition of run2 is the same as that of run1, except for the flow direction. More experimental details can be found in Pis’menny et al. (2006) and Shitsman (1963).
V-momentum:
1 @ @ ðr quv Þ þ ðr qv 2 Þ r @x @r @p 1 @ @ v @u @ @v lv r le þ2 r le þ 2 e2 ¼ þ @x r @x @r @x @r @r r ð3Þ where le is the effective viscosity defined by the turbulent viscosity defined as
lt ¼ qC l fl
k
le ¼ l þ lt and lt is
2
ð4Þ
e
in which fl is a damping function to account for near-wall effects and C l is a constant. Energy:
1 @ @ ðr quhÞ þ ðrqv hÞ r @x @r
¼
1 @ l lT @h r þ r @x Pr rT @x @ l lT @h þ r þ @r Pr rT @r
ð5Þ
where Pr is molecular Prandtl number and rT is the turbulent Prandtl number, the constant 0.9 is used in this study. Six turbulence models have been selected, namely, five k–e type low-Reynolds turbulence models due to Abe et al. (1994) (AKN), Yang and Shih (1993), Chang et al. (1995) (CHC) and Abid (1993) (AB), the k–e–v2–f model of Behnia et al. (1998) (V2F) and one k– x type low-Reynolds turbulence model, i.e., shear stress transport model (SST) developed by Menter (1993), the SST model combines the robustness of the k–e model for the bulk flow with the low-Re treatment of the boundary layer using the modified Wilcox k–x model. The general form of k–e models can be expressed as follows: Turbulent kinetic energy:
3. Governing equations and turbulence models
@ 1 @ ðqukÞ þ ðr qv kÞ @x r @r @ l @k 1 @ l @k lþ t þ r lþ t þ Pk þ Gk qðe DÞ ¼ @x rk @x r @r rk @r
In the FLUENT code, the flow is assumed to be 2D axisymmetry and the governing equations of the continuity, momentum and energy for such a flow written in cylindrical coordinates are as follows: Continuity:
ð6Þ where the production term due to the mean velocity gradient is:
1 @ @ ðrquÞ þ ðr qv Þ r @x @r
¼0
ð1Þ
" (
P k ¼ lt 2
2 2 2 ) 2 # @u @v v @u @ v þ þ þ þ @x @r @x @r r
ð7Þ
The gravitational production term is computed by:
U-momentum:
1 @ @ ðr qu2 Þ þ ðr qv uÞ r @x @r @p 1 @ @u @ @u @ v 2 r le þ r le þ ¼ þ qg þ @x r @x @x @r @r @x
Gk ¼ q0 u0 g x ¼ ch bqgx ð2Þ
k
e
lT
@u @ v @T @u 2 @T þ þ 2lT qk @r @x @r @x 3 @x
ð8Þ
where g x ¼ g for upward flow, g x ¼ g for downward flow and ch ¼ 0:3.
Table 1 Experimental conditions.
a
Cases
P (MPa)
G (kg/m2 s)
qw (kW/m2)
D (mm)
Dir.
Boa 107
Authors
Run1 Run2 Run3
23.5 23.5 23.5
509 509 430
390 390 319.87
6.28 6.28 8
Up Down Up
5.39–5.75 – 9.14–4.57
Pis’menny et al. (2006) Pis’menny et al. (2006) Shitsman (1963)
2
Bo ¼ Gr=ðRe3:425 Pr0:8 Þ, Gr ¼ gbD4 qw =kv .
Q.L. Wen, H.Y. Gu / Annals of Nuclear Energy 37 (2010) 1272–1280
Turbulence dissipation rate:
@ 1 @ ðqueÞ þ ðr qv eÞ @x r @r
physical properties are generated using IAPWS-95 and are incorporated in the CFD code.
@ l @e ¼ lþ t @x re @x 1 @ l @e þ r lþ t r @r re @r 1 qe þ qE þ C e1 f1 ðPk þ Gk Þ C e2 f2 Tt Tt
4. Results and discussion
ð9Þ where T t ¼ k=e. The V2F model has two additional equations for v 2 and f, which are shown below: scale v 2: Turbulent velocity @ 1 @ 2 r qvv 2 quv þ @x r @r " # " # @ l @v 2 1 @ l @v 2 e þ þ kf 6v 2 r lþ t lþ t ¼ @x r @r re @x re @r k
ð10Þ Production f:
0¼
@ @f 1 @ @f 1 ð2=3 v 2 =kÞ þ r 2 f ðC 1 1Þ @x @x r @r @r Tt L þ
C2 1 1 5v 2 =k P þ 2 2 qk k Tt L L
1275
ð11Þ 2
where the turbulent viscosity is defined as: lt ¼ qC l v 2 T t ¼ qC l fl ke , h h 3=2 i qffiffiffiffii 3 where T t ¼ max ke ; 6 qle , L ¼ C L max k e ; C l ððl=eqÞ Þ1=4 , C 1 ¼ 1:4, 2
C 2 ¼ 0:3, C l ¼ 0:22, C L ¼ 0:23, fl ¼ v 2 T t =ðk =eÞ. Other model constants (C l , C e1 , C e2 , rk , re ), damping function (f2 , f2 , fl ) and boundary conditions at wall are specified in Tables 2–4. The complete computational domain, which covered the whole unheated and heated lengths of the test section and ranged from the centre of the tube to the inner wall, was discretized into a mesh of grids, typically, 120 90 (axial radial). The mesh was refined in the radial direction towards the tube wall and in the axial direction in the region where the heating commenced. The grid size is adjusted in each individual run to satisfy the requirement y+ < 0.2 for the first node close to the wall. The sensitivity of calculation results on mesh structure is carefully checked and the nearly grid independent solutions are obtained. The QUICK scheme is used for discretization of momentum and energy equations and the first order upward scheme is used to discretize the turbulence equations to achieve numerical stability. The SIMPLE scheme is used for coupling the pressure and the velocity fields. The convergence criterion for normalized residual of individual equation is set to be less than 106. The effect of temperature is considered, while the effect of pressure on physical properties is neglected since the pressure drop along the ducts is very small. The tables of thermal–
Table 2 Constants in the turbulent models. Models
Cl
C e1
C e2
rk
re
Abe–Kondoh–Nagano (1994) (AKN) Yang–Shih (1993) (YS) Chang–Hsieh–Chen (1995) (CHC) Abid (1993) (AB) Behnia–Parneix–Durbin (1998) (V2F) Menter (1993) (SST)
0.09 0.09 0.09 0.09 0.22 0.09
1.50 1.44 1.5 1.45 1.4 1.55
1.90 1.92 1.90 1.83 1.90 1.83
1.4 1.0 1.4 1.0 1.0 2.0
1.4 1.3 1.4 1.4 1.3 2.0
Figs. 2 and 3 show the comparison between the predictions and the measurements of wall temperature for run1 with upward flow and run2 with downward flow, respectively, where the wall temperature is plotted against the bulk enthalpy. It can be seen that, for run1, the experimental wall temperature increases smoothly at first and then exhibits an abrupt increase when Tw reaches Tpc, which is typical HTD phenomenon. After the peak temperature is achieved, the Tw reduces rapidly along the tube. All turbulent models, except CHC model, are able to predict the wall temperature well in the normal heat transfer region and to capture the general trend of HTD. All k–e type models seriously over-predict the deterioration while the SST model considerably underestimates the peak in wall temperature and delays its appearance, which is in a good agreement with the condition presented by He et al. (2008a). The V2F model performs the best in terms of the prediction of deterioration onset. The deterioration also occurs for run2 with downward flow after Tw > Tpc, but the magnitude of wall temperature is considerably lower than that for run1 as shown in Fig. 3. All models can qualitatively reproduce this heat transfer behavior in downward flow. But CHC, ABID, YS, and AKN models significantly over-estimate the wall temperature in downstream region. The V2F and SST model can predict the wall temperature quite well. Run1 and run2 have same flow rate and wall heat flux. The only difference between the two runs is the buoyancy influence due to the opposite flow direction. Consequently the buoyancy influence must be responsible for the large localized increase of wall temperature observed in the upward flow case. According to Mikielewicz et al. (2002), the effect of buoyancy is negligible for normal fluids such as air or water at atmospheric pressure if the value of buoyancy parameter Bo* based on bulk fluid properties is less than 6.0E7. However, the value of Bo* for run1 is in the range of 5.39E7 to 5.75E7. Clearly, the onset of significant buoyancy influence for supercritical water is very different from that for normal fluid. This is due to the special fluid properties variation and the flow fields in supercritical water. Fig. 4 gives the predictions by V2F model for run1 and run2. The buoyancy effect due to the flow direction on wall temperature is significant. Five axial positions are selected for both runs, as shown in Fig. 4, to present details of the flow and thermal fields based on the modeling results. For run1 with upward flow, P1 is the point selected from the normal heat transfer region and P2 is chosen from the onset of HTD. P3 and P4 are in the deterioration region, and P5 is in the heat transfer recovery region. Fig. 5 presents the predicted radial distributions of fluid temperature, fluid properties variation and flow field at five axial locations obtained using the V2F model for run1 with upward flow. Although the fluid temperature varies sharply near the wall, the sharp variation of fluid properties occurs in the region near the pseudo-critical value, referred to as the large-property-variation (LPV) region (He et al., 2008a). This is a special feature for fluids at supercritical pressure. To specify the region of LPV at different positions, a band of temperature specified by Tpc ± 5 °C is marked in Fig. 5a. The LPV is located within a narrow layer very close to the wall at upstream positions P1 and P2. The LPV moves away from the wall and spreads more widely at downstream locations as P3, P4 and P5. The specific heat distributions at different locations are shown in Fig. 5b. The range of large specific heat at downstream positions is much wider than that at upstream positions due to the spread of LPV region.
1276
Q.L. Wen, H.Y. Gu / Annals of Nuclear Energy 37 (2010) 1272–1280
Table 3 Functions in the turbulent models. Models AKN
f1
f2
fl
1.0
n 2 oh i2 y a t 1 0:3 exp Re 1 exp 3:1 6:5 pffiffiffiffiffiffi t pReffiffiffiffiffi ffi
n
pffiffiffiffiffiffi t pReffiffiffiffiffi ffi
YS
1þ
a
CHC
1.0
ABID
1.0
Ret
1þ
V2F
1 þ 0:045
SST
1.0 2
1=2
Ret ¼ qlke , Rey ¼ yqlk , yþ ¼ ylq
Ret
h i2 1 0:0631 exp 0:0631Rey 1 0:01 exp Re2t h 2 i2 h i Re t 1 29 exp Re 1 exp 12y 36
qffiffiffiffi k
v2
2 oh i2 y Ret 1 þ Re10:75 exp 200 1 exp 1:4 t 2 0 130:5 1:5 104 Rey pffiffiffiffiffiffiffiffi 6 B 5:0 107 Re3 C7 1 þ 1= Ret 41 exp @ y A5 1:0 1010 Re5y 2 1 0:01 exp 0:0215Rey 1 þ 31:66Re1:125 t 3=4
tanhð0:008Rey Þð1 þ 4Ret
1.0
v
1.0
1.0
Þ
2
k
pffiffiffiffiffiffiffiffiffiffiffiffi yq u ffi , sw =q, y ¼ l ue , uþ ¼ pffiffiffiffiffiffiffiffi sw = q
sw ¼ l du , ue ¼ ðle=qÞ0:25 . dy
Table 4 Damping terms D and E and wall boundary conditions for k and e. Models
D
E
AKN
0.0
0.0
Wall BC
YS
0.0
CHC
0.0
ABID
kw ¼ 0,
ew ¼ 2 lq yk2
kw ¼ 0,
ew ¼ 2 lq
0.0
kw ¼ 0,
ew ¼ lq
0.0
0.0
kw ¼ 0,
ew ¼ lq
V2F
0.0
ew ¼ lq @@y2k
0.0
0.0 h i l @k @ e @ e 2 lq þ 2qt @k @y þ @x @y k þ @x k
kw ¼ 0,
SST
ll
2 q2t
2
@ v @x2
2
þ
2
@ v @y2
2
pffiffi 2 @ k @y
@2 k @y2 @2 k @y2
2
2 kw ¼ 0, ew ¼ lq @@y2k
Fig. 4. Comparison of predictions by V2F model with experiment for upward and downward flow (run1 and run2).
Fig. 2. Comparison of predictions with experiment for upward flow (run1).
Fig. 3. Comparison of predictions with experiment for downward flow (run2).
The buoyancy influence in a heated tube is mainly determined by the non-uniformity distribution of density along the radial direction. Fig. 5c shows the density distributions at different axial
locations. At location P1, the density varies sharply in vicinity very close to the wall, which leads to a strong local buoyancy influence. However, the influence is limited to the viscous sub-layer (y+ < 5) and turbulence is not significantly modified. So the buoyancy effect on the heat transfer at P1 is not significant, as supported by wall temperature results shown in Fig. 4. Proceeding to downstream, at location P2, the LPV region moves away from the wall to the buffer layer region. Further downstream, at locations P3, P4 and P5, the LPV region moves further away from the wall and the spread is much wider. The movement of LPV region leads to such a distribution that the lighter water occupy the near-wall region while heavier water stay in the core flow region as shown in Fig. 5c. The density difference between the two regions generates strong buoyancy force to modify the flow field. Fig. 5d gives the velocity profiles at different locations. At location P1, the buoyancy influence is negligible, so the velocity profile is a typical forced convection type. Proceeding to downstream, the velocity profile becomes flattened at locations P2 and P3 where the effect of buoyancy begins to be evident. As the flow proceeds to downstream at location P4, the distortion of the velocity profile due to buoyancy effect is more evident, with the profiles becoming inverted in the core flow region and presenting an M-shape. At P5 of downstream, the velocity profile maintains the M-shape but the extent to which the profile is inverted reduces. The velocity distribution significantly affects the turbulent characteristics and subsequently changes the heat transfer behavior. Fig. 5e–g shows the radial profiles of turbulent shear stress, turbulent shear production and turbulent kinetic energy, respectively, at different axial locations. At location P1, the buoyancy force is limited in the viscous sub-layer, the buoyancy effect is negligible. Significant buoyancy influence can be observed at location P2, where the turbulent shear tress, turbulent shear production and
Q.L. Wen, H.Y. Gu / Annals of Nuclear Energy 37 (2010) 1272–1280
(a) Fluid temperature
1277
(b) Specific heat capability
(c) Density
(d) Velocity
(e) Turbulent shear stress
(f) Turbulent shear production
(g) Turbulent kinetic energy
(h) Turbulent conductive coefficient
Fig. 5. Fluid properties and flow field near-wall region for upward flow (run1).
turbulent kinetic energy are obviously reduced. Subsequently, the onset of deterioration occurs. Stronger buoyancy influence can be observed in following positions P3 and P4, where the shear stress
and turbulent kinetic energy are significantly suppressed. The turbulent conductive coefficient decreases sharply at the locations of P3 and P4 due to the buoyancy effect as shown in Fig. 5h, resulting
1278
Q.L. Wen, H.Y. Gu / Annals of Nuclear Energy 37 (2010) 1272–1280
(a) Fluid temperature
(b) Specific heat capability
(c) Density
(d) Velocity
(e) Turbulent shear stress
(f) Turbulent shear production
(g) Turbulent kinetic energy
(h) Turbulent conductive coefficient
Fig. 6. Fluid properties and flow field near-wall region for downward flow (run2).
in strong HTD. Proceeding further to downstream, at the location P5, turbulent shear stress is negative with the help of even stronger buoyancy effect. The turbulent shear stress, turbulent production
and turbulent kinetic energy start to recover near the wall as seen in Fig. 5. As a result, the turbulent recovery increases the heat transfer and decreases the wall temperature as shown in Fig. 4.
Q.L. Wen, H.Y. Gu / Annals of Nuclear Energy 37 (2010) 1272–1280
Fig. 6 shows the fluid properties and flow field near-wall region for run2 with downward flow. The fluid properties vary sharply near the wall. For the upstream positions of P1, P2 and P3, the LPV regions are limited to a region very close to the wall and subsequently the buoyancy effect on turbulence is negligible. Further downstream, at the positions P4 and P5, the LVP spreads to the buffer layer and the buoyancy effect on turbulence becomes evident as shown in Fig. 6e–g. For downward flow, the direction of the buoyancy force is opposite to flow direction, so the buoyancy force decreases the velocity of lighter water near the wall. On the other hand, the volume expansion of heated water leads to flow accelerations in the near-wall region. The turbulent shear and kinetic energy are reduced, caused by the above two canceling effects. However, the reduction is much smaller than that for upward flow. So the HTD is also observed at downstream position for downward
Fig. 7. Comparison of predictions with experiment for upward flow (run3).
(a) Specific heat capability
(c) Turbulent shear stress
1279
flow, although the deterioration is much weaker than that for upward flow, as shown in Fig. 4. The experimental data of Shitsman (1963) clearly shows two HTD phenomena with two peaks in the wall temperature: one close to the entrance and the other broader. The wall temperature decrease downstream the location of each peak indicates recovery of heat transfer as shown in Fig. 7. It can be seen that each deterioration phenomenon occurs under the condition that the wall temperature is higher than pseudo-critical point temperature Tpc while the bulk fluid temperature is lower than Tpc. All k–e type models used in this study can predict the onset of the first deterioration phenomenon. However, those models under-predict the heat transfer recovery significantly and subsequently seriously overpredict the wall temperature at the downstream of the first peak location. However, only the V2F and YS models show a little deeper valley between the two peaks Remarkably, the SST model can predict interesting heat transfer behavior quantitatively and qualitatively as shown in Fig. 7. The properties and flow field at four axial positions, as marked in Fig. 7, are presented based on the numerical results by SST model. Fig. 8 gives the properties variation and flow fields at different axial locations. It can be seen that the flow behavior due to the buoyancy influence in both HTD phenomena are similar to that in run1 as described above. With the help of the strong buoyancy, the velocity profile with M-shape promotes very effective renewed turbulent production in the near-wall region and the heat transfer recovery. Consequently, the wall temperature decreases significantly from the wall to the first peak and the LVP moves to the near-wall region again as can be seen from the property profile distributions shown in Fig. 8a. At the same time, the velocity in core increases due to the increase in bulk temperature, which makes the velocity gradient in the near-wall region decreases again as
(b) Velocity profile
(d) Turbulent kinetic energy
Fig. 8. Fluid properties and flow field near-wall region for upward flow (run3).
1280
Q.L. Wen, H.Y. Gu / Annals of Nuclear Energy 37 (2010) 1272–1280
shown in Fig. 8b. The turbulent production near the wall decreases, causing the second HTD phenomenon. After the wall temperature reaches the second peak point, the strong buoyancy force recovers the heat transfer again. As the bulk fluid temperature increases, the non-uniformity of density distribution along the radial direction decreases, so the effect of buoyancy on the second heat transfer is much weaker than that on the first one. As a result, the deterioration of heat transfer with low fluid bulk temperature is stronger than that with relatively high fluid bulk temperature. This is also proved by the experimental results by Kim et al. (2005).
5. Conclusions In this study, HTD phenomenon is numerically studied using six low-Reynolds models including AKN, CHC, ABID, YS, V2F and SST to provide detailed information on flow and turbulence. The main conclusions obtained are as follows: (1) All low-Reynolds models considered are to some extent able to reproduce the effect of buoyancy on heat transfer. But all k–e type models seriously over-predict the heat transfer deterioration and do not reproduce the subsequent heat transfer recovery. V2F and SST models perform better than other models in predicting the onset of deterioration due to strong buoyancy force. (2) The second heat transfer deterioration phenomenon may occur if the full recovery of the heat transfer after the first heat transfer deterioration is achieved with the help of very strong buoyancy force. Among the six models used in present study, only the V2F and YS models show a little deeper valley between the two peaks. In particular, the SST model can qualitatively predict this interesting phenomenon. (3) Large property variation in the vicinity of the pseudo-critical temperature causes significant axial density non-uniformity between the near-wall region and the core region in tube, resulting in deformation of velocity profile and impairment of turbulent production. Strong buoyancy effect may occur when the buoyancy parameter Bo* based on bulk flow properties is less than the limiting value of Bo* for onset of buoyancy effects for a normal fluid.
Acknowledgements The authors are grateful to National Natural Science Foundation (No. 50806044) and National Basic Research Program of China (No. 2007CB209804) for providing the financial support for this study. References Abe, K., Kondoh, T., Nagano, Y., 1994. A new turbulence model for predicting fluid flow and heat transfer in separating and reattaching flows – I. Flow field calculations. Int. J. Heat Mass Transfer 37, 139–151. Abid, R., 1993. Evaluation of two-equation turbulence models for predicting transitional flows. Int. J. Eng. Sci. 31, 831–840. Bae, J.H., Yoo, J.Y., Choi, H., 2005. Direct numerical simulation of turbulent supercritical flows with heat transfer. Phys. Fluids 17, 1–24. Behnia, M., Parneix, S., Durbin, P.A., 1998. Prediction of heat transfer in an axisymmetric turbulent jet impinging on a flat plate. Int. J. Heat Mass Transfer 41, 1845–1855. Chang, K.C., Hsieh, W.D., Chen, C.S., 1995. A modified low-Reynolds number turbulence model applicable to recirculating flow in pipe expansion. Trans. ASME J. Fluid Eng. 117, 417–423.
Cheng, X., Schulenburg, T., 2001. Heat transfer at supercritical pressure. Literature review and application to a HPLWR. Wissenschaftliche Berichte FZKA 6609, Forschungszentrum, Karlsruhe, May, 2001. Cheng, X., Kuang, B., Yang, Y.H., 2007. Numerical analysis of heat transfer in supercritical water-cooled flow channels. Nucl. Eng. Des. 237, 240–252. Chien, K.Y., 1982. Predictions of channel and boundary-layer flows with a lowReynolds number turbulence model. AIAA J. 20, 33–38. He, S., Kim, W.S., Jiang, P.X., Jackson, J.D., 2004. of mixed convection heat transfer to carbon dioxide at supercritical pressure. J. Mech. Eng. Sci. 218, 1281–1296. He, S., Kim, W.S., Jackson, J.D., 2008a. A computational study of convective heat transfer to carbon dioxide at a pressure just above the critical value. Appl. Therm. Eng. 28, 1662–1675. He, S., Kim, W.S., Bae, J.H., 2008b. Assessment of performance of turbulence models in predicting supercritical pressure heat transfer in a vertical tube. Int. J. Heat Mass Transfer 51, 4659–4675. Jackson, J.D., 2001. Some striking features of heat transfer with fluids at pressures and temperatures near the critical point. In: International Conference on Energy Conversion and Application (ICECA-2001), Wuhan, China, June 6–12. Jackson, J.D., Hall, W.B., 1979. Forced convection heat transfer to fluids at supercritical pressure. Turbulent Forced Convection in Channels and Bundles, vol. 2. Hemisphere, New York, pp. 563–611. Jones, W.P., Launder, B.E., 1972. The prediction of laminarization with two-equation models of turbulence. Int. J. Heat Mass Transfer 15, 301–314. Kim, J.K., Jeon, H.K., Yoo, J.Y., Lee, J.S., 2005. Experimental study on heat transfer characteristics of turbulent supercritical flow in vertical circular/non-circular tubes. In: Proceedings of the 11th NURETH-11, Avignon, France, October 2–6. Koshizuka, S., Takano, N., Oka, Y., 1995. Numerical analysis of deterioration phenomena in heat transfer to supercritical water. Int. J. Heat Mass Transfer 38, 3077–3084. Lam, C.K.G., Bremhorst, K.A., 1981. Modified form of k–e model for predicting wall turbulence. ASME J. Fluids Eng. 103, 456–460. Licht, J., 2008. Heat Transfer and Fluid Flow Characteristics in Supercritical Water. Doctoral dissertation. University of Wisconsin-Madison, USA. Lien, F.S., Kalitzin, G., Durbin, P.A., 1998. RANS modeling for compressible and transitional flow. In: Proceedings of the Summer Program. Stanford University, USA, p. 267. McEligot, D.M., Jackson, J.D., 2004. ‘‘Deterioration” criteria for convective heat transfer in gas flow through non-circular ducts. Nucl. Eng. Des. 232, 327–333. Menter, F.R., 1993. Zonal two equation k–x turbulence models for aerodynamic flows. In: AIAA 24th Fluid Dynamics Conference, Orlando, Florida USA, July 6–9. Mikielewicz, D.P., Shehata, A.M., Jackson, J.D., McEligot, D.M., 2002. Temperature, velocity and mean turbulence structure in strongly heated internal gas flows: comparison of numerical predictions with data. Int. J. Heat Mass Transfer 45, 4333–4352. Myoung, H.K., Kasagi, N., 1990. A new approach to the improvement of k–e turbulence model for wall bounded shear flows. JSME Int. J. 33, 63–72. Ogata, K., Sato, S., 1972. Measurement of forced convection heat transfer to supercritical helium. In: Proc. of the 4th. Int. Cryogenic Engineering Conference, Eindhoven, Netherlands, May 24–26. Oka, Y., Ishiwatari, Y., Koshizuka, S., 2007. Research and development of super LWR and super fast reactor. In: 3rd Int. Symp. Supercritical Water-cooled Reactors – Design and Technology, Shanghai, China, March 12–15. Palko, D., Anglart, H., 2008. Theoretical and Numerical Study of Heat Transfer Deterioration in High Performance Light Water Reactor, Science and Technology of Nuclear Installations, vol. 2008. Hindawi Publishing Corporation, 5 p. doi:10.1155/2008/405072 (Article ID 405072). Petuhkov, B.S., Kurganov, V.A., 1983. Heat transfer and flow resistance in the turbulent pipe flow of a fluid with near-critical state parameters. Teplofiz. Vys. Temp. 21 (January (1)), 92–100. Pioro, I.L., Duffey, R.B., 2005. Experimental heat transfer in supercritical water flowing inside channels (survey). Nucl. Eng. Des. 235, 2407–2430. Pis’menny, E.N., Razumovskiy, V.G., Maevskiy, A.E., Pioro, I.L., 2006. Heat transfer to supercritical water in gaseous state or affected by mixed convection in vertical tubes. In: Proceedings of the ICONE14, Miami, USA, July 17–20. Schulenberg, T., Starflinger, J., Heinecke, J., 2008. Three pass core design proposal for a high performance light water reactor. Prog. Nucl. Energy 50, 526–531. Sharabi, M., Ambrosini, W., He, S., Jackson, J.D., 2008. Prediction of turbulent convective heat transfer to a fluid at supercritical pressure in square and triangular channels. Ann. Nucl. Energy 35, 993–1005. Shitsman, M.E., 1963. Impairment of the heat transmission at supercritical pressures. High Temp. 1, 237–244. Styrikovich, M.A., Margulova, T.K., Miropolskiy, Z.L., 1967. Problem in the development of designs of supercritical boilers. Teploenergetika 14, 4–7. Yakhot, V., Orszag, S.A., 1986. Renormalization group analysis of turbulence. I. Basic theory. J. Sci. Comput. 1, 1–51. Yamagata, K., Nishikawa, K., Hasegawa, S., Fuji, T., Yoshida, S., 1972. Forced convection heat transfer to supercritical water flowing in tubes. Int. J. Heat Mass Transfer 15, 2575–2593. Yang, Z., Shih, T.H., 1993. New time scale based k–e model for near-wall turbulence. AIAA J. 31, 1191–1198.