diff --git a/docs/api/electrochemistry.md b/docs/api/electrochemistry.md index 5021441..99faf0f 100644 --- a/docs/api/electrochemistry.md +++ b/docs/api/electrochemistry.md @@ -48,6 +48,9 @@ .. autoclass:: MassTransferButlerVolmer :members: + +.. autoclass:: PistonVelocityButlerVolmer + :members: ``` ## Half-cells diff --git a/src/sumfem/electrochemistry/__init__.py b/src/sumfem/electrochemistry/__init__.py index 51853bb..9ad5a7e 100644 --- a/src/sumfem/electrochemistry/__init__.py +++ b/src/sumfem/electrochemistry/__init__.py @@ -6,9 +6,10 @@ a cell filled with electrolyte. - Charge and transport in the electrolyte: :class:`ElectrolyteConduction` and :class:`NernstPlanck`. -- Kinetics: :class:`ExchangeCurrentButlerVolmer`, :class:`ButlerVolmer` - and :class:`MassTransferButlerVolmer`, on :class:`ElectrodeReaction`, - with the Butler-Volmer equation :func:`~sumfem.models.butler_volmer`. +- Kinetics: :class:`ExchangeCurrentButlerVolmer`, :class:`ButlerVolmer`, + :class:`MassTransferButlerVolmer` and :class:`PistonVelocityButlerVolmer`, + 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. @@ -24,6 +25,7 @@ ElectrodeReaction, ExchangeCurrentButlerVolmer, MassTransferButlerVolmer, + PistonVelocityButlerVolmer, ) from .porous import PorousElectrode, Region from .transport import NernstPlanck @@ -40,4 +42,5 @@ "ElectrodeReaction", "ExchangeCurrentButlerVolmer", "MassTransferButlerVolmer", + "PistonVelocityButlerVolmer", ] diff --git a/src/sumfem/electrochemistry/kinetics.py b/src/sumfem/electrochemistry/kinetics.py index ed09660..268fc3c 100644 --- a/src/sumfem/electrochemistry/kinetics.py +++ b/src/sumfem/electrochemistry/kinetics.py @@ -7,6 +7,8 @@ from __future__ import annotations +import typing + import numpy.typing as npt import sympy.physics.units as syu import ufl @@ -21,6 +23,7 @@ "ElectrodeReaction", "ExchangeCurrentButlerVolmer", "MassTransferButlerVolmer", + "PistonVelocityButlerVolmer", ] INPUT = ufl.core.expr.Expr | npt.ArrayLike | None @@ -33,8 +36,9 @@ class ElectrodeReaction(Residual): The volumetric current ``f = a i``, with the specific area ``a`` and the interfacial current density ``i`` of :meth:`current_density`, leaves the solid and enters the electrolyte. It consumes the reduced species and - produces the oxidized one at the rate ``f / (n F)``. Species without a - field get no source. Subclasses give the kinetics. + produces the oxidized one at the rate ``f / (n F)``, and every other + species of the reaction, e.g. protons or water, ``nu`` times as much. + Species without a field get no source. Subclasses give the kinetics. Args: name: Name of the model. @@ -47,6 +51,10 @@ class ElectrodeReaction(Residual): ``".specific_area"`` in 1/m. electrons: Number of electrons ``n`` transferred. markers: Cell markers of the electrode. Empty for the whole domain. + others: Other species of the reaction: concentrations with the moles + ``nu`` of each produced per mole oxidized, negative for one + consumed. E.g. VO2+ + H2O -> VO2+ + 2 H+ + e- produces two protons + and consumes one water. """ _solid_potential: Field @@ -73,6 +81,10 @@ class ElectrodeReaction(Residual): _markers: tuple[str, ...] """Cell markers of the electrode.""" + _others: tuple[tuple[ufl.core.expr.Expr, float], ...] + """Other species of the reaction, with the moles produced per mole + oxidized.""" + def __init__( self, name: str, @@ -84,6 +96,7 @@ def __init__( specific_area: INPUT = None, electrons: int = 1, markers: tuple[str, ...] = (), + others: typing.Sequence[tuple[ufl.core.expr.Expr, float]] = (), ): super().__init__(name) self._solid_potential = solid_potential @@ -100,6 +113,7 @@ def __init__( ) self._electrons = electrons self._markers = tuple(markers) + self._others = tuple((c, float(nu)) for c, nu in others) @property def electrolyte(self) -> Electrolyte: @@ -126,13 +140,15 @@ def residual(self) -> ufl.Form: dx = self._solid_potential.domain.dx(*self._markers) F = f * self._solid_potential.v * dx - f * self._electrolyte_potential.v * dx rate = f / (self._electrons * self._electrolyte.F) + # The residual of a species gains what the reaction consumes species = ( (self._reduced_concentration, 1.0), (self._oxidized_concentration, -1.0), + *((c, -nu) for c, nu in self._others), ) - for c, sign in species: + for c, consumed in species: if isinstance(c, Field): - F += sign * rate * c.v * dx + F += consumed * rate * c.v * dx return F @@ -290,6 +306,13 @@ class ButlerVolmer(ElectrodeReaction): Defaults to the parameter ``".cathodic_transfer"``, 0.5. electrons: See :class:`ElectrodeReaction`. markers: See :class:`ElectrodeReaction`. + stoichiometry: Other ions of the reaction, by name, with the moles + produced per mole of the reduced ion oxidized, see ``others`` of + :class:`ElectrodeReaction`. + + Raises: + ValueError: If an ion of the stoichiometry has no field, as the ion + that follows from electroneutrality. """ _reduced: Ion @@ -325,7 +348,15 @@ def __init__( cathodic_transfer: INPUT = None, electrons: int = 1, markers: tuple[str, ...] = (), + stoichiometry: typing.Mapping[str, float] | None = None, ): + stoichiometry = dict(stoichiometry or {}) + for ion in stoichiometry: + if electrolyte[ion].concentration is None: + raise ValueError( + f"Ion {ion!r} of the reaction follows from electroneutrality, " + "so it has no field to produce or consume." + ) super().__init__( name, solid_potential, @@ -336,6 +367,7 @@ def __init__( specific_area=specific_area, electrons=electrons, markers=markers, + others=[(electrolyte[ion].field, nu) for ion, nu in stoichiometry.items()], ) domain = solid_potential.domain self._reduced = electrolyte[reduced] @@ -522,6 +554,7 @@ def __init__( cathodic_transfer: INPUT = None, electrons: int = 1, markers: tuple[str, ...] = (), + stoichiometry: typing.Mapping[str, float] | None = None, ): super().__init__( name, @@ -537,6 +570,7 @@ def __init__( cathodic_transfer=cathodic_transfer, electrons=electrons, markers=markers, + stoichiometry=stoichiometry, ) domain = solid_potential.domain self._velocity = as_expr(velocity) @@ -609,12 +643,115 @@ def surface_concentrations( self, ) -> tuple[ufl.core.expr.Expr, ufl.core.expr.Expr]: """Surface concentrations balancing reaction and mass transfer.""" - k_a, k_c = self.rate_constants() m_red, m_ox = ( m + self._transfer_floor for m in self.mass_transfer_coefficients() ) - c_red, c_ox = self._bulk() - denominator = 1.0 + k_a / m_red + k_c / m_ox - surface_ox = ((1.0 + k_a / m_red) * c_ox + k_a / m_ox * c_red) / denominator - surface_red = ((1.0 + k_c / m_ox) * c_red + k_c / m_red * c_ox) / denominator - return surface_red, surface_ox + return _balanced_surface(*self.rate_constants(), m_red, m_ox, *self._bulk()) + + +def _balanced_surface( + k_a: ufl.core.expr.Expr, + k_c: ufl.core.expr.Expr, + m_red: ufl.core.expr.Expr, + m_ox: ufl.core.expr.Expr, + c_red: ufl.core.expr.Expr, + c_ox: ufl.core.expr.Expr, +) -> tuple[ufl.core.expr.Expr, ufl.core.expr.Expr]: + """Surface concentrations where mass transfer balances the reaction. + + The solution of ``m_red (c_red - s_red) = k_a s_red - k_c s_ox`` and + ``m_ox (c_ox - s_ox) = -(k_a s_red - k_c s_ox)``, see equation 20 of + `Aghabarari et al., 2026 `_. + """ + denominator = 1.0 + k_a / m_red + k_c / m_ox + surface_ox = ((1.0 + k_a / m_red) * c_ox + k_a / m_ox * c_red) / denominator + surface_red = ((1.0 + k_c / m_ox) * c_red + k_c / m_red * c_ox) / denominator + return surface_red, surface_ox + + +class PistonVelocityButlerVolmer(ButlerVolmer): + """Butler-Volmer kinetics limited by diffusion between the fibres. + + The surface concentrations balance the reaction with mass transfer + ``gamma (c - c_s)`` from the bulk, with the piston velocity + ``gamma = D / delta`` of each ion, its diffusivity in free electrolyte + over the mean distance ``delta`` between the fibres, see equations 23 + and 24 of `Shah et al., 2008 + `_. Unlike + :class:`MassTransferButlerVolmer`, the transfer does not depend on the + flow. + + Args: + name: Name of the model. + solid_potential: See :class:`ButlerVolmer`. + electrolyte_potential: See :class:`ButlerVolmer`. + electrolyte: See :class:`ButlerVolmer`. + reduced: See :class:`ButlerVolmer`. + oxidized: See :class:`ButlerVolmer`. + fiber_distance: Mean distance ``delta`` between the fibres. + Defaults to the parameter ``".fiber_distance"`` in m. + transfer_floor: Regularization added to ``gamma``, in m/s. + **kwargs: The other arguments of :class:`ButlerVolmer`. + """ + + _fiber_distance: ufl.core.expr.Expr + """Mean distance between the fibres.""" + + _transfer_floor: ufl.core.expr.Expr + """Regularization of the piston velocities.""" + + def __init__( + self, + name: str, + solid_potential: Field, + electrolyte_potential: Field, + electrolyte: Electrolyte, + reduced: str, + oxidized: str, + fiber_distance: INPUT = None, + transfer_floor: float = 1e-14, + **kwargs: typing.Any, + ): + super().__init__( + name, + solid_potential, + electrolyte_potential, + electrolyte, + reduced, + oxidized, + **kwargs, + ) + domain = solid_potential.domain + self._fiber_distance = self.input( + "fiber_distance", + fiber_distance, + domain=domain, + unit=syu.m, + doc="Mean distance between the fibres", + ) + self._transfer_floor = self.parameter( + "transfer_floor", + domain, + default=transfer_floor, + unit=syu.m / syu.second, + doc="Regularization of the piston velocities", + ) + + def mass_transfer_coefficients( + self, + ) -> tuple[ufl.core.expr.Expr, ufl.core.expr.Expr]: + """The piston velocities ``D / delta`` of the reduced and oxidized ion.""" + e = self._electrolyte + return ( + e.diffusivity(self._reduced.name) / self._fiber_distance, + e.diffusivity(self._oxidized.name) / self._fiber_distance, + ) + + def surface_concentrations( + self, + ) -> tuple[ufl.core.expr.Expr, ufl.core.expr.Expr]: + """Surface concentrations balancing reaction and mass transfer.""" + m_red, m_ox = ( + m + self._transfer_floor for m in self.mass_transfer_coefficients() + ) + return _balanced_surface(*self.rate_constants(), m_red, m_ox, *self._bulk()) diff --git a/tests/electrochemistry/test_electrochemistry_from_ufl.py b/tests/electrochemistry/test_electrochemistry_from_ufl.py index 3f8ae20..22effa5 100644 --- a/tests/electrochemistry/test_electrochemistry_from_ufl.py +++ b/tests/electrochemistry/test_electrochemistry_from_ufl.py @@ -19,9 +19,10 @@ Ion, MassTransferButlerVolmer, NernstPlanck, + PistonVelocityButlerVolmer, PorousElectrode, ) -from sumfem.symbolic import Domain, Field, UFLResidual +from sumfem.symbolic import Domain, Field, Residual, UFLResidual CELLS = ["interval", "triangle", "quadrilateral", "tetrahedron", "hexahedron"] """The cell types each comparison runs on.""" @@ -239,6 +240,116 @@ def test_kinetics_of_a_redox_couple(cell, mass_transfer, box, compare): ) +@pytest.mark.parametrize("cell", CELLS) +def test_positive_couple_with_protons_water_and_piston_velocities(cell, box, compare): + """VO2+ + H2O <-> VO2+ + 2 H+ + e- of Shah et al. (2008), limited by the + piston velocities ``D / delta`` of each vanadium ion. + + The reaction produces two protons and consumes one water molecule per + vanadium ion oxidized. + """ + gdim = ufl.Cell(cell).topological_dimension + mesh = ufl.Mesh(element("Lagrange", cell, 1, shape=(gdim,))) + V = ufl.FunctionSpace(mesh, element("Lagrange", cell, 1)) + W = ufl.MixedFunctionSpace(V, V, V, V, V, V) + phi_s, phi_l, c4, c5, ch, cw = (ufl.Coefficient(V) for _ in range(6)) + v_s, v_l, w4, w5, wh, ww = ufl.TestFunctions(W) + F, R, T, floor, a, k0, E0 = (ufl.Constant(mesh) for _ in range(7)) + delta, D4, D5 = (ufl.Constant(mesh) for _ in range(3)) + factor = F / (R * T) * (phi_s - phi_l - E0) + k_a, k_c = k0 * ufl.exp(0.5 * factor), k0 * ufl.exp(-0.5 * factor) + c_red, c_ox = positive_part(c4, floor), positive_part(c5, floor) + m_red = D4 / delta + 1e-14 + m_ox = D5 / delta + 1e-14 + denominator = 1.0 + k_a / m_red + k_c / m_ox + c_ox, c_red = ( + ((1.0 + k_a / m_red) * c_ox + k_a / m_ox * c_red) / denominator, + ((1.0 + k_c / m_ox) * c_red + k_c / m_red * c_ox) / denominator, + ) + f = a * F * (k_a * c_red - k_c * c_ox) + rate = f / F + # A mass term for the protons and water, which only the reaction would + # otherwise reach, so that each has a Jacobian block of its own + residual = ( + f * v_s - f * v_l + rate * w4 - rate * w5 - 2 * rate * wh + rate * ww + ) * ufl.dx + (ch * wh + cw * ww) * ufl.dx + names = { + DOMAIN: mesh, + **dict(phi_s=phi_s, phi_l=phi_l, c4=c4, c5=c5, ch=ch, cw=cw), + **dict(F=F, R=R, T=T, floor=floor, a=a, k0=k0, E0=E0), + **dict(delta=delta, D4=D4, D5=D5), + } + unknowns = [phi_s, phi_l, c4, c5, ch, cw] + plain = UFLResidual("plain", residual, unknowns, names=names) + + domain = Domain(DOMAIN, cell) + potentials = [Field(name, domain, ("Lagrange", 1)) for name in ("phi_s", "phi_l")] + fields = [Field(name, domain, ("Lagrange", 1)) for name in ("c4", "c5", "ch", "cw")] + ions = [ + Ion("VO", 2, 3.9e-10, fields[0], positive=True), + Ion("VO2", 1, 3.9e-10, fields[1], positive=True), + Ion("H", 1, 9.3e-9, fields[2]), + Ion("H2O", 0, 2.3e-9, fields[3]), + Ion("SO4", -2, 1.1e-9), + ] + electrolyte = Electrolyte("e", domain, ions, concentration_floor=0.1) + model = PistonVelocityButlerVolmer( + "bv", + *potentials, + electrolyte, + "VO", + "VO2", + fiber_distance=1e-5, + specific_area=1e3, + rate_constant=1e-6, + formal_potential=1.0, + stoichiometry={"H": 2, "H2O": -1}, + ) + + class Mass(Residual): + def residual(self): + ch, cw = fields[2], fields[3] + return (ch * ch.v + cw * cw.v) * domain.dx() + + model = model + Mass("mass") + values = { + "F": electrolyte.F.default_si, + "R": electrolyte.R.default_si, + "T": electrolyte.T.default_si, + "floor": 0.1, + "a": 1e3, + "k0": 1e-6, + "E0": 1.0, + "delta": 1e-5, + "D4": 3.9e-10, + "D5": 3.9e-10, + } + compare( + plain, + plain.unknowns, + model, + [*potentials, *fields], + box(cell), + [(constants(values), {})], + ) + + +def test_the_stoichiometry_needs_ions_with_fields(): + domain = Domain(DOMAIN, "interval") + potentials = [Field(name, domain, ("Lagrange", 1)) for name in ("phi_s", "phi_l")] + c = [Field(name, domain, ("Lagrange", 1)) for name in ("c4", "c5")] + ions = [ + Ion("VO", 2, 3.9e-10, c[0]), + Ion("VO2", 1, 3.9e-10, c[1]), + Ion("H", 1, 9.3e-9), + ] + electrolyte = Electrolyte("e", domain, ions) + with pytest.raises(ValueError, match="electroneutrality"): + ButlerVolmer( + "bv", *potentials, electrolyte, "VO", "VO2", stoichiometry={"H": 2} + ) + + def test_mass_transfer_takes_the_properties_of_the_electrolyte(): """Without a density and viscosity of its own, the mass transfer uses the electrolyte's, so models of one electrolyte share them."""