Add rod banks, xenon kinetics, and document feature set

src/reactor_sim/neutronics.py

@@ -7,7 +7,7 @@ import logging
 
 from . import constants
 from .fuel import fuel_reactivity_penalty
-from .state import CoreState
+from .state import CoreState, clamp
 
 LOGGER = logging.getLogger(__name__)
 
@@ -28,15 +28,29 @@ class NeutronDynamics:
     beta_effective: float = 0.0065
     delayed_neutron_fraction: float = 0.0008
     external_source_coupling: float = 1e-6
-    shutdown_bias: float = -0.012
+    shutdown_bias: float = -0.014
+    iodine_yield: float = 1e-6  # inventory units per MW*s
+    iodine_decay_const: float = 1.0 / 66000.0  # ~18h
+    xenon_decay_const: float = 1.0 / 33000.0  # ~9h
+    xenon_burnout_coeff: float = 1e-13  # per n/cm2
+    xenon_reactivity_coeff: float = 0.002
 
-    def reactivity(self, state: CoreState, control_fraction: float) -> float:
+    def reactivity(self, state: CoreState, control_fraction: float, rod_banks: list[float] | None = None) -> float:
+        if rod_banks:
+            weights = constants.CONTROL_ROD_BANK_WEIGHTS
+            worth = 0.0
+            total = sum(weights)
+            for w, pos in zip(weights, rod_banks):
+                worth += w * (1.0 - clamp(pos, 0.0, 0.95) / 0.95)
+            rod_term = constants.CONTROL_ROD_WORTH * worth / total
+        else:
+            rod_term = constants.CONTROL_ROD_WORTH * (1.0 - control_fraction)
         rho = (
             self.shutdown_bias +
-            constants.CONTROL_ROD_WORTH * (1.0 - control_fraction)
+            rod_term
             + temperature_feedback(state.fuel_temperature)
             - fuel_reactivity_penalty(state.burnup)
-            - xenon_poisoning(state.neutron_flux)
+            - self._xenon_penalty(state)
         )
         return rho
 
@@ -48,8 +62,15 @@ class NeutronDynamics:
         source_term = self.external_source_coupling * external_source_rate
         return ((rho - beta) / generation_time) * state.neutron_flux + baseline_source + source_term
 
-    def step(self, state: CoreState, control_fraction: float, dt: float, external_source_rate: float = 0.0) -> None:
-        rho = self.reactivity(state, control_fraction)
+    def step(
+        self,
+        state: CoreState,
+        control_fraction: float,
+        dt: float,
+        external_source_rate: float = 0.0,
+        rod_banks: list[float] | None = None,
+    ) -> None:
+        rho = self.reactivity(state, control_fraction, rod_banks)
         rho = min(rho, 0.02)
         shutdown = control_fraction >= 0.95
         if shutdown:
@@ -65,3 +86,15 @@ class NeutronDynamics:
             state.neutron_flux,
             d_flux,
         )
+
+    def update_poisons(self, state: CoreState, dt: float) -> None:
+        prod_I = max(0.0, state.power_output_mw) * self.iodine_yield
+        decay_I = state.iodine_inventory * self.iodine_decay_const
+        state.iodine_inventory = max(0.0, state.iodine_inventory + (prod_I - decay_I) * dt)
+        prod_Xe = decay_I
+        burn_Xe = state.neutron_flux * self.xenon_burnout_coeff
+        decay_Xe = state.xenon_inventory * self.xenon_decay_const
+        state.xenon_inventory = max(0.0, state.xenon_inventory + (prod_Xe - decay_Xe - burn_Xe) * dt)
+
+    def _xenon_penalty(self, state: CoreState) -> float:
+        return min(0.03, state.xenon_inventory * self.xenon_reactivity_coeff)