2020
Ind. Eng. Chem. Res. 2009, 48, 2020–2033
Validations of Thermohydraulic Models for Geysering in a Natural Circulation Loop Using an Impedance Needle Probe S. Paruya,† A. K. Saha,‡ and P. Bhattacharya*,§ Department of Chemical Engineering, National Institute of Technology, Durgapur, India-713209, Department of Chemical Engineering, Haldia Institute of Technology, India-721657, and Department of Chemical Engineering, JadaVpur UniVersity, India-700032
This paper describes the validations of thermohydraulic models that predict the geysering instability observed at the riser section of a short-tube natural circulation evaporator (STNCE) at low heat flux and low pressure. The validations have been made by measuring the void fraction in the steam-water system at the riser section of a STNCE with the help of an impedance void-needle probe. The calibration of the probe shows that in the range of low void fractions the probe indicates lower values while in the range of high void fractions it indicates higher values compared to those measured by the manometer. Measured data of void fraction have been analyzed by thermohydraulic models. Both measured data and simulated data of temporal variation of void fractions during geysering have been observed to be periodic pulses of high amplitudes. Effects of riser height and loop resistance on the geysering frequency have been demonstrated in detail. The paper also discusses both the measurement uncertainties and the model uncertainties. 1. Introduction Studies on thermal hydraulics in a boiling natural circulation loop have received significant attention in regard to boiling water reactors (BWRs), simplified boiling water reactors (SBWRs), nuclear steam generators, thermal drum-type boilers, and all the thermosyphon reboiler units of chemical, refining, and petrochemical plants. Natural circulation systems are preferred to forced circulation systems in these applications, as the former eliminates the use of circulation pumps, which can suffer occasional failures and require regular maintenance. Thus, the natural circulation systems offer better operational safety, although they are inherently more prone to flow instabilities because of low driving heads. During startup operations, the steam-water natural circulation (SWNC) loops experience various flow-instability phenomena. Among these, geysering is one of the most important phenomena occurring at low pressures and low heat fluxes. Geysering in the loop is described by regular or irregular cycles of phases due to appreciable change in static head along with its vertical boiling channels. The phases in a cycle consist of (1) refilling through condensation of vapor, increasing static head to stop boiling; (2) incubation to initiate steam formation; and (3) expulsion of vapor, as a result of self-evaporation induced by decrease in static head. Griffith first identified geysering phenomenon in a vertical boiling channel at low heat flux and low pressure.1 He described it as rapid expulsion of boiling liquid and its vapor from a vertical tube and reported a time period of geysering on the order of 0.1-0.01 Hz for the system studied. Murphy carried out experimental investigation to derive correlations between geysers and nongeysers for water, hydrogen, and R113 in a vertical pipe of large diameter.2 Boure et al.3 reported that geysering was an oscillatory phenomenon. Yadigaroglu4 carefully defined three sequences that occur in a geysering cycle, * To whom correspondence should be addressed. E-mail: pinaki_che@ yahoo.com. † National Institute of Technology. ‡ Haldia Institute of Technology. § Jadavpur University.
viz., (1) incubation, (2) vapor generation and expulsion, and (3) refilling. Ozwaw et al.5 clearly described the geysering instability as nonlinear oscillations. Aritomi et al.6 observed geysering phenomenon of periodic nature in an experiment on the startup of a steam-water natural circulation (SWNC) loop. They also observed that it was a condensation-induced thermal nonequilibrium phenomenon under subcooled boiling conditions and it occurred as a result of flow-regime transition from bubbly flow to slug flow. Wang et al.7 made experimental observations of the geysering instability in an SWNC loop with a single boiling channel to conclude that there were four quasiperiodic oscillations of liquid temperature at adiabatic riser (ARI) exit, and the time period of the oscillations of the largest one was about 50.0 s (0.02 Hz). Kyung and Lee8 reported the geysering of periodic nature, whereas Jiang et al.9 experimentally observed that the geysering phenomenon generates high-frequency oscillations with high peaks like pulses with and without regular frequency due to the combined effects of self-evaporation and flashing instability. Janssens-Maenhout et al.10 experimentally investigated the effect of riser width and initial static height on geysering frequency. Subki et al.11 observed geysering phenomenon in a double-channel boiling SWNC loop and confirmed that the startup at low pressures and low heat fluxes results in geysering instability of high frequencies and flashinginduced instability. Zhang et al.12 investigated the geysering of cryogenic propellant and recommended a method for suppressing its severity. The studies on the geysering phenomenon in the chemical, refining, and petrochemical disciplines have also been reported in the literature. Chexal and Bergles13 investigated two-phase flow instabilities in vertical thermosiphon used in the chemical process industry and found nonperiodic geysering oscillations. Geysering can also be found in the storage of cryogenic fluids,14 natural circulation evaporators, and distillation columns. Geysering coupled with manometer oscillations had been identified experimentally in evaporators.15 On the basis of experiments with cryogenic fluids, Kuncoro et al.16 concluded that liquid superheating was the driving force for vapor expulsion and pressure surge during geysering.
10.1021/ie8005462 CCC: $40.75 2009 American Chemical Society Published on Web 12/29/2008
Ind. Eng. Chem. Res., Vol. 48, No. 4, 2009 2021
A survey of the literature reported above indicates that the oscillations of liquid temperatures, loop flow rate, void fraction and static pressure, etc. take place during geysering phenomenon. The oscillations of the loop flow rate are strongly coupled with void fractions in the boiling channel as the void fractions play a key role in estimating the driving force for loop flow rate and two-phase frictional pressure drop. So, the measurement of void fraction in the boiling channel has immense importance for the detailed understanding of void-to-circulation sensitivity of the natural circulation loop of a given configuration and the instability phenomena as well. Various intrusive devices (that disturb the flow structure of a two-phase mixture) for void fraction measurements are based on the measurement of electrical conductance, impedance, and optical properties of a two-phase mixture as mentioned by Hewitt.17 A comprehensive review of impedance-based void fraction measurements has been made by Ceccio et al.18 Extensive uses of impedance probes for the fast transients have been recommended by many investigators19,20 because of the high speed of response and high sensitivity of the probes. Furthermore, they can be suitably operated in either conductance mode or capacitance mode by choosing the proper frequency of the ac signal for the excitation of the probe. Impedance has been measured by Bernier19 using two semicircular electrodes for quantifying void fractions for air-water bubbly flow. One drawback of the impedance probe, as mentioned by Andreussi et al.,20 is that its sensitivity is strongly related to the change of flow regime, and thus, finding a unique calibration curve for the probe is difficult. Several geometries of the probe, such as ring, arc, parallel wire, and parallel plates, are found in the literature. Kok et al.21 measured void fractions in subchannels of natural circulation loop facility DESIRE at both steady and unsteady conditions using γ-transmission tomography and established that the loop flow rate-tovoid fraction transfer function followed a first-order dynamics with a time constant of 0.1 s for the loop. Devis and Fossa22 employed an impedance method to measure void fractions in an air-water mixture flow in a horizontal tube and observed that ring electrodes were more sensitive to small phase discontinuities than half-ring electrodes. Manera23 measured void fractions during flashing-induced instability in natural circulation loop using a wire-mesh sensor and γ-ray technique. The study showed that at high sampling frequencies of 500 Hz and high signal-to-noise ratio, the wire-mesh sensor was capable of capturing the characteristic frequencies of 2-4 Hz, at which water and steam alternate in the test section during flashing. Chanson and Waniewski et al. determined the quantity of air in air-entrainment experiments using impedance needle probes and detected individual bubbles and their velocities and diameters.24,25 It is evident that the geysering phenomenon in a SWNC loop has not been extensively studied using impedance-base needle probes for measuring void fractions in high-frequency instability phenomena in steam-water system in a natural circulation loop. It is also evident that simplified mathematical formulations (HEM or four-equation scheme) for boiling two-phase flow dynamics have been carried out to explain the geysering phenomenon. It is, therefore, felt that the present knowledge on the phenomenon can be improved by making nonlinear characterization of geysering in the loop, such as limit cycles, and by formulating a unified mechanism for the same. For our purposes, both experimental and theoretical investigations on the parametric sensitivity of the phenomenon are required to be carried out. It is also realized that the performance evaluation of the impedance probes of simpler design (needle probe) in
Figure 1. Schematic of the SWNC loop for the investigations (after modifying the natural circulation evaporator).
such a high-frequency phenomenon should be studied in detail. The advantages of the probe are that it is simple in design and has a very high speed of response. It is, therefore, anticipated that the probe will be capable of capturing high-frequency instability phenomena compared to other probes. As it does not seriously affect the local flow structure, it helps in achieving relatively accurate local measurements. In the present investigation, the performance of the impedance needle probe has been analyzed in regard to identifying such high-oscillating instabilities using a complex moving-boundary-based thermohydraulic model developed for a natural circulation evaporator with a short riser.26 Validation of the proposed model has been carried out through the measurements of the void fraction in the riser section and fluid temperatures in various sections of the loop. Parametric sensitivities, such as the effect of riser length and loop resistance on geysering, have also been studied experimentally and theoretically to characterize the nonlinear oscillations. 2. Experimental Setup 2.1. Description of the Setup. A STNCE, supplied by K. C. Engineers, UP, India, was suitably modified to meet the objectives of the present investigation. The schematic of the modified experimental setup is shown in Figure 1. The setup consisted of a steam drum (steam-water separator), down comer, boiling channel (heated section and a short ARI), vapor condenser, condensate tank, and boiler. The dimensions of these components are mentioned in Table 1. In addition, there were rotameter, RTDs, pressure gauges, and other small accessories. The heated section was a vertical shell-and-tube heat exchanger (four tubes) in which the circulating water boiled inside tubes and the heating-steam condensed in the shell side to provide the heat required for boiling the water in the tubes. The steam-water mixture went upward through the tubes to reach the ARI, where liquid-flashing was intense due to the reduction of the static head. The vapor-liquid mixture then entered into the separator (steam drum) from where the vapor left through the top after being separated from its mixture. The water from the separator then passed through the down comer and mixed with fresh feedwater that came from the overhead condensate tank. The vapor from the steam drum went to the condenser, and subsequently, the condensate accumulated in a tank called the condensate tank. The subcooled vapor condensate was recirculated via rotameter to the tube side of the heated section to achieve inlet subcooling to some extent. RTDs were inserted in the equipment at different points to measure temperatures. Another rotameter was installed at the cooling water line to measure the cooling water inlet flow rate. One
2022 Ind. Eng. Chem. Res., Vol. 48, No. 4, 2009 Table 1. Geometry of the Experimental Setup shell tubes
1. Heated Section material shell i.d. length material tube o.d. tube i.d. no. of tubes 2. Riser Section
material i.d. length
SS 304 108.0 mm 500.0 mm SS 304 12.7 mm 9.50 mm 4 SS 304 108.0 mm 250.0 mm
3. Down Comer i.d. length (wrt riser exit) shell tube
25.0 mm 750.0 mm 4. Condenser material shell i.d. length material tube o.d. tube i.d. no. of tubes 5. Boiler
material capacity no. of heaters (3 kW)
SS 304 108.0 mm 500.0 mm SS 304 12.7 mm 9.50 mm 12 SS 304 12.0 L 2
Figure 2. Needle probe for void fraction measurement.
6. Electricity Supply 1-Phase, 220 V ac, 4 kW
manometer was installed for measuring the pressure drop between the top and bottom point of the separator. Another manometer was installed for measuring the pressure drop boiling channel, i.e., between the top point (from where vapor-liquid mixture goes to the separator section) and the bottom point (to which the subcooled vapor condensate is fed). A PT-100 temperature scanner was used for online measurement of temperatures at the heating-steam inlet point, heating-steam condensate point, feed inlet point, the point from where the twophase mixture leaves the riser section, and the point where the vapor leaves the separator. A personal computer (PC) was installed for online data-recording. Additionally, a digital storage oscilloscope (DSO) was installed for recording and storing graphical voltage output signals from the probe that detected the presence of bubbles and measured the void fraction in the ARI.27 2.2. Instrumentations. An impedance-based needle probe was fabricated to measure the void fraction (R) in the ARI in an unsteady two-phase steam-water flow with a maximum liquid velocity of 0.15 m/s. A function generator was used for ac excitation (500 kHz, (5.0 V) of the probe, and the dc signal after the demodulation of the output signal from the probe was analyzed with the help of a DSO (Aplab 36100 DCA, 100 MHz and 1.0 GSa/s). Platinum resistance thermometers (PT-100) with a sensitivity of 0.384 Ω/° C, a range of -75 to +250 °C, and an accuracy of (0.15 at 0 °C have been used. Locations of the sensor installation are shown in Figure 1. T1, T2, T3, T4, T5, T6, and T7 are the temperature sensors. The temperature sensors were interfaced with an online scanner, which were connected to a PC via QMODEM (Telecommunications Program Version 2.0, The Forbin Project Inc.) interfaced with the PT-100s. They were continuously sending the temperature signal to the PC at a frequency of 10 Hz. 3. Impedance Void Needle Probe 3.1. Probe Design. The probe for the study was developed using the design methodology of Chanson,24 who fabricated an
Figure 3. Block diagram for the electronic circuit used in the measurement.
impedance-based needle probe to measure local void fraction and bubble frequencies in an air-entrainment experiment. The advantages of the needle probe are its simple design and easy installation. It consists of two concentric conductors; the schematic is shown in Figure 2. The inner conductor (copper wire, conductor 1) was excited with an ac signal of 500 kHz and (5.0 V, while the outer one (stainless steel needle, conductor 2) was grounded. The gap between the wire and the needle was filled with an insulating material, such as epoxy resin, that could withstand temperatures of 200 °C. It was found that the impedance across the two electrodes, which increased with void fraction, was mainly resistive in the case of excitation frequencies below the order of MHz. Above this, the impedance was predominately capacitive. The output signal from the probe was led to a low-pass filter with a cutoff of 40 kHz and was demodulated to the dc signal being proportional to the local void fraction. The dc signal was then analyzed with the help of a DSO. The analysis is described later with the help of Figures 6 and 7. Figure 3 shows the block diagram of the electronic circuit used in the measurement. The method of Ishii,28 based
Ind. Eng. Chem. Res., Vol. 48, No. 4, 2009 2023
Figure 6. Direct current output from probe Lm stored in DSO.
Referring to Figure 5, one may write Figure 4. Installation scheme of the probe in the riser.
p1 ) pa + FmLlg
(2)
p2 ) pa + FmLrg
(3)
Subtracting eq 2 from eq 3 gives ∆p ) p2 - p1 ) Fm(Lr - L1)g
(4)
By neglecting the contribution of the steam bubbles to the total pressure drop ∆p as per the above assumption, ∆p may be defined as ∆p ) (1 - R)F1Lsg
(5)
L1 + Ls ) Lr + Lm
(6)
Combining eqs 4-6, one may obtain R)1-
Figure 5. Inverted manometer setup for the reference void fraction.
on Eulerian time-averaging, was applied to calculate the experimental void fraction. According to the method, the timeaverage void fraction R(z,τ) at the point z for the duration of τ is defined by
∫
1 R(x,τ) ) 0 τMk(z,t) dt (1) τ Mk(z,t) is state density function of the kth phase at the distance of z from the reference point and time t. τ is the integration time. The value of Mk(z,t) is either 0 or 1, depending on whether the kth phase is present (the conditioned output signal exceeds a chosen threshold voltage Vr) or absent (the conditioned output signal does not exceed Vr). It indicates whether the probe is registering the presence of the kth phase. 3.2. Probe Installation. The installation scheme is shown in Figure 4. First, the probe with its lead wires was inserted into a PVC tube and was tightly fitted to the tube by Teflon tape. Then, the PVC tube containing the probe was led to the center of the riser through a copper tube. The copper tube was a 90° bend and its axis of one section was aligned with the centerline of the riser. The copper tube was inserted into the riser by perforating the wall of the riser. The two pipes and the gap between them were well-sealed using a mechanical seal to prevent any leakage of fluid. 3.3. Probe Calibration. The probe was calibrated with the reference measurement of the local void fraction Rm using an inverted U-tube manometer. Rm was calculated using the principle of hydrostatics and the assumption of a negligible density of the steam phase29 at low pressures.
Fm(Ls-Lm)g FlLsg
(7)
In the case of an inverted manometer with the assumption of the absence of steam bubbles in the liquid column, the mixture density may be expressed as Fm)Fl
(8)
From eqs 7-8, the expression of R is R)
Lm Ls
(9)
In the present experiment, eq 9 was used to calculate the reference value of the local void fraction from the values of manometer reading (Lm) and the spacing of the manometer pressure tapping (Ls). In the present experimental investigation, Ls was 15 cm. The needle probe was placed at the midpoint of the spacing. Manometer reading Lm was recorded at various void fraction levels. During calibration, output signals were recorded at various pseudo-steady-state conditions of evaporator inlet temperatures, loop resistances, and heating rate. ∑k factors for loop resistances were varied from 12.89 to 14.29, heater inlet temperature from 45 to 61.0 °C, and heating rate from 1.2 to 2.0 kW. This is worthwhile to mention that choosing the proper value of Vr is a key to the calibration. It was observed that for a change in Vr of 0.05 V, the change in Rp was found to be approximately 75%. The change was 110% with respect to Rm. Table 2 presents the calculation of void fraction from the variation of dc voltage signal shown in Figure 6 over the time scale of 1.0 s and using eq 1 at a liquid temperature of 61 °C at the inlet of the heated section and heating rate of 2 kW. The dc voltage vs time plots in small time intervals, as shown in Figure 7a-e, reflect some measurement noises (voltage drops are stochastic in nature). In the figures, although the voltage drops indicate the presence of bubbles, all the drops have not
2024 Ind. Eng. Chem. Res., Vol. 48, No. 4, 2009
Figure 7. Output dc voltage over time intervals.
been counted for the calculation of void fraction. Two voltages have been chosen in the present study to determine Vr for which Rp becomes close to Rm. The calculations indicate that increasing the threshold voltage lowers the void fraction because of fewer output signals exceeding the higher threshold voltage. It is also evident from Table 2 that the choice of a proper threshold voltage depends on Rm. Table 2 shows that the probe indicates the void fractions Rp of 0.56 and 0.98 at the threshold voltages of -0.65 and -0.6 V, respectively. Lm, corresponding to the output signal from the probe (Figure 6), is 7.0 cm and the void fraction Rm ()Lm/Ls) becomes 0.467. The corresponding void fraction of 0.98 seems to be unreasonable, whereas the void fraction of 0.56 is close to Rm. Figure 8 shows the calibration of the probe at Vr of -0.65 V. Statistical analysis based on analysis of variance (ANOVA) indicates a strong correlation between Rm and Rp with the coefficient of 0.9742, a confidence limit of above 95%, and a standard error of 0.0438. Table 3 presents the summary of ANOVA analysis. Figure 8 also reveals that at low void fraction, the probe measures the lower values compared to reference values and at high void fraction, the probe gives higher values compared to reference values. This may be explained by the fact that in the range of high void fractions, bubble-slugs are expected to form and the threshold voltage needs to be tuned.
Its magnitude will be slightly higher than 0.65. It was observed from calculations of Rp that using a higher magnitude of the threshold yielded lower values of void fractions. Figure 9 shows the temperature history of a two-phase mixture at the riser exit during a measurement of void fraction. The average amplitude of the temperature oscillation was 1.3 °C. It is required to mention that the calibration data were generated at the parameter constellations of inlet subcooling and heating rate, in which steady or quasi-steady flow conditions were achieved. In the range of low void fractions, the bubble detection by a single probe stationed at the center of the ARI probably will not be accurate. At low bubbling conditions (at high inlet subcooling), the bubbles are localized at the vicinity of the walls of the heated section. It is anticipated that multiple numbers of the probes along axial and radial directions will improve the quality of the measurement. 4. Thermohydraulic Models for Void Fraction Dynamics Basic mathematical formulations of two-phase flow dynamics based on conservation principles called field equations are well described by Wallis30 and Ishii.28 More rigorous formulations are also found in the literature.31,32 Lahey and his research group33-35 extensively carried out modeling work for studying
Ind. Eng. Chem. Res., Vol. 48, No. 4, 2009 2025 Table 2. Calculation of Void Fraction Using the Probe at Different Threshold Voltages at Heating Rate of 2 kW and Heater Inlet Temperature of 61.0 °C threshold voltage ) -0.65 V
threshold voltage ) -0.6 V
time (s), t
∆t
M(x,t)
∆tM(x,t)
M(x,t)
∆tM(x,t)
0.01 0.06 0.08 0.1 0.13 0.15 0.17 0.24 0.3 0.35 0.39 0.44 0.51 0.56 0.62 0.71 0.78 0.82 0.85 0.87 0.91 0.98
0.01 0.05 0.02 0.02 0.03 0.02 0.02 0.07 0.06 0.05 0.04 0.05 0.07 0.05 0.06 0.09 0.07 0.04 0.03 0.02 0.04 0.07
0 1 1 0 0 0 1 0 0 0 1 1 1 0 1 1 1 0 0 1 0 1
0 0.05 0.02 0 0 0 0.02 0 0 0 0.04 0.05 0.07 0 0.06 0.09 0.07 0 0 0.02 0 0.07 0.56 0.56
1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
0.01 0.05 0.02 0.02 0.03 0.02 0.02 0.07 0.06 0.05 0.04 0.05 0.07 0.05 0.06 0.09 0.07 0.04 0.03 0.02 0.04 0.07 0.98 0.98
void fraction ) ∑M(x,t)∆t τ
Figure 8. Calibration of the needle probe at Vr of -0.65 V. Table 3. ANOVA Analysis for Regression regression statistics multiple R R2 adjusted R standard error observations
0.974 291 101 0.949 243 15 0.942 898 544 0.043 829 151 10
the dynamics in a typical boiling channel and BWRs. The unique feature of the model is that it essentially incorporates a scheme of moving boiling boundary, which renders a moving boundary problem as a whole for the systems. The Eulerian-Eulerian approach of Ishii has been adopted in the present study for the discritization of PDEs evolved during the formulation of the problem. According to this approach, both liquid phase and vapor phase behave separately as a continuous phase. The continuity formulation (local one or instantaneous one) of a quantity being conserved is usually done through either spaceor time-averaging of the quantity or averaging both. This leads to the equations for two-phase flow in terms of space-averaged or time-averaged or both space- and time-averaged flow variables. The proposed mathematical analysis of void fraction dynamics has been carried out with moving-boundary formulation of boiling two-phase flow in a natural circulation loop along with due consideration of more realistic dynamics of heated tube wall temperature.26,36 The essential feature of the formulation (five-equation drift flux) is that it eliminates the assumption
Figure 9. Temperature oscillation at the riser exit during a run.
of vapor temperature or liquid temperature at the state of saturation (as is the common practice in the four-equation scheme). Computation of both vapor temperature and liquid temperature in the five-equation scheme has been carried out by considering two energy equations. In addition, the proposed analysis has been carried out with the help of the finite volume method, instead of the finite difference method, since the former technique has been identified to be more befitting in the present case.37 A five-equation model involving conservation equations of mass, momentum, and energy of a mixture, mass conservation equation of vapor, and energy conservation equation of liquid is presented below with the following simplifying assumptions: (1) flow is one-dimensional, (2) the average of products is approximated by the product of averages over cross section, (3) heat and momentum conduction in phases and viscous dissipation are neglected, (4) evaluation of thermophysical properties is based on uniform pressure throughout the selected cross-section, (5) ARI and down comer are adiabatic, (6) heat loss is neglected, (7) dynamic behavior of heated tubes has been represented by the model of a single tube, (8) the pressure at the inlet of the down comer is the pressure in the steam-water separator (system pressure) and it is constant, and (9) there is uniform heat flux throughout the heated section. One may note that the energy balance in the present case considers only the thermal energy equation, which can be obtained by subtracting the equation of kinetic energy from the total energy equation consisting of internal energy and kinetic energy. This was elegantly discussed by Bird et al.38 The contributions of potential energy and kinetic energy have been neglected because the present investigation is focused only on heat transfer phenomenon and their magnitudes are insignificantly low compared to internal energy. 4.1. Conservation Equations. The finite volume forms of the conservation equations applicable for boiling two-phase flow in the heated section presented in Figure 10 are given by eqs 10-14. It may be noted that in Figure 10 the boiling channel consists of a heated section (four tubes) and an ARI section. It may also be noted that a two-phase (2-φ) mixture exists above z ) λ. Liquid phase mass balance:
∫
λ
Lj
∂ [(1 - R)Fl]j dz + ∂t
∫
Lj
λ
∂ [(1 - R)FlVl]j dz ) ∂z
∫
λ
Lj
-Γg dz (10)
Vapor phase mass balance:
∫
∫
∂ ∂ [RFg]j dz + λ Lj [RFgVg]j dz ) ∂t ∂z Two-phase mixture momentum balance: λ
Lj
∫
λ
Lj
Γg dz
(11)
2026 Ind. Eng. Chem. Res., Vol. 48, No. 4, 2009
∫
λ
∫
∂ ∂ [(1 - R)FlVl + RFgVg]j dz + λ Lj [(1 - R)FlVlVl + ∂t ∂z RFgVgVg]j dz ) ∂p - λ Lj + g[(1 - R)Fl + RFg] sin θ + ∆pfT dz (12) ∂z Two-phase mixture internal energy balance: Lj
{
∫
∫
λ
∫
∂ ∂ [(1 - R)Flhl + RFghg]j dz + λ Lj [(1 - R)FlhlVl + ∂t ∂z RFghgVg]j dz ) ∂p ∂pj j Lj dz + λ Lj[(1 - R)Vl + RVg]j dz + λ LjQj dz (13) λ ∂t ∂z Liquid phase internal energy balance: Lj
∫
∫
λ
}
Lj
∫
∫
∫
∂ ∂ [(1 - R)Flhl]j dz + λ Lj [(1 - R)FlhlVl]j dz ) ∂t ∂z ∂pj ∂pj Lj dz + λ Lj[(1 - R)Vl]j dz + λ LjQlj dz (14) (1 - R) λ ∂t ∂z
∫
∫
∫
The exit conditions (temperature, void fraction, phase enthalpies) of a computational node equal to average conditions of the node (eq 18) has been solved implicitly to calculate φ. It is assumed that computational nodes behave as lumped systems. Thus, a set of ODEs obtained from PDEs given by eqs 10-14 for the 2-φ section have been solved simultaneously by Gear’s method. The method is illustrated in Section 5.1. 4.2. Dynamics of Heated Tube-Wall Temperature. The dynamics of the heated tube wall is highly important, particularly during transient conditions such as startup and shutdown to replicate the real-time complicated thermohydraulic phenomena that occur under the conditions. Calculation of the tube wall temperature for flow boiling in a tube has been done on the basis of the balance of thermal energy using the lumped parameter approach. The model has been formulated in reference to Figure 11, which shows the wall heat transfer scheme. The model is represented by the following equation: Fwcpw
qex′′Do 4Diqb′′ 4Dip1p(Tw - Tl) dTw ) 2 - 2 2 2 dt D -D D -D D 2-D2 o
where
i
o
i
o
i
(19) Γg ) Γb - Γc
(15)
Γb is due to the contribution of boiling and Γc is due to the contribution of bubble condensation. Qlj is contributed by bubble condensation and wall heat transfer to the liquid-phase by convection. Total wall heat transfer Qj consists of heat transfer due to boiling and convective heat transfer to the single-phase liquid. λ is the one-phase (1-φ) liquid height (as shown in Figure 10, locating the interface between and 2-φ sections of the boiling channel) above which the liquid starts to boil and is known as the boiling height. It is required to mention that the complexities lie in determining λ because of its variations with loop flow rate, enthalpy of liquid at the inlet of heated section, system pressure, and heat input during the startup of a SWNC. Thus, a moving boundary problem evolves while computing λ and flow variables in both 1-φ and 2-φ sections. Computation of λ has been done by coupling the liquid-energy balance equation and continuity equation in the 1-φ section. The resulting equation has been integrated from z ) 0 to z ) λ based on finite volume method, as given below.
∫
∫
λ(t) ∂h
∫
∂h dz + 0 λ(t)Vl dz ) 0 λ(t)qw dz (16a) ∂t ∂z In the above equation, qw is the volumetric heat transfer rate to the 1-φ liquid zone and it is in W/kg. In the 1-φ zone of N number of computational grids, eq 16a results in a set of ODEs (lumped-parameter model) when the Leibnitz rule is applied.35-37 The general form of the ODEs is given by the following equation for the integration from z ) 0 to z ) λ: 0
d〈h〉 dλ + [〈h〉 - hλ] ) Vlh|0 - Vlh|λ + 〈qw 〉λ (16b) dt dt The derivation of lumped-parameter models from the distributed parameter models in 2-φ section has been made by the above method. The lumped-parameter models are represented by eq 18 obtained with a general distributed-parameter model given by eq 17: λ
∫
L ∂φ
dz +
∫
L ∂(Vφ)
dz )
∫
L
S dz
(17)
d〈φ〉 〈φ〉 - φλ dλ Vφ|λ - V|L〈φ〉 ) + + 〈S〉 dt L - λ dt L-λ
(18)
λ(t)
∂t
λ(t)
∂z
λ(t)
The last term on the right-hand side of eq 19 is due to the contribution of heat convection to a single-phase liquid during flow boiling. It is required to mention that eq 19 has also been used to compute the wall temperature in the 1-φ section by equating qb′′ to zero. Figure 11 shows the wall heat transfer scheme on which the formulation of the problem has been made. The external heat flux qex′′ is due to the steam-heating (steam being condensed outside the tube) of the tube wall. The heat flux for each of four tubes is calculated by the heat exchanged from the steam side divided by four divided by the heated surface area of each tube. qex′′ is given by qex′′ )
Hex 4(πDoLH)
(20)
LH is the heated length. Hex is external heating rate (W). In the present experiment, Hex has been estimated from the enthalpy change of heating steam corresponding to the measured temperature change and its measured condensation rate. The pressure drop on the shell side has been neglected. Substituting the expression of qb′′ (nonlinear function of tube-wall temperature) in eq 19, one obtains a stiff ODE that computes Tw. The stiff ODE has been solved using multistep Gear’s methods.39 4.3. Constitutive Relations. The relations include the correlations for Fanning’s friction factors for single-phase flow, the two-phase multiplier for two-phase flow, heat flux for boiling and bubble condensation, and heat transfer coefficients for single-phase convection. Table 4 presents correlations for friction factors and the two-phase multiplier, while Table 5 presents correlations for heat flux and heat transfer coefficients, as used in the present investigation. The drift-flux correlation of Ishii48 for void distribution coefficient and drift velocity, as used in this work, is presented in Table 6. 5. Solution Scheme for Model Equations The solution schemes of field equations for two-phase flow are well-addressed in the literature, in which implicit and semiimplicit schemes using finite difference method have been detailed.49-51 In the present work of moving boundary analysis, the FVM scheme followed by multistep Gear’s method for stiff ODEs has been proposed. 5.1. Overall Solution Scheme. For the overall solution scheme for a large set of model equations, the modular approach
Ind. Eng. Chem. Res., Vol. 48, No. 4, 2009 2027
hl) for single-phase liquid flow. (Single-phase zone and twophase zone are separately solved using the implicit scheme of Gear’s method. The stiff ODE set is obtained for a computation grid by integrating finite-volume forms of conservation equations using Leibnitz’s rule. For two-phase flow, eqs 10-14 have been solved simultaneously); (6) solution of model equations (equations of single-phase liquid mass, momentum, and energy) of the down comer for (Vl, p, hl); (7) solution of loop momentum balance; and (8) calculation of local mass flow rate and velocities of phases. The following ODE set is obtained for a node volume by integrating conservation equations (eqs 10-14) using Leibnitz’s rule. From liquid-phase mass balance: Figure 10. Schematic of experimental boiling channels.
dp ) fl(hg, Jm, p, R, hl) dt From two-phase momentum balance:
(21)
dJm ) f2(hg, Jm, p, R, hl) dt From two-phase mixture energy balance:
(22)
dhg ) f3(hg, Jm, p, R, hl) dt From liquid-phase energy balance:
(23)
dhl ) f4(hg, Jm, p, R, hl) dt From vapor-phase mass balance:
(24)
dR ) f5(hg, Jm, p, R, hl) (25) dt Within a time-step, the above ODE set (eqs 21-25) have been solved simultaneously iteratively till a desired convergence is attained. The unknown vector Y is updated in an iteration using the implicit Euler method as shown by eq 26. Yk+1 )
[
+ ∆tF ( ∂F ∂Y ) ] ∂F I - ∆t( ) ∂Y
Yk I - ∆t
k
k
(26)
k
Figure 11. Schematic of wall heat transfer in a single-tube boiling channel. Table 4. Correlations for Friction Factors and the Two-Phase Multiplier correlations
friction factor
Hagen-Poiseuille fl ) 16/Re law 40,41 for laminar flow Blasius equation fl ) 0.079/Re for turbulent flow42 Becker et al.43 -
two-phase friction multiplier -
0.25
φ2Lo ) 1.48623 × 108(x/p)0.96 + 1
has been adopted. The scheme has been presented in detail by Paruya and Bhattacharya36 for startup simulation of a nuclear steam generator. The sequential steps are presented below for showing the computational algorithm in new time step: (1) computation of given boundary conditions (heat flux) and given initial conditions at the previous time step; (2) calculation of heated-wall temperature, boiling heat flux, and convective heat flux to a single phase; (3) solution of boiling boundary in a heated tube; (4) calculation of Γg; (5) solution of model equations (equations of mass, momentum, and energy) in boiling channel for (Jm, p, hl, hg, R) for two-phase flow and for (Vl, p,
where Y ) [hg Jm p hl R]T for two-phase steam-water flow through a channel, [Vl p hl]T for single-phase water flow through a channel, [p R]T for a two-phase mixture in a tank (the original model equation is ODE), or [p hl]T for a single liquid in a tank (the original model equation is ODE); F ) [f1 f2 f3 f4 f5]T for two-phase steam-water flow through a channel, [f1 f2 f4]T for single-phase water flow through a channel, [p() R()]T for a twophase mixture in a tank, or [p() hl()]T for a single liquid in a tank; and (∂F/∂Y)k ) the Jacobian matrix after the kth iteration. The methods for automatic time-step selection that compromises solution accuracy, solution stability, and CPU time have been discussed in the literature37 and have been implemented for model computations. 6. Validation of Simulation Results It was found experimentally that it is very difficult to study the geysering instability phenomena in the real-time process because of its metastable nature over narrow ranges of heater inlet temperature, loop resistance, and heating rate, given the system configuration under investigation. Even a slight change in each of loop resistance, heater inlet temperature, and heating rate outside the range forces the unstable conditions to reach
2028 Ind. Eng. Chem. Res., Vol. 48, No. 4, 2009 Table 5. Boiling Heat Transfer Coefficients or Heat Fluxes and Single-Phase Convective Heat Transfer Coefficients correlation
heat transfer regimes
p orq′′
Collier, Re < 2000 Sieder et al.,45 Re > 3000
single-phase convection single-phase convection
Thom et al.46 Lahey47
nucleate boiling (subcooled or bulk) bubble condensation
44
pD/k ) 0.17Re Pr0.43(Pr/Prw)0.35Gr0.1 pD/k ) 0.023Re0.33Prm(µ/µw)0.14 m ) 0.3 for cooling m ) 0.4 for heating qb′′ ) 1.0945 × 103 exp(p/630)(Tw - Tsat)2 qc′′ ) 0.075hfg(Fl-Fg)DhR(Tsat - Tl) 0.33
Table 6. Ishii Correlations for C0 and Vgj flow regime bubbly flow (0 < R
