diff --git a/src/reactor_sim/reactor.py b/src/reactor_sim/reactor.py index 3e26918..b6be387 100644 --- a/src/reactor_sim/reactor.py +++ b/src/reactor_sim/reactor.py @@ -389,10 +389,6 @@ class Reactor: 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: - state.secondary_loop.pressure = max(0.1, saturation_pressure(state.secondary_loop.temperature_out)) if not self.secondary_pump_active or state.secondary_loop.mass_flow_rate <= 1.0: transferred = 0.0 @@ -401,6 +397,15 @@ class Reactor: residual = max(0.0, total_power - transferred) self.thermal.step_core(state.core, state.primary_loop, total_power, dt, residual_power_mw=residual) self.thermal.step_secondary(state.secondary_loop, transferred, dt) + if self.primary_relief_open: + state.primary_loop.pressure = max( + 0.1, min(state.primary_loop.pressure, saturation_pressure(state.primary_loop.temperature_out)) + ) + if self.secondary_relief_open: + state.secondary_loop.pressure = max( + 0.1, min(state.secondary_loop.pressure, saturation_pressure(state.secondary_loop.temperature_out)) + ) + self._vent_secondary_relief(state.secondary_loop, dt) if not self.control.manual_control and not self.shutdown: self.control.safety_backoff(state.core.subcooling_margin, state.core.dnb_margin, dt) self._apply_secondary_boiloff(state, dt) @@ -794,6 +799,32 @@ class Reactor: remaining = max(0.0, remaining - delivered) return total_power + def _vent_secondary_relief(self, loop: CoolantLoopState, dt: float) -> None: + """Model relief valve venting: gradual depressurization with mass/enthalpy loss.""" + target_pressure = 1.0 # MPa surrogate relief sink (not full vacuum) + if loop.inventory_kg <= 0.0: + loop.pressure = max(target_pressure, loop.pressure) + return + # Vent rate scales with overpressure; capped to keep a multi-second depressurization. + overfrac = max(0.0, (loop.pressure - target_pressure) / max(1e-6, loop.pressure)) + vent_rate = min(0.05, 0.01 + 0.1 * overfrac) # fraction of mass per second + vent_mass = min(loop.inventory_kg, loop.inventory_kg * vent_rate * dt) + if vent_mass > 0.0: + specific_enthalpy = ( + loop.steam_quality * constants.STEAM_LATENT_HEAT + + constants.COOLANT_HEAT_CAPACITY * max(loop.temperature_out, constants.ENVIRONMENT_TEMPERATURE) + ) + loop.inventory_kg = max(0.0, loop.inventory_kg - vent_mass) + loop.energy_j = max(0.0, loop.energy_j - vent_mass * specific_enthalpy) + # Pressure ramps toward target with a ~10s time constant. + ramp = min(1.0, dt / 10.0) + loop.pressure = max(target_pressure, loop.pressure - (loop.pressure - target_pressure) * ramp) + # Re-resolve temperature/quality/pressure to reflect the vented state. + try: + self.thermal._resolve_secondary_state(loop) # type: ignore[attr-defined] + except AttributeError: + pass + def _set_turbine_state(self, active: bool, index: int | None = None) -> None: if index is None: for idx in range(len(self.turbine_unit_active)):