- Email: [email protected]
A numerical study of boiling flow instability of a reactor thermosyphon system
Applied Thermal Engineering 26 (2006) 644–653 www.elsevier.com/locate/apthermeng
A numerical study of boiling flow instability of a reactor thermosyphon system A.K. Nayak a
a,*
, D. Lathouwers a, T.H.J.J. van der Hagen a, Frans Schrauwen b, Peter Molenaar b, Andrew Rogers b
Interfaculty Reactor Institute, Delft University of Technology, Mekelweg 15, 2629 JB Delft, The Netherlands b Shell Research and Technology Centre, Badhuisweg 3, 1031 CM Amsterdam, The Netherlands Received 27 July 2004; accepted 22 May 2005 Available online 8 November 2005
Abstract A numerical study has been carried out to investigate the boiling flow instability of a reactor thermosyphon system. The numerical model solves the conservation equations of mass, momentum and energy applicable to a two-fluid and three-field steam–water system using a finite difference technique. The computer code MONA was used for this purpose. The code was applied to the thermosyphon system of an EO (ethylene oxide) chemical reactor in which the heat released by a catalytic reaction is carried by boiling water under natural circulation conditions. The steady-state characteristics of the reactor thermosyphon system were predicted using the MONA code and conventional two-phase flow models in order to understand the model applicability for this type of thermosyphon system. The two-fluid model was found to predict the flow closest to the measured value of the plant. The stability behaviour of the thermosyphon system was investigated for a wide range of operating conditions. The effects of power, subcooling, riser length and riser diameter on the boiling flow instability were determined. The system was found to be unstable at higher power conditions which is typical for a Type II instability. However, with an increase in riser diameter, oscillations at low power were observed as well. These are classified as Type I instabilities. Stability maps were predicted for both Type I and Type II instabilities. Methods of improving the stability of the system are discussed. Ó 2005 Published by Elsevier Ltd. Keywords: Natural circulation; Chemical; Reactor; Two-phase flow; Stability; Type I instability; Type II instability
1. Introduction Use of a thermosyphon is the preferred option for heat removal of many industrial systems, such as nuclear reactors [1] and chemical reactors [2], since it eliminates problems associated with pump maintenance and operation. However, it is well known that this mode of cooling may suffer from instabilities which can hinder the operation and safety of plants. Many types of insta* Corresponding author. Address: Reactor Engineering Division, Bhabha Atomic Research Centre, Trombay, Mumbai 400 085, India. E-mail addresses: nayak_arun@rediffmail.com, arunths@apsara. barc.ernet.in (A.K. Nayak).
1359-4311/$ - see front matter Ó 2005 Published by Elsevier Ltd. doi:10.1016/j.applthermaleng.2005.05.019
bilities are known to occur in boiling two-phase systems; some of which depend on the static or steady-state characteristics of the system, whereas others depend on the dynamic characteristics of the system [3,4]. A densitywave instability is the typical dynamic instability which may occur due to the multiple regenerative feedback between the flow rate, enthalpy, density and pressure drop in the boiling system. The occurrence of the instability depends on the perturbed pressure drop in the two-phase and single-phase regions of the system and the propagation time delay of the void fraction or density in the system. Such an instability can occur at very low power and at high power conditions. This depends on the relative importance of the respective components
A.K. Nayak et al. / Applied Thermal Engineering 26 (2006) 644–653
645
Nomenclature A Cp D f g h hfg K L
cross-sectional area, m2 specific heat, J/kg K hydraulic diameter, m friction coefficient acceleration due to gravity, m/s2 specific enthalpy, J/kg latent heat of vapourisation, J/kg loss coefficient length, m
of pressure drop such as gravity or frictional losses in the system. Fukuda and Kobori [5] have classified the density-wave instability as Type I and Type II for the low power and high power instabilities respectively. The mechanisms can be explained as follows. 1.1. Type I instability For this type of instability to occur, the presence of a long riser plays an important role such as in a boiling two-phase natural circulation loop. Under low quality conditions, a slight change in quality due to any disturbance can cause a large change in void fraction (Fig. 1) and consequently in the driving head. Therefore the flow can oscillate at such low power conditions. But as the power increases, the flow quality increases where the slope of the void fraction vs. quality reduces. This can suppress the fluctuation of the driving head for a small change in quality. Hence, the flow stabilises at higher power. 1.2. Type II instability Unlike the Type I instability, the Type II instability occurs at high power conditions. This instability is driven by the interaction between the single and two-phase frictional component of pressure losses, mass flow, void
α
x Fig. 1. Typical variation of void fraction with quality.
p Q T t w x z a q
pressure, N/m2 volumetric heat addition rate, W/m3 temperature, K time, s mass flow rate, kg/s quality axial distance, m vapour volume fraction density, kg/m3
formation and propagation in the two-phase region. At high power, the flow quality or void fraction in the system is very large. Hence, the two-phase frictional pressure loss may be high owing to the smaller two-phase mixture density. Having a large void fraction will increase the void propagation time delay in the twophase region of the system. Under these conditions, any small fluctuation in flow can cause a larger fluctuation of the two-phase frictional pressure loss due to fluctuation of density and flow, which propagates slowly in the twophase region. On the other hand, the fluctuation of the pressure drop in the single-phase region occurs due to fluctuation of flow alone since the fluctuation of the density is negligible. The pressure drop fluctuation in this region travels much faster due to incompressibility of single-phase region. If the two-phase pressure drop fluctuation is equal in magnitude but opposite in phase with that of the single-phase region, the fluctuation or oscillation is sustained in the system since there are no attenuating mechanisms. Divergent oscillations can occur depending on the magnitude of the pressure loss fluctuation in the two-phase and single-phase regions and the propagation time delay. Both Type I and Type II instabilities can occur in a boiling thermosyphon system depending on its geometry and operating conditions. A literature review suggests that there have been many studies on this phenomenon in the past in two-phase natural circulation systems [6– 11]. Many numerical studies have also been performed in the past to model the density-wave instability [12,13]. However, there are no studies on boiling twophase thermosyphon phenomena inside the vessel of a chemical reactor such as an EO reactor. In these reactors, the heat released by chemical reaction in the catalyst tubes, is carried by boiling steam and water mixture under natural circulation conditions. Flow instability is highly undesirable, since oscillations of flow can cause oscillations of the catalyst temperature which can induce oscillations in heat generation rate in the catalysts. In that case, it is impossible to run the reactor for safety reason. Hence, it is essential to investigate the flow
646
A.K. Nayak et al. / Applied Thermal Engineering 26 (2006) 644–653
stability behaviour of such reactors in order to find out methods of improving the stability characteristics of operating reactor systems and scope for upscaling of future reactors. In EO reactors, the diameter of the vessel in which heat is generated due to a chemical reaction from a large number of catalyst tubes is very large (typically 5 m). Hence, the natural circulation velocity in the reactor is very small. Under these conditions, flow stratification can occur in the vessel. With spacers present to support the catalyst tubes, the two-phase flow can have a very complex distribution. To solve the stability problem of this system, use of a simple homogeneous model or a drift flux model as was mostly used by previous investigators, cannot give accurate results since these models do not consider flow pattern transitions and use empirical correlations for the friction factor and void fraction distribution which may not be applicable for all conditions. Further, treating multi-dimensional two-phase flows using numerical codes is complex. So for studying this phenomenon, the flow is assumed to be one-dimensional. It is required to predict the natural circulation stability characteristics using a robust two-fluid model that considers the transport equations of each phase separately. The objective of this study is to analyse the stability of boiling two-phase thermosyphon phenomena inside an EO (ethylene oxide) chemical reactor using a numerical two-fluid and three-field model. The model comprises conservation equations of mass, momentum and energy applicable to a steam–water flow and is numerically solved using the MONA code. First, the code was applied for simulation of the steady-state natural circulation behaviour of a chemical reactor for which a flow rate measurement at operating power was available. In order to understand the applicability of conventional two-phase flow models for void fraction and friction factor multiplier for this type of thermosyphon system, a theoretical model was developed. The model comprises simplified conservation equations of mass, momentum and energy for steady-state conditions. Different void fraction models are considered to calculate the density in the two-phase region and different two-phase friction factor multiplier models are used to calculate the frictional pressure loss in the two-phase region. The MONA code predictions were compared with that of the simplified model and the measured plant data. The limitations of the conventional models for application to such systems were determined. The stability behaviour of the system was analysed considering the effects of power and subcooling as parameters. The effects of riser diameter and riser length on flow stability were simulated. Stability maps were predicted for different operating conditions and geometries of the reactor. Methods of improving the stability of the system are discussed.
2. Computer code MONA The computer code MONA [14] simulates the flow of steam–water/inert vapour in any pipe line assuming one-dimensional flow. The code is basically an updated version of the RAMONA code used for nuclear reactor thermal hydraulic analysis. MONA contains a set of seven conservation equations, based on the modelling of three-flow fields; the bulk liquid or the liquid film at the wall, the liquid droplets and the vapour or steam. The code solves three mass conservation equations, two momentum conservation equations and two energy conservation equations. Closure relations provide expressions for two-phase friction and comprise a complete set of correlations for heat transfer at the wall and between the respective phases. A phenomenological approach is employed in MONA to model the different flow regimes and the transition between them, thus no separate user specified correlations for void distribution, etc. are required. The exchange of energy and momentum at the interface between vapour and liquid depends on the interfacial area and the topology of the two-phase flow. So the flow regime determines the wetted perimeters and friction factors. The interfacial heat transfer rates are estimated from the interfacial heat transfer coefficient, interfacial area and the temperature difference between the phases. The heat transfer coefficient is determined from using suitable correlations depending on the flow pattern. The MONA code uses a staggered mesh scheme in which the flow variables (velocities, mass flows, etc.) are defined on section boundaries and scalar variables (pressure, void fraction, densities, etc.) are located in the centre of section volumes. Donor cell differencing is adopted. For ease in solution, the mass conservation equations are decoupled from the momentum, energy and volume equations. The time step is limited by the Courant criterion based on the mass propagation velocity. The code has been extensively validated against data for various two-phase thermal hydraulic conditions including pressure wave propagation, U-tube flow oscillations, void fraction, pressure drop, etc. [15,16].
3. Theoretical model for steady-state natural circulation In order to understand the applicability of conventional two-phase flow models for void fraction and friction factor multiplier for the prediction of steady-state natural circulation in a chemical reactor, a theoretical model was developed. A brief description of the model is given in Appendix A. The model predicts the flow rate in a natural circulation loop for a given operating power, subcooling and pressure. The model incorporates various two-phase friction factor multiplier models and void fraction models which can influence on the natural circulation characteristics of the loop.
A.K. Nayak et al. / Applied Thermal Engineering 26 (2006) 644–653
used as described below in order to understand the model applicability for simulation of this type of thermosyphon system:
4. Geometry of the reactor thermosyphon system In the present analysis, we have considered the thermosyphon behaviour inside the vessel of an EO chemical reactor [17]. This is because the measured natural circulation flow rate of the plant at the operating power is available and it is useful for validation of the models. Fig. 2 shows a schematic outline of the reactor considered. The reactor vessel diameter is very large (5 m) which houses many closed tubes containing the EO catalyst. The height of the vessel is about 10 m. The diameter of riser pipes is of the order of 30 mm and the total height of the system is about 15 m. There are 10 spacers in the vessel to hold the catalyst tubes in position and the spacer loss coefficient was considered to be about 3.5. The local losses in the riser section was considered to be about 5.75. The heat released by a chemical reaction is carried by a steam–water mixture through natural circulation. The steam pressure is maintained at about 4.74 MPa in order to maintain the required temperature of the catalyst. Subcooled water from the downcomer enters the vessel and is heated as the fluid flows upward due to buoyancy. Boiling takes place and the steam and water mixture moves to a steam drum through small size riser pipes, wherein the steam and water separation takes place by gravitational settling. The steam is sent to a cooler and an equal amount of feedwater at lower temperature joins the saturated water in the steam drum.
(i) homogeneous model for void fraction and twophase friction factor multiplier; (ii) homogeneous model for void fraction and Martinelli–Nelson model for two-phase friction factor multiplier; (iii) Zuber–Findlay model for void fraction and Martinelli–Nelson model for two-phase friction factor multiplier, and (iv) two-fluid/three-field model as implemented in the MONA code. The results show that (Fig. 3) even though the trend for the variation of the natural circulation flow rate with power is similar, the magnitude of the flow rate as predicted by the various models is quite different. With increasing power, all the models predict that the flow rate initially rises due to an increase in buoyancy force and then reduces due to an increase in two-phase resistance at higher void fraction. The homogeneous model gives the largest flow rate since this model predicts a higher void fraction and smaller two-phase frictional pressure loss as compared to the other models at given power. Of course, when considering slip, the flow rate reduces, but the flow rate was still found to be higher than the measured natural circulation flow rate of the reactor at 18 MW power and subcooling of 2 K. The two-fluid model predicts a flow rate of about 183 kg/s at the operating power, which is closest to the measured value of the reactor at 18 MW power. The reason behind this is that the two-fluid model predicts a higher pressure drop in the two-phase region as compared to the other models.
5. Steady-state characteristics of the system The steady-state behaviour of the reactor thermosyphon system as described above, was predicted at different power keeping the steam drum pressure and downcomer subcooling constant. Several models were
Steam Steam drum Feed water Riser pipes
Reactor vessel
647
Steam and water mixture
Catalyst assemblies
downcomers
Steam + Water
Water
Fig. 2. Schematic of the thermosyphon system considered in the present analysis.
648
A.K. Nayak et al. / Applied Thermal Engineering 26 (2006) 644–653 350 homogenous void fraction model + Martinelli-Nelson Two-phase Multiplier model
homogeneous model
300
flow (kg/s)
250
Zuber-Findlay void-fraction model + Martinelli-Nelson Two-phase multiplier model
200
two-fluid model
150
subcooling = 2K
measurement
pressure = 4.74 MPa
100 0
10
20
30
40
50
60
70
80
power (MW) Fig. 3. Comparison of steady-state flow rate using different models.
Table 1 Comparison of pressure drops in different sections of the thermosyphon as predicted by different models Component
HEM
HEM for void fraction + M–N for two-phase multiplier
Zuber–Findlay for void fraction + M–N for two-phase multiplier
Two-fluid model
Heated vessel Horizontal riser Vertical riser Steam drum Downcomer (vertical)
0.532 bar 0.0127 bar 0.279 bar 0.077 bar 1.129 bar
0.533 bar 0.0343 bar 0.28 bar 0.077 bar 1.129 bar
0.72 bar 0.0345 bar 0.288 bar 0.077 bar 1.129 bar
0.78 bar 0.08 bar 0.346 bar 0.077 bar 1.129 bar
Table 1 shows a comparison of the pressure loss in various components at a flow rate of 183 kg/s for a power of 18 MW and a subcooling of 2 K. The difference in pressure drop among different models is found to be in the heated vessel and in the horizontal piping mostly due to the occurrence of flow stratification in the large diameter vessel and the horizontal riser pipe.
6. Stability characteristics of the system The stability characteristics of the reactor were predicted using the MONA code for different operating conditions and geometric parameters in order to understand their influence on stability. 6.1. Effects of operating conditions on stability Fig. 4 shows the influence of operating power on the stability of the system using the two-fluid model. The thermosyphon flow behaviour at different power levels was predicted at an operating pressure of 4.74 MPa and a reactor inlet subcooling of 3 K. For all the conditions, the computations were started with an assumed
initial flow rate of 283 kg/s, zero void fraction in the heated and adiabatic regions and a uniform pressure (4.74 MPa) throughout the loop. A sensitivity analysis was done to investigate the effects of node sizes and time step on thermosyphon stability phenomena. The code predictions were found to be insensitive to this variation. At low power, the flow is found to be stable. With a rise in power from 10 MW to 18 MW, the thermosyphon flow is found to increase. When the power is further increased to 75 MW, small amplitude flow oscillations appear which die out as time proceeds. However, the corresponding flow rate is found to be lower than that at 18 MW. When the power is further increased to about 95 MW, a limit cycle characterised by constant amplitude and frequency of flow oscillations, is observed. Similar behaviour is also observed at a subcooling of 5 K as shown in Fig. 5. In this case, the computations were started with the same assumed conditions as in the previous case (Fig. 4). However, the starting transients up to 600 s are not shown here. The flow is found to be stable up to a power of 80 MW (Fig. 5(a)). With a further increase in power, the flow becomes unstable as expected (Fig. 5(b)). The frequency and amplitude of oscillations were found to increase with increasing
A.K. Nayak et al. / Applied Thermal Engineering 26 (2006) 644–653
649
450 400 350
flow rate (kg/s)
300
subcooling = 3 K pressure = 4.74 MPa
250
18 MW
10 MW
200 150 100
95 MW
75 MW
50 0 0
200
400
600
800
1000
1200
time (s) Fig. 4. Effect of power rise on natural circulation behaviour at 3 K subcooling.
400 350 300
flow (kg/s)
250 75 MW 200
80 MW
150 100 subcooling = 5 K pressure = 4.74 MPa
2.5 MW 50 600
800
(a)
1000
1200
time (s) 200 85 MW
subcooling = 5 K
flow (kg/s)
150
100
95 MW
90 MW
105 MW
50 600
(b)
800
1000
1200
time (s)
Fig. 5. Effect of power rise on natural circulation behaviour at 5 K subcooling.
650
A.K. Nayak et al. / Applied Thermal Engineering 26 (2006) 644–653 250 power = 95 MW pressure = 4.74 MPa subcooling = 10 K
200
flow (kg/s)
subcooling = 5 K
150
100 subcooling = 15 K subcooling = 25 K
subcooling = 3 K
50 600
700
800
900
1000
time (s)
Fig. 6. Effect of subcooling on natural circulation behaviour.
power. With a rise in power, the flow rate also increases and hence, the time period of fluid to pass through the loop reduces. As a result, the frequency of oscillations which depends on the circulation period, increases. The increase in amplitude is due to the larger flow-void feedback which increases with a rise in power or void fraction in the system. The effect of reactor inlet subcooling on the flow stability is shown in Fig. 6. The operating power was fixed at 95 MW and the pressure at 4.74 MPa. The subcooling was increased from 3 K to 30 K. At low subcooling the flow was found to be unstable and stabilised above a subcooling of 15 K. Fig. 7 shows the stability map plotted in the powersubcooling plane. Since the flow reduces with an increase in power, the indication of increase of threshold power for instability with subcooling, clearly indicates that the type of instability observed in this system is 140 120 Unstable
power (MW)
100 80 60
Stable
40 20 0 0
10
20
30
subcooling (K) Fig. 7. Stability map of the system.
40
50
the so-called Type II instability. This type of instability occurs due to the dominance of frictional losses occurring at higher power when the natural circulation flow rate and two-phase mixture density is small. So, in order to have larger stability margin for the operating system, it is recommended to operate at higher subcooling. 6.2. Effect of riser geometry on stability The riser length and area are two important parameters that play a significant role in the natural circulation characteristics of a system. An increase in riser length or a decrease in riser area increases the two-phase resistance which can reduce the natural circulation flow rate and increases the void fraction for the same heating power. This may affect the performance of the reactor because of low flow and high void fraction conditions which favour Type II density-wave instabilities as discussed earlier. On the other hand, a decrease in riser length or an increase of flow area, even though increases the natural circulation flow rate, reduces the flow quality or void fraction for the same operating power. Under these conditions, slug flow may occur in the riser pipes which can cause undesirable flow induced vibrations. Besides, operation under such low quality conditions may induce the low quality Type I instability. Hence, it is essential to investigate the effects of the riser geometry on the natural circulation characteristics of the reactor. This study is important not only for finding the scope for improving the stability of current reactors but also in finding the effects of increasing the reactor size (up-scaling) on the stability. In this work, the horizontal riser length and the riser area were varied keeping all other parameters constant and their effects on stability were predicted using the MONA code.
A.K. Nayak et al. / Applied Thermal Engineering 26 (2006) 644–653
140
unstable
651
250
length = 5.7 m
1.2 Power = 10 MW subcooling = 2 K
length = 10.7 m
200
length = 15.7 m
80
Area = 0.1459 m2 Area = 0.2918 m2 Area = 0.4377 m2
stable 60
0.8
0.6 150
40
0.4
pressure = 4.74 MPa
20 0 0
10
20
void fraction
1.0
100
flow (kg/s)
threshold power (MW)
120
30
40
100
50
200
400
600
subcooling (K)
800
1000
1200
0.2 1600
1400
time (s)
Fig. 8. Effect of the horizontal riser length on stability.
Fig. 10. Effect of an increase in riser area on stability.
Fig. 8 shows the effect of riser length on flow stability. The length of the horizontal riser was varied from 5.7 m to 15.7 m keeping its diameter constant. All other dimensions of the systems were unchanged as well. A decrease in riser length was found to be stabilising. With decreasing riser length, the two-phase resistance is reduced which stabilises the two-phase flow. Fig. 9 shows the effect of riser flow area on stability. In one case, the actual two riser pipes of flow area 0.0729 m2 each was considered. In the second case, the two riser pipes were combined into one. So even though the total flow area remained the same, the hydraulic diameter increased. This reduces the frictional resistances and hence stabilises the flow. If the riser area is increased further, instabilities can appear at low power condition as seen in Fig. 10. For smaller riser area (0.1459 m2), the flow rate is smaller due to larger resistance in small riser pipes. As the flow area is increased, the flow rate increases, which gives rise to small frequency oscillations, typical of low quality Type I density-wave instability [5]. Such an instability occurs in natural circulation systems having tall risers. Under low power conditions, a slight change in quality
can cause a large change in void fraction or flow rate to induce this instability in the system. The corresponding void fraction oscillations are also shown in Fig. 10. With an increase in riser area, the time period of oscillation reduces due to the increase in flow rate in the system. Fig. 11 shows the instability characteristics of the system for a riser area of 0.2918 m2 (double the nominal 0.5 300
0.4
flow (kg/s)
0.2 200
0.1 0.0
150
flow (5 MW ) flow (10 MW) flow (15 MW)
void fraction (5 MW) void fraction (10 MW) void fraction (15 MW) 400
600
-0.2
(a)
800
1000
1200
1400
time (s) 350
1.0
unstable
140
one pipe - flow area 0.1459 m2
120
void fraction (120 MW) void fraction (140 MW)
80 two pipes - flow area of each 0.0729 m2 60 stable
40
0.8
300
flow rate (120 MW) flow rate (140 MW)
275
void fraction
325
100 flow (kg/s)
threshold power (MW)
-0.1
100 200
riser void fraction
0.3 250
0.6
20 250 200
0 0
10
20
30
40
50
(b)
400
600
800
1000
time (s)
subcooling (K) Fig. 9. Effect of the riser area on stability.
Fig. 11. Effect of an increase in power on stability for a riser area of 0.2918 m2.
652
A.K. Nayak et al. / Applied Thermal Engineering 26 (2006) 644–653
flow models for analyzing the natural circulation behaviour of such reactor systems were verified. The main conclusions of the study are summarised below.
70
threshold power (MW)
60 Stable 50 subcooling = 2 K
40 30
Unstable 20 10 1.5
2.0
2.5
3.0
3.5
4.0
4.5
5.0
5.5
riser area/nominal area Fig. 12. Stability boundary for Type I instability.
area) for different powers. Unlike for the nominal case as seen before (Figs. 4 and 5), with increasing riser area, only low power instability occurs in the system. The high power instability disappears as the riser area is increased. Similar behaviour is observed with further increase in riser area also. Type I instability maps as observed in this case are plotted in Fig. 12 for the threshold power of instability as a function of riser area at a subcooling of 2 K. It is interesting to note that the threshold power for the flow to become stable increases as the riser area is increased.
7. Conclusions This study has focussed on the prediction of the flow stability behaviour in natural circulation cooled chemical reactors, which are highly undesirable since they can induce oscillations in the reaction rate or heat generation rate in the catalysts through catalyst temperature feedback. For example, flow oscillations in the channels will cause changes in the heat transfer rates and consequently of the catalyst temperature in turn causing the power release to oscillate. This will affect the performance of the reactor significantly. Hence, it is essential not only to identify the stable and unstable operating regions of such reactors but also to find out various methods which can improve the stability of the operating as well as the future reactor systems. In this work, for the first time, the natural circulation stability behaviour of an EO (ethylene oxide) chemical reactor was investigated over a wide range of operating conditions and geometric parameters through numerical analysis. While simple two-phase flow models have been extensively used to predict the stability behaviour of natural circulation reactor systems; on the other hand their applicability for the chemical reactor systems has never been studied earlier. In the present work, the limitations of simple two-phase
(a) Conventional homogeneous two-phase flow models with empirical relations for void fraction and two-phase friction factor multiplier overestimate the natural circulation flow of the reactor. (b) The two-fluid model predicts the natural circulation flow closest to the measured value of the reactor. (c) The natural circulation flow rate in the reactor is stable at 18 MW power. The reactor can experience Type II instability if operating at higher power. (d) An increase in power or a decrease in subcooling has a destabilising effect on the natural circulation. (e) The margin to instability increases with decreasing riser length or combining the riser pipes. (f) With increasing riser area, the Type II instabilities disappears but a low-quality Type I instability appears.
Appendix A. Model for steady-state natural circulation The model assumes that in the single-phase region (a) fluid is incompressible, (b) Boussinesq approximation is valid for the variation of density with temperature. In the two-phase region op/ot in energy conservation is neglected; subcooled boiling is neglected; both phases are in thermodynamic equilibrium; no carry-over and carry-under takes place in the steam drum, and (e) heat losses in the loop piping are neglected.
(a) (b) (c) (d)
With these assumptions the conservation equations for one-dimensional two-phase flow for steady-state conditions are given by Mixture mass ow ¼ 0; oz
ðA1Þ
where w is the mass flow rate. Mixture momentum 1 o 2 f op ðw =qÞ þ qg cos / þ ¼0 ðw2 =qÞ þ 2 2 oz oz A 2DA
ðA2Þ
where q is the density of fluid, f is the friction factor, D is the hydraulic diameter and p is the pressure.
A.K. Nayak et al. / Applied Thermal Engineering 26 (2006) 644–653
Mixture energy QA heated region; oh w ¼ oz 0 adiabatic region
653
specific enthalpy conditions for the boundary condition that the closed loop pressure loss is zero. ðA3Þ References
where Q is the volumetric heat addition rate to the mixture. The local pressure losses due to bends and restrictions are estimated as DP k ¼ Kw2 =2qA2
ðA4Þ
where K is the loss coefficient. In the single-phase heated region, the variation of density is obtained using the Boussinesq approximation as given by qðzÞ ¼ q0 ð1 XðT ðzÞ T in ÞÞ ffi q0 ð1 XðhðzÞ hin Þ=C p Þ
ðA5Þ
where X is the volumetric thermal expansion coefficient and Cp is the specific heat of the mixture. In the two-phase heated region, the variation of density is given by qðzÞ ¼ aðzÞqg þ ð1 aðzÞÞqf
ðA6Þ
where a(z) is the void fraction at any axial location. In the present work, the homogeneous model and the Zuber–Findlay model were used in order to understand the applicability of them for the reactor thermosyphon system. In the homogeneous two-phase flow model, the void fraction is given by aðzÞ ¼
1 1 xðzÞ qg 1þ xðzÞ qf
ðA7Þ
Using the Zuber–Findlay drift flux model, the void fraction is given by aðzÞ ¼
jg ðzÞ C 0 ðjf ðzÞ þ jg ðzÞÞ þ V gj
ðA8Þ
where jg is the superficial velocity of steam, jf is the superficial velocity of liquid, C0 is a distribution parameter and Vgj is the drift velocity. In the two-phase region, a two-phase friction factor multiplier (U2lo ) was used to estimate the frictional pressure loss. In this work, the homogeneous model and the the Martinelli–Nelson model [18] were used to estimate the two-phase friction factor multiplier. Eqs. (A1)– (A4) are solved together to obtain the steady-state flow rate for a given power and reactor inlet temperature or
[1] B.S. Shiralkar, Md. Alamgir, J.G.M. Andersen, Thermal hydraulic aspects of the SBWR design, in: Proceedings of the International Conference on Design and Safety of Advanced Nuclear Power Plants (ANP92), AESJ, 1992. [2] A.A. Lammers, JGC Corporation Yokohama HPS Reactor Thermosyphon Cooling System, No. AT171.EB3601, September, 1989. [3] J.A. Boure, A.E. Bergles, L.S. Tong, Review of two-phase flow instability, Nucl. Eng. Des. 25 (1973) 165–192. [4] G. Yadigaroglu, Two-phase flow instabilities and propagation phenomena, Hemisphere Publishing Corporation, 1981, pp. 353– 402 (Chapter 17). [5] K. Fukuda, T. Kobori, Classification of two-phase flow instability by density wave oscillation model, J. Nucl. Sci. Technol. 16 (1979) 95–108. [6] V.K. Chexal, A.E. Bergles, Two phase instabilities in a low pressure natural circulation loop, AIChE Symp. Ser. 69 (1973) 37– 45. [7] S.Y. Lee, M. Ishii, Thermally induced flow oscillation in vertical two-phase natural circulation loop, Nucl. Eng. Des. 122 (1990) 119–132. [8] S.Y. Lee, D.W. Lee, Linear analysis of flow instabilities in an open two-phase natural circulation loop, Nucl. Eng. Des. 128 (1991) 317–330. [9] J.H. Chiang, M. Aritomi, M. Mori, Transient behavior of natural circulation for boiling two-phase flow (flow after pump trip), Nucl. Eng. Des. 152 (1994) 263–275. [10] Z. Zhou, Stability structure of low quality density wave oscillations in a natural circulation system at heating reactor conditions, Kerntechnik 62 (1997) 154–159. [11] S. Guanghui, J. Dounan, K. Fukuda, G. Yujun, Theoretical and experimental study on density wave oscillation of two-phase natural circulation of low equilibrium quality, Nucl. Eng. Des. 215 (2002) 187–198. [12] V. Chatoorgoon, SPORTS—a simple non-linear thermalhydraulic stability code, Nucl. Eng. Des. 93 (1986) 51–67. [13] J. Paniagua, U.S. Rohatgi, V. Prasad, Modeling of thermal hydraulic instabilities in single heated channel loop during startup transients, Nucl. Eng. Des. 193 (1999) 207–226. [14] N. Hoyer, MONA, A 7-equation transient two-phase flow model for LWR dynamics, in: Proceedings of International Conference on New Trends in Nuclear System Thermal Hydraulics, 1994. [15] N. Hoyer, Validation of the transient two-phase flow model MONA against density-wave oscillations, in: Proceedings of the First International Symposium on Two-phase Flow Modelling and Experimentation, Rome, Italy, 1994. [16] R. Zboray, An experimental and modelling study of naturalcirculation boiling water reactor dynamics, Ph.D. Thesis, Delft University of Technology, 2002. [17] A.K. Nayak, D. Lathouwers, T.H.J.J. van der Hagen, A.N.R. Bos, F.J.M. Schrauwen, A numerical study of a closed loop thermosyphon system, in: Proceedings of the Third International Conference on Computational Heat and Mass Transfer (APM2003), 2003, pp. 1–10. [18] R.C. Martinelli, D.B. Nelson, Prediction of pressure drop during forced circulation boiling of water, Trans. ASME 70 (1948) 695.
A numerical study of boiling flow instability of a reactor thermosyphon system A.K. Nayak a
a,*
, D. Lathouwers a, T.H.J.J. van der Hagen a, Frans Schrauwen b, Peter Molenaar b, Andrew Rogers b
Interfaculty Reactor Institute, Delft University of Technology, Mekelweg 15, 2629 JB Delft, The Netherlands b Shell Research and Technology Centre, Badhuisweg 3, 1031 CM Amsterdam, The Netherlands Received 27 July 2004; accepted 22 May 2005 Available online 8 November 2005
Abstract A numerical study has been carried out to investigate the boiling flow instability of a reactor thermosyphon system. The numerical model solves the conservation equations of mass, momentum and energy applicable to a two-fluid and three-field steam–water system using a finite difference technique. The computer code MONA was used for this purpose. The code was applied to the thermosyphon system of an EO (ethylene oxide) chemical reactor in which the heat released by a catalytic reaction is carried by boiling water under natural circulation conditions. The steady-state characteristics of the reactor thermosyphon system were predicted using the MONA code and conventional two-phase flow models in order to understand the model applicability for this type of thermosyphon system. The two-fluid model was found to predict the flow closest to the measured value of the plant. The stability behaviour of the thermosyphon system was investigated for a wide range of operating conditions. The effects of power, subcooling, riser length and riser diameter on the boiling flow instability were determined. The system was found to be unstable at higher power conditions which is typical for a Type II instability. However, with an increase in riser diameter, oscillations at low power were observed as well. These are classified as Type I instabilities. Stability maps were predicted for both Type I and Type II instabilities. Methods of improving the stability of the system are discussed. Ó 2005 Published by Elsevier Ltd. Keywords: Natural circulation; Chemical; Reactor; Two-phase flow; Stability; Type I instability; Type II instability
1. Introduction Use of a thermosyphon is the preferred option for heat removal of many industrial systems, such as nuclear reactors [1] and chemical reactors [2], since it eliminates problems associated with pump maintenance and operation. However, it is well known that this mode of cooling may suffer from instabilities which can hinder the operation and safety of plants. Many types of insta* Corresponding author. Address: Reactor Engineering Division, Bhabha Atomic Research Centre, Trombay, Mumbai 400 085, India. E-mail addresses: nayak_arun@rediffmail.com, arunths@apsara. barc.ernet.in (A.K. Nayak).
1359-4311/$ - see front matter Ó 2005 Published by Elsevier Ltd. doi:10.1016/j.applthermaleng.2005.05.019
bilities are known to occur in boiling two-phase systems; some of which depend on the static or steady-state characteristics of the system, whereas others depend on the dynamic characteristics of the system [3,4]. A densitywave instability is the typical dynamic instability which may occur due to the multiple regenerative feedback between the flow rate, enthalpy, density and pressure drop in the boiling system. The occurrence of the instability depends on the perturbed pressure drop in the two-phase and single-phase regions of the system and the propagation time delay of the void fraction or density in the system. Such an instability can occur at very low power and at high power conditions. This depends on the relative importance of the respective components
A.K. Nayak et al. / Applied Thermal Engineering 26 (2006) 644–653
645
Nomenclature A Cp D f g h hfg K L
cross-sectional area, m2 specific heat, J/kg K hydraulic diameter, m friction coefficient acceleration due to gravity, m/s2 specific enthalpy, J/kg latent heat of vapourisation, J/kg loss coefficient length, m
of pressure drop such as gravity or frictional losses in the system. Fukuda and Kobori [5] have classified the density-wave instability as Type I and Type II for the low power and high power instabilities respectively. The mechanisms can be explained as follows. 1.1. Type I instability For this type of instability to occur, the presence of a long riser plays an important role such as in a boiling two-phase natural circulation loop. Under low quality conditions, a slight change in quality due to any disturbance can cause a large change in void fraction (Fig. 1) and consequently in the driving head. Therefore the flow can oscillate at such low power conditions. But as the power increases, the flow quality increases where the slope of the void fraction vs. quality reduces. This can suppress the fluctuation of the driving head for a small change in quality. Hence, the flow stabilises at higher power. 1.2. Type II instability Unlike the Type I instability, the Type II instability occurs at high power conditions. This instability is driven by the interaction between the single and two-phase frictional component of pressure losses, mass flow, void
α
x Fig. 1. Typical variation of void fraction with quality.
p Q T t w x z a q
pressure, N/m2 volumetric heat addition rate, W/m3 temperature, K time, s mass flow rate, kg/s quality axial distance, m vapour volume fraction density, kg/m3
formation and propagation in the two-phase region. At high power, the flow quality or void fraction in the system is very large. Hence, the two-phase frictional pressure loss may be high owing to the smaller two-phase mixture density. Having a large void fraction will increase the void propagation time delay in the twophase region of the system. Under these conditions, any small fluctuation in flow can cause a larger fluctuation of the two-phase frictional pressure loss due to fluctuation of density and flow, which propagates slowly in the twophase region. On the other hand, the fluctuation of the pressure drop in the single-phase region occurs due to fluctuation of flow alone since the fluctuation of the density is negligible. The pressure drop fluctuation in this region travels much faster due to incompressibility of single-phase region. If the two-phase pressure drop fluctuation is equal in magnitude but opposite in phase with that of the single-phase region, the fluctuation or oscillation is sustained in the system since there are no attenuating mechanisms. Divergent oscillations can occur depending on the magnitude of the pressure loss fluctuation in the two-phase and single-phase regions and the propagation time delay. Both Type I and Type II instabilities can occur in a boiling thermosyphon system depending on its geometry and operating conditions. A literature review suggests that there have been many studies on this phenomenon in the past in two-phase natural circulation systems [6– 11]. Many numerical studies have also been performed in the past to model the density-wave instability [12,13]. However, there are no studies on boiling twophase thermosyphon phenomena inside the vessel of a chemical reactor such as an EO reactor. In these reactors, the heat released by chemical reaction in the catalyst tubes, is carried by boiling steam and water mixture under natural circulation conditions. Flow instability is highly undesirable, since oscillations of flow can cause oscillations of the catalyst temperature which can induce oscillations in heat generation rate in the catalysts. In that case, it is impossible to run the reactor for safety reason. Hence, it is essential to investigate the flow
646
A.K. Nayak et al. / Applied Thermal Engineering 26 (2006) 644–653
stability behaviour of such reactors in order to find out methods of improving the stability characteristics of operating reactor systems and scope for upscaling of future reactors. In EO reactors, the diameter of the vessel in which heat is generated due to a chemical reaction from a large number of catalyst tubes is very large (typically 5 m). Hence, the natural circulation velocity in the reactor is very small. Under these conditions, flow stratification can occur in the vessel. With spacers present to support the catalyst tubes, the two-phase flow can have a very complex distribution. To solve the stability problem of this system, use of a simple homogeneous model or a drift flux model as was mostly used by previous investigators, cannot give accurate results since these models do not consider flow pattern transitions and use empirical correlations for the friction factor and void fraction distribution which may not be applicable for all conditions. Further, treating multi-dimensional two-phase flows using numerical codes is complex. So for studying this phenomenon, the flow is assumed to be one-dimensional. It is required to predict the natural circulation stability characteristics using a robust two-fluid model that considers the transport equations of each phase separately. The objective of this study is to analyse the stability of boiling two-phase thermosyphon phenomena inside an EO (ethylene oxide) chemical reactor using a numerical two-fluid and three-field model. The model comprises conservation equations of mass, momentum and energy applicable to a steam–water flow and is numerically solved using the MONA code. First, the code was applied for simulation of the steady-state natural circulation behaviour of a chemical reactor for which a flow rate measurement at operating power was available. In order to understand the applicability of conventional two-phase flow models for void fraction and friction factor multiplier for this type of thermosyphon system, a theoretical model was developed. The model comprises simplified conservation equations of mass, momentum and energy for steady-state conditions. Different void fraction models are considered to calculate the density in the two-phase region and different two-phase friction factor multiplier models are used to calculate the frictional pressure loss in the two-phase region. The MONA code predictions were compared with that of the simplified model and the measured plant data. The limitations of the conventional models for application to such systems were determined. The stability behaviour of the system was analysed considering the effects of power and subcooling as parameters. The effects of riser diameter and riser length on flow stability were simulated. Stability maps were predicted for different operating conditions and geometries of the reactor. Methods of improving the stability of the system are discussed.
2. Computer code MONA The computer code MONA [14] simulates the flow of steam–water/inert vapour in any pipe line assuming one-dimensional flow. The code is basically an updated version of the RAMONA code used for nuclear reactor thermal hydraulic analysis. MONA contains a set of seven conservation equations, based on the modelling of three-flow fields; the bulk liquid or the liquid film at the wall, the liquid droplets and the vapour or steam. The code solves three mass conservation equations, two momentum conservation equations and two energy conservation equations. Closure relations provide expressions for two-phase friction and comprise a complete set of correlations for heat transfer at the wall and between the respective phases. A phenomenological approach is employed in MONA to model the different flow regimes and the transition between them, thus no separate user specified correlations for void distribution, etc. are required. The exchange of energy and momentum at the interface between vapour and liquid depends on the interfacial area and the topology of the two-phase flow. So the flow regime determines the wetted perimeters and friction factors. The interfacial heat transfer rates are estimated from the interfacial heat transfer coefficient, interfacial area and the temperature difference between the phases. The heat transfer coefficient is determined from using suitable correlations depending on the flow pattern. The MONA code uses a staggered mesh scheme in which the flow variables (velocities, mass flows, etc.) are defined on section boundaries and scalar variables (pressure, void fraction, densities, etc.) are located in the centre of section volumes. Donor cell differencing is adopted. For ease in solution, the mass conservation equations are decoupled from the momentum, energy and volume equations. The time step is limited by the Courant criterion based on the mass propagation velocity. The code has been extensively validated against data for various two-phase thermal hydraulic conditions including pressure wave propagation, U-tube flow oscillations, void fraction, pressure drop, etc. [15,16].
3. Theoretical model for steady-state natural circulation In order to understand the applicability of conventional two-phase flow models for void fraction and friction factor multiplier for the prediction of steady-state natural circulation in a chemical reactor, a theoretical model was developed. A brief description of the model is given in Appendix A. The model predicts the flow rate in a natural circulation loop for a given operating power, subcooling and pressure. The model incorporates various two-phase friction factor multiplier models and void fraction models which can influence on the natural circulation characteristics of the loop.
A.K. Nayak et al. / Applied Thermal Engineering 26 (2006) 644–653
used as described below in order to understand the model applicability for simulation of this type of thermosyphon system:
4. Geometry of the reactor thermosyphon system In the present analysis, we have considered the thermosyphon behaviour inside the vessel of an EO chemical reactor [17]. This is because the measured natural circulation flow rate of the plant at the operating power is available and it is useful for validation of the models. Fig. 2 shows a schematic outline of the reactor considered. The reactor vessel diameter is very large (5 m) which houses many closed tubes containing the EO catalyst. The height of the vessel is about 10 m. The diameter of riser pipes is of the order of 30 mm and the total height of the system is about 15 m. There are 10 spacers in the vessel to hold the catalyst tubes in position and the spacer loss coefficient was considered to be about 3.5. The local losses in the riser section was considered to be about 5.75. The heat released by a chemical reaction is carried by a steam–water mixture through natural circulation. The steam pressure is maintained at about 4.74 MPa in order to maintain the required temperature of the catalyst. Subcooled water from the downcomer enters the vessel and is heated as the fluid flows upward due to buoyancy. Boiling takes place and the steam and water mixture moves to a steam drum through small size riser pipes, wherein the steam and water separation takes place by gravitational settling. The steam is sent to a cooler and an equal amount of feedwater at lower temperature joins the saturated water in the steam drum.
(i) homogeneous model for void fraction and twophase friction factor multiplier; (ii) homogeneous model for void fraction and Martinelli–Nelson model for two-phase friction factor multiplier; (iii) Zuber–Findlay model for void fraction and Martinelli–Nelson model for two-phase friction factor multiplier, and (iv) two-fluid/three-field model as implemented in the MONA code. The results show that (Fig. 3) even though the trend for the variation of the natural circulation flow rate with power is similar, the magnitude of the flow rate as predicted by the various models is quite different. With increasing power, all the models predict that the flow rate initially rises due to an increase in buoyancy force and then reduces due to an increase in two-phase resistance at higher void fraction. The homogeneous model gives the largest flow rate since this model predicts a higher void fraction and smaller two-phase frictional pressure loss as compared to the other models at given power. Of course, when considering slip, the flow rate reduces, but the flow rate was still found to be higher than the measured natural circulation flow rate of the reactor at 18 MW power and subcooling of 2 K. The two-fluid model predicts a flow rate of about 183 kg/s at the operating power, which is closest to the measured value of the reactor at 18 MW power. The reason behind this is that the two-fluid model predicts a higher pressure drop in the two-phase region as compared to the other models.
5. Steady-state characteristics of the system The steady-state behaviour of the reactor thermosyphon system as described above, was predicted at different power keeping the steam drum pressure and downcomer subcooling constant. Several models were
Steam Steam drum Feed water Riser pipes
Reactor vessel
647
Steam and water mixture
Catalyst assemblies
downcomers
Steam + Water
Water
Fig. 2. Schematic of the thermosyphon system considered in the present analysis.
648
A.K. Nayak et al. / Applied Thermal Engineering 26 (2006) 644–653 350 homogenous void fraction model + Martinelli-Nelson Two-phase Multiplier model
homogeneous model
300
flow (kg/s)
250
Zuber-Findlay void-fraction model + Martinelli-Nelson Two-phase multiplier model
200
two-fluid model
150
subcooling = 2K
measurement
pressure = 4.74 MPa
100 0
10
20
30
40
50
60
70
80
power (MW) Fig. 3. Comparison of steady-state flow rate using different models.
Table 1 Comparison of pressure drops in different sections of the thermosyphon as predicted by different models Component
HEM
HEM for void fraction + M–N for two-phase multiplier
Zuber–Findlay for void fraction + M–N for two-phase multiplier
Two-fluid model
Heated vessel Horizontal riser Vertical riser Steam drum Downcomer (vertical)
0.532 bar 0.0127 bar 0.279 bar 0.077 bar 1.129 bar
0.533 bar 0.0343 bar 0.28 bar 0.077 bar 1.129 bar
0.72 bar 0.0345 bar 0.288 bar 0.077 bar 1.129 bar
0.78 bar 0.08 bar 0.346 bar 0.077 bar 1.129 bar
Table 1 shows a comparison of the pressure loss in various components at a flow rate of 183 kg/s for a power of 18 MW and a subcooling of 2 K. The difference in pressure drop among different models is found to be in the heated vessel and in the horizontal piping mostly due to the occurrence of flow stratification in the large diameter vessel and the horizontal riser pipe.
6. Stability characteristics of the system The stability characteristics of the reactor were predicted using the MONA code for different operating conditions and geometric parameters in order to understand their influence on stability. 6.1. Effects of operating conditions on stability Fig. 4 shows the influence of operating power on the stability of the system using the two-fluid model. The thermosyphon flow behaviour at different power levels was predicted at an operating pressure of 4.74 MPa and a reactor inlet subcooling of 3 K. For all the conditions, the computations were started with an assumed
initial flow rate of 283 kg/s, zero void fraction in the heated and adiabatic regions and a uniform pressure (4.74 MPa) throughout the loop. A sensitivity analysis was done to investigate the effects of node sizes and time step on thermosyphon stability phenomena. The code predictions were found to be insensitive to this variation. At low power, the flow is found to be stable. With a rise in power from 10 MW to 18 MW, the thermosyphon flow is found to increase. When the power is further increased to 75 MW, small amplitude flow oscillations appear which die out as time proceeds. However, the corresponding flow rate is found to be lower than that at 18 MW. When the power is further increased to about 95 MW, a limit cycle characterised by constant amplitude and frequency of flow oscillations, is observed. Similar behaviour is also observed at a subcooling of 5 K as shown in Fig. 5. In this case, the computations were started with the same assumed conditions as in the previous case (Fig. 4). However, the starting transients up to 600 s are not shown here. The flow is found to be stable up to a power of 80 MW (Fig. 5(a)). With a further increase in power, the flow becomes unstable as expected (Fig. 5(b)). The frequency and amplitude of oscillations were found to increase with increasing
A.K. Nayak et al. / Applied Thermal Engineering 26 (2006) 644–653
649
450 400 350
flow rate (kg/s)
300
subcooling = 3 K pressure = 4.74 MPa
250
18 MW
10 MW
200 150 100
95 MW
75 MW
50 0 0
200
400
600
800
1000
1200
time (s) Fig. 4. Effect of power rise on natural circulation behaviour at 3 K subcooling.
400 350 300
flow (kg/s)
250 75 MW 200
80 MW
150 100 subcooling = 5 K pressure = 4.74 MPa
2.5 MW 50 600
800
(a)
1000
1200
time (s) 200 85 MW
subcooling = 5 K
flow (kg/s)
150
100
95 MW
90 MW
105 MW
50 600
(b)
800
1000
1200
time (s)
Fig. 5. Effect of power rise on natural circulation behaviour at 5 K subcooling.
650
A.K. Nayak et al. / Applied Thermal Engineering 26 (2006) 644–653 250 power = 95 MW pressure = 4.74 MPa subcooling = 10 K
200
flow (kg/s)
subcooling = 5 K
150
100 subcooling = 15 K subcooling = 25 K
subcooling = 3 K
50 600
700
800
900
1000
time (s)
Fig. 6. Effect of subcooling on natural circulation behaviour.
power. With a rise in power, the flow rate also increases and hence, the time period of fluid to pass through the loop reduces. As a result, the frequency of oscillations which depends on the circulation period, increases. The increase in amplitude is due to the larger flow-void feedback which increases with a rise in power or void fraction in the system. The effect of reactor inlet subcooling on the flow stability is shown in Fig. 6. The operating power was fixed at 95 MW and the pressure at 4.74 MPa. The subcooling was increased from 3 K to 30 K. At low subcooling the flow was found to be unstable and stabilised above a subcooling of 15 K. Fig. 7 shows the stability map plotted in the powersubcooling plane. Since the flow reduces with an increase in power, the indication of increase of threshold power for instability with subcooling, clearly indicates that the type of instability observed in this system is 140 120 Unstable
power (MW)
100 80 60
Stable
40 20 0 0
10
20
30
subcooling (K) Fig. 7. Stability map of the system.
40
50
the so-called Type II instability. This type of instability occurs due to the dominance of frictional losses occurring at higher power when the natural circulation flow rate and two-phase mixture density is small. So, in order to have larger stability margin for the operating system, it is recommended to operate at higher subcooling. 6.2. Effect of riser geometry on stability The riser length and area are two important parameters that play a significant role in the natural circulation characteristics of a system. An increase in riser length or a decrease in riser area increases the two-phase resistance which can reduce the natural circulation flow rate and increases the void fraction for the same heating power. This may affect the performance of the reactor because of low flow and high void fraction conditions which favour Type II density-wave instabilities as discussed earlier. On the other hand, a decrease in riser length or an increase of flow area, even though increases the natural circulation flow rate, reduces the flow quality or void fraction for the same operating power. Under these conditions, slug flow may occur in the riser pipes which can cause undesirable flow induced vibrations. Besides, operation under such low quality conditions may induce the low quality Type I instability. Hence, it is essential to investigate the effects of the riser geometry on the natural circulation characteristics of the reactor. This study is important not only for finding the scope for improving the stability of current reactors but also in finding the effects of increasing the reactor size (up-scaling) on the stability. In this work, the horizontal riser length and the riser area were varied keeping all other parameters constant and their effects on stability were predicted using the MONA code.
A.K. Nayak et al. / Applied Thermal Engineering 26 (2006) 644–653
140
unstable
651
250
length = 5.7 m
1.2 Power = 10 MW subcooling = 2 K
length = 10.7 m
200
length = 15.7 m
80
Area = 0.1459 m2 Area = 0.2918 m2 Area = 0.4377 m2
stable 60
0.8
0.6 150
40
0.4
pressure = 4.74 MPa
20 0 0
10
20
void fraction
1.0
100
flow (kg/s)
threshold power (MW)
120
30
40
100
50
200
400
600
subcooling (K)
800
1000
1200
0.2 1600
1400
time (s)
Fig. 8. Effect of the horizontal riser length on stability.
Fig. 10. Effect of an increase in riser area on stability.
Fig. 8 shows the effect of riser length on flow stability. The length of the horizontal riser was varied from 5.7 m to 15.7 m keeping its diameter constant. All other dimensions of the systems were unchanged as well. A decrease in riser length was found to be stabilising. With decreasing riser length, the two-phase resistance is reduced which stabilises the two-phase flow. Fig. 9 shows the effect of riser flow area on stability. In one case, the actual two riser pipes of flow area 0.0729 m2 each was considered. In the second case, the two riser pipes were combined into one. So even though the total flow area remained the same, the hydraulic diameter increased. This reduces the frictional resistances and hence stabilises the flow. If the riser area is increased further, instabilities can appear at low power condition as seen in Fig. 10. For smaller riser area (0.1459 m2), the flow rate is smaller due to larger resistance in small riser pipes. As the flow area is increased, the flow rate increases, which gives rise to small frequency oscillations, typical of low quality Type I density-wave instability [5]. Such an instability occurs in natural circulation systems having tall risers. Under low power conditions, a slight change in quality
can cause a large change in void fraction or flow rate to induce this instability in the system. The corresponding void fraction oscillations are also shown in Fig. 10. With an increase in riser area, the time period of oscillation reduces due to the increase in flow rate in the system. Fig. 11 shows the instability characteristics of the system for a riser area of 0.2918 m2 (double the nominal 0.5 300
0.4
flow (kg/s)
0.2 200
0.1 0.0
150
flow (5 MW ) flow (10 MW) flow (15 MW)
void fraction (5 MW) void fraction (10 MW) void fraction (15 MW) 400
600
-0.2
(a)
800
1000
1200
1400
time (s) 350
1.0
unstable
140
one pipe - flow area 0.1459 m2
120
void fraction (120 MW) void fraction (140 MW)
80 two pipes - flow area of each 0.0729 m2 60 stable
40
0.8
300
flow rate (120 MW) flow rate (140 MW)
275
void fraction
325
100 flow (kg/s)
threshold power (MW)
-0.1
100 200
riser void fraction
0.3 250
0.6
20 250 200
0 0
10
20
30
40
50
(b)
400
600
800
1000
time (s)
subcooling (K) Fig. 9. Effect of the riser area on stability.
Fig. 11. Effect of an increase in power on stability for a riser area of 0.2918 m2.
652
A.K. Nayak et al. / Applied Thermal Engineering 26 (2006) 644–653
flow models for analyzing the natural circulation behaviour of such reactor systems were verified. The main conclusions of the study are summarised below.
70
threshold power (MW)
60 Stable 50 subcooling = 2 K
40 30
Unstable 20 10 1.5
2.0
2.5
3.0
3.5
4.0
4.5
5.0
5.5
riser area/nominal area Fig. 12. Stability boundary for Type I instability.
area) for different powers. Unlike for the nominal case as seen before (Figs. 4 and 5), with increasing riser area, only low power instability occurs in the system. The high power instability disappears as the riser area is increased. Similar behaviour is observed with further increase in riser area also. Type I instability maps as observed in this case are plotted in Fig. 12 for the threshold power of instability as a function of riser area at a subcooling of 2 K. It is interesting to note that the threshold power for the flow to become stable increases as the riser area is increased.
7. Conclusions This study has focussed on the prediction of the flow stability behaviour in natural circulation cooled chemical reactors, which are highly undesirable since they can induce oscillations in the reaction rate or heat generation rate in the catalysts through catalyst temperature feedback. For example, flow oscillations in the channels will cause changes in the heat transfer rates and consequently of the catalyst temperature in turn causing the power release to oscillate. This will affect the performance of the reactor significantly. Hence, it is essential not only to identify the stable and unstable operating regions of such reactors but also to find out various methods which can improve the stability of the operating as well as the future reactor systems. In this work, for the first time, the natural circulation stability behaviour of an EO (ethylene oxide) chemical reactor was investigated over a wide range of operating conditions and geometric parameters through numerical analysis. While simple two-phase flow models have been extensively used to predict the stability behaviour of natural circulation reactor systems; on the other hand their applicability for the chemical reactor systems has never been studied earlier. In the present work, the limitations of simple two-phase
(a) Conventional homogeneous two-phase flow models with empirical relations for void fraction and two-phase friction factor multiplier overestimate the natural circulation flow of the reactor. (b) The two-fluid model predicts the natural circulation flow closest to the measured value of the reactor. (c) The natural circulation flow rate in the reactor is stable at 18 MW power. The reactor can experience Type II instability if operating at higher power. (d) An increase in power or a decrease in subcooling has a destabilising effect on the natural circulation. (e) The margin to instability increases with decreasing riser length or combining the riser pipes. (f) With increasing riser area, the Type II instabilities disappears but a low-quality Type I instability appears.
Appendix A. Model for steady-state natural circulation The model assumes that in the single-phase region (a) fluid is incompressible, (b) Boussinesq approximation is valid for the variation of density with temperature. In the two-phase region op/ot in energy conservation is neglected; subcooled boiling is neglected; both phases are in thermodynamic equilibrium; no carry-over and carry-under takes place in the steam drum, and (e) heat losses in the loop piping are neglected.
(a) (b) (c) (d)
With these assumptions the conservation equations for one-dimensional two-phase flow for steady-state conditions are given by Mixture mass ow ¼ 0; oz
ðA1Þ
where w is the mass flow rate. Mixture momentum 1 o 2 f op ðw =qÞ þ qg cos / þ ¼0 ðw2 =qÞ þ 2 2 oz oz A 2DA
ðA2Þ
where q is the density of fluid, f is the friction factor, D is the hydraulic diameter and p is the pressure.
A.K. Nayak et al. / Applied Thermal Engineering 26 (2006) 644–653
Mixture energy QA heated region; oh w ¼ oz 0 adiabatic region
653
specific enthalpy conditions for the boundary condition that the closed loop pressure loss is zero. ðA3Þ References
where Q is the volumetric heat addition rate to the mixture. The local pressure losses due to bends and restrictions are estimated as DP k ¼ Kw2 =2qA2
ðA4Þ
where K is the loss coefficient. In the single-phase heated region, the variation of density is obtained using the Boussinesq approximation as given by qðzÞ ¼ q0 ð1 XðT ðzÞ T in ÞÞ ffi q0 ð1 XðhðzÞ hin Þ=C p Þ
ðA5Þ
where X is the volumetric thermal expansion coefficient and Cp is the specific heat of the mixture. In the two-phase heated region, the variation of density is given by qðzÞ ¼ aðzÞqg þ ð1 aðzÞÞqf
ðA6Þ
where a(z) is the void fraction at any axial location. In the present work, the homogeneous model and the Zuber–Findlay model were used in order to understand the applicability of them for the reactor thermosyphon system. In the homogeneous two-phase flow model, the void fraction is given by aðzÞ ¼
1 1 xðzÞ qg 1þ xðzÞ qf
ðA7Þ
Using the Zuber–Findlay drift flux model, the void fraction is given by aðzÞ ¼
jg ðzÞ C 0 ðjf ðzÞ þ jg ðzÞÞ þ V gj
ðA8Þ
where jg is the superficial velocity of steam, jf is the superficial velocity of liquid, C0 is a distribution parameter and Vgj is the drift velocity. In the two-phase region, a two-phase friction factor multiplier (U2lo ) was used to estimate the frictional pressure loss. In this work, the homogeneous model and the the Martinelli–Nelson model [18] were used to estimate the two-phase friction factor multiplier. Eqs. (A1)– (A4) are solved together to obtain the steady-state flow rate for a given power and reactor inlet temperature or
[1] B.S. Shiralkar, Md. Alamgir, J.G.M. Andersen, Thermal hydraulic aspects of the SBWR design, in: Proceedings of the International Conference on Design and Safety of Advanced Nuclear Power Plants (ANP92), AESJ, 1992. [2] A.A. Lammers, JGC Corporation Yokohama HPS Reactor Thermosyphon Cooling System, No. AT171.EB3601, September, 1989. [3] J.A. Boure, A.E. Bergles, L.S. Tong, Review of two-phase flow instability, Nucl. Eng. Des. 25 (1973) 165–192. [4] G. Yadigaroglu, Two-phase flow instabilities and propagation phenomena, Hemisphere Publishing Corporation, 1981, pp. 353– 402 (Chapter 17). [5] K. Fukuda, T. Kobori, Classification of two-phase flow instability by density wave oscillation model, J. Nucl. Sci. Technol. 16 (1979) 95–108. [6] V.K. Chexal, A.E. Bergles, Two phase instabilities in a low pressure natural circulation loop, AIChE Symp. Ser. 69 (1973) 37– 45. [7] S.Y. Lee, M. Ishii, Thermally induced flow oscillation in vertical two-phase natural circulation loop, Nucl. Eng. Des. 122 (1990) 119–132. [8] S.Y. Lee, D.W. Lee, Linear analysis of flow instabilities in an open two-phase natural circulation loop, Nucl. Eng. Des. 128 (1991) 317–330. [9] J.H. Chiang, M. Aritomi, M. Mori, Transient behavior of natural circulation for boiling two-phase flow (flow after pump trip), Nucl. Eng. Des. 152 (1994) 263–275. [10] Z. Zhou, Stability structure of low quality density wave oscillations in a natural circulation system at heating reactor conditions, Kerntechnik 62 (1997) 154–159. [11] S. Guanghui, J. Dounan, K. Fukuda, G. Yujun, Theoretical and experimental study on density wave oscillation of two-phase natural circulation of low equilibrium quality, Nucl. Eng. Des. 215 (2002) 187–198. [12] V. Chatoorgoon, SPORTS—a simple non-linear thermalhydraulic stability code, Nucl. Eng. Des. 93 (1986) 51–67. [13] J. Paniagua, U.S. Rohatgi, V. Prasad, Modeling of thermal hydraulic instabilities in single heated channel loop during startup transients, Nucl. Eng. Des. 193 (1999) 207–226. [14] N. Hoyer, MONA, A 7-equation transient two-phase flow model for LWR dynamics, in: Proceedings of International Conference on New Trends in Nuclear System Thermal Hydraulics, 1994. [15] N. Hoyer, Validation of the transient two-phase flow model MONA against density-wave oscillations, in: Proceedings of the First International Symposium on Two-phase Flow Modelling and Experimentation, Rome, Italy, 1994. [16] R. Zboray, An experimental and modelling study of naturalcirculation boiling water reactor dynamics, Ph.D. Thesis, Delft University of Technology, 2002. [17] A.K. Nayak, D. Lathouwers, T.H.J.J. van der Hagen, A.N.R. Bos, F.J.M. Schrauwen, A numerical study of a closed loop thermosyphon system, in: Proceedings of the Third International Conference on Computational Heat and Mass Transfer (APM2003), 2003, pp. 1–10. [18] R.C. Martinelli, D.B. Nelson, Prediction of pressure drop during forced circulation boiling of water, Trans. ASME 70 (1948) 695.