diff --git a/TODO.md b/TODO.md index b9daa0f..b082c0b 100644 --- a/TODO.md +++ b/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. - - Adjust HX/pressure handling to use stored energy (saturation clamp and pressure rise) and validate steam formation with both pumps at ~3 GW. - - Update turbine power mapping to consume steam enthalpy/quality and align protection trips with real steam presence. - - 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 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. 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; 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 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. diff --git a/src/reactor_sim/reactor.py b/src/reactor_sim/reactor.py index 1bb545c..74ec5e7 100644 --- a/src/reactor_sim/reactor.py +++ b/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 diff --git a/src/reactor_sim/state.py b/src/reactor_sim/state.py index 0625b31..2c03fd3 100644 --- a/src/reactor_sim/state.py +++ b/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"]), - secondary_loop=CoolantLoopState(**data["secondary_loop"]), + primary_loop=CoolantLoopState(**_with_energy(data["primary_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 diff --git a/src/reactor_sim/thermal.py b/src/reactor_sim/thermal.py index 3f3417c..78b6496 100644 --- a/src/reactor_sim/thermal.py +++ b/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.""" - 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 + """Update secondary loop using a stored-energy steam-drum balance.""" cp = constants.COOLANT_HEAT_CAPACITY - sat_temp = saturation_temperature(max(0.05, secondary.pressure)) - energy_j = max(0.0, transferred_mw) * constants.MEGAWATT * dt + mass = max(1e-6, secondary.inventory_kg) + if secondary.energy_j <= 0.0: + secondary.energy_j = mass * cp * secondary.average_temperature() - # Energy needed to heat incoming feed to saturation. - sensible_j = max(0.0, sat_temp - temp_in) * mass_flow * cp * dt - if energy_j <= sensible_j: - delta_t = temperature_rise(transferred_mw, mass_flow) - secondary.temperature_out = temp_in + delta_t - secondary.steam_quality = 0.0 + # Add transferred heat; if no heat, bleed toward ambient. + if transferred_mw > 0.0: + secondary.energy_j += transferred_mw * constants.MEGAWATT * dt else: - energy_left = energy_j - sensible_j - steam_mass = energy_left / constants.STEAM_LATENT_HEAT - produced_fraction = steam_mass / max(1e-6, mass_flow * dt) - secondary.temperature_out = sat_temp - secondary.steam_quality = min(1.0, max(0.0, produced_fraction)) + bleed = 0.01 * (secondary.temperature_out - constants.ENVIRONMENT_TEMPERATURE) + secondary.energy_j = max( + mass * cp * constants.ENVIRONMENT_TEMPERATURE, secondary.energy_j - max(0.0, bleed) * mass * cp * dt + ) + + 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, ) diff --git a/src/reactor_sim/turbine.py b/src/reactor_sim/turbine.py index 1844cd7..494a542 100644 --- a/src/reactor_sim/turbine.py +++ b/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 diff --git a/tests/test_simulation.py b/tests/test_simulation.py index da9fa7f..577bf75 100644 --- a/tests/test_simulation.py +++ b/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 diff --git a/tests/test_thermal.py b/tests/test_thermal.py index b606ddc..1589ed5 100644 --- a/tests/test_thermal.py +++ b/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) - assert loop.temperature_out == pytest.approx(sat_temp, rel=0.02) + solver.step_secondary(loop, transferred_mw=120.0, dt=100.0) + assert loop.temperature_out == pytest.approx(sat_temp, rel=0.05) assert loop.steam_quality > 0.0 assert loop.steam_quality < 1.0