Skip to content

Updated PEM cell modeling schema #75

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
elenya-grant wants to merge 4 commits into
base: feature/BERT
Choose a base branch
from
Open

Updated PEM cell modeling schema #75

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
b7a2ed8
reformatted
elenya-grant Jan 26, 2024
a7e3778
updated PEM schema and removed LCOH tests temporarily
elenya-grant Jan 30, 2024
43808e1
updated alkaline modeling schema and alkaline Urev0 calc
elenya-grant Jan 30, 2024
aae1054
update lcoh glue_code test regression values
bayc Feb 1, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
56 changes: 36 additions & 20 deletions electrolyzer/PEM_cell.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -42,7 +42,17 @@ class PEMCell(FromDictMixin):
cell_area: float
turndown_ratio: float
max_current_density: float

# everything below is new
p_anode: float # bar
p_cathode: float # bar
alpha_a: float
alpha_c: float
i_0_a: float
i_0_c: float
e_m: float
R_ohmic_elec: float
f_1: float
f_2: float
# If we rework this class to be even more generic, we can have these be specified
# as configuration params

Expand All @@ -64,8 +74,8 @@ def calc_open_circuit_voltage(self, temperature):
"""Calculates open circuit voltage using the Nernst equation."""
T_K = convert_temperature([temperature], "C", "K")[0]
E_rev_0 = self.calc_reversible_voltage()
p_anode = P_ATMO # (Pa) assumed atmo
p_cathode = 30 * P_STD # (Pa) 30 bars
p_anode = self.p_anode * P_STD # (Pa)
p_cathode = self.p_cathode * P_STD # (Pa)

# noqa: E501
# Arden Buck equation T=C, https://www.omnicalculator.com/chemistry/vapour-pressure-of-water#vapor-pressure-formulas # noqa
Expand All @@ -83,6 +93,9 @@ def calc_open_circuit_voltage(self, temperature):
p_o2 = p_anode - p_h2O_sat

# General Nernst equation, 10.1016/j.ijhydene.2017.03.046
# E_cell = E_rev_0 + ((R * T_K) / (self.n * F)) * (
# np.log((p_h2 / p_h2O_sat) * np.sqrt(p_o2 / p_h2O_sat))
# )
E_cell = E_rev_0 + ((R * T_K) / (self.n * F)) * (
np.log((p_h2 / P_ATMO) * np.sqrt(p_o2 / P_ATMO))
)
Expand All @@ -105,20 +118,24 @@ def calc_activation_overpotential(self, i, temperature):
T_cathode = T_K # TODO: updated with realistic anode temperature?

# anode charge transfer coefficient TODO: is this a realistic value?
alpha_a = 2
# alpha_a = 2

# cathode charge transfer coefficient TODO: is this a realistic value?
alpha_c = 0.5
# alpha_c = 0.5

# anode exchange current density TODO: update to be f(T)?
i_0_a = 2e-7
# i_0_a = 2e-7

# cathode exchange current density TODO: update to be f(T)?
i_0_c = 2e-3
# i_0_c = 2e-3

# derived from Butler-Volmer eqs
V_act_a = ((R * T_anode) / (alpha_a * F)) * np.arcsinh(i / (2 * i_0_a))
V_act_c = ((R * T_cathode) / (alpha_c * F)) * np.arcsinh(i / (2 * i_0_c))
V_act_a = ((R * T_anode) / (self.alpha_a * F)) * np.arcsinh(
i / (2 * self.i_0_a)
)
V_act_c = ((R * T_cathode) / (self.alpha_c * F)) * np.arcsinh(
i / (2 * self.i_0_c)
)

# alternate equations for Activation overpotential
# Option 2: Dakota: I believe this may be more accurate, found more
Expand All @@ -144,24 +161,24 @@ def calc_ohmic_overpotential(self, i, temperature):
# pulled from https://www.sciencedirect.com/science/article/pii/S0360319917309278?via%3Dihub # noqa
# TODO: pulled from empirical data, is there a better eq?
lambda_nafion = ((-2.89556 + (0.016 * T_K)) + 1.625) / 0.1875
t_nafion = 0.02 # (cm) confirmed that membrane thickness is <0.02.
# t_nafion = 0.02 # (cm) confirmed that membrane thickness is <0.02.

# TODO: confirm with Nel, is there a better eq?
sigma_nafion = ((0.005139 * lambda_nafion) - 0.00326) * np.exp(
1268 * ((1 / 303) - (1 / T_K))
)
R_ohmic_ionic = t_nafion / sigma_nafion
R_ohmic_ionic = self.e_m / sigma_nafion

# TODO: confirm realistic value with Nel https://www.sciencedirect.com/science/article/pii/S0378775315001901 # noqa
R_ohmic_elec = 50e-3
# R_ohmic_elec = 50e-3

# Alternate R_ohmic_elec from https://www.sciencedirect.com/science/article/pii/S0360319918309017 # noqa
# rho = (ohm*m) material resistivity
# l_path = (m) length of electron path
# A_path = (m2) cross-sectional area of conductor path
# R_ohmic_elec = ((rho*l_path)/A_path)

V_ohmic = i * (R_ohmic_elec + R_ohmic_ionic)
V_ohmic = i * (self.R_ohmic_elec + R_ohmic_ionic)

return V_ohmic

Expand Down Expand Up @@ -225,25 +242,24 @@ def calc_faradaic_efficiency(self, T_C, I):
return :: eta_F [-]: Faraday's efficiency
Reference: https://res.mdpi.com/d_attachment/energies/energies-13-04792/article_deploy/energies-13-04792-v2.pdf
""" # noqa
f_1 = 250 # (mA2/cm4)
f_2 = 0.996
# f_1 = 250 # (mA2/cm4)
# f_2 = 0.996
I *= 1000

eta_F = (
((I / self.cell_area) ** 2) / (f_1 + ((I / self.cell_area) ** 2))
) * f_2
((I / self.cell_area) ** 2) / (self.f_1 + ((I / self.cell_area) ** 2))
) * self.f_2

return eta_F

def calc_mass_flow_rate(self, T_C, Idc, dryer_loss=6.5):
def calc_mass_flow_rate(self, T_C, Idc):
"""
Idc [A]: stack current
dryer_loss [%]: loss of drying H2
T_C [C]: cell temperature (currently unused)
return :: mfr [kg/s]: mass flow rate
"""
eta_F = self.calc_faradaic_efficiency(T_C, Idc)
mfr = eta_F * Idc * self.M / (self.n * F) * (1 - dryer_loss / 100.0) # [g/s]
mfr = eta_F * Idc * self.M / (self.n * F) # [g/s]
# mfr = mfr / 1000. * 3600. # [kg/h]
mfr = mfr / 1e3 # [kg/s]
return mfr
115 changes: 76 additions & 39 deletions electrolyzer/alkaline_cell.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -129,11 +129,15 @@ class AlkalineCell(FromDictMixin):
turndown_ratio: float
max_current_density: float

f_1: float
f_2: float

# cell_area: float = field(init=False)
cell_area: float = field(init=False)

# Electrode parameters #
####################
A_electrode: float = field(init=False) # [cm^2]
# A_electrode: float = field(init=False) # [cm^2]
e_a: float = field(init=False) # [cm] anode thickness
e_c: float = field(init=False) # [cm] cathode thickness
d_am: float = field(init=False) # [cm] Anode-membrane gap
Expand All @@ -151,7 +155,7 @@ class AlkalineCell(FromDictMixin):

# Membrane parameters #
####################
A_membrane: float = field(init=False) # [cm^2]
# A_membrane: float = field(init=False) # [cm^2]
e_m: float = field(init=False) # [cm] membrane thickness

# THIS ONE IS PRIMARLY BASED ON
Expand All @@ -166,6 +170,7 @@ class AlkalineCell(FromDictMixin):
M_K: float = 39.0983 # molecular weight of Potassium [g/mol]
lhv: float = 33.33 # lower heating value of H2 [kWh/kg]
hhv: float = 39.41 # higher heating value of H2 [kWh/kg]
gibbs: float = 237.24e3 # Gibbs Energy of global reaction (J/mol)

def __attrs_post_init__(self) -> None:
# Cell parameters #
Expand All @@ -174,11 +179,17 @@ def __attrs_post_init__(self) -> None:

# Electrode parameters #
########################
self.A_electrode = self.electrode["A_electrode"]
self.e_a = self.electrode["e_a"]
self.e_c = self.electrode["e_c"]
self.d_am = self.electrode["d_am"]
self.d_cm = self.electrode["d_cm"]
# self.A_electrode = self.electrode["A_electrode"]
# self.e_a = self.electrode["e_a"]
# self.e_c = self.electrode["e_c"]
self.e_a = self.electrode["e_e"]
self.e_c = self.electrode["e_e"]

# self.d_am = self.electrode["d_am"]
# self.d_cm = self.electrode["d_cm"]

self.d_am = self.electrode["d_em"]
self.d_cm = self.electrode["d_em"]
self.d_ac = self.electrode["d_ac"]

# Electrolyte parameters #
Expand All @@ -187,7 +198,7 @@ def __attrs_post_init__(self) -> None:

# Membrane parameters #
#######################
self.A_membrane = self.membrane["A_membrane"]
# self.A_membrane = self.membrane["A_membrane"]
self.e_m = self.membrane["e_m"]

# calcluate molarity and molality of KOH solution
Expand All @@ -214,7 +225,8 @@ def calculate_bubble_rate_coverage(self, T_C, I):
T_k = convert_temperature([T_C], "C", "K")[0]
J_lim = 30 # [A/cm^2] [Vogt,Balzer 2005]
T_amb = T_k = convert_temperature([25], "C", "K")[0]
j = I / self.A_electrode # [A/cm^2] "nominal current density"
# j = I / self.A_electrode # [A/cm^2] "nominal current density"
j = I / self.cell_area # [A/cm^2] "nominal current density"

# Eqn 19 of [Gambou, Guilbert,et al 2022]
Pv_H20 = np.exp(
Expand Down Expand Up @@ -247,7 +259,8 @@ def calc_current_density(self, T_C, I):
# actual current density reflecting impact of bubble rate coverage
theta_epsilon = self.calculate_bubble_rate_coverage(T_C, I)
theta = theta_epsilon[0]
A_electrode_eff = self.A_electrode * (1 - theta) # [cm^2]
# A_electrode_eff = self.A_electrode * (1 - theta) # [cm^2]
A_electrode_eff = self.cell_area * (1 - theta) # [cm^2]
current_dens = I / A_electrode_eff # [A/cm^2]
return current_dens

Expand Down Expand Up @@ -369,7 +382,8 @@ def calc_Urev(self, T_C, P):
# Eqn 16
Pv_KOH = np.exp((2.302 * a) + (b * np.log(Pv_H20)))

Urev0 = self.calc_open_circuit_voltage(T_C)
Urev0 = self.calc_reversible_voltage()
# Urev0 = self.calc_open_circuit_voltage(T_C)

# Eqn 14
U_rev = Urev0 + ((R * T_K) / (self.z * F)) * np.log(
Expand Down Expand Up @@ -450,15 +464,25 @@ def calc_electrode_resistance(self, T_C):
# Eqn 21 - effective resistance of electrode
rho_nickle_eff = rho_nickle_0 / ((1 - epsilon_Ni) ** 1.5)
# Eqn 20 - resistivity of anode
# Ra = (
# rho_nickle_eff
# * (self.e_a / self.A_electrode)
# * (1 + (temp_coeff * (T_C - tref)))
# )
Ra = (
rho_nickle_eff
* (self.e_a / self.A_electrode)
* (self.e_a / self.cell_area)
* (1 + (temp_coeff * (T_C - tref)))
)
# Eqn 20 - resistivity of cathode
# Rc = (
# rho_nickle_eff
# * (self.e_c / self.A_electrode)
# * (1 + (temp_coeff * (T_C - tref)))
# )
Rc = (
rho_nickle_eff
* (self.e_c / self.A_electrode)
* (self.e_c / self.cell_area)
* (1 + (temp_coeff * (T_C - tref)))
)

Expand Down Expand Up @@ -505,8 +529,11 @@ def calc_electrolyte_resistance(self, T_C, I):
# R_ele_bf: Bubble-free electrolyte resistance
# Eqn 32 of [Gambou, Guilbert,et al 2022] and Eqn 19 of [Henou, Agbossou, 2014]

# R_ele_bf = (100 / sigma_bf) * (
# (self.d_am / self.A_electrode) + (self.d_cm / self.A_electrode)
# )
R_ele_bf = (100 / sigma_bf) * (
(self.d_am / self.A_electrode) + (self.d_cm / self.A_electrode)
(self.d_am / self.cell_area) + (self.d_cm / self.cell_area)
)
# Eqn 2 of [Brauns,2021] says R_eg=(1/(sigma_koh))*((dcs+das)/A_el)
# where A_el is the metal area, not the entire area (which has holes)
Expand Down Expand Up @@ -540,8 +567,11 @@ def calc_membrane_resistance(self, T_C):
# [Gambou, Guilbert,et al 2022]
# S_mem=54.48 # membrane surface area in cm^2

# Rmem = (0.06 + 80 * np.exp(T_C / 50)) / (
# 10000 * self.A_membrane
# ) # Equation 36 - Ohms
Rmem = (0.06 + 80 * np.exp(T_C / 50)) / (
10000 * self.A_membrane
10000 * self.cell_area
) # Equation 36 - Ohms
# ^ Equation 21 of [Henou, Agbossou, 2014]
# output: Rmem=0.23 ohm*cm^2
Expand Down Expand Up @@ -595,29 +625,36 @@ def calc_ohmic_overpotential(self, T_C, I):
V_ohm = I * R_tot # [V/cell]
return V_ohm

def calc_open_circuit_voltage(self, T_C):
def calc_reversible_voltage(self):
"""
I [A]: current
T_C [C]: temperature
return :: E_rev0 [V/cell]: open-circuit voltage

TODO: Are we correcting for temperature twice? U_rev0 should be just 1.229 and
never change (possibly?)

Reference: [Gambou, Guilbert,et al 2022]: Eqn 14
Calculates reversible cell potential at standard state.
"""
# General Nerst Equation
# Eqn 14 of [Gambou, Guilbert,et al 2022]
T_K = convert_temperature([T_C], "C", "K")[0]
E_rev0 = (
1.5184
- (1.5421 * (10 ** (-3)) * T_K)
+ (9.523 * (10 ** (-5)) * T_K * np.log(T_K))
+ (9.84 * (10 ** (-8)) * (T_K**2))
)
# OR should this just be 1.229?
# E_rev_fake = 1.229
return E_rev0
return self.gibbs / (self.z * F)

# def calc_open_circuit_voltage(self, T_C):
# """
# I [A]: current
# T_C [C]: temperature
# return :: E_rev0 [V/cell]: open-circuit voltage

# TODO: Are we correcting for temperature twice? U_rev0 should be just 1.229 and
# never change (possibly?)

# Reference: [Gambou, Guilbert,et al 2022]: Eqn 14
# """
# # General Nerst Equation
# # Eqn 14 of [Gambou, Guilbert,et al 2022]
# T_K = convert_temperature([T_C], "C", "K")[0]
# E_rev0 = (
# 1.5184
# - (1.5421 * (10 ** (-3)) * T_K)
# + (9.523 * (10 ** (-5)) * T_K * np.log(T_K))
# + (9.84 * (10 ** (-8)) * (T_K**2))
# )
# # OR should this just be 1.229?
# # E_rev_fake = 1.229

# return E_rev0

def calc_faradaic_efficiency(self, T_C, I):
"""
Expand All @@ -627,8 +664,8 @@ def calc_faradaic_efficiency(self, T_C, I):
Reference: [Oystein Ulleberg, 2003] Eqn 9, Table 3
"""
# f1 and f2 values from Table 3
f1 = 250 # [mA^2/cm^4]
f2 = 0.96 # [-]
# f1 = 250 # [mA^2/cm^4]
# f2 = 0.96 # [-]
j = self.calc_current_density(T_C, I) # [A/cm^2]
j *= 1000 # [mA/cm^2]

Expand All @@ -640,7 +677,7 @@ def calc_faradaic_efficiency(self, T_C, I):
# f_2=[0.99,0.985,0.98,0.96] #[0-1]

# Eqn 9 from [Oystein Ulleberg, 2003]
eta_F = f2 * (j**2) / (f1 + j**2)
eta_F = self.f_2 * (j**2) / (self.f_1 + j**2)
return eta_F

def calc_mass_flow_rate(self, T_C, I):
Expand Down
Loading