-
Notifications
You must be signed in to change notification settings - Fork 10
/
ORNL-4628.txt
23318 lines (19521 loc) · 417 KB
/
ORNL-4628.txt
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
L
3 445k 0050L48 2 {?. 106
CRIGEN — THE ORNL ISOTOPE GENERATION
AMD DEPLETION CODE
M, 1 EBell
Printed in the United States of America. Available from
National Technical Information Service
U.S. Department of Commerce
5285 Port Roval Road, Springfield, Virginia 22151
Price: Printed Copy $5.45; Microfiche $0.95
This report was prepared as an accouni of work sponsorzd by the United
States Government. Neither the United States nor the United States Atomic
Energy Commission, nor any of their employees, nor any of their contractors,
subcontractors, or their employees, makes any warranty, express or implied, or
assumes any legal liability or respansibility for the accuracy, completeness or
usefulness of any information, appsratus, product or process disclosed, or
represents that its use would not infringe privately owned rights,
. ?
ORNL~-L4628
UC-32 — Mathematics and Computers
Contract No. W-7405-eng-26
CHEMICAL TECHNCOLOGY DIVISION
ORIGEN — THE ORNI, ISOTOPE GENERATION AND DEPLETION CODE
*
M. J. Bell
e A o R H .
MAY 1973
s
o o~ e e A XY e
A Sy B AR 6 e
3 .
2
* :
Present address: USAEC, Washington, D.C., . -..-
A S s L
!.‘..1
.
4
OAK RIDGE NATIONAL LABORATORY
Oak Ridge, Termessee 37830
operated by
UNICN CARBIDE CORPORATICN cxe ]
o.0. o s comessros || INEIGTARAN
3 445k 0050L48 &2
iiil
CONTENTS
Abstract e e e e e e e e e e e
1. TIntroduction
1.1 References for Section 1 . . . .+ & v v v v v « « «
Description of Mathematical Method .
2.1 Computation of the Matrix Exponential Series
2.2 Use of Asymptotic Solutions of the Nuclide Chain
BEquations for Short-Lived Isobtopes . . . . . . . « . .
2.3 Application of the Matrix Exponential Method
for Nonhomogeneous Systems e e e e
2.4 Computation of Neutron Flux and Specific Power .
2.5 Construction of the Transition Matrix . . . . . .
2.6 References for Section 2 . . . v v v v v v v e e e 0.
Description of Nuclear Data Library . . . . . . . . .
3.1 MNuclear Properties of Cladding snd Structural Materials
3.2 Nuclear Properties of the Isotopes of the Actinide
Flements and Their Daughters . . . + « « + « « .
3.3 Nuclear Properties of the Fission Froducts
3.4 Photon Yield Library . . . . . . . . .
3.5 Miscellaneous Nuclear Properties . . . . . . . .
3.6 References for Section 3 . . . « v +v v v v 4 . o« .
Uszer's Manual for the ORIGEN Code
4.1 Description of ORIGEN Code Programs . . . . .
4.2 Description of Card Input to the ORTGEN Program . .
4.3 Input for Sample Problem . . . . .
L,
)i
Additional Programming Considerations . . . . . . .
.5 References for Section % . . . . . ¢ v v v v 0 e e ..
Page
12
13
14
18
19
19
5.
6.
iv
Glogsary of Warning and Error Messages .
AppendiX . . 0 v 4 e e e e e e e e e e
ORIGEN - THE ORNL T8OTOPE GENERATTON AND DEPLETION CODE
M. J. Bell
ABSTRACT
ORIGEN is a versatlle point depletion code which solves the
equations of radiocactive growth and decay for large numbers of
isotopes with arbltrary coupling. The code uses the matrix
exponential method to golve a large system of coupled, linesr,
irst-order ordinary differential eguations with constant coef-
icients, The general nature of the matrix exponential method
ermits the {reatment of complex decay and transmutation schemes.
An extensive library of nuclear data has been compiled, including
alf-lives and decay schemes, neutron absorption cross sections,
fission yields, disintegration energies, and multigroup photon
release data, ORIGEN has been used to compute the compositions
and radicactivity of fission products, cladding materials, and
f'uel materials in LWHs, LMFBRs, MSBRs3, and HTGRs. The applica-
tions are illustrated with calculated inventories and radiation
levels for spent fuel irradiated to a2 burnup of 33,000 MWd/metric
ton in a PWR spectrum.
o Hy
=
1. INTRODUCTION
One of the problems commonly encountered in the field of nuclear
energy is the solution of eguations involving nuclear transmutation and
decsy. To a good approximation, these nuclide chain equations can bhe
represented as a simultanecus system of linear, homogeneocus, Tirst-order
ordinary differential equations with consgtant coefficients. In many
ingtances, the matrix of nuclear transmutetion ccefficients is triangular
and the system of equations can be zgolved using the method of Bateman,l
which has been employed in a number of computer codes. =2 However, each
of these codes has generally suffered from an inasbility to treat more than
a few specific types of transmutations snd from difficulties encountered
- ':é » -
in treating 'feedback," despite some vrogress in these areas.
""Feedback” is a term used to describe a reaction in which an isotope
decays to produce one of its precursors, as is encountered in alpha
decay of the actinides, Mathematically, the appearance of such terms
results in a transition matrix that is not triangular.
N0
An alternative approach to solving the same system of eguations
/’
utilizes the matrix expcnential method, which has been employed by Pease”
and others.TnlO
Although completely general, this method has been limited
in the past by the storage that 1s regquired to generate the malrix expo-
nential series and by computational inaccuracies, Computational difficul-
ties arise because the nuclide chain eguations constitute a classic example
of a set of "stiff" ordinary differential equaticns, that is, one in which
the eigenvalues of the characteristic equation for the system are widely
Separated.ll The matrix exponential method has been employed in the
ORTIGEN code in order to exploit its extreme generality. The limitation
cn the number cf nuclides that can be treated, which results from the
necessity to store large arrays, has been overcome by twe devices:
(1) ORIGEN stores only the nonzero elements of the normally sparse tran-
sition matrix and two vectors that are used to locate them; and {2) the
cxpansion of the matrix exponential function is performed using a recursion
relation which requires storage of only cne vector in addition to the
solution. Computaticnal difficulties arising from eigenfunctions with
very large eigenvalues (corresponding to nuclides with very short half-
lives) are avoided by using asymptotic solutions of the nuclide chain
equations for the conditions of secular and transient equilibrium.
The ORIGEN code has been designed to be extremely versatile in its
applications. It is capable of computing isotopic compositilions of fuel,
fission products, and cladding in both fixed and fluid fuel reactors.
The library of nuclear data that has been compiled for use with the code
is sufficiently extensive to treat 235U and 239Pu fuels in both Tast and
233
thermal spectra, and fission of U in thermal spectra. The library also
contains multigroup photon release rates for the Tission products and the
heavy metals, which permits the calculation of gamma-ray spectra in spent
and refabricated fuels. One other impcortant feature of the code is that
The matrix exponential technigue has been developed to sclve a nonhomo-
geneous systen of equations. This fealture makes it possible for ORIGEN
to be employed in calculating the accumulation of activity in processing
plants, in waste disposal operations, and in the environment.
The remainder of this report consists of three major sections, Section
2 develops the mathematical techniques used in the ORIGEN cowputations.
Section 3 summarizes the information that is contained in the library of
nuclear data. Section 4 examines the function of each subroutine in the
program, discusses programming considerations, describes the input for
the code, and givesg input for a sample problem, Output for a sample
calculation is included in the Appendix. The individual sections are
intended to be independent, so that Section U4 can be employed as a user's
manual without reference to preceding sectlons.
A complete code package has been deposited with the Radiation Shielding
Tnformation Center. Inguirieg or requests for the code may be mailsd to:
Codes Coordinator
Radiation Shielding Information Center
Oak Ridge National Laboratory
- - 5
Oak Ridge, Tennessee 37830 $ c Al
J
P
or telephoned to:
Area code 615, W83-8611, ext. 3-694k4, or
FTS XX - 615-L83-60LL,
Acknowledgments. — The author is deeply Indebted to a number of
individuals who contributed to the development of the ORIGEN code over
a period of several years. The program was initiated by J. P, Nichols,
who supervised the development of the code. Much of the nuclear data was
compiled by E. D, Arnold” with the assistance of H. P, Soard. The author
was assisted in the programming by R. S. Dillon "~ and Mrs. L. G. Knauer,*
and is deeply grateful for the efforts of Mrs. 1. G. Loope, who prepared
the manuscript. Editing was done by C. W. Kee, who now has the task of
data revisions and program modifications. TFinally, the author appreciates
the encoursgement and suggestions of a number of userz of the program who
have patiently awaited the publication of this report.
*
Pregsent address: Osk Ridge Gaseous Diffusion Plant.
Present address: AECOP, Oak Ridge, Tenmn.
TPresent address: Robertsville Junior High School, Cak Ridge, Tenn.
1.
z.
4=
i0.
11,
1.1 References for Section 1
H. Bateman, Proc. Cambridge Phil. Soc. 15, 423 (1910).
M. P. Leitzke and H. C. Claiborne, CRUNCH — An TBM-70L Code for
Calculating N Successive First Order Reactions, ORNL-2958
(Cet. 24, 1960).
D. R. Vondy, Development of a General Method of Explicit Solution
to the Nuelide Chain FMquations for Digital Machine Caleunlations,
ORNL-TM-361 (Oct. 17, 1952).
H. H. Van Tuyl, ISOGEN — A Computer Code for Radioisoctope Generation
Calculations, HW-83755 (1964).
R. 0. Gumprecht, "Computer Code RIBD," Trans. Am. Nucl. Soc. 12(1),
1h1-hb2 (June 1969),
L. Peage, DEEMS, A FORTRAN Program for Solving the First Degree
Coupled Differential Equations by Expansion in Matrix Series,
TSI-49 (October 1963).
H. H. Poynter and J. Suez, "Automatic Digital Set-up and Scaling of
Analog Computers,” Proc. 1963 Joint Automatic Control Conf., . 156.
B. H. Duane, in Physics Research Quarterly Report, Ccl.-Dec. 1963,
HW-30020 (196L).
H. B. Krug, J. E. Olhoeft, and J. Alsina, Chap. 8 in Supplementary
Report on Evaluation of Mass Spectrometric and Radiochemical Analyses
of Yankee Core T Spent Fuel, Including TIsotopes of Elements Thorium
through Curium, WCAP-56086 (August 1969).
S. J. Ball and R. K. Adams, MATEXP, A General Purpose Digital Computer
Program for Solwving Ordinary Differential Equations by the Matrix
Exponential Method, ORNL-TM-1933 (August 1907).
M. ®, Fowler and R. M. Warten, 1BM J. Res. Develop. 11, 537 (1967).
2. DESCRTPTION OF MATHEMATICAL METHOD
A general expression for the formation and disappearance of a nuclide
by nuclear transmutatlon and radioactive decay may be written as Tollows:
N v
I Lij?\jxj + ¢ Z IikUka - (?\i + qfioi).zs.i (1 =1,...N), (1)
J=1 =1,
where Xi ig the atom density of nuclide i, hi is the radioactive disinte-
gration constant for nuclide 1, S; is the spectrum-averaged neutron
absorption c¢ross section of nuclide i, and 4.. and fik are the fractions
of radicactive disintegration and neutron absorption by other nuclides
which lead to the formation of species 1. Also in Eg. (1), ¢ is the
position- and energy-averaged neutron flux, which is also assumed to be
constant over short intervals of time, Rigorously, the system of equations
described by Eq. (1) is nonlinear since the neutron flux will vary with
changes in the composition of the fuel. However, the variation with time
is slow and, if the neutron flux is considered to be constant over short
time intervals, the gystem of Eq. (1) is a homogeneous set of simultaneous
first-order crdinary differential equations with constant coefficients,
which may be written in matrix notation:
»
f-ax. (2)
Equation (2) has the known solution
X - exs (48) £(0) (3)
S
where %(O) is a vector of initial atom densities and A dis a transition
matrix containing the rate coefficients for radiocactive decay and neutron
capture. The function exp (At) in Eg. (3) is the matrix exponential
N
2
function, a matrix of dimension N7, which is defined as
(46)° D (at)”
exp (A) = L+ At + o+ ... = ) o ()
if one can generate this function accurately from the transition matrix,
then the solution of the nuclide chain equations is readily obtained,
2.1 Compuctation of the Matrix Exponential Series
Two principal difficulties are encountered in employing the matrix
exponential technique to solve large systems of equations: (1) a large
amount of memory is required to store the transition matrix and the
matrix exponential function, and (2) computational problems are encoun-
tered in applying the matrix exponential method to systems of equations
with widely separated elgenvalues. The generation and storage of the
transition matrix are explained in Sect. 2.5. 'The computation and storage
of the matrix exponential function have been facilitated by developing a
recursion relation for this function which does not require storage of
the entire matrix, Thus it is possible to derive an expression for one
nuclide in Eg. (3) which is given by
@©
n
Xi(t) = ) Ci 3 (5)
]
n=0Q
n . . .
where Ci 1s generated by use of a recursion relation
0
¢ = x;(0) (62)
N
ntl 0t ‘ n -
Cl Con+l E; aijcj 2
J=1
Here, aij is an element in the transition matrix that is the first-order
rate constant for the formstion of species i from species j. This algorithm
} N n . .y . i
requires storage of only one vector C° in addition to the current value
of the solution.
In performing the summation indicated by Eq. (5) it is necessary to
ensure that precision in the answer will not be lost due to the addition
and subtraction of nearly equal large numbers. In the past, this objective
has been accomplished by scaling the time step by repeatedly dividing
by 2 until the norm of the matrix is less than some acceptable small
value, computing the matrix exponentizl functicn Tor the reduced time
step, and repeatedly squaring the resulting matrix to obtain the desired
—
tine step.
Such a procedure would be impractical for a computation
involving large numbers of nuclides (many of which have short half-lives)
corresponding to large norms of the é} matrix. However, it is ,just for
these short-lived isotopes that the conditions of secular and transient
equilibrium are known to apply. Thus, in the compubtations performed by
ORIGEN, only the compositions of those nuclides whose diagonal matrix
elements are less than a predetermined value are computed by the matrix
eprnential method. The concentrations of the isotopes with large diagonal
matrix elements are computed using an analytical expression for the condi-
tions of secular or transient eguilibrium, as described in Sect, 2.2.
. h 2
Lapidus and Luu53 have shown fthat the accuracy of the computed mabrix
exponential function can be maintained at any desired value by controlling
the time step such that the norm of the matrix At is less than a predeter-
mined value which 1s fized by the word length of the digital computer
used in the calculations. They define a norm of the matrix é, denoted
by [AT, which is given by the smaller of the maximum-row sbsolute sum
, max E:
i 7
J
and the maximum~column absolube sum:
[A] = min {max >: ‘a .
. i
I3
aijl} ’ (7)
where Iaij’ denotes the absolute value of the element ajj' They show that
the maximum term in the summation for any element in the matrix exponential
n
- . n . .
function cannot exceed ai oo where n is the largest integer not larger than
[A7t, Consideration of the word length of the computer used to perform
the calculations will indicate the maximum value of n that can be used
while obtaining a desgired degree of significance in the results. Using
double precision arithmetic, the IBM 360 cperating system can perforn
operations retaining 16 significant decimal figures. In the ORIGEN code
the norm of the transition matrix is restricted to be legs than
[A] < -2 1n 0.001 = 13.8155, so that the maximum term that will be
celculated will be approximately 49,000. Thus, a value as small as
exp (-13.8155) = 10"6 can be computed, while retaining five significant
figures. A sufficient number of terms must be added to the infinite
summation given by FEg. (5) to ensure that the series has converged. The
m
mth term in the series for e[A] is equal to [i? , waich, for large values
-1/2
of m, can be approximated by (Léfifs)m(Bmfl) 1/ using Stirling's approxi-
mation. The value of the norm, [A], is calculated by the code; and m is
set equal to the largest integer in % [A] + 5, which has been determined
as a "rule of thumb" for the number of terms necessary to limit the error
to <0.1%. Thus, for [A] equal 13.8155, 53 terms will be reguired in the
summation., The absolute value of the last term added to the summation
in this case will be less than 6.4 x 10_10, which is sufficiently small
compared with 10"6. It has been observed that the norm is usually less
than its maximum value and, in most cases, 3C or fewer terms are required
Lo evaluate the series.
It has been mentioned that, in previous applicationgsof the malrix
exponential method, the restricticn of the size of the norm of the
transition matrix necessary to treat nuclides with large elgenvalues
was accomplished by repeatedly dividing the matrix by 2, and the final
value of the matrix exponential function was obtained by repeatedly
squaring the resulting intermediate matrix exponential function. In
the present application, the suggestion of Ball and Adamsh that the
transitions involving isotopes with large decay constants be considered
"instantaneous" was adopted; that is, if A » B - C and if the decay
constant for B is large, the matrix is reformulated as if C were formed
from A directly, and the concentration of B is obtained by an alternative
technique. Similarly, if the time constant for A is very large, the
transition matrix is rewritten as if the amcunt of isotope B initially
present were equal to A + B, and only the transition B - C is obtained
by the matrix exponential technique. This reduction of the transition
maltrix and the generation of the solution by the matrix exponential
method are performed by the subrcutine TERM (see Sect. 4.1).
2.2 Use of Asymptotic Solutionz of the Nuclide Chain
Equations for Short-Lived Isotopes
The numerical technigues described in Sect. 2.1 are applied only
to obtain the golutions for isotopes that are sufficiently long-lived
to satisfy the criterion thal the norm of the transition matrix be less
than 2 1ln 1000. Short-lived isotopes are treated by using linear combi-
nations of the homogeneous and particular solutions of the nuclide chain
equations that are computed using alternative procedures. The guantity
of a short-lived nuclide (originally present at the begimning of an
interval) that remains at the end of the interval is computed in subrou-
tine DECAY using a generalized form of the Bateman equations which treats
an arvitrary forward-branching chain. The generalized treatment is
achieved by searching thréugh the transition matrix and forming a queue
of all short~lived precursors of an isotope. The Bateman equation solution
ig then applied to this queue. The queue is terminated when an isotope
having no short-lived precurscrs is encountered. The algorithm also has
provisions for treating two isotopes with equal eigenvalues and for
treating cyelic chains.
Bateman's solution for the ith member in az chain at time t may be
written in the form:5
. ~ds b
Ni(t) = Ni(O)e
i-1 i-1 . -
Py omoy | ) 2yt e () Pt |
k — (d' - d') uj"}l)j ‘1‘ d]:-] - d-' ?
k=1 e 1 J n=k R J
J n#]j
where NJ(O) is the amount of isotope J initially present and the members
of the chain are numbered consecutively for simplicity. This method of
i-1
solution used the conventicn that T[]
n=kK
and that the empty product (x = 1) is equal
an+l,n is equal to the product
"k1,k, Tkie, kel P4,1.17
to unity. The notation a. 5 Lor the first~order rate constant is the
-2
10
il
1
w©
same as that described in Sect. Z.1, and di In the present
spplication, Egq. (8) is recast in the form
~dit
N, (t) = I, (O)e
-1 i1
qvfi i-1 Bil 1l * exp (mdjt) -
+ ZJ N (0) 1 2 E: d
k 9y ] (ds
k=1 n= jak
11
by multiplication and division by N d . The first product in Eg. (9)
n=k
has significance because it is the fraclion of atoms of isolope k that
follow a particular sequence of decays and captures. If this product
becones less than 10—6, contributions from nuclide k and its precursors
to the concentration of nuclide 1 are neglected. The inner summation in
Eq. (9) is performed in double precision arithmetic to preserve accuracy.
This procedure is unnecessary for evaluating the outer summation because
all the terms in this sum are known to be positive. The difficulties
described by Vondy in applying the Bateman equations for small wvalues of
dit do not occur in the present application since, when this condition
occurs, the matrix exponential solution is employed. The matrix exponen-
tial method and the Bateman equations complement each other; that 1s, the
former method is qguite accurate when the magnitude of the characteristic
values of the equations tc be solved 1s small, whereas the Bateman
solution encounters numerical problems in this range. TFor the case
where two isotopes have equal removal constants (d; = dj), the second
summation in Egq. (9) becomes:
i-1 ]
\ -d st i-~1 d,
E; d.t e I —— (10)
:{_ Il i
ik n { J
n#gj
An analogous expression is derived for the case when 4, = dj. These forms
of the Bateman equations are applied when two isotopes in a chain have
the same diagonal element or when a cyclic chain 1s encountered, in which
case & nuclide is considered to be its own precursor.
11
In the situation where a short-lived nuclide has a long-lived
precursor, a second alternative solution is employed. 1In this instance,
the short-lived daughter is assumed to be in secular equilibrium with
ite parent at the end of any time Interval. The concentration of the
parent is obtalined from subroutine TERM, and the concentration of the
daughter is calculated in a subroutine named EQUIL (see Sect. 4.1) by
setting Eq. (3) equal to zero:
[
S
N
. =0 = :
Xy ) E: aij Xj . (11
J=1
Equation (11), which is a set of linear algebraic equations Tor the
concentrations of the short-lived isotopes, is readily solved by the
Gauss-Seidel iterative technique.6 The coefficients in Eq. (11) have
the property that all the diagonal elements of the matrix are negative
: and all off-diagonal elements are positive. The algorithm involves
inverting Eq. (11) and using assumed or previously calculated values
for the unknown concentrations to estimate an improved value, that is,
N
+1 . i
& :....;.Ll . (12)
1 ciii :LJ J
3=1
JH
The iterative procedure has been found to converge very rapidly since,
for these short~lived isotopes, cyclic chains are not usually encountered
and the procedure reduces to a direct sclution.
2.3 Application of the Matrix Exponential Method
for Nonhomogeneous Systems
Certain problems that involve the accumulation of radicactive
materials at a constant rate and are of engineering interest require
the solution of a nonhomogeneous system of first-order linear, ordinary
differential eguations. I[n matrix notation, one writes
%: ,\X_,Jr%' (13)
Q=
This set of equations has the particular solution
x(t) = [exp (88) - 1A B (14)
~
¢
. ~1 . . . s axo s . -
provided that A = exists. Substituting the infinite series representation
for the matrix exponential function, one obtains
At (at)”
X(6) = [L+ Sr+ =+ .1 B (152)
A 2 \
=\ /) T ) B (15b)
m=0
The particular solution may also be expressed as the sum of an infinite
serles:
whose terms are generated by use of a recursion relation
1
Dy =bst, (17a)
N
a+1 t ja
Yoot o 43 DJ (170)
J=1
Once again, the algorithm isg applied only to long-lived nuclides, and the
concentrations of the short-lived nuclides are obtained by an alternative
technique. In this situation, Hg. (11) is modified to the form
N
¥, = 0 = 24 Ay Xs b, (18)
J=1
and is solved by the Gauss-Seidel method. After the homogeneous and
partlicular sclutions have been cobtalned, they are added to obtain the
complete solution of the gsystem of eguations.
2.4 Computation of Neutron Flux and Specific Power
In order to compute changes in fuel composition during irradiation
at constant power, it is necessary to take into account changes in the
neutron flux with vime as the fuel is depleted. At the start of the
computation, the known parameters are the initial fuel composition and
the constant specific power that the fuel must produce during a time
interval. The instantaneous neubtron flux may be related to the fuel
composition at a fixed time by the equation
P 3.20x 107 5, (19)
where P is the specific power, in MW per unit of fuel; ¥ is the macroscopic
fission cross section, in cm? per unit of fuel; and ¢ is the instantaneous
neutron flux, in neutrons per cm?-sec. The constant in Zg. (19) is derived
by assuming a value of 200 MeV per fission. An approximate expression for
the value of the neutron flux as a function of time is obhtained by expansion
in a Taylor series about the start of the interval:
6
. 2 o . 2 . 0_1
@) = 3.125 x 1017 p r|__z:(o)'l L 2O 5 (’*[2(0” - 2(0) MO)) bt (
$(0)7 5(0)°
or
s(t) = 9(0) {1 - t%%—%+§-—(—?ufil) ] (200)
2(0)2 -
14
The average neutron flux during the interval is obtained by integrating
over the interval and dividing by t:
7
]
2(0) . t2 /25(0)” - £(0)%(0)
L..l - ‘ 1(O2 " ‘é_' ( 2(0)2 )> * ...] : (21)
oo et
¢ = ¢(0)
[\/‘
Here, the notation £(0) is used for the macroscopic fission cross section
at the start of the time interwval, and i(O) and i(O) are the first and
second time derivatives evaluated at the start of the interval. The
values of é(O) and f(O) can be evaluated since %(O) = A X(0) and
%(O) = éfi%(o) = é? X(0). FEquation (10) is used in the computer program
in subroutine FLUXO to estimate the average flux during an interval,
based on conditions at the start of the interval. The term involving
the second derivative is only employed for the first time interval where,
for some isotopes, X(O) is zero but %(O) is nonzero. The average power
produced during a time interval for a fuel in a fixed neutron flux is
estimated from the initial composition using a similarly derived equation:
. .o z
P = 3.20 x 1077 px(0) [1+ £(0) £+ Ho) -+ ...] (22)
In order for the above procedures to estimate the average neutron flux
or average specific power correctly, the changes in neutron flux during
thie interval must be relatively small. If the average value of either
of these guantities differs from the initial value by more than 20%, a
message Will be printed out advising the user to employ smaller {ime
increments.
2.5 Construction of the Transition Matrix
An extensive library of nuclear properties of radiocactive isotopes
has been compiled for use with the ORIGEN code. {See Sect. 3 for details.)
The data are in the form of half-lives, fractions of transitions that
produce a given nuclear particle, cross sections, and fractions of
absorptions thatl yield certain particles. These data are read from
15
the library tape and processed into a form for use by the mathematics
routines in subroutine NUDATA (see Sect. 4.1 for description of NUDATA),
nuclides uvsing the present ccde., However, stralghtforward construction
of a generalized transition matrix would require the storage of an
800 x 800 array, which would tax the storage capacity of the largest
computers available today. On the other hand, the transition matrix is
normally very sparse, and storage requlrements can be reduced gubstantially
by storing only the nonzero elements of the matrix and two relatively
small vectors thatl are used to locate the elements. Subroutine NUDATA
is also used to generate the compacted transiftion matrix and the two
storage vechbors.
Subroutine NUDATA processes The library tape by reading a six-digit
identifying nuwber, NUCL(I), for each nuclide, Tollowed by the half-life
and the fraction of each decay that occurs by several competing processes.
A zecond card, which contains neutron absorphtion cross sections for (n,y),
o
(n,@), n,p), (n,2n), (n,3n), and (n,fission) reactions for one of four
reactor spectra, is then read. The six~diglt identifying number is equal
to Z * 10,000 + W * 10 + I3, where 7 is the atomic number, W is the atomic
welght (in integral atomic mass units), and I3 is either O or 1 to indi-
cate a ground state or a metastable state, respectively. This information
is processed into a compacted transition matrix, as described below.
first, the half-life is used to calculate the radiocactive digintegra~-
ct
ion constant, A. IFirst-order rate constants for various competing decay
processes are calculabted by multiplying A by the fraction of transitions
to that final state. The product nueclide resulting from each nuclear
transition is next identified by addition of a suitable constant to the
six-digit identification number for the parent nuelide. {(For example,
for a B decay, 10,000 is added to the parent identifier; for nsutron
capture, 10 is added; or for isomeric transition, the quantity -1 is
added.) Two arrays are constructed: the first, NPROD(J,M), contains
all the products which can be directly formed from any nuclide J by the
transitions considered in the library; and the second, COEFF(J,M), con-
tains the Tirst-order rate constants for the corresponding transition
93]
16
When these arrays have been ccnstructed, a search of the NPROD array is
conducted to identify all the parents of a given nuclide I. [Nuclide J
is a parent of I if NPROD(J,M) equals NUCL(I) for any reaction of type M.]
When a parent of nuclide 1 has been located, the value of the correspond-
ing coefficient aij in the transition matrix is equal to COEFF(J,M).
However, direct storage of aij in a square array would require an exces-
sive amount of storage. Hence, this procedure 1s avoided by incrementing
an index, N, cach time a ccefficient is identified. The coefficients
are stored sequentially in a one-dimensional array, A(N); the value of
J is stored in another one-dimensional array, LOC(N); and the total number
of coefficients for production of nuclide I are stored in a third array,
NANO(T). VWhen all of the coefficients for every nuclide have been stored,
the NfiNO(I) array is converted to indicate the cumulative number of matrix
coefficients for all the isotopes up to and including I [i.e., N¢NO(T + 1)
= NYNO(I) + N@NO(Ll + 1) for all I greater than 1]. After this procedure
has been executed, the NfiNO array is a monoctonically increasing list of
integers whose final value is the number of nonzero, off-dlagonal matrix
elements in the transgition matrix. This final value is preserved separ-
ately as the variable NON, TFor computational convenience, the wvalues of
the diagonal .matrix elements are stored in a separate vector, D(I). To
perform the multiplication of the transition matrix by a vector (e.g.,
N
ii = ¥ aijxj), as is required to execute the algorithm described in
h J=1
Sect. 2.1, the operations described in the flow chart given in Fig. 2.1
are employed.
Two types of data In the nuclear library require a departure from
the procedure Jjust descrived, In the case of neutron-induced reactions,
1t i1s necessary to specify the neutron flux before first-order rate
coefficients can be calcvlated as products of flux and cross sections.
At the time the matrix is generated, the neutron flux is unknown., Also,
to perform a fuel depletion calculation, the flux must be permitied to
vary with time. Thus, when the nonzero, off-diagonal matrix elements
for isotope 1 are stored, all those for formation by radicactive decay
are grouped first and are Tollowed by those for formation of T by neutron
17