Tie loop pressure to saturation when cold

src/reactor_sim/reactor.py

@@ -16,7 +16,7 @@ from .fuel import FuelAssembly, decay_heat_fraction
 from .generator import DieselGenerator, GeneratorState
 from .neutronics import NeutronDynamics
 from .state import CoolantLoopState, CoreState, PlantState, PumpState, TurbineState
-from .thermal import ThermalSolver, heat_transfer
+from .thermal import ThermalSolver, heat_transfer, saturation_pressure
 from .turbine import SteamGenerator, Turbine
 
 LOGGER = logging.getLogger(__name__)
@@ -237,7 +237,7 @@ class Reactor:
             target_pressure = (
                 0.5 + (constants.PRIMARY_NOMINAL_PRESSURE - 0.5) * pump_demand
             ) * power_ratio
-            loop_pressure = 0.5
+            loop_pressure = max(0.1, saturation_pressure(state.primary_loop.temperature_out))
             target_flow = self.primary_pump.flow_rate(pump_demand) * power_ratio
             for idx, pump_state in enumerate(state.primary_pumps):
                 unit_enabled = (
@@ -270,7 +270,7 @@ class Reactor:
             )
         else:
             state.primary_loop.mass_flow_rate = 0.0
-            state.primary_loop.pressure = 0.5
+            state.primary_loop.pressure = max(0.1, saturation_pressure(state.primary_loop.temperature_out))
             for pump_state in state.primary_pumps:
                 pump_state.active = False
                 pump_state.flow_rate = 0.0
@@ -281,7 +281,7 @@ class Reactor:
             target_pressure = (
                 0.5 + (constants.SECONDARY_NOMINAL_PRESSURE - 0.5) * 0.75
             ) * power_ratio
-            loop_pressure = 0.5
+            loop_pressure = max(0.1, saturation_pressure(state.secondary_loop.temperature_out))
             target_flow = self.secondary_pump.flow_rate(0.75) * power_ratio
             for idx, pump_state in enumerate(state.secondary_pumps):
                 unit_enabled = (
@@ -314,7 +314,7 @@ class Reactor:
             )
         else:
             state.secondary_loop.mass_flow_rate = 0.0
-            state.secondary_loop.pressure = 0.5
+            state.secondary_loop.pressure = max(0.1, saturation_pressure(state.secondary_loop.temperature_out))
             for pump_state in state.secondary_pumps:
                 pump_state.active = False
                 pump_state.flow_rate = 0.0

src/reactor_sim/thermal.py

@@ -31,6 +31,16 @@ def temperature_rise(power_mw: float, mass_flow: float) -> float:
     return (power_mw * constants.MEGAWATT) / (mass_flow * constants.COOLANT_HEAT_CAPACITY)
 
 
+def saturation_pressure(temp_k: float) -> float:
+    """Return approximate saturation pressure (MPa) for water at temp_k."""
+    temp_c = max(0.0, temp_k - 273.15)
+    # Antoine equation parameters for water in 1-374C range, pressure in mmHg.
+    a, b, c = 8.14019, 1810.94, 244.485
+    psat_mmHg = 10 ** (a - (b / (c + temp_c)))
+    psat_mpa = psat_mmHg * 133.322e-6
+    return min(constants.MAX_PRESSURE, max(0.01, psat_mpa))
+
+
 @dataclass
 class ThermalSolver:
     primary_volume_m3: float = 300.0
@@ -55,7 +65,9 @@ class ThermalSolver:
         delta_t = temperature_rise(transferred_mw, secondary.mass_flow_rate)
         secondary.temperature_out = secondary.temperature_in + delta_t
         secondary.steam_quality = min(1.0, max(0.0, delta_t / 100.0))
-        secondary.pressure = min(constants.MAX_PRESSURE, 6.0 + delta_t * 0.01)
+        secondary.pressure = min(
+            constants.MAX_PRESSURE, max(secondary.pressure, saturation_pressure(secondary.temperature_out))
+        )
         LOGGER.debug(
             "Secondary loop: transferred=%.1fMW temp_out=%.1fK quality=%.2f",
             transferred_mw,