Neutron Diffusion Equation Alpha-Eigenvalue Calculation

Description: Solves neutron diffusion equation in slab geometry. Determines alpha-eigenvalue using LANL's first algorithm (Hill).

Alpha-Eigenvalue Neutron Diffusion Equation in Slab with Fission Source

The NDE in a slab is given by

$$ -\frac{d}{dx}D(x)\frac{d\phi(x)}{dx} + \bigg ( \Sigma_a + \frac{\alpha}{v} \bigg ) \phi(x) = \frac{1}{k}\nu \Sigma_f \phi(x) $$

where $D(x)$ is the diffusion coefficient, $\Sigma_a$ and $\Sigma_f$ are the absorption and fission macroscopic cross sections, $\nu$ is the average number of neutrons emitted in fission, $k$ is k-effective, and $\alpha$ is the alpha-eigenvalue.

Import Python Libraries



In [16]:

    
%matplotlib inline

import numpy as np
import matplotlib.pyplot as plt

Material Properties



In [17]:

    
D = 0.9
nusigf = 0.70
siga = 0.066

Slab Geometry Width and Discretization



In [18]:

    
#Lx = np.pi*((nusigf-siga)/D)**(-0.5)
Lx = 15.0

N = 55;
h = Lx/(N-1)

x = np.zeros(N)

for i in range(N-1):
    x[i+1] = x[i] + h

Generation of Leakage Matrix



In [19]:

    
L = np.zeros((N,N))
A = np.zeros((N,N))
M = np.zeros((N,N))

for i in range(N):
    L[i][i] = L[i][i] + (-2*(-D/(h**2)))
    
for i in range(1,N):
    L[i][i-1] = L[i][i-1] + (1*(-D/h**2))
    
for i in range(N-1):
    L[i][i+1] = L[i][i+1] + (1*(-D/h**2))

Algorithm

Using Hill's (LANL) algorithm, we calculate the $\alpha$-eigenvalue as follows:

Guess $\alpha_0$ (usually $\alpha_0$ = 0)
Solve for $k_0$
Modify the initial $\alpha_0$ using some modifier ($\alpha_1 = \alpha_0 + \xi$)
Solve for $k_1$
Using the points $(k_i,\alpha_i)$ and $(k_{i-1},\alpha_{i-1})$, linearly extrapolate $\alpha_{i+1}$ such that $k_{i+1} = 0$.
Repeat until $k_i = 1$.

Problems:

1) If the system is subcritical enough

\begin{equation} \frac{\alpha}{v} - \Sigma_{a} < 0 \end{equation}

which causes instabilities (confirm by looking at the condition number of matrix iterations). Need to add parameter to prevent this according to Warsa and Fichtl.

2) Bad guesses of $\alpha_0$ causes algorithm to fail and does not converge.



In [20]:

    
#Generate flux vector

#Tolerance, k-effective and alpha initial guesses, alpha-eigenvalue modifier
tol = 1e-10
k = 1.00
alpha = 0.0
evm = 0.01

#Alpha-eigenvalue outer iteration
for j in range(1,100):
    
    kprev = k
    
    A = np.zeros((N,N))
    phi0 = np.ones((N,1))
    phi0[0] = 0
    phi0[N-1] = 0
    k = 1.0
    
    for i in range(N):
        A[i][i] = A[i][i] + siga + alpha
    
    M = L + A
    M[0][0] = 1
    M[0][1] = 0
    M[N-1][N-1] = 1
    M[N-1][N-2] = 0

    #k-effective inner iteration
    for i in range(100):
    
        kold = k
        psi = np.linalg.solve(M,nusigf*phi0)
    
        k = sum(nusigf*psi)/sum(nusigf*phi0)
        phi0 = (1/k)*psi
        phi0[0] = 0
        phi0[N-1] = 0
    
        residual = np.abs(k-kold)
    
        if residual <= tol:
            break
    
    #Modify alpha-eigenvalue after first iteration, or linearly extrapolate new guess such that k equals 1
    #Reminder: does not calculate both k and alpha
    if j == 2:        
            
        alpha_prev = alpha
        alpha = alpha + evm
        
    elif j > 2:
        
        #print "alpha-alpha_prev = ", alpha-alpha_prev
        #print "alpha = ", alpha
        #print "alpha_prev = ", alpha_prev
        #print "j = ", j
        
        if abs(alpha - alpha_prev) < tol:
            
            break
        
        alpha_new = alpha_prev + (1-kprev)/(k-kprev)*(alpha-alpha_prev)
        alpha_prev = alpha
        alpha = alpha_new
        
    else:
        
        continue
        
    #if abs(k-1) < tol:
    #    print "alpha = ", alpha
    #    print "k-effective = ", k
    #    break
        
#print("abs(k-1) = %.20f" % abs(k-1))        
        
print("alpha = %.15f" % alpha)
print('k-effective = %.15f' % k)

#plt.plot(x,phi0)
#plt.xlabel('Slab (cm)')
#plt.ylabel('Neutron Flux')
#plt.grid()

#print " alpha = ", (k-1)/k * sum(nusigf*phi0)/sum(phi0)









    



alpha = 0.594532716010151
k-effective = 0.999999999999995



In [ ]: