From 1a2cbaa5e6d0ff87e49660358bc0e05509fddb43 Mon Sep 17 00:00:00 2001 From: Jorgen Schartum Dokken Date: Wed, 11 Nov 2026 23:23:00 +0000 Subject: [PATCH] Massive refactor of half-cell to make it easier to work with. Adds additional features, e.g. fiber distance,. Moves the definition of hte redox reaction out of the core part of the library. --- docs/api/electrochemistry.md | 3 + .../electrochemistry/vanadium_half_cell_1d.py | 35 +- .../electrochemistry/vanadium_half_cell_2d.py | 24 +- .../electrochemistry/vanadium_half_cell_3d.py | 18 +- src/sumfem/electrochemistry/__init__.py | 4 +- src/sumfem/electrochemistry/couple.py | 99 ++++ src/sumfem/electrochemistry/half_cell.py | 469 +++++++++++------- src/sumfem/electrochemistry/kinetics.py | 138 +++--- src/sumfem/models/reactions.py | 12 +- tests/electrochemistry/half_cell/conftest.py | 23 +- .../electrochemistry/half_cell/test_models.py | 214 ++++++-- .../test_electrochemistry_from_ufl.py | 14 +- 12 files changed, 736 insertions(+), 317 deletions(-) create mode 100644 src/sumfem/electrochemistry/couple.py diff --git a/docs/api/electrochemistry.md b/docs/api/electrochemistry.md index 99faf0f..d9e3d9a 100644 --- a/docs/api/electrochemistry.md +++ b/docs/api/electrochemistry.md @@ -56,6 +56,9 @@ ## Half-cells ```{eval-rst} +.. autoclass:: RedoxCouple + :members: + .. autoclass:: RedoxHalfCell :members: ``` diff --git a/examples/electrochemistry/vanadium_half_cell_1d.py b/examples/electrochemistry/vanadium_half_cell_1d.py index 9eddc8d..c20f353 100644 --- a/examples/electrochemistry/vanadium_half_cell_1d.py +++ b/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 ` 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 ` 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 ` 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, diff --git a/examples/electrochemistry/vanadium_half_cell_2d.py b/examples/electrochemistry/vanadium_half_cell_2d.py index 95ae2bd..93d7be8 100644 --- a/examples/electrochemistry/vanadium_half_cell_2d.py +++ b/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,23 +226,26 @@ } # - -# 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, @@ -250,6 +253,7 @@ 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 diff --git a/examples/electrochemistry/vanadium_half_cell_3d.py b/examples/electrochemistry/vanadium_half_cell_3d.py index ef79e3c..0a3b99a 100644 --- a/examples/electrochemistry/vanadium_half_cell_3d.py +++ b/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",), diff --git a/src/sumfem/electrochemistry/__init__.py b/src/sumfem/electrochemistry/__init__.py index 9ad5a7e..4850bcb 100644 --- a/src/sumfem/electrochemistry/__init__.py +++ b/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", diff --git a/src/sumfem/electrochemistry/couple.py b/src/sumfem/electrochemistry/couple.py new file mode 100644 index 0000000..eb93950 --- /dev/null +++ b/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}." + ) diff --git a/src/sumfem/electrochemistry/half_cell.py b/src/sumfem/electrochemistry/half_cell.py index 5e42a80..35afe7f 100644 --- a/src/sumfem/electrochemistry/half_cell.py +++ b/src/sumfem/electrochemistry/half_cell.py @@ -29,8 +29,13 @@ ) from .charge import ElectrolyteConduction +from .couple import RedoxCouple from .electrolyte import CONCENTRATION, Electrolyte, Ion -from .kinetics import ButlerVolmer, MassTransferButlerVolmer +from .kinetics import ( + ButlerVolmer, + MassTransferButlerVolmer, + PistonVelocityButlerVolmer, +) from .porous import PorousElectrode, Region from .transport import NernstPlanck @@ -47,36 +52,38 @@ class RedoxHalfCell: The unknowns are the potentials ``phi_s`` of the solid and ``phi_l`` of the electrolyte, in V, and the concentration ``c_`` of every ion with a reservoir concentration, in mol/m^3. The one ion without follows - from electroneutrality. In the porous electrode, the couple - ``reduced <-> oxidized + n e-`` follows Butler-Volmer kinetics, limited - by mass transfer to the fibres when there is flow. In free channels the - electrolyte is not corrected for porosity, and there is no reaction or - solid phase. + from electroneutrality. In the porous electrode, the reaction of the + :class:`~sumfem.electrochemistry.RedoxCouple` follows Butler-Volmer + kinetics, limited by mass transfer to the fibres when there is flow or a + fibre distance. In free channels the electrolyte is not corrected for + porosity, and there is no reaction or solid phase. The model, :attr:`model`, is the sum of: - - ``"solid"``: Ohmic conduction in the solid, + - ``"_solid"``: Ohmic conduction in the solid, :class:`~sumfem.models.Conduction` with Bruggeman's effective conductivity of the electrode; - - ``"reaction"``: the kinetics of the couple, - :class:`~sumfem.electrochemistry.ButlerVolmer`, or + - ``"_reaction"``: the kinetics of the couple, + :class:`~sumfem.electrochemistry.ButlerVolmer`, :class:`~sumfem.electrochemistry.MassTransferButlerVolmer` with mass - transfer, which moves the current between the phases and consumes and - produces the ions of the couple; - - ``"collector_current"`` and ``"membrane_current"``: the current of the - cell into the solid and out of the electrolyte, - :class:`~sumfem.models.BoundaryFlux`; + transfer, or :class:`~sumfem.electrochemistry.PistonVelocityButlerVolmer` + with a fibre distance, which moves the current between the phases and + consumes and produces the ions of the reaction; + - with collector and membrane markers, ``"_collector_current"`` + and ``"_membrane_current"``: the current of the cell into the + solid and out of the electrolyte, :class:`~sumfem.models.BoundaryFlux`; - for the electrode's pores and, with a suffix ``"_channels"``, the free channels: charge conservation in the electrolyte, :class:`~sumfem.electrochemistry.ElectrolyteConduction`, and for every transported ion its storage, :class:`~sumfem.models.Accumulation` weighted with the porosity, and its transport, :class:`~sumfem.electrochemistry.NernstPlanck`, named - ``"storage_"`` and ``"transport_"``; - - with an outlet, ``"outflow_"``: advection out of the cell, + ``"_storage_"`` and ``"_transport_"``; + - with an outlet, ``"_outflow_"``: advection out of the cell, :class:`~sumfem.models.Outflow`; - - with a proton flux, ``"proton_flux"``: the membrane current carried by - protons, a :class:`~sumfem.models.BoundaryFlux` of the ion ``"H"``. + - with a proton flux, ``"_proton_flux"``: the membrane current + carried by protons, a :class:`~sumfem.models.BoundaryFlux` of the ion + ``"H"``. The cell is galvanostatic: the parameter ``".current"`` is the current density through the membrane, in A/m^2, positive on discharge. @@ -94,30 +101,24 @@ class RedoxHalfCell: reservoir: The concentration of each transported ion in the reservoir, in mol/m^3 or a quantity with units, the default of the parameter ``".c0_"``. - reduced: Name of the reduced ion of the couple. With ``oxidized`` - and ``electrons`` it gives the reaction - ``oxidized + n e- <-> reduced``, e.g. V3 + e- <-> V2. - oxidized: Name of the oxidized ion of the couple, of ``electrons`` - more charge than the reduced one. + couple: The redox couple and its reaction. Its parameters are those + of the cell: ``".rate_constant"``, + ``".formal_potential"``, ``".anodic_transfer"`` and + ``".cathodic_transfer"``. electrode: The porous electrode. - oxidation_on_discharge: Whether discharge oxidizes the couple, as in - the negative half-cell. - rate_constant: Standard rate constant of the couple in m/s or with - units, the default of the parameter ``".rate_constant"``. - formal_potential: Formal potential of the couple in V or with - units, the default of the parameter ``".formal_potential"``. - anodic_transfer: Anodic transfer coefficient, the default of the - parameter ``".anodic_transfer"``. - cathodic_transfer: Cathodic transfer coefficient, the default of the - parameter ``".cathodic_transfer"``. - electrons: Number of electrons the couple transfers. channels: Cell markers of free channels. Empty for none. velocity: Superficial velocity of the electrolyte, or ``None`` without flow. mass_transfer: Limit the kinetics by mass transfer to the fibres, which needs a velocity. Defaults to whether there is one. - collector: Facet markers of the current collector. - membrane: Facet markers of the membrane. + fiber_distance: Mean distance between the fibres, in m or with + units. Limits the kinetics by the piston velocities of + :class:`~sumfem.electrochemistry.PistonVelocityButlerVolmer` + instead, with or without flow. + collector: Facet markers of the current collector, or ``None`` for + a cell whose collector is resolved, e.g. in a full cell. + membrane: Facet markers of the membrane, or ``None`` for a cell + whose membrane is resolved. outlet: Facet markers of an outlet, through which ions leave by advection. Empty for none. proton_flux: Carry the membrane current by protons, the ion ``"H"``, @@ -128,19 +129,36 @@ class RedoxHalfCell: viscosity: Dynamic viscosity of the electrolyte, given like the density. degree: Degree of the Lagrange elements of the unknowns. - name: Prefix of the cell's own parameters. + name: Prefix of the cell's own parameters, ``"."``, + and of the names of its models, ``"_"``. + phi_s: The potential of the solid, if it exists, e.g. over the + electrode and its current collector. Created on the domain + otherwise. + phi_l: The potential of the electrolyte, if it exists, e.g. over the + electrodes and the membrane of a full cell. + fields: Existing fields of transported ions, by name, e.g. of water + shared with a membrane. The others are created on the domain. + electrolyte_name: Prefix of the parameters of the electrolyte, see + :class:`~sumfem.electrochemistry.Electrolyte`. Raises: ValueError: For mass transfer or an outlet without a velocity, a - reservoir concentration of an ion that the electrolyte does not - have, a reservoir that lacks not exactly one ion, an ion of the - couple without one, a couple that does not conserve charge, or a - balanced ion without charge. + proton flux without membrane markers, a reservoir concentration + of an ion that the electrolyte does not have, a reservoir that + lacks not exactly one ion, an ion of the reaction without one, a + reaction that does not conserve charge, or a balanced ion + without charge. """ _domain: Domain """The domain of the cell.""" + _name: str + """Prefix of the cell's parameters and models.""" + + _couple: RedoxCouple + """The redox couple and its reaction.""" + _electrode: PorousElectrode """The porous electrode.""" @@ -151,12 +169,6 @@ class RedoxHalfCell: _electrolyte: Electrolyte """The electrolyte.""" - _reduced: str - """Name of the reduced ion of the couple.""" - - _oxidized: str - """Name of the oxidized ion of the couple.""" - _phi_s: Field """Potential of the solid.""" @@ -184,11 +196,11 @@ class RedoxHalfCell: _model: ResidualSum """The sum of all models of the cell.""" - _collector: tuple[str, ...] - """Facet markers of the current collector.""" + _collector: tuple[str, ...] | None + """Facet markers of the current collector, if not resolved.""" - _membrane: tuple[str, ...] - """Facet markers of the membrane.""" + _membrane: tuple[str, ...] | None + """Facet markers of the membrane, if not resolved.""" _functionals: types.MappingProxyType[str, ufl.Form] """The integrals behind the cell's diagnostics, created once so that a @@ -199,90 +211,107 @@ def __init__( domain: Domain, ions: typing.Sequence[Ion], reservoir: typing.Mapping[str, float | sympy.Expr], - reduced: str, - oxidized: str, + couple: RedoxCouple, electrode: PorousElectrode, - oxidation_on_discharge: bool, - rate_constant: float | sympy.Expr, - formal_potential: float | sympy.Expr, - anodic_transfer: float = 0.5, - cathodic_transfer: float = 0.5, - electrons: int = 1, channels: tuple[str, ...] = (), velocity: ufl.core.expr.Expr | None = None, mass_transfer: bool | None = None, - collector: tuple[str, ...] = ("CurrentCollector",), - membrane: tuple[str, ...] = ("Membrane",), + fiber_distance: float | sympy.Expr | None = None, + collector: tuple[str, ...] | None = ("CurrentCollector",), + membrane: tuple[str, ...] | None = ("Membrane",), outlet: tuple[str, ...] = (), proton_flux: bool = False, density: ufl.core.expr.Expr | float = 1399.8, viscosity: ufl.core.expr.Expr | float = 6.963e-3, degree: int = 1, name: str = "cell", + phi_s: Field | None = None, + phi_l: Field | None = None, + fields: typing.Mapping[str, Field] | None = None, + electrolyte_name: str = "electrolyte", ): # Mass transfer to the fibres depends on the speed of the flow, so it - # is on by default whenever there is one - mass_transfer = velocity is not None if mass_transfer is None else mass_transfer + # is on by default whenever there is one; piston velocities take the + # place of the flow's correlation + if mass_transfer is None: + mass_transfer = velocity is not None + if fiber_distance is not None: + mass_transfer = False if velocity is None and (mass_transfer or outlet): raise ValueError("Mass transfer and an outlet need a velocity.") - unknown = sorted(set(reservoir) - {ion.name for ion in ions}) - if unknown: - raise ValueError(f"The electrolyte has no ions {unknown}.") - # The one ion without a reservoir concentration has no field, and - # follows from electroneutrality, which closes the equations - balanced = [ion.name for ion in ions if ion.name not in reservoir] - if len(balanced) != 1: - raise ValueError( - "The reservoir must give the concentrations of all ions but one, " - "whose concentration follows from electroneutrality; it lacks " - f"{balanced}." - ) - # The couple is the reaction oxidized + n e- <-> reduced, which must - # conserve charge - charges = {ion.name: ion.charge for ion in ions} - for ion in (reduced, oxidized): - if ion not in reservoir: - raise ValueError( - f"The ion {ion!r} of the couple needs a reservoir concentration." - ) - if charges[oxidized] - charges[reduced] != electrons: - raise ValueError( - f"{oxidized} + {electrons} e- <-> {reduced} does not conserve " - f"charge: the charges are {charges[oxidized]} and " - f"{charges[reduced]}." - ) + if proton_flux and membrane is None: + raise ValueError("A proton flux through the membrane needs its markers.") + _check_reservoir(ions, reservoir, couple) + self._domain = domain + self._name = name + self._couple = couple self._electrode = electrode - self._reduced, self._oxidized = reduced, oxidized + self._collector = None if collector is None else tuple(collector) + self._membrane = None if membrane is None else tuple(membrane) # The regions filled with electrolyte, with the suffix of their models: # the pores of the electrode, whose porosity weights storage and # corrects transport, and free channels, of porosity one - regions = [("", self._electrode.region)] + regions = [("", electrode.region)] if channels: regions.append(("_channels", Region(tuple(channels)))) self._regions = tuple(regions) - # The unknowns. A concentration is of the order of a mole per litre, - # its scale for the dimensional analysis and the scaled equations + self._create_fields(ions, reservoir, degree, phi_s, phi_l, fields) + self._electrolyte = self._create_electrolyte( + ions, electrolyte_name, density, viscosity + ) + self._create_parameters(reservoir) + self._reaction = self._create_reaction(velocity, mass_transfer, fiber_distance) + self._model = ResidualSum( + self._create_models(velocity, outlet, proton_flux), name=name + ) + self._functionals = self._create_functionals() + + def _create_fields( + self, + ions: typing.Sequence[Ion], + reservoir: typing.Mapping[str, float | sympy.Expr], + degree: int, + phi_s: Field | None, + phi_l: Field | None, + fields: typing.Mapping[str, Field] | None, + ) -> None: + """The unknowns: the potentials, and the concentration of every ion + with a reservoir concentration. Existing fields, e.g. shared with a + membrane, are used as given. A concentration is of the order of a mole + per litre, its scale for the dimensional analysis and the scaled + equations.""" + domain = self._domain element = ("Lagrange", degree) - self._phi_s = Field("phi_s", domain, element, unit=syu.volt) - self._phi_l = Field("phi_l", domain, element, unit=syu.volt) + existing = dict(fields or {}) + self._phi_s = phi_s or Field("phi_s", domain, element, unit=syu.volt) + self._phi_l = phi_l or Field("phi_l", domain, element, unit=syu.volt) self._concentrations = types.MappingProxyType( { - ion.name: Field( + ion.name: existing.get(ion.name) + or Field( f"c_{ion.name}", domain, element, unit=CONCENTRATION, scale=1e3 ) for ion in ions if ion.name in reservoir } ) - # The electrolyte, with the field of every transported ion. The ion - # without a reservoir concentration, and so without a field, follows - # from electroneutrality. Its density and viscosity enter the mass - # transfer, and may be shared with a model of the flow - self._electrolyte = Electrolyte( - "electrolyte", - domain, + + def _create_electrolyte( + self, + ions: typing.Sequence[Ion], + name: str, + density: ufl.core.expr.Expr | float, + viscosity: ufl.core.expr.Expr | float, + ) -> Electrolyte: + """The electrolyte, with the field of every transported ion. The ion + without a reservoir concentration, and so without a field, follows + from electroneutrality. Its density and viscosity enter the mass + transfer, and may be shared with a model of the flow.""" + return Electrolyte( + name, + self._domain, [ dataclasses.replace( ion, concentration=self._concentrations.get(ion.name) @@ -292,8 +321,16 @@ def __init__( density=density, viscosity=viscosity, ) - # The concentrations of the reservoir: the initial state, and the - # values of the conditions at an inlet, see inlet_conditions + + def _create_parameters( + self, reservoir: typing.Mapping[str, float | sympy.Expr] + ) -> None: + """The parameters of the cell: the reservoir concentrations, the + initial state and the values of the conditions at an inlet; the + current density through the membrane, and the ratio of the areas + that makes the current through the collector equal, see + :meth:`balance_currents`; and the formal potential of the couple.""" + domain, name, couple = self._domain, self._name, self._couple self._inlet = types.MappingProxyType( { ion: Parameter( @@ -306,9 +343,6 @@ def __init__( for ion, value in reservoir.items() } ) - # The current density of the cell through the membrane, and the - # ratio of the areas that makes the current through the collector - # equal, see balance_currents self._current = Parameter( f"{name}.current", domain, @@ -322,90 +356,134 @@ def __init__( default=1.0, doc="Membrane area over current collector area", ) - # The kinetics of the couple, parameters of the cell that the reaction - # shares, and the specific area of the electrode self._formal_potential = Parameter( f"{name}.formal_potential", domain, - default=formal_potential, + default=couple.formal_potential, unit=syu.volt, - doc=f"Formal potential of {reduced}/{oxidized}", + doc=f"Formal potential of {couple.reduced}/{couple.oxidized}", ) - kinetics: dict[str, typing.Any] = dict( + + def _create_reaction( + self, + velocity: ufl.core.expr.Expr | None, + mass_transfer: bool, + fiber_distance: float | sympy.Expr | None, + ) -> ButlerVolmer: + """The reaction in the electrode: Butler-Volmer at the concentrations + of the pores, or at the surface of the fibres with mass transfer, by + the flow's correlation or by diffusion between the fibres. Its + parameters are those of the cell.""" + domain, name, couple = self._domain, self._name, self._couple + pair = f"{couple.reduced}/{couple.oxidized}" + common: dict[str, typing.Any] = dict( specific_area=self._electrode.specific_area, formal_potential=self._formal_potential, rate_constant=Parameter( f"{name}.rate_constant", domain, - default=rate_constant, + default=couple.rate_constant, unit=syu.m / syu.second, - doc=f"Standard rate constant of {reduced}/{oxidized}", + doc=f"Standard rate constant of {pair}", ), anodic_transfer=Parameter( f"{name}.anodic_transfer", domain, - default=anodic_transfer, - doc=f"Anodic transfer coefficient of {reduced}/{oxidized}", + default=couple.anodic_transfer, + doc=f"Anodic transfer coefficient of {pair}", ), cathodic_transfer=Parameter( f"{name}.cathodic_transfer", domain, - default=cathodic_transfer, - doc=f"Cathodic transfer coefficient of {reduced}/{oxidized}", + default=couple.cathodic_transfer, + doc=f"Cathodic transfer coefficient of {pair}", ), - electrons=electrons, + electrons=couple.electrons, markers=self._electrode.markers, + reactants=couple.reactants, + products=couple.products, ) - # The reaction in the electrode: Butler-Volmer at the concentrations - # of the pores, or at the surface of the fibres with mass transfer - e = self._electrolyte + arguments = ( + f"{name}_reaction", + self._phi_s, + self._phi_l, + self._electrolyte, + couple.reduced, + couple.oxidized, + ) + if fiber_distance is not None: + distance = Parameter( + f"{name}.fiber_distance", + domain, + default=fiber_distance, + unit=syu.m, + doc="Mean distance between the fibres", + ) + return PistonVelocityButlerVolmer( + *arguments, fiber_distance=distance, **common + ) if mass_transfer: assert velocity is not None - self._reaction = MassTransferButlerVolmer( - "reaction", - self._phi_s, - self._phi_l, - e, - reduced, - oxidized, + return MassTransferButlerVolmer( + *arguments, velocity, fiber_diameter=self._electrode.fiber_diameter, - **kinetics, - ) - else: - self._reaction = ButlerVolmer( - "reaction", self._phi_s, self._phi_l, e, reduced, oxidized, **kinetics + **common, ) + return ButlerVolmer(*arguments, **common) + def _create_models( + self, + velocity: ufl.core.expr.Expr | None, + outlet: tuple[str, ...], + proton_flux: bool, + ) -> list[Residual]: + """The models of the cell, see the class documentation.""" + name, e = self._name, self._electrolyte # The current into the solid at the collector, and out of the # electrolyte at the membrane, on discharge. Oxidation on discharge # releases electrons into the solid, which leave through the # collector, so the current of the cell enters there; reduction on # discharge takes them, so the current leaves the solid there - inflow = self._current if oxidation_on_discharge else -self._current + current = self._current + inflow = current if self._couple.oxidation_on_discharge else -current models: list[Residual] = [ Conduction( - "solid", + f"{name}_solid", self._phi_s, self._electrode.effective_conductivity(), markers=self._electrode.markers, ), self._reaction, - BoundaryFlux( - "collector_current", self._phi_s, inflow * self._area_ratio, collector - ), - BoundaryFlux("membrane_current", self._phi_l, -inflow, membrane), ] + if self._collector is not None: + collector_current = inflow * self._area_ratio + models.append( + BoundaryFlux( + f"{name}_collector_current", + self._phi_s, + collector_current, + self._collector, + ) + ) + if self._membrane is not None: + models.append( + BoundaryFlux( + f"{name}_membrane_current", self._phi_l, -inflow, self._membrane + ) + ) # Every region filled with electrolyte conserves charge, and stores # and transports every ion, with the porosity of the region for suffix, region in self._regions: models.append( - ElectrolyteConduction(f"electrolyte{suffix}", self._phi_l, e, region) + ElectrolyteConduction( + f"{name}_electrolyte{suffix}", self._phi_l, e, region + ) ) for ion, c in self._concentrations.items(): models.append( Accumulation( - f"storage_{ion}{suffix}", + f"{name}_storage_{ion}{suffix}", c, weight=region.porosity, markers=region.markers, @@ -413,7 +491,7 @@ def __init__( ) models.append( NernstPlanck( - f"transport_{ion}{suffix}", + f"{name}_transport_{ion}{suffix}", e[ion], e, self._phi_l, @@ -425,41 +503,51 @@ def __init__( if outlet: assert velocity is not None for ion, c in self._concentrations.items(): - models.append(Outflow(f"outflow_{ion}", c, velocity, outlet)) + models.append(Outflow(f"{name}_outflow_{ion}", c, velocity, outlet)) # The membrane passes the current as protons: their flux is the # current over their charge and Faraday's constant if proton_flux: + assert self._membrane is not None models.append( BoundaryFlux( - "proton_flux", + f"{name}_proton_flux", self._concentrations["H"], -inflow / (e["H"].charge * e.F), - membrane, + self._membrane, ) ) - self._model = ResidualSum(models, name=name) - - # The integrals of the diagnostics: the areas and mean potentials of - # the collector and the membrane, for the voltage, and the amounts of - # the couple in the electrolyte of every region, for the state of - # charge - self._collector = tuple(collector) - self._membrane = tuple(membrane) - ds_collector = domain.ds(*self._collector) - ds_membrane = domain.ds(*self._membrane) - functionals = { - "collector_area": 1.0 * ds_collector, - "membrane_area": 1.0 * ds_membrane, - "collector_potential": self._phi_s * ds_collector, - "membrane_potential": self._phi_l * ds_membrane, - } - for role, ion in (("reduced", reduced), ("oxidized", oxidized)): + return models + + def _create_functionals(self) -> types.MappingProxyType[str, ufl.Form]: + """The integrals of the diagnostics: the areas and mean potentials of + the collector and the membrane, for the voltage, and the amounts of + the couple in the electrolyte of every region, for the state of + charge.""" + domain = self._domain + functionals: dict[str, ufl.Form] = {} + if self._collector is not None: + ds_collector = domain.ds(*self._collector) + functionals["collector_area"] = 1.0 * ds_collector + functionals["collector_potential"] = self._phi_s * ds_collector + if self._membrane is not None: + ds_membrane = domain.ds(*self._membrane) + functionals["membrane_area"] = 1.0 * ds_membrane + functionals["membrane_potential"] = self._phi_l * ds_membrane + for role, ion in ( + ("reduced", self._couple.reduced), + ("oxidized", self._couple.oxidized), + ): c = self._concentrations[ion] functionals[f"amount_{role}"] = sum( region.storage * c * domain.dx(*region.markers) for _, region in self._regions ) - self._functionals = types.MappingProxyType(functionals) + return types.MappingProxyType(functionals) + + @property + def couple(self) -> RedoxCouple: + """The redox couple and its reaction.""" + return self._couple @property def model(self) -> Residual: @@ -532,8 +620,20 @@ def initialize(self, discretization: Discretization) -> None: values[self._phi_l] = -discretization.value(self._formal_potential) discretization.set(values) + def _check_boundaries(self, what: str) -> None: + if self._collector is None or self._membrane is None: + raise ValueError( + f"The {what} needs the facet markers of the current collector and " + "the membrane; this cell resolves them." + ) + def balance_currents(self, discretization: Discretization) -> None: - """Set the area ratio from the areas of membrane and current collector.""" + """Set the area ratio from the areas of membrane and current collector. + + Raises: + ValueError: For a cell without collector or membrane markers. + """ + self._check_boundaries("area ratio") assemble = discretization.assemble_scalar ratio = assemble(self._functionals["membrane_area"]) / assemble( self._functionals["collector_area"] @@ -545,7 +645,11 @@ def voltage(self, discretization: Discretization) -> float: That is the mean solid potential on the current collector minus the mean electrolyte potential on the membrane. + + Raises: + ValueError: For a cell without collector or membrane markers. """ + self._check_boundaries("voltage") assemble = discretization.assemble_scalar f = self._functionals return assemble(f["collector_potential"]) / assemble( @@ -561,6 +665,35 @@ def state_of_charge(self, discretization: Discretization) -> float: def local_state_of_charge(self) -> ufl.core.expr.Expr: """The share of the reduced ion in the couple, ``c_red / (c_red + c_ox)``.""" - c_red = self._concentrations[self._reduced] - c_ox = self._concentrations[self._oxidized] + c_red = self._concentrations[self._couple.reduced] + c_ox = self._concentrations[self._couple.oxidized] return c_red / (c_red + c_ox + self._electrolyte.floor) + + +def _check_reservoir( + ions: typing.Sequence[Ion], + reservoir: typing.Mapping[str, float | sympy.Expr], + couple: RedoxCouple, +) -> None: + """Check the reaction and the reservoir against the ions. + + The reservoir names the transported ions: all but one, which follows + from electroneutrality and closes the equations. Every ion of the + reaction is transported. + """ + couple.check(ions) + unknown = sorted(set(reservoir) - {ion.name for ion in ions}) + if unknown: + raise ValueError(f"The electrolyte has no ions {unknown}.") + balanced = [ion.name for ion in ions if ion.name not in reservoir] + if len(balanced) != 1: + raise ValueError( + "The reservoir must give the concentrations of all ions but one, " + "whose concentration follows from electroneutrality; it lacks " + f"{balanced}." + ) + for ion in couple.species: + if ion not in reservoir: + raise ValueError( + f"The ion {ion!r} of the reaction needs a reservoir concentration." + ) diff --git a/src/sumfem/electrochemistry/kinetics.py b/src/sumfem/electrochemistry/kinetics.py index 268fc3c..11ecf22 100644 --- a/src/sumfem/electrochemistry/kinetics.py +++ b/src/sumfem/electrochemistry/kinetics.py @@ -13,8 +13,8 @@ import sympy.physics.units as syu import ufl -from sumfem.models import butler_volmer -from sumfem.symbolic import Field, Residual, as_expr +from sumfem.models import RateReaction, butler_volmer +from sumfem.symbolic import Field, as_expr from .electrolyte import Electrolyte, Ion @@ -30,33 +30,37 @@ """A model input: an expression, a value, or ``None`` for a parameter.""" -class ElectrodeReaction(Residual): - """The reaction ``reduced -> oxidized + n e-`` in a porous electrode. +class ElectrodeReaction(RateReaction): + """An electrode reaction ``reactants -> products + n e-`` in a porous electrode. - 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)``, 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. + A :class:`~sumfem.models.RateReaction` that also moves charge between + the phases. 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. + The reaction happens ``f / (n F)`` times per unit volume and time, so it + consumes each reactant and produces each product at that rate times its + coefficient, e.g. ``V2+ -> V3+ + e-``, or ``VO2+ + H2O -> VO2+ + 2 H+ + + e-``, which also produces two protons and consumes one water molecule. + Participants without a field get no source. Subclasses give the + kinetics. Args: name: Name of the model. solid_potential: The potential of the solid, ``phi_s``. electrolyte_potential: The potential of the electrolyte, ``phi_l``. electrolyte: The electrolyte, for the constants and temperature. - reduced: Concentration of the reduced species, if any. - oxidized: Concentration of the oxidized species, if any. + reactants: The coefficient of each species that oxidation consumes, + by concentration. None for a reaction that only moves charge. + products: The coefficient of each species that oxidation produces. specific_area: Specific area ``a``. Defaults to the parameter ``".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. """ + _derives_rate = True + _requires_species = False + _solid_potential: Field """The potential of the solid.""" @@ -66,44 +70,34 @@ class ElectrodeReaction(Residual): _electrolyte: Electrolyte """The electrolyte.""" - _reduced_concentration: ufl.core.expr.Expr | None - """Concentration of the reduced species, if any.""" - - _oxidized_concentration: ufl.core.expr.Expr | None - """Concentration of the oxidized species, if any.""" - _specific_area: ufl.core.expr.Expr """Specific area of the electrode.""" _electrons: int """Number of electrons transferred.""" - _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, solid_potential: Field, electrolyte_potential: Field, electrolyte: Electrolyte, - reduced: ufl.core.expr.Expr | None = None, - oxidized: ufl.core.expr.Expr | None = None, + reactants: typing.Mapping[ufl.core.expr.Expr, int] | None = None, + products: typing.Mapping[ufl.core.expr.Expr, int] | None = None, specific_area: INPUT = None, electrons: int = 1, markers: tuple[str, ...] = (), - others: typing.Sequence[tuple[ufl.core.expr.Expr, float]] = (), ): - super().__init__(name) + # The reaction happens in the electrode, where the solid is + super().__init__( + name, + dict(reactants or {}), + dict(products or {}), + measure=solid_potential.domain.dx(*markers), + ) self._solid_potential = solid_potential self._electrolyte_potential = electrolyte_potential self._electrolyte = electrolyte - self._reduced_concentration = reduced - self._oxidized_concentration = oxidized self._specific_area = self.input( "specific_area", specific_area, @@ -112,8 +106,6 @@ def __init__( doc="Specific area of the electrode", ) self._electrons = electrons - self._markers = tuple(markers) - self._others = tuple((c, float(nu)) for c, nu in others) @property def electrolyte(self) -> Electrolyte: @@ -135,21 +127,19 @@ def source(self) -> ufl.core.expr.Expr: """The volumetric current ``f = a i``.""" return self._specific_area * self.current_density() + def rate(self) -> ufl.core.expr.Expr: + """The rate ``f / (n F)`` of the reaction, per unit volume.""" + return self.source() / (self._electrons * self._electrolyte.F) + def residual(self) -> ufl.Form: f = self.source() - 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), + dx = self._measure + current = ( + f * self._solid_potential.v * dx - f * self._electrolyte_potential.v * dx ) - for c, consumed in species: - if isinstance(c, Field): - F += consumed * rate * c.v * dx - return F + if not self._fields(): + return current + return current + super().residual() class ExchangeCurrentButlerVolmer(ElectrodeReaction): @@ -214,8 +204,8 @@ def __init__( solid_potential, electrolyte_potential, electrolyte, - reduced=reduced, - oxidized=oxidized, + reactants=None if reduced is None else {reduced: 1}, + products=None if oxidized is None else {oxidized: 1}, specific_area=specific_area, electrons=electrons, markers=markers, @@ -306,13 +296,17 @@ 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`. + reactants: The coefficient of each ion that oxidation consumes, by + name, see :class:`ElectrodeReaction`. Defaults to the reduced + ion alone, e.g. ``{"VO": 1, "H2O": 1}`` for VO2+ + H2O -> + VO2+ + 2 H+ + e-. + products: The coefficient of each ion that oxidation produces. + Defaults to the oxidized ion alone, e.g. ``{"VO2": 1, "H": 2}``. Raises: - ValueError: If an ion of the stoichiometry has no field, as the ion - that follows from electroneutrality. + ValueError: If the reactants lack the reduced ion or the products the + oxidized one, or if an ion of the reaction other than the couple + has no field, as the ion that follows from electroneutrality. """ _reduced: Ion @@ -348,26 +342,40 @@ def __init__( cathodic_transfer: INPUT = None, electrons: int = 1, markers: tuple[str, ...] = (), - stoichiometry: typing.Mapping[str, float] | None = None, + reactants: typing.Mapping[str, int] | None = None, + products: typing.Mapping[str, int] | None = None, ): - stoichiometry = dict(stoichiometry or {}) - for ion in stoichiometry: + reactants = dict(reactants or {reduced: 1}) + products = dict(products or {oxidized: 1}) + if reduced not in reactants or oxidized not in products: + raise ValueError( + f"The reactants {reactants} must hold the reduced ion {reduced!r}, " + f"and the products {products} the oxidized ion {oxidized!r}." + ) + others = (set(reactants) | set(products)) - {reduced, oxidized} + for ion in sorted(others): 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." ) + + def concentrations(ions: dict[str, int]) -> dict[ufl.core.expr.Expr, int]: + return { + electrolyte.concentration(ion, positive=False): nu + for ion, nu in ions.items() + } + super().__init__( name, solid_potential, electrolyte_potential, electrolyte, - reduced=electrolyte.concentration(reduced, positive=False), - oxidized=electrolyte.concentration(oxidized, positive=False), + reactants=concentrations(reactants), + products=concentrations(products), 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] @@ -554,7 +562,8 @@ def __init__( cathodic_transfer: INPUT = None, electrons: int = 1, markers: tuple[str, ...] = (), - stoichiometry: typing.Mapping[str, float] | None = None, + reactants: typing.Mapping[str, int] | None = None, + products: typing.Mapping[str, int] | None = None, ): super().__init__( name, @@ -570,7 +579,8 @@ def __init__( cathodic_transfer=cathodic_transfer, electrons=electrons, markers=markers, - stoichiometry=stoichiometry, + reactants=reactants, + products=products, ) domain = solid_potential.domain self._velocity = as_expr(velocity) diff --git a/src/sumfem/models/reactions.py b/src/sumfem/models/reactions.py index 59df784..d618db9 100644 --- a/src/sumfem/models/reactions.py +++ b/src/sumfem/models/reactions.py @@ -184,7 +184,9 @@ class RateReaction(Residual): measure: The measure to integrate over, e.g. ``domain.dx("Electrode")``, ``domain.ds("Surface")`` for a reaction of volume species at a surface, or an interface measure. Defaults to the whole domain - of the participants. + of the participants. A measure of cells or facets gives the + domain of the reaction, so its species may live on several + subdomains that overlap there. Example: An enzyme converts a substrate ``S`` into a product ``P`` at the @@ -225,7 +227,7 @@ def __init__( self._reactants = {as_expr(k): int(v) for k, v in reactants.items()} self._products = {as_expr(k): int(v) for k, v in products.items()} fields = self._fields() - if not fields: + if not fields and type(self)._requires_species: raise ValueError(f"Reaction {name!r} has no species that is a field.") for nu in (*self._reactants.values(), *self._products.values()): if nu < 1: @@ -234,6 +236,8 @@ def __init__( sides = self._sides() if sides is not None: self._domain = sides[0].parent + elif measure is not None and isinstance(measure.ufl_domain(), Domain): + self._domain = measure.ufl_domain() else: self._domain = self._reaction_domain(fields) if measure is None: @@ -247,6 +251,10 @@ def __init__( _derives_rate: typing.ClassVar[bool] = False """Whether a subclass derives the rate, instead of taking it.""" + _requires_species: typing.ClassVar[bool] = True + """Whether a reaction needs a species that is a field. A subclass with + other terms, e.g. the current of an electrode reaction, may have none.""" + @property def domain(self) -> Domain: """The domain of the participants, which holds the rate's parameters.""" diff --git a/tests/electrochemistry/half_cell/conftest.py b/tests/electrochemistry/half_cell/conftest.py index 76f9e95..d37e0e7 100644 --- a/tests/electrochemistry/half_cell/conftest.py +++ b/tests/electrochemistry/half_cell/conftest.py @@ -23,7 +23,7 @@ import ufl from sumfem.discrete import MeshData -from sumfem.electrochemistry import Ion, PorousElectrode, RedoxHalfCell +from sumfem.electrochemistry import Ion, PorousElectrode, RedoxCouple, RedoxHalfCell from sumfem.symbolic import Domain W = 2.5e-3 @@ -398,18 +398,24 @@ def electrode(domain: Domain, markers: tuple[str, ...] = ()) -> PorousElectrode: ) -def half_cell(domain: Domain, porous=None, **kwargs) -> RedoxHalfCell: - """The negative half-cell of the examples, V2+ oxidized on discharge.""" +COUPLE = RedoxCouple( + reduced="V2", + oxidized="V3", + rate_constant=2.5e-7 * syu.m / syu.second, + formal_potential=-0.36 * syu.volt, + oxidation_on_discharge=True, +) +"""V2+ <-> V3+ + e- of the examples, oxidized on discharge.""" + + +def half_cell(domain: Domain, porous=None, couple=None, **kwargs) -> RedoxHalfCell: + """The negative half-cell of the examples, by default with their couple.""" return RedoxHalfCell( domain, IONS, RESERVOIR, - "V2", - "V3", + COUPLE if couple is None else couple, electrode(domain) if porous is None else porous, - oxidation_on_discharge=True, - rate_constant=2.5e-7 * syu.m / syu.second, - formal_potential=-0.36 * syu.volt, **kwargs, ) @@ -427,6 +433,7 @@ def reference(): measures=measures, IONS=IONS, RESERVOIR=RESERVOIR, + COUPLE=COUPLE, electrode=electrode, half_cell=half_cell, electrochemistry=electrochemistry, diff --git a/tests/electrochemistry/half_cell/test_models.py b/tests/electrochemistry/half_cell/test_models.py index 045f0f9..9c23b2e 100644 --- a/tests/electrochemistry/half_cell/test_models.py +++ b/tests/electrochemistry/half_cell/test_models.py @@ -1,5 +1,7 @@ """The RFBniCS models: units, structure, schedules and post-processing.""" +import dataclasses + from mpi4py import MPI import dolfinx.fem @@ -9,7 +11,13 @@ import sympy.physics.units as syu from sumfem.discrete import DomainMaps, discretize -from sumfem.electrochemistry import RedoxHalfCell +from sumfem.electrochemistry import ( + Ion, + PistonVelocityButlerVolmer, + PorousElectrode, + RedoxCouple, + RedoxHalfCell, +) from sumfem.models import ( BoundaryFlux, BrinkmanForchheimer, @@ -108,8 +116,9 @@ def test_half_cells_have_consistent_units(reference): ) assert_consistent(cell.model) names = [m.name for m in cell.model.components] - assert "transport_V2_channels" in names and "storage_H_channels" in names - assert "electrolyte_channels" in names + assert "cell_transport_V2_channels" in names + assert "cell_storage_H_channels" in names + assert "cell_electrolyte_channels" in names def test_shared_parameters_are_used_by_flow_and_kinetics(reference): @@ -148,17 +157,8 @@ def test_the_couple_is_oxidized_or_reduced_on_discharge(reference): residuals = [] for oxidation in (True, False): domain = Domain.from_mesh("cell", marked.mesh) - cell = RedoxHalfCell( - domain, - reference.IONS, - reference.RESERVOIR, - "V2", - "V3", - reference.electrode(domain), - oxidation_on_discharge=oxidation, - rate_constant=0.0, - formal_potential=0.0, - ) + couple = RedoxCouple("V2", "V3", 0.0, 0.0, oxidation_on_discharge=oxidation) + cell = reference.half_cell(domain, couple=couple) scheme = BDF(1) system = build_system( cell.model, substitutions=scheme.substitutions(cell.model) @@ -198,12 +198,8 @@ def test_the_reservoir_names_ions_of_the_electrolyte(reference): domain, reference.IONS, {**reference.RESERVOIR, "Cl": 1.0}, - "V2", - "V3", + reference.COUPLE, reference.electrode(domain), - oxidation_on_discharge=True, - rate_constant=2.5e-7, - formal_potential=-0.36, ) @@ -221,21 +217,23 @@ def test_post_processing(reference): def test_the_couple_conserves_charge(reference): - """V3 + 2 e- <-> V2 does not conserve charge, nor V3 + e- <-> SO4.""" + """V2 <-> V3 + 2 e- does not conserve charge. SO4 <-> HSO4 + e- does, but + sulfate follows from electroneutrality, so it has no field to react.""" domain = Domain("cell", "interval") with pytest.raises(ValueError, match="does not conserve charge"): - reference.half_cell(domain, electrons=2) + reference.half_cell( + domain, couple=dataclasses.replace(reference.COUPLE, electrons=2) + ) with pytest.raises(ValueError, match="needs a reservoir"): - RedoxHalfCell( + reference.half_cell( domain, - reference.IONS, - reference.RESERVOIR, - "SO4", - "V3", - reference.electrode(domain), - oxidation_on_discharge=True, - rate_constant=2.5e-7, - formal_potential=-0.36, + couple=dataclasses.replace( + reference.COUPLE, + reduced="SO4", + oxidized="HSO4", + reactants=None, + products=None, + ), ) @@ -250,10 +248,156 @@ def test_the_reservoir_lacks_exactly_one_ion(reference, missing): domain, reference.IONS, reservoir, - "V2", - "V3", + reference.COUPLE, reference.electrode(domain), - oxidation_on_discharge=True, - rate_constant=2.5e-7, - formal_potential=-0.36, ) + + +CONCENTRATION_UNIT = syu.mol / syu.m**3 +DIFFUSIVITY_UNIT = syu.m**2 / syu.second + +POSITIVE_IONS = ( + Ion("VO", 2, 3.9e-10 * DIFFUSIVITY_UNIT, positive=True), + Ion("VO2", 1, 3.9e-10 * DIFFUSIVITY_UNIT, positive=True), + Ion("H", 1, 9.312e-9 * DIFFUSIVITY_UNIT), + Ion("HSO4", -1, 1.33e-9 * DIFFUSIVITY_UNIT), + Ion("H2O", 0, 2.3e-9 * DIFFUSIVITY_UNIT), + Ion("SO4", -2, 1.065e-9 * DIFFUSIVITY_UNIT), +) +"""The positive electrolyte of Shah et al. (2008): VO2+, VO2+, protons, +bisulfate and water, sulfate by electroneutrality.""" + +POSITIVE_RESERVOIR = { + "VO": 1053.0 * CONCENTRATION_UNIT, + "VO2": 27.0 * CONCENTRATION_UNIT, + "H": 1200.0 * CONCENTRATION_UNIT, + "HSO4": 1200.0 * CONCENTRATION_UNIT, + "H2O": 4.2e4 * CONCENTRATION_UNIT, +} + + +POSITIVE_COUPLE = RedoxCouple( + reduced="VO", + oxidized="VO2", + rate_constant=0.68 * 3e-9 * syu.m / syu.second, + formal_potential=1.004 * syu.volt, + oxidation_on_discharge=False, + reactants={"VO": 1, "H2O": 1}, + products={"VO2": 1, "H": 2}, +) +"""VO2+ + H2O <-> VO2+ + 2 H+ + e- of Shah et al., reduced on discharge.""" + + +def positive_half_cell(domain, electrode, couple=POSITIVE_COUPLE, **kwargs): + """The positive half-cell of Shah et al.""" + return RedoxHalfCell( + domain, POSITIVE_IONS, POSITIVE_RESERVOIR, couple, electrode, **kwargs + ) + + +def test_a_positive_half_cell_produces_protons_and_consumes_water(reference): + """VO2+ + H2O <-> VO2+ + 2 H+ + e- conserves charge only with its protons, + and diffusion between the fibres limits it without flow.""" + domain = Domain("cell", "interval") + electrode = reference.electrode(domain) + cell = positive_half_cell(domain, electrode, fiber_distance=10 * syu.micrometer) + assert isinstance(cell.reaction, PistonVelocityButlerVolmer) + assert_consistent(cell.model) + with pytest.raises(ValueError, match="does not conserve charge"): + positive_half_cell( + domain, + electrode, + dataclasses.replace(POSITIVE_COUPLE, reactants=None, products=None), + ) + with pytest.raises(ValueError, match="has no ions"): + # An electrolyte without water + RedoxHalfCell( + domain, + tuple(ion for ion in POSITIVE_IONS if ion.name != "H2O"), + {k: v for k, v in POSITIVE_RESERVOIR.items() if k != "H2O"}, + POSITIVE_COUPLE, + electrode, + ) + + +def test_existing_fields_are_used_as_given(reference): + domain = Domain("cell", "interval") + phi_s = Field("psi", domain, ("Lagrange", 1), unit=syu.volt) + phi_l = Field("phi", domain, ("Lagrange", 1), unit=syu.volt) + protons = Field("c_protons", domain, ("Lagrange", 1), unit=CONCENTRATION_UNIT) + cell = reference.half_cell(domain, phi_s=phi_s, phi_l=phi_l, fields={"H": protons}) + assert cell.phi_s is phi_s and cell.phi_l is phi_l + assert cell.concentrations["H"] is protons + assert cell.electrolyte["H"].field is protons + assert cell.concentrations["V2"].name == "c_V2" + + +def test_a_resolved_collector_and_membrane_carry_no_currents(reference): + """In a full cell the collector and the membrane are domains of their + own, so the half-cell puts no currents on its boundary.""" + domain = Domain("cell", "interval") + cell = reference.half_cell(domain, collector=None, membrane=None) + names = [m.name for m in cell.model.components] + assert not any("current" in name for name in names) + with pytest.raises(ValueError, match="resolves them"): + cell.voltage(None) + with pytest.raises(ValueError, match="needs its markers"): + reference.half_cell(domain, membrane=None, proton_flux=True) + + +def test_two_half_cells_share_the_electrolyte_potential_and_water(reference): + """A negative and a positive half-cell side by side, with one electrolyte + potential and one water field, as in a full cell: their parameters and + models have names of their own, and their sum has consistent units.""" + domain = Domain("cell", "interval") + phi_l = Field("phi_l", domain, ("Lagrange", 1), unit=syu.volt) + water = Field("c_H2O", domain, ("Lagrange", 1), unit=CONCENTRATION_UNIT) + sides = {} + for side in ("negative", "positive"): + fields = { + ion: Field( + f"{side}_c_{ion}", domain, ("Lagrange", 1), unit=CONCENTRATION_UNIT + ) + for ion in ("H", "HSO4") + } + common = dict( + phi_s=Field(f"{side}_phi_s", domain, ("Lagrange", 1), unit=syu.volt), + phi_l=phi_l, + collector=None, + membrane=None, + fiber_distance=10 * syu.micrometer, + name=side, + electrolyte_name=f"{side}_electrolyte", + ) + electrode = PorousElectrode( + f"{side}_electrode", + domain, + (side.capitalize(),), + porosity=0.68, + conductivity=500.0, + specific_area=2e6, + ) + if side == "negative": + ions = (*reference.IONS[:4], Ion("H2O", 0, 2.3e-9), reference.IONS[4]) + reservoir = {**reference.RESERVOIR, "H2O": 4.2e4} + couple = RedoxCouple("V2", "V3", 0.68 * 1.75e-7, -0.255, True) + sides[side] = RedoxHalfCell( + domain, + ions, + reservoir, + couple, + electrode, + fields={**fields, "H2O": water}, + **common, + ) + else: + sides[side] = positive_half_cell( + domain, electrode, fields={**fields, "H2O": water}, **common + ) + model = sides["negative"].model + sides["positive"].model + system = build_system(model, check_units=True) + names = [f.name for f in system.unknowns] + assert names.count("phi_l") == 1 and names.count("c_H2O") == 1 + assert {"negative_c_H", "positive_c_H", "c_VO", "c_V2"} <= set(names) + assert "negative.rate_constant" in system.parameters + assert "positive_electrolyte.D_H" in system.parameters diff --git a/tests/electrochemistry/test_electrochemistry_from_ufl.py b/tests/electrochemistry/test_electrochemistry_from_ufl.py index 22effa5..f32181a 100644 --- a/tests/electrochemistry/test_electrochemistry_from_ufl.py +++ b/tests/electrochemistry/test_electrochemistry_from_ufl.py @@ -303,7 +303,8 @@ def test_positive_couple_with_protons_water_and_piston_velocities(cell, box, com specific_area=1e3, rate_constant=1e-6, formal_potential=1.0, - stoichiometry={"H": 2, "H2O": -1}, + reactants={"VO": 1, "H2O": 1}, + products={"VO2": 1, "H": 2}, ) class Mass(Residual): @@ -334,7 +335,7 @@ def residual(self): ) -def test_the_stoichiometry_needs_ions_with_fields(): +def test_the_reaction_needs_its_couple_and_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")] @@ -346,8 +347,15 @@ def test_the_stoichiometry_needs_ions_with_fields(): electrolyte = Electrolyte("e", domain, ions) with pytest.raises(ValueError, match="electroneutrality"): ButlerVolmer( - "bv", *potentials, electrolyte, "VO", "VO2", stoichiometry={"H": 2} + "bv", + *potentials, + electrolyte, + "VO", + "VO2", + products={"VO2": 1, "H": 2}, ) + with pytest.raises(ValueError, match="must hold the reduced ion"): + ButlerVolmer("bv", *potentials, electrolyte, "VO", "VO2", reactants={"H": 1}) def test_mass_transfer_takes_the_properties_of_the_electrolyte():