diff --git a/src/reactor_sim/constants.py b/src/reactor_sim/constants.py index 5ab5734..6b41202 100644 --- a/src/reactor_sim/constants.py +++ b/src/reactor_sim/constants.py @@ -12,6 +12,9 @@ STEAM_LATENT_HEAT = 2_200_000.0 # J/kg approximate latent heat of vaporization CORE_MELTDOWN_TEMPERATURE = 2_873.0 # K (approx 2600C) threshold for irreversible meltdown MAX_CORE_TEMPERATURE = CORE_MELTDOWN_TEMPERATURE # Allow simulation to approach meltdown temperature MAX_PRESSURE = 15.0 # MPa typical PWR primary loop limit +CLAD_MAX_TEMPERATURE = 1_200.0 # K clad softening / DNB concern +CHF_MASS_FLUX_REF = 1_500.0 # kg/m2-s reference mass flux surrogate +CHF_PRESSURE_REF_MPA = 7.0 # MPa reference pressure for CHF surrogate CONTROL_ROD_SPEED = 0.03 # fraction insertion per second CONTROL_ROD_WORTH = 0.042 # delta rho contribution when fully withdrawn CONTROL_ROD_BANK_WEIGHTS = (0.4, 0.35, 0.25) diff --git a/src/reactor_sim/dashboard.py b/src/reactor_sim/dashboard.py index e28cb75..9a25f19 100644 --- a/src/reactor_sim/dashboard.py +++ b/src/reactor_sim/dashboard.py @@ -541,7 +541,7 @@ class ReactorDashboard: "Toggle turbine units (1/2/3) for staggered maintenance.", "Use m/n/,/. for pump maintenance; B/V for generators.", "Press 'r' to reset/clear state if you want a cold start.", - "Watch component health to avoid automatic trips.", + "Watch component health, DNB margin, and subcooling to avoid automatic trips.", ] for idx, tip in enumerate(tips, start=y + 2): if not _add_safe(idx, 4, f"- {tip}"): diff --git a/src/reactor_sim/reactor.py b/src/reactor_sim/reactor.py index 74ec5e7..e21bff2 100644 --- a/src/reactor_sim/reactor.py +++ b/src/reactor_sim/reactor.py @@ -402,6 +402,12 @@ class Reactor: if (not self.secondary_pump_active or state.secondary_loop.mass_flow_rate <= 1.0) and total_power > 50.0: self._handle_heat_sink_loss(state) + if state.core.dnb_margin is not None and state.core.dnb_margin < 0.3: + LOGGER.critical("DNB margin low: %.2f, initiating SCRAM", state.core.dnb_margin) + self.shutdown = True + self.control.scram() + if state.core.subcooling_margin is not None and state.core.subcooling_margin < 5.0: + LOGGER.warning("Subcooling margin low: %.1fK", state.core.subcooling_margin) failures = self.health_monitor.evaluate( state, diff --git a/src/reactor_sim/state.py b/src/reactor_sim/state.py index 2c03fd3..4d710e9 100644 --- a/src/reactor_sim/state.py +++ b/src/reactor_sim/state.py @@ -21,6 +21,8 @@ class CoreState: power_output_mw: float # MW thermal burnup: float # fraction of fuel consumed clad_temperature: float | None = None # Kelvin + dnb_margin: float | None = None # ratio to critical heat flux surrogate + subcooling_margin: float | None = None # K until boiling xenon_inventory: float = 0.0 iodine_inventory: float = 0.0 delayed_precursors: list[float] = field(default_factory=list) @@ -30,6 +32,10 @@ class CoreState: def __post_init__(self) -> None: if self.clad_temperature is None: self.clad_temperature = self.fuel_temperature + if self.dnb_margin is None: + self.dnb_margin = 1.0 + if self.subcooling_margin is None: + self.subcooling_margin = 0.0 def update_burnup(self, dt: float) -> None: produced_energy_mwh = self.power_output_mw * (dt / 3600.0) diff --git a/src/reactor_sim/thermal.py b/src/reactor_sim/thermal.py index 78b6496..380e8c6 100644 --- a/src/reactor_sim/thermal.py +++ b/src/reactor_sim/thermal.py @@ -95,6 +95,11 @@ class ThermalSolver: core.fuel_temperature += heating - cooling # Keep fuel temperature bounded and never below the coolant outlet temperature. core.fuel_temperature = min(max(primary.temperature_out, core.fuel_temperature), constants.MAX_CORE_TEMPERATURE) + core.clad_temperature = max(primary.temperature_out, core.clad_temperature or primary.temperature_out) + core.subcooling_margin = max(0.0, saturation_temperature(primary.pressure) - primary.temperature_out) + chf = self._critical_heat_flux(primary) + heat_flux = (power_mw * constants.MEGAWATT) / max(1.0, self._core_surface_area()) + core.dnb_margin = max(0.0, chf / max(1e-6, heat_flux)) avg_temp = 0.5 * (primary.temperature_in + primary.temperature_out) primary.energy_j = max( 0.0, primary.inventory_kg * constants.COOLANT_HEAT_CAPACITY * avg_temp @@ -151,3 +156,17 @@ class ThermalSolver: secondary.steam_quality, secondary.energy_j, ) + + def _critical_heat_flux(self, primary: CoolantLoopState) -> float: + """Rough CHF surrogate using mass flux and pressure.""" + # Use a coarse mass-flux and pressure scaling to emulate higher CHF with more flow/pressure. + mass_flux = max(1.0, primary.mass_flow_rate / 50.0) # kg/m2-s surrogate + flux_factor = max(0.5, min(3.0, (mass_flux / 200.0) ** 0.5)) + pressure_factor = 0.5 + 0.5 * min(1.5, primary.pressure / max(0.1, constants.MAX_PRESSURE)) + base_chf = 1.0e7 # W/m2 surrogate + return base_chf * flux_factor * pressure_factor + + def _core_surface_area(self) -> float: + # Simple surrogate: area scaling with volume^(2/3) + volume = self.primary_volume_m3 + return max(1.0, (volume ** (2.0 / 3.0)) * 10.0)