Split into nernst planck and transient nernst planck.

docs/api/electrochemistry.md

@@ -32,6 +32,9 @@
 
 .. autoclass:: NernstPlanck
    :members:
+
+.. autoclass:: TransientNernstPlanck
+   :members:
 ```
 
 ## Kinetics

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",

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.

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 <https://doi.org/10.1016/j.jpowsour.2011.08.024>`_, and the
     surface concentrations those of equation 20 of `Aghabarari et al., 2026
     <https://arxiv.org/abs/2610.01840>`_.
@@ -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 <https://arxiv.org/abs/2610.01840>`_.
+
+    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
     <https://doi.org/10.1016/j.electacta.2008.05.067>`_. Unlike
     :class:`MassTransferButlerVolmer`, the transfer does not depend on the
     flow.

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:

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
+    ``"<name>_storage"``, an :class:`~sumfem.models.Accumulation` weighted
+    with the region's :attr:`~sumfem.electrochemistry.porous.Region.storage`,
+    and ``"<name>_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()

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
 
 

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)