UFDC Home  myUFDC Home  Help 
Title Page  
Dedication  
Acknowledgement  
Table of Contents  
List of Tables  
List of Figures  
Notation  
Abstract  
I. Introduction  
II. Description of method  
III. Application of the hybrid...  
IV. Affects of the characteristics...  
V. Discussion  
Appendixes  
A. Analysis of iterative techn...  
B. Hybrid computer programs  
C. Exact solution to the approximate...  
References  
Biographical sketch 



Table of Contents  
Title Page
Page i Dedication Page ii Acknowledgement Page iii Table of Contents Page iv Page v List of Tables Page vi List of Figures Page vii Page viii Notation Page ix Page x Abstract Page xi Page xii Page xiii I. Introduction Page 1 Page 2 Page 3 Page 4 Page 5 Page 6 Page 7 Page 8 Page 9 Page 10 Page 11 Page 12 Page 13 Page 14 II. Description of method Page 15 Page 16 Page 17 Page 18 Page 19 Page 20 Page 21 Page 22 Page 23 Page 24 Page 25 Page 26 Page 27 Page 28 Page 29 Page 30 Page 31 Page 32 Page 33 Page 34 Page 35 Page 36 Page 37 Page 38 Page 39 Page 40 Page 41 Page 42 Page 43 Page 44 Page 45 Page 46 Page 47 Page 48 Page 49 Page 50 Page 51 Page 52 Page 53 Page 54 Page 55 III. Application of the hybrid method: Examples Page 56 Page 57 Page 58 Page 59 Page 60 Page 61 Page 62 Page 63 Page 64 Page 65 Page 66 Page 67 Page 68 Page 69 Page 70 Page 71 Page 72 Page 73 IV. Affects of the characteristics of nuclear problems upon the hybrid technique Page 74 Page 75 Page 76 Page 77 Page 78 Page 79 Page 80 Page 81 Page 82 Page 83 Page 84 Page 85 Page 86 V. Discussion Page 87 Page 88 Page 89 Page 90 Page 91 Page 92 Page 93 Page 94 Page 95 Page 96 Page 97 Appendixes Page 98 A. Analysis of iterative techniques Page 99 Page 100 Page 101 Page 102 Page 103 Page 104 Page 105 Page 106 Page 107 Page 108 Page 109 Page 110 Page 111 Page 112 Page 113 Page 114 Page 115 Page 116 Page 117 Page 118 Page 119 Page 120 Page 121 B. Hybrid computer programs Page 122 Page 123 Page 124 Page 125 Page 126 Page 127 Page 128 Page 129 Page 130 Page 131 Page 132 Page 133 Page 134 Page 135 Page 136 Page 137 Page 138 Page 139 Page 140 Page 141 Page 142 Page 143 Page 144 Page 145 Page 146 Page 147 Page 148 Page 149 C. Exact solution to the approximate equations Page 150 Page 151 Page 152 Page 153 Page 154 Page 155 Page 156 Page 157 Page 158 Page 159 Page 160 Page 161 Page 162 Page 163 Page 164 Page 165 References Page 166 Page 167 Page 168 Page 169 Page 170 Biographical sketch Page 171 Page 172 Page 173 

Full Text  
HYBRID COMPUTATION OF SPACE AND ENERGY DEPENDENT NUCLEAR REACTOR KINETICS By Kenneth R. Schultz A Dissertation Presented to the Graduate Council of the University of Florida in Partial Fulfillment of the Requirements for the Degree of Doctor of Philosophy UNIVERSITY OF FLORIDA 1971 To my Mother and Father ACKNOWLEDGMENTS The author would like to express his appreciation to the members of his supervisory committee. Dr. T. E. Bul lock, Dr. M. J. Ohanian, Dr. E. E. Carroll, Dr. R. Gaither and Dr. Luehr for their assistance during the course of this research. Dr. M. J. Ohanian and Dr. Thomas Bullock are due special thanks for the many hours they have devoted to directing this study. The entire faculty and staff of the University of Florida Nuclear Engineering Sciences Department were always helpful beyond the call of duty. Dr. A. J. Mockel was of special help in the analysis of the mathematical techniques. Use of the hybrid computer was provided by the De partment of Nuclear Engineering Sciences at the University of Florida. The author's studies at the University of Florida were supported by United States Atomic Energy Commission Traineeships and Fellowships and a Graduate Assistantship under National Science Foundation Grant No. GK2106. Thanks for this support are gratefully given. Above all the author would like to thank his wonderful wife, Mary Lou, for her patience and support during his six years of graduate study. Without her en couragement and help this research could not have been done. iii TABLE OF CONTENTS ACKNOWLEDGMENTS ..... ................................ LIST OF TABLES ....... .............................. LIST OF FIGURES ................................ ... NOTATION .......................................... ABSTRACT ........................................... Chapter I. INTRODUCTION ............................... Space and Energy Dependent Reactor Kinetics .............................. Survey of Kinetics Methods .............. The Use of Hybrid Computers ............. II. DESCRIPTION OF METHOD ...................... Mathematical Model ...................... Description of Hybrid Technique ......... Selection of Iterative Scheme .......... Description of Hybrid System ............ Implementation .......................... III. APPLICATION OF THE HYBRID METHOD: EXAMPLES ................................ Example 1: A Onegroup Model ........... Example 2: A Twogroup Model ........... Example 3: A Twogroup Model with Temperature Feedback .................. IV. AFFECTS OF THE CHARACTERISTICS OF NUCLEAR PROBLEMS UPON THE HYBRID TECHNIQUE ..... ....................... Playback .................... .......... Accuracy ................................ Run Time ................................ Page iii vi vii ix xi 1 1 3 10 15 15 22 31 48 50 56 56 60 65 74 74 77 79 TABLE OF CONTENTSContinued Page DISCUSSION ................................. Recommendations for Further Study ....... Conclusions ............................. Appendix A. ANALYSIS OF ITERATIVE TECHNIQUES .......... B. HYBRID KINETICS PROGRAMS .................. C. EXACT SOLUTION TO THE APPROXIMATE EQUATIONS .............................. REFERENCES ........................................ BIOGRAPHICAL SKETCH ............................... Chapter V. 99 122 150 166 171 LIST OF TABLES Table Page 1. Onegroup Turkey Point Example ............. 60 2. Twogroup Yasinsky and Henry Example ................................. 64 3. Two Groups with Feedback J. B. Yasinsky Example ........................ 73 4. Example 3: Reactor Data ................... 146 5. Example 3: Initial Conditions ............. 148 6. Example 3: Pot List ....................... 149 LIST OF FIGURES Figure Page 1. Hybrid computer ............................. 12 2. Hybrid diagram ............................. 24 3. Iterative procedure ....................... 25 4. Hybrid flow diagram ......................... 28 5. Example of nonconvergence ................... 43 6. Overrelaxation results .................... 46 7. Analog implementation ...................... 53 8. Onegroup exampleFlux vs. distance ........ 57 9. Onegroup exampleFlux vs. time ............ 58 10. Twogroup exampleFlux vs. distance ........ 61 11. Twogroup exampleFlux vs. time ............ 62 12. Example threeThermal shape function ....... 69 13. Example threeReactor power vs. time ....... 70 14. Example threeFuel temperature vs. time .... 71 15. Data playback ............................... 76 16. Sensitivity to offset error ................. 78 17. Computation time vs. number of space points .................................... 83 18. Computation time vs. integration interval ................................. 86 19. HYKIN 3 analog neutronics .................. 126 vii LIST OF FIGURESContinued Figure Page 20. HYKIN 3 analog thermalhydraulics .......... 127 21. Eigenvalue error ........................... 158 viii NOTATION Standard reactor physics notation as defined below is used throughout this dissertation. Neutroni C S ex Z (E 'E s X(E cs = Neutron flux = Delayed neutron precursor density in the Zth group = External neutron source = Total neutron cross section Neutron cross section for absorption = Neutron cross section for absorption a = Neutron cross section for fission ) = Neutron cross section for scattering from group E' to E D = Neutron diffusion coefficient v = Neutron fission yield ;) = Energy distribution of neutrons produced by fission 8 = Total delayed neutron fraction S= Delayed neutron fraction to th precursor group = Delayed neutron precursor decay constant Thermal Hydraulics V = Reactor coolant volume c Volumetric specific heat of coolant Tc T = Coolant temperature AH = Heat transfer area per unit coolant volume U = Conductivity of cladding h = Heat transfer coefficient W = Coolant flow rate T = Fuel temperature C = Specific heat of coolant c T. = Coolant inlet temperature in r = Fraction of fission power appearing in coolant q = Fission power V = Reactor fuel volume C = Conversion factor; fissions to energy Pf = Fuel density Cf = Specific heat of fuel Abstract of a Dissertation Presented to the Graduate Council of the University of Florida in Partial Fulfillment of the Requirements for the Degree of Doctor of Philosophy HYBRID COMPUTATION OF SPACE AND ENERGY DEPENDENT NUCLEAR REACTOR KINETICS By Kenneth R. Schultz March, 1971 Chairman: Dr. M. J. Ohanian CoChairman: Dr. Thomas E. Bullock Major Department: Nuclear Engineering Sciences The neutron flux in a nuclear reactor may generally be described as a function of time, space, and energy. Substantial errors can result if spatial effects are ignored when studying the kinetics of large reactors or if spectral effects are ignored in the study of fast reactor dynamics. Hence, there is substantial current interest in the study of space and energy dependent nuclear reactor kinetics. The digital computer computations used in these studies can be lengthy and hence expensive. An analog com puter, which solves differential equations efficiently, might be used in kinetics studies but solution of realistic time, space, and energy dependent reactor problems is beyond the ability of an analog machine alone. By combining a digital computer with an analog computer, a hybrid machine is obtained that can possess the best features of both digital and analog machines. A method of direct numerical solution of the time, space, and energy dependent neutron diffusion equation has been developed that utilizes the unique capabilities of the hybrid computer. The hybrid method uses an iterative discrete space continuous time algorithm. The space derivatives of the conventional multigroup reactor kinetics equations are represented by finite difference approximations. The re sulting very large set of coupled ordinary differential equations in time can be solved quite readily on the hybrid computer by multiplexing the analog hardware and using an iterative procedure. The flux, precursor, and temperature equations at one space point are set up on the analog com puter. The digital computer supplies the space coupling terms from adjacent space points; it also observes and stores the solutions for later playback. Feedback effects are calculated by the digital machine and are implemented by means of variable coefficients. Automatic digital rescaling of the analog equations allows consideration of realistic reactor dynamics studies involving flux level changes of many orders of magnitude. The hybrid method is well suited for solution of the space and energy dependent reactor kinetics equations and shows strong potential for use in nuclear reactor xii engineering design. It offers reasonable accuracy at modest cost. Hybrid solutions to typical spatial kinetics examples and a description of the HYKIN 3 hybrid kinetics code, which is functionally equivalent to the WIGL 2 digital kinetics program, are presented. xiii CHAPTER I INTRODUCTION During the 28.years since Enrico Fermi's group of physicists brought the first nuclear reactor to critical ity on December 2, 1942, nuclear engineering has come a long way. Reactor power levels have been raised from the original 500 milliwatts of Fermi's Chicago Pile 1 to thousands of megawatts. Temperatures in the more exotic reactors approach 40000 F. This year over 60 per cent of all the orders placed for electric generating stations will be for nuclear power plants. In 28 years nuclear reactors have evolved from primitive physics experiments into sophisticated pieces of production machinery. Nuclear engineering has become an established professional field. Advances in the field of nuclear engineering have brought about the need for more detailed calculational models. Space and Energy Dependent Reactor Kinetics The neutron flux in a reactor may generally be adequately described as a function of time, space and energy. The first reactor kinetics studies, however, made use of the point model equations which ignore the  1   2  space and energy dependence of the flux. In the design of large reactors with high power densities, space de pendent effects that had not been of previous concern began to be noticed. It has been shown[]* that the reactivitiy worth of a local perturbation can be much greater than would be calculated by a space independent analysis. It is necessary, therefore, to consider space dependent effects in the analysis of nuclear reactor kinetics. The need to study the energy variable in detail is peculiar to fast reactors. In these reactors the system multiplication is a strong function of the neutron energy spectrum. The predominant feedback effects, the coolant density effect and the Doppler effect, are also very spectrum dependent. It is, therefore, important to have knowledge of the neutron energy spectrum during a transient if the dynamic behavior of a fast reactor is to be calculated accurately. Thus it can be seen that a definite need exists for a method of nuclear reactor analysis that can consider time, space and energy depend ence and that includes realistic feedback representation. Superior numbers in square brackets refer to the numbered references in the Bibliography.  3  Survey of Reactor Kinetics Methods Before embarking on any program to develop a new method of reactor kinetics analysis, we must first examine what has already been done. The basic equations that describe the behavior of a population of neutrons is the Boltzmann neutron transport equation. The full seven dimensions of the Boltzmann phase space are too much to cope with for practical kinetic cal culations. In a large reactor the flux is essentially isotropic. By ignoring the angular dependence of the neutron flux, the time, space and energy dependent diffusion model is obtained. 3] Neutron Flux Equation STET (r,E,t) = V D(r,E,t)VQ(r,E,t) Zt(r,E,t) x @(r,E,t) + J dE' Z (r,E' + Et) (r,E',t) N + x(E)X C (r,t) + x(E) dE' (1 (E') x vCf(r,E',t)D(rE',t) + Sext(r,E,t) Delayed Neutron Precursor Equation 8C (r,t) t dE'vZ (r,E',t)O(r,E',t) X C (r,t)  4  These are the basic equations from which most reactor kinetics analyses begin. The approximation resul ting from completely ignoring the angular dependence of the neutron flux may seem crude, but it has been shown[2] that the results of reactor kinetics studies using the dif fusion approximation compare quite favorably with the results of much more rigorous treatments. The treatment of this equation may be divided into two general categories: Point model studies.In these studies the time dependence of the neutron flux is assumed to be separable from the other independent variables, space and energy. 4(f,E,t) = f(l,E) T(t) This may be made rigorous by integrating over space and energy rather than simply assuming separability.[48] The rigorous integral parameters of the equations however are then functions of time. This means that the "rigorous point model" equations are still as much work to solve as the time, space and energy dependent diffusion equation itself. Space and energy dependent studies.These studies attempt to account for the fact that the reactor kinetics equation is nonseparable. A great deal of work has been done using the point model equations. This includes analytic studies and  5  numerical methods using digital, analog, and hybrid computers. Extensive work on analytic stability analyses has been done using these equations. Despite the attention the point model equations have enjoyed due to their com parative simplicity, the lack of ability to consider spatial and spectral effects means that point model studies, by themselves, are not sufficient for reactor kinetics analyses. Those studies that do not assume that the kinetics equations are separable but rather attempt to account for the interaction of space and energy dependence as well, may further be subdivided and classified into four categories: 1. Analytic Solutions 2. Numerical Treatments 3. Series Expansion 4. Factorizing Methods Analytic solutions that include time, space and energy dependence are possible for special cases. These cases include the small variation, linearized problems that can be handled by means of Laplace and Fourier transforms of the time and space variables. These lead to space and energy dependent transfer functions.[12] The limitation to small perturbation linear problems however restricts the usefulness of these analytic transformation methods. The reactor kinetics equations are analytically unsolvable when realistic feedback effects are considered. Even with  6  no feedback and very simple, uniform system changes, the analytic study includes expansion in infiniteseries to handle the space and energy dependence and involves a great deal of computation. Analytic methods of stability analysis have been developed that consider the complete set of space and energy dependent partial differential equations that de scribes the flux.[13,14] These are purely stability studies and do not calculate the flux behavior. Complete numerical treatments have been developed such as the WIGLE and FREAK codes. 11 522] The only approxi mations inherent in these approaches are those due to dis cretization of variables and use of finite differences for the derivative terms. The drawback of these methods is cost. Since the required computations are both involved and numerous, they are time consuming. On the large digital computers necessary to handle such a problem, time and money go hand in hand. Since computation time rises at least linearly and possibly with the square of the number of energy groups used, detailed spectral dependence is ex pensive to include and requires still larger computers. There is no theoretical reason why the full seven dimensional Boltzmann equation, complete with delayed neutron equations and all feedback effects, cannot be solved numerically. From a practical point of view, however, this is impossible anduannecessary. The greatest obstacle to  7  direct numerical computation is the small size required for the time step. If the equations are linear, an implicit time differencing scheme may be used. If the equation co efficients change with time, the time steps must be kept quite short to insure accuracy. If an explicit time dif ferencing scheme is used, computational stability restricts the maximum time step size to on the order of the prompt neutron lifetime. In a fast reactor this is on the order of 107 seconds. Because of this the large alldigital codes using direct numerical integration can be slow and costly. Following the advice of Fermi who said, "When in doubt, make a series expansion,"[21] much work has been done along the line of expansion of the flux in a series of time independent modes or nodes with time dependent coefficients 4(F,E,t) = ) mi(F,E) Pi(t) This technique can be cumbersome, particularly in the case of reactors with feedback which causes the flux equation coefficients to vary with time and may cause the modes used at the beginning of the problem to be completely unsuited for the system at a later time.[1,243s] Error is always introduced when the series is truncated after some finite number of terms.  8  The factorizing approach[3l637] separates most of the time dependence out of the kinetics equation: <(r,E,t) = ((r,E,t) (t) where 4 is a weak function of time. This method is further subclassified by the approximations now made. 1. Improved quasistatic method uses large time steps for 2t, small steps for . 2. Quasistatic method ignores . at 3. Adiabatic approximation ignores both and 3. at* 4. Linearized adiabatic method assumes separability and then uses the new value of ( to calculate a new shape function 4 every so often, which in turn is used to continue calculation of (. At one extreme (4.) the factorization approach is the simple point model method. At the other extreme (1.), it is the complete numerical treatment. Accordingly, its accuracy and cost may be adjusted to obtain a compromise. The most successful approximate method of space and energy dependent reactor kinetics analysis developed to date appears to be the quasistatic approximation code QX 1 developed by Ott and Meneley[36,37] at Argonne National Laboratory. Ott's results[37] show that the quasistatic approximation yields results for fast reactors that are almost identical to the "exact" results. This would cause one to question the advisability of attempting to develop  9  a more detailed method of fast reactor kinetics analysis. However, it must be noted that the "exact" standard used for this comparison is the twogroup WIGLE code. Just as Yasinsky and HenryP[] found large discrepancies in certain cases when comparing the results of point model and space dependent kinetics analyses, it is quite possible that detailed analysis of the energy variable will yield new insight on fast reactor kinetics and control. Recent work by Hitchcock, Hansen, and Henry[38] indicates that more detailed spectral analysis may indeed be necessary in reactor kinetics studies. Use of an electronic analog computer, which treats time as a continuous variable, eliminates the time step size problems inherent in a digital computer approach. The lack of storage and retrieval capability on an analog com puter, however, limits its use to one set of differential equations that is solved simultaneously. On the largest analog computers available a maximum of about 100 first order equations (50 secondorder equations, etc.) can be handled simultaneously. Due to this limit, a realistic time, space and energy dependent reactor problem cannot be solved on an analog computer. Despite the fact that much work has been done during recent years in the area of space and energy dependent kinetics, it can be seen from the above discussions that  10  a great deal of interest remains in further development in the field. The Use of Hybrid Computers Until quite recently all of the current methods of space and energy dependent reactor kinetics analysis have been implemented only on digital computers. A number of different methods for solution of partial differential equa tions using hybrid computers have been developed.[3as The hybrid computer is naturally suited for many reactor kinetics problems. Several of the space and energy dependent kinetics methods currently in use could readily be implemented on a hybrid system. Recent applications of hybrid computers to the solution of nuclear engineering problems[[4~1 and the new hybrid facilities at Westinghouse and Babcock and Wilcox among others indicate a growing interest in hybrid computation by the nuclear industry. Use of hybrid computers in the analysis of spatial kinetics problems has been reported in three recent papers. Saphier and Yiftah of the Israel Atomic Energy Commission[41. have attempted a onedimensional kinetics analysis using a discrete spacecontinuous time method. Their method is a simplified version of the technique developed by the present author. Feedback effects are not included and the amplitude scale of the problems considered is fixed and thus of quite limited range. Unfortunately Saphier and Yiftah were unable  11  to make their hybrid technique work due to hardware problems. Wozny, Bartells and Bailey at Purdue[42] used a continuous spacediscrete time procedure to study a onedimensional Xenon oscillation problem. Spatial integration was done on the analog portion of the hybrid computer using Vichnevetsky's decomposition method. [43 Carter et al. at Argonne National Laboratory[4. used a hybrid computer to tackle an ambitious two dimensional kinetics problem including temperature feedback effects. They used a fairly crude nodal representation of the core equations. These papers further indicate growing interest in nuclear appli cation of hybrid computers. A hybrid computer consists of an analog computer and a digital computer connected together with interface devices, as shown in Figure 1, in such a manner that the two may operate as one unit. Analog computers are particu larly well suited for the solution of large sets of simul taneous ordinary differential equations. The equations may be intricately crosscoupled and highly nonlinear without affecting the ease of their solution or the high speed at which the solution may be obtained. By combining the analog computer's ability to solve simultaneous differential equations quickly with a digital computer's prodigious memory and ability to make complex logic decisions rapidly, we can obtain a hybrid system that combines the best features of both analog and digital machines.  12  HYBRID COMPUTER ANALOG MACHINE INTERFACE DIGITAL MACHINE ANALOG TO CONVERTER DIGITAL TO CONVERTER DIGITAL (ADC) ANALOG (DAC) CONTROL SIGNALS Figure 1. Hybrid computer. 4I 'OU_ 1~ K [~ m ,~YY nc ~ ~YIYI   ~~CI~DI~~D ~ar~ypuMsaa~~U^~~~ j~aiwasfwevwaw  13  Hybrid computers have been used extensively in the past, primarily by the aerospace industry, tohelp facilitate what are basically large analog simulations. In recent years interest has shifted from hybrid simulation to hybrid computation. Both the analog and digital machines partici pate fully and equally in the calculations. One of the fields of major interest for hybrid computation is the solution of partial differential equations.["3 9,. The methods used may generally be divided into two categories: 1. Continuous SpaceDiscrete Time (CSDT) Methods 2. Discrete SpaceContinuous Time (DSCT) Methods The names are selfdescriptive. In CSDT methods the partial differential equations are reduced to sets of coupled ordinary differential equations in space by replacing the time de rivatives with finite difference approximations. The space equations are solved on the analog machine at each point in time. This method has been used to study Xenon oscilla tions in nuclear reactors. [42 The Xenon problem has a very slow time scale. For normal kinetics studies the need to use a very small time step makes the CSDT method impractical. The method is further limited to one space dimension. In the DSCT method the space derivatives are replaced by finite difference approximations. This yields a large set of coupled ordinary differential equations in time.  14  These are solved on the analog portion of the computer, treating time as a continuous variable. The technique developed at the University of Florida is a discrete spacecontinuous time method of direct nu merical solution to the time, space, and energy dependent neutron diffusion equation that utilizes the unique capa bilities of the hybrid computer. The method appears to be well suited to the solution of the kinetics equations but certain limitations are imposed by the characteristics of nuclear problems. Specifically, the means of data playback, the computation time, and the method's accuracy are affected by the type of problems to be solved. Auto matic rescaling is imperative due to the large range of the flux encountered in typical reactor problems. CHAPTER II DESCRIPTION OF METHOD Mathematical Model The equations of the diffusion theory model as presented in Chapter I describe the space, time and energy dependent behavior of the neutron flux in a reactor and comprise the basic nuclear model for most reactor kinetics calculations. 1) Flux 1 4 V DVQ E + Z (E' E) (E') dE' v t t oJ s + X(E) I v(E')(1 ((E')) Zf(E') (E') dE' + I Xt(E) XkC (1) 9 2) Delayed Neutron Precursor ac = X (v(E') v(E') Ef(E') D(E')) dE' C (2) where q,D,jfEt and s = f(r,E,t) Ct = C (r,t) and v,v and X = f(E)  15   16  The system feedback effects are contained explicitly in the time dependent coefficients of the above equations. The feedback terms that must be included are density changes in the fuel and coolant due to.temperature and pressure changes, thermal distortion of the core, and microscopic crosssection changes due to the Doppler effect. To cal culate the effects of these feedback phenomena, the temperature must be calculated along with the flux level throughout the reactor. An approximation to the actual temperature distribution is the twocomponent represen tation.[20] Fuel Temperature dTf Vc 1 1 1 pfCf dT ( r)q V AHh+ (T T ) (3) ff dt Vf AHU AHho f c Coolant Temperature S c t1 1 f Tc + V i c + (T TT + rq V. (4) c where q = C I dE Ef(r,E,t) (r,E,t)  17  From the temperature distributions it is possible to calculate the time dependent coefficients of the flux and precursor equations. Each of the parameters may be represented as aT c This completes the feedback loop. Because of this feedback the reactor kinetics equations are nonlinear. Standard reactor notation is used throughout. A standard method of treating the energy dependence of equation (1) is the multigroup approach. The energy range (0 to 10 Mev) is divided into a number of segments or groups. The flux equation, (1), is then integrated over each segment, one by one, obtaining for each group one equation. An approximation obtained from a static calculation is used to represent the energy dependence of the flux in the energy integrals and suitable energy averages are defined for the group constants: Ek E E (r,Et)dE' A k .k E a, (k r,t) (r,E',t) dE' ak (r,t) E E Ek kl k E 'i(r,E',t) dE' Ek1 etc.  18  The resulting time dependent multigroup diffusion equations are the beginning of most space and energy dependent reactor kinetics analyses. kth Prompt Neutron Energy Group: 1 < k < G a4, G 1 kk + G V Xsk ,k Vk D Vk k t k k =1 'k k G + Xk(l 8) v fk' k '=1 k' L + Xk XCQ + Sk (5) =1 th Delayed Neutron Precursor Group: 1 < Z < L 3C G = 8 ~ k~ X C (6) k=l k where k' Dk3 tk s ,k'IkSk and C, = f(r,t) The thermalhydraulic equations remain unchanged except for replacing the energy integral in the fission energy term with a sum over all the energy groups: G q(r,t) = C X Ef (r,t) Dk(r,t) k=l k  19  To handle the space dependence by direct numerical techniques, the space derivatives are replaced by finite difference notation using the central second difference. In cartesian coordinates, for example, d2z j+1 20 + j1 d x. (Ax.) 2 The central difference is used because it is the most ac curate of the threepoint finite difference approximations. Using this approximation the multigroup time dependent finite difference model of the reactor is obtained. (For simplicity of notation slab geometry and one space dimen sion are used.) kth energy group: d, Oj+1 24 + ~j1 G 1 k D k k k X k dt k (Axj 2 k'=l Sk'k k? G j j + k tk k + Xkl k=1 fk k' L + Xk C (7) =i kth precursor group: dC G k= Vfz k X (8) k=l k  20  where the subscripts and superscripts denote: j = space point number 1 < j < N k = energy group number 1 < k < G A = delayed neutron precursor 1 < k < L group number Since the thermalhydraulic equations contain no space derivatives, they remain unchanged. A set of tem perature equations may be defined for each space point: ThermalHydraulic Equations: j dT v PfC d (1 r)qj + AH1 3 f f dt V AHU AHho f c (9) 9p HC dT (i at Ao T) + 2 V fj H H x T. Tj + rq [] (10) in c V where G q3 = C E (t) k(t) (11) k=l k k And feedback a (t) = G(0) + (T (t)) (T(0) k k D f 2 p aj + "kapJ(t) Pi (0)J (12) ap c c c I  21  The standard procedure for solving this large set of coupled ordinary differential equations on a digital computer is to replace the time derivatives with finite difference approximations and use a standard numerical integration routine. If the equations are linear, an implicit time differencing scheme may be used. The calcu lation involves either spatial iterations or the inversion of a large but fairly sparse matrix at each time step. If the equation coefficients change with time, due to feedback effects or an external perturbation, the time steps must be kept quite short to insure accuracy. If an explicit time differencing scheme is used, computational stability requires that the time step size must be smaller than the inverse of the largest eigenvalue of the coefficient matrix, which is related to and on the order of the prompt neutron lifetime. In a fast reactor this is on the order of 107 seconds. Because of this, the large alldigital codes using direct numerical integration can be time consuming and costly. By modeling equations (7) through (11) on the analog portion of a hybrid computer, time is treated as a contin uous variable. Thus time step size limitations are avoided. This set of equations is too large to be modeled simulta neously on even large analog computers. By setting up a block of only a few space points at a time on the analog portion of a hybrid computer and using digital control and  22  memory as the analog integrates each space block sequentially, the equations can be solved quite readily. In this study the space block contained one space point. The algorithm is an iterative discrete spacecontinuous time method with variable coefficients.[46] The iterative discrete spacecontinuous time method has been investigated by Howe and Hsul 46 for nonnuclear application. Their study must be considered a preliminary investigation because the hybrid procedure was simulated using a digital computer. They considered very few space points and simple "model" problems with constant coeffi cients. In so doing this eliminated many of the physical realities that seriously affect the hybrid method. In the present study realistic reactor problems were tackled using actual hybrid hardware. Description of Hybrid Technique The solution proceeds in the following manner: the flux, precursor and temperature equations at one space block are solved on the analog computer for a short time interval. The coupling terms from neighboring space points are supplied by the digital computer which also samples and stores the solution for later playback during the calcu lation of the adjacent space blocks. With the analog computer operating under digital control the equations are  23  solved at each space block, one by one, using coupling terms stored in, and then retrieved from the digital memory. Figure 2 shows a hybrid schematic diagram. Use of the second central difference and sequential solutions means that this procedure must involve iteration. Equation (5) which describes the time behavior of the neu tron flux at space point j requires as input the solutions of 4(t) at the adjacent points j 1 and j + 1. Since the points are handled sequentially, a solution for jl (t) is already available for playback from digital storage. Dj+ (t) however has not yet been obtained. Hence one must use an estimate of what Dj+1(t) will be. Since it is only an estimate, the solution obtained for J (t) will be in error. When this procedure has been done for all space points, a first iteration .estimate of all J (t) will have been ob tained. The entire procedure must then be repeated using the previous solutions as the estimate for j+l1(t), yielding a second iteration estimate of J (t) at all space points j. The procedure is repeated again and again until subsequent iterations agree within some convergence criteria. Figure 3 illustrates the iterative procedure. After the solution has converged at all space points for the first integration interval, the final values of flux, precursor and temperature are used as initial conditions for another time interval. This procedure may be repeated to follow the solution as  24  HYBRID DIAGRAM TEMPERATURE EQUATIONS FEEDBACK COMPUTATION Figure 2. Hybrid diagram. ANALOG FLU( AND PRECURSOR EQUATIONS DIGITAL STORAGE AND CONTROL . ir  25  (p)th iteration of ( +(t) from digital storage (p)th iteration of J (t) (p+l)th iteration of i (t) being generated on the analog (p+l)th iteration of Dj1(t) from digital storage xjI xj Xj Figure 3. Iterative procedure.  26  long as is desired. An initial guess must be used to start the iteration procedure. For the first time interval the flux is assumed to remain constant. For subsequent time intervals linear extrapolation of the previous solution is used as the initial guess. External perturbations, time variant parameters and feedback effects are calculated by the digital computer and are implemented on the analog machine by means of variable coefficients in the differential equations. On an analog computer accuracy is fixed at a percentage of the machine's full scale range by noise. Hence a small signal is less accurate than a large one. Because of this, the useful dynamic range of the dependent variables is limited to about a factor of 20. If the value of a variable falls below 1/20 of full scale, accuracy * becomes very poor. The restriction to a factor of 20 change in dependent variables is a severe limitation for nuclear reactor kinetics studies. In these problems it is often of interest to follow the reactor flux levels through changes of many orders of magnitude. Therefore an auto matic rescaling feature was incorporated in the hybrid method which frees it from the usual restriction of the analog computer dynamic range. If at any time the solu tion at a space point exceeds the machine range of 100 volts, that point may be automatically rescaled and the  27  solution continued. At each point the scale factor is adjusted so that the maximum value of the solution during the time interval being considered is slightly under 100 volts. This allows the best use of the computer's avail able range. Although the automatic rescale feature is quite straightforward, it appears to be a new idea and is particularly applicable to nuclear reactor problems. A flow diagram of the hybrid algorithm is presented in Figure 4. Although the scale factors are automatically adjusted at each space point as the iterative scheme con verges upon the solution, the scale factors are constant during each integration interval. Thus it is important that some foresight be used in the selection of the length of the integration interval. This length must be such that the changes in the dependent variables do not exceed a factor of 20. This in no way restricts the total range that can be handled. Many integration intervals may be stacked end to end, each with up to a factor of 20 change. The automatic rescaling is handled by simply auctioning the scale factor to be used to the highest bidder among the variables to be used at the space point being computed. If any variable exceeds 100 volts, its scale factor is doubled, the auction repeated, and the calculation performed again. This scaling is relative  28  Figure 4. Hybrid flow diagram.  29  scaling. That is, the ratio of any two variables at a space point is not changed although their magnitudes are. Space and time dependent absolute scaling may also be added to the hybrid kinetics method quite easily. In this way the magnitude of all variables at all space points and at all times may be kept within the desired 5 to 100 volt range. Absolute scaling has not been made automatic; the scale factors at various space points and times must be selected and fixed before the solution is performed. With a hybrid computer that had digital pots and digitally con trolled remote patching for input and feedback impedances the absolute scaling could also be performed automatically. It is questionable if this would be desirable. The level of complexity of the hybrid program would be significantly raised and the same effect could be obtained by simply * using more forethought with the nonautomatic scaling before the problem was run. Analog computer scaling can be a particularly nasty problem for reactor kinetics studies because of the dif ferences in behavior of variables. The fluxes and pre cursors can vary by several orders of magnitude at different space points while the temperature distributions are generally fairly flat. The flux values at any space point may change by many orders of magnitude during a reactor transient while the precursors and temperatures will not  30  be expected to change by more than a factor of 10. Because of these differences the scaling that is appropriate at one time and a given space point may be completely wrong at a different time and/or location. By using the space and time dependent absolute scaling and automatic relative rescaling described above the hybrid kinetics method can readily handle these problems. Even the most clever techniques cannot solve all scaling problems. It is always possible that the equation coefficients dictated by the scaling are not obtainable on the analog machine. Gains in the range from 104 to 102 for summers and 102 to 104 for integrators are generally available with reasonable accuracy on analog computers. If scaling of the equations yields coefficients outside of this range, the equations cannot be handled by the analog machine. It should be noted however that the factor of 106 range in possible coefficient values is quite large. If a problem cannot be scaled to within these limits, the ma chinery is probably trying to indicate what close inspection of the physical model should have revealed. That is, that the model is unnecessarily exact, and that a suitable ap proximation may be made that will simplify the problem considerable at virtually no cost in accuracy. Some of the approximations that are often made in digital kinetic studies are particularly useful in the  31  hybrid method. If 1/vk in equation (5) is assumed to be zero for the higher energy groups, the flux equations for these groups are reduced to algebraic relationships. This approximation is useful and valid for studying slow tran sients. If the precursor concentrations in equation (5) are assumed constant, equation (6) may be eliminated. For rapid transients this approximation is valid. Use of one or the other of these approximations, whichever is appli cable, allows the analog integrator gains to be kept within the available 106 range quite easily. Selection of Iteration Scheme Iteration is a standard technique for the solution of secondorder boundary value problems.['17. Iteration is necessary because of the coupling terms from adjacent space points that arise from use of the central difference repre sentation of the second derivative in space. Use of the backwards second difference to represent the space deriva tives of the diffusion equation would eliminate the need for iteration. This scheme, however, is computationally unsta ble[47.] and will not work. The unsuitability of this method has been experimentally verified. A number of different iterative schemes exist. Three schemes were investigated analytically and experimentally for use in the hybrid solution of reactor kinetics problems.  32  These schemes were: Simultaneous Displacements, Successive Displacements and Odd/Even Iteration. Methods of iter ation acceleration were investigated in an attempt to improve upon these schemes. In order to analyze these techniques the model of the reactor was simplified. No feedback effects were included, and one energy group was used. The resulting finite difference flux equation is: 1 dO D (j j+l v a X '(  24) + ij1) _ + (1_)vE f3 + XC3 + S3 Sext dCJ B Xc dt f 2N C3 (t) = 2N i=l (13) (14) 2N  (t) = 2N eYit i(t) = e i i=l 2N C (t) = c eYt i=l Retaining only the one dominant terms of each series S(t) =3 eYit CJ(t) c3 eYkt k Note that Let  33 J(t) y C3 k k c (t) ykc P = Yi eYit = Yij (t) eykt = ykC (t) and A = yk then and so that ACj = 3zEf j ACj f c = VZf j thus equation (13) becomes 1 r< = D3 (j+l v Ax21  20 + ~jl) + ((l8)vZ Ea + vy X' + j sj f X + A ext and Let J (t) A (t) cJ (t) ACj (t)  34  D (j+l 1 r 2D 1 j + A+I v a v Ax2 + Sj =0 (15) D ext This may be represented as j+1 + aj3 + j +j = 0 (16) where aj K2 + D V +a (1(1 )] (17) D3 a V (+A f and _j A Ax2 Sj D ext The aj and B3 are assumed constant in time. The particular iteration scheme used will define a relation ship using equation (16) whereby the (p + 1)th iterate for P may be obtained from the information from the pth iter ation. Equation (16) may be represented as AT + ly = 0 (18) where A is the N x N matrix  35  ca' 1 0 0 1 0 1 a2 1 0 0 1 a3 1 A = 0 0 1 a4 1 1 0 N1 1 0 01 a N (t) p2 = and y = N N N (t) BN N is the number of space points used. To solve equation (18) without directly inverting the matrix A, iterative techniques may be defined that allow a (p + l)th estimate of the solution T to be ob tained given a (p)th estimate. The iteration equation is '(p + 1) = BO(p) + Cy (19) where B is the iteration matrix. In order for the iteration scheme to be successful, it must be consistent and convergent. Consistency means  36  that iteration using the true solution to the original differential equation must return the true solution: Dt = BDt + Cy (20) Convergence means that an error in the initial trial function will die away as the iterations progress. The rate of convergence of an iteration scheme is a measure of its efficiency and hence its relative merit. Let the true solution be defined as the pth iterate plus the pth error: t = (p) + E(p) (21) Substitute equation (21) into equation (20) and subtract equation (19) (p + 1) = BE(p) (22) Similarly with equation (18), AE = 0 Equation (22) shows that the error obeys the homogeneous version of the iteration equation and E(p) = BP~(0) The iteration scheme is convergent if the magnitudes of all the eigenvalues of the matrix B are less than one.  37  Assuming that B has a complete set of eigenvectors and eigenvalues given by Be. = IX., one may expand E(0) in terms of these eigenvectors: N E(0) = ai. 1 1 i=l and s(l) = BE(O) N B a.e i=l = C ai.X.iei 1 then N E(P) = e a.2e. i=l If the eigenvalues are ordered such that l11l > 1x21 > *** > IxI then E(p) = a1X1 e or e(p) = aleP n Alei (23) p Jn A or e where That is, the error decays as e n or e where v = kn X1 is defined as the rate of convergence. Thus the eigenvalues X. of the iteration matrix, B, determine the convergence rate of the iteration scheme. If XA11 < 1,  38  then the scheme is convergent. In the general case the result still holds where ei now includes pseudoeigenvectors of B. The three iterative schemes were examined for convergence and consistency and their rates of convergence were compared. The details are presented in Appendix A. The results are summarized below. A. Simultaneous Displacements The equation for the method of simultaneous displacements is aJj (p + 1) = j+l(p) + 0jl(p) + gj (24) The solutions at the various space points are calculated sequentially. Hence the (p + 1)th solution for j1 is available during the solution of O (p + 1). This information, however, is not used. The complete (p + l)th iteration is performed using only pth information. This scheme is consistent and is convergent as long as a01 > 2(1 ()) for all j (25) where N is the number of space points and a. is defined by equation (17). The rate of convergence is given by v = kn [(1 ) )) (26)  39  B. Successive Displacements The iterative equation for this method is aj (p + 1) = j+1(p) + jl(p + 1) + 0 (27) Again the solutions at the various space points are calculated sequentially. In contrast to the method of simultaneous displacements the information of the (p + l)th iteration is used as soon as it is available. This scheme is also consistent and is convergent if (a2 >4(1 T ) 2)2 for all j (28) The rate of convergence is given by v n 4(1 ( ) (29) C. Odd/Even Iteration In this iterative method the solutions are not obtained sequentially. Rather, the order of solution is 1,3,5,...N; then 2,4,6,...N1. For the "odd" solutions the iterative equation is aJj (p + 1) = j+p) + jl(p) + j (30) For the "even" solutions, aJ (p + 1) = Oj+1(p + 1) + 0jl(p + 1) + J (31)  40  is used. This method is equivalent to successive displace ments. It is consistent and is convergent if (aj)2 > 4(1 () 2 (32) The rate of convergence is given by v) n 4(1 (4 1 )2) (33) a2 N Since Odd/Even Iteration is equivalent to Successive Displacements and slightly more difficult to implement it will not be considered further. Squaring the convergence criteria for Simultaneous Displacements, equation (25), we obtain a > 4(1 )2)2 (34) This is the same as that of Successive Displacements, equation (28). Furthermore the rate of convergence of Successive Displacements is twice that of Simultaneous Displacements: v [n 4(l (I) 2) V [T 2 suc a Vsim n [(l (7 2 n( 2) + Zn( (T) Zn(_) + n(l l(,) 2  41  1 2 1] 2 kn( 2) + I 2) 2 1 2 ( l) + T)2 a 2 S2 Experiment verifies that the method of Successive Dis placements does indeed converge twice as fast as Simul taneous Displacements. Of the three iterative methods examined so far, Successive Displacements is the best. As a consequence it was selected for use in the hybrid method of reactor kinetics analysis. From equation (15) it is evident that the possibility exists that the solution will not converge. If the nuclear parameters of the system being studied are such that 2 1 T 2 (a) < 4(1 () )2 at any point j (35) then the largest eigenvalue of the iteration matrix will be greater than unity and the iteration will diverge. Substitute into equation (25) the definition of a , equation (17). If (2 x1 1 1 D a v A+X 2 N (36) then the iteration will not converge.  42  Several important conclusions may be drawn from equation (36). 1. If the term [Za + (1 1 v a v A+ fl is negative, it is possible that the iterations will not converge. 2. As the space mesh size Ax is made smaller, the lefthand side of equation (36) is decreased, making the inequality more easily met. Thus, the likelihood of nonconvergence is increased. 3. Convergence is improved by increasing Ea and decreasing vZf, that is, by considering a less reactive system. These observations have been experimentally verified. It was further discovered that a nonconvergent calculation could be made to converge by time scaling the equations to run more slowly on the analog computer. The existence of a regime of nonconvergence depends upon the specific parameters of the system being studied. If such a regime exists, it may be avoided by time scaling the problem down, by considering a less reactive perturbation, or by using a coarser space mesh. Figure 5 shows the results of a typical problem which was convergent if calculated for 26 space points, but was nonconvergent if 32 or more space points were used.  43  . N = 32 7 . y. .. N = 26 6 S 4 (t) at reactor center line 0 2 .3 .4 .5 TIMESECONDS Figure 5. Example of nonconvergence. CiVESS  44  To verify that the simultaneous solution of more than one flux energy group would not affect the iterative procedure, the analysis was extended to include more groups. Because of difficulties that later were found to be due to a scaling error, the effect of multigroup equations in a reflector region was of particular concern. The analysis is presented in Appendix A. The results of this analysis show that the addition of more energy groups should not affect the iterative scheme adversely. Iteration Acceleration A number of schemes are available to use in accelerating an iterative method.[47] Two schemes, overrelaxation and asymptotic extrapolation, were ap plied to the method of successive displacements without marked success. Successive overrelaxation is defined by the equation O (p + 1) = y (0+(p) + jl(p + 1) + i3] + (1 y) J(p) where y is the overrelaxation factor. It is theoretically possible to increase the rate of convergence greatly by proper choice of y.[47] The maximum rate of convergence with optimum y is  45  1 4 2 1 / v cos () v = an Ca N 1 + l 4 os2 1 2 cos ( ) aj2 N In practice the "optimum" overrelaxation factor is generally determined experimentally. Unfortunately, the expectation for great improvement in rate of convergence was not borne out. An alldigital simulation of the hybrid algorithm was written to check the entire hybrid procedure initially in the absence of any possible hard ware related problems. This simulation was used to determine the overrelaxation factor. Although this factor will depend upon the characteristics of the particular reactor being modeled, some typical reactor models yielded "optimum" factors in the neighborhood of 1.70. The reduction in number of iterations necessary was considerable as shown in Figure 6. With over relaxation factors greater than 1.70, the solution oscillated. Factors greater than 1.75 caused nondamped oscillations; the solution did not converge. When overrelaxation was used on the true hybrid model rather than the digital simulation, the solution would not converge for factors larger than 1.24. As before, the solution did not diverge until y exceeded 1.75, but for y between 1.24 and 1.75 small, bounded, but  46  .. 16bit digital simulation ...>.. 12bit simulation or hybrid method Results shown for a typical reactor example "a' I I 1.25 1.50 1.75 OVERRELAXATION FACTOR, y Figure 6. Overrelaxation results. 120 100 80 60 0 1.0  47  nondamped oscillations about the true solution occurred. By including some of the physical realities of the hybrid algorithm in the alldigital simulation, it was determined that these oscillations are caused by the 12bit reso lution of the analogtodigital converter. Why the 12bit accuracy of the ADC's should cause this oscillation is not known. These results are summarized in Figure 6. It is possible that a nonstationary technique that allows the overrelaxation factor to be reduced as the solution converges could improve the situation. This possibility was not investigated. Use of the 1.24 overrelaxation factor does result in up to a 30 per cent reduction in number of iterations necessary for convergence. The maximum benefit is ob tained when the initial guess for the solution is very poor. The number of iterations that are necessary may be reduced by using short integration intervals and by extrapolating the solution from the previous interval for use as the initial guess. The nonstationary convergence acceleration method known as asymptotic extrapolation was also tried.[46] The residual of an iteration is defined as r(p + 1) = Q(p + 1) #(p) Using the definition of the iteration error,  48  E(P) = true p) r(p + ) = + + 1) E(p) Thus from equation (23) r(p + 1) = Xr(p) and True = T(p) + r(p) + r(P + 1) + ... + r(m) S(p) + (p) + Xr(p) + X2(p) + XrrCp) S(p) + (p) thus true (p) + F(p) r(p) true F(p) r(p 1) This equation defines the asymptotic extrapolation technique. The transients introduced by the extrapolation took so long to die out that any acceleration of convergence was lost. Description of the Hybrid System The hybrid technique was developed using the Hybrid Computation facility at the University of Florida Department of Nuclear Engineering Sciences. The digital portion of this facility consists of an IBM1800 process  49  control computer with 16K of core storage. The analog portion is an Applied Dynamics80 volt analog computer with 80 amplifiers. The interface between the two machines is provided by 16 direct digital lines to and from the analog control patchboard, a 16 channel, 12bit 50 KHz Analog to Digital Converter, an 8 channel, 15bit 50 KHz Multiplying Digital to Analog Converter and a 17bit BCD digital voltmeter with remote addressing. An extensive library of general purpose hybrid subroutines is available. These subroutines were written at the University of Florida in assembly language and allow FORTRAN control over the hybrid facility. The relatively modest number of digital to analog converter channels available on the facility restricts the scope of the problem that may be attempted. Eight channels will handle at most four energy groups. By doubling up on the DACs, four precursor groups may be treated at the same time. If temperature and other equations requiring DAC channels are desired, the maximum number of energy groups is reduced. The most complicated problem run in the course of this research was a onedimensional model using two energy groups, one delayed neutron precursor group, and one fuel and one coolant temperature at each space point. Only one space point at a time was modeled on the analog. Up to  50  50 space points were used. Fuel and coolant temperature feedback are included in the model. This version of the hybrid kinetics program, which is called HYKIN 3, was modeled after and is equivalent to the WIGL 2 digital spatial kinetics code..[20] HYKIN 3 is described in detail in Appendix C. On a larger hybrid system such as the new facility at Westinghouse with 32 DAC channels up to 16 energy groups with 6 delayed neutron groups and 10 temperature equations could readily be solved. Alternately, larger space blocks could be used and the solution time reduced for a few group problem. Since the equations are solved simultaneously on the analog for each space block the addition of more equations does not increase the computation time. Implementation The hybrid reactor kinetics program is implemented in a very straightforward manner. The differential equations for the flux, precursors, and temperature at each space point as given by equations (7) through (11) are programmed on the analog computer. The feedback effects, equation (12) are calculated on the digital computer.  51  In typical reactor kinetics problems with temperature feedback, the flux behavior is very sensitive to changes in temperature. Because of this care must be used to properly scale the temperature equations to insure good accuracy. To aid in this scaling the slowly varying equations, those for the precursor and temperatures, were divided into two parts. The variables were separated into a constant term and a time dependent change term, e.g., Tf (t) = Tf(0) + AT (t) Only the change term is computed on the analog computer. The constant term is kept in the digital memory and is played back to the analog solution of the change term as required. This allows the change term to be computed with greater accuracy than could be obtained if the entire variable were modeled on the analog. The initial conditions for the analog equations are supplied by the digital machine via DAC channels. The coupling terms from space points (j 1) and (j + 1) in the flux equations are also supplied by the digital machine. The complete solutions in time of 3j+l and ~j1 are available in the digital memory and are fed to the analog as a sequence of points and slopes between these points using two DAC channels per equation. Firstorder interpolation with these points and slopes is performed on the analog computer  52  to reconstruct the coupling terms. Those coefficients of the equations that are affected by temperature feedback, and any others that vary in time, are supplied in a similar manner. The digital computer observes the reactor temperature during the pth iteration and uses this information to calculate the feedback effects and new coefficients for the (p + 1)th iteration. These are output to the analog as points and slopes at the appro priate times to coincide with the flux solution. The coefficients are reconstructed by analog firstorder interpolation from these points and slopes. The coef ficients are then multiplied by the appropriate variable using analog multipliers. Equation coefficients that are space dependent but are not functions of time are implemented using digitally controlled switches and manual pots. A much better method would be to use digital pots but these are not available on the University of Florida Nuclear Engi neering Hybrid facility. Servoset pots would not be adequate as they are generally quite slow, requiring several seconds to set. Their use would increase the run time and hence cost of the hybrid method by many orders of magnitude. A simplified schematic illustrating this implementation is shown in Figure 7.  53  = Di j+1 + j1+ l Ax2 + (1P) VE 2 f2 + (2D + E Ax2 a  (1f)v) } f 1 ( 1 + XcJ j1 i j )1 = D (x,t) 1=) Use analog interpolation with points and ,3 (t) slopes from DACs i = C (x,t) a a S= (x) Use switches or digital pots constant Use manual pot X = constant Use manual pot D +1 + 1) DAC 2 (slope) DAC 3 (point) cJ C DAC 4 (slope) DAC 5 (point) Analog implementation. where Figure 7.  54  Timing is very important in this hybrid method. It is essential that data be played back at the same time during an integration interval as it was taken. An analog clock consisting of two adjustable pulsers wired together to form a digital oscillator was used as the master time base. The digital computer was slaved to this clock and it in turn controlled the analog computer. This timing mechanism proved adequate although a crystal clock would be preferable if one were available. By carefully studying the timing of the hybrid subroutines and by making use of programmed delays and natural delays due to switch closure and amplifier reset times, it was possible to ensure that the data was played back at very close to the same time during an integration internal as it was sampled. All the digital programming was done in FORTRAN using general purpose assembly language hybrid subroutines. This allows easy modification of the programs which was imperative during development of the method. Data output is fairly slow with the present DAC control program. Each call to the routine outputs only one value and requires 274 seconds. The DAC is capable of outputting data at 50 KHz or 20 seconds each. This capability ought to be utilized by making the program output whole arrays at a single call. This would speed up the  55  hybrid program somewhat. Now that the time and other adjustments are determined for the particular pair of analog and digital computers used, it would be possible to make a special purpose hybrid kinetics subroutine. This could be writtenin assembly language and would handle all data input, output, and interpolater control functions. This would offer no advantage beyond elegance, however, since the analog interpolaters limit the speed of the solution, and they are going as fast as they can already. The essential parts of the digital portion of the hybrid program are shown below: Do for each space point Output initial conditions and begin main analog integration Do for each time sample point Do for each flux and coefficient using analog interpolation Output on DACs the point and slope for interpolation Wait for clock Output control to: enable DACs, reset analog interpolators Sample solution for all fluxes, precursors, and temperatures Wait 440 seconds for switch closure and reset Output control to: disable DACs, operate analog interpolators Reset analog CHAPTER III APPLICATION OF THE HYBRID METHOD: EXAMPLES During the period from June, 1968, to the present, work has taken place at the University of Florida on devel opment and investigation of the hybrid method of space and energy dependent reactor kinetics analysis. A large number of one and two energy group calculations have been made. The problems generally considered are those of a reactor, initially critical and at steady state, experiencing a localized perturbation at time zero. Onedimension slab geometry was used in the problems considered. The hybrid method is not restricted to any particular geometry nor is it limited to only one dimension. Howe and Hsu 6] have shown that extension of the technique to more dimensions is straightforward. From 5 to 50 space points were used. Some typical examples are shown below. Example 1: A Onegroup Model Figures 8 and 9 show the results from a hybrid computation of a one energy group model of a large thermal reactor. The reactor is the Westinghouse PWR at Turkey Point, Florida. The example shown used 24 space points across the core, which is 388 cm. thick. The perturbation  56   57  'rEGON * UNPERTURBED REGION 7 S 6 T .2 SEC4 SEC 3 2 * SPACE POINT T 0.0 SEC7, 6 12 18 0 .2 .4 .6 .8 1.0 DISTANCE Onegroup exampleFlux vs. distance. Figure 8.  58  8" SPACE POINT 7 6 THEORY 3  18 Is I ! l2 0 .1 .2 .3 .4 TIME, SECONDS Figure 9. Onegroup exampleFlux vs. time.  59  is a step change of 1.03 per cent in Ea in the leftmost onethird of the core. The perturbation has a reactivity worth of about 50 and results in a delayed supercritical transient which was followed for 400 milliseconds. The computation was performed with unity time scaling and took two integration intervals of 200 milliseconds each. The intervals each contained 20 sample points and required 35 and 16 iterations respectively to converge. The total solution run time was 6.5 minutes on the AD 80/IBM 1800 hybrid system. At $25 per hour the computation cost'was $2.70. A fairly efficient alldigital code using direct Euler integration was used to obtain a "theoretical" solu tion to the same problem using the same nuclear model. This solution is shown by the dotted line in Figure 9. The hybrid solution agrees quite well with the alldigital result. The digital solution required 3.5 minutes on the IBM 1800 and, at $15 per hour, cost $0.90. These figures are summarized in Table 1.  60  TABLE 1 Onegroup Turkey Point Example (24 Space Points) Hybrid Alldigital Euler Integration Time (Minutes) 6.5 3.5 Cost (dollars) 2.70 0.90 System AD 80/IBM 1800 IBM 1800 Example 2: A Twogroup Model The results of a twogroup hybrid computation are shown in Figures 10 and 11. The reactor model is that used by Yasinsky and Henry in their classic paper on spatial kinetics.[1] The example is a 240 cm thick, threeregion slab reactor experiencing a perturbation in the leftmost onequarter of the core. In the example shown eleven space points were used across the core. This is the number used by Yasinsky and Henry. The results are virtually identical to those obtained when the same problem was also solved with up to 50 space points.  61  07 10 105 10O SPACE POINT 0 3 6 9 1l Jjj .2 .4 .6 .8 1.0 DISTA NCE Figure 10. Twogroup exampleFlux vs. distance.  UNPE RTURBED PERTURBED IT I r 1I E TIM 8 10" 10 10 10  62  DELAYED NEUTRONS __ IGNORED Y. Q H. G I 10 uA, A POINTS SPACE 10 ,. POINT? 10 IL 33 t10 i02 3 10 I 0 I I 0 1a+++++ 0 2 4 6 8 10 TIME, MILLISEC ONDS Figure 11. Twogroup exampleFlux vs. time.  63  The reactor is initially critical and at steady state. At time zero the fast and thermal group fission cross sections in the left onequarter of the core are increased by a step of +9.6 per cent. The step increase is followed by a ramp decrease to 9.6 per cent in 10 mil liseconds. The perturbation has a maximum reactivity worth of about $8. The transient is thus highly prompt super critical. The calculation was performed using x 1/100 scaling. Twenty integration intervals of 50 milliseconds were used. The solution was sampled five times during each integration interval. Each interval required on the average nine iterations to converge. Again, an alldigital Euler integration code was used to obtain a "theoretical" solu tion for comparison with the hybrid result. The alldigital result is shown in the figures by the dotted line. The WIGLE code was also used for comparison. All three solu tions agree very closely. The hybrid computation with two prompt neutron groups, two delayed neutron precursor groups and 11 space points used 20 integration intervals with an average of nine iterations for each one. The solution required 9.6 minutes on the AD 80/IBM 1800 hybrid system. Time on the hybrid facility costs $25 per hour so the hybrid computation cost is $3.80. The same problem solved by an alldigital explicit Euler integration scheme required 0.65 minutes on  64  the University of Florida IBM 360/65. At $500 per hour, the solution cost is $5.40. Using the WIGLE code with one delayed neutron group and a fivemicrosecond time step size, the solution required 1.44 minutes on the IBM 360 and thus cost $12.00. The results are summarized in Table 2. TABLE 2 Twogroup Yasinsky and Henry Example (11 Space Points) Hybrid Alldigital Euler WIGLE Time (Minutes) 9.6 0.65 1.44 Cost (dollars) 3.80 5.40 12.00 System AD 80/IBM 1800 IBM 360 IBM 360 In their calculations Yasinsky and Henry ignored the effects of delayed neutrons. This may have been done as an economy measure since their work was done in 1965 when only slower secondgeneration digital computers were avail able and since the run time and hence cost rises at least linearly with the number of equations included at each space point. Since more equations may be added to the hybrid simulation at each space point with no additional computation  65  time cost, hybrid solutions were obtained using two delayed neutron groups as well as with the delayed neutrons ignored. There is a marked difference between the two models. In both cases the hybrid solution agrees well with the digital solution of the same model. Example 3: A Twogroup Model with Temperature Feedback The results of hybrid computation of a model of a control rod ejection accident are shown in Figures 12, 13 and 14. The reactor model contains two prompt neutron energy groups, one delayed neutron precursor and one fuel and one coolant temperature equation'at each space point. Feedback arising from the fuel temperature, Doppler effect, and coolant density changes is included. The reactor is a light water breeder reactor with a fairly fast neutron en ergy spectrum. The core has seven regions in slab geometry and is 300 cm. thick. Thirtytwo space points were used. A control rod is removed from the lefthand side of the core over a period of .5 seconds. This is modeled as a rampwise decrease in a in this region. The total worth of the a2 control rod is $9.34. The solution was followed for 600 milliseconds. The example is one used by J. B. Yasinsky in a recent comparison of spatial and point model kinetics.[481 Time scaling of x 1/10 was used together with 48 integration intervals of 125 milliseconds, each containing five sample  66  points. As many as 58 and as few as 6 iterations were required for the solution within an interval to converge depending upon the severity of the flux changes during that interval. The example is a physically realistic problem. It is typical of one that reactor designers must solve as part of a reactor safety analysis. The problem fully dem onstrates the capabilities of the hybrid system. The hybrid solution differs substantially from the results obtained by Yasinsky using the WIGL 2 spatial ki netics code. The data given in Yasinsky's paper is sparse but should be sufficient, given some ingenuity, to recon struct the complete problem definition and the appropriate initial conditions. Despite this, the solutions do not agree. Although Dr. Yasinsky was interested and sympathetic, he had discarded any further information he once had con cerning the example. Unfortunately the WIGL 2 code is not available at the University of Florida at the present time. To verify that the problem being studied with the hybrid technique was indeed the same as that considered by Dr. Yasinsky, an alldigital spatial kinetics code, ADCAL 3, was written to solve the problem digitally. This code uses the same twogroup nuclear model and the same thermalhy draulics model as that used by the HYKIN 3 hybrid code and the WIGL 2 digital code. It uses a simple Euler integration scheme.  67  The ADCAL 3 results are shown in Figures 12, 13 and 14 along with Yasinsky's results and the results of the hybrid computation. The hybrid and ADCAL 3 results agree reasonably well. They do not agree with Yasinsky's data. This indicates that the problem studied with the hybrid and ADCAL 3 codes is not precisely the same one that Yasinsky solved. One or more of the parameters of the nuclear or thermalhydraulic models differs from those used by Yasinsky. A typographical error in Yasinsky's paper or a misinterpretation of his problem definition could easily cause this difference. All the problem vari ables were initially at steady state and the spatial shapes obtained agree with those presented by Yasinsky. Because of this it is suspected that the error lies somewhere in the feedback terms that dictate the system's time response. The coefficients of the power equation which relates flux level to the time behavior of the fuel temperature, the relationship between fuel and coolant temperature and the feedback coefficients which relate the temperatures back to the time behavior of the flux are all possible culprits. Regardless of the disagreement with Yasinsky's results, the ADCAL 3 code provides a digital solution with which to compare the hybrid results. This comparison shows that the hybrid results are not exceptionally good for this example. The thermal flux shape function is shown in  68  Figure 12. The function obtained by hybrid computation with the HYKIN 3 code agrees quite well with the results from ADCAL 3. Reactor power as a function of time is shown in Figure 13. The hybrid results underestimate the peak power by 16 per cent. The hybrid solution also appears to be about 25 milliseconds early. The average fuel temperature in the reactor center is shown in Figure 14. The hybrid results underestimate the peak fuel temperature by 13 per cent. The cause of the poor accuracy in this example is suspected to be the scaling of the flux and temperature equations. During the course of the problem the flux changes amplitude by over a factor of 100 while fuel temperature changes only a factor of two and coolant temperature changes less than a factor of 1.1. If the scaling is set such that the flux is ten volts while the temperatures are 100 volts at time zero then at the flux peaks the flux is 100 volts and the temperatures are less than ten volts. This scaling is marginal since the flux is very sensitive to temperature through the feedback equations. By separating the temperature equations into a constant term and a time dependent change term and only solving for the change term on the analog the accuracy is improved. Without this scheme the problem is unsolvable. Unfortunately even this does not appear to be enough. Because the change equations contain the constant terms  69  DISTANCE (cm.) 120 180 220 270 300 0 5 10 15 20. 25 30 MESH POINT Figure 12. Example threeThermal shape function. 0 30 2.0 1.5 J 1.0 . H .5  70  20 f E \1 15 II II o1 I i \ o  o I I1 o I o r 1 I / o 0 I I I ji 5k I I 0 .1 .2 .3 .4 .5 .6 TIME SECONDS Figure 13. Example threeReactor power vs. time.  71  1600 1400 1200 1000 800 600 .1 .2 .3 .4 .5 .6 .7 TIME SECONDS Figure 14. Example threeFuel temperature vs. time.  72  the scaling will result in very small voltage values for the temperatures at the flux peaks and very small voltages for the fluxes at other times. Time dependent absolute scaling is needed. The scale factor between flux and temperature ought to be changed as the flux rises and falls. This would allow both flux and temperature voltage values to be kept between 50 and 100 volts during the entire solution which would improve accuracy for both flux and temperature. Unfortunately, implementation of this idea would require digital pots or switchable amplifier impedances, neither of which are available on the University of Florida Nuclear Engineering Sciences hybrid computer. The total hybrid run time was 4.0 hours. At $25 per hour on the hybrid facility,.this cost $100. The ADCAL 3 Euler integration code is simple and straight forward. It is also slow. Run time was 30 minutes on the IBM 360. At $500 per hour, the ADCAL 3 solution cost was $250. The WIGL 2 run time was estimated at 3.65 minutes on a CDC 6600 from data supplied by Dr. Yasinsky. The CDC machine is generally considered to be four times faster than the IBM 360/65 so the estimated WIGL 2 cost is $215. These numbers are summarized in Table 3.  73  TABLE 3 Two Groups with Feedback J. B. Yasinsky Example (32 Space Points) Hybrid WIGL 2 Alldigital ADCAL 3 Run time (Minutes 240 3.65 30 Cost (Dollars) 100 125 250 Machine AD 80/IBM 1800 CDC 6600 IBM 360 As can be seen from the examples above, the hybrid method does work and appears to be well suited to the so lution of reactor kinetics problems. It offers quite reasonable accuracy at costs that are competitive with alldigital methods. The method does have limitations, however. These are discussed in the next chapter. __ _____i__ CHAPTER IV AFFECTS OF THE CHARACTERISTICS OF NUCLEAR PROBLEMS UPON THE HYBRID TECHNIQUE The results of the work at the University of Florida indicate that the hybrid method is well suited to the solution of the kinetics equations but certain limita tions are imposed by the characteristics of nuclear problems. Specifically, the means of data playback, the computation time, and the method's accuracy are affected by the type of problems to be solved. The results are summarized below. Playback A standard method of data reconstruction used in hybrid methods is zerothorder interpolation with the playback onehalf cycle advanced.[46] This form of play back is not adequate for the class of nuclear reactor problems considered in this study. There are two reasons for this. First, the fluxes at adjacent space points are very closely coupled in reactors. Hence, the solution at a point is strongly affected by the shape of the playback at the neighboring points. The problem becomes more severe as the space mesh is made smaller. The second reason for the inadequacy of zerothorder playback is the fact that a  74 .  75  fairly slow sample rate was used. In order to maintain maximum flexibility and ease of modification in the hybrid program the coding was done in FORTRAN using general purpose hybrid subroutines. Because of this there is a minimum "sampleupdate playback" cycle time. It is about 5 milli seconds for one group, 10 milliseconds for 2 group, and 15 milliseconds for 2 group with feedback. The problems were scaled to allow the analog to run as fastas possible in order to minimize the total computation time. This means that the flux level can change appreciably between two samples. The combination of large flux changes and close coupling between space points makes zerothorder playback inadequate for all but very small numbers of space points. This point is illustrated in Figure 15. Fortunately, firstorder interpolation between data points is easily implemented on the analog computer. At each time sample point during an integration interval the digital computer outputs the value and also the slope of the fluxes at the neighboring space points via two DAC channels. An analog integrator uses these values to reconstruct ~J+1(t) and Qjl(t) as series of straight lines between time sample points. Furthermore, since the DAC channels that provide flux and precursor initial condi tions are not used during the integration cycle when the interpolation is needed, firstorder analog playback may be  76  DATA PLAYBACK SAMPLE POINTS ZEROTH ORDER '' EXPECTED SOLUTION ACTUAL SOLUTION SAMPLE POINTS / ORIGINAL CURVE /I / FIRST ORDER PLAYBACK r ... _. .I T 1 Figure 15. Data playback.  77 added without using any extra interface hardware. Use of "open loop" interpolationthat is, of simply supplying the new slope at each sample pointis not sufficiently accurate since errors accumulate. It is necessary to supply both the point and the slope. Accuracy The hybrid method is sensitive to offset or gain differences between the analog computer and the hybrid interface devices. Random error and the coarse (12bit) resolution of the ADC do not affect the accuracy but offset or gain errors which are not random have a strong effect. The reason for this sensitivity is the iterative nature of the solution. A small error repeated many times will accumulate. Figure 16 illustrates this sensitivity. Some possible sources of nonrandom offset and linearity errors are the ADC and the DAC, the analog multipliers and analog amplifiers. By checking the analog components occasionally for correct operation and by adjusting the ADC and DAC offsets and gains frequently or by using the digital computer to compensate for these problems, errors due to these sources can be kept within tolerable limits. It is readily possible to keep the total error in the solution below 5 per cent. On the University of Florida Nuclear Engineering Hybrid  78  50 N r" N U lork f m'. 100 150 OF ITERATIO0NS Figure 16. Sensitivity to offset error. 150 I00 50 0 U~ t ___  79  Computer it was found that the ADC offset would drift up and down by one or two bits from day to day. Although this is only .050 volts out of 100 volts, the high sensi tivity to offset of the iterative method made even this small error unacceptable. By measuring and compensating for the offset before each run, the problem can be elimi nated. Aggravating the sensitivity problem is the fact that the time derivative of the flux in the reactor kinetics equations is defined as a large multiple of the small dif ference between very large terms: S= v(V DVD Za + (1 )VZ ff + XC) A small error in the production, absorbtion, or diffusion terms can make a large difference in the calculated flux derivative. For less sensitive equations the fact that the hybrid computer is not ideal would have much less ef fect upon the solution accuracy. Even with the very sensi tive flux equations, the realities of the hybrid machine may be accounted for quite easily. Run time No definite numbers can be given for the computation time due to the iterative nature of the hybrid method. The run time is a direct function of the number of iterations  80  necessary for convergence. This depends upon a number of variables: the number of space points, the length of the integration interval, the degree of coupling between space points, the severity of the perturbation being studied and the initial guess used to start the iteration procedure. As was pointed out in Chapter II the iteration process may be represented as an exponential decay with the number of iterations done of the error function E(r,t). The initial value of this error is the difference between the true solution and the initial guess used. g (r,t) = qt(r,t) c1 (r,t) and Sp(r,t) = (r;,t)eVP where v is the rate of convergence and p is the iteration number When Ep+l p is smaller than some convergence criteria, p+l is assumed to be small enough and the itera tion is considered to have converged. For a particular set of It and 0 and a given convergence criterion, the error reduction necessary for convergence is eVP. Thus vp is a constant. Computation time is proportional to the number of space points to be iterated, N, multiplied  81  by the number of iterations necessary for convergence. For vp constant the number of iterations is proportional to 1/v. Thus computing time, T, is proportional to N/v: T N/v. From Chapter II, v = n A = 2 An( cos( )) So T cc a N 2 n(Z cos( )) Using cos() 1 (N)2 N 2N 2 An( ) + (n(1 1) Now nn(1 l ) T 2 N So T cc N T (2 2 n( N Jai. For the method of successive displacements I1A = 2 1 l = I 2 +x (E + P (1f)v) =2+ a v   82  L and Ax = L for onedimensional geometry so finally 1 N2 + (E + (lr)vE N2D a v . It is difficult to pick any general functional dependence out of this proportionality. Using the one group nuclear parameters from Example 1 of Chapter III, run times were estimated for various N. From this data, computation time appears to be proportional to N2*6 in the region of interest. Empirically, run time was found to vary with N2'3 for this example. These results are shown in Figure 17. Run time for Example 2 of Chapter III, which is a completely different nuclear problem, was also found experimentally to rise as N2"3, For comparison, the computation time for an alldigital code using an explicit Euler integration scheme will rise in proportion to the cube of the number of space points. Run time for an efficient implicit fewgroup code such as WIGLE may rise as slowly as N10. However, as more energy groups or other equations are added to these codes the coefficient matrix of the computation equation becomes less sparse and the dependence of run time upon N will become stronger. In the worst case it could be proportional to N2.0  83  2,6 COMPUTATION S(ARBITRARY UNI 2o 0 0 O CO 0 5 10 NUM.vR OF 20 50 50 POINTS FIGURE 17 COM PUTAT ION TI ME OF SPACE POINTS vs. NUMBER 00o0 1000 300 200 100 EXPERIMENT 50  84  These computation time estimates hold only for onedimensional calculations. The rate of convergence v depends upon Ax, the number of computations to be performed upon N. Thus the total computation time will be proportional to NLo' and to ( For two or threedimensional sys tems the same T N2*3 functional dependence is expected to hold. As estimate of the computation time for a two or threedimensional problem, however, cannot be made from one dimensional run times by simply multiplying by the ratio of the number of space points to the 2.3 power. Rather, if Ax and the other nuclear parameters are the same, the two or threedimension computation times may be estimated as simply N 1. T = T .(N2D) 2 ID J D N Unfortunately most nuclear problems are rather closely coupled in space and hence are slowly convergent; they require a large number of iterations. As the number of iterations goes up, the solution time rises and the solution accuracy falls. In order to improve upon the hybrid method the number of iterations necessary must be reduced. This is discussed further in Chapter V. The total computation time is reduced as the length of the integration interval is reduced, which agrees with the results of Howe and Hsu.[46] This is explained by the fact that the error in the initial guess made by extrapolation of  85  the previous solution becomes smaller as the extrapolation is shortened. The rate of convergence of the iterative scheme, given by equation (29) in Chapter II, is not affected by the length of the integration interval. The total error in the initial guess, however, does increase as the interval is lengthened. When the initial guess used is P(t) = GD, constant in time, then the total computation time remains about constant regardless of the length of the inte gration interval. When extrapolation is used a limit is reached, however, at which run time increases again as the integration interval is shortened. This is shown in Figure 18. Shorter integration periods cause fewer iterations to be required but also decrease the "duty cycle," the frac tion of time spent getting useful data from the analog. The minimum point will depend upon the particular problem being solved and the particular computer used. The major benefit of the hybrid method is the fact that it allows more equations to be added at each space point, up to the limit of the available hardware, with no computation time penalty. Thus the hybrid procedure allows the addition of more energy and precursor groups or more temperature equations virtually free. In an alldigital analysis the run time will vary at least linearly with the number of equations at each point.  86  500 0 400  z / 0 T 300 300  0 0/ E  00 200 300 400 INTEGRATION INTERVALM4ILLISECONDS Figure 18. Computation time vs.Data is fontegration interval. 0 problem with 26 space points. 100 0 ____ 0 100 200 300 400 INTEGRATION INTERVALM4ILLISECONDS Figure 18. Computation time vs. integration interval. CHAPTER V DISCUSSION The results presented in previous chapters show that the hybrid computer technique may be successfully applied to solution of space and energy dependent reactor kinetics problems. The method works and is easily imple mented on a very modest sized hybrid facility. The method's accuracy is quite acceptable. It has no inherent limita tions; it can readily include feedback effects, it is not limited to any particular geometry or number of dimensions. The number of energy groups, precursors, temperature equa tions and so forth that may be included in the system being modeled depends only upon the amount of hardware available. The cost of the method in its present stage of development is comptative with all digital codes if modest numbers of space points are used. The Discrete SpaceContinuous Time method is a known hybrid technique[39,41,45) and the space and energy dependent reactor kinetics problem is one that has received much attention.[1,2,537, 12,4'44,48] This study is apparently the first successful application of the hybrid technique to the kinetics problem. It also appears to be the first successful application of the iterative DSCT method to  87  