Stabilize secondary steam quality

src/reactor_sim/reactor.py

@@ -358,8 +358,7 @@ class Reactor:
             transferred = 0.0
         else:
             transferred = heat_transfer(state.primary_loop, state.secondary_loop, total_power)
-        net_power = total_power - transferred
-        self.thermal.step_core(state.core, state.primary_loop, net_power, dt)
+        self.thermal.step_core(state.core, state.primary_loop, total_power, dt)
         self.thermal.step_secondary(state.secondary_loop, transferred)
 
         self._step_turbine_bank(state, transferred, dt)

src/reactor_sim/thermal.py

@@ -15,7 +15,7 @@ LOGGER = logging.getLogger(__name__)
 
 def heat_transfer(primary: CoolantLoopState, secondary: CoolantLoopState, core_power_mw: float) -> float:
     """Return MW transferred to the secondary loop."""
-    if secondary.mass_flow_rate <= 0.0:
+    if primary.mass_flow_rate <= 0.0 or secondary.mass_flow_rate <= 0.0:
         return 0.0
     delta_t1 = max(1e-3, primary.temperature_out - secondary.temperature_in)
     delta_t2 = max(1e-3, primary.temperature_in - secondary.temperature_out)
@@ -26,7 +26,20 @@ def heat_transfer(primary: CoolantLoopState, secondary: CoolantLoopState, core_p
     else:
         lmtd = (delta_t1 - delta_t2) / math.log(delta_t1 / delta_t2)
     ua = constants.STEAM_GENERATOR_UA_MW_PER_K
-    transferred = max(0.0, min(core_power_mw, ua * lmtd))
+    ua_limited = ua * lmtd
+
+    # Prevent the heat exchanger from over-transferring and inverting the outlet temperatures.
+    primary_capacity = primary.mass_flow_rate * constants.COOLANT_HEAT_CAPACITY
+    secondary_capacity = secondary.mass_flow_rate * constants.COOLANT_HEAT_CAPACITY
+    approach = 2.0  # K minimum approach between loop outlets
+    pinch_limited = 0.0
+    if primary_capacity > 0.0 and secondary_capacity > 0.0:
+        temp_gap = primary.temperature_out - secondary.temperature_in - approach
+        if temp_gap > 0.0:
+            pinch_watts = temp_gap / ((1.0 / primary_capacity) + (1.0 / secondary_capacity))
+            pinch_limited = max(0.0, pinch_watts / constants.MEGAWATT)
+
+    transferred = max(0.0, min(core_power_mw, ua_limited, pinch_limited))
     LOGGER.debug("Heat transfer %.2f MW with LMTD=%.1fK (ΔT1=%.1f ΔT2=%.1f)", transferred, lmtd, delta_t1, delta_t2)
     return transferred
 
@@ -51,11 +64,20 @@ def saturation_pressure(temp_k: float) -> float:
 class ThermalSolver:
     primary_volume_m3: float = 300.0
 
-    def step_core(self, core: CoreState, primary: CoolantLoopState, power_mw: float, dt: float) -> None:
+    def step_core(
+        self,
+        core: CoreState,
+        primary: CoolantLoopState,
+        power_mw: float,
+        dt: float,
+        residual_power_mw: float | None = None,
+    ) -> None:
+        if residual_power_mw is None:
+            residual_power_mw = power_mw
         temp_rise = temperature_rise(power_mw, primary.mass_flow_rate)
         primary.temperature_out = primary.temperature_in + temp_rise
         # Fuel heats from any power not immediately convected away, and cools toward the primary outlet.
-        heating = 0.005 * max(0.0, power_mw - temp_rise) * dt
+        heating = 0.005 * 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
         # Keep fuel temperature bounded and never below the coolant outlet temperature.