Computer simulations have shown that stability calculations in boiling water reactors (BWRs) are very sensitive to a number of input parameters, modeling assumptions, and numerical integration techniques. Following the 1988 LaSalle instability event, a significant industry-wide effort was invested in identifying these sensitivities. One major conclusion from these studies was that existing time-domain codes could best predict BWR stability by using explicit methods for the energy equation with a Courant number as close to unity as possible. This paper presents a series of sensitivity studies using simplified models, which allow us to determine the effect that different numerical integration techniques have on the results of stability calculations. The present study appears to indicate that, even though using explicit integration with a Courant number of one is adequate for existing codes using time-integration steps of less than 10 ms, second-order solution techniques for the time integration can result in significant improvements in the accuracy of linear (i.e., decay ratio) stability calculations

March-Leuba, J.

Peiro, D. G.

UNT Digital Library

ORNL/CP -94956 _- SENSITIVITY OF BWR STABILITY CALCULATIONS TO NUMERICAL INTEGRATION TECHNIQUES Damian Ginestar Peiro Departamento de Matematica Aplicada; Universidad Politecnica de Valencia Camino de Vera. 14; 46071 Valencia, SPAIN C JosC March-Leuba Instrumentation & Controls Division; Oak Ridge National Laboratory’ P.O. Box 2008: Oak Ridge, TN 37831-6010; U.S.A. ABSTRACT Computer simulations have shown that stability calculations in boiling water reactors (BWRs) are very sensitive to a number of input parameters, modeling assumptions, and numerical integration techniques. Following the 1988 LaSalle instability event, a significant industry-wide effort was invested in identifying these sensitivities. One major conclusion from these studies was that existing time-domain codes could best predict BWR stability by using explicit methods for the energy equation with a Courant number as close to unity as possible. This paper presents a series of sensitivity studies using simplified models, which allow us to determine the effect that different numerical integration techniques have on the results of stability calculations. The present study appears to indicate that, even though using explicit integration with a Courant number of one is adequate for existing codes using time-integration steps of less than 10 ms, second-order solution techniques for the time integration can result in significant improvements in the accuracy of linear (i.e., decay ratio) stability calculations. The problems associated with the numerical solution of oscillatory-type systems can be illustrated by solving the following ordinary differential equation (ODE) For this particular simulation, we have used cs = 0.1 and w = z, which results in an oscillation of frequency 0.5 Hz and DR = 0.82. Following the conventional numerical solution approach, we can convert Eq ( 1 )  into two first order ODE’S and solve them numerically using different time-integration aIgorithms and time step. The results are shown in Fig. 1 and indicate that first order integration methods (Euler’s method) either overestimate (explicit) or underestimate (implicit) the decay ratio significantly even for reasonable time steps. Second order time-integration methods are very accurate even for large integration time steps. All methods converge to the analytic solution as the time step is reduced (see Fig. 3). Managed by Lockheed Martin Energy Research Corp. For the U.S. Department of Energy under contract DE-AC05-960R22464 DISCLAIMER This report was prepared as an account of work spomored by an agency of the United States Government. Neither the United States Government nor any agency thereof, nor any of their employees, make any warranty, express or implied, or assumes any legal liabii- ty or responsibility for the accuracy, completeness, or usefulness of any information, appa- ratus, pmduct, or process disclosed, or represents that its use would not infringe privately owned rights. Reference herein to any specific commercial product, process, or service by trade name, trademark, manufacturer, or dhenvise does not necessarily constitute or imply its endorsement, recommendation, or favoring by the United States Government or any agency thereof. The views and opinions of authors expressed herein do not necessar- ily state or reflect those of the United States Government or any agency thereof. 2 By analytically studying the propagation of errors in the numerical solution of Equation (l), it can be shown that the error in the numerically calculated decay ratio, DRcalc, is DRLalc = DR e’ RWAr I Explicit Runge Kutta Second Order , aac t  where the plus sign applies to explicit methods the minus sign applies to implicit methods. For small time steps, At (in units of seconds), can be $ approximated by n I I 0.5 DRcalc = DRexacr ( 1 10 At ) implicit Euler i (3) Thus, the percent decay ratio error for first-order integration methods is approximately equal to the integration time step in milliseconds. Most time- 0 1  I 0.001 0.01 0 .  1 INTEGRATION TIME STEP (s) domain system codes use integration times of at most 1 to 10 ms; thus, the decay ratio error caused by first-order integration errors should be less than ~i~~~ 1, Decay 10% for these codes. Higher-order time integration numerical integration on schemes (e.g., explicit Runge Kutta of second order) result in negligible errors for integration steps of the order of 100 ms or less. Stability calculations are unlikely to use time steps larger than 100 ms due to Courant number constraints; thus, using time- integration schemes of order higher than second for linear stability (Le. decay ratio) calculations is not likely to result in significant improvements. errors induced by first-order (1) Even though the partial differential equations (PDEs) that represent BWR stability are significantly more complex than Eq. ( l ) ,  the qualitative behavior of the solution is similar and one should expect similar behavior. A more complex model of BWR dynamics can be represented by the following equations where Eq (4) models the convection of the density wave, a, which propagates at constant velocity, V, in this simplified model. The right hand side of Eq ( 5 )  models the reactivity feedback, which is proportional to the core-average density, and an approximation to the power, Q, dynamics. The constant K controls the stability of the system and can be thought of as being proportional to the product of the steady-state power times the density reactivity coefficient. For the analysis presented in this paper, the value of the constants are V = 2 ds, H = 4 m, and K = 3, which result in a decay ratio of 0.41. The boundary conditions are a(z=O, t) = 0; a(z,  r = O )  = 0; Q ( t = O )  = 0 This model, although simple, captures most of the numerical problems associated with solving the linear BWR stability problem in time and space (i.e., decay ratio of the reactor core-wide mode of oscillation). The non-linear response (e.g. large-amplitude limit cycles) is beyond the scope of this study. The space dependence of this model can be discretized by defining node-averaged void fractions ‘--I Integrating Equations (4) and ( 5 )  for each node, we obtain the following nodal equations. dai a(zJ  -a(zj-l> - + v  = Q A z  dt Az i = N  dt i = l  - -  dQ - -K ai We must note that the above equations are exact. The discretization errors occur when we must represent the node boundaries values, a(z,), as a function of the node-average values, ai. This is accomplished by integrating our approximation for the derivative over each of the nodes. For example, for a first-order up-wind differentiation, we approximate the spatial derivative as and we integrate Eq (9) over each node to obtain an approximation for the convective term in Eq (8) a, - ai-1 (10) N Z i )  -a(z,-J Az Az Z,-I V I  Az Thus, to first order in space, our system of partial differential equations becomes da, a, - a,-1 - + v  = Q A z  dt Az Using a similar technique, we can develop a second order up-wind approximation for our equations. The second-order up-wind approximation for the spacial derivative is and the second-order-approximation equations become da, 3a, -4aj-,  +ai-? - + v  = Q A z  dt 2Az i = N  dQ - = -K a; dt i = l  These discretized models. represented by Eqs ( 1  1 )  and (13), allow us to perform a number of analytical 0.60 0.50 0 5 Y a > 0.40 4 n 0.30 Runge Kutta m m --e Euler. 12 nodes -++ RK-2, 12 nodes --E. Euler, 24 nodes -0- RK-2, 24 nodes Euler, 48 nodes +RK-2, 48 nodes ni Euler, 96 nodes -9-RK-2, 96 nodes 0.20 I 0.00 0.02 0.04 0.06 0.08 0.10 INTEGRATION TIME STEP (s) Figure 2. Calculated decay ratios error for Eqs (4) and (5) for different time integration schemes using the first-order explicit space nodalization described by Eq (11). and numerical simulations to evaluate the error introduced by the spatial discretization and by the time integration. For these analyses, we have assumed that the lower plenum has a constant a=O, so negative indexes in Eq (13) are defined. For example, by integrating the first-order spatial approximation of Eq (1 1) using different number of nodes and first and second-order time-integration algorithms, we obtain the results shown in Fig. 2. The numerical-integration results shown in Fig. 2 indicate that space-nodalization errors are not as important as time-integration errors. For example, while the relative error between the 24-node and 96 node case (at constant integration time step) is less than 5%, increasing the time integration step increases the error exponentially to values greater than 50%. It appears from these results that the fact that better results are obtained with larger number of nodes is not because of spatial nodalization errors, but because it forces the integration algorithm to take smaller steps to maintain a Courant limit lower than 1. Figure 2 indicates that using a higher-order time-integration algorithm reduces the decay ratio error significantly, and the error reduction is independent of nodalization scheme. These results confirm our earlier conclusion that decay ratio errors are controlled by errors in the time integration, and not by spatial discretization errors. Most stability calculations with time-domain system codes are performed with at least 24 axial nodes and the average Courant number in the core is from 0.2 to 0.5, because the controlling Courant number occurs at the core exit where the velocity is highest. The integration time step is smaller than a Courant number of one would require because it is limited by either the smallest node or the node with the fastest velocity. Based on this simulation, we conclude that decay ratio errors induced by first-order explicit integration inaccuracies are expected to be in the order of lo%, which is adequate for stability calculations (all NRC-licensed stability codes are assumed accurate to within 20%). ~ . . . . . .. .. .. .. .. .. .... .. .. ....... . .. ..... .. .. . .... .. .. .. .. . . . ........ .... .. . ../ + ..... . 9 i +12 Nodes, 1s t  Order -8-24 Nodes, 1st  Order e 9 6  Nodes,  1st  Order +24 Nodes, 2nd Order I 1 -48 Nodes,  2nd Order +96 Nodes, 2nd Order *48 Nodes,  1 s t  Order -E3- 12  Nodes, 2nd Order ~ 0 0.02 0.04 0.06 0.08 INTEGERATION T IME S T E P  0.1 Figure 3. Decay ratio error dependence with time is similar for first- and second-order space discretization. Using the second-order spatial discretization scheme of Eq (13), we obtain results consistent with the one shown in Figure 2 for the first-order scheme (Eq (1 1)). These results are shown in Fig. 3 for the first- order explicit (i.e, Euler) time integration scheme. We observe that using a higher-order spatial nodalization does not change qualitatively the error dependence with time step integration, which appears to be completely independent of the space nodalization scheme. A counter-intuitive result from these simulations is that using a higher-order space nodalization scheme actually increases the error when a small number of nodes (e.g. 12) is used. These results are confirmed by an analytical study, where the errors induced by the spacial and time nodalization can be treated separately by calculating the dominant eigenvalues of Eqs ( 1  1 )  or (13). The results of this analysis are shown in Fig. 4. We can conclude form these results that higher-order spatial integration results in better convergence (i.e., lower error) only if the number of nodes is sufficiently large (more than 24 nodes). We also note in Figs. 3 and 4 that the space-nodalization errors are at least an order of magnitude smaller than the errors induced by the time integration if relatively large time steps are used. For example, even using a coarse 12-node mesh, the space nodalization error is of the order of lo%, while using a time step of 50 ms results in an error of approximately 50%. These numerical results can also be confirmed analytically by calculating the error induced in the dominant system poles by the time discretization. The results of this analytical study are shown in Figure 5.  For the results shown in this figure, the following time discretization was used which results in explicit first order for 8 = 0, implicit first order for 8 = 1 ,  and semi-implicit second order for 8 = 0.5 0 % C 0 .- c. E cn 0 S V m a v) 0 a 3 U - - a -5% * L 2 li -10% a 0 m % m 0 .- - 8 -15% I Space Integration 0 12 24 36 48 60 72 84 96 Number of Nodes  Figure 4. Decay ratio error due to space nodalization, analytical results. Higher-order methods converge faster, but errors can be larger if large nodes are used. For more than 24 nodes, errors are -5% U K w 2 K >. a n Y 75 yo 50% 25% 0 Yo -25% -50% -f+- Second-Order Semi-Implicit 0 0.02 0.04 0.06 0.08 0.1 TIME INTEGRATION STEP (s) Figure 5. Decay ratio error due to time integration, analytical results. First-order integration errors can be very large. Higher-order time-integration results in negligible errors for time steps up to 100 ms. o.6 1 I 0.3 0.2 0 75- 12 Nodes,  1 st Order +48 Nodes ,  1 s t  Order  I - 1 2  Nodes ,  Semi-Implicit -3-24 Nodes,  Semi-Implicit 1 1  I +24 Nodes ,  1 s t  Order  - 4 9 6  Nodes,  1s t  Order  ___ 5 48 Nodes,  Semi-Implicit +96 Nodes,  Semi-Implicit 11 0.02 0.04 0.06 0.08 INTEGERATION TIME STEP 0.1 Figure 6. Using a semi-implicit scheme, cancellation of errors appears to increase the accuracy. This method is not recommended. Large system-type simulation codes rarely use a consistent integration scheme for all state variables. For example, these codes may use an explicit time integration scheme for the energy equation, an implicit scheme for the neutronics, and a semi-implicit scheme for the mass and momentum equations. Under these conditions, the numerical errors are essentially unpredictable, unless the time steps are so small that the errors become negligible. We have attempted to simulate this situation with our simple model by using an explicit scheme for the energy equation [Eq (4)] and an implicit scheme for Eq (5) .  The numerical implementation consists on simply stepping the a variables in time and then using these new values of a at time t to calculate the power, Q(t). The results of this simulation are shown in Fig. 6. Surprisingly, this unorthodox “mixed” scheme performs better than either fully explicit or fully implicit. The reason for this result is a lucky cancellation of errors. The explicit a integration produces a positive error, which is canceled in part by the negative error introduced by the implicit integration of Q. Although the results shown in Fig. 6 appear to indicate that a semi-implicit method should be the preferred course of action, this method is not recommended unless is thoroughly tested under a wide variety of circumstances. Relying of cancellation of errors to obtain accuracy is extremely dangerous because, under some untested circumstances, the magnitude of errors could change dramatically. CONCLUSIONS This paper documents the results of a series of numerical simulations related to the problem of linear boiling water reactor stability (Le.. decay ratio calculations). These simplified simulations indicate that the numerical techniques used to calculate the linear reactor stability introduce two types of errors: ( 1 )  time-discretization errors, and (2) space-discretization errors, which are independent of each other. The space-discretization errors are always non-conservative (the calculated decay ratio is lower than the actual), and the time-discretization errors could be conservative (for explicit methods) or non- conservative (for non-explicit methods). This simulations also indicate that the time-discretization errors can be significantly large for reasonable time steps (e.g., more than 50% error for time steps of 50 ms,) while space-discretization errors are typically small (e.g. less than 5% for 24 or more axial nodes.) These simulations indicate that, while second-order space discretization errors are more accurate for fine nodalizations, the errors can actually be larger than those of first-order space discretization when coarse nodalizations are used (e.g., 12 or less axial nodes.) For all cases tested, second-order time discretizations decreased the error to insignificant levels (e.g., using a 50 ms time step, the error is reduced from more than 50% with first order to - 1 % error if a second order time discretization scheme is used.) This simulations also indicate that, using mixed semi-implicit methods, one can achieve cancellation of errors that appear to increase the accuracy of the calculation. These methods are not recommended unless they are thoroughly and exhaustively tested, because this cancellation of errors may not be reliable under all circumstances. The preferred course of action is to eliminate the source of the error by using higher-order time-integration and small nodalization schemes. Based on our simplified simulations, we conclude that detailed numerical simulation studies must be performed to evaluate the particular numerical integration scheme used in each code. It is important that these tests include the complete closed-loop dynamics, as it is clear from these simulations that a feedback loop [e.g., Eq (5 ) ]  changes significantly the numerical response. 

Sensitivity of BWR stability calculations to numerical integration techniques

https://digital.library.unt.edu/ark:/67531/metadc694152/m2/1/high_res_d/548854.pdf

Sensitivity of BWR stability calculations to numerical integration techniques

Abstract

Similar works

Full text

Available Versions

UNT Digital Library