SCWR single channel stability analysis using a response matrix method

Nuclear Engineering and Design 241 (2011) 2528–2535

Contents lists available at ScienceDirect

Nuclear Engineering and Design journal homepage: www.elsevier.com/locate/nucengdes

SCWR single channel stability analysis using a response matrix method Jiyun Zhao ∗ , C.P. Tso, K.J. Tseng Division of Power Engineering, Nanyang Technological University, Block S2, Level B2c, Room 126, 50 Nanyang Avenue, Singapore 639798, Singapore

a r t i c l e

i n f o

Article history: Received 7 January 2011 Received in revised form 5 April 2011 Accepted 14 April 2011

a b s t r a c t A system response matrix method, which directly solves the linearized differential equations in the matrix form without Laplace transformation, is introduced for the supercritical fluids flow instability analysis. The model is developed and applied to the single channel or parallel channel type instability analyses of a typical proposed Supercritical Water Reactor (SCWR) design. A uniform axial heat flux is assumed, and the dynamics of the fuel rods and water rods are not considered in this paper. The sensitivity of the decay ratio (DR) to the axial mesh size is analyzed and found that the DR is not sensitive to mesh size once sufficient number of axial nodes is applied. The sensitivity of the stability to inlet orifice coefficient is conducted for the hot channel and found that a higher inlet orifice coefficient will make the system more stable. The susceptibility of stability to operating parameters such as mass flow rate, power and system pressure is also performed. It is found that the SCWR stability sensitivity feature can be improved by carefully choosing the inlet orifice coefficients and operating parameters. The stability feature of the average channel is also analyzed with an equivalent inlet orifice coefficient. Finally, the manufacturing feasibility of the inlet orifices for both the hot channel and average channel is studied and found to be favorable. © 2011 Elsevier B.V. All rights reserved.

1. Introduction Due to the drastic water density change through the core, the Density wave oscillation (DWO) instability is a well known issue during the development of the supercritical water cooled reactors (SCWRs). Flow stability in a nuclear reactor system during DWO can be studied using either a non-linear model (non-linear stability analysis) or a linear model (linear stability analysis) of the onedimensional flow governing equations. In the non-linear stability analysis, two kinds of methods are widely used. One method is using the transient time domain codes such as RAMONA-5B and RELAP5 that performs numerical integration of the non-linear differential equations in the time domain. However, using system codes requires large computational efforts and substantial work to handle the numerical uncertainties. The other method is using reduced-order analytical models such as those developed in Dokhane et al. (2003a,b). For the linear stability analysis, the system models are simplified through linearization of the complex non-linear differential equations. These linearized equations can then be solved through Laplace transform in the frequency domain or directly solved in system response matrix form. However, the system response matrix method needs to collect all of the system information into one

∗ Corresponding author. Tel.: +65 67904508; fax: +65 67933318. E-mail address: [email protected] (J. Zhao). 0029-5493/$ – see front matter © 2011 Elsevier B.V. All rights reserved. doi:10.1016/j.nucengdes.2011.04.026

matrix, which makes the matrix very complicated and necessitates matrix computation techniques such as sparse matrix and matrix partition (Hanggi, 2001). Although some transient information is lost through model linearization, the high computational efficiency and relatively accurate results make the method attractive, especially for prediction of the onset of instability. In fact, boiling water reactors (BWRs), have been experimentally found to behave as linear systems under normal operating conditions (Carmichael and Niemi, 1978). Obviously, the transient codes are expected to represent the system more accurately without major simplifications, but the computational time consumption and numerical stability difficulties often plague the application of this method. Similarly SCWRs could also have three density wave oscillation type instabilities, which are single channel instability, core wide coupled neutronic inphase instability and core wide coupled neutronic out-of-phase instability. During the single channel type oscillations, only one channel or a small fraction of the parallel channels oscillates, while the other channels remain at steady state. This imposes a constant pressure drop boundary condition across the oscillating channel or channels. This type of DWO is also called parallel channel type DWO (Podowski, 2003). Since only a small fraction of the core flow oscillates while the bulk flow remains at steady state during the single channel oscillation, the neutronic feedback due to this small fraction oscillation can be neglected, decoupling the thermalhydraulics dynamics from the neutronic dynamics. This type of

J. Zhao et al. / Nuclear Engineering and Design 241 (2011) 2528–2535

2529

Table 1 Typical SCWR coolant conditions and reactor power. Fluid Operating pressure (MPa) Inlet/outlet temperature (◦ C) Reactor flow rate (kg/s) Total water rods flow (%) Thermal power (MWt)

Water 25 280/500 1843 90 3575

instability could be very detrimental and must be avoided during the design stage. During the core wide in-phase oscillation, the whole core oscillates in the same phase and the fundamental mode of neutronic dynamics is excited. During the core wide out-of-phase instability, half of the core rises in power while the power in the other half decrease to maintain an approximately constant total core power. The first subcritical mode of neutronic dynamics is excited during the out-of-phase oscillations. The single channel stability for the U.S. reference SCWR design was analyzed with linear model in frequency domain by Yang and Zavaljevski (2003), and the Japanese SCWR design was analyzed by Koshizuka et al. (2003) and Yi et al. (2004). By analogy, the analyses developed for the BWR stability, Zhao et al. (2007a,b) developed a three-region model for the supercritical water flow and derived new non-dimensional parameters that govern the supercritical water flow instabilities during the study of the SCWR stabilities. The linear method was also used by Ortega-Gómez et al. (2008) and Gallaway and Podowski (2010) for the supercritical fluids single channel stability analyses. A non-linear method was used to perform supercritical flow stability analysis in horizontal channels (Chatoorgoon, 2008). In this paper, using a linear model, the response matrix method will be introduced and applied to a typical SCWR design single channel type stability analysis, for comparison with the Laplace transform method. 2. Description of a typical SCWR design The proposed typical SCWR design (Buongiorno and MacDonald, 2003) is a thermal neutron spectrum using supercritical water as the coolant. The system pressure is 25 MPa with an inlet water temperature of 280 ◦ C and core average outlet temperature of 500 ◦ C. Due to the small coolant density, especially in the upper part of the core, water rods are required to provide additional moderation. The supercritical water flows into the reactor vessel through the inlet nozzle; then splits and flows partly through the downcomer and partly through the water rods. After mixing in the lower plenum, water flows upward through the core channels to remove the fission energy. For the upward core flow, inlet orificing is expected to be used for improving the power to flow ratio and for stabilizing the SCWR system. The inlet orifice coefficient for the SCWR assemblies must be carefully selected to direct more flow to the higher power assemblies so as to minimize the differences in coolant enthalpy increase among different assemblies. 2.1. Core and fuel assembly descriptions The coolant conditions and core power for the proposed typical SCWR design are summarized in Table 1, as follows. Just like a BWR assembly, the SCWR assembly is also contained in a channel box to avoid coolant mixing among different assemblies. The cross section of a typical SCWR fuel assembly can be seen in Fig. 1 (Buongiorno and MacDonald, 2003). The proposed typical design of the SCWR includes 145 fuel assemblies. The core average power density is 70 kW/L with a target power peaking factor of 2.0. The core pressure drop target is

Fig. 1. SCWR fuel assembly.

Table 2 Typical SCWR core and assembly parameters. Core Number of fuel assembly Axial/radial/local/total peaking factor Average linear power (kW/m) Target core pressure drop (MPa) Fuel assembly Number of fuel pins per assembly Number of water rods per assembly Water rod side (outside dimension) (mm) Assembly duct thickness (mm) Assembly side (outside dimension) (mm) Assembly hydraulic diameter (mm)

145 1.4/1.3/1.1/2.0 (best estimates) 19.2 0.15 300 36 33.6 3 286 3.4

0.15 MPa which is comparable to the traditional PWR core pressure drop. The basic design parameters of the SCWR core and assembly are listed in Table 2. Due to introduction of water rods, the SCWR assembly hydraulic diameter is small, only about 3.4 mm. Although the fuel pin diameter is comparable to that of the traditional LWRs, the fuel length of a typical SCWR design is somewhat longer than that of the traditional LWRs. The basic design parameters of the fuel pin can be found in Table 3. 2.2. Hot channel description From Table 2, a radial power peaking factor of 1.3 has been chosen for the hot channel with an assumption of uniform axial heat flux. Since inlet orifices can be applied for the SCWR and the orifice coefficient can be adjusted to deliver more flow to the hot channel, the hot channel and average channel can have the same enthalpy rise. Table 3 Typical SCWR fuel pin parameters. Fuel pin OD (mm) Fuel pin pitch (mm) Heated length (m) Fission gas plenum length (m) Total fuel pin height (m)

10.2 11.2 4.27 0.6 4.87

2530

J. Zhao et al. / Nuclear Engineering and Design 241 (2011) 2528–2535

It should be noted that, in reality, the same enthalpy rise in the hot and average channels may not be achieved exactly because of two reasons. Firstly, the radial power profile changes with time and, secondly, due to design constraints, only zone-by-zone orificing can be provided (i.e. grouping together assemblies with similar power). However, orifice optimization can be sought to achieve a condition of similar enthalpy rise in all channels.

work shown in Eq. (4), a generalized system response matrix will be constructed and the generalized eigenvalue problem to be solved as shown below. The discussions of the generalized eigenvalue problem can also be found in Golub et al. (1996).

3. Response matrix method and its application to SCWR single channel stability analysis

where Ageneral is the generalized system response matrix. The generalized eigenvalue problem can be described as follows:



ıx 0



 = Ageneral

ıx ıy

 ,

(5)

Ageneral ei = i Bgeneral ei ,

(6)

where matrix Bgeneral has the following form:

0

… 1

…

0

…

0⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ 0⎦

Referring to the state variables

…

⎡1 ⎢0 ⎢ ⎢ Bgeneral = ⎢ ⎢ ⎢ ⎢ ⎣0

…

As shown in Fig. 1, the SCWR fuel assembly is contained in a channel box similar to a BWR assembly to avoid coolant mixing between different assemblies. During the parallel or single channel flow instability, the flow rate fluctuation of the unstable channel will not affect the flow of the remaining channels because the single channel flow is only a small fraction of the whole core flow. Therefore, a constant pressure drop boundary condition can be imposed on that single channel. Due to the small fraction of the single channel flow compared to that of the whole reactor core, the neutronics feedback due to the flow fluctuation of a single channel will not affect the whole core neutronic properties much, and the neutronic is uncoupled (March-Leuba and Rey, 1993). Hence, only the thermal-hydraulics model is implemented here.

d dt

Referring to the constitutive variables

3.1. Analysis methodology For investigating the dynamic features of a thermal-hydraulic system, the three governing equations are the mass conservation equation, the momentum conservation equation and the energy conservation equation. These equations along with the constitutive equations will close the system. Accordingly, the variables can be categorized to be state variables which are described by differential equations and constitutive variables which are described by algebraic equations. Therefore, the set of system equations can be represented as follows: d x(t) = f (x(t), y(t)) dt

and

0 = g(x(t), y(t)),

(1)

where x(t) is the vector of state variables and y(t) is the vector of constitutive variables. Perturbation and linearization of Eq. (1) yields the following first order ordinary differential equations. d ∂f ∂f ıx = ıx + ıy = Aıx + Bıy, dt ∂x ∂y

(2)

∂g ∂g ıx + ıy = Cıx + Dıy, ∂x ∂y

(3)

0=

where A, B, C, and D are coefficient matrices corresponding to respective variable perturbation vectors. If we solve the constitutive variables from Eq. (3) and substitute them into Eq. (2), we obtain a new differential equation set, which includes only the state variables, as: d ıx = (A − BD−1 C)ıx = As ıx. dt

i : ith eigenvalue of system response matrix Ageneral and ei : eigenvector corresponding to ith eigenvalue i . From the dynamic theory, if all of the eigenvalues of a system response matrix Ageneral have negative real parts, the disturbance will asymptotically decay away and the system is stable in that case. After sufficient time, the system oscillation can be determined by the eigenvalue which has the largest real part, called the dominant eigenvalue. For calculating the dominant eigenvalue, an inverse iteration combined with Newton’s method was applied in this work as discussed below. The discussions of this method can also be found in Peters and Wilkinson (1979) and Ipsen (1997). The iteration method consists of the following algorithm:

(4)

The matrix As , in the above Eq. (4) is called the system response matrix. The discussion of the system response matrix can also be found in Hanggi (2001). The system equations considered in the reactor stability analysis are very complicated. To avoid the tedious algebraic derivation

1.  Start withguess  of e0 and 0 2.

ek+1 − ek

k+1 − k

=−

uT1

0

Ageneral − k Bgeneral −Bgeneral ek

 \



uT1 ek − 1 (Ageneral − k Bgeneral )ek

,

uT1

where = [1, 0, . . . , 0] is of the same size as ek and the back slash “\” is a mathematic function in MATLAB, where“A\b” means inverse “A” times “b” or “A−1 b”. 3. ek+1 = ek + e, 4. k+1 = k + , 5. If ||Ageneral ek+1 − k+1 ek+1 || > tolerance go to step 2. Using the above algorithm, the dominant eigenvalue can be calculated very quickly. The dominant eigenvalue has a formula of 1 =  + i , where  is the real part of the dominant eigenvalue, and  the imaginary part. The system response to the dominant eigenvalue is illustrated in Fig. 2. Therefore, the decay ratio, DR, can be defined as the ratio of the amplitudes at time t1 and t2 , i.e. at sequential oscillating periods. DR =



|||| u(t2 ) = exp −2 |||| u(t1 )

 .

(7)

In the following section, this methodology is applied to calculate the DR for a single channel of SCWR at supercritical pressure.

J. Zhao et al. / Nuclear Engineering and Design 241 (2011) 2528–2535

2531

u(t )

t1

t2

t

Fig. 2. Inlet velocity response corresponding to the dominant eigenvalue.

3.2. Decay ratio calculation for a typical SCWR The governing equations for SCWR at steady state for a single channel can be written as:

(a) Mass conservation equation ∂(u) ∂ + = 0. ∂t ∂z

Fig. 3. SCWR single channel nodalization and perturbation.

(8)

(b) Energy conservation equation ∂h ∂(uh) ∂p ∂p P + = +u + h q (z, t). Ac ∂t ∂z ∂t ∂z

(9)

Applying the above constitutive equation to node i and ignoring the fluctuation of surface heat flux, we obtain the following three differential equations: u (h − hi ) i−1 + i−1 ui−1 dıhi = i−1 i−1 ıhi−1 dt zi

(c) Momentum conservation equation ∂(u) ∂(u2 ) ∂p f u2 + =− − g − . De 2 ∂t ∂z ∂z

ui−1 (hi−1 − hi ) i−1  (h − hi ) ıpi−1 + i−1 i−1 ıui−1 zi zi

−

i−1 ui−1 ıhi , zi

(10)

−

(11)

ui−1 (ui − ui−1 ) i−1 ıhi−1 − zi

 After linearization and perturbation of the above four equations, one can obtain three first order differential equations corresponding to three state variables. The state variables were chosen as: pressure P, enthalpy h and coolant velocity u. If the system is divided into n axial nodes, the fluctuation of coolant density at node i can be obtained from the state equation: ıi = i ıhi + i ıpi ,

(12)

where i and i are the derivatives of the coolant density with respect to enthalpy and pressure at node i. The derivatives are calculated using NIST Standard Reference Database (Harvey et al., 2008).



i =

 i =



∂ ∂p

h

i . ∂ ∂h p

i

(14)

1 − ui−1 (ui − ui−1 ) i−1 dıui = ıpi−1 − dt zi

(d) Coolant state (constitutive) equation  = (p, h).

+

(13)

−

dıpi = dt



−



1 + zi



fi u2i g + i 2De i

fi u2i g + i 2De i



i ıhi −

 

i ıpi

i−1 (ui − 2ui−1 ) ıui−1 zi

(15)

ıui ,



i−1 i−1 (hi−1 − hi ) i−1 ui−1 − i ui i − i ıui−1 − ıui + i ıhi z i zi i z i zi i ui ıpi + z

 +

i−1 ui−1 fi ui + De zi









ui−1 i−1 ui−1 (hi−1 − hi ) i−1 + i−1 ui−1 − i ıhi−1 z i zi i



ui−1 i−1 ui−1 (hi−1 − hi ) − i i−1 ıpi−1 . z i zi i

(16)

The single channel nodalization and perturbations of the state variables are illustrated in Fig. 3. It can be seen from Fig. 3 that there are (N + 3) nodes in the single channel consisting of N heated nodes, one inlet orifice node, one lower gas plenum (non-heated) node and one upper gas plenum (non-heated) node. Therefore, there are in total 3(N + 4) state variables, since the variables are defined at node boundaries. Applying Eqs. (14)–(16) at every node, 3(N + 3) equations can be obtained. Thus, 3 additional equations or boundary conditions are required to close the system. The boundary conditions used in the present analysis are as follows:

2532

J. Zhao et al. / Nuclear Engineering and Design 241 (2011) 2528–2535

Axial nodes number Axial mesh size (m) Decay ratio

40 0.1068 0.1066

60 0.0712 0.1107

80 0.0534 0.1132

100 0.0427 0.1143

(a) Ignore the fluctuation of inlet enthalpy, i.e. ıhin = 0. (b) If the inlet and outlet pressures are assumed constant, i.e. ıhin = 0 and ıhexit = 0, this automatically satisfies the constant pressure drop boundary condition across the flow channel. Writing the above mentioned 3(N + 4) equations in matrix form described as Eq. (5), one can obtain the generalized system matrix Ageneral . The geometry of the assembly flow channel together with the coolant state and flow conditions for every axial node can be calculated using the data provided in Tables 1–3. Then, using the decay ratio calculation methodology described in Section 3.1, we calculate the decay ratios for single hot and average channels. 3.3. Results for the typical SCWR design As mentioned earlier, for SCWR, an inlet orifice scheme is expected to be used to adjust the flow distribution among fuel assemblies to ensure that the hot and the average channels have almost the same coolant properties at any height within the core. From Table 2, it can be seen that the radial power peaking factor of the hot channel is 1.3; thus, to obtain the same enthalpy rise in the hot and the average channels, the flow rate of the hot channel should also be 1.3 times that of the average channel. A uniform axial heat flux profile was assumed. The usual BWR thermal-hydraulic stability criterion of decay ratio less than 0.5 was assumed for both the hot and the average channel analyses. 3.3.1. The axial mesh size effects on decay ratio Using a traditional Laplace transformation method, Yang and Zavaljevski (2003) found that the decay ratio is significantly dependent on the axial mesh size. To obtain the decay ratios at zero mesh size, the decay ratios at several axial mesh sizes were calculated and extrapolated to zero mesh size. To determine the axial mesh size effects on decay ratio for the response matrix method, the decay ratios at different mesh sizes were calculated for the hot channel at inlet orifice coefficient 20. The results are presented in Table 4, and graphically shown in Fig. 4. From Table 4 and Fig. 4, it is seen that the decay ratio is not significantly dependent on mesh size, especially at high node number, say above 40. Therefore, eighty (80) nodes were applied axially for the following analysis. 3.3.2. Minimum inlet orifice coefficient for the hot channel The inlet orifice coefficient has a very important role in the system stability characteristics. The decay ratios at different inlet orifice coefficient values for the hot channel are plotted in Fig. 5.

0.15 0.10 0.05 0.00 40

60 80 Axial mesh number

Fig. 4. Decay ratio dependence on axial mesh number.

100

0.6 0.5 0.4 0.3 0.2 0.1 0 0

5

10

15

20

25

Inlet orifice coefficient at Hot Channel Fig. 5. SCWR hot channel decay ratios at different inlet orifice coefficients.

From Fig. 5, it is seen that the decay ratio will be below 0.5 if the inlet orifice coefficient for the hot channel is greater than about 1.0. However, to have a more stable system and to account for uncertainties and off-normal conditions, the designer should select a larger value; say 10–20, for the hot channel inlet orifice coefficient. The choice of the inlet orifice coefficient has dual influences. A higher orifice coefficient will stabilize the system. On the other hand, a higher orifice coefficient will produce a higher core pressure drop, and that means higher coolant pumping power, which, in turn, is harmful from the economics viewpoint. However, an inlet orifice coefficient of 20.0 for the hot channel would produce a core pressure drop of 0.163 MPa which is quite close to the target core pressure drop of 0.15 MPa, and thus should be acceptable. With the inlet orifice coefficient of 20.0, the decay ratio for the hot channel is only 0.11 as seen from Fig. 5. Thus, the hot channel will be very stable, and in all the following calculations, the inlet orifice coefficient of the hot channel has been assumed to be 20.0. For comparison, a typical BWR/4 has an inlet orifice coefficient of 27.8 at the hot channel (Lahey and Moody, 1993). 3.3.3. Sensitivity analysis of hot channel stability with inlet orifice coefficient of 20.0 The decay ratio sensitivity analysis for three key system parameters is addressed in this part: mass flow rate, power and system pressure. Full 100% nominal flow rate and full 100% nominal power at a system pressure of 25 MPa are taken as the reference conditions. 3.3.3.1. Mass flow rate sensitivity analysis. The decay ratios have been calculated for reduced mass flow rate, while keeping the power and the pressure at reference conditions. Therefore, the power to flow ratio is increased as reducing the flow rate while maintaining the power. The results are illustrated in Fig. 6.

Decay ratio

Decay ratio

0.20

Decay Ratio

Table 4 Mesh size effects on decay ratio for the hot channel (inlet orifice coefficient at 20).

0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 0.8 0.85 0.9 0.95 Fraction of reference mass flow rate Fig. 6. SCWR mass flow rate sensitivity analysis.

1

0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 100

110

120

2533

0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0

Decay ratio

Decay ratio

J. Zhao et al. / Nuclear Engineering and Design 241 (2011) 2528–2535

130

Kin=40 Kin=20

0.80 0.85 0.90 0.95 1.00

Percentage of reference power (%)

Fraction of reference mass flow rate

Fig. 7. SCWR power sensitivity analysis.

Fig. 9. Mass flow rate sensitivity at different inlet orifice coefficients.

From Fig. 6, it is seen that the decay ratio will be above 0.5 at flow rates below about 83% of the reference case. This means that a reduction in mass flow rate, while keeping the power constant, has a destabilizing effect.

design is more sensitive to the operating parameters than BWR. In Section 3.3.3 of this paper, the SCWR stability sensitivity feature was investigated with an inlet orifice coefficient of 20 for the hot channel. To find a measure in order to improve the sensitivity feature of the SCWR, the sensitivities are compared for the different inlet orifice coefficients in Figs. 9 and 10 as follows: From Figs. 9 and 10 above, it can be seen that a higher inlet orifice coefficient can improve not only the stability but also the stability sensitivity to operating parameters. Therefore, the SCWR can be designed stable with an improved stability sensitivity feature by carefully choosing the inlet orifice coefficients.

3.3.3.2. Power sensitivity analysis. The decay ratios are calculated by changing the power level, with constant mass flow rate and pressure equal to the reference condition. Therefore, the power to flow ratio is increased as increasing the power level while maintaining the flow rate. The results are illustrated in Fig. 7. From Fig. 7, it is seen that the decay ratio will be above 0.5 as the power rises above 121% of the reference case. This means that an increase in power, while keeping the mass flow rate constant, has a destabilizing effect. 3.3.3.3. Pressure sensitivity analysis. The decay ratios have been calculated for a range of the system pressure, while keeping the mass flow rate and power at the reference condition. The results are illustrated in Fig. 8. From Fig. 8, it is seen that the decay ratio will decrease as the pressure increases, which means that a higher system pressure will stabilize the system. Conversely, a lower system pressure will destabilize the system, although the decay ratio increases only slightly when the system pressure is lowered from 25 MPa to 23 MPa, thus indicating that the effect of system pressure is less important compared to that of flow rate and/or power. These results are similar to what were obtained using the Laplace transformation method and reported in (Zhao et al., 2004). The slight difference in the two sets of results originates from the different reference SCWR designs since the SCWR design is still undergoing revision.

0.16 0.15 0.14 0.13 0.12 0.11 0.10 0.09 0.08

(a) The same power-to-flow rate ratio as the hot channel. (b) The core pressure drop of 0.163 MPa which yields the hot channel inlet orifice coefficient of 20.0. With an inlet orifice coefficient of 115.0, the decay ratio for the average channel is found to be only 0.01. With such a small decay ratio, the average channel will be very stable. 3.3.6. Manufacturing feasibility study for the inlet orifice To check if the large inlet orifice coefficient such as 115.0 is feasible from manufacturing point of view, the following analysis was conducted. The inlet orifice is illustrated in Fig. 11 as follows: For incompressible flow, from the Bernoulli equation, p1 +

u21 2

23

24 25 26 System pressure (MPa)

Fig. 8. SCWR system pressure sensitivity analysis.

27

+ gz1 = p2 +

Then, 1 p1 − p2 = u21 2

Decay ratio

Decay ratio

3.3.4. Measure to improve the SCWR sensitivity feature The sensitivity of the SCWR stability versus operating parameters was compared with that of a typical boiling water reactor (BWR) in Zhao et al. (2007b), it was found that U.S. reference SCWR

3.3.5. Average channel analysis The inlet orifice coefficient for the average channel is calculated to be 115.0 to satisfy both the requirements of:

u22

 4 D 1

D2

2

+ gz2 .

(17)

 −1 ,

(18)

0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0

Kin=40

Kin=20

100

110

120

130

Percentage of reference power (%) Fig. 10. Power sensitivity at different inlet orifice coefficients.

2534

J. Zhao et al. / Nuclear Engineering and Design 241 (2011) 2528–2535

u1 ,P1

D1

D2 u2 ,P2

Fig. 11. Orifice illustration.

where subscript 1 designates the “upstream” or “pipe” condition and subscript 2 designates the “downstream” or “orifice” condition. From Table 7.1 of Holman (2001), it was found that the permanent pressure loss factor of a square-edged orifice is 0.8. Then, ploss = 0.8(p1 − p2 ) = Korifice Therefore, Korifice = 0.8

 4 D 1

D2

u21 2

.

(19)

 −1 .

(20)

For Korifice = 20 (as suggested for the hot channel), D1 /D2 = 2.26. For Korifice = 115 (as suggested for the average channel), D1 /D2 = 3.47. Thus, both of the above diameter ratios (Dpipe /Dorifice ) are quite reasonable. For the typical SCWR assembly, the equivalent inlet pipe diameter would be about 11.8 cm. Therefore, for the hot channel, the orifice diameter would be ∼5.2 cm corresponding to an inlet orifice coefficient of 20. For the average channel, the orifice diameter would be ∼3.4 cm corresponding to the inlet orifice coefficient of 115. Therefore, there should be no difficulty in manufacturing this orifice. From the above calculations, it is seen that the typical SCWR design should have no density-wave instability problem at a core pressure drop of 0.163 MPa with proper inlet orifices for both the hot and average channels for full-power normal operating condition. 4. Discussions and conclusions A response matrix method, which directly solves the linearized system equations without Laplace transformation, is developed and applied to the single channel stability analysis of a proposed SCWR design. The sensitivity of the decay ratio (DR) to the axial mesh size is analyzed and found that the DR is not sensitive to mesh size once the axial node number is high enough such as above 40. To provide an accurate simulation, 80 axial nodes are used in this analysis. It is well known that an inlet orifice improves the channel flow stability in the condition of the subcritical pressure. To investigate the effect of the inlet orificing for the supercritical flow condition, the sensitivity of the flow stability to inlet orifice coefficient is conducted for the hot channel design. It is clearly shown that a higher inlet orifice coefficient will make the channel flow more stable, and it is found that an inlet orifice coefficient of 1.0 will satisfy the criterion of DR below 0.5. However, to have a more stable system and to account for uncertainties and off-normal conditions, the designer should select a larger value; say 10–20, for the hot channel inlet orifice coefficient. With an inlet orifice coefficient of 20.0, it is found that the hot channel DR will be below 0.5 if the power is below 121% rated power or the flow is above 83% rated flow rate. Previous analyses show that the SCWR stability is more sensitive to operating parameters than the traditional BWRs. The measures to improve the SCWR stability sensitivity feature are also investigated. It is found that a higher inlet orifice coefficient can improve the SCWR stability sensitivity feature to the operating parameters. Similar to the subcritical condition, a higher system pressure results in a more stable super-

critical flow. With a much higher equivalent inlet orifice coefficient, it is found that the average channel is very stable. A manufacturing feasibility study shows that it has no difficulties to make the inlet orifices with desired coefficients for both hot channel and average channel. This paper aims at introducing a new method for SCWR stability analysis. To assess the SCWR stability feature more accurately and thoroughly, further work need to be conducted because of following reasons: (1) The SCWR is under developing and the design may subject to change. (2) Two kinds of fuel rods supporting methods are proposed in the literatures for SCWR. One is using spacer grid, the other one is using spiral wire. In this paper, the local pressure drop due to the fuel pin supporters is not considered. (3) The fuel rods and water rods dynamics are not simulated in this paper, although it is shown that ignoring the fuel rods and water rods dynamics provides conservative results. (4) Some clearances between the subchannels in the SCWR fuel assembly are blocked by the water rods boxes, which makes the out-of-phase oscillations possible inside the fuel assembly. It is an interesting issue that needs to be further investigated. However, although the outof-phase oscillations inside the SCWR fuel assembly are possible, the magnitude of the oscillation wave will be reduced through the cross flow between the remaining clearances of the fuel rods. (5) Besides the uniform axial power shape assumed in this paper, the stability features for other power shapes will be investigated in the future. Although the model needs to be refined further, the results obtained in this paper have clearly shown that the SCWR can be designed stable and the stability sensitivity feature can be improved by carefully managing the inlet orifices and choosing proper operating parameters. Acknowledgements The authors want to extend their deepest appreciations to Prof. Kazimi at MIT and Dr. Saha at GE for their invaluable suggestions and help during this work. References Buongiorno, J., MacDonald, P.E., 2003. Supercritical Water Reactor (SCWR) Progress Report for the FY-03 Generation-IV R&D Activities for the Development of the SCWR in the U.S. INEEL/EXT-03-01210. Carmichael, L.A., Niemi, R.O., 1978. Transient and Stability Tests at Peach Bottom Atomic Power Station Unit 2 at End of Cycle 2. EPRI NP-564. Chatoorgoon, V., 2008. Supercritical flow stability in horizontal channels. Nuclear Engineering and Design 238 (8), 1940–1946. Dokhane, A., et al., 2003. Nuclear-coupled Thermal-hydraulic Nonlinear Stability Analysis Using a Novel BWR Reduced Order Model: Part 1 – The Effects of Using Drift Flux Versus Homogeneous Equilibrium Models. ICONE11-36583. Dokhane, A., et al., 2003. Nuclear-coupled Thermal-hydraulic Nonlinear Stability Analysis Using a Novel BWR Reduced Order Model: Part 2 – Stability Limits of In-phase and Out-of-phase Modes. ICONE11-36584. Gallaway, T.L., Podowski, M.Z., 2010. Stability Analysis of Fluid at Supercritical Pressure in a Heated Channel. ICAP P10-10128. Golub, G.F., van Loan, Ch.F., 1996. Matrix Computations. The Johns Hopkins University Press. Hanggi, P., 2001. Investigating BWR stability with a new linear frequency-domain method and detailed 3D neutronics. A Ph.D. dissertation submitted to the Swiss Federal Institute of Technology, Zurich. Holman, J.P., 2001. Experimental Methods for Engineers. McGraw-Hill, Boston. Ipsen, I.C.F., 1997. Computing an eigenvector with inverse iteration. SIAM Review 39 (2), 254–291. Koshizuka, S., Oka, Y., Suhwan, J., 2003. Stability Analysis of a High Temperature Reactor Cooled by Supercritical Light Water. GENES4/A NP2003-1166. Lahey Jr., R.T., Moody, F.J., 1993. The Thermal Hydraulics of a Boiling Water Nuclear Reactor. American Nuclear Society, USA. March-Leuba, J., Rey, J.M., 1993. Coupled thermohydraulic-neutronic instabilities in boiling water nuclear reactors: a review of the state of the art. Nuclear Engineering and Design 145 (1–2), 97–111. Ortega-Gómez, T., Class, A., Lahey Jr., R.T., Schulenberg, T., 2008. Stability analysis of a uniformly heated channel with supercritical water. Nuclear Engineering and Design 238, 1930–1939. Peters, G., Wilkinson, J.H, 1979. Inverse iteration, ill-conditioned equations and Newton’s method. SIAM Review 21 (3), 339–360.

J. Zhao et al. / Nuclear Engineering and Design 241 (2011) 2528–2535 Podowski, M.Z., 2003. Modeling and Analysis of Two-phase Flow Instabilities. NURETH-10. Yang, W.S., Zavaljevski, N., 2003. Effects of water rods on Supercritical Water Reactor stability. Paper 87886. In: Proc. of Global 2003. Yi, T.T., et al., 2004. A linear stability analysis of supercritical water reactors (II): coupled neutronic thermal-hydraulic stability. Journal of Nuclear Science and Technology 41 (12), 1166–1175. Zhao, J., Saha, P., Kazimi, M.S., 2004. One Dimensional Thermal-Hydraulic Stability Analysis of Supercritical Fluid Cooled Reactors. ICONE 12-49075. Zhao, J., Saha, P., Kazimi, M.S., 2007a. Hot channel stability of supercritical watercooled reactors (part I): steady state and sliding pressure startup. Nuclear Technology 158, 158–173. Zhao, J., Saha, P., Kazimi, M.S., 2007b. Hot channel stability of supercritical watercooled reactors (part II): effect of water rods heating and comparison with BWR stability. Nuclear Technology 158, 174–190. Harvey, A.H., Peskin, A.P., Klein, S.A., 2008. NIST/ASME Steam Properties Version 2.22 Users’ Guide.

2535

Jiyun Zhao is currently an assistant professor in Power Engineering Division at Nanyang Technological University in Singapore. He worked as a senior nuclear engineer for five years at AREVA NP Inc. before that he got his Ph.D. degree from Nuclear Science and Engineering Department at Massachusetts Institute of Technology. C.P. Tso is currently a visiting professor in Thermal and Fluids Division at Nanyang Technological University in Singapore. He obtained his Ph.D. degree in nuclear engineering from University of California, Berkeley, and his Master’s degree in nuclear engineering from Massachusetts Institute of Technology. K.J. Tseng is currently an associate professor and division head of Power Engineering Division at Nanyang Technological University (NTU) in Singapore. He received the B.Eng. (first class) and M.Eng. degrees from the National University of Singapore, and the Ph.D. degree from University of Cambridge. He manages a number of research programs in the area of electrical energy and power. He is currently the director of Center for Smart Energy Systems at NTU.