Electrode reactions as rate reactions, and a reusable half-cell

docs/api/electrochemistry.md

@@ -56,6 +56,9 @@
 ## Half-cells
 
 ```{eval-rst}
+.. autoclass:: RedoxCouple
+   :members:
+
 .. autoclass:: RedoxHalfCell
    :members:
 ```

examples/electrochemistry/vanadium_half_cell_1d.py

@@ -82,7 +82,7 @@
 import sympy.physics.units as syu
 
 from sumfem import discrete, solvers, symbolic
-from sumfem.electrochemistry import Ion, PorousElectrode, RedoxHalfCell
+from sumfem.electrochemistry import Ion, PorousElectrode, RedoxCouple, RedoxHalfCell
 
 # -
 
@@ -185,31 +185,30 @@
 }
 # -
 
-# {py:class}`RedoxHalfCell <sumfem.electrochemistry.RedoxHalfCell>` sums the
-# models of {eq}`half-cell` for a redox couple. The reaction
-# $\mathrm{V^{2+} \rightleftharpoons V^{3+} + e^-}$ of {eq}`half-cell-kinetics`
-# is given by three arguments: the ion `reduced`, V²⁺, the ion `oxidized`,
-# V³⁺, and the number of `electrons` the reaction transfers, which the
-# charges of the ions must match. The kinetics take the standard rate
-# constant $k_0$ and the formal potential $E_0$ of the couple.
-# The other ions take no part in the reaction; they carry current and keep
-# the electrolyte neutral. Discharge oxidizes V²⁺ in this negative
+# A {py:class}`RedoxCouple <sumfem.electrochemistry.RedoxCouple>` gives the
+# reaction $\mathrm{V^{2+} \rightleftharpoons V^{3+} + e^-}$ of
+# {eq}`half-cell-kinetics`: the ion `reduced`, V²⁺, the ion `oxidized`, V³⁺,
+# and the number of `electrons` the reaction transfers, which the charges of
+# the ions must match, with the standard rate constant $k_0$ and the formal
+# potential $E_0$ of its kinetics. Discharge oxidizes V²⁺ in this negative
 # half-cell; the positive half-cell of the battery, with the couple
-# VO²⁺/VO₂⁺, reduces its couple on discharge, which reverses the current. The
-# concentrations of the reservoir are Dirichlet conditions at the membrane:
+# VO²⁺/VO₂⁺, reduces its couple on discharge, which reverses the current:
 
-cell = RedoxHalfCell(
-    domain,
-    ions,
-    reservoir,
+couple = RedoxCouple(
     reduced="V2",
     oxidized="V3",
     electrons=1,
-    electrode=electrode,
-    oxidation_on_discharge=True,
     rate_constant=2.5e-7 * syu.m / syu.second,
     formal_potential=-0.36 * syu.volt,
+    oxidation_on_discharge=True,
 )
+
+# {py:class}`RedoxHalfCell <sumfem.electrochemistry.RedoxHalfCell>` sums the
+# models of {eq}`half-cell` for the couple. The other ions take no part in
+# the reaction; they carry current and keep the electrolyte neutral. The
+# concentrations of the reservoir are Dirichlet conditions at the membrane:
+
+cell = RedoxHalfCell(domain, ions, reservoir, couple, electrode)
 inlet = cell.inlet_conditions("Membrane")
 
 # The discretization bounds the quadrature degree of the nonlinear kinetics,

examples/electrochemistry/vanadium_half_cell_2d.py

@@ -87,7 +87,7 @@
 import ufl
 
 from sumfem import discrete, solvers, symbolic
-from sumfem.electrochemistry import Ion, PorousElectrode, RedoxHalfCell
+from sumfem.electrochemistry import Ion, PorousElectrode, RedoxCouple, RedoxHalfCell
 from sumfem.models import BoundaryFlux, Darcy
 from sumfem.plotting import gather_grid
 
@@ -226,30 +226,34 @@
 }
 # -
 
-# The half-cell is the negative one of {doc}`vanadium_half_cell_1d`. The
-# superficial velocity of the flow carries the ions and drives the mass
+# The half-cell is the negative one of {doc}`vanadium_half_cell_1d`, with its
+# couple. The superficial velocity of the flow carries the ions and drives the mass
 # transfer to the fibres, which a velocity switches on. The flow carries the
 # ions out at the `outlet`, the membrane passes the current as protons with
 # `proton_flux`, and the reservoir concentrations hold at the inlet:
 
-cell = RedoxHalfCell(
-    domain,
-    ions,
-    reservoir,
+# +
+couple = RedoxCouple(
     reduced="V2",
     oxidized="V3",
-    electrons=1,
-    electrode=electrode,
-    oxidation_on_discharge=True,
     rate_constant=2.5e-7 * syu.m / syu.second,
     formal_potential=-0.36 * syu.volt,
+    oxidation_on_discharge=True,
+)
+cell = RedoxHalfCell(
+    domain,
+    ions,
+    reservoir,
+    couple,
+    electrode,
     velocity=darcy.superficial_velocity,
     outlet=("Outlet",),
     proton_flux=True,
     density=density,
     viscosity=viscosity,
 )
 inlet = cell.inlet_conditions("Inlet")
+# -
 
 # The flow and the half-cell are discretized together, so they share their
 # parameters, with the quadrature degree of the nonlinear kinetics bounded as

examples/electrochemistry/vanadium_half_cell_3d.py

@@ -77,7 +77,7 @@
 import ufl
 
 from sumfem import discrete, solvers, symbolic
-from sumfem.electrochemistry import Ion, PorousElectrode, RedoxHalfCell
+from sumfem.electrochemistry import Ion, PorousElectrode, RedoxCouple, RedoxHalfCell
 from sumfem.models import Brinkman, BrinkmanForchheimer, GradDiv, NavierStokes, Stokes
 from sumfem.plotting import gather_grid
 
@@ -334,17 +334,19 @@ def flow_through_electrode(sizes: tuple[float, float]) -> None:
     "H": 1000.0 * concentration,
     "HSO4": 1000.0 * concentration,
 }
-cell = RedoxHalfCell(
-    domain,
-    ions,
-    reservoir,
+couple = RedoxCouple(
     reduced="V2",
     oxidized="V3",
-    electrons=1,
-    electrode=electrode,
-    oxidation_on_discharge=True,
     rate_constant=2.5e-7 * syu.m / syu.second,
     formal_potential=-0.36 * syu.volt,
+    oxidation_on_discharge=True,
+)
+cell = RedoxHalfCell(
+    domain,
+    ions,
+    reservoir,
+    couple,
+    electrode,
     channels=channels,
     velocity=u,
     outlet=("Outlet",),

src/sumfem/electrochemistry/__init__.py

@@ -11,13 +11,14 @@
   on :class:`ElectrodeReaction`, with the Butler-Volmer equation
   :func:`~sumfem.models.butler_volmer`.
 - Half-cells of a redox flow battery: :class:`RedoxHalfCell`, negative or
-  positive, of any couple and electrolyte.
+  positive, of any :class:`RedoxCouple` and electrolyte.
 
 The package imports from the core of :mod:`sumfem`, and never the other
 way around.
 """
 
 from .charge import ElectrolyteConduction
+from .couple import RedoxCouple
 from .electrolyte import Electrolyte, Ion
 from .half_cell import RedoxHalfCell
 from .kinetics import (
@@ -37,6 +38,7 @@
     "Region",
     "ElectrolyteConduction",
     "NernstPlanck",
+    "RedoxCouple",
     "RedoxHalfCell",
     "ButlerVolmer",
     "ElectrodeReaction",

src/sumfem/electrochemistry/couple.py

@@ -0,0 +1,99 @@
+"""Redox couples and their reactions."""
+
+from __future__ import annotations
+
+import dataclasses
+import types
+import typing
+
+import sympy
+
+from .electrolyte import Ion
+
+__all__ = ["RedoxCouple"]
+
+
+@dataclasses.dataclass(frozen=True)
+class RedoxCouple:
+    """A redox couple and its reaction at an electrode.
+
+    Oxidation turns the reduced ion into the oxidized one, as the reaction
+    ``reactants -> products + n e-``. The reactants hold the reduced ion and
+    the products the oxidized one, and either may hold other species of the
+    reaction, e.g. VO2+ + H2O -> VO2+ + 2 H+ + e- consumes water and
+    produces protons. The coefficients are those of
+    :class:`~sumfem.models.RateReaction`. The kinetics are those of
+    :class:`~sumfem.electrochemistry.ButlerVolmer`, in the concentrations of
+    the couple.
+
+    Attributes:
+        reduced: Name of the reduced ion, e.g. ``"V2"``.
+        oxidized: Name of the oxidized ion, e.g. ``"V3"``.
+        rate_constant: Standard rate constant, in m/s or with units.
+        formal_potential: Formal potential, in V or with units.
+        oxidation_on_discharge: Whether discharge oxidizes the couple, as in
+            the negative half-cell of a battery; the positive half-cell
+            reduces its couple on discharge.
+        electrons: Number of electrons ``n`` the reaction transfers.
+        anodic_transfer: Anodic transfer coefficient.
+        cathodic_transfer: Cathodic transfer coefficient.
+        reactants: The coefficient of each ion that oxidation consumes, by
+            name. Defaults to the reduced ion alone, e.g.
+            ``{"VO": 1, "H2O": 1}``.
+        products: The coefficient of each ion that oxidation produces.
+            Defaults to the oxidized ion alone, e.g. ``{"VO2": 1, "H": 2}``.
+
+    Raises:
+        ValueError: If the reactants lack the reduced ion or the products the
+            oxidized one.
+    """
+
+    reduced: str
+    oxidized: str
+    rate_constant: float | sympy.Expr
+    formal_potential: float | sympy.Expr
+    oxidation_on_discharge: bool
+    electrons: int = 1
+    anodic_transfer: float = 0.5
+    cathodic_transfer: float = 0.5
+    reactants: typing.Mapping[str, int] | None = None
+    products: typing.Mapping[str, int] | None = None
+
+    def __post_init__(self) -> None:
+        reactants = dict(self.reactants or {self.reduced: 1})
+        products = dict(self.products or {self.oxidized: 1})
+        if self.reduced not in reactants or self.oxidized not in products:
+            raise ValueError(
+                f"The reactants {reactants} must hold the reduced ion "
+                f"{self.reduced!r}, and the products {products} the oxidized "
+                f"ion {self.oxidized!r}."
+            )
+        object.__setattr__(self, "reactants", types.MappingProxyType(reactants))
+        object.__setattr__(self, "products", types.MappingProxyType(products))
+
+    @property
+    def species(self) -> tuple[str, ...]:
+        """Names of the ions of the reaction: the reactants, then the products."""
+        assert self.reactants is not None and self.products is not None
+        return (*self.reactants, *self.products)
+
+    def check(self, ions: typing.Sequence[Ion]) -> None:
+        """Check the reaction against the ions of an electrolyte.
+
+        Raises:
+            ValueError: If an ion of the reaction is not among the ions, or
+                the reaction does not conserve charge.
+        """
+        assert self.reactants is not None and self.products is not None
+        charges = {ion.name: ion.charge for ion in ions}
+        missing = [ion for ion in self.species if ion not in charges]
+        if missing:
+            raise ValueError(f"The electrolyte has no ions {missing} of the reaction.")
+        produced = sum(nu * charges[ion] for ion, nu in self.products.items())
+        consumed = sum(nu * charges[ion] for ion, nu in self.reactants.items())
+        if produced - consumed != self.electrons:
+            raise ValueError(
+                f"{dict(self.reactants)} -> {dict(self.products)} + "
+                f"{self.electrons} e- does not conserve charge: the charge "
+                f"goes from {consumed} to {produced}."
+            )