PARCS to BWR stability analysis

Annals of Nuclear Energy 36 (2009) 317–323

Contents lists available at ScienceDirect

Annals of Nuclear Energy journal homepage: www.elsevier.com/locate/anucene

Application of TRACE/PARCS to BWR stability analysis Yunlin Xu a,*, Thomas Downar a, R. Walls b, K. Ivanov b, J. Staudenmeier c, J. March-Lueba d a

Purdue University, West Lafayette, IN, USA Penn State University, State College, PA, USA US NRC, Rockville, MD, USA d Oak Ridge National Laboratory, Oak Ridge, TN, USA b c

a r t i c l e

i n f o

Article history: Received 10 October 2008 Received in revised form 12 December 2008 Accepted 17 December 2008 Available online 23 February 2008

a b s t r a c t The work described here is the validation of TRACE/PARCS for Boiling Water Reactor stability analysis. A stability methodology was previously developed, verified, and validated using data from the OECD Ringhals stability benchmark. The work performed here describes the application of TRACE/PARCS to all the stability test points from cycle 14 of the Ringhals benchmark. The benchmark points from cycle 14 were performed using a half-core symmetric, 325 channel TRACE model. Several parametric studies are performed on test point 10 of cycle 14. Two temporal difference methods, Semi-Implicit method (SI) and Stability Enhanced Two Step (SETS) method are applied to three different mesh sizes in heated channels with series of time step sizes. The results show that the SI method has a smaller numerical damping than the SETS method. When applying the SI method with adjusted mesh and Courant time step sizes (the largest time step size under the Courant limit), the numerical damping is minimized, and the predicted Decay Ratio (DR) agrees well with the reference values which were obtained from the measured noise signal. The SI method with adjusted mesh and Courant time step size is then applied to all test points of cycle 14 with three types of initiating perturbations, control rod (CR), pressure perturbation, and noise simulation (NS). There is good agreement between the decay ratios and frequencies predicted by TRACE/ PARCS and those from the plant measurements. Sensitivities were also performed to investigate the impact on the decay ratio and natural frequency of the heat conductivity of the gap between fuel and clad, as well as the impact of the pressure loss coefficient of spacers. Ó 2009 Published by Elsevier Ltd.

1. Introduction The US NRC reactor system analysis code TRACE (TRAC RELAP5 Advanced Computational Engine) (Spore et al., 2000) is used to study the reactor coolant system under a wide variety of flow conditions including multi-phase thermal-hydraulics. Multi-dimensional time dependent power distributions are determined using the (Purdue Advanced Reactor Core Simulator) PARCS ( Downar et al., 2006) multi-dimensional reactor kinetics code which has been coupled to TRACE. TRACE coupled to PARCS has been validated previously for Pressurized Water Reactor transient analysis using the OECD/NEA Main Steam Line Break (MSLB) Benchmark (Kozlowski et al., 2004). The work described here is the validation of TRACE for Boiling Water Reactor stability analysis. A stability methodology was previously developed, verified, and validated using data from the OECD Ringhals stability benchmark (Xu et al, 2005a,b). The work performed here describes the application of TRACE to the stability test points from cycle 14 of the Ringhals benchmark and a preliminary

* Corresponding author. Tel.: +1 765 496 3853. E-mail address: [email protected] (Y. Xu). 0306-4549/$ - see front matter Ó 2009 Published by Elsevier Ltd. doi:10.1016/j.anucene.2008.12.022

investigation of the TRACE numerical performance during stability analysis. The Ringhals-1 (RH1) core was well instrumented and the benchmark (Lefvert, 1994) contains a comprehensive and well-defined set of time series data from measurements at beginning of the cycles (BOC) 14, 15, 16 and 17, as well as in the middle of the cycle 16. In this work, we focus on the ‘‘in-phase” instability of test points in BOC 14. The reactor core consists with 648 fuel assemblies. Three types of assemblies, 8  8 assembly with 63 fuel rods, SVEA assembly with 63 or 64 fuel rods, are loaded in cycles 14–17. The number of assemblies for each type loaded in cycle 14 are 145 8  8 assemblies, 183 SVEA(63) assemblies and 320 SVEA(64) assemblies. The nominal core power is 2270 MW. The LPRM and APRM detector signals were sampled for 660 s with 80 ms interval which yielded 8250 data points for each channel at each test condition. The reference decay ratio and frequency (Table 1) are abstracted from measured APRM noise signals using the ARMA process (Lefvert, 1996). This paper is organized as follows. The methodologies for BWR stability analysis will be introduced in Section 2. The development of the TRACE and PARCS models for RH1 will be described in the Section 3. The parametric study of nodalization and time step size are presented in Section 4. The results of the benchmark obtained

318

Y. Xu et al. / Annals of Nuclear Energy 36 (2009) 317–323

Table 1 Stability data for cycle 14 of Ringhals. Points

Core power (%)

Core flow rate (kg/s)

Decay ratio

Frequency (Hz)

P01 P03 P04 P05 P06 P08 P09 P10

65.0 65.0 70.0 70.0 70.2 75.1 72.6 77.7

4105 3666 3657 3868 4126 3884 3694 4104

0.30 0.69 0.79 0.67 0.64 0.78 0.80 0.71

0.43 0.43 0.55 0.51 0.52 0.52 0.56 0.50

with the TRACE/PARCS coupled code are presented in Section 5. The sensitivities of gap conductivity and spacer loss coefficients are studied in Section 6 Summary and conclusions are given in Section 7. 2. BWR stability analysis methodology This section describes the methodology and software tools developed for the application of TRACE/PARCS for BWR stability analysis. The stability simulation calculation with TRACE/PARCS begins with a stand-alone TRACE initialization of the thermalhydraulics field using a fixed power distribution. Upon convergence, a steady state calculation is then performed with the coupled code TRACE/PARCS to achieve a converged coupled thermalhydraulic/neutron flux solution. The SETS numerical method (Borkowski et al., 1996) is used in the TRACE steady state calculations and an explicit coupling is used between the TRACE and PARCS solutions. The stability simulation is then performed using the semi-implicit method in TRACE by introducing one of three excitations to the core to initiate the transient. These perturbations are: (1) Control rod (neutronics) perturbation: one or more control rods are temporarily moved and then replaced to its original position in the PARCS model. (2) Pressure (thermal-hydraulic) perturbation: a short perturbation is applied to the steam line of the TRACE model. (3) Simulation of the in-core noise: the instantaneous moderator density feedback is perturbed using a superposition of the fundamental, first azimuthal, and random noise components. The noise technique numerically simulated the actual stability monitoring procedure utilized during the Ringhals stability tests. The simulation of noise required implementation of new modules and modifications to PARCS. The approach used was similar to that proposed by Hotta et al. (2001) in which a numerically generated band white noise is imposed using the feedback moderator density in each neutronics node prior to the construction of the group cross-sections. Subroutines were added to PARCS in which a white noise is constructed as a superposition of the fundamental, first azimuthal and random components of the flux distribution. For out-of-phase (regional) oscillations all three components are activated, whereas for in-phase (core-wise) oscillations only the fundamental and random components are used. Each component of the noise is expressed as a product of the amplitude and shape functions. The post processing code, DRARMAX, which was developed at Purdue University, is selected to evaluate the DR and NF from the time series output of the TRACE/PARCS calculation. DRARMAX evaluates the decay ratio and the natural frequency for signals obtained from transients initiated by short perturbations as well as those obtained from noise simulation during stationary operation. DRAMAX

is able to calculate the DR and NF taking into consideration also the signal constituting the perturbation given to the system. DRARMAX also determines DRs and NFs for three types of functions: the output signal, its autocorrelation function (ACF) and the impulse response function (IRF). For each of the three types of functions DRARMAX fits the function with a damped sinusoidal form:

xðtÞ ¼ c0 þ ceat cos ðxt þ uÞ

ð1Þ

The first 2 s of the calculated time trace are not used for the fit because they have contamination from higher-order modes and the particular shape of the perturbation used. If least square fitting is successful, the DR and NF of the system are given as

DRf ¼ e2pa=x

ð2Þ

Frf ¼ x=ð2pÞ

ð3Þ

This definition of DR and NF is appropriate for transients with short perturbations, such as control rod and pressure perturbations. For the measured or simulated noise signals, DRARMAX code will first extract the impulse response function (IRF) by an autoregressive and moving-average (ARMA) process (Marple, 1987):

xðiÞ ¼ 

p X

ak xði  kÞ þ

k¼1

q X

bm uði  mÞ

ð4Þ

m¼0

where u, x are input and output noise signals, p, q are orders of the ARMA process, and coefficients ak, bm are obtained by the prediction-error identification method (PEM). An index referred to as Final Prediction Error (FPE) is defined to estimate the quality of the ARMA model:

FPE ¼

Nþqþ1 kx  ^xk N  2p  q  1

ð5Þ

where N is the number of sampled data points and ^x is approximation of the output signal by ARMA. The optimal ARMA orders which minimize the FPE were used for all noise signals. It should be noted that the IRFs were very stable with respect to ARMA orders near the optimal ARMA value. The Lyapunov DR and NF are used for noise signals:

DRL ¼ e2pReðcÞ=ImðcÞ FrL ¼ ImðcÞ=ð2pDtÞ

ð6Þ ð7Þ

where Dt is the time interval of input and output time series, and P c = ln(bm) and bm is the dominant root of 1 þ pk¼1 ak bk (Munoz et al., 1992). 3. TRACE/PARCS models for Ringhals This section provides an overview of the development and validation the TRACE thermal-hydraulics and PARCS neutronics models. The PARCS core neutronics model was developed using the neutronics data provided in the benchmark specifications. The core nodalization is based on one node per assembly in the radial plane (32  32 planar nodes) and 27 axial nodes. The nodalization includes explicit modelling of both the radial and axial reflectors. The PARCS neutronics solution for both the steady state and transient were based on the two-group nodal diffusion option. The benchmark cross-section library specified for cycle 14 of RH1 was based on the CASMO-3/SIMULATE-3 codes and it was reformatted to comply with the cross-section format required in PARCS. The node-wise burnup and history values for all test points from the core simulator SIMULATE-3 are used in PARCS model. The TRACE model includes the reactor pressure vessel and recirculation loops. The boundaries are at the main steam line and feed

319

Y. Xu et al. / Annals of Nuclear Energy 36 (2009) 317–323

water pipes. The separator and downcomers are modelled within the vessel. One of the most important issues for coupled code analysis is the use of sufficient thermal-hydraulics (T–H) channels such that the respective field equations can be solved with sufficient accuracy for the particular application. As depicted in Fig. 1, a 325 channel model with half-core symmetry model was used in the work reported in this paper on the in-phase stability analysis. This model significantly reduces the running time compared with a one to one, full core model which is required for out-of-phase oscillations. In addition, each T–H channel also modelled a leakage path from the channel to the vessel along the bypass region. This leakage simulates the external core flow bypass, which typically is about 10% of the total core coolant flow. All the material properties, such as heat conductivities and heat capacities as well as densities for fuel pellets and claddings from the specifications (Lefvert, 1994) were used directly in the TRACE model. One of the important parameters for stability analysis is the gap conductance which is specified to depend on the fuel pellet average temperature as follows (Table 2). Most of the hydraulic parameters in the specifications were also used in the TRACE model, such as the channel and fuel rod dimensions, flow areas, hydraulic diameters and pressure loss coefficients for channel orifices, lower and upper core tie plates. The TRACE model used the spacer pressure loss coefficient from the Ringahls specification for the 8  8 type of assemblies. The specification, however, defined SVEA spacer coefficient values that were 10% higher than for 8  8 fuel assembly. Other references (see for example Taleyarkan, 1986) define the SVEA spacer coefficients as significantly lower than standard 8  8. A sensitivity calculation was performed using TRACE and a series of SVEA spacer coefficients until a good match was found between the predicted 8  8 and SVEA channel flow distributions. Through this sensitivity analysis, the final SVEA best-estimate spacer coefficient was determined to be approximately 35% loser than the coefficients used for 8  8 bundles. In Section 5 of this paper, the specified spacer coefficients were used for SVEA type assemblies in order to compare the results with our best-estimate value.

Table 2 Gap conductance specifications. Pellet average temperature (°C)

Gap conductance (W/m2 °C)

400 500 600 700 800 900

3800 4000 4300 4800 5400 6000

4. Spatial and temporal discretization Previous studies (Mahaffy, 1983, 1993) have suggested the potential for numerical damping in the SETS method (Mahaffy, 1982), which is the primary numerical solution method in TRACE. Because of these concerns, the semi-implicit (SI) method option in TRACE was used for stability analysis in the work here. However, when the semi-implicit method is used for stability analysis, it is important to achieve a Courant number near unity in order to minimize numerical dissipation:

c ¼ mDt=Dx

ð8Þ

where v is velocity, Dx is mesh size and Dt is time step size. The potential numerical dissipation in the SI fluid solution can be demonstrated using the single phase mass continuity equation:

qt þ mqx ¼ 0

ð9Þ

The explicit and implicit upwind finite differencing of Eq. (9) can be written as Eqs. (10) and (11), respectively.

qjnþ1  qnj

v j1 qnj1  v j qnj

¼0 Dt Dx nþ1 nþ1 n q  qj v j1 qj1  v j qj þ ¼0 Dt Dx þ

nþ1 j

ð10Þ ð11Þ

The modified equation for the explicit and implicit upwind finite difference of mass continuity equation can be written as

Fig. 1. TRACE 325 channel model for in-phase stability analysis.

Y. Xu et al. / Annals of Nuclear Energy 36 (2009) 317–323

 xx =2 þ O½Dx2 ; Dt2  qt þ v qx ¼ v Dxð1þcÞq

ð12Þ

where the right hand side of Eq. (12) is the error introduced by finite differencing. The leading term is the artificial viscosity and the sign  is negative for the explicit scheme and positive for the implicit scheme. The numerical diffusion of the explicit or implicit upwind difference scheme is the result of two sources of error. The spatial differencing error can be written as

v Dxqxx =2

ð13Þ

and the time differencing error can be written as

 v Dxqxx =2 þDtqtt =2 ¼ þc

Table 3 Cell lengths for non-uniform Mesh. Cells indicesa

Length (mm)

1 2–6 7–10 11–14 15–16 17–18 19–20 21–22 23–25

147.2b 58.88 73.6 110.4 147.2 184 220.8 257.6 294.4

a

ð14Þ

In the implicit method, the two errors always enhance each other, and the minimum dissipation occurs at the smallest time step size. However, in the explicit method, the two errors cancel each other. If c = 1, numerical dissipation would disappear. However, if c < 1 then numerical dissipation will increase with decreasing time step, and if c > 1, then numerical dissipation is negative and false viscosity will destabilize the solution. If the upwind scheme is used for spatial finite differencing, then the implicit method is numerically more stable than the explicit method. However, the implicit method has potentially more numerical dissipation then the explicit method for practical applications. For stability analysis in the work here, the explicit method with larger time step sizes was preferred. However, the time step size was chosen to match the Courant limit as closely as possible. In this manner, the potential for numerical instability was minimized which would compromise the accuracy of the solution. It should be noted that the modified Eq. (12) also suggests that numerical diffusion would be reduced if the spatial mesh size is reduced. However, several references (see Dinh et al., 2003) argue that the two fluid model is ill-posed and spatial convergence cannot generally be achieved. Nonetheless, spatial convergence analysis was not the subject of this paper and further mesh refinement was not considered since the computational burden of a full system model with more than 300 channels would have prohibited practical stability analysis with TRACE. The TRACE two fluid model solves six differential equations and there is no pure explicit or implicit method in TRACE for temporal finite differencing. There are two temporal finite difference methods in TRACE: the semi-implicit method and the Stability Enhanced Two Step (SETS) method. In the SI method, the pressure equation is solved implicitly, and the density, velocity and energy equations are solved explicitly. Therefore, the numerical dissipation in SI is similar to the explicit method and the time step size can exceed the acoustic (pressure velocity propagation) Courant limit. Because the velocities in each cell of the TRACE model are generally different due to differences in the cell void fraction, an adjustable axial mesh size was introduced for stability modelling to allow the material Courant numbers in most cells to approach unity simultaneously and therefore to minimize the numerical dissipation. As the transients simulated for the stability analysis are caused by very small perturbations, the velocities of steam vary only slightly during the transient. The approach used here was to determine a mesh size for each cell in the active core region according to the steady state solution which would allow the material Courant numbers in most cells to approach unity simultaneously. In the model used here, the active core is discretized using 25 axial meshes which would result in an axial mesh size of 147.2 mm for uniform meshing. For the case of non-uniform meshing, the results from a previous study (Xu et al., 2005b) were used which approximately matched the channel steam velocity. The chosen axial cell lengths were as follows (Table 3),

b

The cell indices are ordered from bottom to top of the active core. The first cell is in the non-boiling region.

The decay ratios from the transient simulation of test point 10 of Ringhals-1 cycle 14 for the SI and SETS methods are shown in Fig. 2 for six different time step sizes from 1 ms to 37 ms. (Note: 37 ms is the material Courant limit for the adjusted axial mesh. Since the semi-implicit method is unstable when the Courant number exceeds unity, the time step in TRACE is not permitted to exceed 37 ms). For comparison, the decay ratio from the SI method at the Courant limited step size, 19 ms, is also shown in the Figure. The measured decay ratio for this test point of 0.71 is also shown in the Figure. The results in Fig. 2 indicate that for the SI method with adjusted axial mesh size the numerical damping decreases when the time step size is increased. The numerical damping for SI is minimized when the step size reaches material Courant limits (about 37 ms). Conversely, for the SETS method the numerical damping decreases when the step size is decreased and the numerical damping is largest at the material Courant limit. The results from the SETS and the semi-implicit method are similar when the time step size is very small. The results in Fig. 2 are consistent with the modified equation analysis and indicate that the non-uniform mesh help to reduce numerical damping. 5. Stability results for test points of RH1 cycle 14 The methodology for BWR stability analysis described above was applied to all 8 test points of RH1 cycle 14. Detailed results will first be discussed for point 10 for both the steady state and transient and then results will be summarized for all the points. The steady state results for test point 10 are in good agreement with the plant data summarized in the benchmark specifications. PARCS/TRACE results are also consistent with SIMULATE-3 solutions, as the cross-sections of assemblies and core bunrup distributions are provided from CASMO-3/SIMULATE-3 code suite in the benchmark specification. Fig. 3 shows the channel flow rates com-

0.75 0.70

Decay Ratio

320

reference

0.65

SI (adjusted mesh)

0.60 SETS (adjusted mesh)

0.55

SI (uniform mesh)

0.50 0

10

20

30

Step Size (ms) Fig. 2. Decay ratios for various time step sizes.

40

321

Y. Xu et al. / Annals of Nuclear Energy 36 (2009) 317–323

79

7

Sample Data Least-square fitted

8x8 central

78.5

6 SVEA central

SVEA semi-peripheral

5.5

Sample Data

Flow rate kg/s

6.5

8x8 semi-peripheral

5 SVEA peripheral

4.5

8x8 peripherav

4 3.5

0

0.2

0.4

0.6

1

1.2

1.4

1.6

Relative Power

void fraction

0.6 0.5 0.4 0.3

Specification

0.2

Simulate-3

0.1

TRACE

0.0 200

300

height(cm)

Axial Power Profile

Fig. 4. C14P10 axial void fraction.

Reference Simulate-3 TRACE

1

3

5

7

9 11 13 15 17 19 21 23 25

Axial Plane

Reactivity (c) or change(%)

Fig. 5. C14P10 axial power profile.

1.5

Power change Heat flux change Doppler reactivity Coolant reactivity

1.0 0.5

5

10

6

8

10 12 14 16 18 20

6. Sensitivity studies

0.0 0

4

pared to the specifications and Figs. 4 and 5 show the axial void and power profiles compared with specifications and the Simulate-3 solution, respectively. The stability transient results initiated with a control rod perturbation are shown in Fig. 6. This figure shows that the heat flux at fuel rod surface is delayed from the power change by about 0.48 s. The Doppler reactivity feedback which helps to stabilize the transient is in the same phase as the surface heat flux. The coolant reactivity feedback caused by changes in the coolant density is driven by the surface heat flux change with about 0.57 s lag. The power change follows almost immediately changes in the coolant reactivity. The delay in the change in the surface heat flux due to changes in the power is determined by the heat conduction time constant of the fuel rods. And the lag in the change of the coolant density due to the change in the surface heat flux is determined by the flow velocity, the heated length, and the axial power profile. These two time lags make up half of the oscillation period. The core power trace and its least squares fitting is shown in Fig. 7 for the stability transient initiated by a pressure perturbation. And the core power trace predicted by noise simulation and its IRF are shown in Figs. 8 and 9, respectively. As summarized in Table 4, the decay ratios and oscillation frequencies predicted from control rod, pressure perturbation, and noise simulation are very similar. The solutions of all test points at cycle 14 of Ringhals-1 are also summarized in Table 4. The average channel flows for 8  8 type of assemblies are higher by 2.5–3% than the specifications while the average channel flow for SVEA type of fuel assemblies are about same as specifications. The power distributions also agree well with the specification. The RMS of the radial power differences are about 3%. And with the exception of point 9, the RMS of the axial and nodal wise power differences are under 8.5% and 9.8%, respectively. The predicted decay ratios are higher for all points except point 3 and 6, the differences are about 16–13% with the exception of point 1. The predicted natural oscillation frequency are all lower by about 5–9%. In all cases, the decay ratios and natural frequencies from the three methods of initiating the transients are consistent.

0.7

1.6 1.4 1.2 1 0.8 0.6 0.4 0.2 0

2

Fig. 7. Core power after pressure perturbation.

0.8

100

76.5 0

Time (s)

Fig. 3. C14P10 channel flow rates.

0

77.5 77

TRACE flow Refe. flow

0.8

78

15

-0.5 -1.0 -1.5

Time (second) Fig. 6. C14P10 control rod perturbation transient.

20

Sensitivity studies were performed on the pressure loss coefficients of the spacer grids and on the heat conductivity of the gap between the fuel pellet and cladding. It was determined these were two of the most important parameters impacting the stability results. The loss coefficients of the SVEA spacer grids provided in the specifications was substantially different from the values provided in other references (e.g., Taleyarkan, 1986). All of test points in the cycle 14 were reanalyzed with the higher value provided in the

322

Y. Xu et al. / Annals of Nuclear Energy 36 (2009) 317–323

86

Table 5 Stability solutions for RH1 C14 with higher pressure loss coefficients for the SVEA spacer grid.

84

Points

88

Sample Data

Sample Data

82 P01 P03 P04 P05 P06 P08 P09 P10

80 78 76 74 72 0

50

100

150

200

250

Flow error (%)

Decay ratio

88

SVEA

CR

PP

Frequency CR

PP

3.89 3.90 4.11 3.85 3.91 3.86 3.82 3.94

1.72 1.31 1.34 1.49 1.72 1.57 1.30 1.81

0.465 0.673 0.928 0.848 0.749 0.981 1.026 0.825

0.515 0.679 0.915 0.839 0.746 0.968 1.018 0.836

0.406 0.395 0.505 0.485 0.478 0.476 0.523 0.473

0.409 0.395 0.504 0.483 0.477 0.474 0.522 0.475

300

Fig. 8. Core power from noise simulation.

Decay Ratio

1 IRF--ARMA(101, 0)

0.8

IRF

0.6 0.4 0.2 0

0.85

0.85

0.80

0.80

0.75

0.75

0.70

0.70

0.65

0.60

Nature frequency

0.55

0.55

0.50

0.50

0.45

0.45

0.40 0.6

0.8

-0.2 -0.4

0.65

Decay Ratio

0.60

1

1.2

1.4

nature frequency (Hz)

Time (s)

0.40 1.6

Gap conductivity multiplier 0

2

4

6

8

Fig. 10. Sensitivity on gap conductivity.

10 12 14 16 18 20

Time (s) Reactivity (c) or change(%)

Fig. 9. IRF for C14P10 core power.

specifications. The results of this sensitivity analysis are summarized in Table 5. Comparison of the results in Tables 4 and 5, indicate the channel flow rates for the SVEA assemblies are lower with the increased loss coefficients of the specifications, while the flow rates for the SVEA assemblies with our best-estimate value are much closer to the specifications. The flow rates for the 8  8 assemblies also have larger differences from the specifications with the increased SVEA loss coefficients. As indicated in Table 5, the decay ratios with the higher spacer loss coefficients are significantly higher and have large differences from the specifications. The natural frequencies are not significantly affected by changes in the loss coefficients for the in-phase oscillations analyzed here.

1.5

Power change Heat flux change Doppler reactivity Coolant reactivity

1.0 0.5 0.0 0

5

10

15

20

-0.5 -1.0 -1.5

Time (second) Fig. 11. Transient results for 1.44 times of specified gap conductivity.

Table 4 Summary of stability solutions for all test points at cycle 14 of Ringhals-1. Points

P01 P03 P04 P05 P06 P08 P09 P10 a b c d

Flow error (%)a

Power RMS (%)a

88

SVEA

Axial

Radial

Nodal

Ref.

CRb

PPc

NSd

Ref.

CRb

PPc

NSd

2.21 2.23 2.34 2.08 2.14 2.07 2.05 2.13

0.88 0.42 0.35 0.55 0.81 0.60 0.32 0.87

6.77 5.89 8.45 6.73 0.60 4.85 13.5 3.23

2.73 2.90 3.40 3.14 3.02 2.88 2.25 2.66

9.14 8.19 9.80 8.41 7.95 6.72 14.4 5.68

0.30 0.69 0.79 0.67 0.64 0.78 0.80 0.71

0.377 0.580 0.833 0.741 0.619 0.870 0.896 0.733

0.391 0.582 0.813 0.757 0.622 0.854 0.891 0.724

0.444 0.616 0.677 0.736 0.647 0.851 0.803 0.729

0.43 0.43 0.55 0.51 0.52 0.52 0.56 0.50

0.382 0.391 0.505 0.487 0.481 0.477 0.526 0.476

0.382 0.391 0.505 0.483 0.477 0.476 0.524 0.474

0.395 0.391 0.511 0.484 0.474 0.475 0.532 0.475

Compared with reference values. Control rod perturbation. Pressure perturbation. Noise simulation.

Decay ratio

Natural frequency

Oscillation period or heat flux delay time (second)

Y. Xu et al. / Annals of Nuclear Energy 36 (2009) 317–323

2.5 2.0 Oscillation period

1.5

delay of heat flux change

1.0 0.5 0.0 0.6

0.8

1

1.2

1.4

1.6

Gap conductivity multiplier Fig. 12. The impact of gap conductivity to time constants.

Because of potentially large variations in the gap conductance, a sensitivity study was performed to determine its impact on the stability results. The specified fuel rod temperature dependent gap conductivity was multiplied by constant factors of 0.64, 0.8, 1.2 and 1.44 for four of the test points. The decay ratios and natural frequencies for these 4 test points as well as the base cases for test point 10 are shown in Fig. 10. As indicated, the decay ratio and the natural frequency increase with increases in the gap conductivity. The reason for these increases is apparent in Fig. 11 which shows the detailed transient results for the case of the gap conductance multiplied by 1.44 times the specified value. By comparison with the base case in Fig. 6, it is apparent that the response of the fuel rod surface flux and coolant reactivity due to changes in the core power are significantly larger when the gap conductance is increased. Both of these changes destabilize the system and lead to a larger decay ratio. The impact of the gap conductivity on the time delay of the surface heat flux change from changes in the core power is shown in Fig. 12. The oscillation periods are also shown in the figure which shows clearly the correlation between the period and heat flux delay. 7. Summary TRACE/PARCS was applied to the stability analysis of cycle 14 of the Ringhals benchmark. The work here focused on the in-phase (i.e., core-wide) oscillation and used a half-core symmetric, 325 channel TRACE model. A detailed numerical study was performed on the SI and SETS numerical integration methods. Parametrics on the axial nodalization and time step size were performed using test point 10 of cycle 14. When the SI method is applied with a variable axial mesh and the Courant time step size (the largest time step size under the Courant limit), the numerical damping is minimized and the predicted Decay Ratio (DR) agrees well with the measured values. The SI method with an adjusted mesh and the Courant time step size was then applied to all test points of cycle 14. Three types of initiating perturbations were used for each of the test points: control rod (CR), pressure perturbation, and noise simulation (NS). The

323

decay ratios and frequencies are in good agreement with the measured values which were extracted from noise simulation during the plant tests. A sensitivity study has been performed to examine the impact of variations in the gap conductance and the pressure loss coefficient of the spacer grids on the decay ratio and natural oscillation frequency. The results indicate that the spacer loss coefficient has a significant impact on the decay ratio but does not significantly impact the natural frequency. However, the gap conductance was found to have an impact on both the decay ratio and the natural frequency. The results here indicate that time domain stability analysis can be successfully performed with the TRACE code and work is currently being completed on all of the test points of cycles 15, 16 and 17 of the OECD Ringhals benchmark, to include both out-ofphase as well as in-phase oscillations. While the results here are encouraging, the use of problem specific time and spatial discretization methods to minimize numerical damping is not entirely satisfying. An important area of future research will be the development and implementation of robust higher-order numerical methods that simplify the modelling and improve the reliability of BWR time domain stability analysis. References Borkowski, J. et al., 1996. SIMULATE-3K simulation of the Ringhals 1 BWR stability measurements. In: PHYSOR96, Mito, Japan. Dinh, T. et al. 2003.Understanding the ill-posed two-fluid model. In: NURETH-10, October, Seoul, Korea. Downar, T., Xu, Y., Seker, V., Ward, A., 2006. PARCS: Purdue advance reactor core simulator. In: Proceedings of ANS Reactor Physics Topical Meeting PHYSOR 6. Vancouver, BC. Hotta, A., Zhang, M., Ninokata, H., 2001. BWR regional instability analysis by TRAC/ BF1-ENTREE II: application to Ringhals unit stability Test. Nucl. Technol. 135 (1), 17–38. Kozlowski, T., Downar, T., et al., 2004. Analysis of the OECD/NEA PWR main steam line break benchmark with TRAC-M/PARCS. Nucl. Technol.. Lefvert, T., 1994. Ringhals 1 stability benchmark, specifications. NEA/NSC/DOC(94) 15. Lefvert, T., 1996. Ringhals 1 stability benchmark. Final Report. NEA/NSC/DOC(96)22. Mahaffy, J.H., 1982. A stability-enhancing two-step method for fluid flow calculations. J. Comput. Phys. 46, 329–341. Mahaffy J.H., 1983. The advantages and limitations of the SETS method. In: International Conference on Numerical Methods in Nuclear Engineering, Canadian Nuclear Society, Montreal, Quebec, pp28–39. Mahaffy, J.H., 1993. Numerics of codes: stability, diffusion, and convergence. Nucl. Eng. Des. 145, 131–145. Marple, S.L., 1987. Digital Spectral Analysis with Applications. Prentice Hall, Englewood Cliffs (Chapter 8). Munoz, J.L., Verdù, G., Pereira, C., 1992. Dynamic reconstruction and lyapunov exponents from time series data in boiling water reactor. Application to BWR stability analysis. Annu. Nucl. Energy 19, 223. Spore, J. et al., 2000. LA-UR-00–910 TRAC-M/FORTRAN 90. Theory Manual. Los Alamos National Laboratory, Los Alamos, NM. Taleyarkan, R., 1986. MAZDA-NF – A Parallel Channel Stability Analysis Code for BWR Fuel Assemblies. WCAP-11146 Westinghouse Nuclear Energy Systems. Xu, Y., Downar, T., Ivanov, K., Petruzzi, A., Maggini, F., Miró, R., Staudenmeier, J., 2005a. Methodolgies for BWR stability analysis with TRACE/PARCS. ANS Mathematics and Computation Meeting, Avignon, France. Xu, Y., Downar, T., Ivanov, K., Vedovi, J., Petruzzi, A., Staudenmeier, J., 2005b. Analysis of the OECD/NEA Ringhalls instability benchmark with TRACE/PARCS. ANS Mathematics and Computation Meeting, Avignon, France.