diff --git a/src/reactor_sim/reactor.py b/src/reactor_sim/reactor.py index b6be387..d61c36e 100644 --- a/src/reactor_sim/reactor.py +++ b/src/reactor_sim/reactor.py @@ -281,6 +281,8 @@ class Reactor: state.primary_loop.pressure, saturation_pressure(state.primary_loop.temperature_out), 0.1 ) target_pressure = max(0.5, base_head * power_ratio) + if self.primary_relief_open: + target_pressure = min(target_pressure, 1.0) primary_flow_scale = min( self._inventory_flow_scale(state.primary_loop), self._npsh_factor(state.primary_loop) ) @@ -335,6 +337,8 @@ class Reactor: demand = 0.75 base_flow, base_head = self.secondary_pump.performance(demand) target_pressure = max(0.5, base_head * power_ratio) + if self.secondary_relief_open: + target_pressure = min(target_pressure, 1.0) loop_pressure = max( state.secondary_loop.pressure, saturation_pressure(state.secondary_loop.temperature_out), 0.1 ) @@ -398,14 +402,25 @@ class Reactor: 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)) + self._vent_relief( + state.primary_loop, + target_pressure=1.0, + vent_rate_max=0.02, + ramp_time=12.0, + dt=dt, ) + for pump_state in state.primary_pumps: + pump_state.pressure = state.primary_loop.pressure 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_relief( + state.secondary_loop, + target_pressure=1.0, + vent_rate_max=0.05, + ramp_time=10.0, + dt=dt, ) - self._vent_secondary_relief(state.secondary_loop, dt) + for pump_state in state.secondary_pumps: + pump_state.pressure = state.secondary_loop.pressure 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) @@ -799,15 +814,21 @@ class Reactor: remaining = max(0.0, remaining - delivered) return total_power - def _vent_secondary_relief(self, loop: CoolantLoopState, dt: float) -> None: + def _vent_relief( + self, + loop: CoolantLoopState, + target_pressure: float, + vent_rate_max: float, + ramp_time: float, + 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_rate = min(vent_rate_max, 0.01 + vent_rate_max * 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 = ( @@ -816,9 +837,18 @@ class Reactor: ) 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) + # Pressure ramps toward target with the requested time constant. + ramp = min(1.0, dt / max(1e-6, ramp_time)) loop.pressure = max(target_pressure, loop.pressure - (loop.pressure - target_pressure) * ramp) + # Cool toward saturation at the new pressure to avoid re-pressurizing from superheat. + sat_target = saturation_temperature(target_pressure) + if loop.temperature_out > sat_target: + temp_drop = (loop.temperature_out - sat_target) * ramp + loop.temperature_out -= temp_drop + loop.temperature_in = min(loop.temperature_in, loop.temperature_out) + loop.steam_quality = 0.0 + cp = constants.COOLANT_HEAT_CAPACITY + loop.energy_j = max(0.0, loop.inventory_kg * cp * loop.average_temperature()) # Re-resolve temperature/quality/pressure to reflect the vented state. try: self.thermal._resolve_secondary_state(loop) # type: ignore[attr-defined]