Add DNB/subcooling margins and keypad rod control

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)

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

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,

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)

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)