diff --git a/src/reactor_sim/constants.py b/src/reactor_sim/constants.py index ac53b76..bcae0c3 100644 --- a/src/reactor_sim/constants.py +++ b/src/reactor_sim/constants.py @@ -17,6 +17,8 @@ ENVIRONMENT_TEMPERATURE = 295.0 # K AMU_TO_KG = 1.660_539_066_60e-27 MEV_TO_J = 1.602_176_634e-13 ELECTRON_FISSION_CROSS_SECTION = 5e-16 # cm^2, tuned for simulation scale +PUMP_SPOOL_TIME = 5.0 # seconds to reach commanded flow +TURBINE_SPOOL_TIME = 12.0 # seconds to reach steady output # Threshold inventories (event counts) for flagging common poisons in diagnostics. KEY_POISON_THRESHOLDS = { "Xe": 1e20, # xenon diff --git a/src/reactor_sim/coolant.py b/src/reactor_sim/coolant.py index 1894860..ce87b1d 100644 --- a/src/reactor_sim/coolant.py +++ b/src/reactor_sim/coolant.py @@ -5,6 +5,7 @@ from __future__ import annotations from dataclasses import dataclass import logging +from . import constants from .state import CoolantLoopState LOGGER = logging.getLogger(__name__) @@ -14,6 +15,7 @@ LOGGER = logging.getLogger(__name__) class Pump: nominal_flow: float efficiency: float = 0.9 + spool_time: float = constants.PUMP_SPOOL_TIME def flow_rate(self, demand: float) -> float: demand = max(0.0, min(1.0, demand)) diff --git a/src/reactor_sim/reactor.py b/src/reactor_sim/reactor.py index 41cdf24..6d2cd75 100644 --- a/src/reactor_sim/reactor.py +++ b/src/reactor_sim/reactor.py @@ -167,50 +167,78 @@ class Reactor: pump_demand = overrides.get("coolant_demand", self.control.coolant_demand(state.primary_loop)) if self.primary_pump_active: total_flow = 0.0 - pressure = 12.0 * pump_demand + 2.0 + target_pressure = 12.0 * pump_demand + 2.0 + loop_pressure = 0.5 + target_flow = self.primary_pump.flow_rate(pump_demand) for idx, pump_state in enumerate(state.primary_pumps): - if idx < len(self.primary_pump_units) and self.primary_pump_units[idx]: - flow = self.primary_pump.flow_rate(pump_demand) - pump_state.active = True - pump_state.flow_rate = flow - pump_state.pressure = pressure - total_flow += flow - else: - pump_state.active = False - pump_state.flow_rate = 0.0 - pump_state.pressure = state.primary_loop.pressure + unit_enabled = idx < len(self.primary_pump_units) and self.primary_pump_units[idx] + desired_flow = target_flow if unit_enabled else 0.0 + desired_pressure = target_pressure if unit_enabled else 0.5 + pump_state.flow_rate = self._ramp_value( + pump_state.flow_rate, desired_flow, dt, self.primary_pump.spool_time + ) + pump_state.pressure = self._ramp_value( + pump_state.pressure, desired_pressure, dt, self.primary_pump.spool_time + ) + pump_state.active = unit_enabled or pump_state.flow_rate > 1.0 + total_flow += pump_state.flow_rate + loop_pressure = max(loop_pressure, pump_state.pressure) state.primary_loop.mass_flow_rate = total_flow - state.primary_loop.pressure = pressure + state.primary_loop.pressure = loop_pressure if total_flow > 0 else self._ramp_value( + state.primary_loop.pressure, 0.5, dt, self.primary_pump.spool_time + ) else: - state.primary_loop.mass_flow_rate = 0.0 - state.primary_loop.pressure = 0.5 + state.primary_loop.mass_flow_rate = self._ramp_value( + state.primary_loop.mass_flow_rate, 0.0, dt, self.primary_pump.spool_time + ) + state.primary_loop.pressure = self._ramp_value( + state.primary_loop.pressure, 0.5, dt, self.primary_pump.spool_time + ) for pump_state in state.primary_pumps: pump_state.active = False - pump_state.flow_rate = 0.0 - pump_state.pressure = state.primary_loop.pressure + pump_state.flow_rate = self._ramp_value( + pump_state.flow_rate, 0.0, dt, self.primary_pump.spool_time + ) + pump_state.pressure = self._ramp_value( + pump_state.pressure, state.primary_loop.pressure, dt, self.primary_pump.spool_time + ) if self.secondary_pump_active: total_flow = 0.0 - pressure = 12.0 * 0.75 + 2.0 + target_pressure = 12.0 * 0.75 + 2.0 + loop_pressure = 0.5 + target_flow = self.secondary_pump.flow_rate(0.75) for idx, pump_state in enumerate(state.secondary_pumps): - if idx < len(self.secondary_pump_units) and self.secondary_pump_units[idx]: - flow = self.secondary_pump.flow_rate(0.75) - pump_state.active = True - pump_state.flow_rate = flow - pump_state.pressure = pressure - total_flow += flow - else: - pump_state.active = False - pump_state.flow_rate = 0.0 - pump_state.pressure = state.secondary_loop.pressure + unit_enabled = idx < len(self.secondary_pump_units) and self.secondary_pump_units[idx] + desired_flow = target_flow if unit_enabled else 0.0 + desired_pressure = target_pressure if unit_enabled else 0.5 + pump_state.flow_rate = self._ramp_value( + pump_state.flow_rate, desired_flow, dt, self.secondary_pump.spool_time + ) + pump_state.pressure = self._ramp_value( + pump_state.pressure, desired_pressure, dt, self.secondary_pump.spool_time + ) + pump_state.active = unit_enabled or pump_state.flow_rate > 1.0 + total_flow += pump_state.flow_rate + loop_pressure = max(loop_pressure, pump_state.pressure) state.secondary_loop.mass_flow_rate = total_flow - state.secondary_loop.pressure = pressure + state.secondary_loop.pressure = loop_pressure if total_flow > 0 else self._ramp_value( + state.secondary_loop.pressure, 0.5, dt, self.secondary_pump.spool_time + ) else: - state.secondary_loop.mass_flow_rate = 0.0 - state.secondary_loop.pressure = 0.5 + state.secondary_loop.mass_flow_rate = self._ramp_value( + state.secondary_loop.mass_flow_rate, 0.0, dt, self.secondary_pump.spool_time + ) + state.secondary_loop.pressure = self._ramp_value( + state.secondary_loop.pressure, 0.5, dt, self.secondary_pump.spool_time + ) for pump_state in state.secondary_pumps: pump_state.active = False - pump_state.flow_rate = 0.0 - pump_state.pressure = state.secondary_loop.pressure + pump_state.flow_rate = self._ramp_value( + pump_state.flow_rate, 0.0, dt, self.secondary_pump.spool_time + ) + pump_state.pressure = self._ramp_value( + pump_state.pressure, state.secondary_loop.pressure, dt, self.secondary_pump.spool_time + ) self.thermal.step_core(state.core, state.primary_loop, total_power, dt) if not self.secondary_pump_active or state.secondary_loop.mass_flow_rate <= 1.0: @@ -219,7 +247,7 @@ class Reactor: transferred = heat_transfer(state.primary_loop, state.secondary_loop, total_power) self.thermal.step_secondary(state.secondary_loop, transferred) - self._step_turbine_bank(state, transferred) + self._step_turbine_bank(state, transferred, dt) self._maintenance_tick(state, dt) if (not self.secondary_pump_active or state.secondary_loop.mass_flow_rate <= 1.0) and total_power > 50.0: @@ -253,7 +281,7 @@ class Reactor: sum(t.load_demand_mw for t in state.turbines), ) - def _step_turbine_bank(self, state: PlantState, steam_power_mw: float) -> None: + def _step_turbine_bank(self, state: PlantState, steam_power_mw: float, dt: float) -> None: if not state.turbines: return active_indices = [ @@ -265,9 +293,9 @@ class Reactor: break turbine_state = state.turbines[idx] if idx in active_indices: - turbine.step(state.secondary_loop, turbine_state, steam_power_mw=power_per_unit) + turbine.step(state.secondary_loop, turbine_state, steam_power_mw=power_per_unit, dt=dt) else: - self._reset_turbine_state(turbine_state) + self._spin_down_turbine(turbine_state, dt, turbine.spool_time) self._dispatch_consumer_load(state, active_indices) def _reset_turbine_state(self, turbine_state: TurbineState) -> None: @@ -276,6 +304,23 @@ class Reactor: turbine_state.load_demand_mw = 0.0 turbine_state.load_supplied_mw = 0.0 + @staticmethod + def _ramp_value(current: float, target: float, dt: float, time_constant: float) -> float: + if time_constant <= 0.0: + return target + alpha = min(1.0, max(0.0, dt / time_constant)) + return current + (target - current) * alpha + + def _spin_down_turbine(self, turbine_state: TurbineState, dt: float, time_constant: float) -> None: + turbine_state.shaft_power_mw = self._ramp_value(turbine_state.shaft_power_mw, 0.0, dt, time_constant) + turbine_state.electrical_output_mw = self._ramp_value( + turbine_state.electrical_output_mw, 0.0, dt, time_constant + ) + turbine_state.load_demand_mw = 0.0 + turbine_state.load_supplied_mw = self._ramp_value( + turbine_state.load_supplied_mw, 0.0, dt, time_constant + ) + def _dispatch_consumer_load(self, state: PlantState, active_indices: list[int]) -> None: total_electrical = sum(state.turbines[idx].electrical_output_mw for idx in active_indices) if self.consumer: diff --git a/src/reactor_sim/turbine.py b/src/reactor_sim/turbine.py index 266e77d..7c4c442 100644 --- a/src/reactor_sim/turbine.py +++ b/src/reactor_sim/turbine.py @@ -26,12 +26,14 @@ class Turbine: generator_efficiency: float = constants.GENERATOR_EFFICIENCY mechanical_efficiency: float = constants.STEAM_TURBINE_EFFICIENCY rated_output_mw: float = 400.0 # cap per unit electrical output + spool_time: float = constants.TURBINE_SPOOL_TIME def step( self, loop: CoolantLoopState, state: TurbineState, steam_power_mw: float = 0.0, + dt: float = 1.0, ) -> None: enthalpy = 2_700.0 + loop.steam_quality * 600.0 mass_flow = loop.mass_flow_rate * 0.6 @@ -43,8 +45,8 @@ class Turbine: shaft_power_mw = electrical / max(1e-6, self.generator_efficiency) condenser_temp = max(305.0, loop.temperature_in - 20.0) state.steam_enthalpy = enthalpy - state.shaft_power_mw = shaft_power_mw - state.electrical_output_mw = electrical + state.shaft_power_mw = _ramp(state.shaft_power_mw, shaft_power_mw, dt, self.spool_time) + state.electrical_output_mw = _ramp(state.electrical_output_mw, electrical, dt, self.spool_time) state.condenser_temperature = condenser_temp LOGGER.debug( "Turbine output: shaft=%.1fMW electrical=%.1fMW condenser=%.1fK", @@ -52,3 +54,10 @@ class Turbine: electrical, condenser_temp, ) + + +def _ramp(current: float, target: float, dt: float, time_constant: float) -> float: + if time_constant <= 0.0: + return target + alpha = min(1.0, max(0.0, dt / time_constant)) + return current + (target - current) * alpha diff --git a/tests/test_simulation.py b/tests/test_simulation.py index cc3b263..279b597 100644 --- a/tests/test_simulation.py +++ b/tests/test_simulation.py @@ -129,3 +129,18 @@ def test_secondary_pump_unit_toggle_can_restart_pump(): assert reactor.secondary_pump_units == [True, False] assert reactor.secondary_pump_active is True assert state.secondary_loop.mass_flow_rate > 0.0 + + +def test_primary_pumps_spool_up_over_seconds(): + reactor = Reactor.default() + state = reactor.initial_state() + # Enable both pumps and command full flow; spool should take multiple steps. + target_flow = reactor.primary_pump.flow_rate(1.0) * len(reactor.primary_pump_units) + reactor.step(state, dt=1.0, command=ReactorCommand(primary_pump_on=True, coolant_demand=1.0)) + first_flow = state.primary_loop.mass_flow_rate + assert 0.0 < first_flow < target_flow + + for _ in range(10): + reactor.step(state, dt=1.0, command=ReactorCommand(coolant_demand=1.0)) + + assert state.primary_loop.mass_flow_rate == pytest.approx(target_flow, rel=0.1) diff --git a/tests/test_turbine.py b/tests/test_turbine.py new file mode 100644 index 0000000..affa158 --- /dev/null +++ b/tests/test_turbine.py @@ -0,0 +1,33 @@ +import pytest + +from reactor_sim.state import CoolantLoopState, TurbineState +from reactor_sim.turbine import Turbine + + +def test_turbine_spools_toward_target_output(): + turbine = Turbine() + loop = CoolantLoopState( + temperature_in=600.0, + temperature_out=650.0, + pressure=6.0, + mass_flow_rate=20_000.0, + steam_quality=0.9, + ) + state = TurbineState( + steam_enthalpy=0.0, + shaft_power_mw=0.0, + electrical_output_mw=0.0, + condenser_temperature=300.0, + ) + target_electric = min( + turbine.rated_output_mw, 300.0 * turbine.mechanical_efficiency * turbine.generator_efficiency + ) + + dt = 5.0 + turbine.step(loop, state, steam_power_mw=300.0, dt=dt) + assert 0.0 < state.electrical_output_mw < target_electric + + for _ in range(5): + turbine.step(loop, state, steam_power_mw=300.0, dt=dt) + + assert state.electrical_output_mw == pytest.approx(target_electric, rel=0.05)