ChEn-3170: Computational Methods in Chemical Engineering Spring 2020 UMass Lowell; Prof. V. F. de Almeida 29Apr20
$ \newcommand{\Amtrx}{\boldsymbol{\mathsf{A}}} \newcommand{\Bmtrx}{\boldsymbol{\mathsf{B}}} \newcommand{\Mmtrx}{\boldsymbol{\mathsf{M}}} \newcommand{\Imtrx}{\boldsymbol{\mathsf{I}}} \newcommand{\Pmtrx}{\boldsymbol{\mathsf{P}}} \newcommand{\Lmtrx}{\boldsymbol{\mathsf{L}}} \newcommand{\Umtrx}{\boldsymbol{\mathsf{U}}} \newcommand{\Smtrx}{\boldsymbol{\mathsf{S}}} \newcommand{\xvec}{\boldsymbol{\mathsf{x}}} \newcommand{\uvar}{\boldsymbol{u}} \newcommand{\fvar}{\boldsymbol{f}} \newcommand{\avec}{\boldsymbol{\mathsf{a}}} \newcommand{\bvec}{\boldsymbol{\mathsf{b}}} \newcommand{\cvec}{\boldsymbol{\mathsf{c}}} \newcommand{\rvec}{\boldsymbol{\mathsf{r}}} \newcommand{\mvec}{\boldsymbol{\mathsf{m}}} \newcommand{\gvec}{\boldsymbol{\mathsf{g}}} \newcommand{\zerovec}{\boldsymbol{\mathsf{0}}} \newcommand{\norm}[1]{\bigl\lVert{#1}\bigr\rVert} \newcommand{\transpose}[1]{{#1}^\top} \DeclareMathOperator{\rank}{rank} \newcommand{\Power}{\mathcal{P}} $
Background for this topic can be found in the textbook: Dynamics of Nuclear Reactors, David L. Hetrick, 1993, ANS.
A power nuclear reactor is a heat generation device. In most cases, nuclear heat (heat generated by nuclear fission reactions) is then converted in work through a heat-work cycle on the heated fluid passed through the reactor (figure below). A very simplified modeling of all types of nuclear reactors has been developed in the early days of nuclear engineering. The point-reactor model is an introduction to the subject of nuclear reactor dynamics. It disregards space variations of the neutron density (i.e. it is a pointwise approach), it considers dynamics for one-energy group of neutrons, it only applies to conditions when the reactor is near criticality and the nuclear fuel is nearly static. In this model, a pointwise neutron balance is performed using the main components of neutron transport, namely prompt fission, delayed fission, scattering, and absorption (figure below); scattering is taking into account in the diffusion limit.
Given the neutron generation time $\ell$, delayed neutron fraction, $\beta$, decay constants of a six-group delayed neutron emitters, $\lambda_i$, and corresponding yield of delayed neutron fractions for each emitter, $\beta_i$, calculate the pointwise neutron density variation with time for a step change in neutron reactivity, $\rho$.
Name | Parameter | Value | Unit | |
---|---|---|---|---|
neutron generation time | $\ell$ | $1\times 10^{-4}$ | s | |
delayed neutron fraction | $\beta$ | $6.5\times 10^{-3}$ | - | |
reactivity feedback coef. | $\alpha_n$ | $-3\times 10^{-4}$ | cm$^2$ |
For thermal fission of $^{235}$U the following delay neutron data is typically used (Dynamics of Nuclear Reactors, David Hetrick, 1993, ANS textbook).
Delayed neutron emitter group No. | Decay cte ($\lambda_i$,1/sec) | Relative yield ($\beta_i/\beta$) | |
---|---|---|---|
1 | 0.0124 | 0.033 | |
2 | 0.0305 | 0.219 | |
3 | 0.111 | 0.196 | |
4 | 0.301 | 0.395 | |
5 | 1.14 | 0.115 | |
6 | 3.01 | 0.042 |
'''Parameters'''
params = dict()
params['gen_time'] = 1.0e-4 # s
params['beta'] = 6.5e-3 #
params['species_decay'] = [0.0124, 0.0305, 0.111, 0.301, 1.14, 3.01] # 1/sec
params['species_rel_yield'] = [0.033, 0.219, 0.196, 0.395, 0.115, 0.042]
params['alpha_n'] = -3e-4 # cm^2
The space-invariant neutron balance for the point-reactor model without a startup source (or extraneous source) is
\begin{equation*} \frac{\text{d}n}{\text{d}t} = \frac{\rho(t)-\beta}{\ell}\, n + \sum\limits_{i=1}^{6} \lambda_i\,c_i , \end{equation*}where the first term on the right side of the equation represents the net production of neutrons not account for delayed neutrons, and the second term accounts for the source of delayed neutrons considering 6 groups of delayed neutron emitters resulting from the fission of $^{235}$U nuclei. Therefore a balance of neutron emitter species is also necessary
\begin{equation*} \frac{\text{d}c_i}{\text{d}t} = \frac{\beta_i}{\ell}\, n - \lambda_i\,c_i , \ \ \ \ \ \forall \ \ \ \ \ i=1,\ldots,6. \end{equation*}where the first term on the right side of the equation is the source of emitters as a function of the neutron density $n(t)$, and the second term is the consumption rate of the emitter by radioactive decay obtained as a function of the product of the concentration of the emmiter, $c_i(t)$, multiplied by its decay constant $\lambda_i$. Here the concentration of of the $i$th emitter, $c_i$ is considered in terms of delayed neutron density, hence the units are the same as $n$.
The reactivity function with feedback is postulated as
\begin{equation*} \rho(t) = \rho_0 + \alpha_n\,\bigl(n(t) - n(0)\bigr), \end{equation*}where $\alpha_n \le 0$ is the feedback coefficient and the externally applied reactivity is $\rho_0$.
There exists seven equations, the neutron density balance, and six neutron emitting species balance. There exists seven unknown variables, $n(t)$ and $c_i(t)$. Therefore there are as many unknowns as there are equations and the problem is potentially solvable given the initial conditions for the unknowns, and all the values of the parameters in the foregoing table.
A vector notation for the foregoing system of equations greatly improves the generality of the derived computer code. Towards this goal let us define
\begin{equation*} \frac{d\uvar}{dt} = \fvar( \uvar, t ) \end{equation*}where $\uvar(t) = (u_1,u_2,u_3,u_4,u_5,u_6,u_7)$ and we assign
\begin{align*} u_1(t)&=n(t),\\ u_2(t)&=c_1(t),\\ u_3(t)&=c_2(t),\\ u_4(t)&=c_3(t),\\ u_5(t)&=c_4(t),\\ u_6(t)&=c_5(t),\\ u_7(t)&=c_6(t). \end{align*}Also for $\fvar(\uvar,t) = \bigl(f_1(\uvar,t), f_2(\uvar,t), f_3(\uvar,t), f_4(\uvar,t), f_5(\uvar,t), f_6(\uvar,t), f_7(\uvar,t)\bigr)$ we assign
\begin{align*} f_1 & = \frac{\rho(t)-\beta}{\ell}\, u_1 + \sum\limits_{i=2}^{7} \lambda_i\,u_i, \\ f_2 & = \frac{\beta_1}{\ell}\, u_1 - \lambda_1\,u_2, \\ f_3 & = \frac{\beta_2}{\ell}\, u_1 - \lambda_2\,u_3, \\ f_4 & = \frac{\beta_3}{\ell}\, u_1 - \lambda_3\,u_4, \\ f_5 & = \frac{\beta_4}{\ell}\, u_1 - \lambda_4\,u_5, \\ f_6 & = \frac{\beta_5}{\ell}\, u_1 - \lambda_5\,u_6, \\ f_7 & = \frac{\beta_6}{\ell}\, u_1 - \lambda_6\,u_7 . \end{align*}Finally, the initial conditions given are as follows:
\begin{align*} u_1(0)&=n_0,\\ u_2(0)&=c_{1_0},\\ u_3(0)&=c_{2_0},\\ u_4(0)&=c_{3_0},\\ u_5(0)&=c_{4_0},\\ u_6(0)&=c_{5_0},\\ u_7(0)&=c_{6_0} \end{align*}where $n_0$ is a steady-state value of the neutron density, and the corresponding steady-state values for the concentration of delayed neutron emitters is
\begin{equation*} c_{i_0} = \frac{\beta_i\,n_0}{\lambda_i\,\ell} \ \ \ \ \forall \ \ \ \ i=1,\ldots,6 \end{equation*}'''ODE function'''
def f_vec( time, u_vec, params ):
import numpy as np
assert np.all(u_vec >= 0.0)
n_dens = u_vec[0]
gen_time = params['gen_time']
beta = params['beta']
species_decay = params['species_decay']
lambda_vec = np.array(species_decay)
species_rel_yield = params['species_rel_yield']
beta_vec = np.array(species_rel_yield) * beta
rho_0 = params['reactivity_0']
alpha = params['alpha_n']
n_dens_0 = params['n_0']
reactivity = rho_0 + alpha * (n_dens - n_dens_0)
c_vec = u_vec[1:]
#print('c_vec = ',c_vec)
f_tmp = np.zeros(7,dtype=np.float64) # vector for f_vec return
# neutron balance
f_tmp[0] = (reactivity - beta)/gen_time * n_dens + lambda_vec @ c_vec
#print('a=',(reactivity - beta)/gen_time * n_dens)
#print('b=',lambda_vec @ c_vec)
#print('f_tmp[0]',f_tmp[0])
# loop over 6 species balance
#for i in range(6):
#f_tmp[i+1] = species_rel_yield[i]*beta/gen_time * n_dens - lambda_vec[i] * u_vec[i+1]
f_tmp[1:] = beta_vec/gen_time * n_dens - lambda_vec * c_vec
#print('f_tmp=',f_tmp)
#print('time=',time)
#print('')
return f_tmp
'''Create the point-reactor run function'''
def run_point_reactor(f_vec, time_stamps, params):
from scipy.integrate import odeint
max_n_steps_per_time_step = 100 # max number of nonlinear algebraic solver iterations per time step
n_0 = params['n_0']
c_vec_0 = params['c_vec_0']
u_vec_0 = np.zeros(7,dtype=np.float64)
u_vec_0[0] = n_0
u_vec_0[1:] = c_vec_0
(u_vec_history, info_dict) = odeint( f_vec, u_vec_0, time_stamps,
args=( params, ),
rtol=1e-4, atol=1e-4, mxstep=max_n_steps_per_time_step,
full_output=True, tfirst=True )
assert info_dict['message']=='Integration successful.',\
'Fatal: scipy.integrate.odeint failed %r'%info_dict['message']
return u_vec_history
'''Evolve the point-reactor from a steady state; no feedback'''
beta = params['beta']
# Manipulated input variable
params['reactivity_0'] = 10/100 * beta # "10 cents"
params['alpha_n'] = 0.0 # set no feedback
n_0 = 5.0 # arbitrary value at steady state
params['n_0'] = n_0
import numpy as np
c_vec_0 = np.zeros(6,dtype=np.float64)
species_decay = params['species_decay']
lambda_vec = np.array(species_decay)
species_rel_yield = params['species_rel_yield']
beta_vec = np.array(species_rel_yield) * beta
gen_time = params['gen_time']
c_vec_0 = beta_vec/lambda_vec/gen_time * n_0
params['c_vec_0'] = c_vec_0
time_final = 100.0 # s
n_time_steps = 100 # number of solution values in time
import numpy as np
time_stamps = np.linspace(0.0, time_final, num=n_time_steps)
u_vec_history = run_point_reactor( f_vec, time_stamps, params )
'''Plot neutron density and delayed neutron emitter concentration in the reactor'''
def plot_results(u_vec_history):
import matplotlib.pyplot as plt
fig, ax1 = plt.subplots(1, figsize=(12, 6))
ax2 = ax1.twinx()
color_ids = np.linspace(0,1,u_vec_history[:,1:].shape[1])
for (j,color_id) in zip(range(u_vec_history[:,1:].shape[1]),color_ids):
color=plt.cm.nipy_spectral(color_id)
ax2.plot(time_stamps,u_vec_history[:,j+1]/params['c_vec_0'][j],'-.',color=color,label=r'$c_%i$'%(j+1) )
ax2.set_ylabel(r'$c_i/c_{i_0}$',fontsize=16,color='black')
ax2.tick_params(axis='y', labelcolor='black', labelsize=14)
ax2.legend(loc='lower right',fontsize=12)
#ax2.grid(True)
ax1.plot(time_stamps,u_vec_history[:,0]/params['n_0'],'r-',label=r'$n/n_0$' )
ax1.set_xlabel(r'Time [s]',fontsize=16)
ax1.set_ylabel(r'$n/n_0$',fontsize=16,color='black')
ax1.tick_params(axis='y', labelcolor='black', labelsize=14)
ax1.tick_params(axis='x', labelsize=14)
ax1.legend(loc='best',fontsize=12)
ax1.grid(True)
plt.title(r'Point-Reactor Model: $\rho/\beta=$'
+str(params['reactivity_0']/params['beta'])
+r'; $\alpha_n=$'+str(params['alpha_n']),
fontsize=18)
fig.tight_layout() # otherwise the right y-label is slightly clipped
plt.show()
print('')
plot_results(u_vec_history)
'''Evolve the point-reactor from a steady state with negative feedback'''
beta = params['beta']
params['alpha_n'] = -3e-4 # cm^2
# Manipulated input variable
params['reactivity_0'] = 10/100 * beta # "10 cents"
n_0 = 5.0 # arbitrary value at steady state
params['n_0'] = n_0
import numpy as np
c_vec_0 = np.zeros(6,dtype=np.float64)
species_decay = params['species_decay']
lambda_vec = np.array(species_decay)
species_rel_yield = params['species_rel_yield']
beta_vec = np.array(species_rel_yield) * beta
gen_time = params['gen_time']
c_vec_0 = beta_vec/lambda_vec/gen_time * n_0
params['c_vec_0'] = c_vec_0
time_final = 300.0 # s
n_time_steps = 100 # number of solution values in time
import numpy as np
time_stamps = np.linspace(0.0, time_final, num=n_time_steps)
u_vec_history = run_point_reactor( f_vec, time_stamps, params )
plot_results(u_vec_history)
'''Create interactive plots'''
import plotly.offline as py
py.init_notebook_mode()
import plotly.graph_objs as go
layout=go.Layout(title="Point-Reactor Model",
xaxis={'title':'Time [s]'},
yaxis=dict(side='left',title='HELLO'),
yaxis2=dict(overlaying='y',anchor='x',side='right',title='HELLO')
)
fig = go.FigureWidget(layout=layout)
scatt = fig.add_scatter()
fig.layout.titlefont.size = 22
fig.layout.titlefont.family = 'Rockwell'
fig.layout.xaxis.title = 'Time [s]'
fig.layout.yaxis.title = 'n/n_0'
#fig.layout.yaxis2.title = 'T [K]'
'''Vary parameters and display results interactively'''
from ipywidgets import interact
import ipywidgets as widgets
#display(widgets.FloatSlider(readout_format='8.4e'))
@interact(rho_over_beta=widgets.FloatSlider(description=r'$\rho/\beta=$',value=0.01,min=-0.1,max=0.1,step=0.01,readout_format='5.3f'),
gen_time=widgets.FloatSlider(description=r'$\ell=$',value=1e-4,min=1e-4,max=1e-3,step=1e-4,readout_format='5.3e'),
alpha_n=widgets.FloatSlider(description=r'$\alpha_n=$',value=0.0,min=-3e-4,max=1.0e-5,step=1e-5,readout_format='5.3e'),
select=widgets.Dropdown(description='select:',value='neutron_density',options=['neutron_density', 'c_1'])
)
def update(rho_over_beta, gen_time, alpha_n, select):
with fig.batch_update():
scatt.x=time_stamps
beta = params['beta']
params['reactivity_0'] = rho_over_beta * beta
params['alpha_n'] = alpha_n
params['gen_time'] = gen_time
history = run_point_reactor( f_vec, time_stamps, params )
if select == 'neutron_density':
scatt.y=history[:,0]/params['n_0']
scatt.line.color='red'
fig.layout.yaxis.title = 'n/n_0'
elif select == 'c_1':
scatt.y=history[:,1]/params['c_vec_0'][0]
scatt.line.color='blue'
fig.layout.yaxis.title = 'c_1/c_1_0'
else:
pass
fig
interactive(children=(FloatSlider(value=0.01, description='$\\rho/\\beta=$', max=0.1, min=-0.1, readout_format…
FigureWidget({ 'data': [{'line': {'color': 'red'}, 'type': 'scatter', 'uid': '…