Esg/alkaline

.pre-commit-config.yaml

@@ -1,6 +1,6 @@
 repos:
 -   repo: https://github.com/timothycrosley/isort
-    rev: 5.10.1
+    rev: 5.12.0
     hooks:
         - id: isort
           name: isort

CITATION.cff

@@ -0,0 +1,14 @@
+cff-version: 0.0.0
+message: "If you use this software, please cite it as below."
+authors:
+- family-names: "Mona"
+  given-names: "Lisa"
+  orcid: "https://orcid.org/0000-0000-0000-0000"
+- family-names: "Bot"
+  given-names: "Hew"
+  orcid: "https://orcid.org/0000-0000-0000-0000"
+title: "My Research Software"
+version: 2.0.4
+doi: 10.5281/zenodo.1234
+date-released: 2017-12-18
+url: "https://github.com/github-linguist/linguist"

electrolyzer/cell.py → electrolyzer/PEM_cell.py

@@ -1,6 +1,7 @@
 """
 This module defines a Hydrogen Electrolyzer Cell
 """
+
 # TODOs
 # * refine calcCellVoltage(); compare with alkaline models
 # * refine convertACtoDC(); compare with empirical ESIF model
@@ -11,10 +12,10 @@
 from attrs import define
 from scipy.constants import R, physical_constants, convert_temperature
 
-from .type_dec import FromDictMixin
+from electrolyzer.type_dec import FromDictMixin
 
 
-def electrolyzer_model(X, a, b, c, d, e, f):
+def PEM_electrolyzer_model(X, a, b, c, d, e, f):
     """
     Given a power input (kW), temperature (C), and set of coefficients, returns
     current (A).  Coefficients can be determined using non-linear least squares
@@ -30,14 +31,17 @@ def electrolyzer_model(X, a, b, c, d, e, f):
 #############
 F, _, _ = physical_constants["Faraday constant"]  # Faraday's constant [C/mol]
 P_ATMO, _, _ = physical_constants["standard atmosphere"]  # Pa
+P_STD, _, _ = physical_constants["standard-state pressure"]  # Pa (1bar)
 
 
 @define
-class Cell(FromDictMixin):
+class PEMCell(FromDictMixin):
     # Chemical Params #
     ###################
 
     cell_area: float
+    turndown_ratio: float
+    max_current_density: float
 
     # If we rework this class to be even more generic, we can have these be specified
     # as configuration params
@@ -61,10 +65,11 @@ def calc_open_circuit_voltage(self, temperature):
         T_K = convert_temperature([temperature], "C", "K")[0]
         E_rev_0 = self.calc_reversible_voltage()
         p_anode = P_ATMO  # (Pa) assumed atmo
-        p_cathode = P_ATMO
+        p_cathode = 30 * P_STD  # (Pa) 30 bars
 
         # noqa: E501
         # Arden Buck equation T=C, https://www.omnicalculator.com/chemistry/vapour-pressure-of-water#vapor-pressure-formulas # noqa
+        # Reasonable at temperatures between 0-100C
         p_h2O_sat = (
             0.61121
             * np.exp(
@@ -73,12 +78,13 @@ def calc_open_circuit_voltage(self, temperature):
             )
         ) * 1e3  # (Pa)
 
-        # General Nernst equation
+        # Dalton's Law to find partial pressure of reactants at each electrode.
+        p_h2 = p_cathode - p_h2O_sat
+        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_anode - p_h2O_sat) / P_ATMO)
-                * np.sqrt((p_cathode - p_h2O_sat) / P_ATMO)
-            )
+            np.log((p_h2 / P_ATMO) * np.sqrt(p_o2 / P_ATMO))
         )
 
         return E_cell
@@ -138,7 +144,7 @@ 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.03  # (cm) TODO: confirm actual thickness?
+        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(
@@ -212,9 +218,10 @@ def calc_cell_voltage(self, I, temperature):
     # ------------------------------------------------------------
     # Post H2 production
     # ------------------------------------------------------------
-    def calc_faradaic_efficiency(self, I):
+    def calc_faradaic_efficiency(self, T_C, I):
         """
         I [A]: current
+        T_C [C]: cell temperature (currently unused)
         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
@@ -228,13 +235,14 @@ def calc_faradaic_efficiency(self, I):
 
         return eta_F
 
-    def calc_mass_flow_rate(self, Idc, dryer_loss=6.5):
+    def calc_mass_flow_rate(self, T_C, Idc, dryer_loss=6.5):
         """
         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(Idc)
+        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 = mfr / 1000. * 3600. # [kg/h]
         mfr = mfr / 1e3  # [kg/s]

electrolyzer/__init__.py

@@ -6,13 +6,10 @@
 
 # noqa
 
-from .cell import Cell, electrolyzer_model
 from .lcoh import LCOH
 from .stack import Stack
+from .PEM_cell import PEMCell, PEM_electrolyzer_model
 from .supervisor import Supervisor
 from .alkaline_cell import AlkalineCell, ael_electrolyzer_model
-from .alkaline_stack import AlkalineStack
 from .glue_code.run_lcoh import run_lcoh
-from .alkaline_supervisor import AlkalineSupervisor
-from .glue_code.run_alkaline import run_alkaline
 from .glue_code.run_electrolyzer import run_electrolyzer

electrolyzer/alkaline_cell.py

@@ -1,19 +1,20 @@
 """
 This module defines an Alkaline Hydrogen Electrolyzer Cell
 """
+
 import warnings
 
 import numpy as np
 
-# TODOs
+# TODO
 # * refine calcCellVoltage(); compare with alkaline models
 # * refine convertACtoDC(); compare with empirical ESIF model
 # * refine calcFaradaicEfficiency(); compare with other model
 # * add a separate script to show results
-from attrs import define
+from attrs import field, define
 from scipy.constants import R, physical_constants, convert_temperature
 
-from .type_dec import FromDictMixin
+from electrolyzer.type_dec import FromDictMixin
 
 
 warnings.filterwarnings("ignore")
@@ -118,29 +119,40 @@ def ael_electrolyzer_model(X, a, b, c, d, e, f):
 class AlkalineCell(FromDictMixin):
     # Cell parameters #
     ####################
+    model: str
 
-    # Stack parameters #
-    ####################
-    pressure_operating: float  # [bar]
-    n_cells: int  # number of cells
+    electrode: dict
+    electrolyte: dict
+    membrane: dict
+
+    pressure_operating: float
+    turndown_ratio: float
+    max_current_density: float
+
+    cell_area: float = field(init=False)
 
     # Electrode parameters #
     ####################
-    A_electrode: float  # [cm^2]
-    e_a: float  # [cm] anode thickness
-    e_c: float  # [cm] cathode thickness
-    d_am: float  # [cm] Anode-membrane gap
-    d_cm: float  # [cm] Cathode-membrane gap
-    d_ac: float  # [cm] Anode-Cathode gap
+    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
+    d_cm: float = field(init=False)  # [cm] Cathode-membrane gap
+    d_ac: float = field(init=False)  # [cm] Anode-Cathode gap
 
     # Electrolyte parameters #
     ####################
-    w_koh: float
+    w_koh: float = field(init=False)
+    electrolyte_concentration_percent: float = field(init=False)
+    M_KOH: float = field(init=False)
+    M_H2O: float = field(init=False)
+    m: float = field(init=False)
+    M: float = field(init=False)
 
     # Membrane parameters #
     ####################
-    A_membrane: float  # [cm^2]
-    e_m: float  # [cm] membrane thickness
+    A_membrane: float = field(init=False)  # [cm^2]
+    e_m: float = field(init=False)  # [cm] membrane thickness
 
     # THIS ONE IS PRIMARLY BASED ON
     # VALUES FROM [Henou, Agbossou, 2014]
@@ -156,6 +168,28 @@ class AlkalineCell(FromDictMixin):
     hhv: float = 39.41  # higher heating value of H2 [kWh/kg]
 
     def __attrs_post_init__(self) -> None:
+        # Cell parameters #
+        ###################
+        self.cell_area = self.electrode["A_electrode"]
+
+        # 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.d_ac = self.electrode["d_ac"]
+
+        # Electrolyte parameters #
+        ##########################
+        self.w_koh = self.electrolyte["w_koh"]
+
+        # Membrane parameters #
+        #######################
+        self.A_membrane = self.membrane["A_membrane"]
+        self.e_m = self.membrane["e_m"]
+
         # calcluate molarity and molality of KOH solution
         self.electrolyte_concentration_percent = self.w_koh / 100
         self.create_electrolyte()
@@ -237,7 +271,7 @@ def create_electrolyte(self):
         moles_of_solute = grams_of_solute / self.M_KOH  # [mols of solute / solution]
 
         # solvent is water
-        self.M_H20 = 2 * self.M_H + self.M_0  # [g/mol]
+        self.M_H2O = 2 * self.M_H + self.M_0  # [g/mol]
         grams_of_solvent = solution_weight_g * (
             1 - self.electrolyte_concentration_percent
         )
@@ -390,9 +424,9 @@ def calc_activation_overpotential(self, T_C, I):
         # Eqn 13 - Tafel slope for cathode
         bc = (R * T_anode) / (self.z * F * alpha_c)
         # Eqn 11 - anode activation energy
-        V_act_a = ba * np.log(ja / j0a)
+        V_act_a = ba * np.maximum(0, np.log(ja / j0a))
         # Eqn 12 - cathode activation energy
-        V_act_c = bc * np.log(jc / j0c)
+        V_act_c = bc * np.maximum(0, np.log(jc / j0c))
 
         return V_act_a, V_act_c
 
@@ -567,6 +601,9 @@ def calc_open_circuit_voltage(self, T_C):
         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
@@ -616,8 +653,7 @@ def calc_mass_flow_rate(self, T_C, I):
 
         eta_F = self.calc_faradaic_efficiency(T_C, I)
         # Eqn 10 [mol/sec]
-        h2_prod_mol = eta_F * self.n_cells * I / (self.z * F)
-        # n_cells is number of cells in series
+        h2_prod_mol = eta_F * I / (self.z * F)
         mfr = self.M_H * self.z * h2_prod_mol  # [g/sec]
         # z is valency number of electrons transferred per ion
         # for oxygen, z=4