Compare commits

..

3 Commits

Author SHA1 Message Date
Codex Agent
03595b0d12 Queue dashboard polish items 2025-11-24 00:06:57 +01:00
Codex Agent
8ece3efb04 Improve dashboard spacing and steam context 2025-11-24 00:06:13 +01:00
Codex Agent
f6ff6fc618 Add delayed kinetics and steam-drum balance 2025-11-24 00:06:08 +01:00
10 changed files with 185 additions and 21 deletions

View File

@@ -1,9 +1,9 @@
## C.O.R.E. feature set
- **Core physics**: point-kinetics with delayed neutrons, temperature feedback, fuel burnup penalty, xenon/iodine buildup with decay and burn-out, and rod-bank worth curves.
- **Core physics**: point-kinetics with per-bank delayed neutron precursors, temperature feedback, fuel burnup penalty, xenon/iodine buildup with decay and burn-out, and rod-bank worth curves.
- **Rod control**: three rod banks with weighted worth; auto controller chases 3 GW setpoint; manual mode with staged bank motion and SCRAM; state persists across runs.
- **Coolant & hydraulics**: primary/secondary pumps with head/flow curves, power draw scaling, wear tracking; pressure floors tied to saturation; auxiliary power model with generator auto-start.
- **Heat transfer**: steam-generator UA·ΔT_lm model with a pinch cap to keep the primary outlet hotter than the secondary, coolant heating uses total fission power with fuel heating decoupled from exchanger draw, and the secondary thermal solver includes passive cool-down when flow is low.
- **Heat transfer**: steam-generator UA·ΔT_lm model with a pinch cap to keep the primary outlet hotter than the secondary, coolant heating uses total fission power with fuel heating decoupled from exchanger draw, and the secondary thermal solver includes passive cool-down plus steam-drum mass/energy balance with latent heat.
- **Pressurizer & inventory**: primary pressurizer trims toward 7 MPa with level tracking, loop inventories/levels steer flow availability, secondary steam boil-off draws down level with auto makeup, and pumps reduce flow/status to `CAV` when NPSH is insufficient.
- **Steam cycle**: three turbines with spool dynamics, throttle mapping, condenser back-pressure penalty, load dispatch to consumer, steam quality gating for output, generator states with batteries/spool.
- **Protections & failures**: health monitor degrading components under stress, automatic SCRAM on core or heat-sink loss, relief valves per loop, maintenance actions to restore integrity.

12
TODO.md
View File

@@ -1,7 +1,9 @@
## Future realism upgrades
- Steam generator UA·ΔT_lm heat exchange (done) and pump head/flow curves (done); keep validating temps under nominal load.
- Rod banks with worth curves and xenon/samarium buildup (done); extend with delayed-group kinetics per bank.
- Add pressurizer behavior, primary/secondary inventory and level effects, and pump NPSH/cavitation checks.
- Model feedwater/steam-drum mass-energy balance, turbine throttle/efficiency maps, and condenser back-pressure.
- Introduce CHF/DNB margin, clad/fuel split temps, and SCRAM matrix for subcooling loss or SG level/pressure trips.
- [x] Steam generator UA·ΔT_lm heat exchange and pump head/flow curves; keep validating temps under nominal load.
- [x] Rod banks with worth curves, xenon/samarium buildup, and delayed-group kinetics per bank.
- [x] Pressurizer behavior, primary/secondary inventory and level effects, and pump NPSH/cavitation checks.
- [x] Model feedwater/steam-drum mass-energy balance, turbine throttle/efficiency maps, and condenser back-pressure.
- [ ] Introduce CHF/DNB margin, clad/fuel split temps, and SCRAM matrix for subcooling loss or SG level/pressure trips.
- [ ] Flesh out condenser behavior: vacuum pump limits, cooling water temperature coupling, and dynamic back-pressure with fouling.
- [ ] Dashboard polish: compact turbine/generator rows, color critical warnings (SCRAM/heat-sink), and reduce repeated log noise.

View File

@@ -8,6 +8,7 @@ NEUTRON_LIFETIME = 0.1 # seconds, prompt neutron lifetime surrogate
FUEL_ENERGY_DENSITY = 200.0 * MEGAWATT # J/kg released as heat
COOLANT_HEAT_CAPACITY = 4_200.0 # J/(kg*K) for water/steam
COOLANT_DENSITY = 700.0 # kg/m^3 averaged between phases
STEAM_LATENT_HEAT = 2_200_000.0 # J/kg approximate latent heat of vaporization
CORE_MELTDOWN_TEMPERATURE = 2_873.0 # K (approx 2600C) threshold for irreversible meltdown
MAX_CORE_TEMPERATURE = CORE_MELTDOWN_TEMPERATURE # Allow simulation to approach meltdown temperature
MAX_PRESSURE = 15.0 # MPa typical PWR primary loop limit

View File

@@ -308,18 +308,19 @@ class ReactorDashboard:
log_rows = max(2, min(len(self.log_buffer) + 2, 8))
status_height = min(height - 1, max(min_status, min_status + log_rows))
data_height = max(1, height - status_height)
gap = 2
right_width = max(28, width // 3)
left_width = width - right_width
left_width = width - right_width - gap
if left_width < 50:
left_width = min(50, width - 18)
right_width = width - left_width
left_width = min(50, width - (18 + gap))
right_width = width - left_width - gap
data_height = max(1, data_height)
left_width = max(1, left_width)
right_width = max(1, right_width)
data_win = stdscr.derwin(data_height, left_width, 0, 0)
help_win = stdscr.derwin(data_height, right_width, 0, left_width)
help_win = stdscr.derwin(data_height, right_width, 0, left_width + gap)
status_win = stdscr.derwin(status_height, width, data_height, 0)
self._draw_data_panel(data_win, state)
@@ -341,6 +342,8 @@ class ReactorDashboard:
right_win = win.derwin(inner_height, right_width, 1, 1 + left_width)
for row in range(1, height - 1):
win.addch(row, 1 + left_width, curses.ACS_VLINE)
if left_width + 1 < width - 1:
win.addch(row, 1 + left_width + 1, curses.ACS_VLINE)
left_win.erase()
right_win.erase()
@@ -404,7 +407,7 @@ class ReactorDashboard:
f"{state.secondary_loop.mass_flow_rate:7.0f}/{self.reactor.secondary_pump.nominal_flow * len(self.reactor.secondary_pump_units):.0f} kg/s",
),
("Level", f"{state.secondary_loop.level*100:6.1f}%"),
("Inlet Temp", f"{state.secondary_loop.temperature_in:7.1f} K (Target {constants.PRIMARY_OUTLET_TARGET_K:4.0f})"),
("Inlet Temp", f"{state.secondary_loop.temperature_in:7.1f} K (Target {constants.SECONDARY_OUTLET_TARGET_K:4.0f})"),
("Outlet Temp", f"{state.secondary_loop.temperature_out:7.1f} K (Target {constants.SECONDARY_OUTLET_TARGET_K:4.0f})"),
("Pressure", f"{state.secondary_loop.pressure:5.2f}/{constants.MAX_PRESSURE:4.1f} MPa"),
("Steam Quality", f"{state.secondary_loop.steam_quality:5.2f}/1.00"),
@@ -440,6 +443,10 @@ class ReactorDashboard:
("Load", f"{self._total_load_supplied(state):7.1f}/{self._total_load_demand(state):7.1f} MW"),
("Consumer", f"{consumer_status}"),
("Demand", f"{consumer_demand:7.1f} MW"),
(
"Steam",
f"P={state.secondary_loop.pressure:4.2f} MPa q={state.secondary_loop.steam_quality:4.2f} mdot={state.secondary_loop.mass_flow_rate:6.0f} kg/s",
),
],
)
right_y = self._draw_section(right_win, right_y, "Generators", self._generator_lines(state))
@@ -745,7 +752,16 @@ class _DashboardLogHandler(logging.Handler):
def __init__(self, buffer: deque[str]) -> None:
super().__init__()
self.buffer = buffer
self._last_msg: str | None = None
self._repeat_count: int = 0
def emit(self, record: logging.LogRecord) -> None:
msg = self.format(record)
if msg == self._last_msg:
self._repeat_count += 1
if self._repeat_count > 3:
return
else:
self._last_msg = msg
self._repeat_count = 0
self.buffer.append(msg)

View File

@@ -26,7 +26,7 @@ def xenon_poisoning(flux: float) -> float:
@dataclass
class NeutronDynamics:
beta_effective: float = 0.0065
delayed_neutron_fraction: float = 0.0008
delayed_decay_const: float = 0.08 # 1/s effective precursor decay
external_source_coupling: float = 1e-6
shutdown_bias: float = -0.014
iodine_yield: float = 1e-6 # inventory units per MW*s
@@ -55,12 +55,18 @@ class NeutronDynamics:
return rho
def flux_derivative(
self, state: CoreState, rho: float, external_source_rate: float = 0.0, baseline_source: float = 1e5
self,
state: CoreState,
rho: float,
delayed_source: float,
external_source_rate: float = 0.0,
baseline_source: float = 1e5,
) -> float:
generation_time = constants.NEUTRON_LIFETIME
beta = self.beta_effective
source_term = self.external_source_coupling * external_source_rate
return ((rho - beta) / generation_time) * state.neutron_flux + baseline_source + source_term
prompt = ((rho - beta) / generation_time) * state.neutron_flux
return prompt + delayed_source + baseline_source + source_term
def step(
self,
@@ -77,9 +83,22 @@ class NeutronDynamics:
rho = min(rho, -0.04)
baseline = 0.0 if shutdown else 1e5
source = 0.0 if shutdown else external_source_rate
d_flux = self.flux_derivative(state, rho, source, baseline_source=baseline)
rod_positions = rod_banks if rod_banks else [control_fraction] * len(constants.CONTROL_ROD_BANK_WEIGHTS)
self._ensure_precursors(state, len(rod_positions))
bank_factors = self._bank_factors(rod_positions)
bank_betas = self._bank_betas(len(bank_factors))
delayed_source = self._delayed_source(state, bank_factors)
d_flux = self.flux_derivative(
state,
rho,
delayed_source,
external_source_rate=source,
baseline_source=baseline,
)
state.neutron_flux = max(0.0, state.neutron_flux + d_flux * dt)
state.reactivity_margin = rho
self._update_precursors(state, bank_factors, bank_betas, dt)
LOGGER.debug(
"Neutronics: rho=%.5f, flux=%.2e n/cm2/s, d_flux=%.2e",
rho,
@@ -102,3 +121,38 @@ class NeutronDynamics:
def _xenon_penalty(self, state: CoreState) -> float:
return min(0.05, state.xenon_inventory * self.xenon_reactivity_coeff)
def _bank_betas(self, bank_count: int) -> list[float]:
weights = list(constants.CONTROL_ROD_BANK_WEIGHTS)
if bank_count != len(weights):
weights = [1.0 for _ in range(bank_count)]
total = sum(weights) if weights else 1.0
return [self.beta_effective * (w / total) for w in weights]
def _bank_factors(self, positions: list[float]) -> list[float]:
factors: list[float] = []
for pos in positions:
insertion = clamp(pos, 0.0, 0.95)
factors.append(max(0.0, 1.0 - insertion / 0.95))
return factors
def _ensure_precursors(self, state: CoreState, bank_count: int) -> None:
if not state.delayed_precursors or len(state.delayed_precursors) != bank_count:
state.delayed_precursors = [0.0 for _ in range(bank_count)]
def _delayed_source(self, state: CoreState, bank_factors: list[float]) -> float:
decay = self.delayed_decay_const
return sum(decay * precursor * factor for precursor, factor in zip(state.delayed_precursors, bank_factors))
def _update_precursors(
self, state: CoreState, bank_factors: list[float], bank_betas: list[float], dt: float
) -> None:
generation_time = constants.NEUTRON_LIFETIME
decay = self.delayed_decay_const
new_pools: list[float] = []
for precursor, factor, beta in zip(state.delayed_precursors, bank_factors, bank_betas):
production = (beta / generation_time) * state.neutron_flux * factor
loss = decay * precursor
updated = max(0.0, precursor + (production - loss) * dt)
new_pools.append(updated)
state.delayed_precursors = new_pools

View File

@@ -95,6 +95,7 @@ class Reactor:
reactivity_margin=-0.02,
power_output_mw=0.1,
burnup=0.0,
delayed_precursors=[0.0 for _ in constants.CONTROL_ROD_BANK_WEIGHTS],
fission_product_inventory={},
emitted_particles={},
)
@@ -387,7 +388,7 @@ class Reactor:
else:
transferred = heat_transfer(state.primary_loop, state.secondary_loop, total_power)
self.thermal.step_core(state.core, state.primary_loop, total_power, dt)
self.thermal.step_secondary(state.secondary_loop, transferred)
self.thermal.step_secondary(state.secondary_loop, transferred, dt)
self._apply_secondary_boiloff(state, dt)
self._update_loop_inventory(
state.secondary_loop, constants.SECONDARY_LOOP_VOLUME_M3, constants.SECONDARY_INVENTORY_TARGET, dt

View File

@@ -20,6 +20,7 @@ class CoreState:
burnup: float # fraction of fuel consumed
xenon_inventory: float = 0.0
iodine_inventory: float = 0.0
delayed_precursors: list[float] = field(default_factory=list)
fission_product_inventory: dict[str, float] = field(default_factory=dict)
emitted_particles: dict[str, float] = field(default_factory=dict)

View File

@@ -60,6 +60,19 @@ def saturation_pressure(temp_k: float) -> float:
return min(constants.MAX_PRESSURE, max(0.01, psat_mpa))
def saturation_temperature(pressure_mpa: float) -> float:
"""Approximate saturation temperature (K) for water at the given pressure."""
target = max(0.01, min(constants.MAX_PRESSURE, pressure_mpa))
low, high = 273.15, 900.0
for _ in range(40):
mid = 0.5 * (low + high)
if saturation_pressure(mid) < target:
low = mid
else:
high = mid
return high
@dataclass
class ThermalSolver:
primary_volume_m3: float = 300.0
@@ -89,10 +102,37 @@ class ThermalSolver:
core.fuel_temperature,
)
def step_secondary(self, secondary: CoolantLoopState, transferred_mw: float) -> None:
delta_t = temperature_rise(transferred_mw, secondary.mass_flow_rate)
secondary.temperature_out = secondary.temperature_in + delta_t
secondary.steam_quality = min(1.0, max(0.0, delta_t / 100.0))
def step_secondary(self, secondary: CoolantLoopState, transferred_mw: float, dt: float = 1.0) -> None:
"""Update secondary loop using a simple steam-drum mass/energy balance."""
if transferred_mw <= 0.0 or secondary.mass_flow_rate <= 0.0:
secondary.steam_quality = max(0.0, secondary.steam_quality - 0.02 * dt)
secondary.temperature_out = max(
constants.ENVIRONMENT_TEMPERATURE, secondary.temperature_out - 0.5 * dt
)
secondary.pressure = max(
0.1, min(constants.MAX_PRESSURE, saturation_pressure(secondary.temperature_out))
)
return
temp_in = secondary.temperature_in
mass_flow = secondary.mass_flow_rate
cp = constants.COOLANT_HEAT_CAPACITY
sat_temp = saturation_temperature(max(0.05, secondary.pressure))
energy_j = max(0.0, transferred_mw) * constants.MEGAWATT * dt
# Energy needed to heat incoming feed to saturation.
sensible_j = max(0.0, sat_temp - temp_in) * mass_flow * cp * dt
if energy_j <= sensible_j:
delta_t = temperature_rise(transferred_mw, mass_flow)
secondary.temperature_out = temp_in + delta_t
secondary.steam_quality = 0.0
else:
energy_left = energy_j - sensible_j
steam_mass = energy_left / constants.STEAM_LATENT_HEAT
produced_fraction = steam_mass / max(1e-6, mass_flow * dt)
secondary.temperature_out = sat_temp
secondary.steam_quality = min(1.0, max(0.0, produced_fraction))
secondary.pressure = min(
constants.MAX_PRESSURE, max(secondary.pressure, saturation_pressure(secondary.temperature_out))
)

View File

@@ -43,3 +43,15 @@ def test_xenon_penalty_caps():
assert dynamics.xenon_penalty(state) == 0.05
state.xenon_inventory = 0.2
assert dynamics.xenon_penalty(state) == pytest.approx(0.01)
def test_delayed_precursors_follow_rod_banks():
dynamics = NeutronDynamics()
state = _core_state(power=600.0, flux=5e6)
dynamics.step(state, control_fraction=0.2, dt=1.0, rod_banks=[0.2, 0.2, 0.2])
initial_sum = sum(state.delayed_precursors)
assert len(state.delayed_precursors) == 3
assert initial_sum > 0.0
dynamics.step(state, control_fraction=0.95, dt=2.0, rod_banks=[0.95, 0.95, 0.95])
assert sum(state.delayed_precursors) < initial_sum

37
tests/test_thermal.py Normal file
View File

@@ -0,0 +1,37 @@
import pytest
from reactor_sim import constants
from reactor_sim.state import CoolantLoopState
from reactor_sim.thermal import ThermalSolver, saturation_temperature
def _secondary_loop(temp_in: float = 350.0, pressure: float = 0.5, flow: float = 200.0) -> CoolantLoopState:
nominal_mass = constants.SECONDARY_LOOP_VOLUME_M3 * constants.COOLANT_DENSITY
return CoolantLoopState(
temperature_in=temp_in,
temperature_out=temp_in,
pressure=pressure,
mass_flow_rate=flow,
steam_quality=0.0,
inventory_kg=nominal_mass,
level=constants.SECONDARY_INVENTORY_TARGET,
)
def test_secondary_heats_to_saturation_before_boiling():
solver = ThermalSolver()
loop = _secondary_loop(temp_in=330.0, flow=180.0, pressure=0.5)
sat_temp = saturation_temperature(loop.pressure)
solver.step_secondary(loop, transferred_mw=50.0, dt=1.0)
assert loop.steam_quality == pytest.approx(0.0)
assert loop.temperature_out < sat_temp
def test_secondary_generates_steam_when_energy_exceeds_sensible_heat():
solver = ThermalSolver()
loop = _secondary_loop(temp_in=330.0, flow=180.0, pressure=0.5)
sat_temp = saturation_temperature(loop.pressure)
solver.step_secondary(loop, transferred_mw=120.0, dt=1.0)
assert loop.temperature_out == pytest.approx(sat_temp, rel=0.02)
assert loop.steam_quality > 0.0
assert loop.steam_quality < 1.0