Numerical instability of time-discretized one-point kinetic equations

Annals of Nuclear Energy 27 (2000) 791±803 www.elsevier.com/locate/anucene

Numerical instability of time-discretized one-point kinetic equations Kengo Hashimoto a,*, Hideaki Ikeda b, Toshikazu Takeda c a

Atomic Energy Research Institute, Kinki University, 3-4-1 Kowakae, Higashi-Osaka, Osaka 577-8502, Japan b In-Core Fuel Management System Department, Toden Software Inc., 6-19-15 Shinbashi, Minato-ku, Tokyo 105-0004, Japan c Department of Nuclear Engineering, Graduate School of Engineering, Osaka University, 2-1 Yamada-oka, Suita, Osaka 565-0871, Japan Received 24 June 1999; accepted 26 July 1999

Abstract The one-point kinetic equations with numerical errors induced by the explicit, implicit and Crank±Nicolson integration methods are derived. The zero-power transfer functions based on the present equations are demonstrated to investigate the numerical stability of the discretized systems. These demonstrations indicate unconditional stability for the implicit and Crank± Nicolson methods but present the possibility of numerical instability for the explicit method. An upper limit of time mesh spacing for the stability is formulated and several numerical calculations are made to con®rm the validity of this formula. # 2000 Elsevier Science Ltd. All rights reserved. Keywords: Numerical stability; Finite-di€erence equations; One-point kinetic equations; Mesh spacings

1. Introduction In the choice of numerical time-integration scheme for an initial value problem, not only the accuracy but the numerical stability should be considered. When a numerical scheme is applied to an initial value problem, a pseudo oscillation may appear in the solution. This oscillation originates from the numerical instability of the discretized di€erential equations and is physically of no signi®cance. For the avoidance of the instability, numerical stability of the applied scheme must be

* Corresponding author. Fax: +81-6-6721-3743. E-mail address: [email protected] (K. Hashimoto) 0306-4549/00/$ - see front matter # 2000 Elsevier Science Ltd. All rights reserved. PII: S0306-4549(99)00092-4

792

K. Hashimoto et al. / Annals of Nuclear Energy 27 (2000) 791±803

known in advance. By the application of a more stable scheme or suciently ®ne discretizing meshes, a physically reasonable solution can be obtained. In the ®eld of computational ¯uid dynamics, the numerical stability of ®nite-di€erence equations obtained by various time- and space-integration schemes has been investigated (e.g. Roache, 1976). In these investigations, the stability of each scheme was estimated and the stability limit of time-mesh spacing was determined. In the ®eld of reactor physics, however, the numerical stabilities of discretized neutron equations are still not completely clari®ed. Usually, on the basis of experience or a trial-and-error approach, a numerical scheme and a set of mesh spacings are chosen to obtain a stable solution. We have attempted to estimate theoretically the numerical stability of reactor kinetics, neutron di€usion and transport equations discretized by various time- and space-integration schemes. In this paper, as a ®rst step of these studies, the stabilities of one-point kinetic equations discretized by explicit, implicit and Crank±Nicolson (modi®ed Euler) methods are estimated. Even in the investigations of space-dependent equations, the stability of the above equations is fundamental. Time-integration methods employed are described in Section 2, and the one-point equations including the errors induced by these methods are derived in Section 3. The zero-power transfer functions based on these equations are demonstrated in Section 4. The stability limit of time-mesh spacing for the explicit method is given in Section 5. 2. Time integration method In this study, we consider a discretization in time for the following one-point kinetic equations: dP…t† …t† ÿ ˆ P…t† ‡ lC…t† dt 

…1a†

dC…t† ˆ P…t† ÿ lC…t† dt 

…1b†

where P…t† C…t† …t† l 

= = = = = =

power or neutron density; delayed-neutron precursors concentration; reactivity; decay constant of delayed neutron precursors; delayed-neutron fraction; prompt-neutron generation time.

To simplify the following formulation, the above kinetic equations are described under one delayed-group theory. The ®nite-di€erence forms discretized by several time-integration methods are described below.

K. Hashimoto et al. / Annals of Nuclear Energy 27 (2000) 791±803

793

2.1. Explicit-time integration method Pn‡1 ÿ Pn n ÿ n ˆ P ‡ lCn
t 

…2a†

Cn‡1 ÿ Cn n ˆ P ÿ lCn
t 

…2b†

where
t ˆ tn‡1 ÿ tn ; Pn ˆ P…tn †; Cn ˆ C…tn †; n ˆ …tn †; and tn denotes the time at nth mesh point. 2.2. Implicit time-integration method Pn‡1 ÿ Pn n‡1 ÿ n‡1 P ‡ lCn‡1 ˆ 
t

…3a†

Cn‡1 ÿ Cn n‡1 ˆ P ÿ lCn‡1
t 

…3b†

2.3. Crank±Nicolson method Pn‡1 ÿ Pn Pn‡1 ‡ Pn Cn‡1 ‡ Cn n‡1 Pn‡1 ‡ n Pn ˆÿ ‡l ‡
t 2 2 2 

…4a†

Cn‡1 ÿ Cn Pn‡1 ‡ Pn Cn‡1 ‡ Cn ˆ ÿl
t 2 2 

…4b†

3. Kinetic equations with discretizing error In this section, we derive an expression for the one-point equations with discretizing error (i.e. truncation error) according to an approach of Hirt (1968). The derivation is made for each integration method described in the preceding section. 3.1. Expression for the explicit method as

The Taylor-series expansions of unknowns at a time-mesh point can be described

794

K. Hashimoto et al. / Annals of Nuclear Energy 27 (2000) 791±803

Pn‡1 7Pn ‡
tdP…Tn †=dt ‡
t2 =2d2 P…tn †=dt2 ‡
t3 =6d3 P…tn †=dt3

…5a†

Cn‡1 7Cn ‡
tdC…tn †=dt ‡
t2 =2d2 C…tn †=dt2 ‡
t3 =6d3 C…tn †=dt3

…5b†

where the expansion is truncated in the third order. The substitution of the above equations into Eqs. (2a) and (2b) gives dP…tn †=dt ˆ …n ÿ †=Pn ‡ lCn ÿ
t=2d2 P…tn †=dt2 ÿ
t2 =6d3 P…tn †=dt3

…6a†

dC…tn †=dt ˆ =Pn ÿ lCn ÿ
t=2d2 C…tn †=dt2 ÿ
t2 =6d3 C…tn †=dt3 :

…6b†

Replacing an arbitrarily assigned time tn by universal notation t, the above equations can be rewritten as dP…t†=dt ˆ ……t† ÿ †=P…t† ‡ lC…t† ÿ
t=2d2 P…t†=dt2 ÿ
t2 =6d3 P…t†=dt3

…7a†

dC…t†=dt ˆ =P…t† ÿ lC…t† ÿ
t=2d2 C…t†=dt2 ÿ
t2 =6d3 C…t†=dt3 :

…7b†

The discretizing errors are embodied by the last two terms included in these righthand sides. 3.2. Expression for the implicit method In this case, the Taylor-series expansions around an advanced mesh point are made as Pn 7Pn‡1 ÿ
tdP…tn‡1 †=dt ‡
t2 =2d2 P…tn‡1 †=dt2 ÿ
t3 =6d3 P…tn‡1 †=dt3

…8a†

Cn 7Cn‡1 ÿ
tdC…tn‡1 †=dt ‡
t2 =2d2 C…tn‡1 †=dt2 ÿ
t3 =6d3 C…tn‡1 †=dt3 :

…8b†

Substituting the above equations into Eqs. (3a) and (3b), and making a similar replacement of notation as done in the above subsection, we obtain dP…t†=dt ˆ ……t† ÿ †=P…t† ‡ lC…t† ‡
t=2d2 P…t†=dt2 ÿ
t2 =6d3 P…t†=dt3

…9a†

dC…t†=dt ˆ =P…t† ÿ lC…t† ‡
t=2d2 C…t†=dt2 ÿ
t2 =6d3 C…t†=dt3 :

…9b†

The last two terms of the above equations embody a discretizing error induced by the implicit method. Compared with Eqs. (7a) and (7b), these equations have only a di€erence in signs of these third terms. 3.3. Expression for the Crank±Nicolson method

K. Hashimoto et al. / Annals of Nuclear Energy 27 (2000) 791±803

795

The Taylor-series expansions around a central time tn‡1=2 of adjacent meshes tn and tn‡1 are made as Pn‡1 7Pn‡1=2 ‡
t=2

Pn 7Pn‡1=2 ÿ
t=2

dP…tn‡1=2 † d2 P…tn‡1=2 † d3 P…tn‡1=2 † ‡
t2 =8 ÿ
t3 =48 2 dt dt dt3

Cn‡1 7Cn‡1=2 ‡
t=2

Cn 7Cn‡1=2 ÿ
t=2

dC…tn‡1=2 † d2 C…tn‡1=2 † d3 C…tn‡1=2 † ‡
t2 =8 ‡
t3 =48 2 dt dt dt3

…10a†

…10b†

…10c†

dC…tn‡1=2 † d2 C…tn‡1=2 † d3 C…tn‡1=2 † 3 ‡
t2 =8 ÿ
t =48 dt dt2 dt3

…10d†

d…tn‡1=2 † d2 …tn‡1=2 † d3 …tn‡1=2 † 3 ‡
t2 =8 ‡
t =48 dt dt2 dt3

…10e†

n‡1 7n‡1=2 ‡
t=2

n 7n‡1=2 ÿ
t=2

dP…tn‡1=2 † d2 P…tn‡1=2 † d3 P…tn‡1=2 † ‡
t2 =8 ‡
t3 =48 2 dt dt dt3

d…tn‡1=2 † d2 …tn‡1=2 † d3 …tn‡1=2 † 3 ‡
t2 =8 ÿ
t =48 : dt dt2 dt3

…10f†

Substituting the above equations into Eqs. (4a) and (4b), replacing a time tn‡1=2 by notation t and neglecting higher-order error terms than the fourth order, we obtain dP…t†=dt ˆ ……t† ÿ †=P…t† ‡ lC…t† ÿ
t2 =24d3 P…t†=dt3 ÿ =
t2 =8d2 P…t†=dt2 ‡ l
t2 =8d2 C…t†=dt2 ‡
t2 =8‰P…t†;  d2 …t†=dt2 ‡ …t†d2 P…t†=dt2 ‡ 2d…t†=dtdP…t†=dtŠ=

…11a†

dC…t†=dt ˆ =P…t† ÿ lC…t† ÿ
t2 =24d3 C…t†=dt3 ‡ =
t2 =8d2 P…t†=dt2 ÿ l
t2 =8d2 C…t†=dt2 :

…11b†

The above equations have the second-order errors on spacing, while the discretizing errors induced by the explicit or implicit method are of ®rst order, as shown in the above subsections. This feature indicates that the precision of the Crank±Nicolson method should be higher than those of the other methods.

796

K. Hashimoto et al. / Annals of Nuclear Energy 27 (2000) 791±803

4. Transfer function with discretizing error In this section, the zero-power (neutronic) transfer functions with discretizing errors are derived from the above one-point equations and numerical stabilities of the discretized equation systems are investigated from Bode diagrams of the transfer functions. 4.1. Zero-power transfer function Power and precursors concentration are separated into stationary and variational components as follows. P…t† ˆ P0 ‡ p…t†

…12a†

C…t† ˆ C0 ‡ c…t†

…12b†

Substituting the above equations into Eqs. (7a) and (7b), neglecting a nonlinear term and Laplace-transforming the result, the zero-power transfer function for the explicit method can be obtained as Z…s†  p~…s†=‰~…s†= Š=P0 ˆ =Y…s†

…13†

where Y…s† ˆ s ‡ s2 
t=2 ‡ s3 
t2 =6 ‡ ‰s ‡ s2
t=2 ‡ s3
t2 =6Š= ‰s ‡ l ‡ s2
t=2 ‡ s3
t2 =6Š:

…14a†

Substituting Eqs. (12a) and (12b) into Eqs. (9a) and (9b) and making the same operations as above, we obtain the transfer function for the implicit method. The transfer function can be described by substituting the following function into Eq. (13) Y…s† ˆ s ÿ s2 
t=2 ‡ s3 
t2 =6 ‡ ‰s ÿ s2
t=2 ‡ s3
t2 =6Š= ‰s ‡ l ÿ s2
t=2 ‡ s3
t2 =6Š

…14b†

The substitution into Eqs. (11a) and (11b) gives the transfer function for the Crank±Nicolson method, where Y…s† ˆ ‰s ‡ s3 
t2 =24Š=‰1 ‡ s2
t2 =8Š ‡ ‰s ‡ s3
t2 =24Š= ‰s ‡ l ‡ ls2
t2 =8 ‡ s3
t2 =24Š

:

…14c†

As mesh spacing
t approaches zero, the limits of the above transfer functions provide the conventional form of zero-power transfer function.

K. Hashimoto et al. / Annals of Nuclear Energy 27 (2000) 791±803

797

4.2. Bode diagrams Here the Bode diagrams of the above transfer functions with discretizing error are shown to investigate the numerical stability of each di€erent equation system. The discussion should be con®ned to the lower frequency range rather than the Nyquist frequency. The frequency is de®ned as fN ˆ 1=…2
t†

…15†

which speci®es the resolving frequency of a numerical calculation. A perturbation ranging in a frequency higher than the frequency can never be introduced into a discretized equation system. In what follows, the e€ective delayed-neutron fraction of 0.0075, the decay constant of 0.08 sÿ1 and the prompt-neutron generation time of 0.0001 s have been employed. Fig. 1 shows the Bode diagrams of the transfer functions based on various timeintegration methods, where mesh spacing of 0.01 s is speci®ed. Included in this ®gure for comparison is the exact curve, to which the above transfer functions can be reduced in in®nitely ®ne mesh. Compared with the exact curve, the gain of the transfer function for the explicit method is enhanced in the higher frequency range. The enhancement is nevertheless not so signi®cant as to unstabilize a discretized system. In contrast, the gain for the implicit method is diminished. The diminution leads to a stabilization of the system. The gain for the Crank±Nicolson method is distinct from those for the above two methods. The curve has a valley in the vicinity of the Nyquist frequency, which mathematically results from a zero of the transfer function, that is, the pole of the ®rst term of Eq. (14c). Consider a mechanistic reason for the valley. When a reactivity oscillator has the Nyquist frequency, a perturbation of alternating positive and negative reactivity on the time mesh should be introduced into a discretized system. Since the Crank±Nicolson method averages reactivities at adjacent mesh points, these positive and negative values are cancelled out. Consequently the gradient of unknown vanishes and it leads to no response. As observed in Fig. 1, the phase of the transfer function for the explicit method is delayed in the higher frequency range. On the other hand, the phase for the implicit method is in front of the exact curve. This is because the explicit and the implicit methods evaluate the gradient of unknown at backward and forward mesh points, respectively. The phase for the Crank±Nicolson method has a zero in the vicinity of the Nyquist frequency. This crossing results from a zero of the transfer function, as mentioned above. Fig. 2 shows the Bode diagram of the transfer function for a time spacing of 0.1 s. First, we can observe a broad peak in the gain curve for the explicit method. This peak suggests that numerical instability must force all unknowns to oscillate at the Nyquist frequency and hence no solution can be obtained. In the explicit method, such a coarse mesh should never be employed. The phase is distinct from that in Fig. 1. As a result of signi®cant delay, the phase is apparently advanced.

798

K. Hashimoto et al. / Annals of Nuclear Energy 27 (2000) 791±803

For the implicit method, the gain is diminished in the higher frequency range, which highly stabilizes the di€erence-equation system. This method has no limitation of time spacings for numerical stability. In the vicinity of the Nyquist frequency, the gain curve for the Crank±Nicolson method has a similar valley as observed in Fig. 1. This also originates from a zero of the transfer function. In the further higher frequency region, a sharp peak can be observed. This peak mathematically results from a pole of the transfer function, that is, the zero of Eq. (14c). Since the peak appears at a higher frequency than the Nyquist, numerical instability

Fig. 1. Bode diagram of transfer function for time spacing of 0.01 s.

K. Hashimoto et al. / Annals of Nuclear Energy 27 (2000) 791±803

799

based on the peak cannot be observed. The phase curve has two zeroes at the frequencies corresponding to the zero and the pole of the transfer function. 5. Stability limit for the explicit method As con®rmed in the above section, the application of the explicit method to a coarse mesh system may cause numerical instability. In this section, ®rst, the meshspacing limit for the stability is formulated. Next, several numerical calculations are performed to discuss the validity of this formula.

Fig. 2. Bode diagram of transfer function for time spacing of 0.1 s.

800

K. Hashimoto et al. / Annals of Nuclear Energy 27 (2000) 791±803

5.1. Expression for stability limit Reactor power and delayed-neutron precursors concentration are described as P…t† ˆ P0 et and C…t† ˆ C0 et ;

…16†

respectively. In general, the constant  and the two coecients included in the above equations are complex. The real and the imaginary parts of the constant represent the ampli®cation factor and frequency of oscillation, respectively. The substitution of Eq. (16) into Eqs. (7a) and (7b) gives the following eigenvalue problem: P0 ˆ ÿ =P0 ‡ lC0 ÿ 2
t=2P0 ÿ 3
t2 =6P0

…17a†

C0 ˆ =P0 ÿ lC0 ÿ 2
t=2C0 ÿ 3
t2 =6C0

…17b†

where reactivity term is omitted. We solve the above simultaneous equations to obtain a complex eigenvalue . Eliminating the coecients P0 and C0 from Eqs. (17a) and (17b), we obtain …
t2 =62 ‡
t=2 ‡ 1†…
t2 =63 ‡
t=22 ‡  ‡ l ‡ =† ˆ 0:

…18†

The quadratic equation
t2 =62 ‡
t=2 ‡ 1 ˆ 0

…19†

has conjugated complex roots, but the real part of the roots is invariably negative and this mode decays out rapidly. On the other hand, the possibility of positive real part is present in the cubic equation
t2 =63 ‡
t=22 ‡  ‡ l ‡ = ˆ 0:

…20†

Fig. 3 shows the mesh-spacing dependence of the conjugated complex roots. As the mesh is made coarser, the real part (i.e. ampli®cation factor) increases but the imaginary part (i.e. oscillation frequency) decreases. The sign of the factor varies from negative to positive at a mesh width, which is an upper limit of mesh spacings for numerical stability. From Eq. (20) the limit can be obtained as
tL ˆ 3=…l ‡ =†:

…21†

The time-mesh spacings should be determined according to the following criterion:
t <
tL

…22†

K. Hashimoto et al. / Annals of Nuclear Energy 27 (2000) 791±803

801

Fig. 3. Mesh-spacing dependence of temporal eigenvalue  for explicit integration method.

This criterion is based on the assumption that higher-order errors than fourth and reactivity perturbation term are neglected. This assumption may result in some overestimation of the spacing limit. It is preferred to a ®ner mesh than the estimation. On the other hand, the present description based on one delayed-group theory leads to slight error in the estimation. As can be seen from Fig. 2, numerical instability is generated in the particular high frequency region where the in¯uence of the delayed neutron is negligible. Therefore, a detailed description of delayed neutron term is of no signi®cance. 5.2. Numerical calculation Numerical calculations based on the explicit method are illustrated to examine the validity of the present criterion, Eq. (22). An indicial response of power to stepshaped reactivity insertion of +20¢ was calculated for several mesh spacings. This problem is available for the stability test because such a discrete disturbance highly excites all inherent modes. This result is shown in Fig. 4, where we employ the same kinetic parameters as in Section 4. When the mesh spacing of 0.01 s is speci®ed, the numerical solution is quite stable and agrees very closely with the exact solution. When the coarser meshes are speci®ed, apparent oscillations of solution can be observed. In the mesh spacing of 0.03 s, the perturbation forces the solution to oscillate instantly but the oscillatory behavior decays rapidly. This case can be classi®ed into a category of stability. In the spacing of 0.04 s, however, the amplitude of the oscillation continues to increase with time. From these results, the stability limit is expected to exist between two spacings of 0.03 and 0.04 s. Eq. (21) estimates the limit of 0.039 s, which is consistent with the above expectation.

802

K. Hashimoto et al. / Annals of Nuclear Energy 27 (2000) 791±803

Fig. 4. Numerical results of indicial response to reactivity input of 20¢.

6. Conclusions We derived the one-point kinetic equations by considering the numerical errors induced by several time-integration methods. The zero-power transfer functions based on the present equations were demonstrated to study the numerical stability of the discretized systems. The applications of the implicit and the Crank±Nicolson methods give an undertaking to make numerically stable while the explicit method makes the numerical system unstable in the mesh-spacing range broader than a limit. The upper limit of spacings for numerical stability is formulated.

K. Hashimoto et al. / Annals of Nuclear Energy 27 (2000) 791±803

803

This criterion is based on some assumptions but is nevertheless very useful as a measure of mesh spacings for numerical stability. In the near future, we will present several criteria for the stability of space-dependent di€usion and transport equations. Furthermore, the stability of a nonlinear system is a subject of great interest. References Hirt, C.W., 1968. Heuristic stability theory for ®nite-di€erence equations. J. of Computational Physics 2, 339. Roache, P.J., 1976. Computational Fluid Dynamics. Hermosa Publishers Inc.