From 55565994caeae44d3528e5b953684599c34f1ce1 Mon Sep 17 00:00:00 2001 From: Jorgen Schartum Dokken Date: Thu, 12 Nov 2026 09:35:35 +0000 Subject: [PATCH 1/2] Add protonmembane --- docs/api/electrochemistry.md | 10 + src/sumfem/electrochemistry/membrane.py | 395 ++++++++++++++++++++++++ src/sumfem/models/advection.py | 2 +- tests/electrochemistry/test_membrane.py | 205 ++++++++++++ 4 files changed, 611 insertions(+), 1 deletion(-) create mode 100644 src/sumfem/electrochemistry/membrane.py create mode 100644 tests/electrochemistry/test_membrane.py diff --git a/docs/api/electrochemistry.md b/docs/api/electrochemistry.md index 47fbcca..77ec833 100644 --- a/docs/api/electrochemistry.md +++ b/docs/api/electrochemistry.md @@ -56,6 +56,16 @@ :members: ``` +## Membranes + +```{eval-rst} +.. autoclass:: ProtonExchangeMembrane + :members: + +.. autoclass:: SchloeglFlow + :members: +``` + ## Half-cells ```{eval-rst} diff --git a/src/sumfem/electrochemistry/membrane.py b/src/sumfem/electrochemistry/membrane.py new file mode 100644 index 0000000..29d5257 --- /dev/null +++ b/src/sumfem/electrochemistry/membrane.py @@ -0,0 +1,395 @@ +"""Ion-exchange membranes between the half-cells of a redox flow battery.""" + +from __future__ import annotations + +import typing + +import numpy.typing as npt +import sympy +import sympy.physics.units as syu +import ufl + +from sumfem.models import Accumulation, Conduction, ConservativeAdvection, Diffusion +from sumfem.symbolic import ( + DirichletBC, + Domain, + Field, + Parameter, + Residual, + ResidualSum, + as_expr, +) + +from .electrolyte import CONCENTRATION, DIFFUSIVITY, Electrolyte + +__all__ = ["ProtonExchangeMembrane", "SchloeglFlow"] + +INPUT: typing.TypeAlias = ufl.core.expr.Expr | npt.ArrayLike | None +"""A model input: an expression, a value, or ``None`` for a parameter.""" + +PROPERTY: typing.TypeAlias = ufl.core.expr.Expr | float | sympy.Expr | None +"""A property of the membrane: an expression, used as given, or the default +of the parameter the membrane creates, ``None`` for none.""" + +PERMEABILITY = syu.m**2 +"""Unit of permeabilities.""" + +VISCOSITY = syu.pascal * syu.second +"""Unit of dynamic viscosities.""" + + +class SchloeglFlow(Residual): + r"""Flow of water through a charged membrane, by Schlögl's equation. + + The protons that neutralize the fixed charges of the membrane drag the + water with them, so the velocity of the water has an electro-osmotic + part besides the pressure-driven one: + + .. math:: + :label: schloegl-velocity + + \begin{aligned} + \mathbf{v} &= -\frac{k_\phi}{\mu} \rho \nabla \phi + - \frac{k_p}{\mu} \nabla p, \\ + \nabla \cdot \mathbf{v} &= 0, + \end{aligned} + + with the electrokinetic and hydraulic permeabilities :math:`k_\phi` and + :math:`k_p`, the viscosity :math:`\mu` of water, and the charge density + :math:`\rho = F c_\text{H}` of the mobile protons. The residual is + :math:`-(\mathbf{v}, \nabla q)` on the pressure, as for + :class:`~sumfem.models.Darcy` flow, so the normal velocity is continuous + across the faces of the membrane with a Darcy flow beside it. + + Args: + name: Name of the model. + pressure: The pressure :math:`p`. + potential: The electrolyte potential :math:`\phi`. + charge_density: The charge density :math:`\rho` of the mobile ions. + Defaults to the parameter ``".charge_density"`` in C/m^3. + electrokinetic_permeability: :math:`k_\phi`. Defaults to the + parameter ``".electrokinetic_permeability"`` in m^2. + hydraulic_permeability: :math:`k_p`. Defaults to the parameter + ``".hydraulic_permeability"`` in m^2. + viscosity: :math:`\mu`. Defaults to the parameter + ``".viscosity"`` in Pa s. + markers: Cell markers of the membrane. Empty for the whole domain. + """ + + _pressure: Field + """The pressure.""" + + _potential: ufl.core.expr.Expr + """The electrolyte potential.""" + + _charge_density: ufl.core.expr.Expr + """Charge density of the mobile ions.""" + + _electrokinetic_permeability: ufl.core.expr.Expr + """Electrokinetic permeability.""" + + _hydraulic_permeability: ufl.core.expr.Expr + """Hydraulic permeability.""" + + _viscosity: ufl.core.expr.Expr + """Viscosity of water.""" + + _markers: tuple[str, ...] + """Cell markers of the membrane.""" + + def __init__( + self, + name: str, + pressure: Field, + potential: ufl.core.expr.Expr, + charge_density: INPUT = None, + electrokinetic_permeability: INPUT = None, + hydraulic_permeability: INPUT = None, + viscosity: INPUT = None, + markers: tuple[str, ...] = (), + ): + super().__init__(name) + domain = pressure.domain + self._pressure = pressure + self._potential = as_expr(potential) + self._charge_density = self.input( + "charge_density", + charge_density, + domain=domain, + unit=syu.coulomb / syu.m**3, + doc="Charge density of the mobile ions", + ) + self._electrokinetic_permeability = self.input( + "electrokinetic_permeability", + electrokinetic_permeability, + domain=domain, + unit=PERMEABILITY, + doc="Electrokinetic permeability", + ) + self._hydraulic_permeability = self.input( + "hydraulic_permeability", + hydraulic_permeability, + domain=domain, + unit=PERMEABILITY, + doc="Hydraulic permeability", + ) + self._viscosity = self.input( + "viscosity", viscosity, domain=domain, unit=VISCOSITY, doc="Viscosity" + ) + self._markers = tuple(markers) + + @property + def velocity(self) -> ufl.core.expr.Expr: + r"""The velocity :math:`\mathbf{v}` of :eq:`schloegl-velocity`.""" + mu = self._viscosity + return -self._electrokinetic_permeability / mu * self._charge_density * ( + ufl.grad(self._potential) + ) - self._hydraulic_permeability / mu * ufl.grad(self._pressure) + + def residual(self) -> ufl.Form: + p = self._pressure + dx = p.domain.dx(*self._markers) + return ufl.inner(-self.velocity, ufl.grad(p.v)) * dx + + +class ProtonExchangeMembrane(ResidualSum): + r"""A proton-exchange membrane saturated with water. + + The membrane of `Bernardi and Verbrugge, 1991 + `_, as in the full cell of + `Shah et al., 2008 `_, + Section 3.2. Its fixed charges, of concentration :math:`c_f`, are + neutralized by protons, the only mobile ions, so the proton + concentration is :math:`c_\text{H} = c_f` throughout. The membrane is + then an Ohmic conductor for the electrolyte potential :math:`\phi`, and + the water, of concentration :math:`c_w`, moves with the velocity + :math:`\mathbf{v}` of :eq:`schloegl-velocity`: + + .. math:: + :label: proton-exchange-membrane + + \begin{aligned} + -\nabla \cdot (\kappa_m \nabla \phi) &= 0, \\ + \kappa_m &= \frac{F^2}{R T} D_\text{H} c_f, \\ + \frac{\partial c_w}{\partial t} + + \nabla \cdot (-D_w \nabla c_w + \mathbf{v} c_w) &= 0, + \end{aligned} + + with the proton and water diffusivities :math:`D_\text{H}` and + :math:`D_w` of the membrane. The current density + :math:`-\kappa_m \nabla \phi` uses the sign of + :class:`~sumfem.electrochemistry.ElectrolyteConduction`, so with one + potential across the membrane and the electrodes beside it, the current + is continuous across its faces, and so is the flux of water with one + water field. + + The membrane is the sum of ``"_conduction"``, a + :class:`~sumfem.models.Conduction`; with a water field, + ``"_water_storage"``, an :class:`~sumfem.models.Accumulation`, and + ``"_water_diffusion"``, a :class:`~sumfem.models.Diffusion`; and + with a pressure, ``"_flow"``, a :class:`SchloeglFlow`, and + ``"_water_advection"``, a + :class:`~sumfem.models.ConservativeAdvection` of the water. Without a + pressure, the water does not flow, and only diffuses. + + The membrane creates its constants and properties as parameters named + ``"."``: Faraday's constant, the gas constant and the + temperature, see :class:`~sumfem.electrochemistry.Electrolyte`, and one + for each argument from ``fixed_charge`` on, with the argument as its + default value. + + Args: + name: Name of the membrane, and prefix of its parameters and models. + potential: The electrolyte potential :math:`\phi`, e.g. one field + over the membrane and the electrodes. + markers: Cell markers of the membrane. + water: The water concentration :math:`c_w`, or ``None`` without + water transport. + pressure: The pressure :math:`p` of :eq:`schloegl-velocity`, or + ``None`` without flow. + fixed_charge: The concentration :math:`c_f` of fixed charges, in + mol/m^3 or with units. + proton_diffusivity: :math:`D_\text{H}` in the membrane, in m^2/s or + with units. + water_diffusivity: :math:`D_w` in the membrane, in m^2/s or with + units. + electrokinetic_permeability: :math:`k_\phi`, in m^2 or with units. + hydraulic_permeability: :math:`k_p`, in m^2 or with units. + viscosity: The viscosity :math:`\mu` of water, in Pa s or with units. + temperature: The temperature :math:`T`, in K. + """ + + _potential: Field + """The electrolyte potential.""" + + _water: Field | None + """The water concentration, if transported.""" + + _flow: SchloeglFlow | None + """The flow of water, with a pressure.""" + + _constants: Electrolyte + """Faraday's constant, the gas constant and the temperature.""" + + _fixed_charge: ufl.core.expr.Expr + """Concentration of the fixed charges.""" + + _markers: tuple[str, ...] + """Cell markers of the membrane.""" + + def __init__( + self, + name: str, + potential: Field, + markers: tuple[str, ...], + water: Field | None = None, + pressure: Field | None = None, + fixed_charge: PROPERTY = None, + proton_diffusivity: PROPERTY = None, + water_diffusivity: PROPERTY = None, + electrokinetic_permeability: PROPERTY = None, + hydraulic_permeability: PROPERTY = None, + viscosity: PROPERTY = None, + temperature: float = 298.15, + ): + domain = potential.domain + self._potential = potential + self._water = water + self._markers = tuple(markers) + self._constants = Electrolyte(name, domain, temperature=temperature) + e = self._constants + self._fixed_charge = _property( + f"{name}.fixed_charge", + fixed_charge, + domain, + unit=CONCENTRATION, + doc="Concentration of the fixed charges", + ) + proton = _property( + f"{name}.proton_diffusivity", + proton_diffusivity, + domain, + unit=DIFFUSIVITY, + doc="Proton diffusivity in the membrane", + ) + conductivity = e.F**2 / (e.R * e.T) * proton * self._fixed_charge + models: list[Residual] = [ + Conduction(f"{name}_conduction", potential, conductivity, self._markers) + ] + self._flow = None + if pressure is not None: + self._flow = SchloeglFlow( + f"{name}_flow", + pressure, + potential, + charge_density=e.F * self._fixed_charge, + electrokinetic_permeability=_property( + f"{name}.electrokinetic_permeability", + electrokinetic_permeability, + domain, + unit=PERMEABILITY, + doc="Electrokinetic permeability", + ), + hydraulic_permeability=_property( + f"{name}.hydraulic_permeability", + hydraulic_permeability, + domain, + unit=PERMEABILITY, + doc="Hydraulic permeability", + ), + viscosity=_property( + f"{name}.viscosity", + viscosity, + domain, + unit=VISCOSITY, + doc="Viscosity of water", + ), + markers=self._markers, + ) + models.append(self._flow) + if water is not None: + diffusivity = _property( + f"{name}.water_diffusivity", + water_diffusivity, + domain, + unit=DIFFUSIVITY, + doc="Water diffusivity in the membrane", + ) + models.append( + Accumulation(f"{name}_water_storage", water, markers=self._markers) + ) + models.append( + Diffusion(f"{name}_water_diffusion", water, diffusivity, self._markers) + ) + if self._flow is not None: + models.append( + ConservativeAdvection( + f"{name}_water_advection", + water, + self.velocity, + markers=self._markers, + ) + ) + super().__init__(models, name=name) + + @property + def potential(self) -> Field: + """The electrolyte potential.""" + return self._potential + + @property + def fixed_charge(self) -> ufl.core.expr.Expr: + """The concentration of the fixed charges.""" + return self._fixed_charge + + @property + def velocity(self) -> ufl.core.expr.Expr: + """The velocity of the water, zero without a pressure.""" + if self._flow is None: + return ufl.zero(self._potential.domain.geometric_dimension) + return self._flow.velocity + + def proton_conditions( + self, protons: typing.Sequence[Field], faces: str | tuple[str, ...] + ) -> list[DirichletBC]: + """The proton concentration of the membrane at its faces. + + Electroneutrality holds the protons at the concentration of the + fixed charges, also where the electrodes meet the membrane. + + Args: + protons: The proton concentrations of the electrodes beside the + membrane. + faces: Facet markers of the faces of the membrane. + + Returns: + The conditions :math:`c_\\text{H} = c_f` on the faces. + """ + return [DirichletBC(c, faces, self._fixed_charge) for c in protons] + + +def _property( + name: str, + value: PROPERTY, + domain: Domain, + *, + unit: syu.Quantity | sympy.Expr, + doc: str, +) -> ufl.core.expr.Expr: + """A property of the membrane: an expression, or a new parameter. + + Args: + name: The full name of the parameter, ``"."``. + value: An expression, e.g. a shared parameter, used as given, or the + default value of the new parameter, in SI units or with units; + ``None`` for a parameter without a default. + domain: The domain of the parameter. + unit: The unit of the parameter. + doc: A description of the parameter. + + Returns: + The expression or the parameter. + """ + if isinstance(value, ufl.core.expr.Expr): + return value + return Parameter(name, domain, default=value, unit=unit, doc=doc) diff --git a/src/sumfem/models/advection.py b/src/sumfem/models/advection.py index 0f10cb7..209d1c1 100644 --- a/src/sumfem/models/advection.py +++ b/src/sumfem/models/advection.py @@ -11,7 +11,7 @@ class Advection(Residual): - r"""Advection of a field by a velocity, in a form of a subclass. + r"""Advection of a field by a velocity The velocity is any expression, e.g. a field solved by a flow model or data transferred from another mesh. The advected quantity is the field, diff --git a/tests/electrochemistry/test_membrane.py b/tests/electrochemistry/test_membrane.py new file mode 100644 index 0000000..5a37956 --- /dev/null +++ b/tests/electrochemistry/test_membrane.py @@ -0,0 +1,205 @@ +"""The proton-exchange membrane of Bernardi and Verbrugge, with Schlögl's flow.""" + +from mpi4py import MPI + +import dolfinx.mesh +import numpy as np +import pytest +import sympy.physics.units as syu +import ufl +from basix.ufl import element + +from sumfem.discrete import DomainMaps, MeshData, discretize +from sumfem.electrochemistry import ProtonExchangeMembrane +from sumfem.solvers import BDF, Solver +from sumfem.symbolic import ( + FARADAY, + GAS_CONSTANT, + DirichletBC, + Domain, + Field, + ResidualSum, + UFLResidual, + build_system, +) +from sumfem.symbolic.si import to_si + +CELLS = ["interval", "triangle", "quadrilateral", "tetrahedron", "hexahedron"] + +DOMAIN = "box" +"""The name of the symbolic domain of both formulations, see ``compare``.""" + +VALUES = { + "fixed_charge": 1200.0, + "proton_diffusivity": 1.4e-9, + "water_diffusivity": 5.75e-10, + "electrokinetic_permeability": 1.13e-19, + "hydraulic_permeability": 1.58e-18, + "viscosity": 1e-3, +} +"""The membrane of Shah et al. (2008), Tables 4 and 5.""" + +F = float(to_si(FARADAY, syu.coulomb / syu.mol)) +R = float(to_si(GAS_CONSTANT, syu.joule / (syu.mol * syu.kelvin))) +T = 298.15 + + +def steady(model): + """The membrane without the storage of water.""" + return ResidualSum( + [m for m in model.components if not m.name.endswith("storage")], + name="steady", + ) + + +def membrane(domain, water=True, pressure=True): + phi = Field("phi", domain, ("Lagrange", 1), unit=syu.volt) + c_w = Field("c_w", domain, ("Lagrange", 1), unit=syu.mol / syu.m**3, scale=4e4) + p = Field("p", domain, ("Lagrange", 1), unit=syu.pascal, scale=1e3) + model = ProtonExchangeMembrane( + "membrane", + phi, + (), + water=c_w if water else None, + pressure=p if pressure else None, + **VALUES, + ) + return model, phi, c_w, p + + +@pytest.mark.parametrize("cell", CELLS) +def test_the_membrane_against_plain_ufl(cell, box, compare): + """Ohmic conduction with the conductivity of the fixed protons, Schlögl's + flow and the transport of water by it.""" + gdim = ufl.Cell(cell).topological_dimension + mesh = ufl.Mesh(element("Lagrange", cell, 1, shape=(gdim,))) + V = ufl.FunctionSpace(mesh, element("Lagrange", cell, 1)) + W = ufl.MixedFunctionSpace(V, V, V) + phi, c_w, p = (ufl.Coefficient(V) for _ in range(3)) + v_phi, w, q = ufl.TestFunctions(W) + Fc, Rc, Tc, cf, DH, Dw, kphi, kp, mu = (ufl.Constant(mesh) for _ in range(9)) + kappa = Fc**2 / (Rc * Tc) * DH * cf + v = -kphi / mu * Fc * cf * ufl.grad(phi) - kp / mu * ufl.grad(p) + residual = ( + kappa * ufl.inner(ufl.grad(phi), ufl.grad(v_phi)) + - ufl.inner(v, ufl.grad(q)) + - ufl.inner(-Dw * ufl.grad(c_w) + v * c_w, ufl.grad(w)) + ) * ufl.dx + names = { + DOMAIN: mesh, + **dict(phi=phi, c_w=c_w, p=p), + **dict(F=Fc, R=Rc, T=Tc, cf=cf, DH=DH, Dw=Dw, kphi=kphi, kp=kp, mu=mu), + } + plain = UFLResidual("plain", residual, [phi, c_w, p], names=names) + + domain = Domain(DOMAIN, cell) + model, *fields = membrane(domain) + values = { + "plain.F": F, + "plain.R": R, + "plain.T": T, + "plain.cf": VALUES["fixed_charge"], + "plain.DH": VALUES["proton_diffusivity"], + "plain.Dw": VALUES["water_diffusivity"], + "plain.kphi": VALUES["electrokinetic_permeability"], + "plain.kp": VALUES["hydraulic_permeability"], + "plain.mu": VALUES["viscosity"], + } + compare( + plain, + plain.unknowns, + steady(model), + fields, + box(cell), + [(values, {})], + ) + + +def test_the_membrane_has_consistent_units(): + domain = Domain("cell", "interval") + model, *_ = membrane(domain) + build_system(model, check_units=True) + build_system(model, substitutions=BDF(1).substitutions(model), check_units=True) + names = [m.name for m in model.components] + assert names == [ + "membrane_conduction", + "membrane_flow", + "membrane_water_storage", + "membrane_water_diffusion", + "membrane_water_advection", + ] + still, *_ = membrane(domain, pressure=False) + assert [m.name for m in still.components] == [ + "membrane_conduction", + "membrane_water_storage", + "membrane_water_diffusion", + ] + bare, *_ = membrane(domain, water=False, pressure=False) + assert [m.name for m in bare.components] == ["membrane_conduction"] + + +def test_a_membrane_between_two_potentials_and_pressures(): + """Across a membrane of thickness L, with the potential, the pressure and + the water held at both faces, the potential and the pressure are linear, + so the current density and Schlögl's velocity are uniform, and the water + follows, at the nodes, the steady advection-diffusion profile + ``c(x) = c_0 + (c_L - c_0) (exp(v x / D) - 1) / (exp(v L / D) - 1)``. + """ + L, cells = 180e-6, 400 + mesh = dolfinx.mesh.create_interval(MPI.COMM_WORLD, cells, [0.0, L]) + marked = MeshData.from_locators( + mesh, + { + 0: { + "Negative": lambda x: np.isclose(x[0], 0.0), + "Positive": lambda x: np.isclose(x[0], L), + } + }, + boundary=True, + ) + domain = Domain.from_mesh("membrane", mesh, length_scale=L) + model, phi, c_w, p = membrane(domain) + dphi, dp, c0, cL = 0.01, -5e4, 4.2e4, 4.1e4 + bcs = [ + DirichletBC(phi, "Negative", 0.0), + DirichletBC(phi, "Positive", dphi), + DirichletBC(p, "Negative", 0.0), + DirichletBC(p, "Positive", dp), + DirichletBC(c_w, "Negative", c0), + DirichletBC(c_w, "Positive", cL), + ] + # The pressure's equation is some fifteen orders of magnitude smaller + # than the potential's, so the equations are solved in dimensionless form + system = build_system(steady(model), check_units=True, factorize_equations=True) + disc = discretize([system, *bcs], DomainMaps({domain: marked})) + Solver(disc[system], bcs).solve() + + kappa = F**2 / (R * T) * VALUES["proton_diffusivity"] * VALUES["fixed_charge"] + mu = VALUES["viscosity"] + v = ( + -VALUES["electrokinetic_permeability"] / mu * F * VALUES["fixed_charge"] * dphi + - VALUES["hydraulic_permeability"] / mu * dp + ) / L + current = disc.assemble_scalar(-kappa * phi.dx(0) * domain.dx()) / L + assert current == pytest.approx(-kappa * dphi / L, rel=1e-10) + velocity = disc.assemble_scalar(model.velocity[0] * domain.dx()) / L + assert velocity == pytest.approx(v, rel=1e-10) + + # Linear elements are as good as exact at the nodes of this 1D problem; + # compared there, the exponential profile needs no quadrature + D = VALUES["water_diffusivity"] + water = disc.function(c_w) + owned = water.function_space.dofmap.index_map.size_local + x = water.function_space.tabulate_dof_coordinates()[:owned, 0] + exact = c0 + (cL - c0) * np.expm1(v * x / D) / np.expm1(v * L / D) + error = np.abs(water.x.array[:owned] - exact).max(initial=0.0) + assert mesh.comm.allreduce(error, op=MPI.MAX) < 1e-9 * abs(cL - c0) + + +def test_the_membrane_holds_the_protons_at_its_faces(): + domain = Domain("cell", "interval") + model, *_ = membrane(domain) + protons = [Field(f"c_H_{side}", domain, ("Lagrange", 1)) for side in "ab"] + bcs = model.proton_conditions(protons, ("Negative", "Positive")) + assert [bc.field for bc in bcs] == protons + assert all(bc.value is model.fixed_charge for bc in bcs) From e78991033c31b13d37c743c53808f340875132e3 Mon Sep 17 00:00:00 2001 From: Jorgen Schartum Dokken Date: Thu, 12 Nov 2026 09:35:48 +0000 Subject: [PATCH 2/2] Add missing import --- src/sumfem/electrochemistry/__init__.py | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/src/sumfem/electrochemistry/__init__.py b/src/sumfem/electrochemistry/__init__.py index b7f8da7..79fd35a 100644 --- a/src/sumfem/electrochemistry/__init__.py +++ b/src/sumfem/electrochemistry/__init__.py @@ -10,6 +10,8 @@ :class:`MassTransferButlerVolmer` and :class:`PistonVelocityButlerVolmer`, on :class:`ElectrodeReaction`, with the Butler-Volmer equation :func:`~sumfem.models.butler_volmer`. +- Membranes: :class:`ProtonExchangeMembrane`, with the flow of water of + :class:`SchloeglFlow`. - Half-cells of a redox flow battery: :class:`RedoxHalfCell`, negative or positive, of any :class:`RedoxCouple` and electrolyte. @@ -28,6 +30,7 @@ MassTransferButlerVolmer, PistonVelocityButlerVolmer, ) +from .membrane import ProtonExchangeMembrane, SchloeglFlow from .porous import PorousElectrode, Region from .transport import NernstPlanck, TransientNernstPlanck @@ -39,6 +42,8 @@ "ElectrolyteConduction", "NernstPlanck", "TransientNernstPlanck", + "ProtonExchangeMembrane", + "SchloeglFlow", "RedoxCouple", "RedoxHalfCell", "ButlerVolmer",