Skip to content

Commit 16221e3

Browse files
milanofthemilanofthe
authored
Merge pull request #23 from pathsim/feature/point-kinetics
point kinetics
2 parents aade215 + bc5b399 commit 16221e3

File tree

6 files changed

+446
-0
lines changed

6 files changed

+446
-0
lines changed

‎docs/source/examples/point_kinetics.ipynb‎

Lines changed: 152 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,152 @@
1+
{
2+
"cells": [
3+
{
4+
"cell_type": "markdown",
5+
"metadata": {},
6+
"source": [
7+
"# Reactor Point Kinetics\n",
8+
"\n",
9+
"Simulating the transient neutron population in a nuclear reactor using the point kinetics equations (PKE) with six delayed neutron precursor groups."
10+
]
11+
},
12+
{
13+
"cell_type": "markdown",
14+
"metadata": {},
15+
"source": [
16+
"The point kinetics model couples the neutron density $n$ with $G$ delayed neutron precursor groups $C_i$:\n",
17+
"\n",
18+
"$$\\frac{dn}{dt} = \\frac{\\rho - \\beta}{\\Lambda} \\, n + \\sum_{i=1}^{G} \\lambda_i \\, C_i + S$$\n",
19+
"\n",
20+
"$$\\frac{dC_i}{dt} = \\frac{\\beta_i}{\\Lambda} \\, n - \\lambda_i \\, C_i \\qquad i = 1, \\ldots, G$$\n",
21+
"\n",
22+
"where $\\beta = \\sum_i \\beta_i$ is the total delayed neutron fraction, $\\Lambda$ is the prompt neutron generation time, and $\\rho$ is the reactivity."
23+
]
24+
},
25+
{
26+
"cell_type": "code",
27+
"execution_count": null,
28+
"metadata": {},
29+
"outputs": [],
30+
"source": "import matplotlib.pyplot as plt\n\nplt.style.use('../pathsim_docs.mplstyle')\n\nfrom pathsim import Simulation, Connection\nfrom pathsim.blocks import Source, Scope\nfrom pathsim.solvers import GEAR52A\n\nfrom pathsim_chem.neutronics import PointKinetics"
31+
},
32+
{
33+
"cell_type": "markdown",
34+
"metadata": {},
35+
"source": "The point kinetics system is very stiff — the prompt neutron generation time $\\Lambda \\sim 10^{-5}\\,\\text{s}$ creates eigenvalues on the order of $10^5$. The variable-order BDF solver GEAR52A is ideal here: it requires only one implicit solve per step and adapts both step size and order to the smooth exponential dynamics. The default fixed-point tolerance (`1e-9`) is unnecessarily tight for this problem; relaxing it to `1e-6` gives a ~70x speedup with negligible loss in accuracy.\n\n## 1. Delayed Supercritical Step\n\nInsert a step reactivity of $\\rho = 0.003$ (about $0.46\\beta$). Since $\\rho < \\beta$, the reactor is delayed supercritical — the power rises on a slow time scale governed by the delayed neutrons."
36+
},
37+
{
38+
"cell_type": "code",
39+
"execution_count": null,
40+
"metadata": {},
41+
"outputs": [],
42+
"source": "reactor = PointKinetics(n0=1.0)\n\nrho_step = 0.003\nsrc_rho = Source(lambda t: rho_step if t > 0 else 0.0)\nsrc_s = Source(lambda t: 0.0)\nsco = Scope(labels=[\"n\"])\n\nsim = Simulation(\n [src_rho, src_s, reactor, sco],\n [\n Connection(src_rho, reactor), # rho\n Connection(src_s, reactor[1]), # S\n Connection(reactor, sco), # n\n ],\n dt=0.01,\n Solver=GEAR52A,\n tolerance_fpi=1e-6,\n log=True,\n)\n\nsim.run(100)"
43+
},
44+
{
45+
"cell_type": "code",
46+
"execution_count": null,
47+
"metadata": {},
48+
"outputs": [],
49+
"source": [
50+
"sco.plot(lw=1.5)\n",
51+
"plt.yscale('log')\n",
52+
"plt.title('Delayed Supercritical (ρ = 0.003)')\n",
53+
"plt.show()"
54+
]
55+
},
56+
{
57+
"cell_type": "markdown",
58+
"metadata": {},
59+
"source": [
60+
"The neutron density rises exponentially on a time scale of seconds — much slower than the prompt neutron lifetime ($\\Lambda \\sim 10^{-5}$ s) because the delayed neutrons control the dynamics when $\\rho < \\beta$."
61+
]
62+
},
63+
{
64+
"cell_type": "markdown",
65+
"metadata": {},
66+
"source": [
67+
"## 2. Prompt Supercritical\n",
68+
"\n",
69+
"Insert $\\rho = 0.008 > \\beta \\approx 0.0065$. Now the reactor is prompt supercritical — the power rises on the prompt neutron time scale, producing a rapid excursion."
70+
]
71+
},
72+
{
73+
"cell_type": "code",
74+
"execution_count": null,
75+
"metadata": {},
76+
"outputs": [],
77+
"source": "reactor = PointKinetics(n0=1.0)\n\nrho_prompt = 0.008\nsrc_rho = Source(lambda t: rho_prompt if t > 0 else 0.0)\nsrc_s = Source(lambda t: 0.0)\nsco = Scope(labels=[\"n\"])\n\nsim = Simulation(\n [src_rho, src_s, reactor, sco],\n [\n Connection(src_rho, reactor),\n Connection(src_s, reactor[1]),\n Connection(reactor, sco),\n ],\n dt=0.001,\n Solver=GEAR52A,\n tolerance_fpi=1e-6,\n log=True,\n)\n\nsim.run(0.5)"
78+
},
79+
{
80+
"cell_type": "code",
81+
"execution_count": null,
82+
"metadata": {},
83+
"outputs": [],
84+
"source": [
85+
"sco.plot(lw=1.5)\n",
86+
"plt.yscale('log')\n",
87+
"plt.title('Prompt Supercritical (ρ = 0.008 > β)')\n",
88+
"plt.show()"
89+
]
90+
},
91+
{
92+
"cell_type": "markdown",
93+
"metadata": {},
94+
"source": [
95+
"The power rises orders of magnitude within milliseconds. This is why prompt criticality must be avoided in reactor design — the delayed neutrons are the key safety mechanism that keeps power transients manageable."
96+
]
97+
},
98+
{
99+
"cell_type": "markdown",
100+
"metadata": {},
101+
"source": [
102+
"## 3. Subcritical with External Source\n",
103+
"\n",
104+
"A subcritical assembly ($\\rho = -0.05$) with a constant external neutron source. The system reaches an equilibrium where the source multiplication produces a steady neutron population:\n",
105+
"\n",
106+
"$$n_{ss} = \\frac{-S \\, \\Lambda}{\\rho}$$"
107+
]
108+
},
109+
{
110+
"cell_type": "code",
111+
"execution_count": null,
112+
"metadata": {},
113+
"outputs": [],
114+
"source": "rho_sub = -0.05\ns_ext = 1e5\nLambda = 1e-5\n\nreactor = PointKinetics(n0=0.001, Lambda=Lambda)\n\nsrc_rho = Source(lambda t: rho_sub)\nsrc_s = Source(lambda t: s_ext)\nsco = Scope(labels=[\"n\"])\n\nsim = Simulation(\n [src_rho, src_s, reactor, sco],\n [\n Connection(src_rho, reactor),\n Connection(src_s, reactor[1]),\n Connection(reactor, sco),\n ],\n dt=0.01,\n Solver=GEAR52A,\n tolerance_fpi=1e-6,\n log=True,\n)\n\nsim.run(50)"
115+
},
116+
{
117+
"cell_type": "code",
118+
"execution_count": null,
119+
"metadata": {},
120+
"outputs": [],
121+
"source": [
122+
"n_ss = -s_ext * Lambda / rho_sub\n",
123+
"\n",
124+
"sco.plot(lw=1.5)\n",
125+
"plt.axhline(n_ss, color='r', ls='--', lw=1, label=f'$n_{{ss}}$ = {n_ss:.4f}')\n",
126+
"plt.legend()\n",
127+
"plt.title('Subcritical with Source (ρ = −0.05)')\n",
128+
"plt.show()"
129+
]
130+
},
131+
{
132+
"cell_type": "markdown",
133+
"metadata": {},
134+
"source": [
135+
"The neutron density converges to the expected source-multiplied equilibrium value, confirming the subcritical multiplication physics."
136+
]
137+
}
138+
],
139+
"metadata": {
140+
"kernelspec": {
141+
"display_name": "Python 3",
142+
"language": "python",
143+
"name": "python3"
144+
},
145+
"language_info": {
146+
"name": "python",
147+
"version": "3.11.0"
148+
}
149+
},
150+
"nbformat": 4,
151+
"nbformat_minor": 4
152+
}

‎src/pathsim_chem/__init__.py‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -16,3 +16,4 @@
1616
from .tritium import *
1717
from .thermodynamics import *
1818
from .process import *
19+
from .neutronics import *

‎src/pathsim_chem/neutronics/__init__.py‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1 @@
1+
from .point_kinetics import *

‎src/pathsim_chem/neutronics/point_kinetics.py‎

Lines changed: 146 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,146 @@
1+
#########################################################################################
2+
##
3+
## Reactor Point Kinetics Block
4+
##
5+
#########################################################################################
6+
7+
# IMPORTS ===============================================================================
8+
9+
import numpy as np
10+
11+
from pathsim.blocks.dynsys import DynamicalSystem
12+
13+
# CONSTANTS =============================================================================
14+
15+
# Keepin U-235 thermal 6-group delayed neutron data
16+
BETA_U235 = [0.000215, 0.001424, 0.001274, 0.002568, 0.000748, 0.000273]
17+
LAMBDA_U235 = [0.0124, 0.0305, 0.111, 0.301, 1.14, 3.01]
18+
19+
# BLOCKS ================================================================================
20+
21+
class PointKinetics(DynamicalSystem):
22+
"""Reactor point kinetics equations with delayed neutron precursors.
23+
24+
Models the time-dependent neutron population in a nuclear reactor using
25+
the point kinetics equations (PKE). The neutron density responds to
26+
reactivity changes through prompt fission and delayed neutron emission
27+
from G precursor groups.
28+
29+
Mathematical Formulation
30+
-------------------------
31+
The state vector is :math:`[n, C_1, C_2, \\ldots, C_G]` where :math:`n`
32+
is the neutron density (or power) and :math:`C_i` are the delayed neutron
33+
precursor concentrations.
34+
35+
.. math::
36+
37+
\\frac{dn}{dt} = \\frac{\\rho - \\beta}{\\Lambda} \\, n
38+
+ \\sum_{i=1}^{G} \\lambda_i \\, C_i + S
39+
40+
.. math::
41+
42+
\\frac{dC_i}{dt} = \\frac{\\beta_i}{\\Lambda} \\, n
43+
- \\lambda_i \\, C_i \\qquad i = 1, \\ldots, G
44+
45+
where :math:`\\beta = \\sum_i \\beta_i` is the total delayed neutron
46+
fraction and :math:`\\Lambda` is the prompt neutron generation time.
47+
48+
Parameters
49+
----------
50+
n0 : float
51+
Initial neutron density [-]. Default 1.0 (normalized).
52+
Lambda : float
53+
Prompt neutron generation time [s].
54+
beta : array_like
55+
Delayed neutron fractions per group [-].
56+
lam : array_like
57+
Precursor decay constants per group [1/s].
58+
"""
59+
60+
input_port_labels = {
61+
"rho": 0,
62+
"S": 1,
63+
}
64+
65+
output_port_labels = {
66+
"n": 0,
67+
}
68+
69+
def __init__(self, n0=1.0, Lambda=1e-5,
70+
beta=None, lam=None):
71+
72+
# defaults
73+
if beta is None:
74+
beta = BETA_U235
75+
if lam is None:
76+
lam = LAMBDA_U235
77+
78+
beta = np.asarray(beta, dtype=float)
79+
lam = np.asarray(lam, dtype=float)
80+
81+
# input validation
82+
if Lambda <= 0:
83+
raise ValueError(f"'Lambda' must be positive but is {Lambda}")
84+
if len(beta) != len(lam):
85+
raise ValueError(
86+
f"'beta' and 'lam' must have the same length "
87+
f"but got {len(beta)} and {len(lam)}"
88+
)
89+
if len(beta) == 0:
90+
raise ValueError("'beta' and 'lam' must have at least one group")
91+
92+
# store parameters
93+
G = len(beta)
94+
beta_total = float(np.sum(beta))
95+
96+
self.n0 = n0
97+
self.Lambda = Lambda
98+
self.beta = beta
99+
self.lam = lam
100+
self.G = G
101+
self.beta_total = beta_total
102+
103+
# initial conditions: steady state at rho = 0
104+
x0 = np.empty(G + 1)
105+
x0[0] = n0
106+
for i in range(G):
107+
x0[1 + i] = beta[i] / (Lambda * lam[i]) * n0
108+
109+
# rhs of point kinetics ode system
110+
def _fn_d(x, u, t):
111+
n = x[0]
112+
C = x[1:]
113+
114+
rho = u[0]
115+
S = u[1] if len(u) > 1 else 0.0
116+
117+
dn = (rho - beta_total) / Lambda * n + np.dot(lam, C) + S
118+
119+
dC = beta / Lambda * n - lam * C
120+
121+
return np.concatenate(([dn], dC))
122+
123+
# jacobian of rhs wrt state [n, C_1, ..., C_G]
124+
def _jc_d(x, u, t):
125+
rho = u[0]
126+
127+
J = np.zeros((G + 1, G + 1))
128+
129+
J[0, 0] = (rho - beta_total) / Lambda
130+
for i in range(G):
131+
J[0, 1 + i] = lam[i]
132+
J[1 + i, 0] = beta[i] / Lambda
133+
J[1 + i, 1 + i] = -lam[i]
134+
135+
return J
136+
137+
# output function: neutron density only
138+
def _fn_a(x, u, t):
139+
return x[:1].copy()
140+
141+
super().__init__(
142+
func_dyn=_fn_d,
143+
jac_dyn=_jc_d,
144+
func_alg=_fn_a,
145+
initial_value=x0,
146+
)

‎tests/neutronics/__init__.py‎

Whitespace-only changes.

0 commit comments

Comments
 (0)