Nonlinear Analysis of Chaotic Time Series in a Natural Circulation

Article pubs.acs.org/IECR

Nonlinear Analysis of Chaotic Time Series in a Natural Circulation Boiling Loop Arnab Karmakar and Swapan Paruya* Chemical Engineering Department, NIT Durgapur, West Bengal, India-713209 ABSTRACT: The present study reports chaotic flow oscillations observed in a natural circulation boiling loop. The periodic oscillations of wall temperature (thermal oscillations) initiate at a certain condition of heat flux and inlet subcooling in addition to geysering instability and pressure-drop oscillations. We observed that boiling regime changes from nucleate boiling to transition boiling as a result of decrease in inlet subcooling and increase in heat flux. The thermal oscillations are strongly coupled with pressure-drop oscillations. Nonlinear analysis of the time series of loop flow rate at various heater power and inlet subcooling have been carried out using statistical analysis, fast Fourier transform (FFT), time delay embedding for attractor reconstruction, autocorrelation, and correlation dimension. The analysis confirms that the oscillations are more chaotic at relatively low heater power and high inlet subcooling. The complexity of the oscillations strongly depends on boiling heat transfer regime. Our observations and analysis have been supported by other relevant experiments.

1. INTRODUCTION Studies on flow oscillations in the natural circulation boiling loop (NCBL) in boiling water reactors (BWRs), drum-type boilers, thermosyphon reboilers in distillation columns and chemical evaporators, thermosyphon-based heat removal from chemical reactors, etc. have received significant attention in regard to safe and optimal operations. The natural circulation systems offer better operational safety compared to forced circulation systems; although the former are inherently more prone to several nonlinear flow instabilities (flow reversal, geysering, densitywave instabilities, etc.) because of low driving heads. Experimental studies at the conditions of low pressures and low heat fluxes report five major instabilities that occur in NCBL (see the works of Wissler et al.,1 Ledinegg,2 Chexal and Bergles,3 Fukuda and Kobori,4 Aritomi et al.,5 Wu et al.,6 Kyung and Lee,7 Jiang et al.,8 Zboray et al.,9 Hsieh et al.,10 Jain et al.,11 and Paruya et al.12). The instabilities include flow excursion or flow reversal, geysering instability, natural circulation oscillation, flashinginduced instability, and Type I (low steam quality) density-wave oscillation (DWO). Qualitative descriptions of these instabilities can be found in the review of Boure et al.13 Flow excursion or flow reversal is also known as Ledinegg instability.2 This static instability occurs when the sum of frictional pressure drop, gravitational pressure drop, and acceleration pressure drop decreases with mass velocity. For a given heat flux, with increase in mass velocity, the total pressure drop decreases to a minimum value at which the onset of Ledinegg instability occurs. In NCBL, geysering instability is a condensation-induced phenomenon with a periodic flow reversal that takes place as a result of the inception of bubbles due to liquid superheating. This instability was observed by many investigators5,8,14 to occur in a vertical forced-flow boiling channel and NCBL at low heat fluxes and low pressures. As boiling begins, vapor generated in the column transfers heat to the surrounding cold liquid bringing it to its boiling point to enhance bubble formation. The bubbles subsequently push the water in a jet form out of the column (unloading of column). The water jet is known as a geyser. The hydrostatic head decreases due to the water ejection causing a © 2012 American Chemical Society

vigorous boiling and the formation of the geyser enhances. The flashing refers to the vapor generation in an NCBL due to decrease in hydrostatic pressure (the saturation enthalpy decreases along the flow path).1,15 It usually occurs during the startup of NCBL with a long adiabatic riser at high heat fluxes and low pressures. The effect of flashing becomes significant at low pressures. The time period of the flashing-induced oscillations is usually observed to be 1.5−2.0 times of the fluid transit time (time required by the fluid to travel the length of adiabatic riser). DWOs are the results of multiple thermohydraulic feedbacks associated with the flow parameters such as flow, pressure, temperature or enthalpy, and density or void fraction.4,16 The fundamental characteristics of DWOs is that the pressure drop and the mass flow rate oscillate in out-of-phase mode. At a low quality or a low power, the pressure drop due to gravity dominates. A significant change in void fractions occurs for a small change in quality and so is the pressure drop due to the gravity. Consequently, the flow oscillates (low steam quality DWOs). At a high steam quality or a high power, the contribution of the gravity component of the pressure drop to the oscillations (high steam-quality DWOs) becomes less (friction component dominates) because of the fact that the slope of the void fraction versus quality curve reduces. The transportation delay in the channels also greatly contributes to DWOs of the system. The nonlinear analysis of the flow oscillations observed in NCBL is extremely important in order to accurately understand the complex dynamical behavior of the loop. Although the mathematical formulations based on the conservation of mass, momentum, and energy in the loop can predict the nonlinear behavior, the nonlinear analysis of experimental time series is capable of accurately identifying the complex nonlinear behavior Received: Revised: Accepted: Published: 16467

March 25, 2012 October 17, 2012 November 21, 2012 November 21, 2012 dx.doi.org/10.1021/ie300788t | Ind. Eng. Chem. Res. 2012, 51, 16467−16481

Industrial & Engineering Chemistry Research

Article

of NCBL. The analysis of loop flow rate, local pressures, local temperatures, and void fractions includes FFT (fast Fourier transform),17 time delay phase space reconstruction (TDPSR),18 autocorrelation,19 singular value decomposition (SVD),20 average mutual information (AMI),21 the estimation of correlation dimension D2,22 and Lyapunov exponent λ.23,24 A systematic nonlinear analysis of an oscillating time series is found in the literature.25−28 Wu et al.6 experimentally observed chaotic oscillations in a single-channel NCBL at a pressure P of 151.9875 kPa, a heater power Q of 4.0 kW, and inlet subcooling ΔTsub of 75 °C and identified four fundamental frequencies of walltemperature oscillations. Heish et al.10 identified chaotic oscillations of loop flow rates based on FFT and estimated D2 at a specific operating condition in NCBL. Cammarata et al.20 analyzed the experimentally obtained time series for a singlephase natural circulation flow in a rectangular loop using TDPSR, SVD, and AMI. They reported that D2 ranged from 2 to 3 depending on heat flux and it increased with heat flux. The irregular distribution of the trajectories in their reconstructed phase space revealed the existence of chaotic behaviors. Lakshmanan and Pandey29 numerically analyzed the chaotic time series generated from RELAP5 model of a natural circulation boiling system based on FFT, Lyapunov exponent, and Poincare section using the TISEAN package. Qu et al.30 observed chaotic time series of temperature in natural circulation pulsating heat pipe and determined D2 (2.11−4.77) and λ (0.037−0.065) of the temperature oscillations at several Q and filling ratios. It is also noted in their studies that the oscillations in the evaporation section are relatively more chaotic than those in the condensation section. Narayanan et al.31 systematically determined D2 of oscillations of inlet velocity for a forced-flow boiling channel based on their numerical model. In the survey of the literature, we realize that analysis of chaotic time series for a natural circulation boiling loop has not been systemically addressed in details. The parametric sensitivity of fractal dimension, e.g. D2, has not been studied as well. The study will enable us to get more insight into the mechanism of the oscillations and to derive a more realistic formulation. It is also noted that the mechanism of thermal oscillations in NCBL is not well understood. In the present paper, we analyze the chaotic time series of loop flow rate and pressure for an experimental NCBL at different inlet subcooling and heater power. Nonlinear analysis of loop flow rate Wloop has been carried out using statistical analysis, FFT for power spectrum density (PSD), TDPSR, autocorrelation, and the Grassberger−Procaccia22 algorithm for estimating D2. The effect of Q and ΔTsub on the nonlinear behavior of the boiling loop has been investigated. The mechanism of thermal oscillations experienced in our NCBL has been derived. The main contributions of the paper are that the complexity of the chaotic oscillations in NCBL has been quantified based on the estimation of fractal dimensions, and the same has been attempted to relate to boiling heat transfer regime (BHTR). Gradual change of BHTR from nucleate boiling to transition boiling has been demonstrated. We also attempt to relate the thermal oscillations to pressure drop oscillations. It is worthwhile to mention that BHTR has been identified based on PHI signal and Twall signal, visual observation, and picture motion browser (PMB) with the help of a Sony Handy Cam (HRD-XR160, 24 frames/s, 3.3 mega pixels). We now briefly describe the tools for nonlinear analysis used in the present analysis in the next section.

2. NONLINEAR ANALYSIS OF TIME SERIES 2.1. Fast Fourier Transform (FFT) for PSD. Fourier transform operates on a time series to transform it into a frequency spectrum (time domain to frequency domain or Fourier space). It shows how certain frequencies contribute to the oscillating signals. The technique is very powerful in determining the periodic nature of time series. In general, all oscillating piecewise stationary signals may be described by a sum of sine waves, as presented below. For a continuous system, ∞

1 2π

y(ω) =

∫−∞ y(t )e2πiωt dt

(1)

For a discrete system, yk̃ =

1 N

N

∑ sne2πikn/N ∫ n=1

∞

y(t )e 2πiωt dt

−∞

(2)

A power spectrum containing sharp pulses at regular frequencies indicates periodic oscillations while a continuous power spectrum indicates a chaotic or a quasi-periodic orbit. However, for a chaotic time series, Fourier transform can not reveal any deterministic features. It may be the starting point for the time series of interests to determine its aperiodic nature. 2.2. Time Delay Phase Space Reconstruction. One of the important parts of nonlinear analysis is the reconstruction of phase space (generally referred to as attractor) that represents the trajectory of the measured time series. The most interesting part of the attractor reconstruction is that it does not require any prior information about the underlying dynamics of a system. Thus, it resolves the difficulty of representing the attractor for a single measured time series. Such a good representation of the attractor gives an idea of the complexity of a system and eventually participates in the calculation of correlation dimension. Time delay embedding enables us to draw the same. Takens embedding theorem18 is usually applied to reconstruct the attractor. As adapted from the work of Hegger et al.,25 the theorem states that if {Y(i)} is the sequence of scalar measurements of the state of a dynamical system (Y(i) is expressed below), then under certain genericity assumptions, the time delay embedding provides a one-to-one image of the original set, provided embedding dimension m is large enough. Now, the ith embedding vector Y(i) is defined by Y (i) = {y(i), y(i + τ ), y(i + 2τ ), ..., y(i + (m − 1)τ )} (3)

Thus, Y(1), Y(2), ..., Y(i), ..., Y(N − (m − 1)τ) may be defined. τ is the time delay of the reconstructed attractor. In other terms, the attractor can be directly reconstructed from the lagged time series of the system for m ≥ 2D2 + 1. m signifies the minimum attractor dimensions in which the attractors are completely unfolded, and D2 is the minimum number of independent variables (degrees of freedom) one needs to accurately model the chaotic behavior of the system. However, τ and m may be suitably chosen based on the method of Leibert and Schuster32 and the method of false nearest neighbor of Kennel et al.,33 respectively. One may also use autocorrelation function A(τ) expressed below to estimate the time delay:34 N−τ

A(τ ) = 16468

∑i = 1 (δiδi + τ ) N−τ

∑i = 1 δi 2

(4)

dx.doi.org/10.1021/ie300788t | Ind. Eng. Chem. Res. 2012, 51, 16467−16481


Article

Figure 1. Schematic of our experimental natural circulation boiling loop.

selection of τ is essential because of the presence of measurement noise in the experimental time series of finite data set. If τ is chosen small, y(i + τ) does not significantly differ from y(i). On the other hand, very large τ beyond an optimal value may give rise to a weak correlation between y(i + τ) and y(i) because of loss of information. Thus, an optimal τ enables us to achieve a good reconstruction of the attractors while the underlying dynamics of the system is retained. 2.3. Correlation Dimension. The correlation dimension D2 is one of the measures of the fractal dimension of a chaotic system. D2 quantifies the complexity of the chaotic attractors. It basically dictates the number of ordinary differential equations required to describe the chaotic system. Computation of D2 consists of centering a hypersphere about ith point in hyperspace and then letting the radius (r) of the hypersphere grow until all n points are enclosed. D2 essentially requires the computation of correlation sum C(r) expressed by

where δi = yi − ymean and ymean is the mean value of the time series. Generally, the time delay is selected when A(τ) crosses 1/2 of the value for the first time in the A(τ)−τ plot. We now describe the method of reconstruction of phase space. For a measured time series of N number of samples {y(1), y(2), ..., y(N)}, the reconstructed attractors can be generated as Y (1) = {y(1), y(1 + τ ), y(1 + 2τ ), ..., y(1 + (m − 1)τ )} Y (2) = {y(2), y(2 + τ ), y(2 + 2τ ), ..., y(2 + (m − 1)τ )} ⋮ Y (i) = {y(i), y(i + τ ), y(i + 2τ ), ..., y(i + (m − 1)τ )} ⋮ Y (N − (m − 1)τ ) = {y(N − (m − 1)τ ), y(N − (m − 2)τ ), .... , y(N )} (5)

where, i = 1, 2, ..., N − (m − 1)τ and Nm = N − (m − 1)τ is the total number of points in the reconstructed phase space. Although Takens felt that there is no constraint for the selection of τ for an infinite number of noise-free data points, the proper

C(r ) = 16469

1 Nm 2

Nm Nm

∑ ∑ H(r − max|Yi − Yj|) i=1 j=1

(6)



Article

H(·) is the Heaviside function defined as ⎧1 x ≥ 0 H (x ) = ⎨ ⎩0 x < 0

Table 1. Dimensions of the Various Parts of the NCBL (7)

The determination of D2 requires the guess of m. Grassberger and Procaccia22 first proposed a computationally efficient algorithm to estimate it from chaotic data. Basically, a reference point Yi is randomly chosen from all the Nm (= N − (m − 1)τ) points Yi (i = 1, 2, ..., N − (m − 1)τ) as shown in eq 5. The distance of Yj from the reference point is computed. The maximum norm was used in our study. The maximum norm is easy to compute. With it, the scaling region was observed to remain almost constant in a specific range (rmin, rmax) as m was increased. With the Euclidian norm, we observed different scaling regions as m was increased. For deterministic systems, the correlation sum in eq 6 follows a power law:

C(r ) = r D2

ln C(r ) ln r

orientation

heated section (AB) adiabatic riser section (BC) upper plenum (CD) lower plenum (JK) condenser (EH) condenser tubes (4 nos., FG) down-comer (HI) condenser tube pipes including AB, BC, JK, EH, FG, and HI

vertical vertical horizontal horizontal vertical vertical vertical

diameter (mm)

length (mm) 1000 2000 1200 1200 600 450 2350

10.0 (ID) 20.0 (ID) 25.0 (OD)

the exit of the riser. The differential pressure transmitter (DPT) made of SS 310 with the range of 0−10 bar is placed to measure the pressure drop across the distance of 0.6m in the heated section. An electromagnetic flow meter (EFM) with 20 mm polytetrafluoro ethylene lining (the range of 0−1.8 m3/h and uncertainty of less than 0.5%) is placed at the middle of the lower plenum to measure Wloop. It may be noted that the EFM can accurately measure very low flow and reverse flow as experienced in the natural circulation loop. A two-phase valve is placed at the upper plenum to regulate the loop flow rate by changing the loop resistance. The NCBL system also has the cold-leg temperature controller (CLTC) for regulating ΔTsub, pressure controller (PC), and heater-power controller (HPC). Pressure, heater power, and heater inlet temperatures are the most important degrees of freedom of the loop, which have a coupled effect on the complex flow oscillations in the loop. The control of these quantities is very essential to accurately identify the complexity of the flow oscillations and to gain the knowledge on the effect of each of them on the oscillations. In order to maintain the desired temperature at the cold leg, CLTC manipulates the cooling water flow for condenser by changing the opening of the control valve placed on the discharge line of a centrifugal pump (Make: CG, 0.5 HP) that circulates cooling water from a tank of 500 L. CLTC was operated in proportional−integral (PI) mode. A solenoid valve placed at the line connecting the loop and the compressor is manipulated by the proportional (P) mode of a PID controller to regulate the loop pressure. For uniform and safe heating, six segmental heater-coils (Kanthal super: MoSi2) are wound around the aluminum tube of the heater length. The aluminum tube is coaxially placed outside the pipe to be heated. The heat is transferred from the aluminum surface to the pipe wall by the process of radiation and convection. An accurate control of the heater power is achieved thyristor-based HPC in PID mode. Water (prepared by the ultra quartz distiller) was used as the working fluid. The data logger and SCADA (supervisory control and data acquisition) were put online with the sensors and the controllers. It was observed at the atmospheric pressure that twophase flow incepted in the heated section when the cold leg temperature (after the condenser) attained a value of 65 °C. The cooling water pump was then run for the condenser. It was noted that CLTC performs better at higher ΔTsub. The experiments were carried out at different cold leg temperatures of 70, 60, and 50 °C (corresponding to ΔTsub of 30, 40, and 50 °C, respectively) for a given Q (4.0 and 5.0 kW) to experience the effect of inlet subcooling. At lower ΔTsub, the controller parameters need to be tuned repeatedly. It takes longer time

(8)

The exponent D2 is the parameter that is of the greatest interest to us and is generally estimated from the slope of ln C(r) against ln r over a linear region as

D2 =

parts

(9)

For a deterministic chaos, the slope of ln C(r) vs ln r converges to a constant value with increased m for a small r. The constant value is taken as correlation dimension. Eckmann and Ruelle35 derived an interesting criterion D2 ≤ 2log10 N for determining a realistic D2. In determining D2, the choice of τ and m plays a crucial role as the proper selection of those parameters facilitates complete unfolding of the chaotic attractors. The issue of the proper choice of τ and m should therefore be resolved in order to have a strong relation between unfolding of attractor and D2. As mentioned, a small τ makes the axes of m dimensions too close to each other and a large τ causes the loss of information between the axes. On the other hand, if m is too small, the delayed phase space cannot completely unfold the attractor and false nearest neighbors would occur. As m is increased above the required value, D2 will cease to increase attaining a saturated value. Further increase of m will not help us reveal more dimensions of the attractor.

3. EXPERIMENTAL SET UP AND PROCEDURE Figure 1 shows a schematic of our experimental two-phase NCBL. The NCBL mainly consists of six partsheated section, adiabatic riser, condenser, down comer, and upper and lower plenum. The length and diameter of various sections of the loop are mentioned in Table 1. The horizontal plenums (upper and lower plenums) connect the hot leg and cold leg of the loop. The hot fluid from the riser is cooled in the condenser by cooling water flowing on the shell side. The rectangular loop was fabricated of SS 316. Different parts of the loop are manifolded to fit temperature sensors and pressure sensors. Eight resistance temperature detectors (RTDs) T1−T8 are located at heater inlet, heater exit, riser section (three RTDs), condenser exit (cold leg), cooling water inlet, and cooling water outlet. RTDs are PT-100 (sensitivity of 0.384 Ω/°C and uncertainty of less than 0.5%) with the range of 0−400 °C. Three pressure sensors P1−P3 made of SS 310 (WIKA make, Germany, uncertainty of less than 0.25%; response time of 1 ms) are placed at the heater inlet, the heater outlet, and the top of the riser. The locations of T1−T8 and P1−P3 are shown in Figure 1. The loop pressure is controlled at 16470



Article

for tuning as the response of the condenser becomes more oscillatory at lower ΔTsub. The time series data of temperatures from sensors T1−T8, pressures from the sensors P1−P3, and loop flow rate were recorded from the instant of the initiation of the heater. The data-logger software (DAQ) in the computer records the data from all sensors via the SCADA panel at 0.5 Hz. It is required to mention that the two-phase valve controlling the loop resistance is fully open in all experiments. In the next section, we present our experimental results (time series of Wloop, pressure at heat inlet of heated section PHI, pressure drop across the heated section ΔPH).

4. EXPERIMENTAL RESULTS In this section, we give a qualitative explanation to the time series of Wloop, PHI, and ΔPH obtained at varied conditions of Q and ΔTsub. One may note that the time series of 1800 samples are taken up for discussion from the instant when a sustained oscillation was observed. At the initial instant, the peak values of Wloop do not repeat (we do not consider them to be sustained oscillations). Wloop peaks repeat when a fine control of the temperature and the power is achieved. These two things were carefully maintained in order to achieve the sustained oscillations in all experiments. Another aspect is overall loop pressure. Although all the experiments were carried out at around atmospheric pressure (111.325 kPa), slight fluctuations of the same were observed because of the inception of the bubbles. A separate pressure controller (PC) was employed to counter the situation. While experimenting on boiling flow in an NCBL, there must be an uncertainty with respect to instruments, measurement, and overall process. In this study, the instrument uncertainties are negligible; all the instruments including pressure sensors, RTDs, and EFM have uncertainties less than 0.5%. The standard deviations of loop flow rate for a number of sets of experiments at all conditions have been found to differ within ±10%. This uncertainty comes from the controllability of CLTC. At different ΔTsub, the fluctuations of cold leg temperature in CLTC are ±2.62 °C at ΔTsub = 30 °C, ±1.71 °C at ΔTsub = 40 °C, and 1.47 °C at ΔTsub = 50 °C. The relative uncertainties of loop flow rate (calculated based on standard deviations of the flow rates) are not more than ±4.5% while the maximum uncertainty is 10%. It indicates that the controller performs better at high ΔTsub. On the other hand, though the signals of loop flow rate are mostly homogeneous, a bit of nonstationary behavior is observed. For a pure homogeneous signal the standard deviation σn is independent of the length of time series. On the other hand, depending upon nonstationary signal, σn depends on the time series, and the plot of σn vs length of time series grows with certain slope. At Q = 4 kW and ΔTsub = 40 °C, the Wloop signal is observed to be very close to stationary (slope = 0.007) while it slightly deviates from the stationary behavior (slope = 0.018) at Q = 5 kW and ΔTsub = 40 °C. At higher Q, the data become slightly nonstationary. Figures 2 and 3 show oscillations of Wloop for 1800 samples at ΔTsub = 50, 40, and 30 °C and at Q = 4.0 and 5.0 kW, respectively. The oscillations are chaotic pulses of Wloop. They have a continuous spectrum of frequencies (the power spectra are presented later) with some dominant frequencies. Such chaotic pulses in general reduce as ΔTsub (effect of inlet subcooling) is decreased. It is also interesting to note that at relatively high ΔTsub, the reverse flow is more frequent and its extent increases as well. While exploring the reasons for the occurrence of the

Figure 2. Effect of inlet subcooling on loop flow rate at Q = 4.0 kW.

Figure 3. Effect of inlet subcooling on loop flow rate at Q = 5.0 kW.

extended chaotic flows at high ΔTsub, we experienced highly complex and nonlinear bubble dynamics that occurs. The dynamics includes nucleation, bubble growth, bubble departure, and collapse. The bubble growth is mostly the evaporation of microlayer around bubble surface. Thickness of the microlayer is the function of thermal diffusivity α and bubble growth time tg. Under the subcooled condition, the effect of evaporation and condensation acts simultaneously causing the dynamic fluctua16471



Article

Table 2. Statistical Analysis of the Time Series of Wloop, ΔTsat, PHI, and ΔPH Q (kW)

ΔTsub (K)

SD Wloop

SD ΔTsat

SD PHI

SD ΔPH

ρ Wloop−PHI

ρ Wloop−ΔTsat

ρ Wloop−ΔPH

4.0

50 40 30 50 40 30

55.00 59.76 96.74 61.18 84.86 19.91

0.63 0.83 0.80 1.24 1.15 0.27

1.92 1.79 2.04 3.58 3.03 1.40

2.02 1.38 1.44 2.40 1.67 1.29

0.435

0.653

0.093

5.0

temperature in the heated section decreases and bubbles collapse. The bubble collapse takes place in the riser where the liquid is still subcooled and thermal nonequilibrium exists. Video imaging of tracer elements and our direct observations through the view glass (located at the exit of the heated section) confirm that the condensation of bubbles takes place. The stagnancy of the tracer elements and their downward movement from the riser section to the heated section are observed. 4.1. Statistical Analysis of the Time Series. Statistical analysis of the time series of Wloop, ΔTsat, PHI, and ΔPH has been carried out in order to find a strong correlation of Wloop with the rest of the quantities and to find the scale of the corresponding oscillation amplitudes. Their standard deviations (SDs) are presented in Table 2. SD gives the measure of the amplitude of oscillations. At lower Q (4.0 kW), SD of Wloop decreases with ΔTsub while those of PHI and ΔPH seem to pass through a minimum value with increase in ΔTsub. At relatively higher Q (5.0 kW), with increase in ΔTsub, SD of Wloop seems to pass through a maximum value while SDs of PHI and ΔPH increase. The oscillation amplitude of Wloop decreases with ΔTsub for a given heat flux. A sharp decrease in SD of Wloop from 84.86 to 19.91 at Q = 5.0 kW is noted for a change of ΔTsub from 40 to 30 °C. This sharp change may be due to the change of heat transfer regime from nucleate boiling to transition boiling. Such a switchover may attribute to the fact that reverse flow is reduced because of the availability of high driving force (buoyancy force) for natural circulation at high heat flux and low ΔTsub. Direct visualization shows a slug-to-churn and churn-to-annular transition at high Q and low ΔTsub. Detailed nonlinear analysis of the time series presented later will confirm such a change of the heat transfer regime. Estimation of overall correlation coefficients ρ of SDs of Wloop with SDs of PHI, ΔTsat, and ΔPH (listed in Table 2) indicates that the oscillations of Wloop may better be correlated with the oscillations of PHI and ΔTsat (ΔTsat forms much better correlations with Wloop) compared to those with the oscillations of ΔPH. The oscillations mostly originate due to thermodynamic reasons (wall superheat, nucleation, bubble growth, and bubble collapse) rather than due to hydrodynamic reasons. Due to the thermodynamic origins of the oscillations in which the bubble growth and bubble collapse actively participate, the pressure variations are more pronounced. The large deviation of the correlation coefficients from unity is due to the chaotic phenomena. We also noted an improvement of ρ for SDs of Wloop−SDs of ΔPH for the increase of Q to 5.0 from 4.0 kW. It indicates some contribution of ΔPH to the oscillations of Wloop at high Q because of the increased fluid friction at relatively high vapor quality (steam quality). The hydrodynamic factor plays roles in the oscillations of Wloop at high Q, and the mode of flow instability seems to change from geysering type to DWOs or other hydrodynamic instabilities. The correlations of Wloop with PHI and ΔPH can also be observed for a small time scale in Figures 4 and 5 at Q = 4.0 and 5.0 kW. At higher Q, the oscillations of ΔPH influence Wloop oscillations. It is also noted in the figures

tions of bubbles size before departing from the heater surface. This results in bubble-surface oscillation leading to a complex microconvection of the liquid around the bubble. High-speed video and direct observation through view glass identified complex trajectories of the interface motion and the microconvection of liquid around the bubbles which are believed to have a large fractal dimension. On the other hand, we cite some relevant and exciting results of the investigations of Dhir and his research group36,37 on subcooled boiling to understand the mechanism of subcooled boiling oscillations. The increase of ΔTsub reduces the bubble growth rate and increases the bubble growth period. Reduced bubble growth occurs because of loss of vapor due to the condensation at the vapor−liquid interface. The bubble diameter at departure Dbd also reduces with ΔTsub, although the effect is not appreciable. Heat loss due to the condensation of vapor is negligible during the early stages of bubble growth. However, in the later stages, as the interfacial area protruding into the subcooled liquid becomes large, the heat loss due to condensation increases nearly linearly with time. As the bubble base shrinks, the heat loss just before bubble departure becomes almost equal to the heat required for evaporation. They also noted in their experiments that the vapor-production rate (the number of active nucleation sites per unit surface) or the bubble growth rate increases with wall superheat ΔTsat. Consequently, tg decreases and Dbd becomes large. The relative effect of tg on Dbd is significant compared to that of ΔTsat. An experimental study of Kamil et al.38 on boiling incipience in natural circulation reboiler reports that the nucleate boiling is significantly influenced by ΔTsat, and it is associated with the formation of large turbulent eddies. Interestingly, the combined effect of the bubble waiting period tw and bubble growth period tg results in the frequency of bubble release from the heated surface, as given by the following equation. f=

1 tg + tw

(10)

However, tg is the function of both ΔTsat and ΔTsub whereas tw is only a function of ΔTsat given below. tw = 139.1(ΔTsat−4.1)

(11)

We note the high peaks of loop flow and reverse flow in the figures. Direct visualization and the use of high speed camera in the heated region indicate that the bubbles are formed in the region near the exit of heated region. The phenomenon lasts for a relatively long time and the fluid in the heated region then becomes saturated or near-saturated leading to a violent vapor generation. The flow pattern in the heated region is found to be slug or churn flow from bubbly ones at lower ΔTsub. The slugs and churns expel the liquid in an upward direction. As a result, the cold liquid from the down comer (cold leg) is sucked to the heated section, for which the flow meter registers a positive high flow rate. Due to this high flow of cold liquid, the liquid 16472



Article

Figure 4. Effect of inlet subcooling on oscillations of Wloop, PHI, and ΔPH at Q = 4.0 kW (black line for Wloop).

Figure 6. Power spectrum of Wloop at various ΔTsub for Q = 4.0 kW.

Figure 5. Effect of inlet subcooling on oscillations of Wloop, PHI, and ΔPH at Q = 5.0 kW (black line for Wloop).

that as Q increases or ΔTsub decreases, oscillations of PHI change from chaotic ones to the quasi-periodic.

5. NONLINEAR ANALYSIS OF EXPERIMENTAL TIME SERIES In this section we apply the nonlinear analysis tools on the time series of Wloop to characterize the oscillations at varied conditions of Q and ΔTsub. We present the time series in the form of power spectra using FFT and reconstructed attractors using time-delay embedding theorem. We have been able to establish that the oscillations are strongly chaotic with an appreciable D2 and deterministic in nature as the saturation of D2 with m has been noted. 5.1. Power Spectrum of Loop Flow Rates. The PSD of Wloop using FFT has been determined using the function f f t available in the MATLAB toolbox. The spectra are shown for 1800 samples in Figures 6 and 7 for Q = 4.0 and 5.0 kW, respectively. The dominant frequency decreases with ΔTsub for a given Q. A sharp change is noted when Q is increased from 4.0 to 5.0 kW. A distribution of fundamental frequencies (0.019−0.048 Hz at ΔTsub = 40 °C and 0.002−0.063 Hz at ΔTsub = 30 °C) with

Figure 7. Power spectrum of Wloop at various ΔTsub for Q = 5.0 kW.

dominant frequencies is observed at lower Q (4.0 kW). Interestingly, Karmakar et al.39 noted that with less number of samples, the power spectra of PHI and Wloop at Q = 4.0 kW and ΔTsub = 50 °C have more distinct fundamental frequencies and allows us to easily identify the PHI peak and the Wloop peak with identical frequencies at lower ΔTsub, which has some physical significances. However, the presence of many frequencies indicates a wide distribution of different types (size and shape) of bubbles at lower Q. This is because of the fact that low wall superheat results in a wide distribution of bubble diameters at departure as the bubble growth is affected by low superheat and consequently, the low-amplitude oscillations of Wloop are more 16473



Article

while their strength (PSD) at high Q is found to be very high, as evident in Figure 8.

frequently observed between two large-amplitude oscillations at lower Q. The same can be noted in Figures 4 and 5 showing increased numbers of low-amplitude oscillations of Wloop at Q = 4.0 kW. The boiling dynamics seems to be highly complex at relatively low Q. At ΔTsub = 50 °C and Q = 4.0 kW, one dominant frequency is noted. It indicates that the bubbles are of almost same size. It is noted that bubble frequency is less at a relatively high ΔTsub and bubble growth is uniform. At higher Q (5.0 kW) (refer to Figure 7), the fundamental frequencies are more distinct compared to those at Q = 4.0 kW; mostly 2−3 types of bubbles seem to be present because of 2−3 fundamental frequencies with appreciable PSD noted in the power spectra. Analysis by PMB for our high-speed video confirms that the small bubbles near the heated wall travel much faster than the large bubbles at the center. As shown in Table 3, the time period of oscillations based

Figure 8. Fundamental frequency and corresponding PSD at various ΔTsub.

Table 3. Fundamental Frequencies and Their Harmonics at Various Q and ΔTsub Q (kW)

ΔTsub (°C)

4.0 kW

50 40 30 50 40

5.0 kW

30

fundamental frequency, f f (Hz)

dominant frequency (Hz)

time period (s)

boiling delay, tbod + bubble transit time, tbt (s)

0.036 0.019−0.048 0.002−0.063 0.021, 0.042 0.0186, 0.039, 0.0744 0.003, 0.0075, 0.0177, 0.0466, 0.0613

0.036 0.048 0.063 0.042 0.039

27.77 20.83 15.87 23.83 25.64

27.41 23.75 18.90 23.94 20.70

0.061

16.39

16.97

5.2. Attractor Reconstruction. Reconstruction of phasespace trajectories (called attractors) from our experimental time series gives us valuable information on the chaotic character of the data. It evolves during nonlinear analysis of our data based on Takens delay time embedding theorem. As mentioned in Section 2, τ and m may be suitably chosen to unfold the trajectories as far as possible. We apply the method of Leibert and Schuster32 (elegantly described by Narayanan et al.31) implemented in FORTRAN to evaluate the optimal value of τ . The correlation sum C(τ) is calculated based on eq 6 for different values of τ and m with small r using our delay subroutine. C(τ) is then plotted against τ for several m to find the first minimum. Fraser and Swinney21 established based on the idea of mutual information introduced by them that y(i + τ) have the largest amount of additional information with respect to y(i) at the first minimum while they maintain a strong correlation. The corresponding τ at which the first minimum appears is accepted to be the optimal one. We made a reasonable choice of m for which the two consecutive plots look alike. In the present investigations, we varied delay time from 2.0 to 10.0 s and m from 3 to 10. The value of r (radius of hypersphere) was taken to be 1.0−4. The plots of C(τ) vs delay time-to-sampling time ratio are shown in Figure 9. Note that C(τ) is computed based on the normalized Wloop. The figure reveals that the first minimum corresponds to delay-tosampling time ratio of 3.0 (τ = 6.0 s), and the plots at m ≥ 9 are almost similar at all ΔTsub and Q. m ≥ 9 gives the feel that the saturation of D2 requires a maximum value of m equal to 9. Narayanan et al. found the ratio to be 4.0 for r = 0.005. This also indicates the chaos to be deterministic in nature. As the choice of τ is crucial in regard to unfolding the attractor and accurate estimation of D2, it is important to compare the estimated delay time with those evaluated for the related bubbling systems (refer Table 4). It is understood that as the boiling systems are inherently delayed system (a boiling delay is present) compared to air−water bubbling system, they have a higher τ. Although some studies assumed τ for their time series in NCBL,10,29 the method of Leibert and Schuster exercised in this study gives τ very close to theirs. One may note that the method may not necessarily give the required τ in all situations; further manual adjustment of τ may be required in order to achieve a good reconstruction.27 The method guides us to obtain a suitable τ. Figure 10 presents the reconstructed attractors of the original time series based on time delay embedding with τ = 6.0 s at different ΔTsub for Q = 4.0 kW. Note that the reconstructed attractors are not unfolded completely. It may be due to very complex chaotic oscillations as encountered in the experiments. The effect of ΔTsub on the reconstructions is also noted. The

on the dominant frequency is found to decrease with decrease in ΔTsub at 4.0 kW. It is also important to compare the time periods of the oscillations with the sum of boiling delay (tbod) and bubble transit time through riser (tbt) presented in Table 3. In our experiments, tbod and tbt have been estimated based on the following expressions.5

tbod =

tbt =

dlcpl ΔTsubaHL H 1000Q

LR u ̅in + ugj

⎛ g ΔdDR ⎞1/2 ugj = 0.35⎜ ⎟ ⎝ dl ⎠

(12)

(13)

(14)

tbod is the time required for the boiling of subcooled liquid. At Q = 4.0 kW, the time periods of oscillations are almost equal to (or slightly less than) the sum of tbod and tbt, indicating the occurrence of geysering instability. At Q = 5.0 kW, the time periods of oscillations are greater than the sum of tbod and tbt by the factor of 1.98−2.68 for the low-frequency peak (for the largefrequency peak, the factor is very close to unity), revealing the fact that the geysering instability gets appreciably affected by the inception of PDOs or DWOs; it becomes twice the sum tbod + tbt at ΔTsub = 50 °C and Q = 5.0 kW. It is also noted that PSD of lowfrequency peak increases with decrease in ΔT sub . The observations are in close agreement with those of Aritomi et al.5 who compared geysering period and boiling delay. They also noted the period to be twice the boiling delay in PDOs or DWOs. Geysering has high amplitude compared to DWOs. The dominant frequency at low Q is lower than that at high Q, 16474



Article

Figure 9. Plots of C(τ) vs delay time-to-sampling time ratio.

Table 4. Delay Time in Bubbling Systems SI no.

27

1

Vazquez et al.

2

Lin et al.28

3 5

Heish et al.10 Lakshmanan and Pandey29 present study

4 a

researchers

systems and methods air−water bubbling system, AMI air−water bubble column NCBL NCBL NCBL, method of Leibert and Schuster

delay time,a s 1.5 ( fs = 10.0 Hz) 1.0 ( fs = 50.0 Hz) 5.0 ( fs = 0.5 Hz) 6.0 ( fs = 0.5 Hz) 6.0 ( fs = 0.5 Hz)

fs = sampling frequency.

unfolding significantly increases (better reconstruction) with a decrease in ΔTsub. When ΔTsub is decreased, the temperature of liquid at the heater inlet approaches the saturation point and affects the boiling heat transfer regime. At high ΔTsub of 50 °C, the nucleate boiling regime is more prevalent. The complexity of boiling dynamics in the nucleate boiling regime discussed in section 4 is highest for which the reconstructed attractors have very complicated trajectories as shown in Figure 10. The nucleate boiling regime has been established to be significantly more complex so that higher embedding dimensions are required to unfold the attractor and a higher-order model is required to predict its dynamics compared to transition boiling and film boiling regime.40,41 Onset of nucleate boiling (ONB) has a fractal dimension higher than that of fully developed boiling. A better reconstruction obtained at lower ΔTsub = 30 °C can be explained by the fact that the nucleate boiling develops to be much better as a result of bubble merger facilitated along axial and radial direction of the heated section. This has been confirmed by our visual observation through view glass and high-speed video

Figure 10. Attractor reconstruction at Q = 4.0 kW and τ = 6.0s.

(Sony Handy Cam). Fortunately, this matches with the observations of Gaertner.42 Mushrooms of bubbles with short stems and vapor columns have been detected as they develop and 16475



Article

violent boiling has also been observed to occur. We have clearly explained how the bubble growth phenomena and the associated flow oscillations differ at various ΔTsub (50, 40, and 30 °C) at Q = 4.0 kW. The same can also be evidenced by the fluctuations of PHI signal presented in Figure 11. At ΔTsub = 50 °C, bubble-merger

Figure 12. Time evolution of wall temperature at various inlet subcooling at Q = 4.0 kW.

superheat at higher ΔTsub is affected by the initiation of pressuredrop oscillations (PDOs) across the heated section due to increased steam quality at lower ΔTsub. Vapor compressibility of the system with sufficient steam quality plays a key role in PDOs that occur in the section with negative slope of the characteristic Δp−velocity curve. We also note in Figures 4 and 5 that the oscillations of ΔPH become relatively stronger with decrease in ΔTsub. It seems that PDOs induce the thermal oscillations in the present experiments. There is not enough information on such thermal oscillations in NCBL in the existing literatures. In an attempt to relate thermal oscillations with transition boiling, Liu et al.43 reported the thermal oscillations of two types in their experiments on a vertical forced-flow boiling channel. The first one is a local phenomenon in which the high frequencies are attributed to the transition boiling while the second one is influenced by the thermal capacitance of the heater wall, axial conduction, and transition boiling characteristics. The second one is aggravated by the occurrence of the transition boiling in a hydrodynamically unstable system to give rise to sustained periodic oscillations of the heater wall temperature. These thermal oscillations are highly complicated ones due to the combined effect of thermal conduction, wall heat capacitance, and convective transition boiling. Liu and Kakac44 also reported that hydrodynamic instability triggers periodic thermal oscillations of low frequency in their experiments on forced-flow boiling channel. In our experiments, the appearance of PDOs and thermal oscillations also indicates that the boiling regime has switched to the transition boiling regime. This transition occurs between dryout (film boiling) and nucleate boiling. Our video shows that violent boiling occurs at almost regular interval leading to a highsteam quality flow at ΔTsub = 30 °C. Large vapor-slugs capture almost the entire diameter of flow pipe and a dryout condition appears. Microlayer becomes very thin and liquid-deficient at the heated surface because of vapor blanket. Consequently, the heat transfer from the surface to the fluid falls drastically and the surface temperature increases. However, the dryout condition is not stable, because the heated channel is hydrodynamically unstable (PDOs). The channel fills with liquid due to increased driving force for natural circulation. After a while, the nucleate boiling initiates causing the surface temperature to drop. Thus, the transition to dryout condition and nucleate boiling alternates. We obtain a very high period and low amplitude of the cycles in our experiments. High thermal capacitance of our system and less intense PDOs result in low amplitude of thermal oscillations compared to that observed by Liu et al.43 We also noted the periodic oscillations of surface temperature at Q = 5.0 kW and at

Figure 11. Effect of inlet subcooling on oscillations of pressure at heater inlet at Q = 4.0 kW.

leading to the formation of big bubbles (short slugs) and the involved fluid-expulsion phenomena occur with very irregular frequencies due to extended subcooling effect. The short slugs form infrequently to a small extent and the PHI-signal is more chaotic in nature. Note that the reduction of PHI is the indication of slug or big bubble formation. The fluid expulsion in the heated section (indicated by the peaks of loop flow rate in Figure 2) has been noted to last for longer duration. With decrease in ΔTsub, the nucleate boiling heat transfer is facilitated as a result of which the slugs of larger sizes (higher pressure reduction of 5.557 kPa at ΔTsub = 40 °C compared to that of 4.958 kPa at ΔTsub = 30 °C) form more frequently (increased evaporation rate) and more orderly. Consequently, the chaotic fluctuation of PHI signal is gradually reduced (refer the pressure signal at ΔTsub = 40 °C), which is the signature of enhanced nucleate boiling heat transfer. Interestingly, with further decrease in ΔTsub to 30 °C, the fluctuation of PHI signal becomes quasi-periodic with a dominant frequency of 0.002 Hz determined by FFT of the pressure signal. The quasi-periodic oscillations are of great importance as they match with the in-phase oscillations of wall temperature Twall at the exit of the heated section at ΔTsub = 30 °C and Q = 4.0 kW shown in Figure 12. At higher ΔTsub, Twall does not oscillate and remains almost constant with slightly increasing trend suggesting definitely the occurrence of nucleate boiling. With decrease in ΔTsub, pressure-drop oscillation initiates because of increase in compressible volume (due to increased vapor generation). The pressure-drop oscillations cause a premature critical heat flux (CHF) to occur. As a result, the threshold of Twall oscillations decreases considerably, and Figure 12 shows Twall oscillations at ΔTsub = 30 °C. Such Twall oscillations can be regarded as thermal oscillations and confirm that the boiling heat transfer regime has changed to the transition boiling at ΔTsub = 30 °C. It also indicates that the mode of flow instability of geysering type induced by wall 16476



Article

all ΔTsub. This clearly indicates that the transition boiling occurs at Q = 5.0 kW. 5.3. Autocorrelation Function. Another interesting test is carried out based on the autocorrelation function A(τ) for loop flow rate given by eq 4 to gain the knowledge on the chaotic nature of the time series. Figure 13 shows the plots of A(τ) vs τ.

5.4. Estimation of Correlation Dimensions. The analysis of our experimental time series presented in the preceding sections reveals the complexity of the chaos qualitatively at the parameter space of ΔTsub and Q. The analysis can not quantify the degree or level of complexity of the time series. More precisely, the complexity order of mathematical formulation to replicate the chaotic behavior of our NCBL cannot be determined. Estimation of the fractal dimensions, particularly D2, is therefore essential to gain the knowledge on the modeling of the complexity of chaos in our system. On the basis of the method described in section 2.3, we have determined D2. The effect of ΔTsub and Q on D2 has also been studied. Using the optimal value of τ = 6.0s (refer section 5.2), we varied m to obtain a converged value of D2. For a given value of m, D2 was determined by varying r up to 1.0−5. While determining D2, the major concern was the choice of the linear section of the plot of ln C(r) vs ln (r). According to the Grassberger−Procaccia algorithm, the primary criterion is to have the linear section for which the radius of the hypersphere is small enough, i.e., r → 0. We observed the following results in computing D2: 1. We observed that the range rmin ≤ r ≤ rmax, within which D2 is computed, is very small. The same is also observed by Ruelle.45 2. The lower part and the upper part are not useful. We worked on small middle part. Another observation is that if we increase τ, D2 does not converge well as m is increased. 3. For a given τ, the linear part becomes narrower as m is increased. Lai and Lerner46 also observed that for a small τ, a good scaling region is obtained and D2 is accurately obtained. The linear part can be expanded with m, provided τ is decreased. They

Figure 13. Autocorrelation functions varying with delay.

At all conditions of ΔTsub and Q, A(τ) decreases with τ indicating that the oscillations are chaotic in nature. It is also noted in the figure that the complexity reduces with Q as, the slopes are higher at Q = 4.0 kW than those at Q = 5.0 kW.

Figure 14. Correlation dimension at various ΔTsub and Q at τ = 6.0s. 16477



Article

microlayer−bubble, wall−bubble interface, bubbles themselves (coalescence and collision), and bubble−bulk liquid. To test the accuracy of D2, we varied the length of data series N. We computed D2 with the fist N = 1000 in addition to D2 for the total samples N = 1800. We achieve a good saturation of D2 for N = 1800, which are shown in Figure 14. D2 do not significantly differ from those for whole data set of N = 1800. The maximum deviation is 10.4% shown in Table 5. The deviations are computed with the reference D2 for N = 1800. Although the deviations are not significant, the negative deviations (D2 decreases with decrease in N) are more reasonable rather than the positive deviation.34 The positive deviation indicates that the time series of N = 1000 do not attain the sustained oscillations initially and contain some measurement noise. We therefore take D2 = 2.39 at ΔTsub= 30 °C and Q = 4.0 kW. The negative deviations are more at Q = 5.0 kW. The reductions of D2 for last N = 1000 at Q = 5.0 kW may be explained by the fact that boiling channel attains transition boiling through nucleate boiling which occurs initially. As mentioned, we experienced the transition boiling in the heated section at all conditions for Q = 5.0 kW. Computations of D2 reveal that the oscillations are more complex in nucleate boiling regime compared to those in the transition boiling regime. In an excellent review paper of Shoji,47 we found that the boiling regimes from nucleate to film boiling observed for a hot wire using a high-speed video camera are associated with highly complex nonlinear effects. Nonlinear analysis of the time series of surface temperature carried out by them suggests that the chaotic attractors in the nucleate boiling are more complex than those in the transition boiling regime. The power spectra in these regimes are broad-banded. Film boiling has least complexity. It is noted in Figure 14 that the minimum embedding dimensions remain in the range of 6−8 depending on the operating conditions. The computed values of D2 in Table 5 also verify the Takens’s proposal for minimum m ≥ 2D2 + 1 (the number in parentheses in Table 5 is minimum m after which further change in m changes the slope very little). We now compare our D2 with other related experimental or numerical results for boiling flows. Table 6 shows the comparison. D2 values for our experimental series are comparable with those of Heish et al.,10 Cammarata et al.,20 Narayanan et al.,31 Qu et al.,30 and Sathyamurthi et al.50 We note that the nuclear-coupled natural circulation loop has relatively less complex behavior than externally heated natural circulation loop. Interestingly, our D2 (3.09) at Q = 4.0 kW and ΔTsub = 50 °C are very close to that (3.0) estimated by Sathyamurthi et al.50 in their experiments on nucleate boiling. This is a very strong evidence that nucleate boiling occurs in our experiments at Q = 4.0 kW and ΔTsub = 50 °C. As mentioned earlier, D2 is used to quantify the complexity of the attractor by identifying the approximate number of independent variables (degrees of freedom) to describe both

recommended a lower bound and an upper bound of the linear section. The algorithm has been implemented in FORTRAN and verified with the help of the Lorenz problem. In our investigation, we computed D2 at ΔTsub = 30−50 °C and Q = 4.0 and 5.0 kW with 1800 samples. To validate our computation, we plotted D2 against m at all conditions. All the plots are shown in Figure 14. The figure shows a good saturation of D2. The saturated values of D2 at various conditions have been tabulated (see Table 5). The Table 5. Correlation Dimensions of Time Series of Wloop Q (kW)

ΔTsub (°C)

D2 total N = 1800

D2 last N = 1000

deviations wrt N = 1800 (%)

time delay, s

4 kW

50 40 30 50 40 30

3.09(7) 2.75(7) 2.29(6) 2.50 (7) 2.59(8) 2.17(8)

3.09 (7) 2.72(7) 2.39(6) 2.24(7) 2.35(8) 2.14 (8)

0.0 −1.0 +4.1 −10.4 −9.2 −1.3

6.0

5 kW

table shows D2 decreases with decrease in ΔTsub at a constant Q and decreases with Q at a constant ΔTsub. This confirms that the order of complexity of the chaotic oscillations of Wloop decreases with decrease in ΔTsub and with increase in Q. The phenomena have been explained in section 5.2. Decrease in D2 with Q in our experiments can also be verified with computations of D2 by Zboray et al.9 for Wloop time series generated in the experiments on their natural circulation test facility DESIRE. They determined a low-dimensional chaotic behavior in nuclearcoupled natural circulation loop and observed the decrease of D2 from 1.831 to 1.405 for the increase of power from 29.3 to 45.9 kW. In our experimental range, D2 has a maximum value of about 3.0 at ΔTsub = 50 °C and Q = 4.0 kW and a minimum value of about 2.0 at ΔTsub = 30 °C and Q = 5.0 kW. It indicates that the mathematical model to represent the dynamical system for the complexity of the chaos requires more independent variables (high degrees of freedom) at high ΔTsub and low Q at which nucleate boiling favors. The mathematical model requires less number of independent variables at low ΔTsub and high Q at which transition boiling occurs. The boiling heat transfer regime, which changes depending on inlet subcooling and heat flux, has therefore significant impact on the chaotic behavior of our NCBL. In transition boiling, as mentioned earlier, the bubbles coalesce to form large slug-bubbles that cover the entire diameter of the heated channel leaving the heated surface liquid-deficient. Complex microconvection of liquid around the bubbles is not present as much as in the nucleate boiling. In the nucleate boiling, the complex chaotic oscillation results from the complex interactions including those of the nucleation site (wall)− Table 6. Correlation Dimension in Some Relevant Experiments researchers

systems

experiments (E) /numerical model (M)

minimum m

D2

Heish et al.10 present study Cammarata et al.20 Narayanan et al.31 Qu et al.30 Liu et al.48 Lee and Pan49 Sathyamurthi et al.50

NCBL (water boiling) NCBL (water boiling) NCBL (1 − ϕ water) forced flow (water boiling) closed-loop pulsating heat pipe fluidized bed evaporator of natural circulation nuclear-coupled two-phase natural circulation loop nucleate boiling

E E E M E E M E

7 6−8

2.88 2.14−3.09 2.0−3.0 2.45 2.11−4.77 1.5−2.0 1.79 ± 0.01 3.0

16478

3 7 5



Article

the attractor and the dynamical NCBL. The complexity of the attractor increases with D2. As D2 is less than 3 in most of our test cases, use of two controllers (HPC and CLTC) can be justified, except in the case of very high ΔTsub and very low Q. We also test the accuracy of D2 which is carried out by reducing τ to 4.0 s. Although we obtained longer linear scaling section compared to that with τ = 6.0s, the saturation of D2 is not very convincing, as noted in Figure 15, even if m is increased to a high

boiling heat transfer regime has a significant effect on the complexity of the oscillations in the present NCBL. The effect of high subcooling and low heater power (favoring nucleate boiling) results in most complex geysering oscillations of thermodynamic origins while low inlet subcooling and high heater power (favoring transition boiling) induce hydrodynamic instability of reduced complexity such as PDO. We have been able to successfully identify and explain the geysering instability, PDO, and thermal oscillations in our experiments depending on inlet subcooling and heater power. Our observations have been supported by the experiments of Aritomi5 as well. 3. Uncertainty of D2 computations has been evaluated by changing τ and N. The effect of these parameters is present. A value of τ that is too low does not result in good saturation of D2. Optimal choice of delay time is highly required for a good estimation of D2 and is made in the present study. Varying N helps us identify the presence of noise in the time series and find the accurate value of D2. Low volume of data may result in some uncertainties in D2 estimation, as indicated by some test cases with Q = 5.0 kW. Although we have not computed the largest Lyapunov exponent, our study provides sufficient evidence for understanding the complex nonlinear behavior of our NCBL. Because of low sampling frequency, the study reports the knowledge of the overall dynamics of the loop rather than the bubble dynamics, and the former is strongly coupled with the later. Our future study will concentrate on the formulation of nonlinear models for the chaotic oscillations.

■

AUTHOR INFORMATION

Corresponding Author

Figure 15. Correlation dimension at various ΔTsub and Q at τ = 4.0 s.

*E-mail: [email protected]. Notes

The authors declare no competing financial interest.

value. D2 appears to saturate to a higher value when compared with D2 for τ = 6.0 s as noted in Figure 14. These evidence reveal that the choice of τ is very important for a good reconstruction of the attractor. Vazquez et al.27 used the method of delay to reconstruct attractors for pressure signals and acoustic signals in air−water bubbling flow. Using smaller τ, they obtained bad reconstruction of the attractors that does not reflect any physical implications.

■

ACKNOWLEDGMENTS The authors gratefully acknowledge Prof. S. Pushpavanam, IIT Madras, India, and Dr. Pradip Saha for the valuable discussion suggestions. They sincerely thank the peer reviewers of this paper for their valuable comments for improving the technical content of the paper. The help of N. Goswami during the collection of the experimental data is appreciated. The authors also thankfully acknowledge the financial support of Department of Science and Technology, New Delhi, India, under SERC scheme (Sanction no. SR/S3/CE/089/2009).

6. CONCLUSIONS We described nonlinear characteristics of the oscillating time series obtained in our experiments on NCBL in parameter space of inlet subcooling and heater power. Existing theories on bubble growth and bubble collapse have been found to successfully explain the mechanism of the chaotic oscillations. We summarize the conclusions as follows: 1. The effect of both thermodynamics and hydrodynamics has been discussed based on statistical analysis of the time series. The estimation of correlations coefficients reveals that the thermodynamic effect is more pronounced at high inlet subcooling and low heater power, while the hydrodynamic effect is induced and experienced to a great extent as the inlet subcooling is decreased or the heater power is increased. 2. The attractor reconstruction through time delay embedding and the subsequent computations of D2 reveal that

■

16479

NOMENCLATURE A = autocorrelation function a = cross-sectional area of channel (m2) C(r) = correlation sum cp = specific heat capacity (J kg−1 K−1) d = density D2 = correlation dimension D = bubble diameter (m) Δp = pressure drop (kPa) ΔT = temperature difference (°C) ΔTsub = inlet subcooling (°C) ΔTsat = wall superheating (°C) f = frequency (Hz) dx.doi.org/10.1021/ie300788t | Ind. Eng. Chem. Res. 2012, 51, 16467−16481


■

g = acceleration due to gravity H = Heaviside function k = discrete number L = length of channel (m) m = embedding dimension N = length of the time series n = number of points Nm = total number of points in the reconstructed phase space, Nm = N − (m − 1)τ P = pressure (kPa) Q = heater power (kW) r = radius of hypersphere T = temperature (oC) t = time (s) ugj = drift flux velocity (m/s) uin = inlet velocity (m/s) ui̅ n = average inlet velocity (m/s) Wloop = loop mass flow rate (kg/h) y = measured time series element Y = embedding vector

Article

REFERENCES

(1) Wissler, E.; Isbin, H. S.; Amudson, R. Oscillatory behavior of a twophase natural circulation loop. AIChE J. 1956, 2 (2), 157. (2) Ledinegg, M. Instability flow during natural forced circulation. Wärme 1938, 61 (8), 891−898. (3) Chexal, V. K.; Bergles, A. E. Two-phase Instabilities in a Low Pressure Natural Circulations Loop. AIChE Symp. Ser. 1973, 131, 37. (4) Fukuda, K.; Kobori, T. Classification of two-phase flow stability by density wave oscillation model. J. Nucl. Sci. Technol. 1979, 16, 95. (5) Aritomi, M.; Chiang, J. H.; Nakahashi, T. Fundamental study on thermo-hydraulics during start-up in natural circulation boiling water reactor (I), Thermo-hydraulic instabilities. J. Nucl. Sci. Technol. 1992, 29, 631. (6) Wu, C. Y.; Wang, S. B.; Pan, C. Chaotic Oscillations in a Low Pressure Two-Phase Natural Circulation Loop under Low Power and High Inlet Subcooling Conditions. Nucl. Eng. Des. 1996, 162, 223. (7) Kyung, I. S.; Lee, S. Y. Periodic of flow excursion in an open twophase natural circulation loop. Nucl. Eng. Des. 1996, 162, 233. (8) Jiang, S. Y.; Yao, M. H.; Bo, J. H.; Wu, S. R. Experimental simulation study on startup of the 5 MW nuclear heating reactor. Nucl. Eng. Des. 1995, 158, 111. (9) Zboray, R.; De Kruijf, W. J. M.; Van Der Hagen, T.H. J. J. RizwanUddin. Investigating the nonlinear dynamics of natural-circulation, boiling two-phase flows. Nucl. Technol. 2004, 146, 244. (10) Hsieh, C. C.; Wang, S. B.; Pan, C. Dynamic visualization of twophase flow patterns in a natural circulation loop. Int. J. Multiphase Flow 1997, 23, 1147. (11) Jain, V.; Nayak, A. K.; Vijayan, P. K.; Saha, D.; Sinha, R. K. Experimental investigation on the flow instability behavior of a multichannel boiling natural circulation loop at low-pressures. Exp. Therm. Fluid Sci. 2010, 34, 776. (12) Paruya, S.; Saha, A. K.; Bhattacharya, P. Validations of thermohydraulic models for geysering in a natural circulation loop using an impedance needle probe. Ind. Eng. Chem. Res. 2009, 48, 2020. (13) Boure, J. A.; Bergles, A. E.; Tong, L. S. Review of Two-Phase Flow Instability. Nucl. Eng. Des. 1973, 25, 165. (14) Griffith, P. Geysering in Liquid-Filled Lines; ASME: New York, 1963; paper no. 62-HT-39. (15) Furuya, M.; Inada, F.; Yasuo, A. Density Wave Oscillations of a Boiling Natural Circulation Loop Induced By Flashing. Presented at the 7th International Topical Meeting on Reactor Thermal Hydraulics (NURETH-7), Saratoga Springs, NY, Sept. 10−15, 1995, ANS, 923. (16) Saha, P.; Ishii, M.; Zuber, N. An experimental investigation of the thermally induced two-phase flow oscillations two-phase systems. Trans. ASME, J. Heat Transfer 1976, 98, 616. (17) Parker, T. S.; Chua, L. O. Practical Numerical Algorithms for Chaotic Systems; Springer−Verlag: New York, 1989. (18) Takens, F. Detecting strange attractors in turbulence. Dynamical Systems and Turbulence, Lecture Notes in Mathematics; Rand, D. A., Young, L. S., Eds.; Springer−Verlag: Berlin, 1981; p 366. (19) Schuster, H. G. Deterministic Chaos, An Introduction (in Polish); PWN: Warsaw, Poland, 1993. (20) Cammarata, G.; Fichera, A.; Pagano, A. Nonlinear Analysis of a Rectangular Natural Circulation Loop. Int. Commun. Heat Mass Transfer 2000, 27 (8), 1077. (21) Fraser, A. M.; Swinney, H. L. Independent coordinates for strange attractors from mutual information. Phys. Rev. A 1986, 33, 1134. (22) Grassberger, P.; Procaccia, I. Measuring the strangeness of strange attractors. Phys. D 1983, 9, 189. (23) Wolf, A.; Swift, J. B.; Swinney, H. L.; Vastano, J. A. Determining Lyapunov exponents from a time series. Phys. D 1985, 16, 285. (24) Rosenstein, M. T.; Collins, J. J.; De Luca, C. J. A practical method for calculating largest Lyapunov exponents from small data sets. Phys. D 1993, 65 (1, 2), 117. (25) Hegger, R.; Kantza, H.; Schreiber, T. Practical implementation of nonlinear time series methods: The TISEAN package. Chaos 1999, 9 (2), 413. (26) Mosdorf, R.; Shoji, M. Chaos in bubblingnonlinear analysis and modeling. Chem. Eng. Sci. 2003, 58, 3837.

Greek Symbols

λ = Lyapunov exponent δ = δi = yi − ymean and ymean is the mean value of the time series ρ = power spectral density ỹ = piecewise stationary signal ω = frequency spectrum τ = time delay (s) σn = standard deviation of n points Subscripts

bd = bubble departure bod = boiling delay bt = bubble transit g = bubble growth H = heater HI = heater inlet l = liquid mean = mean value R = riser s = sampling sat = saturation sub = inlet subcooling w = bubble waiting wall = heater wall Abbreviations

AMI = average mutual information BHTR = Boiling heat transfer regime CLTC = cold-leg temperature controller DPT = differential pressure transmitter DWO = density-wave oscillation EFM = electromagnetic flow meter FFT = fast Fourier transform HPC = heater power controller NCBL = natural circulation boiling loop PSD = power spectrum density PDO = pressure drop oscillation RTD = resistance temperature detector SD = standard deviation SVD = singular value decomposition TDPSR = time delay phase space reconstruction 16480



Article

(27) Vazquez, A. Manasseh, R.; Sánchez, R.M.; Metcalfe, G. Experimental comparison between acoustic and pressure signals from a bubbling flow. Chem. Eng. Sci. 2008, 63, 5860. (28) Lin, T.-J.; Juang, R.-C.; Chen, C.-C. Characterizations of flow regime transitions in a high-pressure bubble column by chaotic time series analysis of pressure fluctuation signals. Chem. Eng. Sci. 2001, 56, 6241. (29) Lakshmanan, S. P.; Pandey., M. Analysis of startup oscillations in natural circulation boiling systems. Nucl. Eng. Des. 2009, 239, 2391. (30) Qu, J.; Wu, H.; Cheng, P.; Wang, X. Non-linear analyses of temperature oscillations in a closed-loop pulsating heat pipe. Int. J. Heat Mass Transfer 2009, 52, 3481. (31) Narayanan, S.; Srinivas, B.; Pushpavanam, S.; Murty Bhallamudi, S. Non-linear dynamics of a two phase flow system in an evaporator: The effects of (i) a time varying pressure drop (ii) an axially varying heat flux. Nucl. Eng. Des. 1997, 178, 279. (32) Leibert, W.; Schuster, H. G. Proper choice of delay time for the analysis of chaotic time series. Phys. Lett. A 1991, 142, 107. (33) Kennel, M. B.; Brown, R.; Abarbanel, H. D. I. Determining embedding dimension for phase-space reconstruction using a geometrical construction. Phys. Rev. A 1992, 45, 3403. (34) Kantz, H.; Schreiber, T. Nonlinear Time Series Analysis; Cambridge University Press: Cambridge, U.K., 1997. (35) Eckmann, J.-P.; Ruelle, D. Ergodic theory of chaos and strange attractors. Rev. Mod. Phys. 1985, 57, 617. (36) Dhir, V. K. Numerical Simulations of Pool-Boiling Heat Transfer. AIChE J. 2001, 47, 813. (37) Basu, N.; Warrier, G. R.; Dhir, V. K. Wall Heat Flux Partitioning During Subcooled Flow Boiling: Part 1 & 2Model Development. J. Heat Transfer 2005, 127, 131. (38) Kamil, M.; Shamsuzzoha, M.; Abdul Hakeem, M. Analysis of incipience of nucleate boiling in a reboiler tube. AIChE J. 2007, 53 (1), 39. (39) Karmakar, A.; Goswami, N.; Paruya, S.; Subcooled Boiling Oscillations in Natural Circulation Boiling Loop at Low Pressure. AIChE J., submitted for publication. (40) Mosdorf, R.; Shoji, M. Chaos in nucleate boiling-non-linear analysis and modeling. J. Heat Mass Transfer 2004, 47, 1515. (41) Sathyamurthi, V.; Banerjee, D. Non-linear dynamical analyses of transient surface temperature fluctuations during subcooled pool boiling on a horizontal disk. Int. J. Heat Mass Transfer 2009, 52, 5608. (42) Gaertner, R. F. Photographic Study of Nucleate Pool Boiling on a Horizontal Surface. ASME J. Heat Transfer 1965, 8, 17. (43) Liu, H. T.; Kakaç, S.; Mayinger, F. Characteristics of transition boiling and thermal oscillation in an upflow convective boiling system. Exp. Therm. Fluid Sci. 1994, 8, 195. (44) Liu, H. T.; Kakaç, S. An experimental investigation of thermally induced flow instabilities in a convective boiling upflow system. Waerme−Stoffuebertrag. 1991, 26, 365. (45) Ruelle, D. Deterministic chaos: the science and the fiction. The Claude Bernard Lecture, 1989. Proc. R. Soc. London A 1990, 427, 241. (46) Lai, Y.-C.; Lerner, D. Effective scaling regime for computing the correlation dimension from chaotic time series. Phys. D 1998, 115, 1. (47) Shoji, M. Studies of boiling chaos: a review. Int. J. Heat Mass Transfer 2004, 47, 1105. (48) Liu, M-Y; Qiang, A-H; Sun, B.-F. Chaotic characteristics in an evaporator with a vapor−liquid−solid boiling flow. Chem. Eng. Process.− Process Intens. 2006, 45 (1), 73. (49) Lee, J. D.; Pan, C. Nonlinear analysis for a nuclear-coupled twophase natural circulation loop. Nucl. Eng. Des. 2005, 235 (5), 613. (50) Sathyamurthi, V.; Banerjee, D.; Sakamoto, H.; Kim, J. Measurement of the fractal order of wall void fraction during nucleate boiling. Int. J. Heat Fluid Flow 2008, 29, 207.

16481


Nonlinear Analysis of Chaotic Time Series in a Natural Circulation

Recommend Documents