diff --git a/src/reactor_sim/constants.py b/src/reactor_sim/constants.py index b9376e5..7065a12 100644 --- a/src/reactor_sim/constants.py +++ b/src/reactor_sim/constants.py @@ -78,11 +78,17 @@ HX_FOULING_RATE = 1e-5 # fouling increment per second scaled by impurities/temp HX_FOULING_HEAL_RATE = 5e-6 # cleaning/settling when cool/low steam HX_FOULING_MAX_PENALTY = 0.25 # fractional UA loss cap BORON_WORTH_PER_PPM = 8e-6 # delta rho per ppm relative to baseline boron -BORON_TRIM_RATE_PPM_PER_S = 0.02 # slow boron trim toward setpoint when near target +BORON_TRIM_RATE_PPM_PER_S = 0.04 # slow boron trim toward setpoint when near target # Mild thermal/measurement lags FUEL_TO_CLAD_TIME_CONSTANT = 0.3 # seconds (mild lag) CLAD_TO_COOLANT_TIME_CONSTANT = 0.2 # seconds (mild lag) -POWER_MEASUREMENT_TIME_CONSTANT = 6.0 # seconds +POWER_MEASUREMENT_TIME_CONSTANT = 10.0 # seconds +# Passive cooldown +PASSIVE_COOL_RATE_PRIMARY = 0.05 # K/s toward ambient when low/no transfer +PASSIVE_COOL_RATE_SECONDARY = 0.08 # K/s toward ambient when no steam/heat sink +RHR_ACTIVE = True +RHR_CUTOFF_POWER_MW = 5.0 +RHR_COOL_RATE = 0.2 # K/s forced cooldown when power is near zero # Threshold inventories (event counts) for flagging common poisons in diagnostics. KEY_POISON_THRESHOLDS = { "Xe": 1e20, # xenon diff --git a/src/reactor_sim/control.py b/src/reactor_sim/control.py index 50ebc78..7660e13 100644 --- a/src/reactor_sim/control.py +++ b/src/reactor_sim/control.py @@ -25,6 +25,10 @@ class ControlSystem: rod_banks: list[float] = field(default_factory=lambda: [0.5, 0.5, 0.5]) rod_target: float = 0.5 _filtered_power_mw: float = 0.0 + _integral_error: float = 0.0 + _ramp_start_mw: float = 0.0 + _ramp_progress_mw: float = 0.0 + _last_manual: bool = False def update_rods(self, state: CoreState, dt: float) -> float: if not self.rod_banks or len(self.rod_banks) != len(constants.CONTROL_ROD_BANK_WEIGHTS): @@ -35,7 +39,27 @@ class ControlSystem: if abs(self.rod_fraction - self.effective_insertion()) > 1e-6: self.rod_target = clamp(self.rod_fraction, 0.0, 0.95) self._advance_banks(self.rod_target, dt) + self._last_manual = True return self.rod_fraction + + # On transition from manual -> auto, start a gentle setpoint ramp from current power. + if self._last_manual: + self._ramp_start_mw = max(0.0, state.power_output_mw) + self._ramp_progress_mw = 0.0 + self._last_manual = False + + # Setpoint ramp: close the gap at a limited MW/s. + ramp_rate = 5.0 # MW per second toward setpoint + effective_setpoint = self.setpoint_mw + if self._ramp_start_mw > 0.0 and self.setpoint_mw > self._ramp_start_mw: + self._ramp_progress_mw = min( + self.setpoint_mw - self._ramp_start_mw, + self._ramp_progress_mw + ramp_rate * dt, + ) + effective_setpoint = self._ramp_start_mw + self._ramp_progress_mw + else: + effective_setpoint = self.setpoint_mw + raw_power = state.power_output_mw # Begin filtering once we're in the vicinity of the setpoint to avoid chasing noise. if self._filtered_power_mw <= 0.0: @@ -48,14 +72,35 @@ class ControlSystem: else: measured_power = raw_power - error = (measured_power - self.setpoint_mw) / self.setpoint_mw + error = (measured_power - effective_setpoint) / max(1e-6, effective_setpoint) + near = measured_power > 0.9 * effective_setpoint # Deadband near setpoint to prevent dithering; only apply when we're close to target. - if measured_power > 0.9 * self.setpoint_mw and abs(error) < 0.01: + if near and abs(error) < 0.01: self._advance_banks(self.rod_target, dt) return self.rod_fraction - # When power is low (negative error) withdraw rods; when high, insert them. - adjustment = error * 0.35 - adjustment = clamp(adjustment, -constants.CONTROL_ROD_SPEED * dt, constants.CONTROL_ROD_SPEED * dt) + # Integrate a bit to remove steady-state error. + self._integral_error = clamp(self._integral_error + error * dt, -0.05, 0.05) + p_gain = 0.35 if not near else 0.2 + i_gain = 0.015 if not near else 0.01 + adjustment = p_gain * error + i_gain * self._integral_error + speed = constants.CONTROL_ROD_SPEED * dt + # Hard high-band clamp: above 5% over setpoint, freeze withdrawals; above 10%, insert faster. + if measured_power > 1.1 * effective_setpoint: + speed *= 0.4 + adjustment = -speed # force insertion at capped rate + elif measured_power > 1.05 * effective_setpoint: + adjustment = min(adjustment, 0.0) + speed *= 0.6 + # Soft clamp above setpoint: if we're >5% high, enforce insertion at capped rate. + if error > 0.1: + adjustment = min(adjustment, 0.0) + speed *= 0.5 + elif error > 0.05: + adjustment = min(adjustment, adjustment) + speed *= 0.7 + if near: + speed *= 0.5 # slow down when near setpoint to avoid overshoot + adjustment = clamp(adjustment, -speed, speed) self.rod_target = clamp(self.rod_target + adjustment, 0.0, 0.95) self._advance_banks(self.rod_target, dt) LOGGER.debug("Control rod target=%.3f (error=%.3f)", self.rod_target, error) @@ -98,7 +143,7 @@ class ControlSystem: severity = max(severity, max(0.0, 0.7 - dnb_margin) / 0.7) if severity <= 0.0: return - backoff = (0.003 + 0.02 * severity) * dt + backoff = (0.001 + 0.01 * severity) * dt self.rod_target = clamp(self.rod_target + backoff, 0.0, 0.95) self._advance_banks(self.rod_target, dt) LOGGER.debug("Safety backoff applied: target=%.3f severity=%.2f", self.rod_target, severity) diff --git a/src/reactor_sim/dashboard.py b/src/reactor_sim/dashboard.py index 4246f01..e5c56c3 100644 --- a/src/reactor_sim/dashboard.py +++ b/src/reactor_sim/dashboard.py @@ -750,6 +750,10 @@ class ReactorDashboard: lines.append(("SCRAM", "ACTIVE" if self.reactor.shutdown else "CLEAR", curses.color_pair(4) if self.reactor.shutdown else 0)) if self.reactor.meltdown: lines.append(("Meltdown", "IN PROGRESS", curses.color_pair(4) | curses.A_BOLD)) + cooldown_status = "Normal" + if state.core.power_output_mw <= constants.RHR_CUTOFF_POWER_MW and (self.reactor.primary_pump_active or self.reactor.secondary_pump_active): + cooldown_status = "RHR/Passive" + lines.append(("Cooldown", cooldown_status)) sec_flow_low = state.secondary_loop.mass_flow_rate <= 1.0 or not self.reactor.secondary_pump_active heat_sink_risk = sec_flow_low and state.core.power_output_mw > 50.0 if heat_sink_risk: diff --git a/src/reactor_sim/reactor.py b/src/reactor_sim/reactor.py index ab73f94..f82b3b8 100644 --- a/src/reactor_sim/reactor.py +++ b/src/reactor_sim/reactor.py @@ -209,7 +209,7 @@ class Reactor: ) decay_heat = decay_heat_fraction(state.core.burnup) * state.core.power_output_mw total_power = prompt_power + decay_heat + decay_power - total_power = min(total_power, constants.TEST_MAX_POWER_MW * 0.98) + total_power = min(total_power, constants.TEST_MAX_POWER_MW * 0.98, self.control.setpoint_mw * 1.1) state.core.power_output_mw = total_power state.core.update_burnup(dt) # Track fission products and emitted particles for diagnostics. @@ -484,12 +484,20 @@ class Reactor: env = constants.ENVIRONMENT_TEMPERATURE primary_cooling = temperature_rise(transferred, state.primary_loop.mass_flow_rate) if transferred <= 0.0 or state.secondary_loop.mass_flow_rate <= 1.0: - passive = 0.02 * max(0.0, state.primary_loop.temperature_out - env) * dt + passive = constants.PASSIVE_COOL_RATE_PRIMARY * max(0.0, state.primary_loop.temperature_out - env) * dt primary_cooling = max(primary_cooling, passive) state.primary_loop.temperature_in = max(env, state.primary_loop.temperature_out - primary_cooling) + if state.core.power_output_mw <= constants.RHR_CUTOFF_POWER_MW and not self.turbine_active: + bleed = constants.RHR_COOL_RATE * dt + state.primary_loop.temperature_out = max(env, state.primary_loop.temperature_out - bleed) + state.primary_loop.temperature_in = max(env, state.primary_loop.temperature_out - bleed) if state.secondary_loop.mass_flow_rate <= 1.0: - target_temp = env + # Passive cooldown toward ambient when pumps off/low steam. + rhr = 0.0 + if constants.RHR_ACTIVE and state.core.power_output_mw <= constants.RHR_CUTOFF_POWER_MW: + rhr = constants.RHR_COOL_RATE * dt + target_temp = max(env, state.secondary_loop.temperature_out - rhr) state.secondary_loop.temperature_out = self._ramp_value( state.secondary_loop.temperature_out, target_temp, dt, self.secondary_pump.spool_time ) @@ -499,6 +507,10 @@ class Reactor: excess = max(0.0, state.secondary_loop.temperature_out - env) cooling_drop = min(40.0, max(10.0, 0.2 * excess)) state.secondary_loop.temperature_in = max(env, state.secondary_loop.temperature_out - cooling_drop) + if state.core.power_output_mw <= constants.RHR_CUTOFF_POWER_MW and not self.turbine_active: + bleed = constants.RHR_COOL_RATE * dt + state.secondary_loop.temperature_out = max(env, state.secondary_loop.temperature_out - bleed) + state.secondary_loop.temperature_in = max(env, state.secondary_loop.temperature_out - bleed) # Keep stored energies consistent with updated temperatures/quality. cp = constants.COOLANT_HEAT_CAPACITY @@ -675,7 +687,7 @@ class Reactor: if self.control.setpoint_mw <= 0.0: return error = (state.core.power_output_mw - self.control.setpoint_mw) / self.control.setpoint_mw - if abs(error) < 0.02: + if abs(error) < 0.005: return delta = constants.BORON_TRIM_RATE_PPM_PER_S * error * dt state.boron_ppm = min(constants.CHEM_MAX_PPM, max(0.0, state.boron_ppm + delta)) diff --git a/src/reactor_sim/thermal.py b/src/reactor_sim/thermal.py index 0f07f79..bcb43bc 100644 --- a/src/reactor_sim/thermal.py +++ b/src/reactor_sim/thermal.py @@ -127,7 +127,7 @@ class ThermalSolver: temp_rise = temperature_rise(power_mw, primary.mass_flow_rate) primary.temperature_out = primary.temperature_in + temp_rise # Fuel heats from total fission power (even when most is convected) plus any residual left in the coolant. - heating = (0.002 * max(0.0, power_mw) + 0.01 * max(0.0, residual_power_mw)) * dt + heating = (0.003 * max(0.0, power_mw) + 0.01 * max(0.0, residual_power_mw)) * dt cooling = 0.025 * max(0.0, core.fuel_temperature - primary.temperature_out) * dt core.fuel_temperature += heating - cooling # Simple radial split: fuel -> clad conduction with burnup-driven conductivity drop and gap penalty. diff --git a/tests/test_simulation.py b/tests/test_simulation.py index d524df6..098ecd5 100644 --- a/tests/test_simulation.py +++ b/tests/test_simulation.py @@ -341,3 +341,38 @@ def test_full_power_reaches_steam_and_turbine_output(): assert results[3600]["electric"] > 150.0 # No runaway core temperatures. assert results[3600]["core_temp"] < constants.CORE_MELTDOWN_TEMPERATURE + + +def test_cooldown_reaches_ambient_without_runaway(): + """Shutdown with pumps running should cool loops toward ambient, no runaway.""" + reactor = Reactor.default() + reactor.health_monitor.disable_degradation = True + reactor.allow_external_aux = True + reactor.relaxed_npsh = True + state = reactor.initial_state() + # Start hot. + reactor.step( + state, + dt=1.0, + command=ReactorCommand( + generator_units={1: True, 2: True}, + primary_pumps={1: True, 2: True}, + secondary_pumps={1: True, 2: True}, + rod_manual=True, + rod_position=0.55, + ), + ) + turbines_started = False + for i in range(1800): + cmd = None + if not turbines_started and state.secondary_loop.steam_quality > 0.02 and state.secondary_loop.pressure > 1.0: + cmd = ReactorCommand(turbine_on=True, turbine_units={1: True, 2: True, 3: True}) + turbines_started = True + if i == 900: + cmd = ReactorCommand(rod_position=0.95, turbine_on=False, turbine_units={1: False, 2: False, 3: False}) + reactor.step(state, dt=1.0, command=cmd) + + assert not reactor.meltdown + assert state.core.power_output_mw < 1.0 + assert state.primary_loop.temperature_out < 320.0 + assert state.secondary_loop.temperature_out < 320.0