Add enthalpy-based secondary boil-off and turbine mapping

TODO.md

@@ -8,7 +8,7 @@
 - [ ] Flesh out condenser behavior: vacuum pump limits, cooling water temperature coupling, and dynamic back-pressure with fouling.
 - [ ] Dashboard polish: compact turbine/generator rows, color critical warnings (SCRAM/heat-sink), and reduce repeated log noise.
 - [ ] Incremental realism plan:
-  - Add stored enthalpy for primary/secondary loops and a steam-drum mass/energy balance (sensible + latent) while keeping existing pump logic and tests passing.
+  - Add stored enthalpy for primary/secondary loops and a steam-drum mass/energy balance (sensible + latent) while keeping existing pump logic and tests passing. Target representative PWR conditions: primary 15–16 MPa, 290–320 °C inlet/320–330 °C outlet, secondary saturation ~6–7 MPa with boil at ~490–510 K.
-  - Adjust HX/pressure handling to use stored energy (saturation clamp and pressure rise) and validate steam formation with both pumps at ~3 GW.
+  - Adjust HX/pressure handling to use stored energy (saturation clamp and pressure rise) and validate steam formation with both pumps at ~3 GW. Use realistic tube-side material assumptions (Inconel 690/SS cladding) and clamp steam quality to phase-equilibrium enthalpy.
-  - Update turbine power mapping to consume steam enthalpy/quality and align protection trips with real steam presence.
+  - Update turbine power mapping to consume steam enthalpy/quality and align protection trips with real steam presence; drive inlet steam around 6–7 MPa, quality/enthalpy-based flow to ~550–600 MW(e) per machine class if steam is available.
-  - Add integration test: cold start → gens/pumps 2/2 → ramp to ~3 GW → confirm steam quality threshold → enable all turbines and require electrical output.
+  - Add integration test: cold start → gens/pumps 2/2 → ramp to ~3 GW → confirm steam quality threshold at the secondary drum → enable all turbines and require electrical output. Include a step that tolerates one secondary pump off for a period to prove redundancy still yields steam.

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, saturation_pressure, temperature_rise
+from .thermal import ThermalSolver, heat_transfer, saturation_pressure, saturation_temperature, temperature_rise
 from .turbine import SteamGenerator, Turbine
 
 LOGGER = logging.getLogger(__name__)
@@ -126,6 +126,7 @@ class Reactor:
             steam_quality=0.0,
             inventory_kg=primary_nominal_mass * constants.PRIMARY_INVENTORY_TARGET,
             level=constants.PRIMARY_INVENTORY_TARGET,
+            energy_j=primary_nominal_mass * constants.PRIMARY_INVENTORY_TARGET * constants.COOLANT_HEAT_CAPACITY * ambient,
         )
         secondary = CoolantLoopState(
             temperature_in=ambient,
@@ -135,6 +136,7 @@ class Reactor:
             steam_quality=0.0,
             inventory_kg=secondary_nominal_mass * constants.SECONDARY_INVENTORY_TARGET,
             level=constants.SECONDARY_INVENTORY_TARGET,
+            energy_j=secondary_nominal_mass * constants.SECONDARY_INVENTORY_TARGET * constants.COOLANT_HEAT_CAPACITY * ambient,
         )
         primary_pumps = [
             PumpState(active=self.primary_pump_active and self.primary_pump_units[idx], flow_rate=0.0, pressure=0.5)
@@ -434,6 +436,17 @@ class Reactor:
             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)
 
+        # Keep stored energies consistent with updated temperatures/quality.
+        cp = constants.COOLANT_HEAT_CAPACITY
+        primary_avg = 0.5 * (state.primary_loop.temperature_in + state.primary_loop.temperature_out)
+        state.primary_loop.energy_j = max(0.0, state.primary_loop.inventory_kg * cp * primary_avg)
+        sat_temp_sec = saturation_temperature(max(0.05, state.secondary_loop.pressure))
+        sec_liquid_energy = state.secondary_loop.inventory_kg * cp * min(state.secondary_loop.temperature_out, sat_temp_sec)
+        sec_latent = state.secondary_loop.inventory_kg * state.secondary_loop.steam_quality * constants.STEAM_LATENT_HEAT
+        superheat = max(0.0, state.secondary_loop.temperature_out - sat_temp_sec)
+        sec_superheat = state.secondary_loop.inventory_kg * cp * superheat if state.secondary_loop.steam_quality >= 1.0 else 0.0
+        state.secondary_loop.energy_j = max(0.0, sec_liquid_energy + sec_latent + sec_superheat)
+
         state.primary_to_secondary_delta_t = max(0.0, state.primary_loop.temperature_out - state.secondary_loop.temperature_in)
         state.heat_exchanger_efficiency = 0.0 if total_power <= 0 else min(1.0, max(0.0, transferred / max(1e-6, total_power)))
 
@@ -578,7 +591,13 @@ class Reactor:
         if loop.mass_flow_rate <= 0.0 or loop.steam_quality <= 0.0:
             return
         steam_mass = loop.mass_flow_rate * loop.steam_quality * constants.SECONDARY_STEAM_LOSS_FRACTION * dt
+        if steam_mass <= 0.0:
+            return
+        prev_mass = max(1e-6, loop.inventory_kg)
         loop.inventory_kg = max(0.0, loop.inventory_kg - steam_mass)
+        # Scale stored energy with the remaining mass to keep specific enthalpy consistent.
+        ratio = max(0.0, loop.inventory_kg) / prev_mass
+        loop.energy_j *= ratio
         nominal = self._nominal_inventory(constants.SECONDARY_LOOP_VOLUME_M3)
         loop.level = min(1.2, max(0.0, loop.inventory_kg / nominal)) if nominal > 0 else 0.0
 

src/reactor_sim/state.py

@@ -4,6 +4,8 @@ from __future__ import annotations
 
 from dataclasses import dataclass, field, asdict
 
+from . import constants
+
 from .generator import GeneratorState
 
 
@@ -51,6 +53,7 @@ class CoolantLoopState:
     steam_quality: float  # fraction of vapor
     inventory_kg: float = 0.0  # bulk mass of coolant
     level: float = 1.0  # fraction full relative to nominal volume
+    energy_j: float = 0.0  # stored thermal/latent energy for the loop
 
     def average_temperature(self) -> float:
         return 0.5 * (self.temperature_in + self.temperature_out)
@@ -135,8 +138,8 @@ class PlantState:
         delta_t = data.get("primary_to_secondary_delta_t", 0.0)
         return cls(
             core=CoreState(**core_blob, fission_product_inventory=inventory, emitted_particles=particles),
-            primary_loop=CoolantLoopState(**data["primary_loop"]),
+            primary_loop=CoolantLoopState(**_with_energy(data["primary_loop"])),
-            secondary_loop=CoolantLoopState(**data["secondary_loop"]),
+            secondary_loop=CoolantLoopState(**_with_energy(data["secondary_loop"])),
             turbines=turbines,
             primary_pumps=[PumpState(**p) for p in prim_pumps_blob],
             secondary_pumps=[PumpState(**p) for p in sec_pumps_blob],
@@ -146,3 +149,14 @@ class PlantState:
             primary_to_secondary_delta_t=delta_t,
             time_elapsed=data.get("time_elapsed", 0.0),
         )
+
+
+def _with_energy(loop_blob: dict) -> dict:
+    """Backwards compatibility: derive energy if missing."""
+    if "energy_j" in loop_blob:
+        return loop_blob
+    energy = 0.5 * (loop_blob.get("temperature_in", 295.0) + loop_blob.get("temperature_out", 295.0))
+    energy *= loop_blob.get("inventory_kg", 0.0) * constants.COOLANT_HEAT_CAPACITY
+    out = dict(loop_blob)
+    out["energy_j"] = energy
+    return out

src/reactor_sim/thermal.py

@@ -95,6 +95,10 @@ 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)
+        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
+        )
         LOGGER.debug(
             "Primary loop: flow=%.0f kg/s temp_out=%.1fK core_temp=%.1fK",
             primary.mass_flow_rate,
@@ -103,42 +107,47 @@ class ThermalSolver:
         )
 
     def step_secondary(self, secondary: CoolantLoopState, transferred_mw: float, dt: float = 1.0) -> None:
-        """Update secondary loop using a simple steam-drum mass/energy balance."""
+        """Update secondary loop using a stored-energy steam-drum balance."""
-        if transferred_mw <= 0.0 or secondary.mass_flow_rate <= 0.0:
-            secondary.steam_quality = max(0.0, secondary.steam_quality - 0.02 * dt)
-            secondary.temperature_out = max(
-                constants.ENVIRONMENT_TEMPERATURE, secondary.temperature_out - 0.5 * dt
-            )
-            secondary.pressure = max(
-                0.1, min(constants.MAX_PRESSURE, saturation_pressure(secondary.temperature_out))
-            )
-            return
-
-        temp_in = secondary.temperature_in
-        mass_flow = secondary.mass_flow_rate
         cp = constants.COOLANT_HEAT_CAPACITY
-        sat_temp = saturation_temperature(max(0.05, secondary.pressure))
+        mass = max(1e-6, secondary.inventory_kg)
-        energy_j = max(0.0, transferred_mw) * constants.MEGAWATT * dt
+        if secondary.energy_j <= 0.0:
+            secondary.energy_j = mass * cp * secondary.average_temperature()
 
-        # Energy needed to heat incoming feed to saturation.
+        # Add transferred heat; if no heat, bleed toward ambient.
-        sensible_j = max(0.0, sat_temp - temp_in) * mass_flow * cp * dt
+        if transferred_mw > 0.0:
-        if energy_j <= sensible_j:
+            secondary.energy_j += transferred_mw * constants.MEGAWATT * dt
-            delta_t = temperature_rise(transferred_mw, mass_flow)
-            secondary.temperature_out = temp_in + delta_t
-            secondary.steam_quality = 0.0
         else:
-            energy_left = energy_j - sensible_j
+            bleed = 0.01 * (secondary.temperature_out - constants.ENVIRONMENT_TEMPERATURE)
-            steam_mass = energy_left / constants.STEAM_LATENT_HEAT
+            secondary.energy_j = max(
-            produced_fraction = steam_mass / max(1e-6, mass_flow * dt)
+                mass * cp * constants.ENVIRONMENT_TEMPERATURE, secondary.energy_j - max(0.0, bleed) * mass * cp * dt
-            secondary.temperature_out = sat_temp
+            )
-            secondary.steam_quality = min(1.0, max(0.0, produced_fraction))
+
+        sat_temp = saturation_temperature(max(0.05, secondary.pressure))
+        liquid_energy = mass * cp * sat_temp
+        available = secondary.energy_j
+
+        if available <= liquid_energy:
+            # Subcooled or saturated liquid.
+            temp = available / (mass * cp)
+            secondary.temperature_out = max(temp, constants.ENVIRONMENT_TEMPERATURE)
+            secondary.steam_quality = max(0.0, secondary.steam_quality - 0.01 * dt)
+        else:
+            excess = available - liquid_energy
+            quality = min(1.0, excess / (mass * constants.STEAM_LATENT_HEAT))
+            superheat_energy = max(0.0, excess - quality * mass * constants.STEAM_LATENT_HEAT)
+            superheat_temp = superheat_energy / (mass * cp) if quality >= 1.0 else 0.0
+            secondary.temperature_out = sat_temp + superheat_temp
+            secondary.steam_quality = quality
+            # Re-normalize stored energy to the realized state.
+            secondary.energy_j = liquid_energy + quality * mass * constants.STEAM_LATENT_HEAT + superheat_energy
 
         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",
+            "Secondary loop: transferred=%.1fMW temp_out=%.1fK quality=%.2f energy=%.1eJ",
             transferred_mw,
             secondary.temperature_out,
             secondary.steam_quality,
+            secondary.energy_j,
         )

src/reactor_sim/turbine.py

@@ -6,6 +6,7 @@ from dataclasses import dataclass
 import logging
 
 from . import constants
+from .thermal import saturation_temperature
 from .state import CoolantLoopState, TurbineState
 
 LOGGER = logging.getLogger(__name__)
@@ -50,10 +51,13 @@ class Turbine:
         throttle = min(constants.TURBINE_THROTTLE_MAX, max(constants.TURBINE_THROTTLE_MIN, self.throttle))
         throttle_eff = 1.0 - constants.TURBINE_THROTTLE_EFFICIENCY_DROP * (constants.TURBINE_THROTTLE_MAX - throttle)
 
-        enthalpy = 2_700.0 + loop.steam_quality * 600.0
+        sat_temp = saturation_temperature(max(0.05, loop.pressure))
+        superheat = max(0.0, loop.temperature_out - sat_temp)
+        enthalpy = (constants.STEAM_LATENT_HEAT / 1_000.0) + (superheat * constants.COOLANT_HEAT_CAPACITY / 1_000.0)
         mass_flow = effective_mass_flow * 0.6 * throttle
-        computed_power = (enthalpy * mass_flow / 1_000.0) / 1_000.0
+        computed_power = (enthalpy * mass_flow) / 1_000.0  # MW from enthalpy flow
         available_power = steam_power_mw if steam_power_mw > 0 else computed_power
+        available_power = min(available_power, computed_power)
         backpressure_loss = 1.0 - _backpressure_penalty(loop)
         shaft_power_mw = available_power * self.mechanical_efficiency * throttle_eff * backpressure_loss
         electrical = shaft_power_mw * self.generator_efficiency

tests/test_simulation.py

@@ -266,3 +266,31 @@ def test_auto_control_resets_shutdown_and_moves_rods():
     assert reactor.shutdown is False
     assert reactor.control.manual_control is False
     assert reactor.control.rod_fraction < 0.95
+
+
+def test_full_power_reaches_steam_and_turbine_output():
+    """Integration: cold start -> pumps/gens on -> ramp to ~3 GW -> steam -> turbines online."""
+    reactor = Reactor.default()
+    state = reactor.initial_state()
+    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=False,
+        ),
+    )
+    for i in range(600):
+        cmd = None
+        if i == 200:
+            cmd = ReactorCommand(secondary_pumps={2: False})
+        if i == 300:
+            cmd = ReactorCommand(secondary_pumps={2: True})
+        if i == 400:
+            cmd = ReactorCommand(turbine_on=True, turbine_units={1: True, 2: True, 3: True})
+        reactor.step(state, dt=1.0, command=cmd)
+
+    assert state.secondary_loop.steam_quality > 0.02
+    assert state.total_electrical_output() > 50.0

tests/test_thermal.py

@@ -30,8 +30,9 @@ def test_secondary_heats_to_saturation_before_boiling():
 def test_secondary_generates_steam_when_energy_exceeds_sensible_heat():
     solver = ThermalSolver()
     loop = _secondary_loop(temp_in=330.0, flow=180.0, pressure=0.5)
+    loop.inventory_kg *= 0.1  # reduce mass to let boil-up happen quickly
     sat_temp = saturation_temperature(loop.pressure)
-    solver.step_secondary(loop, transferred_mw=120.0, dt=1.0)
+    solver.step_secondary(loop, transferred_mw=120.0, dt=100.0)
-    assert loop.temperature_out == pytest.approx(sat_temp, rel=0.02)
+    assert loop.temperature_out == pytest.approx(sat_temp, rel=0.05)
     assert loop.steam_quality > 0.0
     assert loop.steam_quality < 1.0