Add feedwater controller and turbine governor

src/reactor_sim/reactor.py

@@ -61,6 +61,9 @@ class Reactor:
         self.turbine_active = any(self.turbine_unit_active)
         if not self.generators:
             self.generators = [DieselGenerator() for _ in range(2)]
+        # Balance-of-plant controls
+        self.feedwater_valve = 0.5
+        self._last_steam_out_kg_s = 0.0
         if not self.primary_pump_units or len(self.primary_pump_units) != 2:
             self.primary_pump_units = [True, True]
         if not self.secondary_pump_units or len(self.secondary_pump_units) != 2:
@@ -218,9 +221,6 @@ class Reactor:
         self._update_loop_inventory(
             state.primary_loop, constants.PRIMARY_LOOP_VOLUME_M3, constants.PRIMARY_INVENTORY_TARGET, dt
         )
-        self._update_loop_inventory(
-            state.secondary_loop, constants.SECONDARY_LOOP_VOLUME_M3, constants.SECONDARY_INVENTORY_TARGET, dt
-        )
 
         pump_demand = overrides.get(
             "coolant_demand",
@@ -424,9 +424,7 @@ class Reactor:
         if not self.control.manual_control and not self.shutdown:
             self.control.safety_backoff(state.core.subcooling_margin, state.core.dnb_margin, dt)
         self._apply_secondary_boiloff(state, dt)
-        self._update_loop_inventory(
-            state.secondary_loop, constants.SECONDARY_LOOP_VOLUME_M3, constants.SECONDARY_INVENTORY_TARGET, dt
-        )
+        self._update_secondary_level(state, dt)
 
         steam_draw = self._step_turbine_bank(state, transferred, dt)
         if steam_draw > 0.0:
@@ -533,9 +531,17 @@ class Reactor:
             if idx in active_indices:
                 # Simple throttle map: reduce throttle when electrical demand is low, open as demand rises.
                 demand = turbine_state.load_demand_mw
-                throttle = 0.4 if demand <= 0 else min(1.0, 0.4 + demand / max(1e-6, turbine.rated_output_mw))
-                turbine.throttle = throttle
+                desired = 0.4 if demand <= 0 else min(1.0, 0.4 + demand / max(1e-6, turbine.rated_output_mw))
+                # Governor: nudge throttle toward desired based on electrical error.
+                error = (demand - turbine_state.electrical_output_mw) / max(1.0, turbine.rated_output_mw)
+                turbine.throttle = max(0.3, min(1.0, turbine.throttle + (desired - turbine.throttle) * 0.5 + 0.2 * error * dt))
                 turbine.step(state.secondary_loop, turbine_state, steam_power_mw=power_per_unit, dt=dt)
+                if turbine_state.electrical_output_mw > 1.05 * turbine.rated_output_mw:
+                    LOGGER.critical("Turbine %d overspeed/overload, tripping unit", idx + 1)
+                    self._spin_down_turbine(turbine_state, dt, turbine.spool_time)
+                    turbine_state.status = "TRIPPED"
+                    self.turbine_unit_active[idx] = False
+                    continue
                 if power_per_unit <= 0.0 and turbine_state.electrical_output_mw < 0.1:
                     turbine_state.status = "OFF"
                 elif turbine_state.electrical_output_mw < max(0.1 * turbine.rated_output_mw, 1.0):
@@ -613,6 +619,39 @@ class Reactor:
         loop.inventory_kg = max(0.0, loop.inventory_kg + correction * nominal_mass * dt)
         loop.level = min(1.2, max(0.0, loop.inventory_kg / nominal_mass))
 
+    def _update_secondary_level(self, state: PlantState, dt: float) -> None:
+        """Steam drum level controller with shrink/swell and feedwater valve."""
+        loop = state.secondary_loop
+        nominal_mass = self._nominal_inventory(constants.SECONDARY_LOOP_VOLUME_M3)
+        if nominal_mass <= 0.0:
+            loop.level = 0.0
+            return
+        if loop.inventory_kg <= 0.0:
+            loop.inventory_kg = nominal_mass * constants.SECONDARY_INVENTORY_TARGET
+        current_level = loop.inventory_kg / nominal_mass
+        steam_out = loop.mass_flow_rate * max(0.0, loop.steam_quality)
+        # Shrink/swell: apparent level drops when steam draw surges.
+        swell = -0.02 * (steam_out - self._last_steam_out_kg_s) / max(1.0, nominal_mass)
+        sensed_level = current_level + swell
+        # PI-ish valve adjustment toward target level.
+        error = constants.SECONDARY_INVENTORY_TARGET - sensed_level
+        valve_delta = 0.3 * error * dt
+        self.feedwater_valve = max(0.0, min(1.0, self.feedwater_valve + valve_delta))
+        # Feedwater adds mass and energy at inlet temperature.
+        steam_factor = min(1.0, max(0.1, steam_out / max(1.0, nominal_mass * 0.1)))
+        feed_rate = (
+            self.feedwater_valve
+            * nominal_mass
+            * 0.002
+            * steam_factor
+        )  # up to ~0.2% of nominal mass per second, scaled by steam draw
+        added_mass = feed_rate * dt
+        loop.inventory_kg = max(0.0, loop.inventory_kg + added_mass)
+        cp = constants.COOLANT_HEAT_CAPACITY
+        loop.energy_j += added_mass * cp * loop.temperature_in
+        loop.level = min(1.2, max(0.0, loop.inventory_kg / nominal_mass))
+        self._last_steam_out_kg_s = steam_out
+
     def _inventory_flow_scale(self, loop: CoolantLoopState) -> float:
         if loop.level <= constants.LOW_LEVEL_FLOW_FLOOR:
             return 0.0