diff --git a/Battery.py b/Battery.py index c6bc905..db2103f 100644 --- a/Battery.py +++ b/Battery.py @@ -1,8 +1,54 @@ from functools import reduce import numpy as np - +import pde import cantera as ct +from scipy import optimize + +class CombustionChamber: + 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 + self.gas.TPX = burned_state # Set equilibrium state + T0, P0, X0 = self.gas.TPX + self.T0 = T0 # K, adiabatic flame temperature + self.P0 = P0 # Pa, pressure + self.X0 = X0 # Composition in mole fractions, Fuel + Air + self.h0 = self.gas.enthalpy_mass # inlet enthalpy + self.hA = hA # HTC x Area + self.Twall0 = 1100 + 273.15 + 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_Twall (self, Twall0=None, Twall1=None): + if Twall0: self.Twall0 = Twall0 + if Twall1: self.Twall1 = Twall1 + + 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, ): + """ Iteratively solve for outlet temperature that balance with heat loss to walls """ + f_found = optimize.root_scalar(self.energy_balance_equation, + bracket=[max(self.Twall0, self.Twall1), self.T0]) + self.T1 = f_found.root + + return f_found.root + + def heat (self, Tout=None): + ''' Heat(W) to walls ''' + 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: @@ -18,10 +64,111 @@ class CokeCharge: def end_baking (self, t): self.t_push = t +class CokeOvenBrickHeatEqn(pde.PDEBase): + """Implementation of the Heat equation""" + + def __init__ (self, bc="auto_periodic_neumann"): + self.bc = bc + self.rho = 1900 # kg / m3 + self.kCoef0 = 0.93 # W / m / K + self.kCoef1 = 0.698e-3 # W / m / K2 + self.cpCoef0 = 837.2 # J / kg / K + self.cpCoef1 = 251.2e-3 # J / kg / K2 + + def k (self, T): + return T * self.kCoef1 + self.kCoef0 + + def cp (self, T): + return T * self.cpCoef1 + self.cpCoef0 + + def update_bc (self, gradT_chamber=None, T_oven=None): + bc0, bc1 = self.bc + if gradT_chamber: + self.bc[0] = {"derivative": gradT_chamber} + if T_oven: + self.bc[1] = {"value": T_oven} + + def evolution_rate(self, state, t=0): + """implement the python version of the evolution equation""" + state_lap = state.laplace(bc=self.bc) + state_grad = state.gradient(bc=self.bc) + + k = self.kCoef1 * state + self.kCoef0 + cp = self.cpCoef1 * state + self.cpCoef0 + + k_grad = self.kCoef1 * state_grad + + return (state_grad.dot(k_grad) + k * state_lap) / self.rho / cp + +brick_thickness = 0.14 # m, +n_grid_brick = 32 # Number of Grid points inside +wall_grid = pde.CartesianGrid([[0, brick_thickness]], n_grid_brick, periodic=False) +wall_area = 6.7 * 16.7 # m^2 , Oven cross section area + +class RefractoryWall: + 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}]) + + 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) + if Q: + gradT = Q / wall_area / k0 + else: + gradT = None + self.eqn.update_bc(gradT, T_oven) + + def solve (self, dt): + # solution = self.eqn.solve (eqn, bc) + self.T_internal = self.eqn.solve(self.T_internal, t_range=dt, dt=1., tracker='consistency') + self.T_chamber = self.T_internal.get_boundary_values(axis=0, upper=False, bc=self.eqn.bc) + + def heat_to_oven (self): + """ NOT YET IMPLEMENTED """ + return 0.0 + + +Twall_table = np.loadtxt('./CokeOvenWallTemperature.csv', delimiter=',').T +Twall_table[0] *= ((66/80) / (100/80)) +Twall_table[1] += 273.15 + + +def Twall_model(x): + ''' + Coke oven wall temperature vs time after charging + Temperature (K) vs Elapsed time (hour) + ''' + return np.interp(x, Twall_table[0], Twall_table[1]) + +class OvenChamber: + def __init__ (self): + self.content = None + + 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 + else: + elapsed_time = 0. + return Twall_model(elapsed_time) + + def bake (self, q): + """ Add transferred heat to coal charge content """ + if self.content: + self.content.bake(q) + + def charge (self, coal_charge): + """ Update content with fresh coal is charged """ + self.content = coal_charge class Battery: - def __init__ (self, name, size, heat_program, charge_program, T_combustion_0): + def __init__ (self, name, size, heat_program, charge_program, burned_gas_state, hv): self.name = name # Battery name self.size = size # Size of battery, number of ovens self.heat_program = heat_program # Heat program or schedule object @@ -30,17 +177,42 @@ class Battery: 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) - self.T0 = T_combustion_0 # Burned gas temperature + 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 + self.hv = hv # Base unit heat J/kg + self.normal_heat = self.heat_program.f(-1) # GJ / rev + + Q0 = self.normal_heat * 1e9 * 3 / 3600 # GJ/rev => J/s (W) + mdot0 = Q0 / hv # (J/s) / (J/kg) => kg/s + self.mdot0 = mdot0 # kg / s + # chambers[0] - walls_0[0] - ovens[0] - walls_1[0] - chambers[1] - walls_0[1] - ... # ... walls_1[i-1] - chambers[i] - walls_0[i] - ovens[i] - walls_1[i] - chambers[i+1] - walls_0[i+1] - ... # ... walls_1[size-2] - chambers[size-1] - walls_0[size-1] - ovens[size-1] - walls_1[size-1] - chambers[size] - self.chambers = np.zeros(self.size+1) - self.ovens = np.zeros(self.size) - self.walls_0 = np.zeros(self.size) - self.walls_1 = np.zeros(self.size) + self.chambers = [ + CombustionChamber(self.mdot0, self.gas, + (self.T0, self.P0, self.X0), hA=700) + for ichamber in range(self.size+1) + ] + self.ovens = [ + OvenChamber() + for ioven in range(self.size) + ] + self.walls_0 = [ + RefractoryWall(Twall_model(0)) + for ioven in range(self.size) + ] + self.walls_1 = [ + RefractoryWall(Twall_model(0)) + for ioven in range(self.size) + ] # For 1~4 Coke Ovens with n+5 P/C sequence start_indices = [1, 3, 5, 2, 4] @@ -60,6 +232,9 @@ class Battery: """ Fill battety with normal charge rate """ self.update(dt) # Time adavancement + def mdot (self, t): + return self.mdot0 * self.heat_program.f(t) / self.normal_heat + 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] @@ -86,13 +261,14 @@ class Battery: except IndexError: wall_upper = None - chmbr.update_mdot() + chmbr.update_mdot(self.mdot(self.t)) chmbr.update_Twall( wall_lower.T_chamber if wall_lower else wall_upper.T_chamber, wall_upper.T_chamber if wall_upper else wall_lower.T_chamber, ) + print(f"t={self.t} - C{i_chamber} with {chmbr.Twall0} K and {chmbr.Twall1} K ") chmbr.solve() - Q1, Q2 = chmbr.heat() + 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) @@ -101,13 +277,13 @@ class Battery: # solve heat equations of all walls # bake charge in oven for i_oven, (oven, wall_lower, wall_upper) in enumerate(zip(self.ovens, self.walls_0, self.walls_1)): - oven.get_charge_temperature() + T_oven = oven.get_charge_temperature(self.t) - wall_lower.update_bc(Toven=Q1) - wall_upper.update_bc(Toven=Q2) + wall_lower.update_bc(T_oven=T_oven) + wall_upper.update_bc(T_oven=T_oven) - wall_lower.solve(dt) - wall_upper.solve(dt) + wall_lower.solve(dt * 60 * 60) + wall_upper.solve(dt * 60 * 60) ql = wall_lower.heat_to_oven() qu = wall_upper.heat_to_oven() @@ -120,10 +296,11 @@ class Battery: # integrate heat to oven # 오븐 벽면 온도 우선 시간 함수로 - dQ = self.dQ(dt) # array, dQ pairs of all oven taking from both walls + '''dQ = self.dQ(dt) # array, dQ pairs of all oven taking from both walls for cc in self.processing: cc.bake(dQ) # bake.(dQ[cc.idx_oven]) + ''' def push_and_charge (self, coke_charge): if len(self.processing) >= self.size: @@ -137,14 +314,15 @@ class Battery: self.product.append(coke) def charge (self, coke_charge): - self.ovens[coke_charge.i_oven].charge(coke_charge) + self.ovens[coke_charge.idx_oven].charge(coke_charge) self.processing.append(coke_charge) def dQ (self, dt): return self.heat_program.dQ(self.t, self.t+dt) def is_pc_time (self, dt): - ''' P/C should be done in this time step ''' + ''' 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): @@ -161,9 +339,9 @@ class Battery: # t_last + period 가 t, t + dt 사이에 들어오는 것 검사 # t + dt 가 다음 추출/장입 시각 이후일 때 => 이번 time step 에 추출/장입을 실행해야함 if self.t + dt >= period + self.t_last : - print(f"Push timing within [ {self.t} , {self.t + dt} ].", - f"{self.t + dt - latest_coke_charge} since last P/C. ", - f"P/C period = {self.charge_program.period(self.t)}",) + print(f"P/C within [ {self.t:7.3} , {self.t + dt:7.3} ].", + f"{self.t + dt - latest_coke_charge:7.3} since last P/C. ", + f"period = {self.charge_program.period(self.t):7.3}",) # 마지막 장입 시각 + 장입 시간 간격 이 이번 time step 에 포함됨 # 일정한 간격으로 장입 진행 중, 마지막 장입 시간 += 장입 간격 @@ -181,4 +359,132 @@ class Battery: self.push_and_charge(fresh_coal) # 시뮬레이션 시간 업데이트 - self.t += dt \ No newline at end of file + self.t += dt + +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:0.21,N2:0.79" + coke_oven_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" + + mix = ct.Solution('gri30.xml') + + 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} + + exhaust = ct.Solution('gri30.xml') + exhaust.TPX = (25+273.15, ct.one_atm, + { + "CO2": element_X['C'], + "H2O": element_X['H']/2, + "O2": (element_X['O'] - 2*element_X['C'] - element_X['H']/2)/2, + "N2": element_X['N']/2, + } + ) + + if return_unburned: + return mix.mole_fraction_dict(threshold=-1), exhaust.mole_fraction_dict(threshold=-1) + else: + return exhaust.mole_fraction_dict(threshold=-1) + +class HeatSchedule: + def __init__ (self, xp, fp): + self.xp = xp + self.fp = fp + self.f = lambda x: np.interp(x, self.xp, self.fp) + + def dQ(self, t0, t1): + 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): + 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): + self.f(t0) + return np.trapz(self.f(x), x) + + 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 + air = "O2:0.21,N2:0.79" + coke_oven_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" + + f_found = optimize.root_scalar(lambda x: coke_oven_exhaust_stoichiometry(x)["O2"] - 0.045, + bracket=[1e-300, 1]) + + # equivalence ratio for O2 4.5 % in exhaust gas (stoichiometric) + phi_O2_045 = f_found.root + + # unburned and burned gas compositions for O2 4.5 % in exhaust gas (stoichiometric) + Xu, Xb = coke_oven_exhaust_stoichiometry(phi_O2_045, return_unburned=True) + + gas = ct.Solution('gri30.xml') + + # Heating value of unburned premixed gas + gas.TPX = 25 + 273.15, ct.one_atm, Xu + hu = gas.enthalpy_mass + gas.TPX = None, None, Xb + hb = gas.enthalpy_mass + hv = hu - hb + print(f'{hu*1e-6} - {hb*1e-6} = {hv*1e-6} MJ/kg') + + # burned premixed gas state (chemical equilibrium with HP constraint) + gas.TP = 600+273.15, ct.one_atm + gas.set_equivalence_ratio(phi_O2_045, fuel=coke_oven_fuel, oxidizer=air) + gas.equilibrate('HP') + gas_in_state = gas.TPX + + # Time(Hours) - GJ/rev + sample_program = np.array('''\ +-3 81 +0 81 +0 69.61621622 +3.5 69.61621622 +3.5 58.83157895 +7 58.83157895 +7 48.6 +10.5 48.6 +10.5 42.039 +14.5 42.039 +14.5 38.88 +24.5 38.88 +24.5 42.039 +27.75 42.039 +27.75 48.6 +30.35 48.6 +30.35 58.83157895 +32.43 58.83157895 +32.43 69.61621622 +34.094 69.61621622 +34.094 81 +36.46688136 81 +36.46688136 72.9 +39.46688136 72.9 +39.46688136 67.23 +42.46688136 67.23 +42.46688136 62.37 +47.46688136 62.37 +47.46688136 67.23 +50.46688136 67.23 +50.46688136 72.9 +53.46688136 72.9 +53.46688136 81 +56.46688136 81 +'''.split(), dtype=np.double).reshape((-1,2)) + + heating_plan = HeatSchedule(*sample_program.T) + charging_plan = ChargeSchedule( 81, 9, 9, 1e-12, 24+13, 3, 1e-12 ) + n_doors = 66 + + bat3A = Battery("3A", n_doors, heating_plan, charging_plan, gas_in_state, hv) \ No newline at end of file