diff --git a/FEATURES.md b/FEATURES.md index fb452f0..95a94de 100644 --- a/FEATURES.md +++ b/FEATURES.md @@ -4,6 +4,7 @@ - **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. +- **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, 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. - **Persistence & ops**: snapshots auto-save/load to `artifacts/last_state.json`; dashboard with live metrics, protections/warnings, heat-exchanger telemetry, component health, and control shortcuts. diff --git a/src/reactor_sim/constants.py b/src/reactor_sim/constants.py index 5478617..7ce78b4 100644 --- a/src/reactor_sim/constants.py +++ b/src/reactor_sim/constants.py @@ -35,6 +35,21 @@ SECONDARY_OUTLET_TARGET_K = 520.0 PRIMARY_NOMINAL_PRESSURE = 7.0 # MPa typical RBMK channel header pressure SECONDARY_NOMINAL_PRESSURE = 7.0 # MPa steam drum/steam line pressure surrogate STEAM_GENERATOR_UA_MW_PER_K = 25.0 # overall UA for steam generator (MW/K) +# Loop volume / inventory assumptions +PRIMARY_LOOP_VOLUME_M3 = 350.0 +SECONDARY_LOOP_VOLUME_M3 = 320.0 +PRIMARY_PRESSURIZER_SETPOINT_MPA = 7.0 +PRIMARY_PRESSURIZER_DEADBAND_MPA = 0.15 +PRIMARY_PRESSURIZER_HEAT_RATE_MPA_PER_S = 0.08 +PRIMARY_PRESSURIZER_SPRAY_RATE_MPA_PER_S = 0.12 +PRIMARY_PRESSURIZER_LEVEL_DRAW_PER_S = 0.002 +PRIMARY_PRESSURIZER_LEVEL_FILL_PER_S = 0.001 +LOOP_INVENTORY_CORRECTION_RATE = 0.08 # fraction of nominal mass restored per second toward target +PRIMARY_INVENTORY_TARGET = 0.95 +SECONDARY_INVENTORY_TARGET = 0.9 +SECONDARY_STEAM_LOSS_FRACTION = 0.02 # fraction of steam mass that leaves the loop each second +NPSH_REQUIRED_MPA = 0.25 +LOW_LEVEL_FLOW_FLOOR = 0.05 # Threshold inventories (event counts) for flagging common poisons in diagnostics. KEY_POISON_THRESHOLDS = { "Xe": 1e20, # xenon diff --git a/src/reactor_sim/dashboard.py b/src/reactor_sim/dashboard.py index a07f706..abfd4b9 100644 --- a/src/reactor_sim/dashboard.py +++ b/src/reactor_sim/dashboard.py @@ -384,9 +384,11 @@ class ReactorDashboard: "Flow", f"{state.primary_loop.mass_flow_rate:7.0f}/{self.reactor.primary_pump.nominal_flow * len(self.reactor.primary_pump_units):.0f} kg/s", ), + ("Level", f"{state.primary_loop.level*100:6.1f}%"), ("Inlet Temp", f"{state.primary_loop.temperature_in:7.1f} K"), ("Outlet Temp", f"{state.primary_loop.temperature_out:7.1f} K (Target {constants.PRIMARY_OUTLET_TARGET_K:4.0f})"), ("Pressure", f"{state.primary_loop.pressure:5.2f}/{constants.MAX_PRESSURE:4.1f} MPa"), + ("Pressurizer", f"{self.reactor.pressurizer_level*100:6.1f}% @ {constants.PRIMARY_PRESSURIZER_SETPOINT_MPA:4.1f} MPa"), ("Relief", "OPEN" if self.reactor.primary_relief_open else "CLOSED"), ], ) @@ -401,6 +403,7 @@ class ReactorDashboard: "Flow", 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})"), ("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"), diff --git a/src/reactor_sim/reactor.py b/src/reactor_sim/reactor.py index fb6c51c..c7cc23d 100644 --- a/src/reactor_sim/reactor.py +++ b/src/reactor_sim/reactor.py @@ -36,6 +36,7 @@ class Reactor: atomic_model: AtomicPhysics consumer: ElectricalConsumer | None = None health_monitor: HealthMonitor = field(default_factory=HealthMonitor) + pressurizer_level: float = 0.6 primary_pump_active: bool = True secondary_pump_active: bool = True primary_pump_units: list[bool] = field(default_factory=lambda: [True, True]) @@ -85,6 +86,9 @@ class Reactor: def initial_state(self) -> PlantState: ambient = constants.ENVIRONMENT_TEMPERATURE + primary_nominal_mass = constants.PRIMARY_LOOP_VOLUME_M3 * constants.COOLANT_DENSITY + secondary_nominal_mass = constants.SECONDARY_LOOP_VOLUME_M3 * constants.COOLANT_DENSITY + self.pressurizer_level = 0.6 core = CoreState( fuel_temperature=ambient, neutron_flux=1e5, @@ -119,6 +123,8 @@ class Reactor: pressure=0.5, mass_flow_rate=0.0, steam_quality=0.0, + inventory_kg=primary_nominal_mass * constants.PRIMARY_INVENTORY_TARGET, + level=constants.PRIMARY_INVENTORY_TARGET, ) secondary = CoolantLoopState( temperature_in=ambient, @@ -126,6 +132,8 @@ class Reactor: pressure=0.5, mass_flow_rate=0.0, steam_quality=0.0, + inventory_kg=secondary_nominal_mass * constants.SECONDARY_INVENTORY_TARGET, + level=constants.SECONDARY_INVENTORY_TARGET, ) primary_pumps = [ PumpState(active=self.primary_pump_active and self.primary_pump_units[idx], flow_rate=0.0, pressure=0.5) @@ -202,6 +210,13 @@ class Reactor: state.core.add_emitted_particles(particles) self._check_poison_alerts(state) + self._update_loop_inventory( + state.primary_loop, constants.PRIMARY_LOOP_VOLUME_M3, constants.PRIMARY_INVENTORY_TARGET, dt + ) + self._update_loop_inventory( + state.secondary_loop, constants.SECONDARY_LOOP_VOLUME_M3, constants.SECONDARY_INVENTORY_TARGET, dt + ) + pump_demand = overrides.get( "coolant_demand", self.control.coolant_demand( @@ -256,12 +271,15 @@ class Reactor: target_flow = base_flow * power_ratio loop_pressure = max(0.1, saturation_pressure(state.primary_loop.temperature_out)) target_pressure = max(0.5, base_head * power_ratio) + primary_flow_scale = min( + self._inventory_flow_scale(state.primary_loop), self._npsh_factor(state.primary_loop) + ) for idx, pump_state in enumerate(state.primary_pumps): unit_enabled = ( self.primary_pump_active and idx < len(self.primary_pump_units) and self.primary_pump_units[idx] ) powered = power_ratio > 0.1 - desired_flow = target_flow if unit_enabled else 0.0 + desired_flow = target_flow * primary_flow_scale if unit_enabled else 0.0 desired_pressure = target_pressure if unit_enabled else 0.5 if not powered: desired_flow = 0.0 @@ -275,6 +293,8 @@ class Reactor: pump_state.active = unit_enabled and powered and pump_state.flow_rate > 1.0 if not powered or not unit_enabled: pump_state.status = "STOPPING" if pump_state.flow_rate > 1.0 else "OFF" + elif primary_flow_scale < 0.99: + pump_state.status = "CAV" elif pump_state.flow_rate < max(1.0, desired_flow * 0.8): pump_state.status = "STARTING" else: @@ -289,8 +309,10 @@ class Reactor: state.primary_loop.mass_flow_rate = self._ramp_value( state.primary_loop.mass_flow_rate, 0.0, dt, self.primary_pump.spool_time ) + pressurizer_floor = constants.PRIMARY_PRESSURIZER_SETPOINT_MPA - 0.5 if self.pressurizer_level > 0.05 else 0.5 + target_pressure = max(pressurizer_floor, saturation_pressure(state.primary_loop.temperature_out)) state.primary_loop.pressure = self._ramp_value( - state.primary_loop.pressure, max(0.1, saturation_pressure(state.primary_loop.temperature_out)), dt, self.primary_pump.spool_time + state.primary_loop.pressure, target_pressure, dt, self.primary_pump.spool_time ) for pump_state in state.primary_pumps: pump_state.active = False @@ -305,12 +327,15 @@ class Reactor: target_pressure = max(0.5, base_head * power_ratio) loop_pressure = max(0.1, saturation_pressure(state.secondary_loop.temperature_out)) target_flow = base_flow * power_ratio + secondary_flow_scale = min( + self._inventory_flow_scale(state.secondary_loop), self._npsh_factor(state.secondary_loop) + ) for idx, pump_state in enumerate(state.secondary_pumps): unit_enabled = ( self.secondary_pump_active and idx < len(self.secondary_pump_units) and self.secondary_pump_units[idx] ) powered = power_ratio > 0.1 - desired_flow = target_flow if unit_enabled else 0.0 + desired_flow = target_flow * secondary_flow_scale if unit_enabled else 0.0 desired_pressure = target_pressure if unit_enabled else 0.5 if not powered: desired_flow = 0.0 @@ -324,6 +349,8 @@ class Reactor: pump_state.active = unit_enabled and powered and pump_state.flow_rate > 1.0 if not powered or not unit_enabled: pump_state.status = "STOPPING" if pump_state.flow_rate > 1.0 else "OFF" + elif secondary_flow_scale < 0.99: + pump_state.status = "CAV" elif pump_state.flow_rate < max(1.0, desired_flow * 0.8): pump_state.status = "STARTING" else: @@ -349,6 +376,7 @@ class Reactor: pump_state.pressure = state.secondary_loop.pressure pump_state.status = "STOPPING" if pump_state.flow_rate > 0.1 else "OFF" + self._apply_pressurizer(state.primary_loop, dt) if self.primary_relief_open: state.primary_loop.pressure = max(0.1, saturation_pressure(state.primary_loop.temperature_out)) if self.secondary_relief_open: @@ -360,6 +388,10 @@ class Reactor: 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._apply_secondary_boiloff(state, dt) + self._update_loop_inventory( + state.secondary_loop, constants.SECONDARY_LOOP_VOLUME_M3, constants.SECONDARY_INVENTORY_TARGET, dt + ) self._step_turbine_bank(state, transferred, dt) self._maintenance_tick(state, dt) @@ -487,6 +519,59 @@ class Reactor: turbine_state.load_demand_mw = demand_per_unit turbine_state.load_supplied_mw = share if demand_per_unit <= 0 else min(share, demand_per_unit) + def _nominal_inventory(self, volume_m3: float) -> float: + return volume_m3 * constants.COOLANT_DENSITY + + def _update_loop_inventory( + self, loop: CoolantLoopState, volume_m3: float, target_level: float, dt: float + ) -> None: + nominal_mass = self._nominal_inventory(volume_m3) + if nominal_mass <= 0.0: + loop.level = 0.0 + return + if loop.inventory_kg <= 0.0: + loop.inventory_kg = nominal_mass * target_level + current_level = loop.inventory_kg / nominal_mass + correction = (target_level - current_level) * constants.LOOP_INVENTORY_CORRECTION_RATE + loop.inventory_kg = max(0.0, loop.inventory_kg + correction * nominal_mass * dt) + loop.level = min(1.2, max(0.0, loop.inventory_kg / nominal_mass)) + + def _inventory_flow_scale(self, loop: CoolantLoopState) -> float: + if loop.level <= constants.LOW_LEVEL_FLOW_FLOOR: + return 0.0 + if loop.level <= 0.25: + return max(0.0, (loop.level - constants.LOW_LEVEL_FLOW_FLOOR) / (0.25 - constants.LOW_LEVEL_FLOW_FLOOR)) + return 1.0 + + def _npsh_factor(self, loop: CoolantLoopState) -> float: + vapor_pressure = saturation_pressure(loop.temperature_in) + available = max(0.0, loop.pressure - vapor_pressure) + if available <= 0.0: + return 0.0 + return max(0.0, min(1.0, available / constants.NPSH_REQUIRED_MPA)) + + def _apply_pressurizer(self, primary: CoolantLoopState, dt: float) -> None: + target = constants.PRIMARY_PRESSURIZER_SETPOINT_MPA + band = constants.PRIMARY_PRESSURIZER_DEADBAND_MPA + heat_rate = constants.PRIMARY_PRESSURIZER_HEAT_RATE_MPA_PER_S + spray_rate = constants.PRIMARY_PRESSURIZER_SPRAY_RATE_MPA_PER_S + if primary.pressure < target - band and self.pressurizer_level > 0.05: + primary.pressure = min(target, primary.pressure + heat_rate * dt) + self.pressurizer_level = max(0.0, self.pressurizer_level - constants.PRIMARY_PRESSURIZER_LEVEL_DRAW_PER_S * dt) + elif primary.pressure > target + band: + primary.pressure = max(target - band, primary.pressure - spray_rate * dt) + self.pressurizer_level = min(1.0, self.pressurizer_level + constants.PRIMARY_PRESSURIZER_LEVEL_FILL_PER_S * dt) + primary.pressure = min(constants.MAX_PRESSURE, max(saturation_pressure(primary.temperature_out), primary.pressure)) + + def _apply_secondary_boiloff(self, state: PlantState, dt: float) -> None: + loop = state.secondary_loop + if loop.mass_flow_rate <= 0.0 or loop.steam_quality <= 0.0: + return + steam_mass = loop.mass_flow_rate * loop.steam_quality * constants.SECONDARY_STEAM_LOSS_FRACTION * dt + loop.inventory_kg = max(0.0, loop.inventory_kg - steam_mass) + nominal = self._nominal_inventory(constants.SECONDARY_LOOP_VOLUME_M3) + loop.level = min(1.2, max(0.0, loop.inventory_kg / nominal)) if nominal > 0 else 0.0 + def _handle_failure(self, component: str) -> None: if component == "core": LOGGER.critical("Core failure detected. Initiating SCRAM.") @@ -765,6 +850,7 @@ class Reactor: "generator_auto": self.generator_auto, "primary_relief_open": self.primary_relief_open, "secondary_relief_open": self.secondary_relief_open, + "pressurizer_level": self.pressurizer_level, "maintenance_active": list(self.maintenance_active), "generators": [ { @@ -800,6 +886,7 @@ class Reactor: self.generator_auto = metadata.get("generator_auto", self.generator_auto) self.primary_relief_open = metadata.get("primary_relief_open", self.primary_relief_open) self.secondary_relief_open = metadata.get("secondary_relief_open", self.secondary_relief_open) + self.pressurizer_level = metadata.get("pressurizer_level", self.pressurizer_level) maint = metadata.get("maintenance_active") if maint is not None: self.maintenance_active = set(maint) diff --git a/src/reactor_sim/state.py b/src/reactor_sim/state.py index 717a095..89014b4 100644 --- a/src/reactor_sim/state.py +++ b/src/reactor_sim/state.py @@ -43,6 +43,8 @@ class CoolantLoopState: pressure: float # MPa mass_flow_rate: float # kg/s steam_quality: float # fraction of vapor + inventory_kg: float = 0.0 # bulk mass of coolant + level: float = 1.0 # fraction full relative to nominal volume def average_temperature(self) -> float: return 0.5 * (self.temperature_in + self.temperature_out) diff --git a/tests/test_pressurizer.py b/tests/test_pressurizer.py new file mode 100644 index 0000000..96950ca --- /dev/null +++ b/tests/test_pressurizer.py @@ -0,0 +1,35 @@ +from reactor_sim.reactor import Reactor +from reactor_sim.commands import ReactorCommand + + +def test_pressurizer_raises_pressure_with_level(): + reactor = Reactor.default() + state = reactor.initial_state() + state.primary_loop.pressure = 5.0 + reactor.pressurizer_level = 0.8 + reactor.primary_pump_active = False + reactor.secondary_pump_active = False + + reactor.step(state, dt=1.0, command=ReactorCommand.scram_all()) + + assert state.primary_loop.pressure > 5.0 + assert reactor.pressurizer_level < 0.8 + + +def test_low_npsh_limits_primary_flow(): + reactor = Reactor.default() + state = reactor.initial_state() + reactor.shutdown = False + reactor.control.manual_control = True + reactor.control.rod_fraction = 0.0 + reactor.primary_pump_active = True + reactor.primary_pump_units = [True, True] + reactor.secondary_pump_active = False + state.primary_loop.pressure = 0.05 # near-vacuum to force cavitation + state.primary_loop.temperature_in = 400.0 + state.primary_loop.temperature_out = 600.0 + + reactor.step(state, dt=1.0, command=ReactorCommand(generator_units={1: True})) + + assert state.primary_pumps[0].status == "CAV" + assert state.primary_loop.mass_flow_rate < 100.0