Theory & Method

Point Kinetics Simulator · AntraSIM · version 1.0.0 · build 2026-09-19. This document describes the model as implemented in this application. Where the two source notebooks (Euler and RK4, v1.0.0) disagree with each other or with their own written text, the executable model, physical reasoning and the references below have been used, and the resolution is documented in §17.

How to read this document. Each topic follows the same progression: Intuition → Physical mechanism → Mathematical formulation → Engineering interpretation → Assumptions and limitations. A newcomer can stop after the first two layers of each section; an experienced engineer will find the equations, units, numerical implementation and limits in the last three.

Model status. Established point-kinetics theory (§1–§5, §14) is combined with a deliberately simplified first-order feedback model (§7–§8) that is not a fuel heat-conduction or two-phase thermal-hydraulic model. The application is numerically verified (§16). It is not validated against plant or experimental transient data and is not safety-grade.

1. Physical Background

Intuition A chain reaction is a population of neutrons in which every neutron causes, on average, a certain number of further fissions. If exactly one new neutron survives per neutron lost the population stays constant (critical); if slightly more, it grows; if slightly fewer, it decays.

Physical mechanism The balance is expressed by the effective multiplication factor $k$ , the ratio of neutrons produced in one generation to neutrons lost in the preceding generation. The departure from criticality is the reactivity $ρ$ [3]:

ρ = \frac{k - 1}{k}

Engineering interpretation Reactivity is dimensionless ( $Δ k / k$ ). Because it is a small number in normal operation it is usually quoted in pcm (per cent mille): $1 pcm = 10^{- 5} Δ k / k$ . Reactivity is also compared with the delayed-neutron fraction $β$ : the ratio $ρ / β$ is expressed in “dollars” (not used in the interface, which shows pcm). All reactivity readouts in this application are in pcm; all internal calculations use $Δ k / k$ .

The causal chain implemented by the simulator is:

External reactivity ρ_ext(t)→neutron flux φ and precursors C_i→fuel temperature T_f and void v (first-order lag)→Doppler ρ_D and void ρ_v→total reactivity ρ→back to the flux

2. Point-Reactor Approximation

Intuition A real core has a flux shape that varies from place to place. The point-reactor model treats the whole core as a single “point” whose neutron population rises and falls as one body, with a fixed shape.

Physical mechanism The time-dependent neutron flux is assumed to separate into a fixed spatial shape and a time-dependent amplitude, $ϕ (𝐫, E, t) \approx n (t) ψ (𝐫, E)$ . Weighting the space- and energy-dependent balance equations with the adjoint flux and integrating removes the spatial and energy dependence and leaves ordinary differential equations for the amplitude and the precursor populations.

Engineering interpretation The state variable $ϕ$ is a measure of the neutron population. In the point-kinetics literature the same equations are written for the neutron density $n (t)$ [neutrons/cm³] with precursor concentrations $C_{i}$ [nuclei/cm³] [11, 12], or for power-like quantities with all states in MW [13]; the forms differ only by constant factors. This application solves the power-normalised form and converts to flux and density for display:

Quantity	Relation	Default at 100 %
Fission power $P = ϕ$ [MWth] (1 a.u. = 1 MWth)	$ϕ_{0} = P_{rated} P (0) / 100 %$	3000 MWth
Neutron flux $Φ$ [n/(cm²·s)]	$Φ = ϕ Φ_{rated} / P_{rated}$	3.0·10¹³
Neutron density $n$ [n/cm³]	$n = Φ / v_{0}$ , $v_{0} = 2.2 \cdot 10^{5}$ cm/s	1.364e+08
Precursor concentration $C_{i}$ [atoms/cm³]	$C_{i}^{atoms} = C_{i}^{MW} Φ_{rated} / (P_{rated} v_{0})$	see §5

Here $Φ$ is the conventional thermal flux $Φ = n v_{0}$ with the 2200 m/s reference speed of thermalised neutrons (the value 2200 m/s at 20 °C is quoted in [13a]); $Φ_{rated} = 3.0 \cdot 10^{13}$ n/(cm²·s) is the core-average thermal flux at rated power. Core-average thermal fluxes of a light-water power reactor of $2.2 \cdot 10^{13}$ (typical PWR value) and up to $3.6 \cdot 10^{13}$ n/(cm²·s) (in-core facility described as comparable to a commercial LWR) were found in [9, 10]; 3·10¹³ lies within that range. The normalisation is a pure scaling of the state vector: the kinetics equations are linear and homogeneous in $(ϕ, C_{i})$ and the feedback gains are defined per percent of rated power (§7, §8), so the relative behaviour does not depend on $P_{rated}$ , $Φ_{rated}$ or $v_{0}$ . This was checked by rescaling the default case from the 1000 a.u. of the source notebooks to 3000 MWth: fuel temperature, void and reactivities are reproduced to about 10⁻¹⁴ and the power is three times larger to round-off.

Consistency of the scaling constants. The prompt lifetime $T_{n} = 4 \cdot 10^{- 5}$ s and $n (0)$ imply a neutron production rate $n / T_{n} = 3.4 \cdot 10^{12}$ neutrons/(cm³·s). With about 2.4 neutrons and about 200 MeV (3.2·10⁻¹¹ J) per fission (general values, not looked up in this release) this corresponds to a volumetric power of about 45 W/cm³ = 45 kW/L, which is of the order of the core power density of a boiling-water reactor (a few tens of kW/L; general knowledge, not separately verified). The three constants $Φ_{rated}$ , $v_{0}$ and $T_{n}$ are therefore mutually consistent to within the accuracy of these general values. Changing one of them without the others breaks this consistency; the application does not check it.

Assumptions and limitations of the normalisation $ϕ$ represents the fission (neutron-generated) power only. Heat from the decay of fission products (decay heat) and the delay between fission power and heat transfer to the coolant are not modelled; after a SCRAM the real thermal power therefore falls more slowly than the $ϕ$ shown. The single conversion constant between power, flux and density ignores spatial peaking and the neutron energy spectrum. The 3000 MWth reference core is BWR-like (core-average void 45 %, §8) and is a set of representative values, not data of a specific plant.

3. Prompt and Delayed Neutrons

Intuition Nearly all fission neutrons appear immediately (“prompt”). A small fraction, of order 0.65 % for thermal fission of U-235 [3], is emitted seconds to a minute later by radioactive decay of certain fission products (“precursors”). This small delayed fraction slows the reactor response from milliseconds to seconds and is what makes reactor control possible.

Physical mechanism Some fission products (e.g. Br-87 with a half-life of about 56 s) decay by β-emission to a nucleus excited above its neutron binding energy, which then emits a neutron. The many precursors are grouped by half-life into six effective groups [1, 3]. Each group is characterised by a delayed fraction $β_{i}$ and a decay constant $λ_{i}$ ; the total is $β = \sum_{i} β_{i}$ .

Engineering interpretation Time scales in the default data set:

Quantity	Value	Meaning
Prompt neutron lifetime $T_{n}$	4.0·10⁻⁵ s	Time scale of the prompt response (default retuned, §15)
Precursor mean lives $1 / λ_{i}$	80.6, 32.8, 9.01, 3.32, 0.877, 0.332 s	Time scale of the delayed groups
Mean delayed-neutron lifetime $\sum β_{i} / λ_{i} / β$	13.04 s	Average delay of a delayed neutron
Effective generation time incl. delayed neutrons $T_{n} + \sum β_{i} / λ_{i}$	0.0848 s	Why the reactor responds in seconds rather than milliseconds

Assumptions and limitations Group constants depend on the fissioning nuclide and on neutron energy; the default set is for thermal fission of U-235 (§15). Different evaluations of the six-group parameters exist and differ noticeably, especially in the shortest-lived groups [6]. The simulator treats one nuclide only; it does not model a fuel mixture or burn-up dependence of β.

4. Six-Group Point Kinetics Equations

Intuition Two things happen in every instant: the neutron population is multiplied (or lost) at a rate set by how far reactivity is from the delayed-neutron fraction, and the precursor “reservoirs” are filled by fissions and drained by decay, feeding the population back with a delay.

Mathematical formulation The model implemented is

\frac{d ϕ}{d t} = \frac{ρ - β}{T_{n}} ϕ + \sum_{i = 1}^{6} λ_{i} C_{i} + S

\frac{d C_{i}}{d t} = \frac{β_{i}}{T_{n}} ϕ - λ_{i} C_{i}, i = 1, \dots, 6

Symbol	Meaning	Unit
$ϕ$	neutron flux, normalised to fission power (§2)	MWth (1 a.u. = 1 MWth)
$C_{i}$	precursor group state; same scaling as $ϕ$ (MWth-equivalent), shown as concentration in atoms/cm³ (§2)	MWth-eq. (atoms/cm³)
$ρ$	total reactivity, $ρ = ρ_{ext} + ρ_{D} + ρ_{v}$ (§9)	Δk/k (shown in pcm)
$β_{i}, β$	delayed-neutron fraction of group $i$ , and its sum	–
$λ_{i}$	precursor decay constant	s⁻¹
$T_{n}$	prompt neutron lifetime (generation time)	s
$S$	external neutron source	MWth/s

Engineering interpretation In the first equation the term $(ρ - β) ϕ / T_{n}$ is the net prompt production: fission neutrons that are prompt, minus losses. The sum $\sum λ_{i} C_{i}$ is the delayed-neutron source, and $S$ an independent source. In the second equation each group is fed in proportion to the flux and drains with its own decay constant. Because $β / T_{n}$ is a large rate (163 s⁻¹ for the default data) the prompt response is fast; when $ρ < β$ the delayed neutrons hold the population back and the response is governed by the slow precursor dynamics. When $ρ \geq β$ prompt neutrons alone sustain the chain reaction (§13, §14).

Assumptions and limitations The precursor states are number densities of the delayed-neutron precursor nuclei of group $i$ : in the density form of the equations $d C_{i} / d t = β_{i} n / T_{n} - λ_{i} C_{i}$ the term $β_{i} n / T_{n}$ is the production rate of precursors per cm³ [11, 12], so that $C_{i}$ is in atoms/cm³ when $n$ is in neutrons/cm³. In the power-normalised form solved here the same $C_{i}$ carries the scaling of $ϕ$ (MWth-equivalent); §2 gives the conversion. The equations are linear in $(ϕ, C_{i})$ for a given $ρ$ ; the nonlinearity of the model enters only through the feedback (§7–§9).

5. Initial Equilibrium Conditions

Intuition The simulation starts from a reactor that has been running steadily: the precursors are in balance with the flux, so nothing changes until an external disturbance occurs.

Mathematical formulation Setting $d C_{i} / d t = 0$ gives the equilibrium precursor state. The application initialises

C_{i} (0) = \frac{β_{i}}{λ_{i} T_{n}} ϕ_{0}

which is the density-form expression $C_{i} (0) = β_{i} n_{0} / (λ_{i} Λ)$ of [12]. For the default case, $n_{0} = 1.364 e + 08$ n/cm³ gives

Group	β_i	λ_i [s⁻¹]	C_i(0) [atoms/cm³]	C_i(0) [MWth-eq.]
1	0.000215	0.0124	5.911e+10	1.3e+06
2	0.001424	0.0305	1.592e+11	3.502e+06
3	0.001274	0.111	3.913e+10	8.608e+05
4	0.002568	0.301	2.908e+10	6.399e+05
5	0.000748	1.14	2.237e+09	4.921e+04
6	0.000273	3.01	3.092e+08	6802

Inserting the equilibrium state into the flux equation gives $d ϕ / d t = ρ ϕ_{0} / T_{n} + S$ . For the default case ( $ρ (0) = 0$ , $S = 0$ ) the derivative vanishes and the initial state is a true steady state; the value is computed from the equilibrium expression and never stored as a constant.

Non-default initial conditions. This initialisation is an equilibrium only if $ρ (0) = 0$ and $S = 0$ . With $S > 0$ the flux initially grows as $d ϕ / d t = S$ ; a true source-driven subcritical steady state ( $ρ < 0$ ) would require $ϕ = - S T_{n} / ρ$ with the same precursor relation, which is not what is initialised. Likewise, an initial fuel temperature or void different from the reference values gives $ρ (0) \neq 0$ through the feedback terms. The application reports such cases as informational messages; the resulting initial transient is a correct solution of the equations for those initial conditions, not a numerical artefact.

6. External Reactivity Model

Intuition The user “operates” the reactor by prescribing an external reactivity history: two steps and a SCRAM (rapid shutdown).

Mathematical formulation The external reactivity is a stateless, piecewise-constant, right-continuous function of time:

ρ_{ext} (t) = {\begin{matrix} 0 & t < t_{1} \\ ρ_{1} & t_{1} \leq t < t_{2} \\ ρ_{1} + ρ_{2} & t_{2} \leq t < t_{SCRAM} \\ ρ_{1} + ρ_{2} + ρ_{SCRAM} & t \geq t_{SCRAM} \end{matrix}

with the ordering constraint $0 \leq t_{1} \leq t_{2} \leq t_{SCRAM} \leq t_{end}$ . Defaults: $ρ_{1} = ρ_{2} = + 100$ pcm at 10 s and 20 s, $ρ_{SCRAM} = - 6000$ pcm (insertion) at 40 s, $t_{end} = 60$ s.

Choice of the SCRAM value. A SCRAM must overcome the power defect: when power and temperatures fall, the feedback reactivities change sign and add positive reactivity (for the parameter set of §7–§8 about +3700 pcm: void 70 pcm/% × 45 % = 3150 pcm, Doppler 1.6 pcm/°C × 340 °C ≈ 540 pcm). The total worth of control rods is design-specific and is quoted as about 6000–9000 pcm, with a PWR power defect of about 2500 pcm [9a]. The default insertion of −6000 pcm (lower end of that range; net −5800 pcm after the +200 pcm of the steps) leaves a shutdown margin of about 2100 pcm in the model, above the 0.38 % Δk/k stuck-rod margin quoted as a BWR licensing standard [9b]. The −1000 pcm net level of the notebooks is smaller than the power defect of a boiling-water core and would let the reactor return to criticality during the void collapse. The steps $ρ_{1} = ρ_{2} = + 100$ pcm are small compared with $β = 650$ pcm (sub-prompt-critical) and are retained from the notebooks. The SCRAM is a step: a real rod insertion takes seconds, and the prompt drop shown is therefore idealised.

Engineering interpretation All three reactivity values are increments that are added to the preceding level, so the SCRAM is handled exactly like the two steps: after $t_{SCRAM}$ the external reactivity is $ρ_{1} + ρ_{2} + ρ_{SCRAM}$ . A SCRAM value of 0 pcm therefore changes nothing, and $- 6000$ pcm represents the worth of all control rods being inserted on top of whatever the earlier steps have done (net $- 5800$ pcm with the default steps). The function is evaluated from the model parameters alone. The source notebooks instead treated the SCRAM as an absolute level (−1000 pcm replacing the accumulated value); with the present convention the notebook net value corresponds to $ρ_{SCRAM} = - 1200$ pcm. The source notebooks stored the first step in a variable that was assigned only inside one branch; with $t_{1} = t_{2}$ that variable was never set, and a re-run silently used a stale value. The production model has no such state: coincident events are merged deterministically (the value after the last coincident event applies), and an event at $t = 0$ is in force from the first instant.

Assumptions and limitations Ideal instantaneous reactivity steps; no rod-motion dynamics, no insertion rate limits, no protection logic. All reactivity values are user inputs in pcm, converted with $1 pcm = 10^{- 5} Δ k / k$ .

7. Fuel-Temperature Feedback

Intuition When power rises the fuel heats up. Hotter fuel absorbs more neutrons in the U-238 resonances (Doppler broadening), which lowers reactivity and pushes the power back down: an inherent, fast-acting negative feedback.

Physical mechanism Increased fuel temperature broadens the capture resonances and raises resonance absorption without a corresponding rise in fission. In a real core the fuel temperature follows the power through heat conduction and transfer to the coolant, on time scales of seconds. The temperature that matters is an effective (weighted) average over the pellet; this application uses one average fuel temperature, in °C (temperature differences in °C and K are identical).

Mathematical formulation — model actually implemented The simulator does not solve heat conduction. It uses a first-order lag toward a power-dependent target temperature:

T_{f, target} = T_{ref} + K_{T} (ϕ - ϕ_{0}), \frac{d T_{f}}{d t} = \frac{T_{f, target} - T_{f}}{τ_{T}}

ρ_{D} = α_{D} (T_{f} - T_{ref})

Symbol	Meaning	Default	Unit
$T_{ref}, T_{f} (0)$	average fuel temperature at rated power (reference and initial value; zero of $ρ_{D}$ )	627 (= 900 K)	°C
$K_{T}$	steady-state rise of the average fuel temperature per percent of rated power above the reference power $ϕ_{0}$ (interface value $K_{T}^{%}$ ; model value $K_{T} = K_{T}^{%} / (P_{rated} / 100)$ in °C/MWth)	3.4 (0.113 °C/MWth)	°C per % rated
$τ_{T}$	time constant of the lag	5.0	s
$α_{D}$	Doppler reactivity coefficient	−1.6 pcm/°C (−1.6·10⁻⁵ Δk/k/°C)	pcm/°C

The interface shows $α_{D}$ in pcm/°C and displays the equivalent Δk/k per °C beneath the field; internally $α_{D} = (pcm/°C) \cdot 10^{- 5}$ .

Reasonableness of the constants (retuned in this release).

$T_{ref} = 627$ °C. An average fuel temperature of 900 K is used for a typical PWR in the neutronic study [15]; a similar value is assumed here for a BWR-like core. It is a representative assumption, not plant data.
$α_{D} = - 1.6$ pcm/°C. The background section of the patent [7] quotes $1.6 \cdot 10^{- 5} Δ k / k /^{\circ} C$ as the Doppler coefficient of boiling-water reactors in operation. The notebook value −1 pcm/K was of the right order but somewhat low. (The patent value is the only core-level BWR figure checked.)
$K_{T} = 3.4$ °C per % rated. In a boiling-water core the coolant stays at the saturation temperature (about 286 °C at roughly 7 MPa; a steam-table value that was not separately verified in this release), so the fuel-to-coolant temperature difference is what changes with power: $K_{T} \approx (627 - 286) / 100 \approx 3.4$ °C per percent. This is a derived assumption (linearised through the rated point), not a measured value; the notebook value (1 K per %) would give an unrealistically small fuel temperature change.
$τ_{T} = 5$ s. No numerical value of the fuel time constant was found in the checked sources (the BWR stability patent [7a] only states that the void response follows the power after a delay set by the fuel thermal time constant). An own order-of-magnitude estimate from generic UO₂ properties (conductivity about 3 W/m·K, volumetric heat capacity about 3·10⁶ J/m³·K, pellet radius about 4 mm; general values, not looked up) gives about 3 s for the fundamental conduction mode plus 1–2 s for gap, cladding and film, i.e. 4–5 s. The value 5 s (unchanged from the notebooks) is therefore plausible but not verified.

Engineering interpretation At a sustained power $ϕ$ , the fuel settles at $T_{f} = T_{ref} + K_{T} (ϕ - ϕ_{0})$ with time constant $τ_{T}$ , so the steady Doppler feedback per percent of rated power above $ϕ_{0}$ is $α_{D} K_{T}^{%} = - 1.6 \cdot 3.4 = - 5.4$ pcm per percent. With the default power peak of about 3626 MWth (120.9 % of rated) the fuel rises by about 19 °C and $ρ_{D}$ reaches about −31 pcm. $K_{T}$ and $τ_{T}$ are modelling parameters that represent, in one number each, the combined heat capacity, heat-transfer and power-to-temperature relation; they do not correspond to a specific fuel design.

Correction relative to the notebooks. The notebook text writes $d T_{f} / d t = k_{T} (ϕ - ϕ_{0})$ , a pure integrator. The executable code of both notebooks instead implements the lag toward the target above. The lag form is used here because it is what the simulations actually computed, and because it has a physical equilibrium: with the integrator form a power that stays above $ϕ_{0}$ would heat the fuel without bound.

Assumptions and limitations A single lumped fuel temperature; no fuel-to-coolant coupling, no heat capacity or gap conductance, no temperature-dependent Doppler coefficient. The temperature has no meaning beyond the stated proportionality; at low power it extrapolates linearly toward $T_{ref} - 100 K_{T}^{%}$ (287 °C at zero power).

8. Moderator-Void Feedback

Intuition In a boiling-water reactor more power means more steam bubbles in the coolant. Steam is a much poorer moderator than water, so fewer neutrons are slowed to fission energies and reactivity drops: a strong, fairly quick negative feedback. At rated power the core-average void fraction of a BWR is of the order of 40–45 %.

Mathematical formulation — model actually implemented With void fraction $v$ (a fraction between 0 and 1):

v_{target} = clip [v_{ref} + K_{v} (ϕ - ϕ_{0}), v_{min}, v_{max}], \frac{d v}{d t} = \frac{v_{target} - v}{τ_{v}}

ρ_{v} = α_{v} (v - v_{ref})

Defaults: $v (0) = v_{ref} = 0.45$ (45 % core-average void at rated power), $K_{v}^{%} = 0.45$ % void per percent of rated power ( $K_{v} = K_{v}^{%} / P_{rated} = 1.5 \cdot 10^{- 4}$ void fraction per MWth), $τ_{v} = 6.0$ s, $α_{v} = - 0.7 Δ k / k$ per unit void fraction (= −70 pcm per % void), $v_{min} = 0$ , $v_{max} = 0.7$ .

Units — handle with care. The interface uses percent void for all void quantities and pcm per percentage point for the coefficient. The conversions are:

Quantity	Interface unit	Internal unit	Default, both forms
Void fraction $v$	% void	fraction	v = 0.01 ⇔ 1 % void; v(0) = 0.45 ⇔ 45 %
Coefficient $α_{v}$	pcm per % void	Δk/k per unit void fraction	−70 pcm/% = −0.7 Δk/k per unit
Gain $K_{v}$	% void per % rated power	void fraction per MWth	0.45 %/% = 1.5·10⁻⁴ per MWth
Limits $v_{min}, v_{max}$	% void	fraction	0 … 70 % = 0 … 0.7

The derivation of the default coefficient: $- 0.7 Δ k / k$ per unit void fraction × 0.01 (one percentage point) = $- 0.007 Δ k / k = - 70$ pcm per percent void.

Reasonableness of the constants (retuned in this release).

$v (0) = 45$ %. Specified for the 3000 MWth reference core; a mean void fraction of about 40 % at design flow (rising to 56 % at reduced flow) is described for a BWR in the patent [8].
$α_{v} = - 70$ pcm/% void. The patent [7] gives $- 7.0 \cdot 10^{- 4} Δ k / k$ per % void for boiling-water reactors in operation; a search hit for a benchmark study quotes about 30 pcm/% void for a typical BWR [6a] (abstract only). The notebook value of −200 pcm/% void was 3–7 times larger than these BWR figures (the range −30 to −250 pcm/% found in [9c] is a PWR statement and not applicable). −70 pcm/% is used; a value near −30 pcm/% is a sensible sensitivity case.
$K_{v} = 0.45$ % per %. The void response to power is nonlinear (boiling starts above the subcooled region); the default is the secant slope 45 %/100 % of a void fraction that vanishes at zero power and reaches 45 % at rated power, linearised through the operating point. It is a modelling assumption. A consequence is that the void target reaches its lower limit only at zero power, so the clip is inactive in normal transients (§11).
$τ_{v} = 6$ s. Boiling follows the heat flux, which lags the fission power by the fuel time constant; the void then responds with the coolant transit time (of the order of 1 s). The estimate $τ_{v} \approx τ_{T} + 1$ s is an own assumption with qualitative support in [7a]; it is not a verified value.
$v_{max} = 70$ %. The 0.5 limit of the notebook manual leaves only 5 percentage points above 45 %, which a 10 % power increase would exhaust ( $45 + 0.45 \cdot 10 > 50$ ). The limit is therefore raised to 70 %, above the 56 % mean void mentioned in [8]. This is a modelling bound, not a sourced physical limit.

Engineering interpretation The steady void feedback per percent of rated power above $ϕ_{0}$ is $α_{v} K_{v}^{%} = - 70 \cdot 0.45 = - 31.5$ pcm, about six times the Doppler value (−5.4 pcm per percent); the total power coefficient of the model is −36.9 pcm per percent. In the default transient the void rises to 47.4 % ( $ρ_{v} \approx - 170$ pcm at the peak). After the SCRAM the void collapses toward 4 %, giving $ρ_{v} \approx + 2900$ pcm, which together with the Doppler term (+500 pcm) is the power defect that the SCRAM value must overcome (§6).

Corrections relative to the notebooks. (a) The notebook text describes an integrator $d v / d t = k_{v} (ϕ - ϕ_{0})$ ; the code implements the lag above, which is used here. (b) The notebook manual states the model limit $0 \leq v \leq 0.5$ , while the executable code clamps the target to $[0, 0.5]$ but the state to $[0, 1]$ . The production model applies one range $[v_{min}, v_{max}]$ (default 0–70 %, see above) consistently: with the target inside the range and $Δ t \leq τ_{v}$ a first-order lag cannot leave it; if the numerical solution nevertheless leaves it (for instance with a very large step or user-modified limits) the application does not clamp the state silently but reports “Model left the defined void-fraction range”. (c) The notebooks did not subtract a reference void in $ρ_{v}$ ; this is now $α_{v} (v - v_{ref})$ , identical for $v_{ref} = 0$ .

Assumptions and limitations The void is a single lumped variable driven by power, not by coolant flow, pressure, subcooling or channel geometry; there is no transport delay and no two-phase flow model, and no feedback from the void to the flow. The clip makes the right-hand side non-smooth wherever the target reaches $v_{min}$ or $v_{max}$ (for example with $K_{v}^{%} > 0.45$ after a SCRAM); this limits the accuracy of higher-order integration (§11, §16).

9. Coupled State Model

Intuition Neutrons, precursors, fuel temperature and void all influence one another; the simulator therefore treats them as one coupled system and advances them together.

Coupling structure of the implemented model. The flux φ drives two lagged states (fuel temperature, void). Their reactivity contributions ρ_D and ρ_v are added to the external reactivity in the summing junction Σ and act back on the point-kinetics equations. All three signals enter with a plus sign; the feedback is negative because α_D and α_v are negative for the default data.

Mathematical formulation One canonical state vector and one derivative function serve both solvers:

𝐲 = {(ϕ, C_{1}, \dots, C_{6}, T_{f}, v)}^{T}, \frac{d 𝐲}{d t} = 𝐟 (t, 𝐲)

ρ (t, 𝐲) = ρ_{ext} (t) + α_{D} (T_{f} - T_{ref}) + α_{v} (v - v_{ref})

The components of $𝐟$ are the right-hand sides of §4, §7 and §8 with this $ρ$ . The reactivity is evaluated from the state of each solver stage, so the feedback is part of the numerical scheme and not frozen during a step.

Engineering interpretation The total reactivity, its three contributions and the derived quantities (§14) are computed from the stored state; no quantity is integrated twice or separately from this system.

10. Explicit Euler Method

Intuition Advance the solution along the current slope for one time step.

Mathematical formulation

𝐲_{n + 1} = 𝐲_{n} + h 𝐟 (t_{n}, 𝐲_{n})

Engineering interpretation The method is first order: halving the step halves the error (verified in §16, observed order 1.0 in the asymptotic range). The reactivity used in the flux equation is that of the state $𝐲_{n}$ itself. The Euler notebook computed the new reactivity but used the previous array element in the flux derivative, an accidental one-step delay of the feedback; that delay is not reproduced.

Assumptions and limitations Stability on the fast prompt mode requires $h | ρ - β | / T_{n} ≲ 2$ (§13). Euler is provided as a comparison and educational solver; RK4 is the production default.

11. Runge–Kutta Fourth-Order Method

Intuition Sample the slope at the start, twice in the middle and at the end of the step, and combine them with weights 1, 2, 2, 1.

Mathematical formulation The classical scheme applied to the complete coupled state vector:

𝐤_{1} = 𝐟 (t_{n}, 𝐲_{n}), 𝐤_{2} = 𝐟 (t_{n} + \frac{h}{2}, 𝐲_{n} + \frac{h}{2} 𝐤_{1}), 𝐤_{3} = 𝐟 (t_{n} + \frac{h}{2}, 𝐲_{n} + \frac{h}{2} 𝐤_{2}), 𝐤_{4} = 𝐟 (t_{n} + h, 𝐲_{n} + h 𝐤_{3})

𝐲_{n + 1} = 𝐲_{n} + \frac{h}{6} (𝐤_{1} + 2 𝐤_{2} + 2 𝐤_{3} + 𝐤_{4})

Engineering interpretation In every stage the flux, precursors, fuel temperature and void are those of the stage, and the feedback reactivities are recomputed from them. The RK notebook instead advanced temperature and void with a separate Euler step first and then held the total reactivity constant in all four stages, so only the linear kinetics were fourth-order; that is corrected here.

Assumptions and limitations The order 4 holds where the right-hand side is smooth. Two features of the model are not smooth: the reactivity discontinuities (handled exactly, §12) and the clip of the void target (a slope discontinuity where the target reaches $v_{min}$ or $v_{max}$ ). In the default parameter set the clip is never active, and the default transient converges with order about 4 (§16, tests C, C2). With a steeper void gain ( $K_{v}^{%} = 0.9$ , test C3) the target is clipped after the SCRAM and the observed accuracy stalls at about $10^{- 7}$ relative.

12. Treatment of Discontinuities

Intuition A step change in reactivity occurs at one instant. A numerical step must never straddle it; otherwise the event would be smeared or shifted by up to one time step.

Mathematical formulation The external reactivity is piecewise constant with breakpoints ${t_{1}, t_{2}, t_{SCRAM}}$ . The time grid consists of the major nodes $t_{n} = n Δ t$ (computed as a product, so no drift accumulates) together with the breakpoints and the end time. The driver advances from the current time $t$ to

t^{+} = min {(n + 1) Δ t, t_{next event}, t_{end}}

so a step is shortened whenever an event lies inside it. Within such a step $ρ_{ext}$ is constant and equal to the value of the segment being integrated, which is exact for all stage times. At each breakpoint two rows are stored with the same time: the state with the value before the event, and the state with the value after it (the flux, precursors, temperature and void are continuous; only $ρ_{ext}$ and $ρ$ jump). This draws the reactivity steps as exactly vertical lines and shows event times unambiguously in the table and CSV (event rows are marked).

Engineering interpretation If $t_{end}$ is not a multiple of $Δ t$ the last step is shortened to end exactly at $t_{end}$ . The reported time stamps are exactly the times at which states were computed; the notebook grid (linspace with int(t_max/dt) points) differed from the integration step. Coincident events are merged into one breakpoint and labelled accordingly. The event-splitting is tested with breakpoints deliberately placed off the time grid (§16, tests D1–D4).

13. Numerical Stability and Timestep

Intuition The model contains very different speeds: a prompt response in milliseconds and a precursor and feedback response in seconds to tens of seconds. An explicit method must resolve the fastest one, even when it is not the one of interest.

Mathematical formulation For $ρ < β$ the fastest mode of the linearised system is approximately the prompt mode with rate

Λ_{fast} \approx \frac{ρ - β}{T_{n}} < 0

An explicit method is stable if $z = h | Λ_{fast} |$ lies inside its stability region on the negative real axis: $z \leq 2$ for explicit Euler and $z \leq 2.785$ for RK4. For the default data the rate is $β / T_{n} = 163$ s⁻¹ at $ρ = 0$ and about 1660 s⁻¹ after the SCRAM (net $ρ_{ext} = - 5800$ pcm, minimum total reactivity −5996 pcm), so the stability limits are approximately $Δ t < 1.2$ ms for Euler and $Δ t < 1.7$ ms for RK4. The default time step is therefore $Δ t = 0.2$ ms ( $z \approx 0.33$ ); the value 1 ms of the notebooks is too large for the stronger SCRAM and shorter prompt lifetime of this parameter set (it would give $z \approx 1.7$ ). For $ρ > β$ the fast mode becomes unstable physically (growth rate $(ρ - β) / T_{n}$ ) and accuracy requires $h (ρ - β) / T_{n} ≪ 1$ .

Engineering interpretation The application evaluates $z$ from the external reactivity levels and issues a warning when $z > 0.5$ (Euler) or $z > 1$ (RK4) and a stability warning at the limits. It also warns if $Δ t > 0.2 min (τ_{T}, τ_{v})$ or is large compared with the shortest precursor life. The user’s time step is never changed silently, also not when the solver is switched. If the state becomes non-finite the run stops cleanly and reports the time and cause; a negative population or precursor, a void outside its limits and the prompt-critical threshold are reported when they occur. Test M (§16) demonstrates the advisory and the clean stop for RK4 at $Δ t = 0.05$ s.

Assumptions and limitations The estimate of $z$ uses external reactivity only (feedback is not known before the run). Total reactivity extremes are shown after the run; a stability warning does not appear if feedback alone causes a large step-size violation, in which case the non-finite check stops the run.

14. Reactor Period and the Inhour Relation

Intuition The reactor period is the time in which the flux changes by a factor $e$ . It is not simply the prompt lifetime divided by how far reactivity is above $β$ : for most reactivities the delayed neutrons determine it.

Mathematical formulation For constant reactivity and no source, the linear system (§4) has solutions $e^{ω t}$ . Substituting $ϕ, C_{i} \propto e^{ω t}$ gives the inhour equation [4, 5]

ρ = ω T_{n} + \sum_{i = 1}^{6} \frac{β_{i} ω}{ω + λ_{i}}

which has seven roots $ω_{j}$ for a given $ρ$ . The largest root $ω_{0}$ dominates after transients have decayed; the stable period is $τ = 1 / ω_{0}$ , positive for $ρ > 0$ (growth) and negative for $ρ < 0$ (decay). The right-hand side is strictly increasing for $ω > - λ_{min}$ , so $ω_{0}$ is unique on that interval and is found by bisection (bracket $[0, ρ / T_{n}]$ for $ρ > 0$ , $[- λ_{min}, 0]$ for $ρ < 0$ ) to machine precision. The prompt-only expression

τ_{prompt} \approx \frac{T_{n}}{ρ - β}

is the limit of the inhour equation for $ω ≫ λ_{i}$ , i.e. when the reactor is well above prompt critical. It is not a general reactor period and is not displayed by the application as one.

Stable period for the default six-group data and $T_{n} = 10^{- 4}$ s, computed with the application’s inhour solver, compared with the prompt-only expression.
ρ [pcm]	ρ/β	ω₀ [s⁻¹]	Stable period τ = 1/ω₀ [s]	Prompt-only T_n/(ρ−β) [s]
50	0.077	0.0073034	136.9	n/a (ρ < β)
100	0.154	0.018199	54.95	n/a (ρ < β)
200	0.308	0.057652	17.35	n/a (ρ < β)
300	0.461	0.14093	7.096	n/a (ρ < β)
500	0.769	0.82913	1.206	n/a (ρ < β)
650.2	1.000	7.5089	0.1332	diverges (ρ = β)
800	1.230	39.079	0.02559	0.0267
-100	-0.154	-0.0076929	-130	n/a (ρ < β)
-300	-0.461	-0.011156	-89.64	n/a (ρ < β)
-1000	-1.538	-0.012104	-82.62	n/a (ρ < β)
-6000	-9.228	-0.012355	-80.94	n/a (ρ < β)

Engineering interpretation The Reactor period (inhour) panel shows $τ$ for the external reactivity in force after each event; the cursor readout shows it for the instantaneous total reactivity. Both assume that $ρ$ is frozen; with feedback active the actual period changes continuously. The inhour root finder was checked against the dominant eigenvalue of the kinetics matrix computed by an independent power iteration (test F2, §16).

Prompt-critical threshold. When $ρ \geq β$ (650.2 pcm for the default data) the prompt neutrons alone make the chain reaction divergent; the flux then rises on the time scale of $T_{n}$ rather than of the precursors, and the delayed-neutron mechanism that normally makes the reactor controllable no longer limits the rate. The application shows “Prompt-critical threshold exceeded” as a message (and marks $β$ in the reactivity plot); it does not alter the calculation.

Assumptions and limitations The inhour relation assumes constant reactivity, the six-group data of the model and a fixed shape. Only the dominant root is used.

15. Model Parameters and Default Values

Provenance of the six-group data. The decay constants $λ_{i}$ in the notebooks (0.0124, 0.0305, 0.111, 0.301, 1.14, 3.01 s⁻¹) and the relative abundances $a_{i} = β_{i} / β$ implied by the notebook $β_{i}$ (0.033, 0.219, 0.196, 0.395, 0.115, 0.042) are identical to the tabulated six-group set for thermal fission of U-235, with half-lives 55.72, 22.72, 6.22, 2.30, 0.61 and 0.23 s [2]. The six-group representation originates from the experiments of Keepin, Wimett and Zeigler [1]. The notebook values $β_{i}$ equal $β a_{i}$ with $β = 0.0065$ [3] rounded to six decimals, which is why the sum is 0.006502 (650.2 pcm) rather than 0.0065. The original publication [1] was not re-read in full for this release, so the attribution of each individual number to it is not independently confirmed; the numbers were verified against [2] and [3]. Other evaluations of six-group data exist and differ [6].

The remaining defaults were retuned in this release to a 3000 MWth boiling-water-type reference core with 45 % average void. The table below lists the values, the second table the check against literature.

Symbol	Description	Default value	Unit	Source / basis	Editable
Delayed-neutron data (thermal fission, U-235) — external data
$β_{1} \dots β_{6}$	delayed fractions of groups 1–6	0.000215, 0.001424, 0.001274, 0.002568, 0.000748, 0.000273	–	$β a_{i}$ , $a_{i} = 0.033, 0.219, 0.196, 0.395, 0.115, 0.042$ , $β = 0.0065$ [2, 3]	yes
$β$	total delayed fraction $\sum β_{i}$	0.006502 (650.2 pcm)	–	calculated; standard value 0.0065 [3]	calculated
$λ_{1} \dots λ_{6}$	decay constants (half-lives 55.72, 22.72, 6.22, 2.30, 0.61, 0.23 s)	0.0124, 0.0305, 0.111, 0.301, 1.14, 3.01	s⁻¹	[2]	yes
Kinetics and normalisation
$T_{n}$	prompt neutron lifetime (used as generation time)	4.0·10⁻⁵	s	within the literature range 2·10⁻⁵ – 10⁻⁴ s [14]; consistent with §2	yes
$S$	external source	0	MWth/s	modelling assumption	yes
$P_{rated}$	rated thermal power (100 %); 1 a.u. = 1 MWth	3000	MWth	user-specified scaling (§2)	yes
$P (0)$	initial power; $ϕ_{0} = P_{rated} P (0) / 100 %$ is also the feedback reference power	100 % (3000 MWth)	% rated	user-specified	yes
$Φ_{rated}$	core-average thermal neutron flux at rated power	3.0·10¹³	n/(cm²·s)	2.2·10¹³ – 3.6·10¹³ for LWRs [9, 10]	yes
$v_{0}$	reference neutron speed for $Φ = n v_{0}$	2.2·10⁵ (2200 m/s)	cm/s	[13a]	yes (advanced)
External reactivity history
$ρ_{1}, t_{1}$	first step	+100 pcm at 10 s	pcm, s	[A1, A2]	yes
$ρ_{2}, t_{2}$	second step (added to $ρ_{1}$ )	+100 pcm at 20 s	pcm, s	[A1, A2]	yes
$ρ_{SCRAM}, t_{SCRAM}$	SCRAM reactivity insertion (added to $ρ_{1} + ρ_{2}$ ) and time	−6000 pcm at 40 s	pcm, s	rod worth 6000–9000 pcm, power defect [9a]; §6 (was −1000 pcm, hard-coded)	yes
$t_{end}$	simulation end	60	s	[A1, A2]	yes
Feedback (§7, §8)
$T_{f} (0), T_{ref}$	average fuel temperature at rated power	627	°C	900 K in a PWR study [15]; assumption	yes (T_ref advanced)
$α_{D}$	Doppler coefficient	−1.6 (−1.6·10⁻⁵ Δk/k/°C)	pcm/°C	BWR in operation [7]	yes
$K_{T}$	fuel target gain	3.4 (0.113 °C/MWth)	°C per % rated	derived: (627 − 286)/100 (§7)	yes
$τ_{T}$	fuel time constant	5.0	s	own estimate 4–5 s (§7); not verified	yes
$v (0), v_{ref}$	core-average void fraction at rated power	45 %	%	user-specified; about 40 % [8]	yes (v_ref advanced)
$α_{v}$	void coefficient	−70 pcm per % void (−0.7 Δk/k per unit)	pcm/% void	BWR in operation [7]; ~−30 [6a]	yes
$K_{v}$	void target gain	0.45 % per % rated (1.5·10⁻⁴ per MWth)	% void per % rated	secant 45 %/100 % (§8); assumption	yes
$τ_{v}$	void time constant	6.0	s	$τ_{T}$ +1 s (§8); estimate; not verified	yes
$v_{min}, v_{max}$	void limits (model validity range)	0 %, 70 %	%	modelling bound (§8; notebooks: 0–50 %)	yes (advanced)
Numerical parameters
Solver	RK4 (default) or explicit Euler	RK4	–	§10, §11	yes
$Δ t$	time step	0.0002	s	stability with $T_{n}$ and SCRAM value (§13); notebooks 0.01 / 0.001 s	yes
—	output / export interval (rounded to a multiple of $Δ t$ )	0.1	s	[A1, A2]	yes

Reasonableness check of the retuned constants

Constant	Notebooks (rescaled)	New default	Literature / reasoning	Assessment
Doppler $α_{D}$	−1 pcm/K	−1.6 pcm/°C	BWR in operation −1.6·10⁻⁵ Δk/k/°C [7]	consistent with the source
Void $α_{v}$	−200 pcm/%	−70 pcm/%	−70 pcm/% [7]; about −30 pcm/% for a typical BWR [6a]; PWR values do not apply	notebook value 3–7× too strong; corrected
Average void at rated power	0 %	45 %	about 40 % at design flow [8]; 45 % specified	consistent
Average fuel temperature	600 K (327 °C)	627 °C (900 K)	900 K used for a typical PWR [15]	consistent (assumption)
$K_{T}$	1 K per %	3.4 °C per %	(T_avg − T_sat)/100 %, with T_sat ≈ 286 °C (steam table, not verified)	derived; notebook value too small
$K_{v}$	0.2 % per %	0.45 % per %	secant of 0 → 45 % void over 0 → 100 % power	assumption
$τ_{T}$	5 s	5 s	own estimate 4–5 s from generic UO₂ properties; no checked source	plausible, not verified
$τ_{v}$	5 s	6 s	fuel lag + about 1 s transit; qualitative support in [7a]	plausible, not verified
Total power coefficient	−41 pcm/% (1 + 40)	−36.9 pcm/% (5.4 + 31.5)	model power defect 3.7 % Δk/k (PWR: about 2.5 % Δk/k quoted [9a]); no BWR-specific value checked	plausible order
$T_{n}$	1·10⁻⁴ s	4·10⁻⁵ s	about 2·10⁻⁵ s [14a] to about 10⁻⁴ s [14b]; consistent with $Φ_{rated}$ and power density (§2)	within range
SCRAM $ρ$	−1000 pcm (absolute level)	−6000 pcm (increment; net −5800)	total rod worth 6000–9000 pcm [9a]; must exceed the 3.7 % power defect; margin about 2.1 % Δk/k (net −5.8 %) vs 0.38 % licensing [9b]	retuned; now an increment
$Φ_{rated}$	—	3·10¹³ n/(cm²·s)	2.2·10¹³ [10] to 3.6·10¹³ [9]	consistent

If any value differs from the reference set, the header shows “Modified configuration” and the field is marked; exports (CSV, IC file) contain the complete parameter set.

16. Verification and Numerical Accuracy

Scope. This section documents software verification (the code solves the stated equations correctly) and numerical accuracy. It is not model validation: agreement between two solvers or with an analytical solution of the same equations says nothing about how well the simplified feedback model represents a real reactor. No validation against plant or experimental data has been performed, and the retuned constants of §15 are representative literature-based assumptions, not calibrated to any plant.

The suite below is executed in the application (footer → Verification) with the same code as the simulator. The table shows the results recorded for this release (v1.0.0) on the build machine; values that depend on floating-point arithmetic only are reproducible bit-for-bit in the same JavaScript engine.

Independent references. With both feedbacks disabled ( $α_{D} = α_{v} = 0$ ) and piecewise-constant $ρ_{ext}$ , the system is linear on every segment and its exact solution is a product of matrix exponentials $𝐲 (t_{k + 1}) = e^{𝐌_{k} Δ_{k}} 𝐲 (t_{k})$ , computed by scaling-and-squaring with an independent implementation. With feedback active, the reference is RK4 at $Δ t = 0.05$ ms.

ID	Test	Acceptance criterion	Result	Measured
A-RK4	Zero-reactivity equilibrium (RK4, Δt = 0.01 s, 60 s)	max \|Δφ/φ0\|, \|ΔCi/Ci0\| < 1e-9; T, v unchanged	PASS	φ 0, C 0, T 0 K, v 0
A-Euler	Zero-reactivity equilibrium (Euler, Δt = 0.01 s, 60 s)	max \|Δφ/φ0\|, \|ΔCi/Ci0\| < 1e-9; T, v unchanged	PASS	φ 0, C 0, T 0 K, v 0
B	Euler → RK4 with decreasing Δt (default transient, feedback on)	error decreases monotonically; observed order ≈ 1	PASS	Δt 0.0008: 5.24e-5; Δt 0.0004: 1.97e-5; Δt 0.0002: 9.85e-6; Δt 0.0001: 4.92e-6; order 1.41, 1.00, 1.00
C	RK4 refinement Δt, Δt/2, Δt/4 (default transient, feedback on)	error decreases with each halving; observed order > 3.0 for the first and between 3.5 and 4.6 for the second halving (first pair still slightly pre-asymptotic: Δt·\|ρ−β\|/Tn ≈ 1.3)	PASS	Δt 0.0008: 8.13e-11; Δt 0.0004: 8.93e-12; Δt 0.0002: 5.42e-13; observed order 3.19, 4.04
C2	RK4 refinement, smooth sub-case (steps at 10 s, 20 s; no SCRAM)	observed order between 3.5 and 4.8 for both halvings	PASS	Δt 0.008: 7.36e-11; Δt 0.004: 2.97e-12; Δt 0.002: 1.47e-13; order 4.63, 4.33
C3	RK4 refinement with the void-target clip active after SCRAM (K_v = 0.9 % void per % rated; informational)	error decreases from the coarsest step and is < 1e-6 at Δt/4; the order is not required to be 4 and the error stalls near 1e-7 (slope discontinuity of the clip)	PASS	Δt 0.0008: 6.33e-6; Δt 0.0004: 9.85e-8; Δt 0.0002: 1.11e-7; observed order 6.01, -0.17
C4	RK4 refinement with void dynamics removed (K_v = 0), full transient incl. SCRAM	observed order between 3.5 and 4.6 for both halvings	PASS	Δt 0.0008: 3.65e-11; Δt 0.0004: 3.10e-12; Δt 0.0002: 1.48e-13; order 3.56, 4.39
D1	Event times off the Δt grid (t1 = 2.0004, t2 = 5.0007, tSCRAM = 9.0001 s; Δt = 0.003 s, T_n = 1 ms)	each breakpoint is an exact time node with left/right ρext; no step straddles a breakpoint	PASS	exact nodes true, left/right rows true, no straddling true, ρext = 0.0/100.0/200.0/-5800.0 pcm
D2	Coincident events (t1 = t2 = tSCRAM; t1 = t2; t1 = 0)	deterministic, stateless profile: all three levels accumulate (ρ1 + ρ2 + ρSCRAM), t = 0 event applies at once	PASS	all-coincident true, t1 = t2 → 200 pcm true, t1 = 0 true
D3	Non-multiple Δt: tEnd = 1.0 s, Δt = 0.3 s (final shortened step)	last node equals tEnd exactly	PASS	t_last = 1
D4	Time grid t_n = n·Δt (no linspace mismatch)	max \|t − nΔt\| < 1e-12 s	PASS	max deviation 0 s
H1	Feedback disabled (αD = αV = 0): RK4 vs exact matrix-exponential solution, default transient	max relative φ error < 1e-8	PASS	RK4 Δt = 0.001: 1.21e-11
H2	Feedback disabled: Explicit Euler vs exact solution	first-order error, < 2e-2 at Δt = 0.001	PASS	Euler Δt = 0.001: 1.17e-4
E	Positive sub-prompt-critical step (+300 pcm, feedback off)	matches exact solution (<1e-8); late growth rate within 0.5 % of inhour root ω0	PASS	vs exact 2.68e-13; ω_sim = 0.14105 s⁻¹, ω0 = 0.14093 s⁻¹ (τ = 7.10 s)
F	Negative step (−300 pcm, feedback off)	matches exact solution (<1e-8); monotonic decay after the prompt drop; decay rate bounded by the inhour root (ω_sim ≤ ω0, since sub-dominant modes decay faster)	PASS	vs exact 8.45e-12; prompt-drop ratio φ(10.02)/φ(9.99) = 0.6866; ω_sim(80–100 s) = -0.01577 s⁻¹, ω0 = -0.01116 s⁻¹
F2	Inhour root ω0 vs dominant eigenvalue of the point-kinetics matrix (independent power iteration on exp(MΔτ))	relative difference < 1e-8 for ρ = −1000 … +800 pcm (both signs, beyond prompt critical)	PASS	max rel. diff 2.59e-11; stable periods: -1000 pcm: -82.618 s, -300 pcm: -89.635 s, -50 pcm: -209.882 s, 50 pcm: 136.922 s, 100 pcm: 54.950 s, 200 pcm: 17.345 s, 300 pcm: 7.096 s, 500 pcm: 1.206 s, 650.2 pcm: 0.133 s, 800 pcm: 0.026 s
G	SCRAM transition (ρSCRAM added to ρ1 + ρ2 like the other steps)	ρext = ρ1 + ρ2 + ρSCRAM for t ≥ tSCRAM (default 200 − 6000 = −5800 pcm; −1500 pcm case); a 0 pcm SCRAM leaves ρext unchanged; prompt drop reproduced by exact solution (<1e-8)	PASS	ρext(40.5 s) = -5800.000 pcm; vs exact 4.01e-12; φ(40.05)/φ(39.99) = 0.227
G2	Zero-worth SCRAM (ρSCRAM = 0 pcm) changes nothing	ρext stays at ρ1 + ρ2 = 200 pcm for t ≥ tSCRAM and the flux shows no discontinuity at tSCRAM	PASS	ρext(t ≥ 41 s) = 200 pcm: true; relative flux change over 39.9 → 40.0 s: 2.27e-6
H3	Feedback disabled: φ independent of KT, Kv, τT, τv	bit-identical flux histories	PASS	identical (4011 nodes)
I	Feedback isolation: none / Doppler only / void only / both	ρD = 0 when αD = 0, ρv = 0 when αv = 0; ρ = ρext + ρD + ρv; φmax(none) > φmax(single) > φmax(both)	PASS	φmax: 22497.3 / 4162.8 / 3645.0 / 3626.3 MWth; sum residual 1.82e-12 pcm
J	IC round trip (non-default values in every group)	parsed configuration identical to source; identical calculated results	PASS	config identical true, results identical true
K	Invalid IC files rejected (11 malformed variants)	parseIC returns errors and no configuration for every variant (state untouched)	PASS	invalid JSON: rejected; wrong format: rejected; unsupported schema: rejected; wrong application: rejected; missing array: rejected; short array: rejected; non-numeric value: rejected; out-of-range value: rejected; event order violated: rejected; unknown solver: rejected; missing group: rejected
L-RK4	Validity and determinism, default case (RK4, default Δt = 0.0002 s)	all finite; φ, Ci ≥ 0; 0 ≤ v ≤ 0.5; repeated run bit-identical	PASS	finite true, min C 3.13e+2, min φ 136.697, v in range true, checksum 870aac55 (repeatable)
L-Euler	Validity and determinism, default case (Euler, default Δt = 0.0002 s)	all finite; φ, Ci ≥ 0; 0 ≤ v ≤ 0.5; repeated run bit-identical	PASS	finite true, min C 3.13e+2, min φ 136.697, v in range true, checksum fd090484 (repeatable)
R-RK4	Default-case regression values (RK4, default Δt = 0.0002 s)	Pmax, t(Pmax), P(60 s), Tmax, T(60 s), vmax, v(60 s), ρmax, ρmin reproduce stored release values (rel. dev. < 1e-9)	PASS	φmax 3626.254 MWth at 20.0368 s; φ(60) 136.6971; Tmax 646.059 °C; vmax 47.3884 %; ρ -5996.15 … 102.83 pcm; max dev 0
R-Euler	Default-case regression values (Euler, default Δt = 0.0002 s)	Pmax, t(Pmax), P(60 s), Tmax, T(60 s), vmax, v(60 s), ρmax, ρmin reproduce stored release values (rel. dev. < 1e-9)	PASS	φmax 3626.512 MWth at 20.0364 s; φ(60) 136.6968; Tmax 646.059 °C; vmax 47.3884 %; ρ -5996.15 … 102.83 pcm; max dev 0
M	Unstable time step (RK4, Δt = 0.05 s) is advised and stops cleanly	advisory issued before the run; run stops with a non-finite-state message; no NaN in stored data	PASS	advisory true, stopped at t = 10.15 s, stored data finite true

Default-case regression values (feedback active, $Δ t = 0.0002$ s)

Quantity (default configuration)	RK4	Explicit Euler
Maximum flux / power [MWth]	3626.2545	3626.5118
Time of maximum flux [s]	20.037	20.036
Flux / power at 60 s [MWth]	136.6971	136.6968
Maximum average fuel temperature [°C]	646.0590	646.0591
Average fuel temperature at 60 s [°C]	310.8765	310.8759
Maximum void fraction [%]	47.38843	47.38843
Void fraction at 60 s [%]	3.98614	3.98605
Maximum total reactivity [pcm]	102.830	102.830
Minimum total reactivity [pcm]	-5996.154	-5996.154
Integration steps	300000	300000

Interpretation. (1) Without feedback RK4 reproduces the exact solution to about $10^{- 11}$ relative (H1), and Euler shows the expected first-order error of about $10^{- 4}$ at $Δ t = 1$ ms (H2). (2) Euler converges toward RK4 with observed order 1.0 as the step is reduced (B). (3) RK4 with feedback converges with observed order about 3.2 for the first halving (still slightly pre-asymptotic, $z = h | ρ - β | / T_{n} \approx 1.3$ ) and about 4.0 for the second (C); the smooth sub-case gives 4.3–4.6 (C2) and the run without void dynamics 3.6–4.4 (C4). The retuning changed this picture relative to the notebook parameters: with $K_{v}^{%} = 0.45$ the void target is never clipped in the default case, so the default transient is smooth; with a steeper gain (0.9 % per %, test C3) the clip is active after the SCRAM and the accuracy stalls at about $10^{- 7}$ relative. (4) Event handling is exact for breakpoints off the grid (D1–D4). (5) The default $Δ t = 0.2$ ms is in the asymptotic range for both solvers; the difference between the two solvers on the default case is 0.007 % in maximum power (3626.51 vs 3626.25 MWth), as expected for a first-order method.

Not covered by this suite. Behaviour for very stiff or prompt-supercritical cases beyond the stability warnings; extremely small time steps (accumulation of round-off); browsers other than the one listed in §17.4.

17. Scope and Limitations

17.1 What the model represents

Point-reactor neutron kinetics with six delayed-neutron precursor groups, an optional external source, and a user-defined external reactivity history (two steps and a SCRAM).
Fuel-temperature (Doppler) and moderator-void reactivity feedback via first-order lags toward flux-dependent targets.
Numerical solution by explicit Euler or classical RK4 of one common coupled state model, with exact treatment of the reactivity breakpoints.

17.2 What the model does not represent

Spatial neutron-flux distribution and local power peaking; xenon and spatial oscillations; detailed control-rod worth distributions; thermal-hydraulic channel dynamics; fuel heat conduction and heat capacity; coolant transport delay; detailed boiling dynamics and flow; core pressure dynamics; changes of the neutron spectrum; multidimensional feedback; burn-up; plant protection logic and reactor-specific safety systems. The thermal-hydraulic part of this application is not a mechanistic model; $K_{T}, K_{v}, τ_{T}, τ_{v}, α_{D}, α_{v}$ are lumped, illustrative parameters. The model is suitable for point-kinetics education, engineering understanding, transient exploration, solver comparison and qualitative feedback studies. It must not be used as evidence of plant behaviour, and no claim of plant validation or safety-grade applicability is made.

17.3 Source-notebook audit and resolutions

Finding in the source notebooks	Resolution in this application
Text states $d T_{f} / d t = k_{T} (ϕ - ϕ_{0})$ and $d v / d t = k_{v} (ϕ - ϕ_{0})$ ; code implements first-order lags to a flux-dependent target	Executable lag model implemented and documented (§7, §8).
Euler notebook: reactivity of the new step is computed but the flux derivative uses the previous element (one-step feedback delay)	No delay: Euler and RK4 call the same derivative on the same state (§10).
RK notebook: temperature and void advanced separately by Euler, total reactivity frozen in all RK stages; text still describes Euler	RK4 integrates the full coupled 9-component state; reactivity is re-evaluated in every stage (§11).
Time grid `linspace(0,t_max,int(t_max/dt))` does not have spacing dt	Grid built from $t_{n} = n Δ t$ ; final step shortened if needed (§12).
No event handling: a step may jump over a reactivity change	Steps are split at every breakpoint (§12).
First-step reactivity kept in a variable set only in one branch (undefined or stale if $t_{1} = t_{2}$ )	Stateless piecewise function (§6).
Reactor-specific constants (fuel 600 K, void 0 %, α_v = −200 pcm/%, K_T, K_v, T_n, τ) without literature basis	Retuned to a 3000 MWth BWR-type core with 45 % void, temperature in °C, constants checked against literature where available; unverified ones flagged (§15).
Flux in arbitrary units (φ₀ = 1000 a.u.), feedback gains per a.u.	Normalised to fission power (1 a.u. = 1 MWth, 3000 MWth at 100 %); gains defined per percent of rated power so that the dynamics are unchanged (§2, §7, §8).
SCRAM reactivity hard-coded to −1000 pcm, replacing the accumulated external reactivity	Editable input, handled as an increment like the two steps (0 pcm = no change); default retuned to −6000 pcm (§6).
Void clamped to [0, 0.5] for the target but [0, 1] for the state; manual states 0–0.5	One range (editable limits, default 0–70 %) applied consistently; excursions are reported, not hidden (§8).
Void reactivity without reference subtraction; fuel target based on the initial temperature; $T_{ref}$ hard-coded to 600 K	$ρ_{v} = α_{v} (v - v_{ref})$ , $ρ_{D} = α_{D} (T_{f} - T_{ref})$ with explicit reference parameters (identical results at the defaults).
Different default time steps (Euler 0.001 s, RK 0.01 s) and slider ranges (e.g. reactivity ±1000 vs ±500 pcm, Doppler range)	One default $Δ t = 0.0002$ s (§13), wider validated ranges for both solvers; RK4 at 0.01 s is unstable for the retuned parameter set.
“Effective initial reactor period” $T_{n} / (ρ_{1} - β)$ , which is negative and not a reactor period for $ρ_{1} < β$	Replaced by the inhour stable period (§14).
Excel export with charts (openpyxl)	CSV with metadata plus IC JSON; no external library or network dependency.

17.4 Software and environment

Single-file HTML, vanilla JavaScript, no external resources; calculation runs in the browser, in chunks so that the interface remains responsive. Integration is limited to 10⁷ steps. All results are deterministic for identical inputs. The release was exercised automatically in headless Chromium (Playwright 1.56): default run, both solvers, hover/zoom/keyboard cursor, IC save/load round trip, rejection of malformed IC files, CSV and PNG export, dialogs, the separate theory window, and the verification suite; Firefox and Edge were not part of the automated tests. Provided for engineering guidance, analysis and training; results requiring safety, regulatory, contractual or operational significance shall be independently verified.

18. References

External references (consulted)

[1] G. R. Keepin, T. F. Wimett, R. K. Zeigler, “Delayed neutrons from fissionable isotopes of uranium, plutonium, and thorium”, Physical Review 107, 1044–1049 (1957), doi:10.1103/PhysRev.107.1044. (Bibliographic record and abstract checked; six-group representation.)
[2] Los Alamos National Laboratory, LA-UR-13-20680, lecture slides on passive neutron analysis, slide “Delayed Neutron Data (U-235 for Thermal Incident Neutron)”, OSTI record 1062155. (Checked: half-lives, decay constants, relative yields.)
[3] U.S. Department of Energy, DOE Fundamentals Handbook: Nuclear Physics and Reactor Theory, DOE-HDBK-1019/1-93 and DOE-HDBK-1019/2-93, January 1993. (Checked: six precursor groups for thermal fission of U-235, delayed-neutron fraction of U-235 = 0.0065, definition of reactivity.)
[4] MIT OpenCourseWare 22.05 Reactor Physics (Fall 2006), lecture notes “Part Twenty-Eight: Dynamic Period and Inhour Equations”. (Checked: inhour equation, roots, limitations of constant-reactivity solution.)
[5] Nuclear Engineering course notes EP4D3, “Reactor Physics: Point Kinetics”, hosted at nuceng.ca. (Checked: seven roots of the inhour equation, stable period, prompt-critical limit.)
[6] Measurements of effective delayed neutron fraction in a fast neutron reactor using the perturbation method, arXiv:1510.06276. (Checked only for its comparison table showing that six-group parameters differ between Keepin and evaluated libraries.)

Additional references for the retuned constants (consulted at abstract/excerpt level via search)

[6a] Article record “Void reactivity coefficient benchmark results for a 10 × 10 BWR assembly in the full 0–100 % void fraction range” (ScienceDirect, S0306454909000693); the excerpt seen quotes about 30 pcm/% void for a typical BWR (attributed there to Jatuff et al., 2009). Abstract-level only; not read in full.
[7] US patent 5,940,461 (also US 5,812,621), “Reactor core for a light water cooled reactor, fuel assembly and control rod”, background section: Doppler coefficient 1.6·10⁻⁵ Δk/k/°C and void coefficient −7.0·10⁻⁴ Δk/k/% void “of BWR now under operation”. [7a] US patent 5,524,128, “Boiling water reactor stability control”: void response follows the flux after a delay set by the fuel thermal time constant (qualitative).
[8] US patent 4,292,128, “Boiling water type nuclear reactor”: mean void fraction across the core of about 40 % at design flow, rising to 56 % at reduced flow.
[9] MIT Nuclear Reactor Laboratory, “In-Core Experiments” (nrl.mit.edu): in-core irradiation facilities with flux levels comparable to a commercial LWR, thermal flux up to 3.6·10¹³ n/(cm²·s), fast flux (>0.1 MeV) up to 1.2·10¹⁴. [9a] nuclear-power.com, “Reactor Shutdown | Condition & SCRAM”: total control-rod worth design-specific, e.g. 6000–9000 pcm; PWR power defect about 2500 pcm. [9b] US patent 4,629,599, “Burnable absorber arrangement for fuel bundle”: BWR licensing shutdown margin 0.38 % with one rod stuck out, 1 % design basis. [9c] U.S. NRC, Westinghouse Technology Systems Manual, Section 2.1 (PWR void coefficient −30 to −250 pcm/% void; PWR only).
[10] ScienceDirect Topics, “Pressurized Water Reactor”: average neutron flux 2.16·10¹⁷ (thermal) and 3·10¹⁸ (fast) n/(m²·s) for a typical PWR.
[11] ScienceDirect Topics, “Prompt Neutron”: precursor concentration $c_{i}$ in nuclei/cm³.
[12] arXiv:1603.00195, “A nonstiff solution for the stochastic neutron point kinetics equations”: point kinetics in neutron density $n$ and precursor density $C_{i}$ , $C_{i} (0) = β_{i} n_{0} / (Λ λ_{i})$ .
[13] arXiv:1306.1119, “Analysis of Reactivity Induced Accident for Control Rods Ejection with Loss of Cooling”: point kinetics with all states in MW. [13a] EOLSS chapter “Nuclear Reactor Kinetics”: neutron speed 2200 m/s at 20 °C; lifetimes in LWR below 1 ms.
[14] Prompt neutron lifetime: [14a] nuclear-power.com, “Prompt Neutron Lifetime”: about 2·10⁻⁵ s in the worked example; [14b] IAEA TCS-70, “Introduction to Water Cooled Reactor Theory with the Micro…”: neutron lifetime timespan of about 10⁻⁴ s.
[15] arXiv:1810.06431, “Neutronic Study of UO2-BeO Fuel with Various Claddings”: an average fuel temperature of 900 K is used for a typical PWR in the simulations.

Antranova model basis

[A1] Antranova Engineering AB, Point Kinetics Simulator Including TH-Feedbacks, Solver: Euler, notebook Antranova_PointKineticsFB_EU, version 1.0.0.
[A2] Antranova Engineering AB, Point Kinetics Simulator Including TH-Feedbacks, Solver: Runge-Kutta (RK4), notebook Antranova_PointKineticsFB_RK, version 1.0.0.

Point Kinetics Simulator

Theory & Method

Theory & Method

1. Physical Background

2. Point-Reactor Approximation

3. Prompt and Delayed Neutrons

4. Six-Group Point Kinetics Equations

5. Initial Equilibrium Conditions

6. External Reactivity Model

7. Fuel-Temperature Feedback

8. Moderator-Void Feedback

9. Coupled State Model

10. Explicit Euler Method

11. Runge–Kutta Fourth-Order Method

12. Treatment of Discontinuities

13. Numerical Stability and Timestep

14. Reactor Period and the Inhour Relation

15. Model Parameters and Default Values

Reasonableness check of the retuned constants

16. Verification and Numerical Accuracy

Default-case regression values (feedback active, $Δ t = 0.0002$ s)

17. Scope and Limitations

17.1 What the model represents

17.2 What the model does not represent

17.3 Source-notebook audit and resolutions

17.4 Software and environment

18. References

External references (consulted)

Additional references for the retuned constants (consulted at abstract/excerpt level via search)

Further reading (bibliographic entries given for orientation; not consulted in preparing this release)

Antranova model basis

About & Licence

Point Kinetics Simulator

Verification suite