Slow secondary relief venting and drop to 1 MPa

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)):