formatting Battery.py, added 2022-12-23 TMS data

Battery.py

@@ -10,8 +10,9 @@ import pde
 import cantera as ct
 from scipy import optimize
 
+
 class CombustionChamber:
-    def __init__ (self, mdot, ct_object, burned_state, hA=700):
+    def __init__(self, mdot, ct_object, burned_state, hA=700):
         self.mdot = mdot                        # kg/s
         self.gas = ct_object                    # gas object
         self.eq_state = burned_state            # HP equilibrium state
@@ -27,21 +28,24 @@ class CombustionChamber:
         self.Twall1 = 1100 + 273.15
         self.Area = 6.7 * 16.7
 
-    def update_mdot (self, mdot_new):
-        if mdot_new : self.mdot = mdot_new
+    def update_mdot(self, mdot_new):
+        if mdot_new:
+            self.mdot = mdot_new
 
-    def update_Twall (self, Twall0=None, Twall1=None):
-        if Twall0: self.Twall0 = Twall0
-        if Twall1: self.Twall1 = Twall1
+    def update_Twall(self, Twall0=None, Twall1=None):
+        if Twall0:
+            self.Twall0 = Twall0
+        if Twall1:
+            self.Twall1 = Twall1
 
-    def energy_balance_equation (self, Tout):
+    def energy_balance_equation(self, Tout):
         self.gas.TP = Tout, None
         h1 = self.gas.enthalpy_mass
         q1, q2 = self.heat(Tout)
 
         return (self.mdot * (self.h0 - h1) - q1 - q2)
 
-    def solve (self, ):
+    def solve(self, ):
         """ Iteratively solve for outlet temperature that balance with heat loss to walls """
 
         meanTwall = (self.Twall0 + self.Twall1) / 2
@@ -55,30 +59,34 @@ class CombustionChamber:
 
         return f_found.root
 
-    def heat (self, Tout=None):
+    def heat(self, Tout=None):
         ''' Heat(W) to walls '''
-        if Tout is None: Tout = self.T1
+        if Tout is None:
+            Tout = self.T1
         Tgas = (self.T0 + Tout) / 2
 
         return self.hA * (Tgas - self.Twall0), self.hA * (Tgas - self.Twall1)
 
+
 class CokeCharge:
 
-    def __init__ (self, t_charge, idx_oven):
+    def __init__(self, t_charge, idx_oven):
         self.t_charge = t_charge
         self.t_push = None
         self.idx_oven = idx_oven
         self.Q = 0
 
-    def bake (self, dQ):
+    def bake(self, dQ):
         self.Q += dQ
 
-    def end_baking (self, t):
+    def end_baking(self, t):
         self.t_push = t
 
+
 brick_thickness = 0.14                            # m,
 n_grid_brick = 16                                 # Number of Grid points inside
-wall_grid = pde.CartesianGrid([[0, brick_thickness]], n_grid_brick, periodic=False)
+wall_grid = pde.CartesianGrid(
+    [[0, brick_thickness]], n_grid_brick, periodic=False)
 wall_area = 6.7 * 16.7  # m^2 , Oven cross section area
 
 # op_grad2 = wall_grid.make_operator_no_bc('gradient_squared', backend='scipy')
@@ -89,10 +97,11 @@ wall_area = 6.7 * 16.7 # m^2 , Oven cross section area
 # op_info_grad = wall_grid._get_operator_info('gradient')
 # op_info_lap = wall_grid._get_operator_info('laplace')
 
+
 class CokeOvenBrickHeatEqn(pde.PDEBase):
     """Implementation of the Heat equation"""
 
-    def __init__ (self, bc="auto_periodic_neumann"):
+    def __init__(self, bc="auto_periodic_neumann"):
         self.bc = bc
         self.rho = 1900            # kg / m3
         self.kCoef0 = 0.93         # W / m / K
@@ -100,13 +109,13 @@ class CokeOvenBrickHeatEqn(pde.PDEBase):
         self.cpCoef0 = 837.2       # J / kg / K
         self.cpCoef1 = 251.2e-3    # J / kg / K2
 
-    def k (self, T):
+    def k(self, T):
         return T * self.kCoef1 + self.kCoef0
 
-    def cp (self, T):
+    def cp(self, T):
         return T * self.cpCoef1 + self.cpCoef0
 
-    def update_bc (self, gradT_chamber=None, T_oven=None):
+    def update_bc(self, gradT_chamber=None, T_oven=None):
         bc0, bc1 = self.bc
         if gradT_chamber:
             self.bc[0] = {"derivative": gradT_chamber}
@@ -134,7 +143,8 @@ class CokeOvenBrickHeatEqn(pde.PDEBase):
         k = self.kCoef1 * state + self.kCoef0
         cp = self.cpCoef1 * state + self.cpCoef0
 
-        state_grad_k_grad = self.kCoef1 * state_grad2 # state_grad.dot(state_grad)
+        state_grad_k_grad = self.kCoef1 * \
+            state_grad2  # state_grad.dot(state_grad)
 
         return (state_grad_k_grad + k * state_lap) / cp / self.rho
 
@@ -156,15 +166,17 @@ class CokeOvenBrickHeatEqn(pde.PDEBase):
         return pde_rhs
 '''
 
+
 class RefractoryWall:
-    def __init__ (self, T0):
+    def __init__(self, T0):
         self.T_oven = T0
         self.T_chamber = T0
         self.q_chamber = 0.
         self.T_internal = pde.ScalarField(wall_grid, T0)
-        self.eqn = CokeOvenBrickHeatEqn(bc=[{"derivative": 0}, {"value": self.T_oven}])
+        self.eqn = CokeOvenBrickHeatEqn(
+            bc=[{"derivative": 0}, {"value": self.T_oven}])
 
-    def update_bc (self, Q=None, T_oven=None):
+    def update_bc(self, Q=None, T_oven=None):
         # Q = - k(T) gradT
         # T_chamber = self.T_internal.get_boundary_values(axis=0, upper=False, bc=self.eqn.bc)
         k0 = self.eqn.k(self.T_chamber)
@@ -174,13 +186,14 @@ class RefractoryWall:
             gradT = None
         self.eqn.update_bc(gradT, T_oven)
 
-    def solve (self, dt):
+    def solve(self, dt):
         # solution = self.eqn.solve (eqn, bc)
         self.T_internal = self.eqn.solve(
             self.T_internal, t_range=dt, dt=30., tracker='consistency', backend="numpy")
-        self.T_chamber = self.T_internal.get_boundary_values(axis=0, upper=False, bc=self.eqn.bc)
+        self.T_chamber = self.T_internal.get_boundary_values(
+            axis=0, upper=False, bc=self.eqn.bc)
 
-    def heat_to_oven (self):
+    def heat_to_oven(self):
         """ NOT YET IMPLEMENTED """
         return 0.0
 
@@ -197,11 +210,12 @@ def Twall_model(x):
     '''
     return np.interp(x, Twall_table[0], Twall_table[1])
 
+
 class OvenChamber:
-    def __init__ (self):
+    def __init__(self):
         self.content = None
 
-    def get_charge_temperature (self, t):
+    def get_charge_temperature(self, t):
         """ Return temperature of coal charge content at oven wall """
         if self.content:
             elapsed_time = t - self.content.t_charge
@@ -209,19 +223,21 @@ class OvenChamber:
             elapsed_time = 0.
         return Twall_model(elapsed_time)
 
-    def bake (self, q):
+    def bake(self, q):
         """ Add transferred heat to coal charge content """
         if self.content:
             self.content.bake(q)
 
-    def charge (self, coal_charge):
+    def charge(self, coal_charge):
         """ Update content with fresh coal is charged """
         self.content = coal_charge
 
+
 def wall_solve_wrapper(t_range, wall):
     wall.solve(t_range)
     return wall.T_internal, wall.T_chamber
 
+
 class Battery:
 
     def load_state(self):
@@ -237,22 +253,25 @@ class Battery:
         with open('oven.state', 'rb') as coke_state_file:
             self.processing = pickle.load(coke_state_file)
 
-    def __init__ (self, name, size, heat_program, charge_program, burned_gas_state, hv, init_from_file=False):
+    def __init__(self, name, size, heat_program, charge_program, burned_gas_state, hv, init_from_file=False):
         self.name = name                        # Battery name
         self.size = size                        # Size of battery, number of ovens
         self.heat_program = heat_program        # Heat program or schedule object
         self.charge_program = charge_program    # Charge program of schedule object
         self.t = 0                              # Battery time
         self.t_last = 0                         # Time of last Push/Charge
-        self.processing = []                    # List of Coke charges under processing(drying)
-        self.product = []                       # List of Coke charges done(completed)
+        # List of Coke charges under processing(drying)
+        self.processing = []
+        # List of Coke charges done(completed)
+        self.product = []
         self.gas = ct.Solution('gri30.xml')
         self.gas.TPX = burned_gas_state         # Burned gas T, P, X
         T0, P0, X0 = self.gas.TPX
         self.T0 = T0
         self.P0 = P0
         self.X0 = X0
-        self.sequence_idx = 0                   # Integer, 0 ~ (size-1), progress index for oven sequence array
+        # Integer, 0 ~ (size-1), progress index for oven sequence array
+        self.sequence_idx = 0
 
         self.wall_t_history = []
         self.gas_t_history = []
@@ -290,7 +309,8 @@ class Battery:
 
         # For 1~4 Coke Ovens with n+5 P/C sequence
         start_indices = [1, 3, 5, 2, 4]
-        self.oven_idx_order = np.concatenate([np.array(range(i0 - 1, self.size, 5)) for i0 in start_indices])
+        self.oven_idx_order = np.concatenate(
+            [np.array(range(i0 - 1, self.size, 5)) for i0 in start_indices])
 
         if init_from_file:
             print("Initializaton from file")
@@ -323,7 +343,8 @@ class Battery:
             normal_period = self.charge_program.period(-1)  # 감산 전 장입 간격 (주기)
             dt = normal_period / period_over_dt            # Simulation Time Step
 
-            self.t = - normal_period * self.size * n_cycle  # 정상상태 생성 모사 시간 = 장입 간격 * 총 장입 횟수
+            self.t = - normal_period * self.size * \
+                n_cycle  # 정상상태 생성 모사 시간 = 장입 간격 * 총 장입 횟수
             self.t_last = self.t                            # 마지막 장입을 정상상태 시뮬레이션 시작 시각으로 설정
 
             # initialization time loop
@@ -332,16 +353,16 @@ class Battery:
                 """ Fill battety with normal charge rate """
                 self.update(dt)     # Time adavancement
 
-    def mdot (self, t):
+    def mdot(self, t):
         return self.mdot0 * self.heat_program.f(t) / self.normal_heat
 
-    def next_oven (self):
+    def next_oven(self):
         ''' Index of the oven to which apply push and charge '''
         next_oven_id = self.oven_idx_order[self.sequence_idx % self.size]
         self.sequence_idx += 1
         return next_oven_id
 
-    def bake (self, dt):
+    def bake(self, dt):
         # update combustion chamber equilibrium temperature
         # Tad = 연료 조성과 공연비로 결정
         # m_dot = 연료 발열량과 공급열량 공연비로 결정
@@ -363,14 +384,19 @@ class Battery:
 
             chmbr.update_mdot(self.mdot(self.t)/self.size)
             chmbr.update_Twall(
-                Twall0=(wall_lower.T_chamber if wall_lower else wall_upper.T_chamber), 
-                Twall1=(wall_upper.T_chamber if wall_upper else wall_lower.T_chamber), 
+                Twall0=(
+                    wall_lower.T_chamber if wall_lower else wall_upper.T_chamber),
+                Twall1=(
+                    wall_upper.T_chamber if wall_upper else wall_lower.T_chamber),
             )
-            print(f"t={self.t:6.2} : {chmbr.Twall0} K | Chamber {i_chamber} | {chmbr.Twall1} K ")
+            print(
+                f"t={self.t:6.2} : {chmbr.Twall0} K | Chamber {i_chamber} | {chmbr.Twall1} K ")
             chmbr.solve()
             Q1, Q2 = chmbr.heat()  # W (J/s)
-            if wall_lower: wall_lower.update_bc(Q=Q1)
-            if wall_upper: wall_upper.update_bc(Q=Q2)
+            if wall_lower:
+                wall_lower.update_bc(Q=Q1)
+            if wall_upper:
+                wall_upper.update_bc(Q=Q2)
 
         # Loop all ovens
         # update oven wall temperatures using coke charge age
@@ -383,7 +409,8 @@ class Battery:
             wall_upper.update_bc(T_oven=T_oven)
 
         with Pool(12) as pool:
-            wall_sln = pool.starmap(wall_solve_wrapper, [((dt*60*60), w) for w in self.walls_0+self.walls_1])
+            wall_sln = pool.starmap(wall_solve_wrapper, [(
+                (dt*60*60), w) for w in self.walls_0+self.walls_1])
 
             # wall_lower.solve(dt * 60 * 60) # convert hours to seconds
             # wall_upper.solve(dt * 60 * 60) # convert hours to seconds
@@ -405,30 +432,30 @@ class Battery:
 
         # integrate heat to oven # 오븐 벽면 온도 우선 시간 함수로
 
-    def push_and_charge (self, coke_charge):
+    def push_and_charge(self, coke_charge):
         if len(self.processing) >= self.size:
             self.push(coke_charge.t_charge)
         self.charge(coke_charge)
 
-    def push (self, t):
+    def push(self, t):
         """ Push complete coke out of oven """
         coke = self.processing.pop(0)
         coke.end_baking(t)
         self.product.append(coke)
 
-    def charge (self, coke_charge):
+    def charge(self, coke_charge):
         self.ovens[coke_charge.idx_oven].charge(coke_charge)
         self.processing.append(coke_charge)
 
-    def dQ (self, dt):
+    def dQ(self, dt):
         return self.heat_program.dQ(self.t, self.t+dt)
 
-    def is_pc_time (self, dt):
+    def is_pc_time(self, dt):
         ''' Whether P/C should be done in this time step '''
         period = self.charge_program.period(self.t)
         return self.t + dt >= period + self.t_last
 
-    def update (self, dt):
+    def update(self, dt):
         # dQ = self.heat_program.dQ(self.t, self.t+dt) # t, t+dt 사이 공급하는 열량, array 로 대체 필요
 
         # t 에서 t+dt 까지 탄화실 가열
@@ -437,11 +464,12 @@ class Battery:
         period = self.charge_program.period(self.t)  # 현재 장입 시간 간격
 
         # 마지막 장입탄 장입 시각
-        latest_coke_charge = self.processing[-1].t_charge if len(self.processing) > 0 else self.t_last
+        latest_coke_charge = self.processing[-1].t_charge if len(
+            self.processing) > 0 else self.t_last
 
         # t_last + period 가 t, t + dt 사이에 들어오는 것 검사
         # t + dt 가 다음 추출/장입 시각 이후일 때 => 이번 time step 에 추출/장입을 실행해야함
-        if self.t + dt >= period + self.t_last :
+        if self.t + dt >= period + self.t_last:
             # 마지막 장입 시각 + 장입 시간 간격 이 이번 time step 에 포함됨
             # 일정한 간격으로 장입 진행 중, 마지막 장입 시간 += 장입 간격
             if self.t < self.t_last + period:
@@ -461,14 +489,16 @@ class Battery:
                   f"{self.t + dt - latest_coke_charge:7.3} since last P/C. ",
                   f"period = {self.charge_program.period(self.t):7.3}",)
 
-        
         # 시뮬레이션 시간 업데이트
         self.t += dt
 
-        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.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)]))
 
-def coke_oven_exhaust_stoichiometry (phi=1.0, return_unburned=False):
+
+def coke_oven_exhaust_stoichiometry(phi=1.0, return_unburned=False):
 
     # Define the oxidizer composition, here air with 21 mol-% O2 and 79 mol-% N2
     air = "O2:1,N2:3.762"
@@ -479,7 +509,8 @@ def coke_oven_exhaust_stoichiometry (phi=1.0, return_unburned=False):
     mix.TP = 25+273.15, ct.one_atm
     mix.set_equivalence_ratio(phi=phi, fuel=coke_oven_fuel, oxidizer=air)
 
-    element_X = {ename: mix.elemental_mole_fraction(ename) for ename in mix.element_names}
+    element_X = {ename: mix.elemental_mole_fraction(
+        ename) for ename in mix.element_names}
 
     exhaust = ct.Solution('gri30.xml')
     exhaust.TPX = (25+273.15, ct.one_atm,
@@ -496,8 +527,9 @@ def coke_oven_exhaust_stoichiometry (phi=1.0, return_unburned=False):
     else:
         return exhaust.mole_fraction_dict(threshold=-1)
 
+
 class HeatSchedule:
-    def __init__ (self, xp, fp):
+    def __init__(self, xp, fp):
         self.xp = xp
         self.fp = fp
         self.f = lambda x: np.interp(x, self.xp, self.fp)
@@ -506,21 +538,23 @@ class HeatSchedule:
         x = np.linspace(t0, t1, 31)
         return np.trapz(self.f(x), x)
 
+
 class ChargeSchedule:
-    def __init__ (self, normal_load, service_start, service_time, service_load, aux_start, aux_time, aux_load):
+    def __init__(self, normal_load, service_start, service_time, service_load, aux_start, aux_time, aux_load):
         self.xp = np.array([service_start, service_start, service_start+service_time, service_start+service_time,
                             aux_start, aux_start, aux_start+aux_time, aux_start+aux_time, ])
         self.fp = np.array([normal_load, service_load, service_load, normal_load,
                             normal_load, aux_load, aux_load, normal_load])
         self.f = lambda x: np.interp(x, self.xp, self.fp)
 
-    def to_charge (self, t0, t1):
+    def to_charge(self, t0, t1):
         self.f(t0)
         return np.trapz(self.f(x), x)
 
-    def period (self, t):
+    def period(self, t):
         return 24 / self.f(t)
 
+
 if __name__ == "__main__":
 
     # Define the oxidizer composition, here air with 21 mol-% O2 and 79 mol-% N2
@@ -553,14 +587,16 @@ if __name__ == "__main__":
     gas_in_state = gas.TPX
 
     # Time(Hours) - GJ/rev
-    sample_program = np.array(open('sample_heat_221128_3A.txt').read().split(), dtype=np.double).reshape((-1,2))
+    sample_program = np.array(open(
+        'sample_heat_221128_3A-Plan2.txt').read().split(), dtype=np.double).reshape((-1, 2))
 
     heating_plan = HeatSchedule(*sample_program.T)
-    charging_plan = ChargeSchedule( 82, 9, 13, 1e-12, 9+13+18, 4, 1e-12 )
+    charging_plan = ChargeSchedule(82, 9, 12, 1e-12, 9+13+18, 4, 1e-12)
     n_doors = 66
 
-    load_state = False
-    bat3A = Battery("3A", n_doors, heating_plan, charging_plan, gas_in_state, hv, init_from_file=load_state)
+    load_state = True
+    bat3A = Battery("3A", n_doors, heating_plan, charging_plan,
+                    gas_in_state, hv, init_from_file=load_state)
 
     if not load_state:
         with open('gas.history', 'wb') as gas_history_file:
@@ -576,7 +612,7 @@ if __name__ == "__main__":
             pickle.dump(bat3A.processing, wall_history_file)
 
     dt = 5. * 1./60.  # 5 min
-    for it in range (int(60/dt)): # 시뮬레이션 시간 도메인 = 60시간
+    for it in range(int(60/dt)):  # 시뮬레이션 시간 도메인 = 60시간
         bat3A.update(dt)
 
     with open('gas.history2', 'wb') as gas_history_file: