Prototype Done
This commit is contained in:
Yeongdo Park
parent
0f7e6bec97
commit
6b7251ffcd
1 changed files with 327 additions and 21 deletions
348
Battery.py
348
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
|
||||
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)
|
||||
Loading…
Add table
Reference in a new issue