$\newcommand{\Amtrx}{\boldsymbol{\mathsf{A}}} \newcommand{\Bmtrx}{\boldsymbol{\mathsf{B}}} \newcommand{\Cmtrx}{\boldsymbol{\mathsf{C}}} \newcommand{\Mmtrx}{\boldsymbol{\mathsf{M}}} \newcommand{\Imtrx}{\boldsymbol{\mathsf{I}}} \newcommand{\Pmtrx}{\boldsymbol{\mathsf{P}}} \newcommand{\Qmtrx}{\boldsymbol{\mathsf{Q}}} \newcommand{\Lmtrx}{\boldsymbol{\mathsf{L}}} \newcommand{\Umtrx}{\boldsymbol{\mathsf{U}}} \newcommand{\xvec}{\boldsymbol{\mathsf{x}}} \newcommand{\yvec}{\boldsymbol{\mathsf{y}}} \newcommand{\zvec}{\boldsymbol{\mathsf{z}}} \newcommand{\avec}{\boldsymbol{\mathsf{a}}} \newcommand{\bvec}{\boldsymbol{\mathsf{b}}} \newcommand{\cvec}{\boldsymbol{\mathsf{c}}} \newcommand{\rvec}{\boldsymbol{\mathsf{r}}} \newcommand{\uvar}{\boldsymbol{u}} \newcommand{\fvar}{\boldsymbol{f}} \newcommand{\norm}[1]{\bigl\lVert{#1}\bigr\rVert} \DeclareMathOperator{\rank}{rank} \DeclareMathOperator{\abs}{abs}$ SciPy 2020: A. Rotker and V. F. de Almeida 06Jul20

# Balance of Plant Network Simulation for a Boiling Water Nuclear Reactor¶

Austin Rotker and Valmor F. de Almeida
Dept. of Chemical Engineering (Nuclear Program)
University of Massachusetts Lowell, USA

## Objectives¶

• Use one-group energy, two-temperature, point-reactor dynamics modeling (six-group delayed neutron emitters) to obtain insight on neutron density behavior with time-dependent reactivity of a Boiling Water Reactor (BWR) during startup, steady-state operation and shutdown.
• Explore and understand intrinsic temperature reactivity coefficient feedback from nuclear sources.
• Develop and apply computational components based on the Python ecosystem. This application uses Cortix as the framework for developing a dynamical network for basic elements of a BWR balance of plant. The BoP BWR application is developed as an example of a Cortix application.

## Problem Statement¶

Consider a boiling water reactor, with UO$_2$ fuel and a light water moderator and coolant. Given the neutron generation time $\ell$, delayed neutron fraction, $\beta$, decay constants of 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 varying neutron reactivity, $\rho(t)$.

The neutron reactivity variation with time will depend indirectly on the effects of the balance of plant around the reactor including startup and shutdown transients. This will include a high pressure turbine, two low pressure turbines and a condenser (Fig. 1).

Fig. 1: Balance of plant layout: Reactor, 1; high-pressure turbine, 2; two low pressure turbines, 3 and 3*, and condenser is 4.

Steam is created from nuclear heat generated in the reactor and passed to the high-pressure turbine, which expands the steam to about 0.5 MPA (or 5 bar). The runoff from the high-pressure turbine is then split in half and fed to the two low-pressure turbines, which expand the runoff to 0.5 bar. The runoff is then fed to a condenser which produces subcooled liquid water, which is then fed back into the bottom of the reactor.

Modeling the entire network is accomplished through the use of the network simulation library Cortix.. The individual plant components (reactor, turbines and condenser) are implemented as Cortix modules, with data exchange handled by the Cortix library.

Additionally, a number of the calculations used in this model rely on the IAPWS library developed by Juan José Gómez Romera (jjgomera). IAPWS is the International Association for the Properties of Water and Steam, which develops power series that can be used to calculate properties of water and steam (such as specific enthalpy and entropy) at a variety of temperatures, pressures and qualities. The IAPWS library devloped by jjgomera implements these series as python functions which can be easily accessed and called to calculate any property required for thermodynamic calculations, drastically simplifying the process of modelling the two-phase systems in use by the turbine and condenser modules. This project specifically uses the IAPWS97 module within the library, which uses the formulations published in 1997 by IAPWS foundation as an update to ones developed in 1967.

## Boiling Water Reactor Model¶

Our model for the BWR module is one-group (no diffusion within the moderator assumed for leakage), tracks two temperatures (the uranium fuel and the water which acts as both a coolant and a moderator), and also behaves like a point reactor. This means that all of the kinetics, neutron density and heat generation is collapsed to a single point. Temperatures for the fuel and coolant are therefore averages of what they would be in a real reactor system.

A fission event releases heat and several neutrons; these neutrons go on to create more fission events, while the heat generated from the fissions is absorbed by the fuel. As the temperature of the fuel rises, it releases heat to the lower-temperature coolant, causing it to boil and produce steam. This steam then leaves off the top of the reactor, removing some heat from the system to be absorbed by the turbines and condensers later on down the line.

### One-Group Energy Neutron Balance¶

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 (with $\rho(t)$ being the reactivity, $\beta$ being the delayed neutron fraction and $\ell$ being the prompt neutron lifetime), 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 (with $\lambda$ being the decay constant of the specific delayed neutron precursor group, $c$ being its concentration and $q(t)$ being the magnitude of the neutron source, if any). 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)$ (with $\beta$ being the yield of each delayed neutron precursor group from fission), 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 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_m(t) - T^{(0)}\bigr). \end{equation*}

Where $\rho(0)$ is the initial reactivity insertion into the reactor, $\alpha_{T_N}$ is the overall temperature reactivity coefficient as a function of current temperature of the moderator $T_m(t)$, and $T^{(0)}$ is the reference temperature (298 $^{\circ}K$).

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 $T_f$ is the temperature of the fuel, $\rho_f$ and $c_{pf}$ is its density and heat capacity, $q^{'''}(t)$ is the heat generated per unit volume within the fuel, $\dot Q_f$ is the rate of heat absorption from the fuel by the moderator and $V_f$ is the total volume of the fuel.

The rate of heat absorption by the moderator is calculated by:

\begin{equation*} {\dot{Q}_\text{f}} = -UA\bigl(T_f-T_c\bigr), \end{equation*}

Where $U$ and $A$ is the overall heat transfer coefficient and heat transfer area, $T_f$ is the temperature of the fuel and $T_c$ is the temperature of the moderator.

The rate of heat generation per unit volume within the fuel is equal to:

\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*}

Where $G_f$ is the energy released per fission, $\sigma_\text{fo}$ is the microscopic fission cross section within the fuel at a reference temperature $T^\text{(o)}$ of 298 $^{\circ} K$, $T$ is the current temperature of the fuel, $\epsilon$ is the fast fission factor, $\rho_\text{fm}$ is the mass density of the fuel, $M_\text{fn}$ is the molar mass of the fuel, $N_A$ is the number density of the fuel, $n(t)$ is the neutron density at that point in time and $v_\text{(o)}$ is the velocity of a thermal neutron.

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*}

$T_c$ is the temperature of the coolant leaving the top of the reactor, $\tau_c$ is the coolant residence time within the reactor, $T_{in}$ is the temperature of the subcooled coolant entering the bottom of the reactor, $\rho_c$ is the density of the coolant, $c_{pc}$ is the heat capacity of the coolant and $V_c$ is the volume of the coolant. The heat rate source is the negative of the heat rate sink in the fuel/moderator.

Neutron density and delayed neutron precursor concentrations 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*}

### Vector ODE System¶

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)$ is the state vector 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),\\ u_8(t)&=T_f(t),\\ u_9(t)&=T_c(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 + q(t), \\ 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, \\ f_8 & = - \frac{1}{\rho_f\,c_{pf}}\biggl(q^{'''}(t) - \frac{\dot{Q}_\text{f}}{V_f} \biggr), \\ f_9 & = -\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{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},\\ u_8(0)&=T_{f_0},\\ u_9(0)&=T_{c_0} \end{align*}

where $n_0$ and $c_{i_0}$ are initial values for the unknowns given by the problem specification.

### Reactivity Feedback Model¶

The moderator temperature reactivity coefficient is the major contributor to the power coefficient in a boiling water reactor. The power coefficient determines the feedback effect of increasing the power generated by the reactor (and therefore neutron density as well as heat generation). Generally, we see that increasing power will lead to an increase in the temperature of the water used to moderate the reactor. This will reduce the density of the water, and therefore its effectiveness as a moderator. We would expect that this would lead to lower reactivity, and the equation developed below reflects this. We can estimate the numerical effects of moderator temperature on reactivity by using a linearization of the effect of increasing temperature on the neutron multiplecation factor (K), which is used to derive reactivity.

\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*}

## Turbine Model¶

A boiling water reactor will typically use several different turbines, either designed to expand high-pressure steam drawn directly off the top of the reactor, or to expand lower-pressure runoff steam drawn from the high-pressure turbines. This is typically done to improve thermodynamic efficiency and avoid condensation of the steam within the turbines themselves. Our model contains a single high-pressure turbine which expands steam coming off the top of the reactor to around 5 bar; There are also two low-pressure turbines which expand the runoff from the high-pressure turbine to 0.05 bar. Both of these different types of turbine use the same model but different parameters to determine runoff temperature and the amount of power generated.

We start with a basic heat balance on the turbine:

\begin{equation*} Q + W_\text{s} = \bigl(\Delta H + \Delta K + \Delta Z\bigr)\,\dot m \end{equation*}

With Q being the heat into or out of the system, 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 on a per-kilogram basis, respectively. $\dot m$ is the mass flowrate of fluid through the turbine. 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 \dot m \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 \dot m = \bigl(H_\text{outlet} - H_\text{inlet}\bigr)\,\dot m \end{equation*}

Which can then be re-arranged as:

\begin{equation*} W_\text{s} = \Delta H \dot m = \bigl(H_\text{inlet} - H_\text{outlet}\bigr)\,\dot m \end{equation*}

If conditions of the inlet steam are given, either being the same temperature and pressure of the steam leaving the top of the reactor (in the case of a high-pressure turbine), or having the runoff temperature and pressure of steam leaving a high-pressure turbine (in the case of a low-pressure turbine). The exit pressure of the steam leaving the turbine is similarly fixed by the design pressure of the turbine ($P_\text{turbine}$, 5 bar in the case of a high-pressure turbine, 0.05 bar in the case of a low-pressure turbine). Therefore, the above equation can be rewritten as:

\begin{equation*} W_\text{s} = \bigl(H_\text{inlet}(T_\text{inlet}, P_\text{inlet}) - H_\text{outlet}(T_\text{runoff}, P_\text{turbine})\bigr)\,\dot m \end{equation*}

It should be noted that the turbines are not assumed to operate isentropically, but instead at around 80% of the maximum isentropic work for any given steam inlet properties.

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 $P_\text{turbine}$ and a temperature which would give the same enthalpy as steam of temperature $T_\text{inlet}$ and pressure $P_\text{inlet}$.

\begin{equation*} W_\text{isentropic} = \bigl(H_\text{inlet}(T_\text{inlet}, P_\text{inlet}) - H_\text{outlet}(T_\text{isentropic} ,P_\text{turbine} )\bigr)\dot m \end{equation*}

Actual work performed by the turbine can then be determined by simply multiplying the known efficiency by the calculated isentropic work.

\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 $P_\text{turbine}$, $H_\text{real runoff}$:

\begin{equation*} H_\text{real runoff}(T_\text{real}, P_\text{turbine}) = \frac{W_\text{real}}{\dot m} - H_\text{inlet}(T_\text{inlet}, P_\text{inlet}) \end{equation*}

With the ending enthalpy of the liquid runoff known, and its pressure fixed by the output pressure of the turbine, the runoff temperature and quality of the liquid from the turbine can be calculated using the IAPWS97 backwards solver functions. To do this, $H_\text{real runoff}$ is first compared to the dew and bubble point enthalpies of water at the runoff pressure to arrive at the runoff quality:

\begin{equation*} X_\text{runoff} = \frac{H_\text{real runoff} - H_\text{bubl pt}}{H_\text{dew pt} - H_\text{bubl pt}} \end{equation*}

If $X_\text{runoff}$ is greater than 1, then the runoff will be a superheated vapor. If it is less than 0, then it will be a subcooled liquid. In either case, the runoff temperature can then be calculated by finding the superheated vapor or liquid with the same enthalpy $H_\text{real runoff}$ and the turbine runoff pressure, using an IAPWS backwards solver function. If $X_\text{runoff}$ is between 0 and 1, this indicates that the runoff is a two-phase mixture with a temperature fixed by the boiling point of the liquid at the specified runoff pressure.

## Condenser Model¶

The condenser takes the runoff from the low-pressure turbines and condenses it into a subcooled liquid to be fed back into the bottom of the reactor. In theory, the runoff from the turbines can be a two-phase mixture, a bubble-point liquid, a dew-point vapor or even superheated vapor. In practice, however, the runoff from the low-pressure turbines is only a two-phase mixture of varying composition for the duration of the simulation.

The first thing the condenser model calculates is the amount of energy which must be withdrawn from the turbine runoff to cool it to its dew point. This is calculated by:

\begin{equation*} Q = \bigl(H_\text{real runoff} - H_\text{bubl pt}\bigr)\,\dot m \end{equation*}

This will be used to determine how much of the condenser's cooling capacity is used to bring the incoming liquid down to the bubble point temperature, and how much of this capacity will be left to subcool it. With the actual condensation duty known, an accurate number for cooling capacity at the current input composition and pressure must be determined. We start by determining the log-mean temperature difference for the condensation, which is calculated by:

\begin{equation*} \Delta T_\text{lm} = \frac{T_\text{c out} - T_\text{c in}}{\ln\Bigl(\frac{T_\text{critical} - T_\text{c in}}{T_\text{boiling pt} - T_\text{c out}}\Bigr)} \end{equation*}

Where $T_\text{c out}$ is the temperature of the cooling water leaving the condensation stage, $T_\text{c in}$ is the temperature of the cooling water entering the condensation stage and $T_\text{boiling pt}$ is the boiling point temperature of the turbine runoff at the turbine runoff pressure. With $T_\text{c in}$ and $Q$ known, the runoff temperature of the cooling water can be calculated from the following equation:

\begin{equation*} T_\text{c out} = T_\text{c in } + \frac{Q}{\dot m \, c_\text{p}}, \end{equation*}

where here $\dot m$ is the cooling water mass flowrate, not the mass flowrate leaving the turbines. $c_\text{p}$ is the heat capacity of the coolant in this case. The resulting $T_\text{c out}$ may then be plugged into the expression for log mean temperature difference. In order to arrive at an accurate overall heat transfer coefficient for the condensation operation, the McNaught expression is used

\begin{equation*} \alpha_\text{sh} = \alpha_\text{l} \Bigl(\frac{1}{x_\text{tt}}\Bigr)^{0.78}, \end{equation*}

where $\alpha_\text{sh}$ is the overall heat transfer coefficient, $\alpha_\text{l}$ is the heat transfer coefficient for a liquid water mixture at the same temperature and pressure as the turbine runoff flowing over a tube (calculated using the Churchill-Bernstein correlation), and $x_\text{tt}$ is the Martinelli parameter. The Martinelli parameter can be calculated using the following expression:

\begin{equation*} x_\text{tt}^2 = \Bigl(\frac{1 - x}{x}\Bigr)^{1.8}\Bigl(\frac{\rho_\text{dew pt}}{\rho_\text{bubbl pt}}\Bigr)\Bigl(\frac{\eta_\text{bubbl pt}}{\eta_\text{dew pt}}\Bigr)^{0.5}, \end{equation*}

where x is the quality of the steam entering the turbine, $\rho$ is the mass density of either the dew point or bubble point mixture at the turbine runoff pressure, and $\eta$ is the viscocity of either the dew point or bubble point mixture at the turbine runoff pressure. With the overall heat transfer coefficient, $Q$ and the temperature driving force $\Delta T_\text{lm}$ known, we can calculate the amount of cooling capacity left to subcool the condensate by determining the amount of heat transfer area that hasn't been used:

\begin{equation*} A = \frac{Q}{\alpha_\text{sh} \Delta T_\text{lm}} \end{equation*}

The total area (and therefore capacity) remaining for subcooling may then be found by subtracting the above result from the total heat transfer area of the condenser.

Subcooling proceeds slightly differently than the algorithm for determining remaining capacity after condensation. It is an iterative process which begins with a guess for the temperature of the cooling water leaving the condenser. With this guess, a preliminary value for the runoff temperature of the process fluid from the condenser may be generated by performing a heat balance between the process fluid (turbine runoff) and the cooling water (note that for this section "c in" refers to the temperature of the coolant leaving the condensation stage and "c out" refers to the runoff temperature of the coolant leaving the condenser itself):

\begin{equation*} T_\text{process out} = T_\text{boiling pt} - \frac{\dot m_\text{coolant}\, c_\text{p coolant} \, \Delta T_\text{coolant}}{\dot m_\text{process}\, c_\text{p process}} \end{equation*}

Where $\Delta T_\text{coolant}$ is the difference between the temperature of the coolant after the condensation step and the guessed value for the temperature of the coolant leaving the condenser (that is, $T_\text{c runoff} - T_\text{c out}$. With $T_\text{process out}$ known, a log-mean temperature difference for this cooling operation may be calculated (assuming parallel flow):

\begin{equation*} \Delta T_\text{lm} = \frac{\bigl(T_\text{boiling pt} - T_\text{c out}\bigr) - \bigl(T_\text{process out} - T_\text{c in}\bigr)}{\ln\Bigl(\frac{T\text{boiling pt} - T_\text{c out}}{T_\text{process out} - T_\text{c in}}\Bigr)} \end{equation*}

An overall heat transfer coefficient for this iteration may then be obtained using the Churchill-Bernstein correlation. Taking into account the remaining area, the overall heat transfer coefficient and the temperature driving force, a Q may then be calculated:

\begin{equation*} Q = \alpha \, A_\text{remaining} \Delta T_\text{lm} \end{equation*}

With this Q, we may loop back around and calculate a new 'guess' value for the runoff temperature of the coolant by another energy balance:

\begin{equation*} T_\text{c out} = T_\text{c in} + \frac{Q}{\dot m_\text{coolant} \, c_\text{p coolant}} \end{equation*}

With a new value for $T_\text{c out}$ obtained that is closer to the actual value, the algorithm may be repeated until it converges on a suitably accurate value (for the purposes of our code, within +/- 0.01 $^{\circ}$C. When the value for $T_\text{c out}$ converges, we will be left with a suitable estimation for the temperature of the subcooled process fluid leaving the condenser and being fed back into the bottom of the reactor.

## Cortix Balance of Plant Network¶

In [2]:
'''Balance of Plant Header'''

# Leave this block here for Azure
try:
import cortix, iapws
except ImportError:
print('Installing missing packages...')
print('')
!pip install cortix iapws
import cortix, iapws

from cortix import Cortix
from cortix import Network

from cortix.examples.bwr.reactor import BWR
from cortix.examples.bwr.turbine import Turbine
from cortix.examples.bwr.condenser import Condenser
from cortix.examples.bwr.params import get_params

import scipy.constants as unit
import matplotlib.pyplot as plt

end_time = 1 * unit.hour
unit.second = 1.0
time_step = 30.0 * unit.second
show_time = (True,5*unit.minute)

use_mpi = False  # True for MPI; False for Python multiprocessing

plant = Cortix(use_mpi=use_mpi, splash=True) # System top level

plant_net = plant.network = Network() # Network

[8128] 2020-07-06 01:44:00,348 - cortix - INFO - Created Cortix object
_____________________________________________________________________________
L A U N C H I N G
_____________________________________________________________________________
...                                        s       .     (TAAG Fraktur)
xH88"~ .x8X                                 :8      @88>
:8888   .f"8888Hf        u.      .u    .      .88      %8P      uL   ..
:8888>  X8L  ^""   ...ue888b   .d88B :@8c    :888ooo    .     [email protected]  @88R
X8888  X888h        888R Y888r ="8888f8888r -*8888888  [email protected]  ""Y888k/"*P
88888  !88888.      888R I888>   4888>"88"    8888    888E    Y888L
88888   %88888      888R I888>   4888> "      8888      888E      8888
88888 > 8888>     888R I888>   4888>        8888      888E      888N
8888L %  ?888   ! u8888cJ888   .d888L .+    .8888Lu=   888E   .u./"888&
8888  -*""   /   "*888*P"    ^"8888*"     ^%888*     888&  d888" Y888*"
"888.      :"      "Y"          "Y"         "Y"      R888"  "Y   Y"
""***~"                                           ""
https://cortix.org
_____________________________________________________________________________

In [3]:
'''Reactor'''

params = get_params()

reactor = BWR(params)  # Create reactor module

reactor.name = 'BWR'
reactor.save = True
reactor.time_step = time_step
reactor.end_time = end_time
reactor.show_time = show_time

plant_net.module(reactor)  # Add reactor module to network

In [4]:
'''Turbine 1'''

params['turbine_inlet_pressure'] = 2
params['turbine_outlet_pressure'] = 0.5
params['high_pressure_turbine'] = True

turbine1 = Turbine(params)  # turbine 1 module

turbine1.name = 'High Pressure Turbine'
turbine1.save = True
turbine1.time_step = time_step
turbine1.end_time = end_time

plant_net.module(turbine1)  # Add turbine 1 module to network

In [5]:
'''Turbine 2'''

params['turbine_inlet_pressure'] = 0.5
params['turbine_outlet_pressure'] = 0.005
params['high_pressure_turbine'] = False
params['steam flowrate'] = params['steam flowrate'] / 2

turbine2 = Turbine(params)  # turbine 2 module

turbine2.name = 'Low Pressure Turbine 1'
turbine2.save = True
turbine2.time_step = time_step
turbine2.end_time = end_time

plant_net.module(turbine2)  # Add turbine 2 module to network

In [6]:
'''Turbine 3'''

params['turbine_inlet_pressure'] = 0.5
params['turbine_outlet_pressure'] = 0.005
params['high_pressure_turbine'] = False

turbine3 = Turbine(params)  # turbine 3 module

turbine3.name = 'Low Pressure Turbine 2'
turbine3.save = True
turbine3.time_step = time_step
turbine3.end_time = end_time

plant_net.module(turbine3)

In [7]:
'''Condenser'''

params['steam flowrate'] = params['steam flowrate'] * 2

condenser = Condenser(params)

condenser.name = 'Condenser'
condenser.save = True
condenser.time_step = time_step
condenser.end_time = end_time

plant_net.module(condenser)

In [8]:
'''Balance of Plant Network Connectivity'''

plant_net.connect([reactor, 'coolant-outflow'], [turbine1, 'inflow'])
plant_net.connect([turbine1, 'outflow-1'], [turbine2, 'inflow'])
plant_net.connect([turbine1, 'outflow-2'], [turbine3, 'inflow'])
plant_net.connect([turbine2, 'outflow-1'], [condenser, 'inflow-1'])
plant_net.connect([turbine3, 'outflow-1'], [condenser, 'inflow-2'])
plant_net.connect([condenser, 'outflow'], [reactor, 'coolant-inflow'])

plant_net.draw()

Out[8]:
In [9]:
'''Run'''

plant.run()  # Run network dynamics simulation

plant.close()  # Properly shutdow plant

[8128] 2020-07-06 01:45:09,487 - cortix - INFO - Launching Module <cortix.examples.bwr.reactor.BWR object at 0x7f58c0186b90>
[8128] 2020-07-06 01:45:09,490 - cortix - INFO - Launching Module <cortix.examples.bwr.turbine.Turbine object at 0x7f58c01ccd50>
[8128] 2020-07-06 01:45:09,494 - cortix - INFO - Launching Module <cortix.examples.bwr.turbine.Turbine object at 0x7f58a0fa3850>
[8128] 2020-07-06 01:45:09,500 - cortix - INFO - Launching Module <cortix.examples.bwr.turbine.Turbine object at 0x7f58c0186690>
[8128] 2020-07-06 01:45:09,491 - cortix - INFO - BWR::run():time[m]=0.0
[8128] 2020-07-06 01:45:09,506 - cortix - INFO - Launching Module <cortix.examples.bwr.condenser.Condenser object at 0x7f58c0186750>
[8128] 2020-07-06 01:45:45,095 - cortix - INFO - BWR::run():time[m]=5.0
[8128] 2020-07-06 01:46:22,822 - cortix - INFO - BWR::run():time[m]=10.0
[8128] 2020-07-06 01:46:54,866 - cortix - INFO - BWR::run():time[m]=15.0
[8128] 2020-07-06 01:47:26,538 - cortix - INFO - BWR::run():time[m]=20.0
[8128] 2020-07-06 01:48:04,010 - cortix - INFO - BWR::run():time[m]=25.0
[8128] 2020-07-06 01:48:54,637 - cortix - INFO - BWR::run():time[m]=30.0
[8128] 2020-07-06 01:49:31,600 - cortix - INFO - BWR::run():time[m]=35.0
[8128] 2020-07-06 01:49:37,574 - cortix - INFO - BWR::run():time[m]=40.0
[8128] 2020-07-06 01:49:37,819 - cortix - INFO - BWR::run():time[m]=45.0
[8128] 2020-07-06 01:49:38,061 - cortix - INFO - BWR::run():time[m]=50.0
[8128] 2020-07-06 01:49:38,306 - cortix - INFO - BWR::run():time[m]=55.0
[8128] 2020-07-06 01:49:38,559 - cortix - INFO - BWR::run():time[m]=60.0
[8128] 2020-07-06 01:49:38,595 - cortix - INFO - run()::Elapsed wall clock time [s]: 338.25
[8128] 2020-07-06 01:49:38,595 - cortix - INFO - Closed Cortix object.
_____________________________________________________________________________
T E R M I N A T I N G
_____________________________________________________________________________
...                                        s       .     (TAAG Fraktur)
xH88"~ .x8X                                 :8      @88>
:8888   .f"8888Hf        u.      .u    .      .88      %8P      uL   ..
:8888>  X8L  ^""   ...ue888b   .d88B :@8c    :888ooo    .     [email protected]  @88R
X8888  X888h        888R Y888r ="8888f8888r -*8888888  [email protected]  ""Y888k/"*P
88888  !88888.      888R I888>   4888>"88"    8888    888E    Y888L
88888   %88888      888R I888>   4888> "      8888      888E      8888
88888 > 8888>     888R I888>   4888>        8888      888E      888N
8888L %  ?888   ! u8888cJ888   .d888L .+    .8888Lu=   888E   .u./"888&
8888  -*""   /   "*888*P"    ^"8888*"     ^%888*     888&  d888" Y888*"
"888.      :"      "Y"          "Y"         "Y"      R888"  "Y   Y"
""***~"                                           ""
https://cortix.org
_____________________________________________________________________________
[8128] 2020-07-06 01:49:38,596 - cortix - INFO - close()::Elapsed wall clock time [s]: 338.25


## Results¶

During startup, the neutron density follows a predictable sigmoid shape (Fig. 2), rapidly shooting up and then slowing down until steady state is reached (1.0 being the reactor steady state neutron density). Steady-state is maintained for about 20 minutes until the reactor enters shutdown and neutron density gradually drops to zero.

Neutron density also drops slightly about 15 minutes after steady state is reached. High moderator temperature results in negative reactivity insertion and a corresponding drop in neutron density.

Fig. 2: Neutron density time-variation during startup and shutdown.

The six delayed neutron emitter groups follow a similar pattern to neutron density during startup and shutdown (Fig. 3), the difference in levels being due to the different $\beta_i$ yields given in the one-group neutron energy balance. This plot (Fig. 3) shows the relative abundance of each of the delayed neutron precursor groups, in relation to steady state reactor neutron density. The small drop in density observed in the neutron density around the 25 minute mark can also be observed here.

Fig. 3: Delayed emitter concentrations variation with time.

Fuel temperature (Fig. 4) follows the neutron density curve fairly closely, which is to be expected since higher neutron density will generate more nuclear heat. The starting and ending temperatures are both room temperature at 297 $^{o}K$.

Fig. 4:Average fuel temperature variation with time.

Coolant temperature (Fig. 5) follows the neutron density curve fairly closely, which is to be expected since more nuclear heat generation will result in a higher runoff coolant temperature. The starting and ending temperatures are both room temperature at 297 $^{o}K$.

Fig. 5: Average temperature of the steam drawn off the top of the reactor, varying with time.

The power produced by the high-pressure turbine (Fig. 6) follows the same general curve seen with outflow coolant temperature (Fig. 5), with the exception of the sharp drop observed in power generation almost immeditately during shutdown. The turbine is tripped as soon as inlet temperature drops below a certain threshold, causing the quick drop in power within the span of a few minutes during shutdown.

Fig. 6: Power produced by the high-pressure turbine, varying with time.

The low-pressure turbines (Fig. 7) follow the same power curves as the high-pressure turbines, with the exception that they produce slightly lower amounts of power at steady state.

Fig. 7: Power produced by a single low-pressure turbine, varying with time.

Water exits the condenser at room temperature during startup, and the runoff temperature gradually increases to around 303 $^{o}K$ as reactor temperature increases and higher-energy steam enters the balance of plant. The observed decrease in outflow temperature occurs at the same time as the drop in neutron density (Fig. 2).

Fig. 8: Temperature of the subcooled water leaving the condenser and being fed back into the bottom of the reactor, varying with time.

## Conclusions¶

The evolution of the plant from startup to shutdown is carried out by message passing (either using the Message Passing Interface for Python or native Python multiprocessing pipes) of time-dependent data specified by the connectivity of the Cortix network. Each module on the network runs on its own process or thread and data is sychronized by message passing through ports of the modules.

The plant network model results follow a sigmoid shape up to the steady state point which would be expected in real plant operation. Some of the individual steady state values (i.e. the power produced by the turbines) are qualitatively similar to real operation, but improvement of the parameters used could make the model quantitatively enhanced.

The Cortix library is instrumental in modeling dynamic systems with multiple modules. Further work could be done to scale the system into a much larger more complex plant, including additional turbines and unit operations such as pumping or safety systems.

## Acknowledgements¶

This work was partially funded by the Francis College of Engineering at UMass Lowell, Department of Chemical Engineering (Nuclear Program), and the Cortix group.