diff --git a/docs/api/electrochemistry.md b/docs/api/electrochemistry.md index d9e3d9a..47fbcca 100644 --- a/docs/api/electrochemistry.md +++ b/docs/api/electrochemistry.md @@ -32,6 +32,9 @@ .. autoclass:: NernstPlanck :members: + +.. autoclass:: TransientNernstPlanck + :members: ``` ## Kinetics diff --git a/src/sumfem/electrochemistry/__init__.py b/src/sumfem/electrochemistry/__init__.py index 4850bcb..b7f8da7 100644 --- a/src/sumfem/electrochemistry/__init__.py +++ b/src/sumfem/electrochemistry/__init__.py @@ -5,7 +5,7 @@ - Porous electrodes: :class:`PorousElectrode`, and the :class:`Region` of a cell filled with electrolyte. - Charge and transport in the electrolyte: :class:`ElectrolyteConduction` - and :class:`NernstPlanck`. + and :class:`NernstPlanck`, with storage :class:`TransientNernstPlanck`. - Kinetics: :class:`ExchangeCurrentButlerVolmer`, :class:`ButlerVolmer`, :class:`MassTransferButlerVolmer` and :class:`PistonVelocityButlerVolmer`, on :class:`ElectrodeReaction`, with the Butler-Volmer equation @@ -29,7 +29,7 @@ PistonVelocityButlerVolmer, ) from .porous import PorousElectrode, Region -from .transport import NernstPlanck +from .transport import NernstPlanck, TransientNernstPlanck __all__ = [ "Electrolyte", @@ -38,6 +38,7 @@ "Region", "ElectrolyteConduction", "NernstPlanck", + "TransientNernstPlanck", "RedoxCouple", "RedoxHalfCell", "ButlerVolmer", diff --git a/src/sumfem/electrochemistry/charge.py b/src/sumfem/electrochemistry/charge.py index a32c477..a19dad6 100644 --- a/src/sumfem/electrochemistry/charge.py +++ b/src/sumfem/electrochemistry/charge.py @@ -13,13 +13,26 @@ class ElectrolyteConduction(Residual): - """Charge conservation in the electrolyte, ``-(i, grad(v))``. + r"""Charge conservation in the electrolyte, :math:`-(\mathbf{i}_l, \nabla v)`. - The current density ``i`` of the electrolyte, + The current density of the electrolyte, :meth:`~sumfem.electrochemistry.electrolyte.Electrolyte.current_density`, - includes migration with the concentration-dependent conductivity and - diffusion of every ion. Diffusivities take the transport correction of - the region. + carries charge by migration, with the concentration-dependent + conductivity :math:`\kappa_\text{eff}`, and by diffusion of every ion + :math:`k`: + + .. math:: + :label: electrolyte-current + + \begin{aligned} + \mathbf{i}_l &= -\kappa_\text{eff} \nabla \phi + - F \sum_k z_k D_{k,\text{eff}} \nabla c_k, \\ + \kappa_\text{eff} &= \frac{F^2}{R T} + \sum_k z_k^2 D_{k,\text{eff}} c_k, + \end{aligned} + + with the diffusivities corrected for the region as in + :eq:`nernst-planck-flux`. Args: name: Name of the model. diff --git a/src/sumfem/electrochemistry/half_cell.py b/src/sumfem/electrochemistry/half_cell.py index 35afe7f..504e632 100644 --- a/src/sumfem/electrochemistry/half_cell.py +++ b/src/sumfem/electrochemistry/half_cell.py @@ -18,7 +18,7 @@ import ufl from sumfem.discrete import Discretization -from sumfem.models import Accumulation, BoundaryFlux, Conduction, Outflow +from sumfem.models import BoundaryFlux, Conduction, Outflow from sumfem.symbolic import ( DirichletBC, Domain, @@ -37,105 +37,163 @@ PistonVelocityButlerVolmer, ) from .porous import PorousElectrode, Region -from .transport import NernstPlanck +from .transport import TransientNernstPlanck __all__ = ["RedoxHalfCell"] class RedoxHalfCell: - """A half-cell of a redox flow battery, of a couple of dissolved ions. + r"""A half-cell of a redox flow battery, of a couple of dissolved ions. The macro-homogeneous porous-electrode model of RFBniCS, `Aghabarari et al., 2026 `_, Section II, for any couple, electrolyte and sign of the half-cell. - 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 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 unknowns are the potential :math:`\phi_s` of the solid and + :math:`\phi_l` of the electrolyte, and the concentration :math:`c_k` of + every ion with a reservoir concentration. In the porous electrode, of + porosity :math:`\varepsilon` and specific area :math:`a`, charge is + conserved in the solid and the electrolyte, which exchange the current + :math:`a i` of the reaction, and each ion is stored, transported and + produced: + + .. math:: + :label: redox-half-cell + + \begin{aligned} + -\nabla \cdot (\sigma_\text{eff} \nabla \phi_s) + a i &= 0, \\ + \nabla \cdot \mathbf{i}_l - a i &= 0, \\ + \varepsilon \frac{\partial c_k}{\partial t} + + \nabla \cdot \mathbf{N}_k &= \nu_k \frac{a i}{n F}, \\ + \sum_k z_k c_k &= 0. + \end{aligned} + + The terms are those of the components of the cell: + + - :math:`\sigma_\text{eff}`, the effective conductivity of the solid, + :meth:`PorousElectrode.effective_conductivity + `; + - :math:`\mathbf{i}_l`, the current density of the electrolyte, + :eq:`electrolyte-current`; + - :math:`\mathbf{N}_k`, the Nernst-Planck flux of the ion of charge + :math:`z_k`, :eq:`nernst-planck-flux`, with the superficial velocity + :math:`\mathbf{u}` of the electrolyte; + - :math:`i`, the interfacial current density of the couple, + :eq:`butler-volmer-couple`, at surface concentrations limited by mass + transfer, :eq:`mass-transfer-balance`, with flow or a fibre distance; + - :math:`\nu_k` and :math:`n`, the coefficient of the ion in the + reaction of the couple, negative for a reactant, positive for a + product and zero for an ion that takes no part, and the number of + electrons, see :class:`~sumfem.electrochemistry.RedoxCouple`. + + The one ion without a reservoir concentration follows from + electroneutrality, the last equation. In free channels, + :math:`\varepsilon = 1` and there is no reaction or solid. + + The cell is galvanostatic, with the current density :math:`j` of the + cell through the membrane, positive on discharge. On discharge the + couple is oxidized in the negative half-cell, :math:`s = 1`, and reduced + in the positive one, :math:`s = -1`. With the outward normal + :math:`\mathbf{n}`, the current enters the solid at the current + collector and leaves the electrolyte at the membrane, + + .. math:: + :label: redox-half-cell-boundary + + \begin{aligned} + \sigma_{\text{eff}} \nabla \phi_s \cdot \mathbf{n} &= s j r + && \text{at the current collector}, \\ + \mathbf{i}_l \cdot \mathbf{n} &= s j + && \text{at the membrane}, \\ + \mathbf{N}_\text{H} \cdot \mathbf{n} &= \frac{s j}{z_\text{H} F} + && \text{at the membrane, with a proton flux}, \\ + (\mathbf{N}_k - \mathbf{u} c_k) \cdot \mathbf{n} &= 0 + && \text{at an outlet}, + \end{aligned} + + where the ratio :math:`r` of the membrane area to the collector area + makes the total currents equal, see :meth:`balance_currents`. Other + boundaries carry no flux. The model, :attr:`model`, is the sum of: - - ``"_solid"``: Ohmic conduction in the solid, - :class:`~sumfem.models.Conduction` with Bruggeman's effective - conductivity of the electrode; - - ``"_reaction"``: the kinetics of the couple, + - ``"_solid"``: the first equation of :eq:`redox-half-cell` + without the reaction, + :class:`~sumfem.models.Conduction`; + - ``"_reaction"``: the kinetics :eq:`butler-volmer-couple`, with + :eq:`mass-transfer-balance` under mass transfer, :class:`~sumfem.electrochemistry.ButlerVolmer`, :class:`~sumfem.electrochemistry.MassTransferButlerVolmer` with mass 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 a fibre distance, which adds the reaction to + :eq:`redox-half-cell`; - 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`; + and ``"_membrane_current"``: the first two conditions of + :eq:`redox-half-cell-boundary`, :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, + channels: the charge in the electrolyte, + :class:`~sumfem.electrochemistry.ElectrolyteConduction`, and the + balance :eq:`nernst-planck-balance` of every ion, + :class:`~sumfem.electrochemistry.TransientNernstPlanck`, named + ``"_"``, whose summands are ``"__storage"`` and + ``"__transport"``; + - with an outlet, ``"_outflow_"``: the outlet condition, :class:`~sumfem.models.Outflow`; - - with a proton flux, ``"_proton_flux"``: the membrane current - carried by protons, a :class:`~sumfem.models.BoundaryFlux` of the ion + - with a proton flux, ``"_proton_flux"``: the proton condition at + the membrane, 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. - On discharge the couple is oxidized in the negative half-cell and - reduced in the positive one, so the current enters the solid at the - current collector of the first, and leaves it at that of the second. - The current collector carries the current density ``current * - area_ratio``, so that the total currents match when the two areas - differ. - Args: domain: The domain of the cell. - ions: The ions of the electrolyte, without fields: the cell creates + ions: The ions of the electrolyte, with their charges :math:`z_k` + and diffusivities :math:`D_k`, without fields: the cell creates the field of every ion that has a reservoir concentration. 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_"``. - 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. + couple: The redox couple and its reaction: the coefficients + :math:`\nu_k`, :math:`n` and :math:`s`, and the kinetics + :eq:`butler-volmer-couple`, whose parameters are those of the + cell: ``".rate_constant"`` :math:`k_0`, + ``".formal_potential"`` :math:`E_0`, + ``".anodic_transfer"`` :math:`\alpha_a` and + ``".cathodic_transfer"`` :math:`\alpha_c`. + electrode: The porous electrode: :math:`\varepsilon`, :math:`a` and + :math:`\sigma`. channels: Cell markers of free channels. Empty for none. - velocity: Superficial velocity of the electrolyte, or ``None`` + velocity: The superficial velocity :math:`\mathbf{u}`, 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. - 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. + :eq:`mass-transfer-balance`, with the coefficients of the flow's + correlation, :eq:`mass-transfer-correlation`, which needs a + velocity. Defaults to whether there is one. + fiber_distance: The mean distance :math:`\delta` between the + fibres, in m or with units. Limits the kinetics by + :eq:`mass-transfer-balance` with the piston velocities, + :eq:`piston-velocity`, 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. + outlet: Facet markers of an outlet. Empty for none. proton_flux: Carry the membrane current by protons, the ion ``"H"``, - through the membrane. - density: Density of the electrolyte, see + through the membrane, the third condition of + :eq:`redox-half-cell-boundary`. + density: Density of the electrolyte, for the mass transfer, see :class:`~sumfem.electrochemistry.Electrolyte`. Pass the parameter of a flow model to share it. 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, ``"."``, - 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. + and of the names of its models, ``"_"``. The + current density :math:`j` is the parameter ``".current"``, + in A/m^2, and the area ratio :math:`r` ``".area_ratio"``. + phi_s: The potential :math:`\phi_s` of the solid, if it exists, e.g. + over the electrode and its current collector. Created on the + domain otherwise. + phi_l: The potential :math:`\phi_l` 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 @@ -190,6 +248,19 @@ class RedoxHalfCell: _formal_potential: Parameter """Formal potential of the couple.""" + _rate_constant: Parameter + """Standard rate constant of the couple.""" + + _anodic_transfer: Parameter + """Anodic transfer coefficient of the couple.""" + + _cathodic_transfer: Parameter + """Cathodic transfer coefficient of the couple.""" + + _fiber_distance: Parameter | None + """Mean distance between the fibres, if the piston velocities limit the + kinetics.""" + _reaction: ButlerVolmer """The reaction in the electrode.""" @@ -261,8 +332,8 @@ def __init__( 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._create_parameters(reservoir, fiber_distance) + self._reaction = self._create_reaction(velocity, mass_transfer) self._model = ResidualSum( self._create_models(velocity, outlet, proton_flux), name=name ) @@ -277,11 +348,22 @@ def _create_fields( 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.""" + """Create the unknowns of the cell. + + Concentrations get the scale of a mole per litre, for the + dimensional analysis and the scaled equations. + + Args: + ions: The ions of the electrolyte. + reservoir: The reservoir concentrations. Their ions are the + transported ones, each with a field; the one ion without + follows from electroneutrality. + degree: Degree of the Lagrange elements of new fields. + phi_s: The potential of the solid, or ``None`` to create it. + phi_l: The potential of the electrolyte, or ``None`` to create it. + fields: Existing concentration fields, by ion name, used as + given. + """ domain = self._domain element = ("Lagrange", degree) existing = dict(fields or {}) @@ -305,10 +387,20 @@ def _create_electrolyte( 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.""" + """Create the electrolyte, with the fields of the transported ions. + + Args: + ions: The ions of the electrolyte. The one without a field + follows from electroneutrality. + name: Prefix of the electrolyte's parameters. + density: Density of the electrolyte, a value or a shared + expression. + viscosity: Dynamic viscosity of the electrolyte, a value or a + shared expression. + + Returns: + The electrolyte. + """ return Electrolyte( name, self._domain, @@ -323,14 +415,25 @@ def _create_electrolyte( ) def _create_parameters( - self, reservoir: typing.Mapping[str, float | sympy.Expr] + self, + reservoir: typing.Mapping[str, float | sympy.Expr], + fiber_distance: float | sympy.Expr | None, ) -> 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.""" + """Create the parameters of the cell. + + These are the reservoir concentrations, the current density, the + area ratio, the kinetic parameters of the couple, and the fibre + distance if one is given. + + Args: + reservoir: The default reservoir concentration of each + transported ion. + fiber_distance: The default mean distance between the fibres, or + ``None`` for kinetics without piston velocities, which then + has no such parameter. + """ domain, name, couple = self._domain, self._name, self._couple + pair = f"{couple.reduced}/{couple.oxidized}" self._inlet = types.MappingProxyType( { ion: Parameter( @@ -361,43 +464,62 @@ def _create_parameters( domain, default=couple.formal_potential, unit=syu.volt, - doc=f"Formal potential of {couple.reduced}/{couple.oxidized}", + doc=f"Formal potential of {pair}", + ) + self._rate_constant = Parameter( + f"{name}.rate_constant", + domain, + default=couple.rate_constant, + unit=syu.m / syu.second, + doc=f"Standard rate constant of {pair}", + ) + self._anodic_transfer = Parameter( + f"{name}.anodic_transfer", + domain, + default=couple.anodic_transfer, + doc=f"Anodic transfer coefficient of {pair}", + ) + self._cathodic_transfer = Parameter( + f"{name}.cathodic_transfer", + domain, + default=couple.cathodic_transfer, + doc=f"Cathodic transfer coefficient of {pair}", + ) + self._fiber_distance = ( + None + if fiber_distance is None + else Parameter( + f"{name}.fiber_distance", + domain, + default=fiber_distance, + unit=syu.m, + doc="Mean distance between the fibres", + ) ) def _create_reaction( - self, - velocity: ufl.core.expr.Expr | None, - mass_transfer: bool, - fiber_distance: float | sympy.Expr | None, + self, velocity: ufl.core.expr.Expr | None, mass_transfer: bool ) -> 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}" + """Create the reaction in the electrode, from the cell's parameters. + + Args: + velocity: The superficial velocity, or ``None`` without flow. + mass_transfer: Limit the kinetics by the flow's mass-transfer + correlation, unless the cell has a fibre distance. + + Returns: + The reaction: :class:`~sumfem.electrochemistry.PistonVelocityButlerVolmer` + with a fibre distance, else + :class:`~sumfem.electrochemistry.MassTransferButlerVolmer` with + mass transfer, else :class:`~sumfem.electrochemistry.ButlerVolmer`. + """ + name, couple = self._name, self._couple 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=couple.rate_constant, - unit=syu.m / syu.second, - doc=f"Standard rate constant of {pair}", - ), - anodic_transfer=Parameter( - f"{name}.anodic_transfer", - domain, - default=couple.anodic_transfer, - doc=f"Anodic transfer coefficient of {pair}", - ), - cathodic_transfer=Parameter( - f"{name}.cathodic_transfer", - domain, - default=couple.cathodic_transfer, - doc=f"Cathodic transfer coefficient of {pair}", - ), + rate_constant=self._rate_constant, + anodic_transfer=self._anodic_transfer, + cathodic_transfer=self._cathodic_transfer, electrons=couple.electrons, markers=self._electrode.markers, reactants=couple.reactants, @@ -411,16 +533,9 @@ def _create_reaction( 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", - ) + if self._fiber_distance is not None: return PistonVelocityButlerVolmer( - *arguments, fiber_distance=distance, **common + *arguments, fiber_distance=self._fiber_distance, **common ) if mass_transfer: assert velocity is not None @@ -438,7 +553,16 @@ def _create_models( outlet: tuple[str, ...], proton_flux: bool, ) -> list[Residual]: - """The models of the cell, see the class documentation.""" + """Create the models of the cell, listed in the class documentation. + + Args: + velocity: The superficial velocity, or ``None`` without flow. + outlet: Facet markers of the outlet. Empty for none. + proton_flux: Carry the membrane current by protons. + + Returns: + The models, to be summed. + """ 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 @@ -480,18 +604,10 @@ def _create_models( f"{name}_electrolyte{suffix}", self._phi_l, e, region ) ) - for ion, c in self._concentrations.items(): - models.append( - Accumulation( - f"{name}_storage_{ion}{suffix}", - c, - weight=region.porosity, - markers=region.markers, - ) - ) + for ion in self._concentrations: models.append( - NernstPlanck( - f"{name}_transport_{ion}{suffix}", + TransientNernstPlanck( + f"{name}_{ion}{suffix}", e[ion], e, self._phi_l, @@ -519,10 +635,13 @@ def _create_models( 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.""" + """Create the integrals of the voltage and the state of charge. + + Returns: + The areas and potential integrals of the collector and the + membrane, when they have markers, and the amounts of the reduced + and oxidized ion, by name. + """ domain = self._domain functionals: dict[str, ufl.Form] = {} if self._collector is not None: @@ -621,6 +740,14 @@ def initialize(self, discretization: Discretization) -> None: discretization.set(values) def _check_boundaries(self, what: str) -> None: + """Check that the collector and the membrane have facet markers. + + Args: + what: The quantity that needs them, for the error message. + + Raises: + ValueError: If either has 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 " @@ -677,9 +804,15 @@ def _check_reservoir( ) -> 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. + Args: + ions: The ions of the electrolyte. + reservoir: The reservoir concentrations. They must name all ions + but one, which follows from electroneutrality, and every ion of + the reaction. + couple: The redox couple. + + Raises: + ValueError: If the reaction or the reservoir does not match the ions. """ couple.check(ions) unknown = sorted(set(reservoir) - {ion.name for ion in ions}) diff --git a/src/sumfem/electrochemistry/kinetics.py b/src/sumfem/electrochemistry/kinetics.py index 11ecf22..bfda662 100644 --- a/src/sumfem/electrochemistry/kinetics.py +++ b/src/sumfem/electrochemistry/kinetics.py @@ -259,20 +259,30 @@ def current_density(self) -> ufl.core.expr.Expr: class ButlerVolmer(ElectrodeReaction): - """Butler-Volmer kinetics of a redox couple of an electrolyte. + r"""Butler-Volmer kinetics of a redox couple of an electrolyte. The Butler-Volmer equation, see :func:`~sumfem.models.butler_volmer`, with an exchange current density from the concentrations of the couple. The interfacial current density is - ``i = n F (k_a c_red - k_c c_ox)``, - ``k_a = k0 exp(alpha_a n F eta / (R T))``, - ``k_c = k0 exp(-alpha_c n F eta / (R T))``, - - with ``eta = phi_s - phi_l - E0`` relative to the formal potential - ``E0``, and the concentrations at the electrode surface, by default the - positive parts of the bulk concentrations. The reduced and oxidized ions - are ions in the electrolyte, given by name. + .. math:: + :label: butler-volmer-couple + + \begin{aligned} + i &= n F \left( k_a c_\text{red}^s - k_c c_\text{ox}^s \right), \\ + k_a &= k_0 e^{\alpha_a n F \eta / (R T)}, \\ + k_c &= k_0 e^{-\alpha_c n F \eta / (R T)}, \\ + \eta &= \phi_s - \phi_l - E_0, + \end{aligned} + + with the standard rate constant :math:`k_0`, the transfer coefficients + :math:`\alpha_a` and :math:`\alpha_c`, the formal potential :math:`E_0`, + and the concentrations :math:`c_\text{red}^s` and :math:`c_\text{ox}^s` + of the reduced and oxidized ion at the electrode surface. These are, by + default, the smooth positive parts of the bulk concentrations, + :math:`c^+ = (c + \sqrt{c^2 + c_\text{floor}^2}) / 2`, see + :meth:`~sumfem.electrochemistry.electrolyte.Electrolyte.positive_part`. + The reduced and oxidized ions are ions in the electrolyte, given by name. Subclasses change the kinetics by overriding :meth:`surface_concentrations`, :meth:`rate_constants` or @@ -476,14 +486,40 @@ def overpotential(self) -> ufl.core.expr.Expr: class MassTransferButlerVolmer(ButlerVolmer): - """Butler-Volmer kinetics limited by mass transfer to the fibers. - - The surface concentrations balance the reaction with mass transfer - ``k_m (c - c_s)`` from the bulk. The mass-transfer coefficient of each - ion is ``k_m = Sh D / d`` with the Sherwood number ``Sh = a Re^b`` and - the Reynolds number ``Re = rho |u| d / mu`` of the flow past fibers of - diameter ``d``. ``D`` is the diffusivity in free electrolyte. The - defaults ``a = 7`` and ``b = 0.4`` are the correlation of `Krishnamurthy + r"""Butler-Volmer kinetics limited by mass transfer to the fibers. + + The surface concentrations of :eq:`butler-volmer-couple` are those at + which mass transfer from the bulk, with the coefficients + :math:`m_\text{red}` and :math:`m_\text{ox}`, balances the reaction: + + .. math:: + :label: mass-transfer-balance + + \begin{aligned} + m_\text{red} \left( c_\text{red}^+ - c_\text{red}^s \right) + &= k_a c_\text{red}^s - k_c c_\text{ox}^s, \\ + m_\text{ox} \left( c_\text{ox}^s - c_\text{ox}^+ \right) + &= k_a c_\text{red}^s - k_c c_\text{ox}^s, + \end{aligned} + + two linear equations, solved in closed form. Each coefficient has a + small floor added. The coefficients follow from the Sherwood number + :math:`\mathrm{Sh} = a \, \mathrm{Re}^b` of the flow past fibers of + diameter :math:`d`: + + .. math:: + :label: mass-transfer-correlation + + \begin{aligned} + m_\text{red} &= \mathrm{Sh} \frac{D_\text{red}}{d}, \\ + m_\text{ox} &= \mathrm{Sh} \frac{D_\text{ox}}{d}, \\ + \mathrm{Re} &= \frac{\rho |\mathbf{u}| d}{\mu}, + \end{aligned} + + with the diffusivities :math:`D_\text{red}` and :math:`D_\text{ox}` of + the reduced and oxidized ion in free electrolyte, and the density + :math:`\rho` and viscosity :math:`\mu` of the electrolyte. The + defaults :math:`a = 7` and :math:`b = 0.4` are the correlation of `Krishnamurthy et al., 2011 `_, and the surface concentrations those of equation 20 of `Aghabarari et al., 2026 `_. @@ -669,9 +705,20 @@ def _balanced_surface( ) -> 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 + Solves ``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 `_. + + Args: + k_a: Anodic rate constant. + k_c: Cathodic rate constant. + m_red: Mass-transfer coefficient of the reduced ion. + m_ox: Mass-transfer coefficient of the oxidized ion. + c_red: Bulk concentration of the reduced ion. + c_ox: Bulk concentration of the oxidized ion. + + Returns: + The surface concentrations ``s_red`` and ``s_ox``. """ 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 @@ -680,13 +727,22 @@ def _balanced_surface( class PistonVelocityButlerVolmer(ButlerVolmer): - """Butler-Volmer kinetics limited by diffusion between the fibres. + r"""Butler-Volmer kinetics limited by diffusion between the fibres. + + The surface concentrations balance the reaction with mass transfer as in + :eq:`mass-transfer-balance`, with the piston velocities + + .. math:: + :label: piston-velocity + + \begin{aligned} + m_\text{red} &= \frac{D_\text{red}}{\delta}, \\ + m_\text{ox} &= \frac{D_\text{ox}}{\delta}, + \end{aligned} - 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 + the diffusivities of the reduced and oxidized ion in free electrolyte + over the mean distance :math:`\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. diff --git a/src/sumfem/electrochemistry/porous.py b/src/sumfem/electrochemistry/porous.py index 3386866..ef29544 100644 --- a/src/sumfem/electrochemistry/porous.py +++ b/src/sumfem/electrochemistry/porous.py @@ -37,7 +37,10 @@ def storage(self) -> ufl.core.expr.Expr | float: @property def transport(self) -> ufl.core.expr.Expr | float: - """Correction of diffusivities, Bruggeman's: one in free electrolyte.""" + r"""Correction of diffusivities, Bruggeman's :math:`\varepsilon^b`. + + It is one in free electrolyte, with porosity one. + """ return bruggeman(1.0, self.porosity, self.bruggeman_exponent) @@ -163,7 +166,12 @@ def fiber_diameter(self) -> Parameter: return self._fiber_diameter def effective_conductivity(self) -> ufl.core.expr.Expr: - """Bruggeman's conductivity of the solid, ``(1 - eps)^1.5 sigma``.""" + r"""Bruggeman's conductivity of the solid. + + That is :math:`\sigma_\text{eff} = (1 - \varepsilon)^{3/2} \sigma`, with + the conductivity :math:`\sigma` of the fibres and the porosity + :math:`\varepsilon`. + """ return bruggeman(self._conductivity, 1.0 - self._porosity) def permeability(self) -> ufl.core.expr.Expr: diff --git a/src/sumfem/electrochemistry/transport.py b/src/sumfem/electrochemistry/transport.py index 3e234ef..694bd87 100644 --- a/src/sumfem/electrochemistry/transport.py +++ b/src/sumfem/electrochemistry/transport.py @@ -4,22 +4,34 @@ import ufl -from sumfem.symbolic import Field, Residual, as_expr +from sumfem.models import Accumulation +from sumfem.symbolic import Field, Residual, ResidualSum, as_expr from .electrolyte import Electrolyte, Ion from .porous import Region -__all__ = ["NernstPlanck"] +__all__ = ["NernstPlanck", "TransientNernstPlanck"] class NernstPlanck(Residual): - """Nernst-Planck transport of one ion, ``-(N, grad(w))``. + r"""Nernst-Planck transport of one ion, :math:`-(\mathbf{N}, \nabla w)`. - The molar flux is ``N = -D grad(c) - D z F / (R T) c grad(phi) + u c``, - with the diffusivity ``D`` corrected for the region. Add a - :class:`~sumfem.models.Accumulation` weighted with the region's - :attr:`~sumfem.electrochemistry.porous.Region.storage` for transient - transport. + The molar flux of the ion of charge :math:`z` and concentration + :math:`c` is + + .. math:: + :label: nernst-planck-flux + + \begin{aligned} + \mathbf{N} &= -D_\text{eff} \nabla c + - \frac{z F}{R T} D_\text{eff} c \nabla \phi + \mathbf{u} c, + \end{aligned} + + with its diffusivity :math:`D` corrected for the region, + :math:`D_\text{eff} = t D` with :math:`t` the region's + :attr:`~sumfem.electrochemistry.porous.Region.transport`, the electrolyte + potential :math:`\phi` and the superficial velocity :math:`\mathbf{u}`. + For transient transport, see :class:`TransientNernstPlanck`. Args: name: Name of the model. @@ -94,3 +106,72 @@ def residual(self) -> ufl.Form: c = self._ion.field dx = c.domain.dx(*self._region.markers) return -ufl.inner(self.flux(), ufl.grad(c.v)) * dx + + +class TransientNernstPlanck(ResidualSum): + r"""Storage and Nernst-Planck transport of one ion. + + The conservation law of the ion in a region of porosity + :math:`\varepsilon`, + + .. math:: + :label: nernst-planck-balance + + \begin{aligned} + \varepsilon \frac{\partial c}{\partial t} + \nabla \cdot \mathbf{N} + &= 0, + \end{aligned} + + with the flux :math:`\mathbf{N}` of :eq:`nernst-planck-flux`. Sources, + e.g. of an :class:`~sumfem.electrochemistry.ElectrodeReaction`, are + added to it as models of their own. It is the sum of + ``"_storage"``, an :class:`~sumfem.models.Accumulation` weighted + with the region's :attr:`~sumfem.electrochemistry.porous.Region.storage`, + and ``"_transport"``, a :class:`NernstPlanck`. + + Args: + name: Name of the sum, and prefix of its summands' names. + ion: The ion, one of the electrolyte's with a field. + electrolyte: The electrolyte of the ion. + potential: The electrolyte potential :math:`\phi`. + region: The region filled with electrolyte. + velocity: The superficial velocity :math:`\mathbf{u}`. Defaults to + zero, without advection. + + Raises: + ValueError: If the ion is not in the electrolyte. + TypeError: If the ion follows from electroneutrality. + """ + + _transport: NernstPlanck + """The transport of the ion.""" + + def __init__( + self, + name: str, + ion: Ion, + electrolyte: Electrolyte, + potential: Field, + region: Region = Region(), + velocity: ufl.core.expr.Expr | None = None, + ): + transport = NernstPlanck( + f"{name}_transport", ion, electrolyte, potential, region, velocity + ) + storage = Accumulation( + f"{name}_storage", + ion.field, + weight=region.storage, + markers=region.markers, + ) + super().__init__([storage, transport], name=name) + self._transport = transport + + @property + def ion(self) -> Ion: + """The transported ion.""" + return self._transport.ion + + def flux(self) -> ufl.core.expr.Expr: + r"""The molar flux :math:`\mathbf{N}`, see :meth:`NernstPlanck.flux`.""" + return self._transport.flux() diff --git a/tests/electrochemistry/half_cell/test_models.py b/tests/electrochemistry/half_cell/test_models.py index 9c23b2e..6bd4020 100644 --- a/tests/electrochemistry/half_cell/test_models.py +++ b/tests/electrochemistry/half_cell/test_models.py @@ -116,8 +116,8 @@ def test_half_cells_have_consistent_units(reference): ) assert_consistent(cell.model) names = [m.name for m in cell.model.components] - assert "cell_transport_V2_channels" in names - assert "cell_storage_H_channels" in names + assert "cell_V2_channels_transport" in names + assert "cell_H_channels_storage" in names assert "cell_electrolyte_channels" in names diff --git a/tests/electrochemistry/test_electrolyte.py b/tests/electrochemistry/test_electrolyte.py index 0990f15..40fb145 100644 --- a/tests/electrochemistry/test_electrolyte.py +++ b/tests/electrochemistry/test_electrolyte.py @@ -12,7 +12,10 @@ MassTransferButlerVolmer, NernstPlanck, PorousElectrode, + Region, + TransientNernstPlanck, ) +from sumfem.models import Accumulation from sumfem.symbolic import Domain, Field, build_system CONCENTRATION = syu.mol / syu.m**3 @@ -127,6 +130,22 @@ def test_transport_needs_a_transported_ion_of_the_electrolyte(domain): NernstPlanck("V2", stranger, electrolyte, phi_l) +def test_transient_transport_is_storage_and_transport(domain): + """The balance of an ion is its storage, weighted with the porosity, and + its Nernst-Planck transport, and needs a transported ion.""" + ions = vanadium(domain) + electrolyte = Electrolyte("e", domain, ions) + _, phi_l = potentials(domain) + region = Region(("Electrode",), porosity=0.7) + balance = TransientNernstPlanck("V2", ions[0], electrolyte, phi_l, region) + storage, transport = balance.components + assert [m.name for m in balance.components] == ["V2_storage", "V2_transport"] + assert isinstance(storage, Accumulation) and isinstance(transport, NernstPlanck) + assert balance.ion is ions[0] + with pytest.raises(TypeError, match="electroneutrality"): + TransientNernstPlanck("SO4", ions[3], electrolyte, phi_l) + + def test_a_vanadium_electrode_has_consistent_units(domain): ions = vanadium(domain) electrolyte = Electrolyte("e", domain, ions, density=1e3, viscosity=1e-3)