Development of a nodal method for the solution of the neutron diffusion equation in cylindrical geometry

Development Of A Nodal Method

For The Solution Of The Neutron

Diffusion Equation In Cylindrical

Geometry

Mehmet MERCIMEK

Istanbul Technical University Energy Institute

NODAL FORMALISM

1. The multigroup neutron balance equation

 

   

   

   

r r , g 1,2, ,G k r r r r r J G 1 g g g f g 1 g 1 g g g g g g r g                        





          2. Fick’s Law

 

r D

 

r

 

r J g _g g         

Basis of nodal formalism:

 

r

   

r r Q

 

r J  _r           

Nodal Mesh

 Decompose the reactor core into relatively large

subregions (nodes) in which the material composition and flux are assumed uniform.

Cell and Edge Averaged

Quantities

E

dge-based

unknowns



  











ir 1/2 1/2 ir 2 1/2 -i 2 1/2 i i

r

rdr

)

r

(r

2



  





ir 1/2 1/2 -ir r 2 1/2 -i 2 1/2 i i

Q(r)rdr

)

r

r

(

2

Q

The

cell-based

unknowns

Φ_i+1/2, Φ_i-1/2 and J_i+1/2, J_i-1/2

Relationships between

edge-based quantities

Under P₁ Approximation

 

 

 

2 r J n 4 r r j u u       







   







_{i 1/2}

2

j

_{i 1/2}

j

_{i 1/2}







   







_{i 1/ 2}

2

j

_{i 1/2}

j

_{i 1/2} -1/2 i 1/2 i 1/2 i

j

j

J

_



_



_     1/2



i 1/2



i 1/2 i

j

j

J

-r r r r

n

J.

j

j

J









Discrete Nodal Balance

Equation

i i

2

r

S





2 r r r i 1/2 i 1/2 i     2 S S S i 1/2 i 1/2 i    



i 1/2 i 1/2



i 1/2



i 1/2 i 1/2



r i i i i i i 1/2 i j j S j j S Q S S _ _  _  _ _  _      





_

_



           i 1/2 1/2 -i 1/2 i 1/2 i i r 1/2 i 1/2 i r r Q(r)rdr 2 (r)rdr r r Σ 2 rdr rJ(r) dr d r r r 1 2

Nodal Expansion Method (NEM)



In this lowest order form, NEM considers a

quadratic expansion of the transverse

averaged flux on each cell



The expansion coefficients are determined by

applying Fick’s law in combination with

continuity of normal current.

Polynomial Basis 1

1/2 i 1/2 -i l lP(r), r r r N 0 l a r   _    



i i r -r  

 where ξ=±1/2 when r=r_i+1/2

1/2 ), ( P N 0 l al l         



P₀()=1

 









2 / 1 2 / 1 l

ξ

dξ

0

,

l

0

P

Polynomial Basis 2



        1/2 1/2 -i i i 2 1/2 -i 2 1/2 i i r ) d ) r (r 2 ) 2 / 1 ( P a ) 2 / 1 ( P a ) 2 / 1 ( P a 1/2) ( P 2 0a 0 0 1 1 2 2 1/2 i       _



_{l l} l ) 2 / 1 ( P a ) 2 / 1 ( P a ) 2 / 1 ( P a 1/2) ( P 2 0 a _{0 0 1 1 2 2} 1/2 i           _



_{l l} l

 

   1 P

 

4 1 ξ 3 ξ P₂  2 

Polynomial Basis 3



 _  _



 _{i i} 1 a i i i i i i i 0 r r a                                               _{ i} i i i i i i i i 2 r r 2 r r a                                                        _{ }   i 2 i i i i i 2 i i i i i ξ ) / r ( 2 1 3 4 1 ) / r ( 8 1 ξ ) / r ( 2 1 3 4 1 ) / r ( 8 1 6 2 3 1/2 ), ( P N 0 l al l     

Fick’s Law 1













d

)

(

d

D

-)

J(

i i                                      _{ }  i 1/2 i i 1/2 i i i i i i 1/2 i ) / r ) / r D J                                     _{ }  i 1/2 i i 1/2 i i i i i i 1/2 i ) / r ) / r D J

Fick’s Law 2

i i i i i i 1/2 i i 2 i 2 i i i 1/2 i i i i i 1/2 i D 4 1 D 6 j D 48 r D 2 1 j 4 r 1 D 2 j                                        i i i i i i 1/2 i i 2 i 2 i i i 1/2 i i i i i 1/2 i 4D D 6 j D 48 r D 2 1 j 4 r 1 D 2 j _                                                    _  2 i i i i D 48 D 16 1

Equations per Node

i i 1/2 i i 1/2 i i 1/2 i

m

j

n

j

p

j

_



_



_





i i 1/2 i i 1/2 i i 1/2 i

t

j

o

j

p

j

_



_



_







 







_                         1/2 i 1/2 i 1/2 i 1/2 i 1/2 i 1/2 i 1/2 i 1/2 i i i i r i S j S j S j S j S 1 Q 1

Matrix Form

S

J

F

k

1

J

A

eff





J

unknown vector A 3Nx3N band matrix,

where N is the total number of nodes

F

3Nx3N diagonal matrix (fission source term)_{only (3i-1)th elements contain nonzero}

Matrix Equation 1

N N eff N N _k F J 1 J A  Nx3N 3 N N N a 1/2 N N a 1/2 -N N a 1/2 -N 4)) -(3x(3N N N 1 -N 1 -N a 1/2 -N 1 -N 3)) -N 3 ( 3) -N 3 (( 1 -N )) 3 ) 6 3 (( N 1 p 0 m S S 1 S S S S -0 0 p 1 o 0 0 n 0 0 S S 0 0 t A 0 A                                                    x x N

Matrix Equation 2







       _{ }  _{1/2 3/2 3/2 N-1/2 N 1/2} T N j j j ... j j J Nx3N 3 (3x3) 1 ) 3 3 ( ) 3 3 ( ) 3 3 ( (3x3) 1 ) 3 3 ( ) 3 3 ( ) 3 3 ( (3x3) 1 N F ... 0 0 0 ... F 0 0 ... 0 F F                  x x x x x x

Iteration

J



is known from initial estimates



New

J

vector is found with a linear system solver



New k_{eff estimate is found after fission source iteration}



Iteration continues until the difference between two successive k_eff estimates drops below the convergence criterion. N N eff N N _k F J 1 J A 

One-group, Bare,

Homogeneous Reactor

D=0.65cm, a=0.12cm-1 and f=0.185cm-1.            r D k 1 D r dr r d r dr d r 1 _f eff a        r B dr r d r dr d r 1 ₂

One-Group Reflected Reactor

fuel refl.

R2=5cm

 Fuel is the same material as previous

problem with same radius.

 Reflector is a graphite material with

thickness 1.25cm.

 Absorption cross-section and diffusion

coefficient of the graphite are taken to be 0.00032cm-1 and 0.84cm

Analytical Solution



B.C.



   _{  }  R₁ R₁

|

|

dr R d D dr R d D(1)  1  _ (2)  1         ₂ (2) R₂ dr d 2 D R 4 1 )) R ( K ) R ( (-LI C ) R ( J C R ₁ 0 1 R 0 4 1 c 0 1      )) R ( K ) R ( I (-L D C ) R ( J D C - ₁ cc ₁ c ₁  ₄ R R ₁ R ₁  R ₁ R ₁ )) R ( K ) R ( I (L D ) R ( K ) R ( LI -) R ( J D ) R ( J 1 R 1 R 1 R 1 R R 1 R 0 1 R 0 1 c 1 c c 1 c 0             k =0.768077605

Two-group, Bare, Reactor

First, critical radius is calculated for zero incoming current boundary condition. Then, QFEMR and NEMR results are compared to see how close they can calculate k_eff to critical value.      S





₀ D S dr d r dr d r 1 1 R1 f2 f1               





0 D -S dr d r dr d r 1 2 a2 2 S,1   _{ }      

Analytical Solution

R)

(B

J

D

2B

R)

(B

J

_{0 1}



_{1 1 1 1}

R)

(B

J

D

2B

R)

(B

J

_{0 2}



_{2 2 1 2} R=41.8931096cm _{S=0.116309825}

TRIGA

 Fast and thermal cross-sections are obtained after

the ring homogenizations

 NEMR finds the effective multiplication factor as

k

_eff

 By multiplying the number of annular regions in the

basic mesh by an integer (degree of refinement), finer meshes may be produced.

 QFEMR with 320 quadratic elements or 16 degree of

refinement (finest mesh) gives keff=1.21051196.

% 00255016 . 0 21051196 . 1 % 100 | 21054283 . 1 21051196 . 1 | Error(%)    

CONCLUSION

 NEM and quadratic FEM were shown to be of better accuracy with respect to linear

FEM. It also appears that NEM is a practical method for the problems in which the mesh is very coarse (1, 2, 3 nodes etc.).