`Cortix`

2019 **03Aug2019**

- This is part of the Cortix-on-Jupyter-Notebook examples.
- You must be in a Jupyter Notebook server to run this notebook.
- Select each of the cells below and run them sequentially (use the run button,
`>|`

on the tool bar or use the`Cell`

option on the menu bar). - Alternatively, on the menu bar run all cells:
`Cell -> Run All`

.

$ \newcommand{\Amtrx}{\boldsymbol{\mathsf{A}}} \newcommand{\Bmtrx}{\boldsymbol{\mathsf{B}}} \newcommand{\Smtrx}{\boldsymbol{\mathsf{S}}} \newcommand{\xvec}{\boldsymbol{\mathsf{x}}} \newcommand{\vvar}{\boldsymbol{v}} \newcommand{\fvar}{\boldsymbol{f}} \newcommand{\Power}{\mathcal{P}} \newcommand{\bm}[1]{{\boldsymbol{#1}}} $

This notebook uses Cortix to simulate the startup and shutdown of a boiling water reactor using a point-reactor model, This work is in progress, authors: Prof. Valmor F. de Almeida, UMass Lowell, Dept. of Chemical Engineering; And Austin Rotker. Cortix Group.

The following equations are developed for use in the point-reactor BWR model.

The space-invariant neutron balance for the point-reactor model 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 + q(t), \end{equation*}where the first term on the right side of the equation represents the net production of neutrons not accounting 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 number 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 number concentration of the emmiter, $c_i(t)$, multiplied by its decay constant $\lambda_i$. Here the number concentration of of the $i$th emitter, $c_i$ is considered in terms of delayed neutron number density, hence the units are the same as $n$.

The current model considers the temperature reactivity coefficient, $\alpha_{T_n}$, that is,

\begin{equation*} \rho(t) = \rho_0 + \alpha_{T_n}(T_f)\,\bigl(T_f(t) - T^{(0)}\bigr), \end{equation*}A heat balance in the static fuel dispersed into the moderator gives

\begin{equation*} \frac{\text{d}T_f}{\text{d}t} = - \frac{1}{\rho_f\,c_{pf}}\biggl(q^{'''}(t) - \frac{\dot{Q}_\text{f}}{V_f} \biggr) \end{equation*}where the heat rate sink in the fuel/moderator is

\begin{equation*} {\dot{Q}_\text{f}} = -UA\bigl(T_f-T_c\bigr), \end{equation*}and the nuclear heating power is given by

\begin{equation*} q^{'''}(t) = G_f\frac{\sqrt{\pi}}{2}\,\sigma_\text{fo}\,\sqrt{\frac{T^\text{(o)}}{T}}\,\biggl(\frac{\epsilon\,w_\text{fs}\,\rho_\text{fm}}{M_\text{fn}}\,i\,N_A\biggr)\,n(t)\,v_\text{(o)} . \end{equation*}A heat balance in the flowing coolant fluid produces

\begin{equation*} \frac{\text{d}T_c}{\text{d}t} = -\frac{1}{\tau_c}\biggl(T_c-T_{in}\biggr) - \frac{1}{\rho_c\,c_{pc}}\biggl( \frac{\dot{Q}_\text{f}}{V_c} \biggr) \end{equation*}where the heat rate source is the negative of the heat rate sink in the fuel/moderator.

Neutron density and delayed neutron precursor concentrations at steady-state are related as follows:

\begin{align*} n_\text{ss} &= -\frac{q_\text{ss}\,\ell }{\rho_\text{ss} } \\ c_{i_\text{ss}} &= \frac{\beta_i}{\lambda_i\,\ell}\, n_\text{ss} \ \ \ \ \forall \ \ \ \ i=1,\ldots,6 . \end{align*}The moderator temperature reactivity coefficient is the major contributor to the power coefficient in a boiling water reactor. It can be calculated by the following equation:

\begin{equation*} \alpha_\text{m} = \frac{dK}{dT_\text{m}}\frac{1}{K}\ \end{equation*}Where $\alpha_\text{m}$ is the moderator temperature feedback coefficient, K is the neutron multiplecation factor for a finite medium, and $T_\text{m}$ is the temperature of the moderator. K is calculated as follows:

\begin{equation*} K = \frac{\eta \epsilon P F}{1 + (\tau + L^2)B^2 + \tau L^2 B^4} \end{equation*}Plugging this in to the original moderator temperature reactivity coefficient equation leaves the following differential equation:

\begin{equation*} \alpha_\text{m} = \frac{d}{dT_\text{m}}\frac{1 + (\tau + L^2)B^2 + \tau L^2 B^4}{\eta \epsilon P F}\ \end{equation*}With $\tau$ being the fermi age or slowing down area, $L$ being the neutron diffusion length, $B$ being buckling, $\eta$ being neutron yield, $\epsilon$ being fast fission factor, $P$ being resonance escape probability and $F$ being thermal utilization. Only $\tau$, $P$ and $F$ vary significantly with moderator temperature, with $L$ and $B$ being assumed constants over the range of temperatures.

$\tau$ is found as follows:

\begin{equation*} \tau = \frac{D_\text{f}}{\Sigma_\text{1}} \end{equation*}with $D_\text{f}$ being the diffusion coefficient and equal to:

\begin{equation*} D_\text{f} = \frac{1}{3 \Sigma_\text{s}(1 - \mu_\text{0})} \end{equation*}With $\Sigma_\text{s}$ being the macroscopic scattering cross section, $\sigma_\text{s} N_\text{m}$, and $\mu_\text{0}$ being a constant dependent on the material used as the moderator.

$\Sigma_\text{1}$ is the neutron removal cross section, or the probability per unit length that a neutron becomes thermalized. It is calculated as:

\begin{equation*} \Sigma_\text{1} = \frac{\xi \Sigma_\text{s}}{ln \frac{E_\text{0}}{E_\text{th}}} \end{equation*}With $\xi$ being the average logarithmic energy decrement for a collision involving a neutron and the moderator (in the case of a light water reactor, $\xi$ = 1), $E_\text{0}$ being the energy of a neutron produced during fission (assumed to be constant at around 2 MeV), and $E_\text{th}$ being the energy of the medium that the neutron is diffusing through, which is dependent on the temperature of the moderator by the following equation:

\begin{equation*} E_\text{th} = (T_\text{m})\frac{0.0862 \ MeV}{K} \end{equation*}With $T_\text{m}$ being in kelvin.

Assuming that only $N_text{m}$ and $E_\text{th}$ vary with the temperature of the moderator, then the differential of $\tau$ with respect to $T_\text{m}$ may be calculated as follows:

\begin{equation*} \frac{d \tau}{d T_\text{m}} = \frac{d}{d T_\text{m}} \frac{\frac{1}{3 \sigma_\text{s} N_\text{m}(1 - \mu_\text{0})}}{\frac{\xi \sigma_\text{s} N_\text{m}}{ln \frac{E_\text{0}}{E_\text{th}}}} \end{equation*}\begin{equation*} \frac{d \tau}{d T_\text{m}} = \frac{((0.0862\frac{E_\text{th}}{E_\text{0}})(3 N_\text{m} )-(ln \frac{E_\text{0}}{E_\text{th}})(6 \frac{d N_\text{m}}{d T_\text{m}})}{(3 N_\text{m})^2 (1 - \mu_\text{0})} \end{equation*}Next, the diffusion area in the moderator before absorption, $L^2$, may be evaluated. It can be calculated through the following relation:

\begin{equation*} L^2 = \frac{1}{3 \Sigma_\text{s} \Sigma_\text{a} (1 - \mu_\text{0})} \end{equation*}And therefore:

\begin{equation*} L = \sqrt{\frac{1}{3 \Sigma_\text{s} \Sigma_\text{a} (1 - \mu_\text{0})}} \end{equation*}Note that L is referred to as the diffsion length whereas L^2 is the diffusion area.

Taking the derivative of the diffusion length, L, with respect to time, and noting that only the number density of the moderator varies with time, and that the macroscopic cross sections $\Sigma$ can be rewritten as $\sigma N_\text{m}$, we obtain the differential of diffusion length with respect to time:

\begin{equation*} \frac{dL}{dT_\text{m}} = \frac{1}{2 \sqrt{\frac{-2 \frac{dN_\text{m}}{dT_\text{m}}}{3 \sigma_\text{s} \sigma_\text{a} N_\text{m}^3 (1 - \mu_\text{0})}}} \end{equation*}Next, the terms in the four factor formula, $\eta \ , \ \epsilon \ , \ P, \ F$ must be developed and the derivatives of the ones which vary with moderator temperature must be calculated.

$\eta$ and $\epsilon$ do not vary with moderator temperature, and are assumed to have constant values of 2.02 and 1.03.

$P$, the resonance escape, is calculated via the following equation:

\begin{equation*} P = e^{- \frac{N_\text{f} V_\text{F} I}{\xi \Sigma_\text{s M} V_\text{M}}} \end{equation*}In the above equation, $N_\text{f}$ is the number density of the fuel, $V_\text{F}$ is the volume of the fuel, $I$ is the resonance escape integral, and $V_\text{M}$ is the volume occupied by the moderator.

With only the macroscopic scattering cross section of the moderator varrying with temperature of the moderator (volume of coolant in the core is assumed constant, only the density changes), the derivative of the resonance escape probability with respect to the temperature of the moderator can be calculated as:

\begin{equation*} \frac{dP}{dT_\text{m}} = \frac{-N_\text{f} V_\text{F} I \sigma_\text{s} \frac{dN_\text{m}}{dT_\text{m}}}{(\xi \Sigma_\text{s M} V_\text{M})^2}e^{- \frac{N_\text{f} V_\text{F} I}{\xi \Sigma_\text{s M} V_\text{M}}} \end{equation*}The thermal utilization, $F$, is calculated as follows:

\begin{equation*} F = \frac{\Sigma_\text{aF}}{\Sigma_\text{aF} + \Sigma_\text{aM}} \end{equation*}With $\Sigma_\text{aF}$ being the macroscopic absorption cross section of the fuel and $\Sigma_\text{aM}$ being the macroscopic absorption cross section of the moderator. As only the macroscopic absorption cross section of the moderator will vary with the temperature of the moderator, the derivative of thermal utilization with respect to the temperature of the moderator can be written as follows:

\begin{equation*} \frac{dF}{dT_\text{m}} = \frac{-(\Sigma_\text{aF})(\sigma_\text{aM} \frac{dN_\text{m}}{dT_\text{m}})}{(\Sigma_\text{aF} + \Sigma_\text{aM})^2} \end{equation*}The final preliminary differential equation is the one that must be developed for the differential of the moderator number density with respect to the temperature of the moderator. The change in number density of the moderator with respect to temperature will be proportional to the change in the mass density of the moderator with respect to temperature, by the molar mass of the moderator (in this case, 18 g/mol) and avogadro's number. The derivative of the mass density with respect to the temperature of the moderator of water is a function of a forty term polynomial used by the IAPWS 97 standards to approximate the behavior of water at different temperatures and pressures, that is:

\begin{equation*} \frac{d \rho}{dT_\text{m}} = \frac{d}{dT_\text{m}}(IAPWS(T)) \end{equation*}The derivative can be approximated by using scipy's inbuilt differentiation function. From there, the differential of mass density with respect to moderator temperature may be related to the differential of number density with respect to moderator temperature by the molar mass and avogadro's number:

\begin{equation*} \frac{dN_\text{m}}{dT_\text{m}} = (18 g/mol)(N_\text{A})(\frac{d \rho}{dT_\text{m}}) \end{equation*}Now that all previous relationships have been developed, the moderator temperature feedback coefficient can be calculated. Recall the original equation for the moderator temperature feedback coefficient:

\begin{equation*} \alpha_\text{m} = \frac{d}{dT_\text{m}}\frac{1 + (\tau + L^2)B^2 + \tau L^2 B^4}{\eta \epsilon P F}\ \end{equation*}Taking the derivatives with respect to moderator temperature yields the following differential:

\begin{equation*} \alpha_\text{m} = \frac{(\eta \epsilon P F)((\frac{d \tau}{dT_\text{m}})(B^2 + L^2 B^4)+(\frac{dL}{dT_\text{m}})(2L B^2 + 2L \tau B^4))-(1 + (\tau + L^2)B^2 + \tau L^2 B^4)(\eta \epsilon)(F \frac{dP}{dT_\text{m}} + P \frac{dF}{dT_\text{m}})} {(\eta \epsilon P F)^2} \end{equation*}Temperature and pressure of the coolant leaving the reactor is given by temp_c and pressure; these two variables are used with the IAPWS package (under the IAPWS-97 standards) to determine the specific enthalpy of the steam leaving the reactor. The turbine expands any steam passing through it to 0.75 kPa. This process is NOT assumed to be isentropic, and instead operates at 80% of maximum isentropic efficiency.

We start with a basic heat balance on the turbine:

\begin{equation*} Q + W_\text{s} = \Delta H + \Delta K + \Delta Z \end{equation*}With Q being the heat extracted from the working fluid, W being the work done by the working fluid (and hence being positive if the turbine is producing useful work), and with $\Delta$H, $\Delta$K and $\Delta$Z representing the changes in enthalpy, kinetic energy and gravitational potential energy of the fluid, respectively. It is assumed that the turbine operates adiabatically and with no change in the kinetic or gravitational potential energy of the fluid as it moves through the turbine. Hence, the equation reduces to:

\begin{equation*} W_\text{s} = \Delta H \end{equation*}If the turbine is producing useful work, then W is positive, and hence $\Delta$H must be positive as well. The enthalpy of the fluid will decrease as it passes through the turbine and is expanded to a lower temperature and pressure, hence the equation can be rewritten as:

\begin{equation*} W_\text{s} = \Delta H = H_\text{inlet} - H_\text{outlet} \end{equation*}The inlet steam is at the same state of the steam leaving the top of the reactor, I.E. it has a pressure $P_\text{reactor}$ and the same temperature of the coolant, $T_\text{coolant}$. The turbine expands whatever is taken into the inlet to a final pressure of 0.75 kPa and a temperature $T_\text{runoff}$ which varies based on the properties of the steam that was fed to the turbine. Hence, the work balance equation can be re-written as:

\begin{equation*} W_\text{s} = H_\text{inlet}(T_\text{coolant}, P_\text{reactor}) - H_\text{outlet}(T_\text{runoff}, 0.75 kPa) \end{equation*}Ideally, final temperature, $T_\text{runoff}$, would be the bubble point temperature of the runoff at 0.75 kPa (~ 3 C$^{\circ}$). However, the turbine does not operate isentropically, and instead is assumed to operate at 80% of its maximum isentropic efficiency.

We denote the maximum possible isentropic work of the turbine, $W_\text{isentropic}$, as the amount of work that would be done if the steam entering the turbine was expanded to a pressure of 0.75 kPa and a temperature of 3 C$^{\circ}$.

\begin{equation*} W_\text{isentropic} = H_\text{steam}(T_\text{coolant}, P_\text{reactor}) - H_\text{water}(0.75 kPa, 3 ^oC) \end{equation*}The efficiency of the turbine is then given in terms of % of maximum isentropic efficiency. For example, with an isentropic efficiency of 80%:

\begin{equation*} W_\text{real} = 0.8 \, W_\text{isentropic} \end{equation*}And, by an energy balance, we can calculate the ending enthalpy of the liquid runoff at 0.75 kPa, $H_\text{real liquid}$:

\begin{equation*} W_\text{real} = H_\text{steam}(T_\text{coolant}, P_\text{reactor}) - H_\text{real liquid}(T_\text{runoff}, 0.75 kPa) \end{equation*}With the ending enthalpy of the liquid runoff known, and its pressure fixed by the output pressure of the turbine (0.75 kPa), the runoff temperature of the liquid from the turbine can be calculated using the IAPWS97 backwards solver functions.

The condenser takes the runoff from the turbine and condenses it to a subcooled liquid at the current operating pressure of the reactor. The work required to do this is calculated as follows:

\begin{equation*} W_\text{condenser} = \frac{H_\text{turbine runoff}(T_\text{runoff}, 0.0035 kPa) - H_\text{subcooled liquid out}(T_\text{subcooled}, P_\text{reactor})}{efficiency}) \end{equation*}Where the enthalpy of the subcooled liquid out is given by:

\begin{equation*} H_\text{subcooled liquid out}(T_\text{subcooled}, P_\text{reactor}) = H_\text{saturated liquid}(T_\text{saturated}, P_\text{reactor}) * (subcooling) \end{equation*}With subcooling being the % by which the enthalpy of the subcooled liquid is lowered from the enthalpy it would have if it were at its boiling point. This % is constant, and is based on the % of subcooling that exists at steady state with a reactor pressure of 7 mPa and a condenser runoff temperature of 220 degrees celcius:

\begin{equation*} subcooling = (1 -\frac{H(T_{saturated}, 7 mPa, x=0) - H(220 C, 7 mPa}{H(T_{saturated}, 7 mPa, x=0}) \end{equation*}In [1]:

```
# Import various packages; must have the Cortix repository installed
import matplotlib.pyplot as plt
import scipy.constants as const
# Leave this block here for Azure
try:
import cortix
except ImportError:
print('Installing the "cortix" package...')
print('')
!pip install cortix
import cortix
from cortix.src.cortix_main import Cortix
from cortix.src.network import Network
from cortix.examples.bwr.bwr import BWR
from cortix.examples.bwr.turbine import Turbine
from cortix.examples.bwr.condenser import Condenser
```

In [2]:

```
# Setup parameters
```

In [3]:

```
# Create the application network
community = Community(n_groups=n_groups, non_offender_adult_population=100, offender_pool_size=0)
city_net.module(community)
community.end_time = end_time
community.time_step = time_step
community.show_time = (True,10*const.day)
community.save = True
prison = Prison(n_groups=n_groups, pool_size=0)
city_net.module(prison)
prison.end_time = end_time
prison.time_step = time_step
prison.save = True
parole = Parole(n_groups=n_groups, pool_size=0)
city_net.module(parole)
parole.end_time = end_time
parole.time_step = time_step
parole.save = True
adjudication = Adjudication(n_groups=n_groups, pool_size=0)
city_net.module(adjudication)
adjudication.end_time = end_time
adjudication.time_step = time_step
adjudication.save = True
jail = Jail(n_groups=n_groups, pool_size=0)
city_net.module(jail)
jail.end_time = end_time
jail.time_step = time_step
jail.save = True
arrested = Arrested(n_groups=n_groups, pool_size=0)
city_net.module(arrested)
arrested.end_time = end_time
arrested.time_step = time_step
arrested.save = True
probation = Probation(n_groups=n_groups, pool_size=0)
city_net.module(probation)
probation.end_time = end_time
probation.time_step = time_step
probation.save = True
city_net.connect( prison, parole,'bidirectional' )
city_net.connect( adjudication, prison )
city_net.connect( jail, prison )
city_net.connect( adjudication, jail )
city_net.connect( arrested, jail )
city_net.connect( arrested, adjudication )
city_net.connect( arrested, probation )
city_net.connect( probation, jail )
city_net.connect( adjudication, probation )
city.save = True
city_net.connect( prison, community )
city_net.connect( jail, community )
city_net.connect( arrested, community, 'bidirectional' )
city_net.connect( parole, community )
city_net.connect( adjudication, community )
city_net.connect( probation, community )
```

In [4]:

```
# View the Cortix network created
city_net.draw()
```

Out[4]:

In [5]:

```
# Run the simulation!
city.run()
```

In [6]:

```
community = city_net.modules[0]
prison = city_net.modules[1]
parole = city_net.modules[2]
adjudication = city_net.modules[3]
jail = city_net.modules[4]
arrested = city_net.modules[5]
probation = city_net.modules[6]
quant_names = {'Prison':'fpg','Parole':'feg','Adjudication':'fag',
'Jail':'fjg','Arrested':'frg','Probation':'fbg','Community':'f0g'}
total_num_params = 0
```

In [7]:

```
'''Inspect Data Function'''
def inspect_module_data(module,quant_names):
population_phase = module.population_phase
if isinstance(quant_names,dict):
(fxg_quant, time_unit) = population_phase.get_quantity_history(quant_names[module.name])
elif isinstance(quant_names,str):
(fxg_quant, time_unit) = population_phase.get_quantity_history(quant_names)
fxg_quant.plot( x_scaling=1/const.day, x_label='Time [day]', y_label=fxg_quant.name+' ['+fxg_quant.unit+' %]')
plt.grid()
plt.show()
# number of parameters in the prison model
n_params = (len(population_phase.GetActors())-1)*n_groups
print('# parameters = ',n_params)
return n_params
```

In [8]:

```
total_num_params += inspect_module_data(prison,quant_names)
```

In [9]:

```
total_num_params += inspect_module_data(parole,quant_names)
```

In [10]:

```
total_num_params += inspect_module_data(adjudication,quant_names)
```

In [11]:

```
total_num_params += inspect_module_data(jail,quant_names)
```

In [12]:

```
total_num_params += inspect_module_data(arrested,quant_names)
```

In [13]:

```
total_num_params += inspect_module_data(probation,quant_names)
```

Offenders removed (negative sign) from the general non-offender community.

In [14]:

```
total_num_params += inspect_module_data(community,quant_names)
```

In [15]:

```
total_num_params += inspect_module_data(community,'f0g_free')
```

In [16]:

```
'''Total number of unknowns and parameters'''
print('total number of unknowns =', n_groups*len(city_net.modules))
print('total number of parameters =', total_num_params)
```