ISTANBUL TECHNICAL UNIVERSITY ★ ENERGY INSTITUTE
DEVELOPMENT OF A NODAL METHOD FOR THE SOLUTION OF THE NEUTRON DIFFUSION EQUATION
IN CYLINDRICAL GEOMETRY
M. Sc. Thesis by
Mehmet MERCİMEK, Nuclear E. Eng.
Department: Energy Science and Technology Programme: Energy Science and Technology
ISTANBUL TECHNICAL UNIVERSITY ★ ENERGY INSTITUTE
DEVELOPMENT OF A NODAL METHOD FOR THE SOLUTION OF THE NEUTRON DIFFUSION EQUATION
IN CYLINDRICAL GEOMETRY
M. Sc. Thesis by
Mehmet MERCİMEK, Nuclear E. Eng. (301061015)
Date of submission: 5 May 2008 Date of defence examination: 10 June 2008
Supervisor (Chairman): Prof. Dr. Atilla ÖZGENER
Members of the Examining Committee: Prof. Dr. Melih GEÇKİNLİ (I.T.U.) Prof. Dr. Cemal YILDIZ (I.T.U.)
İSTANBUL TEKNİK ÜNİVERSİTESİ ★ ENERJİ ENSTİTÜSÜ
SİLİNDİRİR GEOMETRİDE NÖTRON DİFÜZYON DENKLEMİNİN ÇÖZÜMÜ İÇİN NODAL BİR
YÜKSEK LİSANS TEZİ
Nükleer E. Müh. Mehmet MERCİMEK (301061015)
Tezin Enstitüye Verildiği Tarih: 5 Mayıs 2008 Tezin Savunulduğu Tarih: 10 Haziran 2008
Tez Danışmanı: Prof. Dr. Atilla ÖZGENER
Diğer Jüri Üyeleri: Prof. Dr. Melih GEÇKİNLİ (İ.T.Ü.) Prof. Dr. Cemal YILDIZ (İ.T.Ü.)
Nodal m ethods are fast and accurate m ethods which com bine attractive features of the finite elem ent m ethod as well as of the finite difference method. In this work, a nodal m ethod has been developed in cylindrical geometry which gives acceptably accurate results for realistic problems.
I w ould like to express my sincere gratitude to my supervisor Prof. Dr. Atilla ÖZGENER. His invaluable knowledge, expert guidance and care allowed this work to be completed.
Special thanks and appreciation go to Prof. Dr. Bilge Ö ZG EN ER and Prof. Dr. A kif ATALA Y for their great support, advice and time.
I also w ish to thank my colleagues from Istanbul Technical University Institute of Energy; especially physics engineer and ITU TRIGA M A R K II reactor operator M ehm et GENCELİ for his guidance, support and encouragement.
I w ould also like to recognize my dear friends A bdülham it ERTUNGA, Erdal KEKLİK, Ram azan A K Y Ü Z for helping me in their ow n unique ways.
Finally I would like to thank my parents, M ehm et and Yüksel M ERCİM EK, my dear brother M ustafa M ERCİM EK and my sister Serpil M ERC İM EK for encouraging me throughout my life.
TABLE OF CONTENTS
TABLE OF CONTENTS iv
LIST OF ABBREVIATIONS vi
LIST OF TABLES vii
LIST OF FIGURES viii
LIST OF SYMBOLS ix
1. INTRODUCTION 1
1.1 Diffusion Theory 1
1.2 Num erical M ethods for Solving the Neutron Diffusion Equation 2
1.3 Nodal M ethods 3
1.4 Finite Elem ent M ethod 4
1.5 Objectives of the W ork 4
2. NODAL FORMALISM IN CYLINDRICAL GEOMETRY 6
2.1 Cell and Edge A veraged Quantities 6
2.2 Nodal Balance Equation 8
2.3 Transverse Integration 10
2.4 Nodal Expansion M ethod 11
2.4.1 Construction of Polynomial Basis 11
2.4.2 Fick’s Law 15
2.5 Iterative Solution o f the Nodal Equations 20
2.5.1 M atrix Equation 21
2.5.2 One-Node Form ulation 22
2.5.3 Two-Node Form ulation 23
2.5.4 Three-Node Form ulation and General M atrix Form 24 2.5.5 General M atrix Equation with Reflective Boundary Condition 26 2.5.6 Fission Source Iteration in M ultigroup Diffusion Equations 27
3. NUMERICAL APPLICATIONS 29
3.1 One-Group, Bare, Hom ogeneous Reactor 29
3.1.1 Analytical Solution 29
3.1.2 Solution with Nodal M ethod 32
3.1.3 NEM R and QFEM R Results 34
3.2 One-Group, Reflected Reactor 37
3.2.1 Analytical Solution 37
3.3.2 NEM R and QFEM R Results 48
3.4 TRIGA M A R K II Reactor 51
4. CONCLUSION 56
APPENDiX A A MANUAL FOR NEMR 59
A1 Input List 59
A2 Description of NEM R Subprograms 61
APPENDIX B COMPUTER PROGRAMS 64
LIST OF ABBREVIATIONS
FEM : Finite Elem ent M ethod NEM : Nodal Expansion M ethod ODE : Ordinary Differential Equation RHS : Right Hand Side
LIST OF TABLES
Table 3.1 Iteration Steps o f N ew ton’s M ethod for the Solution of (3.7)... 31
Table 3.2 Calculated Values o f the Variables in the Two-Node E qu atio n s... 32
Table 3.3 Radii and Related Surface A reas...33
Table 3.4 keff Results of QFEM R and N EM R Program s...34
Table 3.5 Average Fluxes and Respective E rro rs ...36
T able 3.6 Iteration Results of Newton M ethod ... 39
Table 3.7 Effective M ultiplication Factors Calculated by 3 M ethods...41
Table 3.8 Average Fluxes w ith Their Errors in the Fuel and the R eflector... 44
Table 3.9 Iteration Steps in Two-Group P ro b lem ...47
Table 3.10 Results of QFEM R and NEM R Programs for Two-Group R eactor...48
Table 3.11 H om ogenized Fast Group Cross-Sections for TRIGA R eactor...52
Table 3.12 H om ogenized Thermal Group Cross-Sections for TRIG A R eactor...52
LIST OF FIGURES
Figure 2.1 Nodal M esh Im posed on One-dim ensional Cylindrical D o m a in ...7
Figure 3.1 One-dim ensional, Bare, H om ogeneous, Cylindrical R e a c to r... 29
Figure 3.2 Cylindrical Reactor with T wo N odes...32
Figure 3.3 keff Results of FEM and N E M ... 35
Figure 3.4 Flux Distributions Along the Radial D istance... 36
Figure 3.5 V ariation of keff with Respect to Total N um ber of N o des...42
Figure 3.6 keff versus N um ber of Nodes in The Fuel R eg io n ... 43
Figure 3.7 Flux Distribution in the Reflected R eactor... 44
Figure 3.8 Effective M ultiplication Factor Obtained by Three M ethod... 49
Figure 3.9 Fast and Thermal Flux D istributions... 49
Figure 3.10 Thermal to Fast Flux Ratios Along the R a d iu s ...50
Figure 3.11 ITU TRIGA M A R K II Reactor Core D ia g ra m ...51
Figure 3.12 TRIGA Reactor Flux Distribution (N EM R)...54
LIST OF SYMBOLS
J g ( Î ) : N eutron current density vector for group g D g : Diffusion coefficient for group g
: N eutron flux for group g
: M acroscopic rem oval cross-section : M acroscopic absorption cross-section
^s,g^g’ : M acroscopic scattering cross-section from group g to group g ’ : M acroscopic fission cross-section
keff : Effective m ultiplication factor
X : Fission spectrum
V : Average num ber of neutrons released per fission
Q : Neutron source
ri-1/2, ri+1/2, ri : Radii at node edges and center respectively
Ai : ith node width
j + : Right-going partial current j- : Left-going partial current
Si+1/2, Si-1/2, Si : Surface areas at node edges and center respectively
Pl : lth order polynomial % : Local variable
ai : Expansion coefficients for polynom ial of lth order A (3*3) : 3 by 3 m atrix constituted from one node form ulation
: 3 by 3 m atrix with all elem ents are zero £ : Convergence param eter
B : Buckling term
Jo, J1 : Bessel functions of first kind, zeroth and first order respectively
wf : Recoverable energy per fission : Average flux
X : M axim um eigenvalue, keff
SİLİNDİRİR GEOMETRİDE NÖTRON DİFÜZYON DENKLEMİNİN ÇÖZÜMÜ İÇİN NODAL BİR YÖNTEM GELİŞTİRME
N ükleer reaktörlerin birçok fiziksel özelliği nötron difüzyon teorisi ile anlaşılm aktadır. Difüzyon teorisinin geçerli olabilm esi için reaktör ortamını oluşturan binlerce küçük malzem e, ortalam a tesir kesitleri ve difüzyon katsayıları kullanılarak hom ojenleştirilir.
Bu hom ojenleştirm e işlem ine rağm en reaktör kalbi yine de oldukça heterojen bir ortam oluşturur. B u heterojenlik yakıt demetleri arasındaki yakıt m iktarları farkından, yanıcı zehirlerden, kontrol çubuklarından, su kanallarından, yapısal m alzem elerden vs. kaynaklanır.
G eleneksel sonlu farklar yöntem inde ağ aralığı iki gereksinim i karşılayacak şekilde seçilm elidir: (a) kalan heterojenliği gösterebilm eli (b) term al difüzyon uzunluğundan daha kısa olmalı.
Böyle bir sonlu farklar m odeli 100.000 - 1.000.000 kadar bilinm eyen içerir. Bu ise bilgisayar donanım ında ki gelişmeye rağm en ürkütücü bir problemdir.
Bunun yerine reaktör kalplerinde nötron akı dağılım ını ve etkin çoğaltm a katsayısını bulm ak için çok sayıda yaklaşım yöntem i geliştirilmiştir. Bunlar nodal, kaba ağ ve sentez yöntem leri olarak sınıflandırılır.
N odal yöntemlerde, reaktör kalbi nod denilen büyük hom ojenleştirilm iş alanlara bölünür. Genellikle bir yakıt topluluğu (asem ble) ya da toplulukları bir nod olarak tanımlanır. Böylece bilgisayar zam anından ve depolam a alanından kazanılır. Nodal hesaplam alar sonucu b ir yakıt topluluğu için güç ya da ortalam a akı ve reaktör için etkin çoğaltm a katsayısı bulunur.
N odal yöntem lerin tem el fikri iki nod arasındaki yüzeyde nötron akımları ve bu nodlarda ortalam a nötron akılan arasında ilişki kurmaktır. Bu ilişkiyi sağlayan bir katsayı matrisi oluşturulur.
G eleneksel ve dik yönde integre edilmiş nodal yöntem ler olm ak üzere birbirinden oldukça farklı iki sınıf nodal yöntem geliştirilmiştir. H er ikisi de aynı nodal denge denklem ini kullanm alarına rağmen ayrık sistemi çözecek ek denklem leri farklı şekilde türetirler. Bu tezin teorik tem elini dik yönde integrasyon yaparak elde edilen nodal açılım yöntem i oluşturur.
Bu çalışm ada polinom açılım yöntem lerinden biri olan nodal açılım yöntem lerinden en düşük dereceden olanı kullanılm ıştır. Sistem geometrisi olarak bir boyutlu silindir alınmıştır. Açılım katsayılarının bulunm asında Fick Yasasından, ayrık nodal denge denklem inden ve normal akım ın sürekliliğinden yararlanılm ıştır. H er bir nod için ikisi Fick Yasasından biri ayrık nodal denge denklem inden olm ak üzere üç denklem ya da b ir başka ifadeyle üç vektör elde edilmiştir. Bu denklem ler bir katsayı matrisini oluştururlar.
Çok gruplu difüzyon teorisi için yetkinlik-özdeğer hesaplam aları yapabilen bir bilgisayar programı bu matris form undan yararlanılarak geliştirilmiştir. Bu program FO RTRAN 90 dilinde yazılmış ve W IN DO W S işletim sistem inde koşulmuştur. Derleyici olarak FO RTRAN Pow er Station 4.0 kullanılmıştır. Bu program çok gruplu nötron difüzyon denklem ini çok bölgeli bir sistem için çözerek etkin çoğaltma katsayısını, akı ve akım dağılım ını ve ortalam a akımları bulm a yeteneğine sahiptir. Bu FO RTRAN program ının ismi olarak, R yönünde nodal açılım yöntemi kelim elerinin İngilizce baş harflerinden oluşan N EM R seçilmiştir.
NEM R program ını doğrulam ak için b ir gruplu, b ir grup iki bölgeli, iki gruplu problem lerin analitik çözümleri bulunm uş, bu sonuçlar hem NEM R program ının sonuçları ile hem de lineer ve kuadratik sonlu elem anlar yöntemi ile karşılaştırılmıştır. Sonlu elem anlar yöntem i için QFEM R programı kullanılmıştır. Son olarak iki grup çok bölgeli b ir reaktör olan TRIG A reaktörü için program koşulm uş ve bütün bu problem lerde N EM R program ının tutarlı ve doğru sonuçlar verdiği gözlenmiştir.
Bu tezin amacı sonsuz silindirik bir ortam için nodal yöntem programı geliştirm ek ve nodal yöntem ler ile sonlu elem anlar yöntem ini karşılaştırmak, hangi yöntem in hangi durum larda daha iyi sonuç verdiğini gözlem lem ek olmuştur. Test problem lerinden görüleceği gibi bilgisayar programı doğrulanm ış ve geliştirilen nodal yöntem in sonlu elem anlar yöntem lerine göre nod sayısının oldukça az olduğu kaba ağlarda daha iyi sonuçlar verdiği görülmüştür.
DEVELOPMENT OF A NODAL METHOD FOR THE SOLUTION OF THE NEUTRON DIFFUSION EQUATION IN CYLINDRICAL GEOMETRY
Diffusion theory is sufficiently accurate to provide a quantitative understanding of m any physics features o f nuclear reactors. A hom ogenized m ixture with average cross-sections and diffusion coefficients are used to create a com putational m odel for w hich diffusion theory is valid.
Even after assem bly hom ogenization, the reactor core remains a highly heterogeneous m edium because of assembly-to assem bly variation in fuel com position, burnable poisons, control rods, w ater channels, and structure and so on. The m esh spacing in a conventional few-group finite difference m odel of the core is constrained by two requirem ents: (a) Representation o f the rem aining heterogeneity (b) it m ust be shorter than therm al diffusion length.
Such finite difference m odel requires 100,000 to 1,000,000 unknowns. This is a form idable problem even today despite trem endous advances in com puter hardware. A large num ber o f approxim ation m ethods has been developed to enable a more com putationally tractable solution for the effective m ultiplication constant and neutron flux distribution in reactor cores. These m ethods can be classified as nodal, coarse-m esh and synthesis methods.
In nodal methods, the reactor core is subdivided into large hom ogenized regions and each such region constitutes a node. Usually a subassem bly (or sometimes a group of assem blies) is defined as a node. N odal calculations are carried out to determ ine the effective m ultiplication factor and assembly powers (or assem bly average fluxes). The essential idea of nodal methods is to relate neutron currents across an interface betw een two nodes to the average flux levels in those two nodes; that o f the coefficient m atrix is to relate the fluxes and currents on the nodal surfaces directly to each other.
Two rather distinct classes of nodal m ethods have evolved; conventional and transverse integrated. They have a com m on foundation in the discrete nodal balance equation and are characterized by the techniques they em ploy to derive the additional equations necessary to solve the discrete system. Transverse integrated nodal expansion m ethod constitutes the theoretical basis of this thesis.
In this work, a lowest order nodal expansion m ethod has been developed in one dim ensional cylindrical geometry. The expansion coefficients are determ ined by applying F ick’s law in com bination with discrete nodal balance equation and continuity of normal current. Three equations for each node are obtained. Nodal balance equation constitutes one equation and other two equations are derived from F ick’s law. These equations are used to form a coefficient matrix.
written in FO RTRAN 90 and run in W INDOW S operating system. It is com piled in FO RTRAN Pow er Station Version 4.0. It is capable of calculating effective m ultiplication factor, flux and current distribution, average fluxes for m ulti-group and m ulti-region problems. It has been nam ed N EM R (Nodal Expansion M ethod in R direction)
In order to validate this program , it was run for the problem s with know n analytical solutions. Results were com pared with calculated values and finite elem ent m ethod solutions which were obtained using a FO RTR A N program called QFEM R. (Quadratic Finite Element M ethod in R direction)
Four problem s extend from a sim ple bare, one group reactor to tw o-group, seven- region TRIGA reactor were cosidered. Effective m ultiplication factors, flux distributions and average fluxes were com pared with analytical solutions.
The objective of this thesis has been to develop a nodal m ethod program for infinite cylinder and com pare nodal and finite elem ents m ethods. It has been seen from the benchm arking problem s that this aim has been accomplished. N odal expansion m ethod and quadratic finite elem ent m ethod were shown to be of better accuracy with respect to linear finite elem ent method. It also appears that nodal expansion m ethod is a practical m ethod for the problems in w hich the m esh is very coarse.
1 INTRODUCTION 1.1 Diffusion Theory
In order to design a nuclear reactor properly, it is necessary to be able to predict how the neutrons will be distributed throughout the system. The approxim ate value of the neutron distribution can be found by solving the diffusion equation, essentially the same equation as is used to describe diffusion phenom ena in other branches of engineering. This procedure was used for the design of m ost early reactors and it is still widely used to provide first estimates o f reactor properties .
D iffusion theory provides a strictly valid m athem atical description of the neutron flux when the assumptions m ade in its derivation - absorption m uch less likely than scattering, smooth spatial variation of the neutron distribution, isotropic or linearly anisotropic scattering - are satisfied. The first condition is satisfied for m ost of the m oderating and structural m aterials found in a nuclear reactor, but not for the fuel and control elements. The second condition is satisfied a few m ean free paths away from the boundary o f large (relative to the m ean free path) hom ogeneous m edia with relatively uniform source distributions. The third condition is satisfied for scattering from m ost nuclei.
A m odern nuclear reactor consists of thousands of small elem ents, m any of them highly absorbing with dim ensions on the order of a few m ean free paths or less. Y et diffusion theory is widely used in nuclear reactor analysis and makes accurate predictions. The many small elem ents in a large region are replaced by a hom ogenized m ixture with effective averaged cross sections and diffusion coefficients, thus creating a com putational m odel for w hich diffusion theory is valid. Even after the local fuel pin, clad, coolant, and so on, heterogeneity is replaced by a hom ogenized representation; a reactor core rem ains a highly heterogeneous m edium because of the intra-assem bly and assem bly-to-assem bly variation in fuel composition, burnable poisons, control rods, w ater channels, structure and so on .
coolant densities due to nonuniform core power densities and tem perature distributions as well. Such complexities im m ediately force one to discard analytical m ethods in favor o f a direct num erical solution o f the diffusion equation. In fact even w hen an analytical solution of the diffusion equation is possible, it is frequently more convenient to bypass this in favor o f a num erical solution, particularly w hen the analytical solution may involve numerous functions that have to be evaluated num erically in any event, or w hen param eter studies are required that may involve a great m any such solutions .
1.2 Numerical Methods for Solving the Neutron Diffusion Equation
The general procedure is to rew rite the differential diffusion equation in finite difference form and then solve the resulting system of difference equations on a digital computer. The m esh spacing in a conventional few group finite-difference m odel of a reactor core is constrained by two requirem ents:
• It m ust be sufficiently fine to represent the rem aining spatial heterogeneity
• It m ust be no larger than the shortest (therm al) group diffusion length in order to avoid num erical inaccuracy
A few group finite difference m odel that could adequately describe such a core might w ell have 1 0 0,0 0 0-1,0 0 0 , 0 0 0 unknowns (the fluxes in each group at each m esh
point). The direct solution of such a problem, even in diffusion theory, rem ains a form idable computation. For calculations such as fuel burnup or transient analysis, in w hich many full-core spatial solutions are needed, direct few group finite-difference solutions rem ain impractical.
A large num ber of approxim ation methods have been developed to enable a more com putationally tractable solution for the effective m ultiplication constant and neutron flux distribution in nuclear reactor cores. Follow ing historical precedent, these methods can generally be classified as nodal, coarse-m esh and synthesis methods, although distinction among categories may be largely a m atter of perspective and sequencing o f calculational steps .
1.3 Nodal Methods
In order to avoid the large storage and execution tim e requirem ents of a direct finite difference treatm ent o f the diffusion equation a schem e is provided by so called nodal methods. The general idea is to decom pose the reactor core into relatively large subregions or node cells in w hich the m aterial com position and flux are assum ed uniform . One then attempts to determ ine the coupling coefficients characterizing node cell leakage and then to determ ine the node cell fluxes them selves .
If the num ber o f nodal cells N is large, the nodal m ethod becomes equivalent to the finite difference scheme and hence loses any calculational advantages. The real power o f the nodal approach is realized only w hen the num ber o f node cells N is sm all, since then the cells are large enough that they becom e coupled via neutron diffusion only to nearby cells, that is the transfer m atrix is sparse .
Usually a subassem bly (or sometimes a group of assem blies) is defined as a node so typical dim ension of a node cell is in the order of ~20 cm in a nuclear reactor. Nodal calculations are carried out to determ ine the effective m ultiplication factor and assem bly powers (or assem bly average fluxes). N odal m ethods are classified further into two distinct categories: conventional nodal m ethods and transverse integrated nodal methods.
Conventional nodal methods are based on the calculation of a nodal average flux or fission rate for each hom ogenized subassem bly. N eutron diffusion betw een adjacent nodes is represented by coupling coefficients w hich are usually determ ined by em pirical means. Conventional nodal m ethods are also called sim ulators and can achieve im pressive accuracy for the particular reactor type they were intended. Due to the lack of theoretical foundation in conventional nodal methods, researchers have sought to develop nodal methods w hich have sound theoretical basis and w hich converge to the correct solution as the node sizes are taken smaller. The transverse integrated nodal m ethod is based on integrating the three (or tw o) dim ensional diffusion equation over the transverse direction(s) and reducing it to a one dim ensional diffusion equation, with transverse leakage term s. The m ost popular of the transverse integrated nodal methods is based on a polynomial expansion of the one dim ensional flux and can obtain higher orders of accuracy than conventional
1.4 Finite Element Method
The finite elem ent m ethod (FEM) provides a system atic approach for developing solutions with higher order accuracy than the conventional finite difference equations. Although FEM could be based on a weighed residual approach, FEM applications to the neutron diffusion depend usually on the equivalence of the solution o f the diffusion equation with the m inim ization of a functional. System of interest is divided into hom ogeneous regions o f specified geometric shape which are called finite elements. The spatial dependence of the neutron flux (and sometimes current) is assum ed to be a polynom ial of a certain degree within the element. The degree of the polynom ial dependence determines the nam e of a particular FEM application (i.e. linear, quadratic, cubic etc. finite elements) .
FEM is a coarse-m esh method. Like nodal methods, coarse-m esh m ethods generally require detailed regional heterogeneous flux distributions in order to construct hom ogenized param eters and to com bine with the coarse-m esh solution to construct a detailed heterogeneous flux solution .
1.5 Objectives of the Work
The objective o f this thesis is to develop a nodal m ethod which gives acceptably accurate results for realistic problem s at a considerable savings over finite difference methods. This m ethod is capable of calculating m ultiregional and m ultigroup problem s in a one dim ensional cylindrical geometry. Com parison this m ethod with FEM is another aim of this work. These are accom plished by using a com puter program written in FO RTRAN 90. It is nam ed N EM R as the abbreviation of nodal expansion m ethod in r dimension for the cylinder.
Chapter 2 is about the developm ent of a transverse integrated nodal method. A nodal form alism is obtained using lowest order nodal expansion method. The expansion coefficients are determ ined by applying Fick’s law in com bination with continuity of normal current. Three equations for each node are obtained. N odal balance equation constitutes one equation and other two equations are derived from F ick’s law. These equations are used to form a m atrix equation. Then, a com puter program which can carry out criticality-eigenvalue calculations in m ultigroup diffusion theory is developed using this m atrix form. This program can find effective m ultiplication
factor, flux and current distribution and average fluxes for m ultigroup and m ultiregion problems.
Chapter 3 is devoted to validation of this program. It is run for the problems which have been solved analytically. Results are com pared with analytical and FEM solutions. Linear and quadratic finite elem ent m ethod results are obtained using a FO RTRAN 90 program called QFEM R. (Quadratic Finite Elem ent M ethod in R dim ension)
Chapter 4 includes the conclusion and recom mendations for future work. Finally, appendixes are given for the description of N EM R program.
2 NODAL FORMALISM IN CYLINDRICAL GEOMETRY 2.1 Cell and Edge Averaged Quantities
The m ultigroup neutron balance equation
^ V /Tg W , ^
g=1 k g=1
w ith no group-to-group upscatter assum ption, and Fick’s Law Jg (r ) = - D g ( r ) V^g ( r )
constitute the basis of the nodal formalism.
To avoid com plexity in notation, scattering and fission sources will be defined as
Since group index is of no im m ediate concern in nodal developm ent, the group indices will be suppressed and (2.4) and (2.2) will be written as
A lthough nodal m ethods are generally devised for three dim ensional whole core calculations, the transverse integration procedure fails in (r,0) or (r,0,z) cylindrical
geometry. Because the transverse integration over r (in 2-d) or z and r (in 3-d) leads to an impasse. The difficulty arises w ith the azim uthal term . This im passe may be circum vented using analytical methods but it is not within the scope of this work. Here a one dim ensional developm ent is presented. However, extension from 1 -d (r) to 2-d (r-z) is actually trivial for further studies.
g G X uZg'(r > g ( r ) , g = 1,2,l ,g g=1 (2.3) Thus, (2.1) becomes V -Jg (r ) + Zg ( r ) ^ (r ) = Qg ( r ) (2.4) V - J (r ) + £ r ( r ) <K? ) = Q (?) (2.5) J (r ) = - D ( r ) V ^ (r) (2.6)
Figure 2.1 illustrates a cylindrical m esh having nodes Ai=(ri_i/2, ri+i/2). The radii at node centers are similarly denoted as (ri) and it is convenient to define m esh spacing Ari=ri+1/2-ri_1/2. 0 A2 ♦ --- ► ri —I- - ^ R / *-r1/2 r3/2 r5/2 ri-1/2 ri+i/2 rN+1/2
Figure 2.1 N odal M esh Im posed on One-dim ensional Cylindrical Dom ain Ari is replaced w ith Ai for simplicity
A i _ ri +1/2 ri - 1/2 (2.7)
Inherited from the nodal perspective and com m on to all nodal discretizations is the choice of cell and edge-based unknowns. The cell-based unknowns are defined by
$i = Qi = ( r 2 — r 2 ) VAi+1/2 Ai-1/2 ' r. +1/2
j^ (r )rdr i -1/2 ri+1/2 ( r 2 - r 2 ) VAi+1/2 Ai-1/2 '
r=r-jQ(r)rdr r=ri-1/2 2 2 (2.8) (2.9)
w hich are the cell average flux and source. W hile the edged-based unknowns, namely edge average fluxes and currents are just point fluxes and currents at node boundaries in one dimension, because there is no other dim ension over which to find averages. They are shown as ^ i+1/2, ^ i-1/2 and Ji+1/2, Ji-1/2
In some nodal methods, edge-averaged partial currents are also needed. To avoid notational complexity, the partial currents will be denoted by the lower case letter j. j++1/2 and j——1/2 denote the outgoing partial currents, while j—+1/2 and j+-1/2are the
incom ing partial currents at the right (i+1/2) and left (i-1/2) edges respectively.
U nder P1 approxim ation
j ± (r ) fr(r) ± n u •J ( r )
4 2 (
where u is an arbitrary direction. If diffusion theory is valid,
j± (r )= 1 * (r) + D(r ) d ^ ( r ) JuW 4 VW 2 3u (2.11) from (2.10), obviously: J i+1/2 = j i+1/2 — j i+1/2 (2.12) J i-1/2 = j i—1/2 — j i—1/2 (2.13) and ^i+1/2 = 2 (ji+1/2 + j i+1/2 ) (2.14)
^i—1/2 = 2(ji—1/2 + j i—1/2 ) (2.15) 2.2 Nodal Balance Equation
The nodal view of the first order system (2.5) and (2.6) suggests that a natural starting point is to integrate the exact balance equation (2.5) over an arbitrary cell
JV- J (r ) dV +
JI r (r) $ (r) dV =
JQ (r) dV (2.16) V- m V- m V.n In cylindrical coordinates, i+1/2 2
n J1 — [rJ(r)]rdr
n JQ(r)rdr i+1/2 i+1/2 *1—1/2 r dr ri-1/2 ri-1/2 (2.17)
Since the node is already hom ogenized, m acroscopic rem oval cross section is constant for a node, I r(r) = I r and it is also assum ed constant for all nodes of the same material, I r = I r . W hen the integrations are carried out with the help of (2.8) and (2.9), (2.17) becomes:
2n[ri+1/2J i+1/2 ' ri—1/2J i—1/2 ]
]+ I rfoi2n(ri+1/2 ri-1/2 ) = Q i 2n(riri2+1/2 ri21/2 )
where and Q ; represent the node averaged flux and source respectively. Surface areas can be defined as
S i+1/2 = 2nri+1/2 , and Si—1/2 = 2nri—1/2 (2. 19)
E ^ .A. , x Q.A. , x
Si+1/2Ji+1/2 - Si-1 /2Ji-1/2 + 1 2 1 (Si+1/2 + Si " 1/2 ^ = ~ 2 ± (Si+1/2 + Si_1/2 )
Surface areas at nodes are
s. = 2nri
where r. can be written as
(2.2 0) (2.2 1) ri+1/2 + ri-1/2 ri = 2 (2.2 2) Using (2.21) and (2.22) S + S S = ° i+1/2 T ° i-1/2 2 (2.23) Substituting (2.13), (2.14) and (2.23) into (2.20), finally, the discrete nodal balance equation is
Si+1/2 (ji+1/2 j i+1/2 ) Si—1/2 (ji—1/2 j i-1/2 ) + E r ^ i AiSi = Q iAiSi (2.24)
Various nodal m ethods have a com m on foundation in the discrete nodal balance equation and are characterized by the techniques they em ploy to derive the additional equations necessary to solve the discrete system.
Two rather distinct classes of nodal m ethods have evolved. The first class, often referred to as conventional or sim ulation models, makes use o f detailed calculations or reactor operating experience to evaluate the edge averaged currents in terms of differences in cell averaged fluxes for adjacent nodes, w ith empirically adjusted coupling coefficients. So these methods lack theoretical foundation.
The second class, sometimes referred to as consistently form ulated m odels, makes use of the concept of transverse integration and of higher order (than ordinary finite difference) approxim ations to evaluate the edge averaged currents and the internodal coupling terms in order to derive nodal equations that can be expected to converge to the exact solution in the lim it of small m esh spacing.
2.3 Transverse Integration
Finneman, et al. developed a popular discretization procedure which utilizes the partial current directly. This is transverse integration procedure w hich has becom e a cornerstone o f m odern nodal methods .
The usual strategy for solving the neutron diffusion equation in two or three dimensions by nodal methods is to reduce the m ultidim ensional partial differential equation to a set of ordinary differential equations (O D E’s) in the separate spatial coordinates. This reduction is accom plished by transverse integration o f the equations. In cylindrical coordinates, the three-dim ensional equation is first integrated over r and 0 to obtain an ODE in z, then over r and z to obtain an ODE in 0, and finally over 0 and z to obtain an OD E in r. Then these O D E ’s are solved to obtain one-dim ensional solutions for the neutron fluxes averaged over the other two dim ensions. Because the solution in each node is an exact analytical solution, the nodes can be m uch larger than the m esh elem ents used in finite-difference solutions. Then the solutions in the different nodes are coupled by applying interface conditions, ultim ately fixing the solutions to the external boundary condition .
Transverse integrated equations contain transverse leakage terms. They m ust be satisfied in an integral sense. That is, these equations are m ultiplied with the weight function and integrated over the node. The integrals o f the weighted residue m ust vanish. This is called as weighted residual procedure. In the lowest order nodal expansion m ethod (NEM), weight function is chosen as one. This procedure gives (2.24), discrete nodal balance equation. In the low est order case which is treated here, (2.24) is enough to constitute one equation for a node.
There are two distinct applications of the transverse integrated equations. The first is the weighted residual procedure o f the N EM by Finnem an et al. In one dim ensional analysis here, transverse integration is not necessary to reduce dimension. The second variety depends on the analytical solutions o f the transverse integrated equations. The second approach has been presented in  for cylindrical geometry
2.4 Nodal Expansion Method
The developm ent of m odem consistent nodal discretizations began in the m id 7 0 ’s. These methods were based on local polynom ial expansions. The first polynomial m ethod was the nodal expansion m ethod (NEM). In fact, although some variations and im provem ents have been considered, the NEM ideology still dominates the polynom ial class of nodal methods.
In this low est order form, NEM considers a quadratic expansion of the averaged flux on each cell. The expansion coefficients are determ ined by applying F ick’s law in com bination with discrete nodal balance equation (2.24) and continuity o f normal current.
Considerable effort has been m ade to utilize higher order polynom ial expansion within NEM . The difficulty this creates is centered around the evaluation of the higher order expansion coefficients. In particular, the weighted residual procedure that is typically used relies on transverse-integrated equations and as a result an approxim ation of the transverse normal currents (i.e. transverse leakage) is also required.
2.4.1 Construction of Polynomial Basis
The N EM treatm ent of the transverse integrated O D E’s is based on a low order polynom ial expansion of the transverse integrated flux. In one dim ension this is ju st r dependent flux
Wr ) = Z ^ l^ ri-1/2 < r < ri+1/2 (2.25)
A t this point a local variable £ is defined as r - r
4 = — L (2.26)
where £=±1 / 2 w hen r=ri±1/2.
£ and r are the same order. Then, (2.25) can be written in local variable N
<K4) = Z a i P (4), - 1 / 2 < £ < 1 / 2 (2.27)
For simplicity, the first polynom ial is chosen as Po© = 1
and the higher order polynomials are required to be orthogonal to P0(£)
( 0 d£, = 0 ,
Transform ation of integration operator gives dr d Ş = s : Using (2.30) and (2.26) in (2.8) 2 (2.30) = ( r 2 - r 2 ) V*i+1/2 *1-1/2 2 -1/2
f ^ ) ( ^ 4 + ri)A id^(2.31) Substituting (2.22) into (2.31) 1 1/2 fc = A
- f+ p)Aid^ (2.32) Airi -1/2
Low est order NEM uses quadratic expansion, so N=2. Inserting (2.27) into (2.32) i 1/2 N=2 fc =T —
SarPl(^))(^Ai + p) ^ (2.33) Airi -1/2 l= 0 Thus (2.34) -I 1/2
fc = - f(a0P0
+ri -1/2 Using (2.28) 1 1/2
fc = - f(a0
^ ^ +a2P2
The orthogonality requirem ent results in last two terms to be vanished in the integrand. Also integration of the first term is zero. Thus
fc = 1 1/2 1/2 1/2 l1
© ^ +a2
© ^ +a0ri
^-1/2 -1/2 -1/2 (2.36) -1 /2 r
All odd polynomials satisfy the requirem ent (2.29). Thus
Pife) = a £ (2.37)
where a1 is a constant yet to be determined. If an even polynom ial has only one term
(one monomial), it can never satisfy (2.29), since
f4 nd£, = .r
--f/2 (n +1)2' Specifically 1/2
1«2d« = Tî ■ , n even (2.38) -1/2 12 (2.39)
Thus any quadratic polynom ial of the form
P 2 fe) = « 2 f ^ - ^
w ould satisfy (2.29).
Substituting (2.37) and (2.40) into (2.36)
(2.40) = 1 1/2 1/2 1/2 y 1/2 a
fÇ2^ + a2A itt2
fÇMÇ- a2A ia2
fd^ + aori
fd^ -1/2 -1/2 -1/21 2 -1/2 (2.41)
Second and third integrals have to be vanished since they contain odd polynomials
fci =1 a 1A i a ! 1 2 + a 0ri = a o +
If 4=+1/2, this describes node outer boundary
fci+1/2 = £ a ,Pt (1/2) = aoPo(1/2) + a1P1(1/2) + a2P2(1/2)
, a . a2
fci+1/2 = a0 + a ^ T + a ^ ^
If £=-1/2, this describes node inner boundary
(2.43) (2.44) fci-1/2 = £ a ıPi ( -1/2) = a0P0 ( -1/ 2) + a1P1( - 1 /2) + a2P2( -1/ 2) l= 0 a a (2.45) r
Subtracting (2.46) from (2.44) a i — foi +1/2 Ai-1/2 ) a! Adding (2.46) and (2.44), ^i+1/2 + A i—1/2 2 a 0 + a 2 O i 3 Substituting (2.47) into (2.42), A i (A i+1/2 — Ai—1/2 ) Ai = a0 +■ Thus a0 = ^ i +
12rSubstituting (2.50) into (2.48), (2.47) (2.48) (2.49) (2.50) A i+1/2 + A i—1/2 - 2^ i + ' A i A i—1/2 A i^ii+1/2 6r. 6r. + a2
a 2(2 .5 i)
The last coefficient a2 is found from (2.51)
( ( a 6r + A f 6r f A i+1/2 2^i + i J 6ri - A i V 6ri J ■A A A i—1/ (2.52)
To m ake (2.47) and (2.52) as simple as possible, a 1=1 and a 2=3 are chosen. Polynom ial basis are constructed for the lowest order N EM as
3 a 2
p, (^) - 3 42 —4
Using these polynomials and coefficients in (2.27),
12r1 + K^i+1/2 ^i—1/2 )]3
+6r + A . v 6ri J ^i+ 1/2 — 2^i + ^ 6ri — A i 'a A V 6ri j Ti—1/2 342— 1 4 (2.53) (2.54) (2.55)
Finally, after some arrangements in (2.55) W
4) _ 3 _ 652) ' + - 1 _ 1 v8(ri / Ai) 4+
4+ 3 + -1___
7i+1/2 2 4 + r _ j _____1 v 8(ri / Ai) 4
_ 4+ (
3- - 1 A
72(ri / Ai)
^i-2 4 1/2
2.4.2 Fick’s Law F ick’s law states that
J(r) = -Dd^(r) dr Note that d^(r) _ d^ d4 _ 1 d ^ (r ) dr d4 dr Ai d4 (2.57) m ay be written as J(
4) _ -d
, M )Ai d
4(2.56) (2.57) (2.58) (2.59) Using (2.56) d<K
4) _ _ ( ( 1 2^^i + 1 + 2 4 d 4 l 3 + - 1 ■U 2(ri / Ai)
Substituting (2.60) into (2.59), for 4=+1/2
^i+1/2 + _ 1 + 24 ( 3 _ - 1 2(ri / Ai)
J_ _ D l i+1/2 A. ( _ 6 ^ i + 4 + - 1 l 2(ri / A i) ( ^i+1/2 + 2 _ -1 l 2(ri / Ai) ^ i—1/2 Similarly, for 4=-1/2 j _ _
Di_1/2 A . ( ( 6 ^ i _ l 2 + 1
2(ri / Ai) ^i.
2(ri / A i) ^i_1/2
(2.61) may be written in term s of edge averaged partial currents, using (2.14) and (2.15)
J = « D A _ 2 D 4 + ! '][;+ j- ] Ai A i V 2( ri/ Ai) )
A. 2 2(r. / A i)i a \ | [ji—1/2 + Ji—1/2 ]
Using (2.12) in (2.63) and after arrangements
^ 8D. d. N 1 + ----L + — V A i ri )ji+1/2 \ —8 D —D ı N V Ai ri )ji+1/2 — ( 4D. i iD N V Ai ri ) ji —1/2 ( 4 D i iD N V Ai ri ) — 6D. Ji—1/2 + . ^i — A
Similarly, using (2.13), (2.14), and (2.15), in (2.62) + •— 6D iA 2D. Ji—1/2 Ji—1/2 _ A ^i + A A. A. ( 2 + 1 1 V 2(r. / A.) +2Di ( Ai 4 - 1 \ V 2( r. / A.) + 2Di Ji—i-1/2 1/2 + .Ai ( 4 - 1 V 2(r. / A.)_ Ji+1/2 + Ji -1/ 2D ( Ai 2 + 1 V 2(r. / A.) (2.64) (2.65)
A fter some arrangements N Ji—1/2 \ + 8 D _ D N V A i ri ) ( 4D + D iN V Ai ri )Ji+1/2 + 4 D + D V Ai r. Ji+1/2 ( 1—8 D + Dl Ai r 1 J+ 6Di Ji—1/2 . $i A (2.66)
(2.64) and (2.66) are the expressions for the outgoing partial currents. Substituting the outgoing current j “ 1/2 defined by (2.66) into (2.64) gives
^ 8Dj d. ^ 1 + ---- L + — v A i rı JJi+1/2 _ 8D. D ^ ( A. rı JJi+1/2 —v
4Di Di A. r Ji—1/2 ( 4D i iD V + A; r 4D 1 I D1 v i i jv i + Ai ri ( 4D 1 D V1 Ai ri 4D 1 _l_ D1 ( 1+ 4 D i - D v Ai r Y D ^ Ji+1/2 + ( 1 + rı J v 8D i D ; ^ ( 4D; 1 + 1 1 A. ri J + v Aı v i i jv i + Ai ri A 8Di Di Y - Y Ji+1/2 6Di^i (, 8D. d. ^ Ji-1/2 - ( 8D d. ^ Ai 1 + i V Aı ri J v Ai rı J 6D ; . + — ^ ^ Ai Collecting terms Ji+1/2 ^ 8 D; d . ^ 1 + ---- L + —-v A i ri J 16D
2d Ai2 r i2 ^ 8D ; D ; ^ 1 + ---- L + —-v A i ri J ' Ji—1/2 ( 4 D D; 4 D i iD Y ' i —__^ + V i A; r; A r i JV 1—8Dl+ Dl Ai rı J ( 1 + 8Di — Di > v A i r j + Ji+1/2 (1 — 8D; — D. > v Ai ri +
"ÂT T"i 1 +8Di i iD A ı r. J + ^ı 6Di( 4D ; D ; ^ 6D; \ Aj r; j A; . ( . 8D ; D ^ A 1 + A ı ı ri J 1 Simplifying (2.68) yields Ji+1/2 1 1 + ---L + - 16D; 48D
2' Ai Ai2 ' j i—1/2 2D. r A ı + j i+1/2 2D r 1 4 8 D 2 a2 + ^ ı 6D Ai ^ + 2 4 D 2 a2
Sim ilarly, substituting the outgoing current Ji++1/2 defined by (2.64) into (2.66)
( Q oc D y +
!v Ai ri y 4 D i D Y Ji—1/2 v A i r 1 —8Di Di i Ai ri D 8D i D ^ 1 + — L + Ji+1/2 v A ri y f
2^ v Ai2 ri2 f
2d2 y i f , 8D i D ' ^ ^ 1 + — L + — v A i ri y v Ai2
4D + d, Ai r 6Di A 1 + — L + A r D i ^Ji—1/2 + oc Q D i y 1 + i + i ri y v Ai ri y D *i (2.70) + f 4D i D ^ v A i ri y Ji+1'2 f1—8 D + D i y A r + 6Di Ji—1/2 + A ^i Collecting terms in (2.70) Ji—1/2 f
1 6d2 d2 v A2 ri 2 y f oc Q D y f OC cf D i y 1 + 1 + + i v Ai ri y v Ai ri y f oc cf D iy 1 + + —^ v Ai ri y Ji+1/2 f
4D i iD ^ + v Ai ri y D 8Di Di ^ v Ai ri y Ji—1/2 f _ 8 D + d, yf v Ai ri yv 8Di Di 1 + — L + — - A r i y +
2_ Dl A
2r 2 D 8D i D i ^ v Ai ri y (2.71) + ^i 6D 4D Ai v Ai 1 y f 8D i D iy 1 + i + v Ai ri y Simplifying (2.71) Ji—1/2 1 16Di 48D2 1 + ---L + - i + J+i—1/2 Ai Ai2 1 2D i 48D2 1 + — i - i
_ Ji+1/2 2Di 8Di ri Ai Ai2 + ^i 6D i 24D i2 - + - i Ai Ai2 (2.72)
If a dimensionless variable is defined as 2 r
( 16D. 48D? 1 1 + - v A + ■ A 22i y (2.69) becomes ( 1 4 1 2D. Ji+1/2 v ri A i y ;+ T: J i—1/2 + ( 2d. 48D2 1 v xi A i y ; 6 D A T Ji+1/2 + i v 1+4 D Ai y T and (2.72) becomes - 2D j i—1/2 _ ' ( 1 4 1 —+ — v ri Ai y ( , 2D, 48D 2 1 1 + — ’ T J i+1/2 + ' A2i y T Ji—1/2 + 6Di A i Ti ı + 4 D i v Ai y T = r (2.73) (2.74) (2.75)
(2.74) and (2.75) constitute 2 equations per node. The num ber of unknowns per node is three. The outgoing partial currents ( j++1/2, j-_1/2 ) and the cell flux ^ constitute the three unknowns. The incom ing partial currents ( j+—1/2, j -+1/2) can be considered known quantities, since they are either equal to the outgoing partial currents o f the neighboring cells or are known from boundary condition.
The remarks above are valid if conventional hom ogenization theory is used and, thus continuity o f partial currents is assumed. If equivalence hom ogenization theory is used; the incom ing partial currents are not equal to the outgoing partial currents of the neighboring cells. The incom ing partial current can be written in term s of the outgoing partial current o f the neighboring cell, outgoing partial current of the same cell and the flux discontinuity factors.
If we w ish to use a higher order flux expansion, say N=4, the higher order flux m om ents appear in their counterpart of the two equations (2.74) and (2.75). Thus, extra equations are needed. W eighted residual procedure equates the num ber of unknowns to the num ber of equations provided the treatm ent o f transverse leakage term does not introduce any extra unknowns. The m ost successful approxim ation for the treatm ent o f the transverse leakage term has been an approxim ation called quadratic approxim ation .
(2.74) and (2.75) can be written in shorter forms, if the following variables are defined
2.5 Iterative Solution of the Nodal Equations
2Di m. r
14 2 V ri a i
2Ti r _ 2D i 48D 2 ^ r A2
26D i Ai 1 + 4D. pi = V A i 2 Ti T ( 1+ 2D i 48D 2 7 V A 2 o i = T t i = _ 2Dir 1 .v i T i 4 1 + — A i 2 Therefore (2.75) becomes j i + 1 / 2 = mi ji _ 1 / 2 + niJ( + 1 /2 + Pi^i j__1/2 = ti Ji+1/2 + oij+_1/2 + Pi 6i (2.76) (2.77) (2.78) (2.79) (2.80) (2.81) (2.82) Here J i+1/2 and j+_1/2 are know n quantities. Butj++1/2, j i—1/2 and 6 are unknowns. So a third equation is necessary for the determ ination of the three nodal unknowns. It’s the discrete nodal balance equation. From (2.24)
^ = [S i+1/2 (ji+1/2 j i+1/2 ) Si—1/2 Ui—1/2 (j i+1/2 E r A iSi
■jr- 1 /2 )] + Qi AiSi
E r A iSi
(2.83) can be written as
(2.81) and (2.82) can be put into m atrix form as 2.5.1 Matrix Equation Ji + 1 / 2 J i —1 /2 m i n i . o i t i j i—1 /2 j i + 1 / 2 + ^ i pi pi
(2.85) can be written shortly • o u t
j = R j in + ^ p
and (2.84) can be written as
< fc = 2 - [q, — s t (jo u t— ji n)] ^ r where ST 1 a,s,[si + 1 / 2 S i - 1 / 2 ] R = m i n i - oi ti _ j+ j + • o u t J i + 1 / 2 and ji n = Ji —1 /2 J = _ Ji —1 /2 _ _ Ji + 1 / 2 _ (2.85) (2.86) (2.87) (2.88) (2.89) (2.90)
(2.86) and (2.87) are used for iterative process. One m ethod is to use two m atrix equations separately. Q, w ould be known from estimates keff and ^ of the previous outer iteration or initial guess. Inner iteration provides estimates of the incom ing partial currents and cell average fluxes. A m arching procedure takes places throughout the core in an inner iteration. Inner iteration ends when the per cent difference betw een the consecutive estimates of the cell averaged fluxes drops below a certain predeterm ined convergence criteria.
A fter the convergence of the inner iteration, next outer iteration begins and a new keff estim ate is found. The outer iterations stop when the percent error betw een keff estimates drops a certain predeterm ined convergence criterion.
Alternatively, (2.86) and (2.87) are written in one m atrix equation form as
A J = 1
Here, J is the unknow n vector. It contains outgoing partial currents and cell average fluxes. A is 3Nx3N band matrix, where N is the total num ber of nodes. F is again 3Nx3N diagonal m atrix and only (3i-1)th elem ents contain nonzero fission source term. (i= 1,2,...N ) S is the scattering source vector. In the first iteration, right hand side (RHS) of (2.91) is known. J is known from initial estimates of outgoing partial currents and cell averaged fluxes.
A J (n+1) = b (n) (2.92)
N ew J vector is found with a linear system solver w hich has two m ain subroutines, first makes LU factorization of the m atrix A and second solves the system. N ew keff estim ate is found after fission source iteration. Iteration continues until the difference betw een two successive keff estimates drops below the convergence criterion. Converged J vector contains both fluxes and currents as it is shown in (2.93) and they are separated into flux and current vectors in the next step.
J = [jl/2 $1 j 3/2 j 3/2 $2 j 5/2... j (2N-1)/2 $N j (2N+1)/2 ] (2.93) 2.5.2 One-Node Formulation
In one-group form ulation, rem oval cross section is equal to the absorption cross section. Node thickness is equal to the radius of the cylinder. Only fission source exists.
Then, (2.82), (2.84) and (2.81) become
j 1/2 = t 1 J 3/2 + o 1j 1/2 + p 1$1 (2.94) $1 = VEf $ _ S3/2 j 3/2 Eakeff $ 1 EaAS1 S1/2j 1/2 + S3/2j 3/2 + S1/2j 1/2 E a AS1 E a AS1 EaAS1 (2.95) j 3/2 = m 1 j 1/2 + n 1j 3/2 + p 1$1 (2.96) Boundary conditions are J1/2 = j+/2 _ j1/2 = 0 ^ j+/2 = j1/2 (reflective) a n d j3/2 = 0
(vacuum). A lso S1/2=0, no surface area exists at the center o f the cylinder. Thus, after the application of boundary conditions, the equations (2.94), (2.95) and (2.96) describe the whole system:
S 3/2j 3/2 _ V E f A + ^ aAS: E ak eff j 3/2 — m i j 1/2 — Pl^1 = 0 In m atricial form A lil = ^ - El il eff where 1 _ o l _ Pl 0 !> II 0 1 S °3/2 Ea ASl _ m l - P1 1 iT = [Jl/2 A1 j 3/2 ] E ı = 0 0 0 vE 0 0 E a 0 0 0 2.5.3 Two-Node Formulation
W hole system is hom ogeneous and the nodes have equal thickness. If i=1, (1 - °1 )j+/2 - Pl^l - j = 0
A 1 +
S 3/2j 3/2 S 3/2j 3/2 V E f - $ 1
EaASl EaASl Eakeff j 3/2 — m lj 1/2 — n 1J 3/2 " M l = 0
For i=2, j_/2 = 0 (boundary condition) j 3/2 " o 2 j 3/2 — p 2^2 = 0 A + S5/2j 5/2 + S3/2j 3/2 _ S3/2j 3/2 = VEf A ^ 2 E a AS2 E a AS2 E a AS2 Eakeff ^ j 5/2 m 2j 3/2 " P2A2 = 0 (2.98) (2.99) (2.100) (2.101) (2.102) (2.103) (2.104) (2.105) (2.106) (2.107) (2.108) (2.109)
(2.104) through (2.109) can be written in m atrix form A 2J2 = k L F2 J2 k eff where A2 A1 (3x3) 0 0 - o 0 0 - 1 1 0 S S3/2 0 0 0 Z a AS1 - n1 0 0 - o2 1 - P2 0 S ° 3 /2 ° 3 /2S 1 S5/2 ^a AS2 Z a AS2 Za AS: - m2 0 - P2 1 J 2 [j1/2 j 3/2 j 3/2 ^ 2 j 5/2 ]
£2= p (3x3) 0 (3x3) 0 (3x3) £ (3x3) .= =1 J 6x6
2.5.4 Three-Node Formulation and General Matrix Form Three-node m atrix form is
A J = k L F , J k eff where ( 2.110) (2.111) (2.112) (2.113) (2.114)
0(3x3) As A (6x6) t2 0 0 S5/2 0 0 (2.115) S a AS2 - n2 0 0 - o3 1 - p3 0 II S5/2 S5/2 1 ° 7 /2S S a AS3 S a AS3 Sa AS3 - m3 0 - p3 1 9 x9 jT = .J1/2 ^ 1 J3/2 J3/2 ^ 2 J5/2 j 5/2 ^3 j 7/2 ] (2.116) "-,-^(3x3) F1 0 (3x3) 0(3x3) _ £3 = 0(3x3) (3x3) £1 0(3x3) (2.117) 'ii q^ 0 (3x3) (3x3) F1 J 9 x9
The im plem entation of A m atrix in the com puter code can be accom plished by com paring A m atrices in one-node (2.101), tw o-node (2.111) and three-node (2.115) form ulations. General m atrix form
AnJn = kL Fn Jn (2.118)
An = An-1 ((3N-3) x(3N-3)) q(
4)) -,((3N - 6 ) x3)) t N-l Q Q S °N-l
/2Q Q E a ASn-1 — n N-l Q Q - o N l - Pn Q S °N-l
/2Ea ASn E a ASn E a ASn - m N Q - Pn l J N [jl/2 j 3/2 j 3/2 ... Jn-1/2 ^N Jn+1/2 ] En „(
3) El q
3) q (3x3) q (3x3) q
3) ■•Ei (
3) J3Nx3N J3Nx3N
2.5.5 General Matrix Equation with Reflective Boundary Condition Reflective boundary condition is
J N+l/2 = j N+l/2 — j N+l/2 = Q
( 2 .ll9 )
( 2 .l2 l)
(2.122) And (2.122) states that jN+l/2 is equal t o j N+l/2 at boundary. So, only the equations related w ith the last node will be changed. This corresponds to last colum n o f A , i.e. A(N,N), A (N -l,N ) and A(N-2,N). A fter some simplification, these equations become
j N-l/2 - o N j N-l/2 t Nj N+l/2 PN Q (2.122) ^ N + SN-l/2 j N-l/2 e a s.t SN-l/2 j N-l/2 e a s.t vE f E k " (2.123) (l - n N)jN+l/2 m Nj N-l/2 - P N^N _ Q (2.124) and A is
An = An-1 ((3N-3) *(3N-3)) q(
4)) H O % -6) *
3)) t N-l Q Q s °N-l
/2Q Q E a ASn-1 — n N-l Q Q - on l -pn " t 2 S °N-l
/2l Q Ea ASn E a ASn - m N Q -pn l - n N J 3Nx3N
F and J are the same as before.
2.5.6 Fission S o urce I te ra tio n in M u ltig ro u p D iffusion E q u a tio n s
The multigroup diffusion equations can be written as
V Di V ^ + E ri^ = - 1 - %!Q k. f f ■ V.D 2 V^2 + E R2^2 _ . %2Q + E S1^2^1 k eff V D G V^G + E RG^G _ , %GQ + E Sl^G$1 + ...+ E S(G-l)^G^G-l k eff
It is assumed that there is no upscattering and fission source is defined as
Q(r) = X V g 'M g '(r)
The spatial dependence of the fission source is identical in each group diffusion
The initial estimates of Q(r) and multiplication eigenvalue keff are made before first
Q(r) ~ Q(Q
keff ~ k(Q
- V.DjV^j + E R1^(1)= %ıQ(0)(r) k eff
The flux in the first group is calculated and the diffusion equation for the next lowest energy group is solved
- V .D2V ^ + £ R2^ = ^ o f X2Q(0 )(r) + ^ 1^
and flux of this group is determ ined for every nodes. All of the group fluxes are found using this procedure. Diffusion equation corresponds to neutron balance equation, so only RHS of (2.123) is changed with additional scattering source term. A new fission source can be calculated since ^ (1), ^ , ... ^ are known.
Q (1)(r) = X v g-E fg 'C (r) g'=1
And a new value of keff:
k(1) _ eff J Q(1 )(r)d3r ^ 1 QI0 I(r)d3r keff (2.131) (2.132)
There will always exist a m axim um eigenvalue, keff that is real and positive. The corresponding eigenfunction, cell average fluxes and outgoing partial currents, is unique and nonnegative everywhere within the reactor. It can be dem onstrated that this fission source iteration will converge to this positive dom inant eigenvalue keff and the corresponding eigenfunction .
At this point one tests the source iteration for convergence, such as by comparing keff(n+1) - k(n)eff
k (n+1) ?
< e (2.133)
If the changes in k (enff) are sufficiently small, one assum es that convergence has been achieved, and the iterative procedure is ended. If not, a new fission source is calculated and the iteration continues.
3 N U M E R IC A L A P P L IC A T IO N S
The form ulations derived in the previous chapter have been im plem ented in the N EM R code which is a com puter program written in FO RTRAN 90. This code and the num erical results obtained from it will be described in this chapter. Also com parison with FEM is m ade by using the results of a com puter program, QFEMR.
3.1 O n e-g ro u p , B are, H om ogeneous R e a c to r
In this problem, a bare, cylindrical reactor of diam eter 7.5cm is considered (Figure 3.1). Zero incom ing current boundary condition ( j_ = 0 ) is assum ed at the surface of this cylinder. Effective m ultiplication factor of this system is determ ined using the one group cross sections D=0.65cm, Ea=0.12cm -1 and uEf=0.185cm -1.
= 0 I--- 1--- ► r(cm)
F ig u re 3.1 One-dim ensional, Bare, Hom ogeneous, Cylindrical Reactor 3.1.1 A n aly tical S olution
One-group diffusion equation can be written as 1 ^ r # ( r ) _ E ^ ( r ) = _ v E f , ( r . r dr dr D k eff D (3.1) Sim plifying (3.1) 1 d d ^ (r) 2 r + B 2^( r ) = 0 r dr dr where B is defined as (3.2) B 2 v Ef k eff ■Ea D Solution of (3.2) is (3.3)
W( r) = AJo(Br) (3.4)
Zero incom ing current boundary condition is
= 0 j-(R ) = * R ) + D d t 4 2 dr r=R Using (3.4) in (3.5) gives J0(BR) = 2D B J1(BR) (3.6) can be written as J0(x) = —z ~ Ji(x) R
where x=BR. Substituting the num erical values of D and R in (3.7) gives 2.884615382J0(x) - x J1(x) = 0
(3.8) can be solved by using N ew ton’s M ethod  f(x (t)) x (t+1) = x (t) -f (x (t)) (3.5) (3.6) (3.7) (3.8) (3.9)
where, t is the iteration number. Here f(x) is the left hand side of (3.8). f ’(x) can be found using recurrence relation for Bessel Function of first kind,
x Jn = n Jn - x Jn+1 = - n J n + x Jn-1 (3. 10) Thus, (3.9) becomes
x(t+1) = x(t) 2.884615382J0(x (t)) - x (t)J t (x (t))
+ 2.884615382J1(x (t)) + x (t)J 0(x (t)) (3.11) Initial estim ate of x (0) can be found by assuming critical reactor, keff=1. In this case, (3.3) gives B=0.316227766 and x (0)=BR=1.185854123. Table 3.1 shows the results o f N ew ton’s method.
Table 3.1 I t e r a t i o n S t e p s o f N e w t o n ’ s M e t h o d f o r t h e S o l u t i o n o f ( 3 . 7 ) t x ® e (t) ( % ) 0 1 . 1 8 5 8 5 4 1 2 3 -1 1 . 7 9 9 7 7 3 4 4 5 9 4 . 9 8 3 7 2 1 . 7 7 1 1 7 3 5 3 5 1 . 6 1 4 5 3 1 . 7 7 1 2 8 5 9 8 9 0 . 0 0 6 3 4 1 . 7 7 1 2 8 5 9 9 1 9 . 1 1 4 x 1 0 -8 T h e s e c a l c u l a t i o n s a r e m a d e u s i n g M A T H E M A T I C A 5 . 2 . F i n a l l y , x = 1 . 7 7 1 2 8 5 9 9 1 a n d B = x / R = 0 . 4 7 2 3 4 2 9 3 c m -1 . E f f e c t i v e m u l t i p l i c a t i o n f a c t o r c a n b e c a l c u l a t e d a s v £ f k eff = --- - f — = 0 . 6 9 8 0 6 0 2 6 4 ( 3 . 1 2 ) eff D B 2 + £ a A v e r a g e f l u x i s d e f i n e d a s R
j^ ( r ) 2 n r d r n R 2 A v e r a g e f l u x c a n b e c a l c u l a t e d f r o m ( 3 . 1 3 ) - 2 R ^
p j2 J A J 0 ( 0 . 4 7 2 3 4 2 9 3 r ) r d r = 0 . 5 8 0 8 6 1 3 7 5 4 6 A0 0 ( 3 . 1 3 ) ( 3 . 1 4 ) I n o r d e r t o f i n d a n e x p r e s s i o n f o r A , i t i s n e c e s s a r y t o m a k e a s e p a r a t e c a l c u l a t i o n o f t h e r e a c t o r p o w e r . I n p a r t i c u l a r , t h e r e a r e £ f ^ ( r ) f i s s i o n s p e r c m 3/ s e c a t t h e p o i n t o f r , a n d i f t h e r e c o v e r a b l e e n e r g y i s w f j o u l e s p e r f i s s i o n ( w f = 3 . 2 x 1 0 -11 j o u l e s ) , t h e n t h e t o t a l p o w e r p e r a x i a l d i s t a n c e , i n w a t t s / c m , i s
RP = w f £ f
j^ ( r ) 2 n r d r ( 3 . 1 5 ) 0 P e r f o r m i n g t h e i n t e g r a t i o n g i v e s P = w f £ f n R 2 0 . 5 8 0 8 6 1 3 7 5 4 6 A ( 3 . 1 6 ) I f t h e r e a c t o r p o w e r i s g i v e n a s P = 2 0 0 0 w a t t / c m , a n d t h e m a c r o s c o p i c f i s s i o n c r o s s s e c t i o n i s £ f = 0 . 0 7 6 4 c m -1 , t h e n A i s c a l c u l a t e d a s A = _________________P _________________ w f £ f n R 2 0 . 5 8 0 8 6 1 3 7 5 4 6 3 . 1 8 9 4 9 5 4 x 1 0 13 ( 3 . 1 6 )
F l u x d i s t r i b u t i o n c a n b e f o u n d b y m u l t i p l y i n g c o n s t a n t A w i t h c e l l a v e r a g e f l u x e s .
F i n a l l y , a v e r a g e f l u x i s c a l c u l a t e d u s i n g ( 3 . 1 6 ) i n ( 3 . 1 4 ) a s
$ = 1 . 8 5 2 6 5 4 6 8 4 x 1 0 13 n 0/ c m 2 s e c . ( 3 . 1 7 )
3.1.2 Solution with Nodal Method
H e r e , f o r m u l a t i o n s d e r i v e d i n p r e v i o u s s e c t i o n a r e t e s t e d w i t h t w o n o d e c a l c u l a t i o n s b e f o r e t h e i r i m p l e m e n t a t i o n i n t o t h e c o m p u t e r p r o g r a m . F i g u r e 3 . 2 s h o w s t h e n o d e s o f t h i s s y s t e m . S i n c e A a n d D a r e t h e s a m e f o r t h e s a m e m a t e r i a l i n o n e - g r o u p c a l c u l a t i o n , t i s t h e s a m e f o r t w o n o d e s . F r o m ( 2 . 7 3 ) , t= 1 2 . 3 1 5 2 A i s c a l c u l a t e d b y d i v i d i n g r a d i u s , R b y t h e n u m b e r o f n o d e s , 2 . T h e r e f o r e A = 1 . 8 7 5 Figure 3.2 C y l i n d r i c a l R e a c t o r w i t h T w o N o d e s C a l c u l a t e d v a l u e s o f t h e v a r i a b l e s d e s c r i b e s i n ( 2 . 7 6 ) - ( 2 . 8 0 ) a r e g i v e n i n t a b l e 3 . 2 . R a d i i a n d c o r r e s p o n d i n g s u r f a c e a r e a s a r e g i v e n i n T a b l e 3 . 3 . T h e y a r e u s e d i n ( 2 . 7 6 ) t h r o u g h ( 2 . 8 0 ) . Table 3.2 C a l c u l a t e d V a l u e s o f t h e V a r i a b l e s i n t h e T w o - N o d e E q u a t i o n s i 1 2 m - 0 . 1 1 2 5 9 8 - 0 . 1 8 7 6 6 3 n - 0 . 4 9 9 8 0 5 - 0 . 4 2 4 7 4 t - 0 . 3 3 7 7 9 4 - 0 . 2 6 2 7 2 9 o - 0 . 2 7 4 6 0 9 - 0 . 3 4 9 6 7 4 p 0 . 4 0 3 1 0 1 0 . 4 0 3 1 0 1
Table 3.3 Radii and Related Surface Areas i r i ( c m ) S i ( c m 2 ) 1 / 2 0 0 1 0 . 9 3 7 5 5 . 8 8 7 5 3 / 2 1 . 8 7 5 0 1 1 . 7 7 5 0 2 2 . 8 1 2 5 1 7 . 6 6 2 5 5 / 2 3 . 7 5 0 0 2 3 . 5 5 0 0 T a b l e 3 . 2 c a n b e u s e d t o t e s t f u r t h e r c o m p u t e r p r o g r a m s w h i c h w i l l b e d e v e l o p e d u s i n g n o d a l m e t h o d s i n c y l i n d r i c a l g e o m e t r y . E q u a t i o n s d e r i v e d i n s e c t i o n ( 2 . 5 . 3 ) a r e f o r m e d . T h e s e e q u a t i o n s c o n s t i t u t e t h e m a t r i x , A , w i t h n u m e r i c a l v a l u e s " 1 . 2 7 4 6 1 - 0 . 4 0 3 1 0 1 0 0 . 3 3 7 7 9 4 0 0 ' 0 1 8 . 8 8 8 5 4 - 8 . 8 8 8 5 4 0 0 0 . 1 1 2 5 9 8 - 0 . 4 0 3 1 0 1 1 0 . 4 9 9 8 0 5 0 0 A = 0 0 0 . 3 4 9 6 7 4 1 - 0 . 4 0 3 1 0 1 0 ( 3 . 1 8 ) 0 0 - 2 . 9 6 2 8 4 2 . 9 6 2 8 4 1 5 . 9 2 5 6 9 0 0 0 . 1 8 6 6 3 0 - 0 . 4 0 3 1 0 1 1 6 x6 T h i s m a t r i x h a s 3 u p p e r d i a g o n a l s a n d 3 l o w e r d i a g o n a l s . I t i s a b a n d m a t r i x . A l l o t h e r A m a t r i c e s w i t h i n c r e a s e d n o d e n u m b e r h a v e t h e s a m e c h a r a c t e r . T h e r e f o r e , i n o r d e r t o u s e c o m p u t e r m e m o r y e c o n o m i c a l l y , o n l y 7 d i a g o n a l s o f t h e s e m a t r i c e s a r e s t o r e d . T h e y a r e t r a n s f e r r e d i n t o a n e w m a t r i x i n t h e N E M R c o m p u t e r p r o g r a m . A l l t h e e l e m e n t s o f t h e m a t r i x F a r e z e r o e x c e p t F ( 2 , 2 ) = F ( 5 , 5 ) = 1 . 5 4 1 6 0 5 . H e n c e , o n l y t h e s e e l e m e n t s a r e u s e d i n c a l c u l a t i o n s o f N E M R . I f a n e w m a t r i x d e f i n e d a s B = A -1F ( 3 . 1 9 ) ( 2 . 1 1 0 ) b e c o m e s A -1
FJ = B J = k e f f J ( 3 . 2 0 ) k e f f i s t h e m a x i m u m e i g e n v a l u e o f t h e m a t r i x , B . M A T H E M A T I C A 5 . 2 f i r s t
keff=0.689617 (3.21) and the error
(Analytic keff - N odal keff Error% =
J---Analytic k eff
This is a reasonable error as will be show n in the follow ing section. 3.1.3 NEMR and QFEMR Results
Linear FEM, quadratic FEM and N EM results of m ultiplication factor are given in Table 3.4
Lx 100% = 1.2095% (3.22)
Table 3.4 keff Results of QFEM R and N EM R Programs
M ethod N um ber o f Elem ents(FEM ) or Nodes(NEM ) keff Error (%) Linear FEM 2 0.746411132 6.92646 10 0.699714822 0.23702 50 0.698125364 0.00933 100 0.698076062 0.00226 Quadratic FEM 1 0.720881392 3.26922 2 0.699073117 0.14509 10 0.698061073 0.0001158 50 0.698059636 0.00009 Low est Order NEM 3 0.694169485 0.55737 5 0.696634106 0.204303 11 0.697764156 0.042418 21 0.697979689 0.011543
Quadratic FEM gives better results than linear FEM. In the FEM term inology, the points in the m esh w here the unknowns are introduced are called nodes. In the nodal method, a node is not a point, it describes a cell. A linear finite elem ent has two nodes; both are located at the endpoints o f the element. Therefore, N elem ents correspond to N+1 node in the linear FEM. In case of quadratic finite elements, the num ber of nodes per elem ent increases to three. N quadratic elem ents contain 2N+1 nodes.
NEM is better than linear FEM . NEM also seems the best m ethod for 3 nodes which corresponds to 1 quadratic elem ent in the quadratic FEM , but then, quadratic FEM gives more accurate results than other methods as shown in figure 3.3.
0,720 n , 0,715 0,710 0,705 0,700 __ = ___ ~ 0,695
-— ySrTTZ----0,690 - / ---0,685 - 4 ---0,680 1 6 11 16 21 26 31 Number of Nodes
Figure 3.3 keff Results of FEM and NEM
N EM R code with two nodes gives keff =0.68963090 with an error 1.20754102%. They are very sim ilar to tw o-node calculations of keff (3.21) and its error (3.22). This shows that nodal equations are im plem ented into the code correctly.
If the num ber of nodal cells is large, the nodal m ethod loses its calculational advantages with respect to finite elem ent methods. But if the num ber o f nodes are small as 1, 2, 3 nodes, nodal m ethod seems more advantageous than finite elem ent m ethods as seen in figure 3.3.
(3.13) can be discretized as --- nodal analytic ---linear FEM --- quadratic FEM ^ = N I f