From f6ff6fc618a6a8451ce56b37b705a5a29d4165de Mon Sep 17 00:00:00 2001 From: Codex Agent Date: Mon, 24 Nov 2025 00:06:08 +0100 Subject: [PATCH] Add delayed kinetics and steam-drum balance --- FEATURES.md | 4 +-- TODO.md | 11 ++++--- src/reactor_sim/constants.py | 1 + src/reactor_sim/neutronics.py | 62 ++++++++++++++++++++++++++++++++--- src/reactor_sim/reactor.py | 3 +- src/reactor_sim/state.py | 1 + src/reactor_sim/thermal.py | 48 ++++++++++++++++++++++++--- tests/test_neutronics.py | 12 +++++++ tests/test_thermal.py | 37 +++++++++++++++++++++ 9 files changed, 163 insertions(+), 16 deletions(-) create mode 100644 tests/test_thermal.py diff --git a/FEATURES.md b/FEATURES.md index c74c632..f6ab9f0 100644 --- a/FEATURES.md +++ b/FEATURES.md @@ -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. diff --git a/TODO.md b/TODO.md index f2d0b10..c197b41 100644 --- a/TODO.md +++ b/TODO.md @@ -1,7 +1,8 @@ ## 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. diff --git a/src/reactor_sim/constants.py b/src/reactor_sim/constants.py index d8df13d..d0a0efe 100644 --- a/src/reactor_sim/constants.py +++ b/src/reactor_sim/constants.py @@ -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 diff --git a/src/reactor_sim/neutronics.py b/src/reactor_sim/neutronics.py index 45d335f..e69f443 100644 --- a/src/reactor_sim/neutronics.py +++ b/src/reactor_sim/neutronics.py @@ -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 diff --git a/src/reactor_sim/reactor.py b/src/reactor_sim/reactor.py index 06ec045..c9a9b34 100644 --- a/src/reactor_sim/reactor.py +++ b/src/reactor_sim/reactor.py @@ -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 diff --git a/src/reactor_sim/state.py b/src/reactor_sim/state.py index 89014b4..3238f23 100644 --- a/src/reactor_sim/state.py +++ b/src/reactor_sim/state.py @@ -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) diff --git a/src/reactor_sim/thermal.py b/src/reactor_sim/thermal.py index 867fa13..3f3417c 100644 --- a/src/reactor_sim/thermal.py +++ b/src/reactor_sim/thermal.py @@ -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)) ) diff --git a/tests/test_neutronics.py b/tests/test_neutronics.py index 2c975d2..8826e20 100644 --- a/tests/test_neutronics.py +++ b/tests/test_neutronics.py @@ -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 diff --git a/tests/test_thermal.py b/tests/test_thermal.py new file mode 100644 index 0000000..b606ddc --- /dev/null +++ b/tests/test_thermal.py @@ -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