Compare commits
3 Commits
9dc4ca7733
...
03595b0d12
| Author | SHA1 | Date | |
|---|---|---|---|
|
|
03595b0d12 | ||
|
|
8ece3efb04 | ||
|
|
f6ff6fc618 |
@@ -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
12
TODO.md
@@ -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.
|
||||
|
||||
@@ -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
|
||||
|
||||
@@ -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)
|
||||
|
||||
@@ -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
|
||||
|
||||
@@ -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
|
||||
|
||||
@@ -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)
|
||||
|
||||
|
||||
@@ -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))
|
||||
)
|
||||
|
||||
@@ -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
37
tests/test_thermal.py
Normal 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
|
||||
Reference in New Issue
Block a user