diff --git a/CHANGELOG.md b/CHANGELOG.md index c99486239e..bb902bdfbf 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -7,6 +7,7 @@ ## Features +- Added `"coupled ion-electron transfer"` as an `"intercalation kinetics"` option, following Fraggedakis et al. (2021). It uses the existing reorganization energy and a new temperature-dependent `"... electrode coupled ion-electron transfer rate constant [A.m-2]"` per electrode, and its current saturates at large overpotential. ([#5993](https://github.com/pybamm-team/PyBaMM/pull/5993)) - Replaced Codacy with open checks that also run locally. Ruff's `flake8-bandit` (`S`) rules replace Bandit, and CodeQL scans the Python and C++ code. `nox -s type-check` runs Pyrefly over `pybamm`, which now type-checks cleanly, and a new type error fails CI. clang-tidy reports solver findings to code scanning as an advisory check. ([#5974](https://github.com/pybamm-team/PyBaMM/pull/5974)) ## Bug fixes diff --git a/docs/source/api/models/submodels/interface/kinetics/coupled_ion_electron_transfer.rst b/docs/source/api/models/submodels/interface/kinetics/coupled_ion_electron_transfer.rst new file mode 100644 index 0000000000..ba01c52ebe --- /dev/null +++ b/docs/source/api/models/submodels/interface/kinetics/coupled_ion_electron_transfer.rst @@ -0,0 +1,7 @@ +Coupled ion-electron transfer +============================= + +.. autoclass:: pybamm.kinetics.CoupledIonElectronTransfer + :members: + +.. footbibliography:: diff --git a/docs/source/api/models/submodels/interface/kinetics/index.rst b/docs/source/api/models/submodels/interface/kinetics/index.rst index efb8be4d30..fd61b78fb1 100644 --- a/docs/source/api/models/submodels/interface/kinetics/index.rst +++ b/docs/source/api/models/submodels/interface/kinetics/index.rst @@ -5,6 +5,7 @@ Kinetics base_kinetics butler_volmer + coupled_ion_electron_transfer diffusion_limited linear marcus diff --git a/packages/pybamm/src/pybamm/CITATIONS.bib b/packages/pybamm/src/pybamm/CITATIONS.bib index 2e5c10c125..7753061bfe 100644 --- a/packages/pybamm/src/pybamm/CITATIONS.bib +++ b/packages/pybamm/src/pybamm/CITATIONS.bib @@ -251,6 +251,17 @@ @article{Fan2022 publisher = {Elsevier} } +@article{Fraggedakis2021, + title = {Theory of coupled ion-electron transfer kinetics}, + author = {Fraggedakis, Dimitrios and McEldrew, Michael and Smith, Raymond B and Krishnan, Yamini and Zhang, Yirui and Bai, Peng and Chueh, William C and Shao-Horn, Yang and Bazant, Martin Z}, + journal = {Electrochimica Acta}, + volume = {367}, + pages = {137432}, + year = {2021}, + publisher = {Elsevier}, + doi = {10.1016/j.electacta.2020.137432} +} + @article{Gustafsson2020, doi = {10.21105/joss.02369}, year = {2020}, @@ -904,6 +915,17 @@ @article{Yang2017 doi = {10.1016/j.jpowsour.2017.05.110} } +@article{Zeng2014, + title = {Simple formula for {Marcus--Hush--Chidsey} kinetics}, + author = {Zeng, Yi and Smith, Raymond B and Bai, Peng and Bazant, Martin Z}, + journal = {Journal of Electroanalytical Chemistry}, + volume = {735}, + pages = {77--83}, + year = {2014}, + publisher = {Elsevier}, + doi = {10.1016/j.jelechem.2014.09.038} +} + @article{Zhao2018, title = {Modeling the effects of thermal gradients induced by tab and surface cooling on lithium ion cell performance}, author = {Zhao, Yan and Patel, Yatish and Zhang, Teng and Offer, Gregory J}, diff --git a/packages/pybamm/src/pybamm/models/full_battery_models/base_battery_model.py b/packages/pybamm/src/pybamm/models/full_battery_models/base_battery_model.py index 15d8d347b0..58004502f5 100644 --- a/packages/pybamm/src/pybamm/models/full_battery_models/base_battery_model.py +++ b/packages/pybamm/src/pybamm/models/full_battery_models/base_battery_model.py @@ -160,6 +160,8 @@ class BatteryModelOptions(pybamm.FuzzyDict): Model for intercalation kinetics. Can be "symmetric Butler-Volmer" (default), "asymmetric Butler-Volmer", "linear", "Marcus", "Marcus-Hush-Chidsey" (which uses the asymptotic form from Zeng 2014), + "coupled ion-electron transfer" (which uses the form from Fraggedakis + 2021), or "MSMR" (which uses the form from Baker 2018). A 2-tuple can be provided for different behaviour in negative and positive electrodes. * "interface utilisation": str @@ -353,6 +355,7 @@ def __init__(self, extra_options): "linear", "Marcus", "Marcus-Hush-Chidsey", + "coupled ion-electron transfer", "MSMR", ], "interface utilisation": ["full", "constant", "current-driven"], @@ -1388,6 +1391,8 @@ def get_intercalation_kinetics(self, domain): return pybamm.kinetics.Marcus elif options["intercalation kinetics"] == "Marcus-Hush-Chidsey": return pybamm.kinetics.MarcusHushChidsey + elif options["intercalation kinetics"] == "coupled ion-electron transfer": + return pybamm.kinetics.CoupledIonElectronTransfer elif options["intercalation kinetics"] == "MSMR": return pybamm.kinetics.MSMRButlerVolmer diff --git a/packages/pybamm/src/pybamm/models/submodels/interface/base_interface.py b/packages/pybamm/src/pybamm/models/submodels/interface/base_interface.py index 6624b646f8..759b0198fa 100644 --- a/packages/pybamm/src/pybamm/models/submodels/interface/base_interface.py +++ b/packages/pybamm/src/pybamm/models/submodels/interface/base_interface.py @@ -49,6 +49,65 @@ def __init__(self, param, domain, reaction, options=None, phase="primary"): domain_options = getattr(self.options, domain) self.size_distribution = domain_options["particle size"] == "distribution" + def _get_intercalation_concentrations(self, variables): + """ + Get the electrolyte concentration, particle surface concentration and + temperature seen by the main intercalation reaction. + + Parameters + ---------- + variables : dict + The variables in the full model. + + Returns + ------- + c_e, c_s_surf, T : :class:`pybamm.Symbol` + Electrolyte concentration, particle surface concentration and + temperature, on matching domains. + """ + domain, Domain = self.domain_Domain + phase_name = self.phase_name + c_e = variables[f"{Domain} electrolyte concentration [mol.m-3]"] + T = variables[f"{Domain} electrode temperature [K]"] + + # For "particle-size distribution" submodels, take distribution version + # of c_s_surf that depends on particle size. + if self.size_distribution: + c_s_surf = variables[ + f"{Domain} {phase_name}particle surface " + "concentration distribution [mol.m-3]" + ] + # If all variables were broadcast (in "x"), take only the orphans, + # then re-broadcast c_e + if ( + isinstance(c_s_surf, pybamm.Broadcast) + and isinstance(c_e, pybamm.Broadcast) + and isinstance(T, pybamm.Broadcast) + ): + c_s_surf = c_s_surf.orphans[0] + c_e = c_e.orphans[0] + T = T.orphans[0] + + # as c_e must now be a scalar, re-broadcast to "current collector" + c_e = pybamm.PrimaryBroadcast(c_e, ["current collector"]) + # broadcast c_e, T onto "particle size" + c_e = pybamm.PrimaryBroadcast(c_e, [f"{domain} {phase_name}particle size"]) + T = pybamm.PrimaryBroadcast(T, [f"{domain} {phase_name}particle size"]) + else: + c_s_surf = variables[ + f"{Domain} {phase_name}particle surface concentration [mol.m-3]" + ] + # If all variables were broadcast, take only the orphans + if ( + isinstance(c_s_surf, pybamm.Broadcast) + and isinstance(c_e, pybamm.Broadcast) + and isinstance(T, pybamm.Broadcast) + ): + c_s_surf = c_s_surf.orphans[0] + c_e = c_e.orphans[0] + T = T.orphans[0] + return c_e, c_s_surf, T + def _get_exchange_current_density(self, variables): """ A private function to obtain the exchange current density @@ -65,54 +124,13 @@ def _get_exchange_current_density(self, variables): """ phase_param = self.phase_param domain, Domain = self.domain_Domain - phase_name = self.phase_name domain_options = getattr(self.options, domain) c_e = variables[f"{Domain} electrolyte concentration [mol.m-3]"] T = variables[f"{Domain} electrode temperature [K]"] if self.reaction == "lithium-ion main": - # For "particle-size distribution" submodels, take distribution version - # of c_s_surf that depends on particle size. - domain_options = getattr(self.options, domain) - if domain_options["particle size"] == "distribution": - c_s_surf = variables[ - f"{Domain} {phase_name}particle surface " - "concentration distribution [mol.m-3]" - ] - # If all variables were broadcast (in "x"), take only the orphans, - # then re-broadcast c_e - if ( - isinstance(c_s_surf, pybamm.Broadcast) - and isinstance(c_e, pybamm.Broadcast) - and isinstance(T, pybamm.Broadcast) - ): - c_s_surf = c_s_surf.orphans[0] - c_e = c_e.orphans[0] - T = T.orphans[0] - - # as c_e must now be a scalar, re-broadcast to - # "current collector" - c_e = pybamm.PrimaryBroadcast(c_e, ["current collector"]) - # broadcast c_e, T onto "particle size" - c_e = pybamm.PrimaryBroadcast( - c_e, [f"{domain} {phase_name}particle size"] - ) - T = pybamm.PrimaryBroadcast(T, [f"{domain} {phase_name}particle size"]) - - else: - c_s_surf = variables[ - f"{Domain} {phase_name}particle surface concentration [mol.m-3]" - ] - # If all variables were broadcast, take only the orphans - if ( - isinstance(c_s_surf, pybamm.Broadcast) - and isinstance(c_e, pybamm.Broadcast) - and isinstance(T, pybamm.Broadcast) - ): - c_s_surf = c_s_surf.orphans[0] - c_e = c_e.orphans[0] - T = T.orphans[0] + c_e, c_s_surf, T = self._get_intercalation_concentrations(variables) # Get main reaction exchange-current density (may have empirical hysteresis) j0_option = getattr(domain_options, self.phase)["exchange-current density"] if j0_option == "single": diff --git a/packages/pybamm/src/pybamm/models/submodels/interface/kinetics/__init__.py b/packages/pybamm/src/pybamm/models/submodels/interface/kinetics/__init__.py index 21bf2da7f1..f91b8271de 100644 --- a/packages/pybamm/src/pybamm/models/submodels/interface/kinetics/__init__.py +++ b/packages/pybamm/src/pybamm/models/submodels/interface/kinetics/__init__.py @@ -3,6 +3,7 @@ from .butler_volmer import SymmetricButlerVolmer, AsymmetricButlerVolmer from .linear import Linear from .marcus import Marcus, MarcusHushChidsey +from .coupled_ion_electron_transfer import CoupledIonElectronTransfer from .tafel import ForwardTafel # , BackwardTafel from .no_reaction import NoReaction from .msmr_butler_volmer import MSMRButlerVolmer @@ -15,6 +16,7 @@ from .inverse_kinetics.inverse_butler_volmer import InverseButlerVolmer from .inverse_kinetics.inverse_linear import InverseLinear -__all__ = ['base_kinetics', 'butler_volmer', 'diffusion_limited', +__all__ = ['base_kinetics', 'butler_volmer', 'coupled_ion_electron_transfer', + 'diffusion_limited', 'inverse_kinetics', 'linear', 'marcus', 'msmr_butler_volmer', 'no_reaction', 'tafel', 'total_main_kinetics'] diff --git a/packages/pybamm/src/pybamm/models/submodels/interface/kinetics/coupled_ion_electron_transfer.py b/packages/pybamm/src/pybamm/models/submodels/interface/kinetics/coupled_ion_electron_transfer.py new file mode 100644 index 0000000000..f613e2e3ae --- /dev/null +++ b/packages/pybamm/src/pybamm/models/submodels/interface/kinetics/coupled_ion_electron_transfer.py @@ -0,0 +1,123 @@ +# +# Coupled ion-electron transfer kinetics class +# +from __future__ import annotations + +import numpy as np + +import pybamm + +from .base_kinetics import BaseKinetics + + +class CoupledIonElectronTransfer(BaseKinetics): + """ + Intercalation rate law in which the ion and electron transfer in one concerted + step :footcite:t:`Fraggedakis2021`, using the asymptotic Marcus-Hush-Chidsey + integral of :footcite:t:`Zeng2014`. In units of :math:`RT/F`, + + .. math:: + j = k_0 \\sqrt{\\pi\\lambda}\\,\\mathrm{erfc}\\left(\\frac{\\lambda - \\sqrt{1 + \\sqrt{\\lambda} + \\eta_f^2}}{2\\sqrt{\\lambda}}\\right) \\left[\\frac{\\theta}{1 + e^{-\\eta_f}} - \\frac{\\tilde{c}_e (1 - \\theta)}{1 + e^{\\eta_f}}\\right], + + where :math:`\\eta_f = \\eta_r + \\ln(\\tilde{c}_e (1 - \\theta) / \\theta)`, + :math:`\\theta = c_{s,surf}/c_{max}` and :math:`\\tilde{c}_e = c_e/c_{e,init}`. + Not available for a lithium metal electrode, which has no host lattice. + + Parameters + ---------- + param : parameter class + model parameters + domain : str + The domain to implement the model, either: 'Negative' or 'Positive'. + reaction : str + The name of the reaction being implemented + options: dict + A dictionary of options to be passed to the model. + See :class:`pybamm.BaseBatteryModel` + phase : str, optional + Phase of the particle (default is "primary") + """ + + def __init__(self, param, domain, reaction, options, phase="primary"): + if reaction == "lithium metal plating": + raise pybamm.OptionError( + "Coupled ion-electron transfer kinetics need a host lattice, so " + "they cannot be used for the lithium metal electrode of a half-cell. " + "Set 'intercalation kinetics' per electrode instead, e.g. " + "('Marcus-Hush-Chidsey', 'coupled ion-electron transfer') for a " + "positive working electrode" + ) + if reaction != "lithium-ion main": + raise pybamm.OptionError( + "Coupled ion-electron transfer kinetics are only implemented for " + f"lithium-ion intercalation, not '{reaction}'" + ) + super().__init__(param, domain, reaction, options, phase) + domain_options = getattr(self.options, self.domain) + if getattr(domain_options, self.phase)["exchange-current density"] != "single": + raise pybamm.OptionError( + "Coupled ion-electron transfer kinetics do not use the " + "exchange-current density parameter, so 'exchange-current density' " + "must be 'single'" + ) + pybamm.citations.register("Fraggedakis2021") + pybamm.citations.register("Zeng2014") + + def _get_exchange_current_density(self, variables): + c_e, c_s_surf, T = self._get_intercalation_concentrations(variables) + phase_param = self.phase_param + tol_c_s = pybamm.settings.tolerances["j0__c_s"] + tol_c_e = pybamm.settings.tolerances["j0__c_e"] + stoichiometry = pybamm.maximum( + pybamm.minimum(c_s_surf / phase_param.c_max, 1 - tol_c_s), tol_c_s + ) + relative_c_e = pybamm.maximum(c_e, tol_c_e) / phase_param.c_e_init + + # Kept for the rate law, which only receives j0 and the overpotential + self._stoichiometry = stoichiometry + self._relative_c_e = relative_c_e + + # Linear-response slope, so j ~ j0 F eta_r / RT as for Butler-Volmer + F_RT = self.param.F / (self.param.R * T) + eta_f_equilibrium = self._formal_shift(stoichiometry, relative_c_e) + return ( + relative_c_e + * (1 - stoichiometry) + * self._rate_factor(eta_f_equilibrium, F_RT, T) + * self._fermi(-eta_f_equilibrium) + ) + + @staticmethod + def _formal_shift(stoichiometry, relative_c_e): + return pybamm.log(relative_c_e * (1 - stoichiometry) / stoichiometry) + + @staticmethod + def _fermi(x): + """Return 1 / (1 + exp(-x)) without overflow.""" + return (1 + pybamm.tanh(x / 2)) / 2 + + def _rate_factor(self, eta_f, F_RT, T): + """Return k0 times the Fermi-weighted Marcus integral prefactor.""" + lambda_T = F_RT * self.phase_param.mhc_lambda + k0 = self.phase_param.k0_coupled_ion_electron_transfer(T) + a = 1 + pybamm.sqrt(lambda_T) + arg = (lambda_T - pybamm.sqrt(a + eta_f**2)) / (2 * pybamm.sqrt(lambda_T)) + return k0 * pybamm.sqrt(np.pi * lambda_T) * pybamm.erfc(arg) + + def _match_domains(self, symbol, eta_r): + domain = self.domain + if symbol.secondary_domain == [ + "current collector" + ] and eta_r.secondary_domain == [f"{domain} electrode"]: + symbol = pybamm.SecondaryBroadcast(symbol, [f"{domain} electrode"]) + return symbol + + def _get_kinetics(self, j0, ne, eta_r, T, u): + stoichiometry = self._match_domains(self._stoichiometry, eta_r) + relative_c_e = self._match_domains(self._relative_c_e, eta_r) + + F_RT = self.param.F / (self.param.R * T) + eta_f = F_RT * eta_r + self._formal_shift(stoichiometry, relative_c_e) + oxidation = stoichiometry * self._fermi(eta_f) + reduction = relative_c_e * (1 - stoichiometry) * self._fermi(-eta_f) + return u * self._rate_factor(eta_f, F_RT, T) * (oxidation - reduction) diff --git a/packages/pybamm/src/pybamm/parameters/lithium_ion_parameters.py b/packages/pybamm/src/pybamm/parameters/lithium_ion_parameters.py index 84c27f2fd1..c54cda276e 100644 --- a/packages/pybamm/src/pybamm/parameters/lithium_ion_parameters.py +++ b/packages/pybamm/src/pybamm/parameters/lithium_ion_parameters.py @@ -661,6 +661,16 @@ def j0(self, c_e, c_s_surf, T, lithiation=None): post_processor=regulariser, ) + def k0_coupled_ion_electron_transfer(self, T): + """Coupled ion-electron transfer rate constant [A.m-2]""" + Domain = self.domain.capitalize() + inputs = {"Temperature [K]": T} + return pybamm.FunctionParameter( + f"{self.phase_prefactor}{Domain} electrode coupled ion-electron " + "transfer rate constant [A.m-2]", + inputs, + ) + def j0_stripping(self, c_e, c_Li, T): """Dimensional exchange-current density for stripping [A.m-2]""" tol = pybamm.settings.tolerances["j0__c_e"] diff --git a/packages/pybamm/tests/integration/test_models/test_full_battery_models/test_lithium_ion/base_lithium_ion_tests.py b/packages/pybamm/tests/integration/test_models/test_full_battery_models/test_lithium_ion/base_lithium_ion_tests.py index 19d54ac51e..7e119acbcf 100644 --- a/packages/pybamm/tests/integration/test_models/test_full_battery_models/test_lithium_ion/base_lithium_ion_tests.py +++ b/packages/pybamm/tests/integration/test_models/test_full_battery_models/test_lithium_ion/base_lithium_ion_tests.py @@ -159,6 +159,22 @@ def test_kinetics_mhc(self): ) self.run_basic_processing_test(options, parameter_values=parameter_values) + def test_kinetics_coupled_ion_electron_transfer(self): + options = {"intercalation kinetics": "coupled ion-electron transfer"} + parameter_values = pybamm.ParameterValues("Marquis2019") + parameter_values.update( + { + "Negative electrode reorganization energy [eV]": 0.2, + "Positive electrode reorganization energy [eV]": 0.2, + "Negative electrode coupled ion-electron transfer rate constant " + "[A.m-2]": 10, + "Positive electrode coupled ion-electron transfer rate constant " + "[A.m-2]": 10, + }, + check_already_exists=False, + ) + self.run_basic_processing_test(options, parameter_values=parameter_values) + def test_kinetics_marcus(self): options = {"intercalation kinetics": "Marcus"} parameter_values = pybamm.ParameterValues("Marquis2019") diff --git a/packages/pybamm/tests/unit/test_models/test_full_battery_models/test_base_battery_model.py b/packages/pybamm/tests/unit/test_models/test_full_battery_models/test_base_battery_model.py index 1d573d95bf..e76c39b290 100644 --- a/packages/pybamm/tests/unit/test_models/test_full_battery_models/test_base_battery_model.py +++ b/packages/pybamm/tests/unit/test_models/test_full_battery_models/test_base_battery_model.py @@ -30,7 +30,7 @@ 'exchange-current density': 'single' (possible: ['single', 'current sigmoid']) 'heat of mixing': 'false' (possible: ['false', 'true']) 'hydrolysis': 'false' (possible: ['false', 'true']) -'intercalation kinetics': 'symmetric Butler-Volmer' (possible: ['symmetric Butler-Volmer', 'asymmetric Butler-Volmer', 'linear', 'Marcus', 'Marcus-Hush-Chidsey', 'MSMR']) +'intercalation kinetics': 'symmetric Butler-Volmer' (possible: ['symmetric Butler-Volmer', 'asymmetric Butler-Volmer', 'linear', 'Marcus', 'Marcus-Hush-Chidsey', 'coupled ion-electron transfer', 'MSMR']) 'interface utilisation': 'full' (possible: ['full', 'constant', 'current-driven']) 'lithium plating': 'none' (possible: ['none', 'reversible', 'partially reversible', 'irreversible']) 'lithium plating porosity change': 'false' (possible: ['false', 'true']) diff --git a/packages/pybamm/tests/unit/test_models/test_full_battery_models/test_lithium_ion/base_lithium_ion_half_cell_tests.py b/packages/pybamm/tests/unit/test_models/test_full_battery_models/test_lithium_ion/base_lithium_ion_half_cell_tests.py index e8fb3d9eb6..f70b28a85a 100644 --- a/packages/pybamm/tests/unit/test_models/test_full_battery_models/test_lithium_ion/base_lithium_ion_half_cell_tests.py +++ b/packages/pybamm/tests/unit/test_models/test_full_battery_models/test_lithium_ion/base_lithium_ion_half_cell_tests.py @@ -50,6 +50,21 @@ def test_well_posed_kinetics_mhc(self): options = {"intercalation kinetics": "Marcus-Hush-Chidsey"} self.check_well_posedness(options) + def test_well_posed_kinetics_coupled_ion_electron_transfer(self): + options = { + "intercalation kinetics": ( + "Marcus-Hush-Chidsey", + "coupled ion-electron transfer", + ) + } + self.check_well_posedness(options) + + def test_coupled_ion_electron_transfer_lithium_metal_not_supported(self): + options = {"intercalation kinetics": "coupled ion-electron transfer"} + if self.model is not None: + with pytest.raises(pybamm.OptionError, match=r"lithium metal electrode"): + self.check_well_posedness(options) + def test_well_posed_constant_sei(self): options = {"SEI": "constant"} self.check_well_posedness(options) diff --git a/packages/pybamm/tests/unit/test_models/test_full_battery_models/test_lithium_ion/base_lithium_ion_tests.py b/packages/pybamm/tests/unit/test_models/test_full_battery_models/test_lithium_ion/base_lithium_ion_tests.py index 965f5b8949..3d86628111 100644 --- a/packages/pybamm/tests/unit/test_models/test_full_battery_models/test_lithium_ion/base_lithium_ion_tests.py +++ b/packages/pybamm/tests/unit/test_models/test_full_battery_models/test_lithium_ion/base_lithium_ion_tests.py @@ -270,6 +270,10 @@ def test_well_posed_kinetics_mhc(self): options = {"intercalation kinetics": "Marcus-Hush-Chidsey"} self.check_well_posedness(options) + def test_well_posed_kinetics_coupled_ion_electron_transfer(self): + options = {"intercalation kinetics": "coupled ion-electron transfer"} + self.check_well_posedness(options) + def test_well_posed_sei_constant(self): options = {"SEI": "constant"} self.check_well_posedness(options) diff --git a/packages/pybamm/tests/unit/test_models/test_submodels/test_interface/test_kinetics/test_coupled_ion_electron_transfer.py b/packages/pybamm/tests/unit/test_models/test_submodels/test_interface/test_kinetics/test_coupled_ion_electron_transfer.py new file mode 100644 index 0000000000..dc2e4d6d23 --- /dev/null +++ b/packages/pybamm/tests/unit/test_models/test_submodels/test_interface/test_kinetics/test_coupled_ion_electron_transfer.py @@ -0,0 +1,155 @@ +import numpy as np +import pytest +from scipy.integrate import quad + +import pybamm + +C_MAX = 30000.0 +C_E_INIT = 1000.0 +T = 298.15 +REORGANIZATION_ENERGY = 0.2 +RATE_CONSTANT = 10.0 + + +def mhc_rate(eta): + """Marcus-Hush-Chidsey rate integral, with eta and lambda in units of kT/e.""" + F_RT = pybamm.constants.F.value / (pybamm.constants.R.value * T) + lambda_T = F_RT * REORGANIZATION_ENERGY + return quad( + lambda x: np.exp( + -((x - lambda_T + eta) ** 2) / (4 * lambda_T) - np.logaddexp(0, x) + ), + -np.inf, + np.inf, + )[0] + + +def expected_current(eta_r, stoichiometry, relative_c_e): + F_RT = pybamm.constants.F.value / (pybamm.constants.R.value * T) + eta_f = F_RT * eta_r + np.log(relative_c_e * (1 - stoichiometry) / stoichiometry) + oxidation = stoichiometry * mhc_rate(eta_f) + reduction = relative_c_e * (1 - stoichiometry) * mhc_rate(-eta_f) + return RATE_CONSTANT * (oxidation - reduction) + + +class TestCoupledIonElectronTransfer: + def get_current(self, eta_r, stoichiometry, relative_c_e): + param = pybamm.LithiumIonParameters() + options = pybamm.BatteryModelOptions( + {"intercalation kinetics": "coupled ion-electron transfer"} + ) + submodel = pybamm.kinetics.CoupledIonElectronTransfer( + param, "negative", "lithium-ion main", options, "primary" + ) + variables = { + "Negative electrolyte concentration [mol.m-3]": pybamm.Scalar( + relative_c_e * C_E_INIT + ), + "Negative electrode temperature [K]": pybamm.Scalar(T), + "Negative particle surface concentration [mol.m-3]": pybamm.Scalar( + stoichiometry * C_MAX + ), + } + j0 = submodel._get_exchange_current_density(variables) + j = submodel._get_kinetics( + j0, pybamm.Scalar(1), pybamm.Scalar(eta_r), pybamm.Scalar(T), 1 + ) + parameter_values = pybamm.ParameterValues( + { + "Maximum concentration in negative electrode [mol.m-3]": C_MAX, + "Initial concentration in electrolyte [mol.m-3]": C_E_INIT, + "Negative electrode reorganization energy [eV]": ( + REORGANIZATION_ENERGY + ), + "Negative electrode coupled ion-electron transfer rate constant " + "[A.m-2]": RATE_CONSTANT, + } + ) + return ( + parameter_values.process_symbol(j0).evaluate().item(), + parameter_values.process_symbol(j).evaluate().item(), + ) + + @pytest.mark.parametrize("stoichiometry", [0.05, 0.5, 0.95]) + @pytest.mark.parametrize("relative_c_e", [0.5, 1.0, 2.0]) + def test_no_current_at_equilibrium(self, stoichiometry, relative_c_e): + _, j = self.get_current(0.0, stoichiometry, relative_c_e) + assert j == pytest.approx(0, abs=1e-12) + + @pytest.mark.parametrize("eta_r", [-0.4, -0.2, -0.02, 0.02, 0.2, 0.4]) + @pytest.mark.parametrize("stoichiometry", [0.1, 0.6]) + def test_current_matches_marcus_integral(self, eta_r, stoichiometry): + # The Zeng asymptotic form is within 5% of the integral at this lambda + _, j = self.get_current(eta_r, stoichiometry, 1.2) + np.testing.assert_allclose( + j, expected_current(eta_r, stoichiometry, 1.2), rtol=0.05 + ) + assert np.sign(j) == np.sign(eta_r) + + def test_current_saturates(self): + F_RT = pybamm.constants.F.value / (pybamm.constants.R.value * T) + plateau = 2 * RATE_CONSTANT * np.sqrt(np.pi * F_RT * REORGANIZATION_ENERGY) + stoichiometry, relative_c_e = 0.3, 0.8 + _, j_oxidation = self.get_current(5.0, stoichiometry, relative_c_e) + _, j_reduction = self.get_current(-5.0, stoichiometry, relative_c_e) + np.testing.assert_allclose(j_oxidation, plateau * stoichiometry, rtol=1e-6) + np.testing.assert_allclose( + j_reduction, -plateau * relative_c_e * (1 - stoichiometry), rtol=1e-6 + ) + + @pytest.mark.parametrize("stoichiometry", [0.2, 0.8]) + def test_exchange_current_is_linear_response(self, stoichiometry): + eta_r = 1e-6 + j0, j = self.get_current(eta_r, stoichiometry, 1.0) + F_RT = pybamm.constants.F.value / (pybamm.constants.R.value * T) + np.testing.assert_allclose(j, j0 * F_RT * eta_r, rtol=1e-4) + + def test_lead_acid_not_supported(self): + with pytest.raises(pybamm.OptionError, match=r"lithium-ion intercalation"): + pybamm.kinetics.CoupledIonElectronTransfer( + pybamm.LeadAcidParameters(), + "negative", + "lead-acid main", + pybamm.BatteryModelOptions({}), + ) + + def test_current_sigmoid_not_supported(self): + options = pybamm.BatteryModelOptions( + { + "intercalation kinetics": "coupled ion-electron transfer", + "exchange-current density": "current sigmoid", + } + ) + with pytest.raises(pybamm.OptionError, match=r"exchange-current density"): + pybamm.kinetics.CoupledIonElectronTransfer( + pybamm.LithiumIonParameters(), "negative", "lithium-ion main", options + ) + + def test_citations(self): + citations = pybamm.citations + citations._reset() + pybamm.kinetics.CoupledIonElectronTransfer( + pybamm.LithiumIonParameters(), + "negative", + "lithium-ion main", + pybamm.BatteryModelOptions({}), + ) + assert "Fraggedakis2021" in citations._papers_to_cite + assert "Zeng2014" in citations._papers_to_cite + assert "Sripad2020" not in citations._papers_to_cite + + def test_x_averaged_concentrations_broadcast_to_electrode(self): + submodel = pybamm.kinetics.CoupledIonElectronTransfer( + pybamm.LithiumIonParameters(), + "negative", + "lithium-ion main", + pybamm.BatteryModelOptions({}), + ) + x_averaged = pybamm.PrimaryBroadcast( + pybamm.Variable("c", domain="current collector"), "negative particle" + ) + eta_r = pybamm.PrimaryBroadcast( + pybamm.Variable("eta", domain="negative electrode"), "negative particle" + ) + matched = submodel._match_domains(x_averaged, eta_r) + assert matched.secondary_domain == ["negative electrode"]