- Email: [email protected]
CFD investigating flow and heat transfer characteristics in a natural circulation loop
Annals of Nuclear Energy 58 (2013) 65–71
Contents lists available at SciVerse ScienceDirect
Annals of Nuclear Energy journal homepage: www.elsevier.com/locate/anucene
CFD investigating flow and heat transfer characteristics in a natural circulation loop J.Y. Wang, T.J. Chuang, Y.M. Ferng ⇑ Department of Engineering and System Science, Institute of Nuclear Engineering and Science, National Tsing Hua University, 101, Sec. 2 Kuang-Fu Rd., Hsingchu 30013, Taiwan, ROC
a r t i c l e
i n f o
Article history: Received 28 August 2012 Received in revised form 4 January 2013 Accepted 10 January 2013 Available online 3 April 2013 Keywords: Natural circulation Compressible CFD model Passive safety system
a b s t r a c t Natural circulation is the most important heat removal mechanism for passive safety systems in newdesign nuclear power plants. It is important to investigate the flow and heat transfer characteristics related to the natural circulation mechanism. In this paper, a compressible three-dimensional (3-D) CFD model is proposed to investigate these phenomena in a natural circulation loop. The flow and heat transfer behaviors in a loop can be reasonably captured by the present model, which includes the secondary flow within the elbow and the thermal stratification in the horizontal pipe, etc. This model is also validated with the existing experimental data. Compared with the measured results including time histories of local maximum fluid temperature and temperature difference across the heater, and the relationship between Ress and Grm/NG, the present predicted results show good agreement. These comparisons reveal that the present CFD methodology can be applied in simulating the thermal–hydraulic characteristics related to the natural circulation mechanism in confidence. Ó 2013 Elsevier Ltd. All rights reserved.
1. Introduction With advantage of simplification, economy, and no need of power supply, passive safety systems are generally adopted for new-design nuclear power plants. After occurrence of the Fusikama accident, natural circulation mechanism then plays a significant role in the safety related systems since it can be free of electrical demand. However, driving force for the natural circulation mechanism may be not strong enough to deliver the coolant and to remove the heat required for the accident. Therefore, it is crucial to develop analytical models to investigate the flow and heat transfer characteristics related to this heat removal mechanism. Regarding about the researches on natural circulation, watercooled loops had been generally adopted due to their importance in the industry applications, including nuclear power plants, solar water heaters, geothermal processes, machinery cooling, etc. In a natural circulation loop, water is driven by the buoyancy force that is essentially generated due to the temperature-driven density gradient. Heated in the lower portion of loop, water becomes less dense and moves upwards, and subsequently it can be cooled at the higher portion, becomes denser, and flows downwards, which establishes natural circulation. Vijayan (2002) developed a correlation between two nondimensional groups (Ress ; Gr m =N G ) for single-phase steady state ⇑ Corresponding author. E-mail address: [email protected] (Y.M. Ferng). 0306-4549/$ - see front matter Ó 2013 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.anucene.2013.01.015
natural circulation loops with variable or constant flow cross-sections. Using one-dimensional (1-D) and 3-D CFD code, Pilkhwal et al. (2007) predicted the dynamic behavior related to a singlephase natural circulation loop. Misale et al. (2007) presented experimental investigations for a rectangular single-phase natural circulation mini-loop, which consists of two horizontal copper tubes (heat transfer sections) and two vertical tubes (legs). The effect of power transferred to the fluid and loop inclination had been systematically investigated. Vijayan et al. (2007) studied effects of heater and cooler orientations on the single-phase natural circulation in a rectangular loop. Depending on heating condition, a hysteresis regime in the stable or unstable flow was also observed. Kumar et al. (2011) studied single-phase natural circulation loops under natural and mixed convections. Effects of wall friction factor on loop dynamics had also been presented. Angelo et al. (2012) performed a 3-D numerical analysis of a steady state rectangular natural circulation loop located at the IPEN/CNEN-SP. Good agreement between steady-state numerical and measured results was shown in their work. In this paper, a compressible 3-D CFD model is proposed to investigate the thermal–hydraulic characteristics related to the natural circulation mechanism. The experiments conducted by Misale et al. (2007) are adopted to validate this CFD model, including the temperature transient behaviors and the relationship of Ress on the parameter of Gr m =N G , which demonstrates applicability of the CFD methodology on safety analysis of nuclear reactors with the passive safety features.
66
J.Y. Wang et al. / Annals of Nuclear Energy 58 (2013) 65–71
Nomenclature flow area, m2 specific heat, J/kg K loop diameter, m gravitational acceleration, m2/s generation of turbulence, kg/m s3 modified Grashof number, D3 q2 bDT r =l2 loop height, m turbulent kinetic energy, m2/s2 total length in a loop, m geometric parameter, Lt/D (Vijayan, 2002) pressure, N/m2 total heat source, W steady-state Reynolds number, DW/Al mass flowrate, kg/s time, s
A Cp D g G Grm H k Lt NG p Q Ress W t
T u x
fluid temperature, K fluid velocity, m/s x coordinate, m
Greek symbols fluid density, kg/m3 b temperature expansion coefficient, 1/K e turbulent dissipation rate, m2/s3 di;j delta function l viscosity, kg/m s lt turbulent viscosity, kg/m s s stress, N/m2 k thermal conductivity, W/m K rt turbulent Prandtl number (=0.85) DTr reference temperature difference, QH/AlCp, K
q
where
2. Simulation model Natural circulation essentially belongs to the density-driven phenomenon in a gravitational field, which is caused by the temperature difference. A compressible CFD model is adopted to consider the variation of density on the temperature, instead of the Boussinesq approximation. All of the simulation equations can be described as follows. 2.1. Governing equations
k
lt ¼ C l fl q P ¼ st;ij
e
@ui ; @xj
f1 ¼ 1; ¼ exp
2
mt ¼ lt =q
;
st;ij ¼ lt
ð9Þ
@ui @uj 2=3qkdij þ @xj @xi
f 2 ¼ 1 0:3 expðRe2t Þ; " # 3:4
ð10Þ
fl ð11Þ
ð1 þ Ret =50Þ2
Continuity equation:
@q @ ðqui Þ ¼ 0 þ @t @xi
ð1Þ
¼ 0 for i–j
Momentum equation:
@ðqui Þ @ @p @ sij þ ðqui uj Þ ¼ þ þ qg i @t @xj @xi @xj
dij ¼ 1 for i ¼ j
k
2
ð2Þ
Ret ¼
sij ¼ l
@ui @xj
ð3Þ
Dk ¼ 2m
l ¼ l for laminar flow ¼ l þ lt for turbulent flow
ð4Þ
Energy equation:
@ðqC p TÞ @ @ @T þ ðqC p ui TÞ ¼ k þ ui sij @t @xi @xi @xi
ð5Þ
where
k ¼ k
for laminar flow
¼ k þ Cplt =Prt
for turbulent flow
ð6Þ
In the present work, the low-Reynolds number k e model (Launder and Sharma, 1974) is adopted as the turbulence model if the natural circulation flow becomes the turbulent flow. The turbulent kinetic energy k, and energy dissipation rate e, can be obtained by solving the following transport equations:
@ðqkÞ @ l @k ¼ P qe qDk þ qkui l þ t @t @xi rk @xi
ð7Þ
ð13Þ
me
where
ð12Þ
Ee ¼ 2mmt
pffiffiffi!2 @ k @y
ð14Þ
!2 @2u @y2
ð15Þ
C e1 , C e2 , C l , rk and re are the model constants. The damping functions fl, f1 and f2 and the extra source terms D and E are only active close to solid walls, which makes it possible to solve k and e in the viscous sublayer. The values of these constants are listed in Table 1 (Launder and Sharma, 1974). 2.2. Boundary conditions As shown in Fig. 1a, a rectangular loop (Misale et al., 2007) is simulated in this paper. Its dimensions and the thermocouple locations are also indicated in this plot. A heater is located at the lower section and a cooler is put on the upper section. According to the Table 1 Constants in the low-Re k e turbulence model. C e1 C e2 Cl
@ðqeÞ @ l @e e ¼ ðC e1 f1 P C e2 f2 qeÞ þ qEe þ qeui l þ t @t @xi re @xi k ð8Þ
rk re
1.44 1.92 0.09 1.3 1.0
J.Y. Wang et al. / Annals of Nuclear Energy 58 (2013) 65–71
67
Fig. 1. Schematic of simulation loop (a) and typical mesh model (b).
experiments, the constant heat flux is imposed on the heating section and the constant temperature of 0 °C is maintained in the cooler, which are set as the thermal boundary conditions in the simulations. Adiabatic thermal and non-slip boundary conditions are set on the wall of the loop. In addition, the standard wall functions (Launder and Spalding, 1973) for both the momentum and the energy equations are applied near the wall to capture proper turbulent behaviors.
3. Numerical treatment The equations presented in the above section essentially belong to the non-linear partial differential equations (PDEs), which should be discretized into the finite differencing forms in order to perform numerical calculations. The pressure-based solver is adopted to the compressible flow simulations. The second-order upwind scheme is used to treat the convection terms in the PDEs. The pressure-implicit with splitting of operators (PISOs) scheme (Issa, 1986) is adopted to solve the velocity–pressure coupled equations. Existing in the body force term of momentum equation (Eq. (2)), the density at the mesh surface is treated using the second-order-upwind interpolation scheme. The second-order implicit formulation is used to treat the time-dependent term in the transient simulations. Fig. 1b shows the typical mesh distributions on the cross-section, straight line, and elbow, respectively. Detailed descriptions
about the experiments simulated herein can be referred to the work of Misale et al. (2007). As illustrated in Fig. 1b, 120 uniform nodes and 576 structured nodes are typically set along the axial direction and on the cross-section, respectively. Mesh independent calculations with a grid refinement ratio of 1.4 are performed. Three kinds of mesh distributions are adopted, including the typical nodes of 120 576, the finer nodes of 160 864, and the coarser nodes of 84 384. The maximum deviations in the values of Ress and Grm for the simulation case with input power of 25 W are less than 0.5%. These comparison results indicate that the simulation results with typical mesh are mesh-independent. In addition, the time step of 0.01 s is typically used in the transient simulations. Different values of time step are also adopted in the sensitivity calculations All the simulations presented in this paper are performed using the FLUENT code (ANSYS, 2009) on an Intel_Core i7 3.4G Processor PC. In every time step, the convergence criteria for all of the governing equations are set to be 105.
4. Results and discussion Fig. 2 compares the time histories of temperature at T4 point between the experimental measurements (blue1) (Misale et al., 2007) and the CFD predictions (red). In this simulation, the input 1 For interpretation of color in Figs. 2 and 3, the reader is referred to the web version of this article.
68
J.Y. Wang et al. / Annals of Nuclear Energy 58 (2013) 65–71
Fig. 2. Comparison of temperature histories at point T4 in the loop under input power of 25 kW.
heating power is 25 W. Similar trends of temperature histories between both results are revealed in this figure, except that more oscillations in the initial transient stage are predicted by the present CFD model. The maximum temperature in the T4 point is predicted to be 360 K that is higher than the measured value of 349 K by about 3%. The predicted occurrence time (36 s)at which the peak water temperature at the T4 point is reached is slightly earlier than the measured one (41 s). After the transient time exceeds about 200 s, the temperature at T4 point would approach a steady-state value, which is clearly demonstrated in both the measured and the predicted results. The steady-state temperature predicted by the present CFD model is 337 K slightly higher than the measured one of 336 K. The time histories of temperature difference (DT, DT = T4 T1) cross the heater in the lower section for the measurements (blue) and the predictions (red) are illustrated in Fig. 3. This figure clearly shows that the predicted characteristics of DT are similar to the measured ones. Detailed observation of Figs. 2 and 3 reveals that it would take about 200 s for the test loop to establish the steady-state natural circulation condition. Both experimental data
Fig. 3. Comparison of temperature difference (DT) histories across the heater under input power of 25 kW.
and model predictions confirm this result. The maximum DT in the transient period is predicted to be about 68 K, higher than the measured result of 45 K. However, the predicted steady-state DT is about 24.2 that is close to the measured one of 23.3. The DT is the essential force to drive the water in a natural circulation loop. This agreement implies that the steady-state natural circulation flowrate in the loop predicted by the present CFD model may reproduce the measured one. In order to compare with the steady-state measured data of Misale et al. (2007), six different power levels and three different values of loop inclination (a) are simulated. The effect of loop inclination can be simply considered in the gravitation acceleration vector for the present model. The power level and the loop inclination are listed in Table 2. Fig. 4 compares the steady-state temperature at T4 point for different input heating power levels between the CFD predictions and the experimental measurements. The maximum deviation for the predictions is less than 3%, as clearly shown in this figure. Comparison of the measured and the predicted DT across the heater is also illustrated in Fig. 5. The maximum deviation between the predictions and the measurements is less than 12%. Comparison of Figs. 4 and 5 reveals that the deviation of DT is larger than that of T4 temperature, implying that the present CFD model under-predicts the temperature at the heater inlet (i.e. T1 point). The underprediction may be caused by the ideal cooling condition set at the cooler in the CFD simulations. This ideal assumption resulting in the more cooling for the cooler and consequently the lower water temperature at the cooler outlet. Table 2 Power level and loop inclination (a) in the simulations. Power (W)
a (°)
5 15 25 2.5 7.5 10 2.5 7.5 10
0 0 0 30 30 30 75 75 75
Fig. 4. Comparison of steady-state temperature at T4 point in the loop under different values of input power.
69
J.Y. Wang et al. / Annals of Nuclear Energy 58 (2013) 65–71
vortex are also revealed in the secondary flow pattern on the cross-section of elbow. As the aforementioned description, the fluid flow in a natural circulation loop is essentially resulted from the buoyancy force generated by the temperature difference. Therefore, the amount of driven natural circulation flowrate is proportional to the magnitude of temperature difference. The general correlation for the steady-state flowrate and the reference temperature difference in a natural circulation loop can be expressed as (Vijayan and Austregesilo, 1994):
Ress ¼ C
r Grm NG
ð16Þ
In addition, based on the research results of Swapnalee and Vijayan (2011), this correlation for the laminar, transition and turbulent flow loops, respectively, has the following forms.
Fig. 5. Comparison of steady-state temperature difference (DT) across the heater under different values of input power.
Fig. 6 shows the predicted steady-state flow and heat transfer characteristics in a natural loop with the heating power of 25 kW. The shaded contours represent the calculated temperature and main stream vectors distributions on a 2-D plane near the heater. In order to clearly illustrate the thermal stratification and the secondary flow pattern, their distribution on the cross-section of elbow is enlarged and shown on the inner corner of elbow. As clearly revealed in Fig. 6, the thermal stratification is established in the heater portion, which is similar to the predicted results from Pilkhwal et al. (2007). The hotter fluid at the top region would be driven from the heater towards the left-side section, passes through the elbow, and flows upwards in the vertical pipe, resulting in the natural circulation flow. In addition, two symmetrical
0:5 Grm Ress ¼ 0:1768 NG
for a laminar flow loop
ð17Þ
0:387 Grm Ress ¼ 1:216 NG
for a transition flow loop
ð18Þ
0:364 Grm Ress ¼ 1:956 NG
for a turbulent flow loop
ð19Þ
Fig. 7 shows the relationship of Ress and Grm for the experimental measurements, the CFD predictions, and the correlation calculations. In addition to the experimental data from Misale et al. (2007), the comparison also includes the data from Kumar et al. (2011). In this figure, hollow symbols are the present CFD predicted results; lines are the calculation results from correlations of Eqs. (17) and (18). Two distinct regions: laminar region and transition region, respectively, are clearly revealed in Fig. 7. The temperature difference would increase with the increasing heating power, consequently inducing the higher circulation flowrate. This relationship can be precisely captured in the present CFD predictions. In addition, as the heating power is continuously increased, the induced flowrate is also increased, rendering that the highly turbulent circulation flow condition can be reached. However, in
Fig. 6. 2-D distributions of temperature contours and flow vectors on the axial plane of heater and the cross-section of elbow.
70
J.Y. Wang et al. / Annals of Nuclear Energy 58 (2013) 65–71
Fig. 7. Ress function of Grm/NG for experimental measurements, CFD model predictions, and correlation calculated results.
this condition, the localized boiling phenomenon may be occurred near the heater surface. This localized flow boiling condition is not considered in the present CFD model. Therefore, only the laminar and transition turbulent flow conditions are simulated in this study. It can be clearly seen in Fig. 7 that the present CFD predictions agree well with the measurements. Based on the present predicted results for a laminar flow loop, the relationship between Ress and Grm can be regressed as follows:
0:493 Grm Ress ¼ 0:1653 NG
ð20Þ
Compared with the correlation of Eq. (17) for a laminar flow loop, the deviation for the constant and the exponent in above equation is about 6.5% and 1.4%, respectively. In the transition region, the predicted relationship between Ress and Grm can be regressed as:
0:403 Grm Ress ¼ 0:9833 NG
ð21Þ
Deviation in the values of constant and exponent between the CFD regression formula (Eq. (21)) and the correlation (Eq. (18)) are about 19% and 4%, respectively. Larger discrepancy is revealed in the transition region. However, Eq. (21) regressed from the present CFD predictions is similar to that from the simulation results of Pilkhwal et al. (2007) in the region of Ress 103 . That is,
0:3962 Grm Ress ¼ 0:8422 NG
ð22Þ
These comparisons reveal that the present CFD model can +simulate the thermal–hydraulic characteristics in a single-phase natural circulation loop with enough accuracy.
(1) The thermal–hydraulic characteristics related to the natural circulation in a loop can be reasonably captured by the present CFD model. These phenomena include the flow driven by the density difference, the thermal stratification, and the secondary flow pattern, etc. (2) Compared with the experimental data from Misale et al. (2007), the predicted results from the present CFD model show good agreement. These measured data include the time histories of both the local temperature and the DT across the heater. The predicted steady-state value of maximum temperature in a natural circulation loop corresponds well with the measured one. However, some discrepancy between the measurements and the predictions is also revealed in the DT across the heater, which may be caused by the idea cooling condition set in the cooler for the simulations. (3) Using the relationship between Ress and Grm, the present CFD methodology with laminar model and low-Re k e turbulence model has been validated with the experimental data obtained from Misale et al. (2007) and Kumar et al. (2011). The comparisons reveal that the present predicted relationship between Ress and Grm agrees qualitatively and quantitatively with both the experimental data as well as the appropriate correlations proposed by Swapnalee and Vijayan (2011). The heat removal capability of natural circulation is strongly related to the flowrate. It is important for the design of natural circulation in the passive safety systems to predict the circulation flowrate. Therefore, based on the simulation and the comparison results, the present CFD model can be applied in simulating the natural circulation phenomenon in confidence.
References 5. Conclusions The main objective of this paper is to investigate the CFD capability in simulating the flow and heat transfer characteristics of natural circulation since the passive safety systems are generally implemented in new design of nuclear power plants and the natural circulation is the essential passive cooling mechanism. Several important conclusions can be drawn from the simulation results.
Angelo, G., Andradea, D.A., Angelo, E., Torresa, W.M., Sabundjiana, G., Macedoa, L.A., Silva, A.F., 2012. A numerical and three-dimensional analysis of steady state rectangular natural circulation loop. Nucl. Eng. Des.. http://dx.doi.org/10.1016/ j.nucengdes.2011.12.020. ANSYS Inc. 2009. Fluent R12 User’s Guide. Issa, R.I., 1986. Solution of the implicitly discretized fluid flow equation by operator splitting. J. Comput. Phys. 62, 40–65. Kumar, N., Doshi, J.B., Vijayan, P.K., 2011. Investigations on the role of mixed convection and wall friction factor in single-phase natural circulation loop dynamics. Ann. Nucl. Energy 38, 2247–2270.
J.Y. Wang et al. / Annals of Nuclear Energy 58 (2013) 65–71 Launder, B.E., Sharma, B.I., 1974. Application of the energy-dissipation model of turbulence to the calculation of flow near a spinning disc. Lett. Heat Mass Transfer 1, 131–138. Launder, B.E., Spalding, D.B., 1973. The numerical computational of turbulent flows. Comput. Meth. Appl. Mech. Eng. 3, 269. Misale, M., Garibaldi, P., Passos, J.C., Ghisi de Bitencourt, G., 2007. Experiments in a single-phase natural circulation mini-loop. Exp. Therm. Fluid Sci. 31, 1111– 1120. Pilkhwal, D.S., Ambrosini, W., Forgione, N., Vijayan, P.K., Saha, D., Ferreri, J.C., 2007. Analysis of the unstable behaviour of a single-phase natural circulation loop with one-dimensional and computational fluid-dynamic models. Ann. Nucl. Energy 34, 339–355.
71
Swapnalee, B.T., Vijayan, P.K., 2011. A generalized flow equation for single phase natural circulation loops obeying multiple friction factor. Int. J. Heat Mass Transfer 54, 2618–2629. Vijayan, P.K., 2002. Experimental observations on the general trends of the steady state and stability behaviour of single-phase natural circulation loops. Nucl. Eng. Des. 215, 139–152. Vijayan, P.K., Austregesilo, H., 1994. Scaling laws for single-phase natural circulation loops. Nucl. Eng. Des. 152, 331–347. Vijayan, P.K., Sharma, M., Saha, D., 2007. Steady state and stability characteristics of single-phase natural circulation in a rectangular loop with different heater and cooler orientations. Exp. Therm. Fluid Sci. 31, 925–945.
Contents lists available at SciVerse ScienceDirect
Annals of Nuclear Energy journal homepage: www.elsevier.com/locate/anucene
CFD investigating flow and heat transfer characteristics in a natural circulation loop J.Y. Wang, T.J. Chuang, Y.M. Ferng ⇑ Department of Engineering and System Science, Institute of Nuclear Engineering and Science, National Tsing Hua University, 101, Sec. 2 Kuang-Fu Rd., Hsingchu 30013, Taiwan, ROC
a r t i c l e
i n f o
Article history: Received 28 August 2012 Received in revised form 4 January 2013 Accepted 10 January 2013 Available online 3 April 2013 Keywords: Natural circulation Compressible CFD model Passive safety system
a b s t r a c t Natural circulation is the most important heat removal mechanism for passive safety systems in newdesign nuclear power plants. It is important to investigate the flow and heat transfer characteristics related to the natural circulation mechanism. In this paper, a compressible three-dimensional (3-D) CFD model is proposed to investigate these phenomena in a natural circulation loop. The flow and heat transfer behaviors in a loop can be reasonably captured by the present model, which includes the secondary flow within the elbow and the thermal stratification in the horizontal pipe, etc. This model is also validated with the existing experimental data. Compared with the measured results including time histories of local maximum fluid temperature and temperature difference across the heater, and the relationship between Ress and Grm/NG, the present predicted results show good agreement. These comparisons reveal that the present CFD methodology can be applied in simulating the thermal–hydraulic characteristics related to the natural circulation mechanism in confidence. Ó 2013 Elsevier Ltd. All rights reserved.
1. Introduction With advantage of simplification, economy, and no need of power supply, passive safety systems are generally adopted for new-design nuclear power plants. After occurrence of the Fusikama accident, natural circulation mechanism then plays a significant role in the safety related systems since it can be free of electrical demand. However, driving force for the natural circulation mechanism may be not strong enough to deliver the coolant and to remove the heat required for the accident. Therefore, it is crucial to develop analytical models to investigate the flow and heat transfer characteristics related to this heat removal mechanism. Regarding about the researches on natural circulation, watercooled loops had been generally adopted due to their importance in the industry applications, including nuclear power plants, solar water heaters, geothermal processes, machinery cooling, etc. In a natural circulation loop, water is driven by the buoyancy force that is essentially generated due to the temperature-driven density gradient. Heated in the lower portion of loop, water becomes less dense and moves upwards, and subsequently it can be cooled at the higher portion, becomes denser, and flows downwards, which establishes natural circulation. Vijayan (2002) developed a correlation between two nondimensional groups (Ress ; Gr m =N G ) for single-phase steady state ⇑ Corresponding author. E-mail address: [email protected] (Y.M. Ferng). 0306-4549/$ - see front matter Ó 2013 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.anucene.2013.01.015
natural circulation loops with variable or constant flow cross-sections. Using one-dimensional (1-D) and 3-D CFD code, Pilkhwal et al. (2007) predicted the dynamic behavior related to a singlephase natural circulation loop. Misale et al. (2007) presented experimental investigations for a rectangular single-phase natural circulation mini-loop, which consists of two horizontal copper tubes (heat transfer sections) and two vertical tubes (legs). The effect of power transferred to the fluid and loop inclination had been systematically investigated. Vijayan et al. (2007) studied effects of heater and cooler orientations on the single-phase natural circulation in a rectangular loop. Depending on heating condition, a hysteresis regime in the stable or unstable flow was also observed. Kumar et al. (2011) studied single-phase natural circulation loops under natural and mixed convections. Effects of wall friction factor on loop dynamics had also been presented. Angelo et al. (2012) performed a 3-D numerical analysis of a steady state rectangular natural circulation loop located at the IPEN/CNEN-SP. Good agreement between steady-state numerical and measured results was shown in their work. In this paper, a compressible 3-D CFD model is proposed to investigate the thermal–hydraulic characteristics related to the natural circulation mechanism. The experiments conducted by Misale et al. (2007) are adopted to validate this CFD model, including the temperature transient behaviors and the relationship of Ress on the parameter of Gr m =N G , which demonstrates applicability of the CFD methodology on safety analysis of nuclear reactors with the passive safety features.
66
J.Y. Wang et al. / Annals of Nuclear Energy 58 (2013) 65–71
Nomenclature flow area, m2 specific heat, J/kg K loop diameter, m gravitational acceleration, m2/s generation of turbulence, kg/m s3 modified Grashof number, D3 q2 bDT r =l2 loop height, m turbulent kinetic energy, m2/s2 total length in a loop, m geometric parameter, Lt/D (Vijayan, 2002) pressure, N/m2 total heat source, W steady-state Reynolds number, DW/Al mass flowrate, kg/s time, s
A Cp D g G Grm H k Lt NG p Q Ress W t
T u x
fluid temperature, K fluid velocity, m/s x coordinate, m
Greek symbols fluid density, kg/m3 b temperature expansion coefficient, 1/K e turbulent dissipation rate, m2/s3 di;j delta function l viscosity, kg/m s lt turbulent viscosity, kg/m s s stress, N/m2 k thermal conductivity, W/m K rt turbulent Prandtl number (=0.85) DTr reference temperature difference, QH/AlCp, K
q
where
2. Simulation model Natural circulation essentially belongs to the density-driven phenomenon in a gravitational field, which is caused by the temperature difference. A compressible CFD model is adopted to consider the variation of density on the temperature, instead of the Boussinesq approximation. All of the simulation equations can be described as follows. 2.1. Governing equations
k
lt ¼ C l fl q P ¼ st;ij
e
@ui ; @xj
f1 ¼ 1; ¼ exp
2
mt ¼ lt =q
;
st;ij ¼ lt
ð9Þ
@ui @uj 2=3qkdij þ @xj @xi
f 2 ¼ 1 0:3 expðRe2t Þ; " # 3:4
ð10Þ
fl ð11Þ
ð1 þ Ret =50Þ2
Continuity equation:
@q @ ðqui Þ ¼ 0 þ @t @xi
ð1Þ
¼ 0 for i–j
Momentum equation:
@ðqui Þ @ @p @ sij þ ðqui uj Þ ¼ þ þ qg i @t @xj @xi @xj
dij ¼ 1 for i ¼ j
k
2
ð2Þ
Ret ¼
sij ¼ l
@ui @xj
ð3Þ
Dk ¼ 2m
l ¼ l for laminar flow ¼ l þ lt for turbulent flow
ð4Þ
Energy equation:
@ðqC p TÞ @ @ @T þ ðqC p ui TÞ ¼ k þ ui sij @t @xi @xi @xi
ð5Þ
where
k ¼ k
for laminar flow
¼ k þ Cplt =Prt
for turbulent flow
ð6Þ
In the present work, the low-Reynolds number k e model (Launder and Sharma, 1974) is adopted as the turbulence model if the natural circulation flow becomes the turbulent flow. The turbulent kinetic energy k, and energy dissipation rate e, can be obtained by solving the following transport equations:
@ðqkÞ @ l @k ¼ P qe qDk þ qkui l þ t @t @xi rk @xi
ð7Þ
ð13Þ
me
where
ð12Þ
Ee ¼ 2mmt
pffiffiffi!2 @ k @y
ð14Þ
!2 @2u @y2
ð15Þ
C e1 , C e2 , C l , rk and re are the model constants. The damping functions fl, f1 and f2 and the extra source terms D and E are only active close to solid walls, which makes it possible to solve k and e in the viscous sublayer. The values of these constants are listed in Table 1 (Launder and Sharma, 1974). 2.2. Boundary conditions As shown in Fig. 1a, a rectangular loop (Misale et al., 2007) is simulated in this paper. Its dimensions and the thermocouple locations are also indicated in this plot. A heater is located at the lower section and a cooler is put on the upper section. According to the Table 1 Constants in the low-Re k e turbulence model. C e1 C e2 Cl
@ðqeÞ @ l @e e ¼ ðC e1 f1 P C e2 f2 qeÞ þ qEe þ qeui l þ t @t @xi re @xi k ð8Þ
rk re
1.44 1.92 0.09 1.3 1.0
J.Y. Wang et al. / Annals of Nuclear Energy 58 (2013) 65–71
67
Fig. 1. Schematic of simulation loop (a) and typical mesh model (b).
experiments, the constant heat flux is imposed on the heating section and the constant temperature of 0 °C is maintained in the cooler, which are set as the thermal boundary conditions in the simulations. Adiabatic thermal and non-slip boundary conditions are set on the wall of the loop. In addition, the standard wall functions (Launder and Spalding, 1973) for both the momentum and the energy equations are applied near the wall to capture proper turbulent behaviors.
3. Numerical treatment The equations presented in the above section essentially belong to the non-linear partial differential equations (PDEs), which should be discretized into the finite differencing forms in order to perform numerical calculations. The pressure-based solver is adopted to the compressible flow simulations. The second-order upwind scheme is used to treat the convection terms in the PDEs. The pressure-implicit with splitting of operators (PISOs) scheme (Issa, 1986) is adopted to solve the velocity–pressure coupled equations. Existing in the body force term of momentum equation (Eq. (2)), the density at the mesh surface is treated using the second-order-upwind interpolation scheme. The second-order implicit formulation is used to treat the time-dependent term in the transient simulations. Fig. 1b shows the typical mesh distributions on the cross-section, straight line, and elbow, respectively. Detailed descriptions
about the experiments simulated herein can be referred to the work of Misale et al. (2007). As illustrated in Fig. 1b, 120 uniform nodes and 576 structured nodes are typically set along the axial direction and on the cross-section, respectively. Mesh independent calculations with a grid refinement ratio of 1.4 are performed. Three kinds of mesh distributions are adopted, including the typical nodes of 120 576, the finer nodes of 160 864, and the coarser nodes of 84 384. The maximum deviations in the values of Ress and Grm for the simulation case with input power of 25 W are less than 0.5%. These comparison results indicate that the simulation results with typical mesh are mesh-independent. In addition, the time step of 0.01 s is typically used in the transient simulations. Different values of time step are also adopted in the sensitivity calculations All the simulations presented in this paper are performed using the FLUENT code (ANSYS, 2009) on an Intel_Core i7 3.4G Processor PC. In every time step, the convergence criteria for all of the governing equations are set to be 105.
4. Results and discussion Fig. 2 compares the time histories of temperature at T4 point between the experimental measurements (blue1) (Misale et al., 2007) and the CFD predictions (red). In this simulation, the input 1 For interpretation of color in Figs. 2 and 3, the reader is referred to the web version of this article.
68
J.Y. Wang et al. / Annals of Nuclear Energy 58 (2013) 65–71
Fig. 2. Comparison of temperature histories at point T4 in the loop under input power of 25 kW.
heating power is 25 W. Similar trends of temperature histories between both results are revealed in this figure, except that more oscillations in the initial transient stage are predicted by the present CFD model. The maximum temperature in the T4 point is predicted to be 360 K that is higher than the measured value of 349 K by about 3%. The predicted occurrence time (36 s)at which the peak water temperature at the T4 point is reached is slightly earlier than the measured one (41 s). After the transient time exceeds about 200 s, the temperature at T4 point would approach a steady-state value, which is clearly demonstrated in both the measured and the predicted results. The steady-state temperature predicted by the present CFD model is 337 K slightly higher than the measured one of 336 K. The time histories of temperature difference (DT, DT = T4 T1) cross the heater in the lower section for the measurements (blue) and the predictions (red) are illustrated in Fig. 3. This figure clearly shows that the predicted characteristics of DT are similar to the measured ones. Detailed observation of Figs. 2 and 3 reveals that it would take about 200 s for the test loop to establish the steady-state natural circulation condition. Both experimental data
Fig. 3. Comparison of temperature difference (DT) histories across the heater under input power of 25 kW.
and model predictions confirm this result. The maximum DT in the transient period is predicted to be about 68 K, higher than the measured result of 45 K. However, the predicted steady-state DT is about 24.2 that is close to the measured one of 23.3. The DT is the essential force to drive the water in a natural circulation loop. This agreement implies that the steady-state natural circulation flowrate in the loop predicted by the present CFD model may reproduce the measured one. In order to compare with the steady-state measured data of Misale et al. (2007), six different power levels and three different values of loop inclination (a) are simulated. The effect of loop inclination can be simply considered in the gravitation acceleration vector for the present model. The power level and the loop inclination are listed in Table 2. Fig. 4 compares the steady-state temperature at T4 point for different input heating power levels between the CFD predictions and the experimental measurements. The maximum deviation for the predictions is less than 3%, as clearly shown in this figure. Comparison of the measured and the predicted DT across the heater is also illustrated in Fig. 5. The maximum deviation between the predictions and the measurements is less than 12%. Comparison of Figs. 4 and 5 reveals that the deviation of DT is larger than that of T4 temperature, implying that the present CFD model under-predicts the temperature at the heater inlet (i.e. T1 point). The underprediction may be caused by the ideal cooling condition set at the cooler in the CFD simulations. This ideal assumption resulting in the more cooling for the cooler and consequently the lower water temperature at the cooler outlet. Table 2 Power level and loop inclination (a) in the simulations. Power (W)
a (°)
5 15 25 2.5 7.5 10 2.5 7.5 10
0 0 0 30 30 30 75 75 75
Fig. 4. Comparison of steady-state temperature at T4 point in the loop under different values of input power.
69
J.Y. Wang et al. / Annals of Nuclear Energy 58 (2013) 65–71
vortex are also revealed in the secondary flow pattern on the cross-section of elbow. As the aforementioned description, the fluid flow in a natural circulation loop is essentially resulted from the buoyancy force generated by the temperature difference. Therefore, the amount of driven natural circulation flowrate is proportional to the magnitude of temperature difference. The general correlation for the steady-state flowrate and the reference temperature difference in a natural circulation loop can be expressed as (Vijayan and Austregesilo, 1994):
Ress ¼ C
r Grm NG
ð16Þ
In addition, based on the research results of Swapnalee and Vijayan (2011), this correlation for the laminar, transition and turbulent flow loops, respectively, has the following forms.
Fig. 5. Comparison of steady-state temperature difference (DT) across the heater under different values of input power.
Fig. 6 shows the predicted steady-state flow and heat transfer characteristics in a natural loop with the heating power of 25 kW. The shaded contours represent the calculated temperature and main stream vectors distributions on a 2-D plane near the heater. In order to clearly illustrate the thermal stratification and the secondary flow pattern, their distribution on the cross-section of elbow is enlarged and shown on the inner corner of elbow. As clearly revealed in Fig. 6, the thermal stratification is established in the heater portion, which is similar to the predicted results from Pilkhwal et al. (2007). The hotter fluid at the top region would be driven from the heater towards the left-side section, passes through the elbow, and flows upwards in the vertical pipe, resulting in the natural circulation flow. In addition, two symmetrical
0:5 Grm Ress ¼ 0:1768 NG
for a laminar flow loop
ð17Þ
0:387 Grm Ress ¼ 1:216 NG
for a transition flow loop
ð18Þ
0:364 Grm Ress ¼ 1:956 NG
for a turbulent flow loop
ð19Þ
Fig. 7 shows the relationship of Ress and Grm for the experimental measurements, the CFD predictions, and the correlation calculations. In addition to the experimental data from Misale et al. (2007), the comparison also includes the data from Kumar et al. (2011). In this figure, hollow symbols are the present CFD predicted results; lines are the calculation results from correlations of Eqs. (17) and (18). Two distinct regions: laminar region and transition region, respectively, are clearly revealed in Fig. 7. The temperature difference would increase with the increasing heating power, consequently inducing the higher circulation flowrate. This relationship can be precisely captured in the present CFD predictions. In addition, as the heating power is continuously increased, the induced flowrate is also increased, rendering that the highly turbulent circulation flow condition can be reached. However, in
Fig. 6. 2-D distributions of temperature contours and flow vectors on the axial plane of heater and the cross-section of elbow.
70
J.Y. Wang et al. / Annals of Nuclear Energy 58 (2013) 65–71
Fig. 7. Ress function of Grm/NG for experimental measurements, CFD model predictions, and correlation calculated results.
this condition, the localized boiling phenomenon may be occurred near the heater surface. This localized flow boiling condition is not considered in the present CFD model. Therefore, only the laminar and transition turbulent flow conditions are simulated in this study. It can be clearly seen in Fig. 7 that the present CFD predictions agree well with the measurements. Based on the present predicted results for a laminar flow loop, the relationship between Ress and Grm can be regressed as follows:
0:493 Grm Ress ¼ 0:1653 NG
ð20Þ
Compared with the correlation of Eq. (17) for a laminar flow loop, the deviation for the constant and the exponent in above equation is about 6.5% and 1.4%, respectively. In the transition region, the predicted relationship between Ress and Grm can be regressed as:
0:403 Grm Ress ¼ 0:9833 NG
ð21Þ
Deviation in the values of constant and exponent between the CFD regression formula (Eq. (21)) and the correlation (Eq. (18)) are about 19% and 4%, respectively. Larger discrepancy is revealed in the transition region. However, Eq. (21) regressed from the present CFD predictions is similar to that from the simulation results of Pilkhwal et al. (2007) in the region of Ress 103 . That is,
0:3962 Grm Ress ¼ 0:8422 NG
ð22Þ
These comparisons reveal that the present CFD model can +simulate the thermal–hydraulic characteristics in a single-phase natural circulation loop with enough accuracy.
(1) The thermal–hydraulic characteristics related to the natural circulation in a loop can be reasonably captured by the present CFD model. These phenomena include the flow driven by the density difference, the thermal stratification, and the secondary flow pattern, etc. (2) Compared with the experimental data from Misale et al. (2007), the predicted results from the present CFD model show good agreement. These measured data include the time histories of both the local temperature and the DT across the heater. The predicted steady-state value of maximum temperature in a natural circulation loop corresponds well with the measured one. However, some discrepancy between the measurements and the predictions is also revealed in the DT across the heater, which may be caused by the idea cooling condition set in the cooler for the simulations. (3) Using the relationship between Ress and Grm, the present CFD methodology with laminar model and low-Re k e turbulence model has been validated with the experimental data obtained from Misale et al. (2007) and Kumar et al. (2011). The comparisons reveal that the present predicted relationship between Ress and Grm agrees qualitatively and quantitatively with both the experimental data as well as the appropriate correlations proposed by Swapnalee and Vijayan (2011). The heat removal capability of natural circulation is strongly related to the flowrate. It is important for the design of natural circulation in the passive safety systems to predict the circulation flowrate. Therefore, based on the simulation and the comparison results, the present CFD model can be applied in simulating the natural circulation phenomenon in confidence.
References 5. Conclusions The main objective of this paper is to investigate the CFD capability in simulating the flow and heat transfer characteristics of natural circulation since the passive safety systems are generally implemented in new design of nuclear power plants and the natural circulation is the essential passive cooling mechanism. Several important conclusions can be drawn from the simulation results.
Angelo, G., Andradea, D.A., Angelo, E., Torresa, W.M., Sabundjiana, G., Macedoa, L.A., Silva, A.F., 2012. A numerical and three-dimensional analysis of steady state rectangular natural circulation loop. Nucl. Eng. Des.. http://dx.doi.org/10.1016/ j.nucengdes.2011.12.020. ANSYS Inc. 2009. Fluent R12 User’s Guide. Issa, R.I., 1986. Solution of the implicitly discretized fluid flow equation by operator splitting. J. Comput. Phys. 62, 40–65. Kumar, N., Doshi, J.B., Vijayan, P.K., 2011. Investigations on the role of mixed convection and wall friction factor in single-phase natural circulation loop dynamics. Ann. Nucl. Energy 38, 2247–2270.
J.Y. Wang et al. / Annals of Nuclear Energy 58 (2013) 65–71 Launder, B.E., Sharma, B.I., 1974. Application of the energy-dissipation model of turbulence to the calculation of flow near a spinning disc. Lett. Heat Mass Transfer 1, 131–138. Launder, B.E., Spalding, D.B., 1973. The numerical computational of turbulent flows. Comput. Meth. Appl. Mech. Eng. 3, 269. Misale, M., Garibaldi, P., Passos, J.C., Ghisi de Bitencourt, G., 2007. Experiments in a single-phase natural circulation mini-loop. Exp. Therm. Fluid Sci. 31, 1111– 1120. Pilkhwal, D.S., Ambrosini, W., Forgione, N., Vijayan, P.K., Saha, D., Ferreri, J.C., 2007. Analysis of the unstable behaviour of a single-phase natural circulation loop with one-dimensional and computational fluid-dynamic models. Ann. Nucl. Energy 34, 339–355.
71
Swapnalee, B.T., Vijayan, P.K., 2011. A generalized flow equation for single phase natural circulation loops obeying multiple friction factor. Int. J. Heat Mass Transfer 54, 2618–2629. Vijayan, P.K., 2002. Experimental observations on the general trends of the steady state and stability behaviour of single-phase natural circulation loops. Nucl. Eng. Des. 215, 139–152. Vijayan, P.K., Austregesilo, H., 1994. Scaling laws for single-phase natural circulation loops. Nucl. Eng. Des. 152, 331–347. Vijayan, P.K., Sharma, M., Saha, D., 2007. Steady state and stability characteristics of single-phase natural circulation in a rectangular loop with different heater and cooler orientations. Exp. Therm. Fluid Sci. 31, 925–945.