From 1502d0bf357b0c9f9863175cd975670ec9b27449 Mon Sep 17 00:00:00 2001 From: Eric Fell Date: Thu, 25 Jan 2024 13:47:20 -0500 Subject: [PATCH] refactor names for # of electrons, total cell overpotential. Simplify error messages --- src/rfbzero/experiment.py | 39 ++++++++++------ src/rfbzero/redox_flow_cell.py | 83 ++++++++++++++++++---------------- tests/test_experiment.py | 16 +++---- tests/test_redox_flow_cell.py | 12 ++--- 4 files changed, 84 insertions(+), 66 deletions(-) diff --git a/src/rfbzero/experiment.py b/src/rfbzero/experiment.py index 8404ff0..15db800 100644 --- a/src/rfbzero/experiment.py +++ b/src/rfbzero/experiment.py @@ -74,7 +74,7 @@ def __init__(self, duration: float, time_increment: float, charge_first: bool = #: The combined (CLS+NCLS) mass transport overpotential (V), at each time step. self.mt: list[float] = [0.0] * self.max_steps #: The total cell overpotential (V), at each time step. - self.loss: list[float] = [0.0] * self.max_steps + self.total_overpotential: list[float] = [0.0] * self.max_steps # Total number of cycles is unknown at start, thus list sizes are undetermined self.capacity = 0.0 @@ -107,7 +107,7 @@ def _record_step( ocv: float, n_act: float = 0.0, n_mt: float = 0.0, - losses: float = 0.0 + total_overpotential: float = 0.0 ) -> None: """Records simulation data at valid time steps.""" # Update capacity @@ -119,10 +119,10 @@ def _record_step( self.ocv[self.step] = ocv self.step_is_charge[self.step] = charge - # Record overpotentials and total loss + # Record overpotentials and total total_overpotential self.act[self.step] = n_act self.mt[self.step] = n_mt - self.loss[self.step] = losses + self.total_overpotential[self.step] = total_overpotential # Record species concentrations cls_ox = cell_model.c_ox_cls @@ -185,7 +185,7 @@ def _finalize(self) -> None: self.act = self.act[:self.step] self.mt = self.mt[:self.step] - self.loss = self.loss[:self.step] + self.total_overpotential = self.total_overpotential[:self.step] class CycleStatus(str, Enum): @@ -317,9 +317,13 @@ def validate(self) -> CycleStatus: if abs(self.current) >= min(self.current_lim_cls, self.current_lim_ncls): return CycleStatus.LIMITING_CURRENT_REACHED - losses, *_ = self.cell_model._total_overpotential(self.current, self.current_lim_cls, self.current_lim_ncls) + total_overpotential, *_ = self.cell_model._total_overpotential( + self.current, + self.current_lim_cls, + self.current_lim_ncls + ) ocv = self.cell_model._open_circuit_voltage() - cell_v = self.cell_model._cell_voltage(ocv, losses, self.charge) + cell_v = self.cell_model._cell_voltage(ocv, total_overpotential, self.charge) if self.charge and cell_v >= self.voltage_limit or not self.charge and cell_v <= self.voltage_limit: return CycleStatus.VOLTAGE_LIMIT_REACHED @@ -344,10 +348,10 @@ def cycle_step(self) -> CycleStatus: return self._check_capacity(CycleStatus.NEGATIVE_CONCENTRATIONS) # Calculate overpotentials and the resulting cell voltage - losses, n_act, n_mt = self.cell_model._total_overpotential( + total_overpotential, n_act, n_mt = self.cell_model._total_overpotential( self.current, self.current_lim_cls, self.current_lim_ncls) ocv = self.cell_model._open_circuit_voltage() - cell_v = self.cell_model._cell_voltage(ocv, losses, self.charge) + cell_v = self.cell_model._cell_voltage(ocv, total_overpotential, self.charge) # Check if the voltage limit is reached if self.charge and cell_v >= self.voltage_limit or not self.charge and cell_v <= self.voltage_limit: @@ -356,7 +360,16 @@ def cycle_step(self) -> CycleStatus: cycle_status = self._check_capacity(cycle_status) # Update results - self.results._record_step(self.cell_model, self.charge, self.current, cell_v, ocv, n_act, n_mt, losses) + self.results._record_step( + self.cell_model, + self.charge, + self.current, + cell_v, + ocv, + n_act, + n_mt, + total_overpotential + ) return self._check_time(cycle_status) @@ -442,7 +455,7 @@ def __current_direction(self) -> int: def __find_min_current(self, ocv: float) -> None: """ Solves the current at a given timestep of constant voltage cycling. - Attempts to minimize the difference of voltage, OCV, and losses (function of current). + Attempts to minimize the difference of voltage, OCV, and total overpotential (function of current). """ @@ -600,7 +613,7 @@ class ConstantCurrent(CyclingProtocol): voltage_limit_discharge : float Voltage below which cell will switch to charge (V). current : float - Instantaneous current flowing (A). + Current (A) value used for charging. Negative of this value used for discharging current. current_charge : float Desired charging current for CC cycling (A). current_discharge : float @@ -823,7 +836,7 @@ class ConstantCurrentConstantVoltage(CyclingProtocol): current_cutoff_discharge : float Current above which CV discharging will switch to CC portion of CCCV charge (A). current : float - Instantaneous current flowing (A). + Current (A) value used for charging. Negative of this value used for discharging current. current_charge : float Desired charging current for CC cycling (A). current_discharge : float diff --git a/src/rfbzero/redox_flow_cell.py b/src/rfbzero/redox_flow_cell.py index afa3284..1a5a5ef 100644 --- a/src/rfbzero/redox_flow_cell.py +++ b/src/rfbzero/redox_flow_cell.py @@ -55,6 +55,7 @@ class ZeroDModel: Default is 5.0, a typical lab-scale cell size. cls_negolyte : bool True if negolyte is the CLS, False if posolyte is the CLS. + Default is True. time_increment : float Simulation time step (s). Default is 0.01, providing adequate balance of accuracy vs compute time. @@ -65,10 +66,12 @@ class ZeroDModel: Roughness factor, dimensionless. Total surface area divided by geometric surface area. Default is 26.0, as used in [1]. - n_cls : int + num_electrons_cls : int Number of electrons transferred per active species molecule in the CLS. - n_ncls : int + Default is 1. + num_electrons_ncls : int Number of electrons transferred per active species molecule in the NCLS. + Default is 1. temperature : float Temperature of battery and electrolytes (K). Default is 298 K (25C). @@ -101,8 +104,8 @@ def __init__( time_increment: float = 0.01, k_mt: float = 0.8, roughness_factor: float = 26.0, - n_cls: int = 1, - n_ncls: int = 1, + num_electrons_cls: int = 1, + num_electrons_ncls: int = 1, temperature: float = 298.0 ) -> None: self.cls_volume = cls_volume @@ -122,8 +125,8 @@ def __init__( self.time_increment = time_increment self.k_mt = k_mt self.const_i_ex = F * roughness_factor * self.geometric_area - self.n_cls = n_cls - self.n_ncls = n_ncls + self.num_electrons_cls = num_electrons_cls + self.num_electrons_ncls = num_electrons_ncls self.prev_c_ox_cls = self.c_ox_cls self.prev_c_red_cls = self.c_red_cls @@ -134,13 +137,14 @@ def __init__( self.nernst_const = (R * temperature) / F - for key, value in {'cls_volume': self.cls_volume, 'ncls_volume': self.ncls_volume, 'k_0_cls': self.k_0_cls, + for key, value in {'cls_volume': self.cls_volume, 'ncls_volume': self.ncls_volume, 'cls_start_c_ox': self.c_ox_cls, 'cls_start_c_red': self.c_red_cls, 'ncls_start_c_ox': self.c_ox_ncls, 'ncls_start_c_red': self.c_red_ncls, + 'ocv_50_soc': self.ocv_50_soc, 'resistance': self.resistance, 'k_0_cls': self.k_0_cls, 'k_0_ncls': self.k_0_ncls, 'geometric_area': self.geometric_area, 'time_increment': self.time_increment, 'k_mt': self.k_mt, 'const_i_ex': self.const_i_ex, - 'ocv_50_soc': self.ocv_50_soc, 'resistance': self.resistance, 'n_cls': self.n_cls, - 'n_ncls': self.n_ncls, 'temperature': temperature}.items(): + 'num_electrons_cls': self.num_electrons_cls, 'num_electrons_ncls': self.num_electrons_ncls, + 'temperature': temperature}.items(): if key not in ['ocv_50_soc', 'resistance', 'cls_start_c_ox', 'cls_start_c_red', @@ -153,17 +157,18 @@ def __init__( if not 0.0 < self.alpha_cls < 1.0 or not 0.0 < self.alpha_ncls < 1.0: raise ValueError("Alpha parameters must be between 0.0 and 1.0") - if not isinstance(self.n_cls, int) or not isinstance(self.n_ncls, int): - raise ValueError("'n_cls' and 'n_ncls' must be >= 1") + if not isinstance(self.num_electrons_cls, int) or not isinstance(self.num_electrons_ncls, int): + raise ValueError("'num_electrons_cls' and 'num_electrons_ncls' must be >= 1") if self.ocv_50_soc == 0.0 and self.cls_volume >= self.ncls_volume: raise ValueError("'cls_volume' must be < 'ncls_volume' in a symmetric cell") - if self.ocv_50_soc == 0.0 and self.n_cls != self.n_ncls: - raise ValueError("Symmetric cell (0 volt OCV) requires 'n_cls' and 'n_ncls' to be equal (same species)") + if self.ocv_50_soc == 0.0 and self.num_electrons_cls != self.num_electrons_ncls: + raise ValueError("Symmetric cell ('ocv_50_soc' = 0) requires 'num_electrons_cls' and 'num_electrons_ncls' " + "to be equal (same species)") - self.init_cls_capacity = self.cls_volume * self.n_cls * (self.c_ox_cls + self.c_red_cls) - init_ncls_capacity = self.ncls_volume * self.n_ncls * (self.c_ox_ncls + self.c_red_ncls) + self.init_cls_capacity = self.cls_volume * self.num_electrons_cls * (self.c_ox_cls + self.c_red_cls) + init_ncls_capacity = self.ncls_volume * self.num_electrons_ncls * (self.c_ox_ncls + self.c_red_ncls) if self.init_cls_capacity >= init_ncls_capacity: raise ValueError("Initial capacity of CLS must be less than initial capacity of NCLS") @@ -174,7 +179,6 @@ def __init__( def _exchange_current(self) -> tuple[float, float]: """ Calculates exchange current (i_0) of redox couples in the CLS and NCLS. - Value returned is in Amps. Returns ------- @@ -185,16 +189,16 @@ def _exchange_current(self) -> tuple[float, float]: """ # division by 1000 for conversion from L to cm^3 - i_0_cls = (self.n_cls * self.const_i_ex * self.k_0_cls * (self.c_red_cls ** self.alpha_cls) + i_0_cls = (self.num_electrons_cls * self.const_i_ex * self.k_0_cls * (self.c_red_cls ** self.alpha_cls) * (self.c_ox_cls ** (1 - self.alpha_cls)) * 0.001) - i_0_ncls = (self.n_ncls * self.const_i_ex * self.k_0_ncls * (self.c_red_ncls ** self.alpha_ncls) + i_0_ncls = (self.num_electrons_ncls * self.const_i_ex * self.k_0_ncls * (self.c_red_ncls ** self.alpha_ncls) * (self.c_ox_ncls ** (1 - self.alpha_ncls)) * 0.001) return i_0_cls, i_0_ncls def _limiting_current(self, c_lim: float) -> float: """ Calculates limiting current (i_lim) for a single reservoir. - Value returned is in Amps. This is equation 6 of [1]. + This is equation 6 of [1]. """ # div by 1000 for conversion from L to cm^3 @@ -219,11 +223,11 @@ def _limiting_concentration(self, charge: bool) -> tuple[float, float]: """ if self.cls_negolyte == charge: - i_lim_cls = self._limiting_current(self.c_ox_cls) * self.n_cls - i_lim_ncls = self._limiting_current(self.c_red_ncls) * self.n_ncls + i_lim_cls = self._limiting_current(self.c_ox_cls) * self.num_electrons_cls + i_lim_ncls = self._limiting_current(self.c_red_ncls) * self.num_electrons_ncls else: - i_lim_cls = self._limiting_current(self.c_red_cls) * self.n_cls - i_lim_ncls = self._limiting_current(self.c_ox_ncls) * self.n_ncls + i_lim_cls = self._limiting_current(self.c_red_cls) * self.num_electrons_cls + i_lim_ncls = self._limiting_current(self.c_ox_ncls) * self.num_electrons_ncls return i_lim_cls, i_lim_ncls @@ -235,7 +239,7 @@ def _activation_overpotential(self, current: float, i_0_cls: float, i_0_ncls: fl Parameters ---------- current : float - Instantaneous current flowing (A). + Instantaneous current flowing (A). Positive if charging, negative if discharging. i_0_cls : float Exchange current of CLS redox couple at a given timestep (A). i_0_ncls : float @@ -250,8 +254,8 @@ def _activation_overpotential(self, current: float, i_0_cls: float, i_0_ncls: fl z_cls = abs(current) / (2 * i_0_cls) z_ncls = abs(current) / (2 * i_0_ncls) - n_act = self.nernst_const * ((log(z_cls + ((z_cls ** 2) + 1) ** 0.5) / self.n_cls) - + (log(z_ncls + ((z_ncls ** 2) + 1) ** 0.5) / self.n_ncls)) + n_act = self.nernst_const * ((log(z_cls + ((z_cls ** 2) + 1) ** 0.5) / self.num_electrons_cls) + + (log(z_ncls + ((z_ncls ** 2) + 1) ** 0.5) / self.num_electrons_ncls)) return n_act def _negative_concentrations(self) -> bool: @@ -285,7 +289,7 @@ def __mass_transport_overpotential(self, current: float, i_lim_cls: float, i_lim c_tot_cls = self.c_red_cls + self.c_ox_cls c_tot_ncls = self.c_red_ncls + self.c_ox_ncls - charge = current >= 0 + charge = current >= 0.0 if self.cls_negolyte == charge: c1_cls = self.c_ox_cls c2_cls = self.c_red_cls @@ -300,8 +304,8 @@ def __mass_transport_overpotential(self, current: float, i_lim_cls: float, i_lim i = abs(current) n_mt = self.nernst_const * ( - (log(1 - ((c_tot_cls * i) / ((c1_cls * i_lim_cls) + (c2_cls * i)))) / self.n_cls) - + (log(1 - ((c_tot_ncls * i) / ((c1_ncls * i_lim_ncls) + (c2_ncls * i)))) / self.n_ncls) + (log(1 - ((c_tot_cls * i) / ((c1_cls * i_lim_cls) + (c2_cls * i)))) / self.num_electrons_cls) + + (log(1 - ((c_tot_ncls * i) / ((c1_ncls * i_lim_ncls) + (c2_ncls * i)))) / self.num_electrons_ncls) ) return n_mt * -1 @@ -321,7 +325,7 @@ def _total_overpotential(self, current: float, i_lim_cls: float, i_lim_ncls: flo Returns ------- - n_loss : float + total_overpotential : float Total cell overpotential (V). n_act : float Total activation overpotential (V). @@ -336,9 +340,9 @@ def _total_overpotential(self, current: float, i_lim_cls: float, i_lim_ncls: flo n_act = self._activation_overpotential(current, i_0_cls, i_0_ncls) n_mt = self.__mass_transport_overpotential(current, i_lim_cls, i_lim_ncls) - n_loss = n_ohmic + n_act + n_mt + total_overpotential = n_ohmic + n_act + n_mt - return n_loss, n_act, n_mt + return total_overpotential, n_act, n_mt def _open_circuit_voltage(self) -> float: """ @@ -358,14 +362,15 @@ def _open_circuit_voltage(self) -> float: direction = 1 if self.cls_negolyte else -1 ocv = (self.ocv_50_soc - + direction * (((self.nernst_const / self.n_cls) * log(self.c_red_cls / self.c_ox_cls)) - + ((self.nernst_const / self.n_ncls) * log(self.c_ox_ncls / self.c_red_ncls)))) + + direction * (((self.nernst_const / self.num_electrons_cls) * log(self.c_red_cls / self.c_ox_cls)) + + ((self.nernst_const / self.num_electrons_ncls) * log( + self.c_ox_ncls / self.c_red_ncls)))) return ocv @staticmethod - def _cell_voltage(ocv: float, losses: float, charge: bool) -> float: + def _cell_voltage(ocv: float, total_overpotential: float, charge: bool) -> float: """If charging, add overpotentials to OCV, else subtract them.""" - return ocv + losses if charge else ocv - losses + return ocv + total_overpotential if charge else ocv - total_overpotential def _coulomb_counter( self, @@ -381,7 +386,7 @@ def _coulomb_counter( Parameters ---------- current : float - Instantaneous current flowing (A). + Instantaneous current flowing (A). Positive if charging, negative if discharging. cls_degradation : DegradationMechanism, optional Degradation mechanism for CLS. ncls_degradation: DegradationMechanism, optional @@ -400,8 +405,8 @@ def _coulomb_counter( # Change in concentration from coulomb counting based solely on current direction = 1 if self.cls_negolyte else -1 - delta_cls = ((self.time_increment * current) / (F * self.n_cls * self.cls_volume)) * direction - delta_ncls = ((self.time_increment * current) / (F * self.n_ncls * self.ncls_volume)) * direction + delta_cls = ((self.time_increment * current) / (F * self.num_electrons_cls * self.cls_volume)) * direction + delta_ncls = ((self.time_increment * current) / (F * self.num_electrons_ncls * self.ncls_volume)) * direction # update CLS and NCLS concentrations c_ox_cls = self.c_ox_cls - delta_cls diff --git a/tests/test_experiment.py b/tests/test_experiment.py index 437e535..0b64ba4 100644 --- a/tests/test_experiment.py +++ b/tests/test_experiment.py @@ -26,8 +26,8 @@ def test_cc(self): resistance=1.0, # ohms k_0_cls=1e-3, # cm/s k_0_ncls=1e-3, # cm/s - n_cls=1, # electrons - n_ncls=1, # electrons + num_electrons_cls=1, # electrons + num_electrons_ncls=1, # electrons ) # define degradation mechanisms @@ -62,8 +62,8 @@ def test_cv1(self): resistance=1.0, # ohms k_0_cls=1e-3, # cm/s k_0_ncls=1e-3, # cm/s - n_cls=1, # electrons - n_ncls=1, # electrons + num_electrons_cls=1, # electrons + num_electrons_ncls=1, # electrons ) # define degradation mechanisms @@ -104,8 +104,8 @@ def test_cccv_full_cell(self): resistance=0.8, # ohms k_0_cls=1e-3, # cm/s k_0_ncls=1e-3, # cm/s - n_cls=1, # electrons - n_ncls=1, # electrons + num_electrons_cls=1, # electrons + num_electrons_ncls=1, # electrons ) # define degradation mechanisms @@ -143,8 +143,8 @@ def test_cccv_symmetric_cell(self): resistance=0.8, # ohms k_0_cls=1e-3, # cm/s k_0_ncls=1e-3, # cm/s - n_cls=1, # electrons - n_ncls=1, # electrons + num_electrons_cls=1, # electrons + num_electrons_ncls=1, # electrons ) # define degradation mechanisms diff --git a/tests/test_redox_flow_cell.py b/tests/test_redox_flow_cell.py index 772ba51..ce50bb2 100644 --- a/tests/test_redox_flow_cell.py +++ b/tests/test_redox_flow_cell.py @@ -40,8 +40,8 @@ def test_class_init(self, v_cls, v_ncls, ox_cls, red_cls, ox_ncls, red_ncls, ocv k_0_ncls=k_n, alpha_cls=a_c, alpha_ncls=a_n, - n_cls=n_c, - n_ncls=n_n, + num_electrons_cls=n_c, + num_electrons_ncls=n_n, ) @pytest.mark.parametrize("v_cls,v_ncls,ox_cls,red_cls,ox_ncls,red_ncls,ocv,res,k_c,k_n,time_i", @@ -71,7 +71,7 @@ def test_timestep_warning(self,v_cls,v_ncls,ox_cls,red_cls,ox_ncls,red_ncls,ocv, def test_exchange_current(self): cell = ZeroDModel(cls_volume=0.005, ncls_volume=0.01, cls_start_c_ox=0.01, cls_start_c_red=0.01, ncls_start_c_ox=0.01, ncls_start_c_red=0.01, ocv_50_soc=1.0, resistance=1, k_0_cls=1e-3, - k_0_ncls=1e-3, n_ncls=2) + k_0_ncls=1e-3, num_electrons_ncls=2) i_0_cls, i_0_ncls = cell._exchange_current() assert np.isclose(i_0_cls, 0.12543093175) assert np.isclose(i_0_ncls, 0.25086186351) @@ -80,7 +80,7 @@ def test_limiting_current(self): limiting_c = 0.2 cell = ZeroDModel(cls_volume=0.005, ncls_volume=0.01, cls_start_c_ox=0.01, cls_start_c_red=0.01, ncls_start_c_ox=0.01, ncls_start_c_red=0.01, ocv_50_soc=1.0, resistance=1, k_0_cls=1e-3, - k_0_ncls=1e-3, n_ncls=2) + k_0_ncls=1e-3, num_electrons_ncls=2) i_lim = cell._limiting_current(limiting_c) assert np.isclose(i_lim, 77.1882656) @@ -88,7 +88,7 @@ def test_limiting_current(self): def test_limiting_concentration(self): cell = ZeroDModel(cls_volume=0.005, ncls_volume=0.01, cls_start_c_ox=0.01, cls_start_c_red=0.02, ncls_start_c_ox=0.02, ncls_start_c_red=0.01, ocv_50_soc=1.0, resistance=1, k_0_cls=1e-3, - k_0_ncls=1e-3, n_ncls=2) + k_0_ncls=1e-3, num_electrons_ncls=2) i_lim_cls, i_lim_ncls = cell._limiting_concentration(True) assert np.isclose(i_lim_cls, 3.85941328) @@ -101,7 +101,7 @@ def test_limiting_concentration(self): def test_activation_overpotential(self): cell = ZeroDModel(cls_volume=0.005, ncls_volume=0.01, cls_start_c_ox=0.01, cls_start_c_red=0.01, ncls_start_c_ox=0.01, ncls_start_c_red=0.01, ocv_50_soc=1.0, resistance=1, k_0_cls=1e-3, - k_0_ncls=1e-3, n_ncls=2) + k_0_ncls=1e-3, num_electrons_ncls=2) current = 1 i_0_cls = 0.01 i_0_ncls = 0.01