feat: implement 1D regenerator lumped thermal mass model and reversing cycle logic

Battery.py

@@ -140,6 +140,62 @@ class CombustionChamber:
         return self.hA * (Tgas - self.Twall0), self.hA * (Tgas - self.Twall1)
 
 
+class Regenerator:
+    """Models a single regenerator chamber in the coke oven battery.
+
+    Uses a 1D lumped thermal mass model to track temperature dynamics during
+    heat storage (charging) and release (discharging) cycles.
+    """
+
+    def __init__(self, idx, T_init=1016.9, M_reg=1e9, Cp_reg=1000.0, eff=0.8):
+        """Initializes the Regenerator.
+
+        Args:
+            idx (int): Index identifier of the regenerator.
+            T_init (float, optional): Initial brick temperature (K). Default is 1016.9 K
+                to replicate the 600°C boundary condition.
+            M_reg (float, optional): Mass of checker bricks (kg). Default is 1e9
+                (option A: infinite thermal mass to suppress dynamics for verification).
+            Cp_reg (float, optional): Specific heat of checker bricks (J/kg/K). Default 1000.
+            eff (float, optional): Heat exchange effectiveness (0 to 1). Default 0.8.
+        """
+        self.idx = idx
+        self.T = T_init
+        self.M_reg = M_reg
+        self.Cp_reg = Cp_reg
+        self.eff = eff
+        self.Cp_gas = 1000.0  # J/kg/K, average heat capacity of gas mixture
+
+    def heat_exchange(self, dt_sec, mdot, T_gas_in, is_discharging):
+        """Calculates gas outlet temperature and updates regenerator brick temperature.
+
+        Args:
+            dt_sec (float): Time step in seconds.
+            mdot (float): Mass flow rate of air/flue gas (kg/s).
+            T_gas_in (float): Temperature of incoming gas (K).
+            is_discharging (bool): True if air/gas is being preheated (discharging),
+                False if hot flue gas is reheating the bricks (charging).
+
+        Returns:
+            float: Temperature of the outgoing gas (K).
+        """
+        if mdot <= 1e-8:
+            return T_gas_in
+
+        if is_discharging:
+            # Heat release cycle: incoming cold air at T_gas_in is preheated to T_gas_out
+            T_gas_out = T_gas_in + self.eff * (self.T - T_gas_in)
+            Q_dot = mdot * self.Cp_gas * (T_gas_out - T_gas_in)  # W (J/s)
+            self.T -= (Q_dot * dt_sec) / (self.M_reg * self.Cp_reg)
+            return T_gas_out
+        else:
+            # Heat storage cycle: incoming hot flue gas at T_gas_in heats the bricks
+            T_gas_out = T_gas_in - self.eff * (T_gas_in - self.T)
+            Q_dot = mdot * self.Cp_gas * (T_gas_in - T_gas_out)  # W (J/s)
+            self.T += (Q_dot * dt_sec) / (self.M_reg * self.Cp_reg)
+            return T_gas_out
+
+
 class CokeCharge:
     """Represents a single coal/coke charge inside an oven.
 
@@ -627,6 +683,28 @@ class Battery:
             for ioven in range(self.size)
         ]
 
+        # Regenerators and fuel valves initialization (size + 2 elements)
+        # Option A is default: M_reg=1e9 to replicate the legacy 600°C boundary condition.
+        self.regenerators = [
+            Regenerator(ireg, T_init=1016.9, M_reg=1e9, Cp_reg=1000.0, eff=0.8)
+            for ireg in range(self.size + 2)
+        ]
+        self.fuel_valves = np.zeros(self.size + 2)
+        self.reversing_period = 20.0 / 60.0  # 20 minutes in hours
+        self.control_active = False
+
+        # Cantera parameters for dynamic preheating and equilibrium calculations
+        self.fuel = "H2:6.42, O2:0.39, N2:47.28, CH4:1.79, CO:24.25, CO2:19.72, C2H4:0.13, C2H6:0.04"
+        self.oxidizer = "O2:0.21,N2:0.79"
+        try:
+            f_found = optimize.root_scalar(
+                lambda x: coke_oven_exhaust_stoichiometry(x)["O2"] - 0.045,
+                bracket=[1e-300, 1]
+            )
+            self.phi = f_found.root
+        except Exception:
+            self.phi = 0.814238515
+
         # For 1~4 Coke Ovens with n+5 P/C sequence
         start_indices = [1, 3, 5, 2, 4]
         self.oven_idx_order = np.concatenate(
@@ -647,9 +725,14 @@ class Battery:
             for chmbr, T1 in zip(self.chambers, latest_chamber[1]):
                 chmbr.T1 = T1
 
-            # Recover Wall State
+            # Recover Wall State with copy() to break reference linkage with history lists
             for wl, wu, wallT in zip(self.walls_0, self.walls_1, latest_wall[1]):
-                wl.T_chamber, wl.T_internal.data, wl.T_oven, wu.T_oven, wu.T_internal.data, wu.T_chamber = wallT
+                wl.T_chamber = wallT[0]
+                wl.T_internal.data = wallT[1].copy()
+                wl.T_oven = wallT[2]
+                wu.T_oven = wallT[3]
+                wu.T_internal.data = wallT[4].copy()
+                wu.T_chamber = wallT[5]
 
             # Recover Oven State
             for coal in self.processing:
@@ -684,6 +767,95 @@ class Battery:
         """
         return self.mdot0 * self.heat_program.f(t) / self.normal_heat
 
+    def is_cycle_A(self, t):
+        """Checks if the battery is currently in Cycle A (odd-numbered regenerators discharging).
+
+        Args:
+            t (float): Simulation time (hours).
+
+        Returns:
+            bool: True if in Cycle A, False if in Cycle B.
+        """
+        cycle_num = int(np.floor(t / self.reversing_period))
+        return (cycle_num % 2) == 0
+
+    def get_chamber_inlets(self, t):
+        """Computes mass flow rate and maps regenerators for each combustion chamber.
+
+        Args:
+            t (float): Current simulation time (hours).
+
+        Returns:
+            tuple: (chamber_mdots, inlet_reg_indices, outlet_reg_indices)
+                - chamber_mdots (np.ndarray): Fuel-air flow rate for each chamber (kg/s).
+                - inlet_reg_indices (list): Regenerator index supplying to each chamber.
+                - outlet_reg_indices (list): Regenerator index receiving exhaust from each chamber.
+        """
+        total_mdot = self.mdot(t)
+        size = self.size
+
+        # If external control is inactive, dynamically reset valves to legacy distribution
+        if not getattr(self, "control_active", False):
+            self.fuel_valves = np.zeros(size + 2)
+            if self.is_cycle_A(t):
+                # Cycle A: Odd regenerators (0-indexed even j) discharge
+                self.fuel_valves[0] = total_mdot / size
+                self.fuel_valves[2:size+1:2] = 2.0 * total_mdot / size
+            else:
+                # Cycle B: Even regenerators (0-indexed odd j) discharge
+                self.fuel_valves[1:size+1:2] = 2.0 * total_mdot / size
+                self.fuel_valves[size+1] = total_mdot / size
+
+        # Distribute valve flows to individual combustion chambers
+        chamber_mdots = np.zeros(size + 1)
+        inlet_reg_indices = [-1] * (size + 1)
+        outlet_reg_indices = [-1] * (size + 1)
+
+        is_a = self.is_cycle_A(t)
+
+        if is_a:
+            # Cycle A: Even regenerators (j = 0, 2, ..., size) are INLETS
+            # Odd regenerators (j = 1, 3, ..., size+1) are OUTLETS
+            for j in range(0, size + 1, 2):
+                val = self.fuel_valves[j]
+                if j == 0:
+                    chamber_mdots[0] += val
+                    inlet_reg_indices[0] = 0
+                else:
+                    chamber_mdots[j-1] += val / 2.0
+                    chamber_mdots[j] += val / 2.0
+                    inlet_reg_indices[j-1] = j
+                    inlet_reg_indices[j] = j
+
+            for j in range(1, size + 2, 2):
+                if j == size + 1:
+                    outlet_reg_indices[size] = size + 1
+                else:
+                    outlet_reg_indices[j-1] = j
+                    outlet_reg_indices[j] = j
+        else:
+            # Cycle B: Odd regenerators (j = 1, 3, ..., size+1) are INLETS
+            # Even regenerators (j = 0, 2, ..., size) are OUTLETS
+            for j in range(1, size + 2, 2):
+                val = self.fuel_valves[j]
+                if j == size + 1:
+                    chamber_mdots[size] += val
+                    inlet_reg_indices[size] = size + 1
+                else:
+                    chamber_mdots[j-1] += val / 2.0
+                    chamber_mdots[j] += val / 2.0
+                    inlet_reg_indices[j-1] = j
+                    inlet_reg_indices[j] = j
+
+            for j in range(0, size + 1, 2):
+                if j == 0:
+                    outlet_reg_indices[0] = 0
+                else:
+                    outlet_reg_indices[j-1] = j
+                    outlet_reg_indices[j] = j
+
+        return chamber_mdots, inlet_reg_indices, outlet_reg_indices
+
     def next_oven(self):
         """Returns the index of the next oven to be pushed and charged.
 
@@ -700,14 +872,36 @@ class Battery:
         Args:
             dt (float): Simulation time step (hours).
         """
-        # update combustion chamber equilibrium temperature
-        # Tad = 연료 조성과 공연비로 결정
-        # m_dot = 연료 발열량과 공급열량 공연비로 결정
-        # m(h1 - h0) = hA(Tgas - Twall) => solve with initial T0 = Tad
+        dt_sec = dt * 3600.0
+        size = self.size
+        is_a = self.is_cycle_A(self.t)
 
-        # Loop all combustion chambers
-        # update chamber wall temperatures and mass flow rates
-        # solve for equilibrium heat to walls
+        # Get dynamic mass flows and regenerator mappings for this step
+        chamber_mdots, inlet_reg_indices, outlet_reg_indices = self.get_chamber_inlets(self.t)
+
+        # 1. Aggregate flow rates for the discharging (inlet) regenerators
+        inlet_mdots = np.zeros(size + 2)
+        for i_chamber, mdot in enumerate(chamber_mdots):
+            inlet_reg_idx = inlet_reg_indices[i_chamber]
+            inlet_mdots[inlet_reg_idx] += mdot
+
+        # 2. Perform heat exchange for discharging regenerators once to find preheat temperatures
+        reg_preheat_temps = np.zeros(size + 2)
+        active_inlets = range(0, size + 1, 2) if is_a else range(1, size + 2, 2)
+        for j in active_inlets:
+            reg_preheat_temps[j] = self.regenerators[j].heat_exchange(
+                dt_sec, inlet_mdots[j], 298.15, is_discharging=True
+            )
+
+        # Cache for Cantera HP equilibrium calculation to save CPU overhead.
+        # Maps inlet_reg_idx -> (T0, burned_state, h0)
+        cantera_cache = {}
+
+        # Arrays to aggregate exhaust properties for the charging (outlet) regenerators
+        outlet_mdots = np.zeros(size + 2)
+        outlet_T_weighted = np.zeros(size + 2)
+
+        # 3. Loop all combustion chambers to solve thermal equations
         for i_chamber, chmbr in enumerate(self.chambers):
             if i_chamber > 0:
                 wall_lower = self.walls_1[i_chamber-1]
@@ -719,7 +913,29 @@ class Battery:
             else:
                 wall_upper = None
 
-            chmbr.update_mdot(self.mdot(self.t)/self.size)
+            mdot = chamber_mdots[i_chamber]
+            inlet_reg_idx = inlet_reg_indices[i_chamber]
+            T_pre = reg_preheat_temps[inlet_reg_idx]
+
+            # Compute new combustion state (HP equilibrium) based on T_pre
+            if inlet_reg_idx not in cantera_cache:
+                self.gas.TP = T_pre, self.P0
+                self.gas.set_equivalence_ratio(self.phi, self.fuel, self.oxidizer)
+                self.gas.equilibrate('HP')
+                burned_state = self.gas.TPX
+                h0 = self.gas.enthalpy_mass
+                T0 = self.gas.T
+                cantera_cache[inlet_reg_idx] = (T0, burned_state, h0)
+
+            T0, burned_state, h0 = cantera_cache[inlet_reg_idx]
+
+            # Update chamber state variables
+            chmbr.update_mdot(mdot)
+            chmbr.gas.TPX = burned_state
+            chmbr.T0 = T0
+            chmbr.h0 = h0
+            chmbr.T1 = T0  # reset outlet guess to Tad
+
             chmbr.update_Twall(
                 Twall0=(
                     wall_lower.T_chamber if wall_lower else wall_upper.T_chamber),
@@ -735,6 +951,21 @@ class Battery:
             if wall_upper:
                 wall_upper.update_bc(Q=Q2)
 
+            # Accumulate exhaust properties for charging (outlet) regenerator
+            outlet_reg_idx = outlet_reg_indices[i_chamber]
+            outlet_mdots[outlet_reg_idx] += mdot
+            outlet_T_weighted[outlet_reg_idx] += mdot * chmbr.T1
+
+        # 4. Perform heat exchange for charging regenerators once using mass-weighted exhaust temps
+        active_outlets = range(1, size + 2, 2) if is_a else range(0, size + 1, 2)
+        for j in active_outlets:
+            mdot_tot = outlet_mdots[j]
+            if mdot_tot > 1e-8:
+                T_exh_avg = outlet_T_weighted[j] / mdot_tot
+                self.regenerators[j].heat_exchange(
+                    dt_sec, mdot_tot, T_exh_avg, is_discharging=False
+                )
+
         # Loop all ovens
         # update oven wall temperatures using coke charge age
         # solve heat equations of all walls
@@ -867,8 +1098,8 @@ class Battery:
 
         self.gas_t_history.append(
             (self.t, [chmbr.T1 for chmbr in self.chambers]))
-        self.wall_t_history.append((self.t, [(wl.T_chamber, wl.T_internal.data, wl.T_oven, wu.T_oven,
-                                   wu.T_internal.data, wu.T_chamber) for wl, wu in zip(self.walls_0, self.walls_1)]))
+        self.wall_t_history.append((self.t, [(wl.T_chamber, wl.T_internal.data.copy(), wl.T_oven, wu.T_oven,
+                                   wu.T_internal.data.copy(), wu.T_chamber) for wl, wu in zip(self.walls_0, self.walls_1)]))
 
 
 def coke_oven_exhaust_stoichiometry(phi=1.0, return_unburned=False):