From 31b589dc481c73ed94df5e6fae49895d77cb2936 Mon Sep 17 00:00:00 2001 From: Jorgen Schartum Dokken Date: Thu, 12 Nov 2026 09:42:30 +0000 Subject: [PATCH 1/3] Resolve https://github.com/scientificcomputing/sumfem/issues/20 --- .../electrochemistry/vanadium_half_cell_1d.py | 8 +-- .../electrochemistry/vanadium_half_cell_2d.py | 4 +- .../electrochemistry/vanadium_half_cell_3d.py | 4 +- src/sumfem/electrochemistry/electrolyte.py | 42 ++++++----- src/sumfem/electrochemistry/half_cell.py | 2 +- src/sumfem/electrochemistry/kinetics.py | 9 +-- tests/electrochemistry/half_cell/conftest.py | 4 +- .../electrochemistry/half_cell/test_models.py | 4 +- .../test_electrochemistry_from_ufl.py | 28 +++----- tests/electrochemistry/test_electrolyte.py | 70 ++++++++++++++++--- 10 files changed, 110 insertions(+), 65 deletions(-) diff --git a/examples/electrochemistry/vanadium_half_cell_1d.py b/examples/electrochemistry/vanadium_half_cell_1d.py index 34f5513..18960b9 100644 --- a/examples/electrochemistry/vanadium_half_cell_1d.py +++ b/examples/electrochemistry/vanadium_half_cell_1d.py @@ -139,8 +139,8 @@ # # The ions of the electrolyte are given with their charges and diffusivities, # in their units, which sumfem converts to SI, see {ref}`units-checking`. The -# vanadium ions take the positive part of their concentrations in the -# kinetics and the conductivity, as they may be depleted. The concentrations +# kinetics take the positive part of the vanadium concentrations, as they may +# be depleted. The concentrations # of the reservoir name the ions that are transported, each with a field and # the equations of {eq}`half-cell`. Sulfate has none: electroneutrality, # the last equation of {eq}`half-cell`, with the charges of the ions, @@ -171,8 +171,8 @@ diffusivity = syu.m**2 / syu.second concentration = syu.mol / syu.m**3 ions = ( - Ion("V2", 2, 8.10e-12 * diffusivity, positive=True), - Ion("V3", 3, 1.65e-11 * diffusivity, positive=True), + Ion("V2", 2, 8.10e-12 * diffusivity), + Ion("V3", 3, 1.65e-11 * diffusivity), Ion("H", 1, 9.312e-9 * diffusivity), Ion("HSO4", -1, 1.33e-9 * diffusivity), Ion("SO4", -2, 1.065e-9 * diffusivity), diff --git a/examples/electrochemistry/vanadium_half_cell_2d.py b/examples/electrochemistry/vanadium_half_cell_2d.py index 1c7609e..fb592b9 100644 --- a/examples/electrochemistry/vanadium_half_cell_2d.py +++ b/examples/electrochemistry/vanadium_half_cell_2d.py @@ -212,8 +212,8 @@ diffusivity = syu.m**2 / syu.second concentration = syu.mol / syu.m**3 ions = ( - Ion("V2", 2, 8.10e-12 * diffusivity, positive=True), - Ion("V3", 3, 1.65e-11 * diffusivity, positive=True), + Ion("V2", 2, 8.10e-12 * diffusivity), + Ion("V3", 3, 1.65e-11 * diffusivity), Ion("H", 1, 9.312e-9 * diffusivity), Ion("HSO4", -1, 1.33e-9 * diffusivity), Ion("SO4", -2, 1.065e-9 * diffusivity), diff --git a/examples/electrochemistry/vanadium_half_cell_3d.py b/examples/electrochemistry/vanadium_half_cell_3d.py index e088b40..35edf7b 100644 --- a/examples/electrochemistry/vanadium_half_cell_3d.py +++ b/examples/electrochemistry/vanadium_half_cell_3d.py @@ -322,8 +322,8 @@ def flow_through_electrode(sizes: tuple[float, float]) -> None: diffusivity = syu.m**2 / syu.second concentration = syu.mol / syu.m**3 ions = ( - Ion("V2", 2, 8.10e-12 * diffusivity, positive=True), - Ion("V3", 3, 1.65e-11 * diffusivity, positive=True), + Ion("V2", 2, 8.10e-12 * diffusivity), + Ion("V3", 3, 1.65e-11 * diffusivity), Ion("H", 1, 9.312e-9 * diffusivity), Ion("HSO4", -1, 1.33e-9 * diffusivity), Ion("SO4", -2, 1.065e-9 * diffusivity), diff --git a/src/sumfem/electrochemistry/electrolyte.py b/src/sumfem/electrochemistry/electrolyte.py index ff3929a..32bd24b 100644 --- a/src/sumfem/electrochemistry/electrolyte.py +++ b/src/sumfem/electrochemistry/electrolyte.py @@ -39,16 +39,12 @@ class Ion: default value of the parameter the electrolyte creates. concentration: The field of a transported ion. ``None`` for the one ion whose concentration follows from electroneutrality. - positive: Use a smooth positive part of the concentration in the - conductivity and in reaction kinetics. Use it for ions that can - be depleted, such as the active species. """ name: str charge: float diffusivity: float | sympy.Expr concentration: Field | None = None - positive: bool = False @property def field(self) -> Field: @@ -278,26 +274,30 @@ def __getitem__(self, name: str) -> Ion: return self._ions[name] def positive_part(self, c: ufl.core.expr.Expr) -> ufl.core.expr.Expr: - """The smooth positive part ``(c + sqrt(c^2 + floor^2)) / 2``.""" + """The smooth positive part ``(c + sqrt(c^2 + floor^2)) / 2``. + + Reaction kinetics take it of the concentrations they react, see + :class:`~sumfem.electrochemistry.ButlerVolmer`, so that a + concentration that undershoots zero neither reverses a reaction nor + enters a logarithm. + + Args: + c: A concentration. + """ return 0.5 * (c + ufl.sqrt(c * c + self._floor * self._floor)) - def concentration( - self, name: str, positive: bool | None = None - ) -> ufl.core.expr.Expr: - """Concentration of an ion, as used in conductivity and kinetics. + def concentration(self, name: str) -> ufl.core.expr.Expr: + """Concentration of an ion, with its sign. That is the field of a transported ion, and the expression given by - electroneutrality for the ion without one. + electroneutrality for the ion without one. The conductivity uses it + as it is, as Nernst-Planck migration does, so the current is the + charge the ions carry, see :meth:`current_density`. Args: name: Name of the ion. - positive: Take the positive part of the concentration. Defaults - to :attr:`Ion.positive` of the ion. """ - concentration = self._concentrations[name] - if self._ions[name].positive if positive is None else positive: - return self.positive_part(concentration) - return concentration + return self._concentrations[name] def diffusivity(self, name: str) -> Parameter: """Diffusivity of an ion in free electrolyte, the parameter ``D_``. @@ -331,9 +331,7 @@ def diffusion_current( correction: See :meth:`conductivity`. """ total = sum( - s.charge - * self.diffusivity(s.name) - * ufl.grad(self.concentration(s.name, positive=False)) + s.charge * self.diffusivity(s.name) * ufl.grad(self.concentration(s.name)) for s in self._ions.values() ) return self._F * correction * total @@ -345,6 +343,12 @@ def current_density( ) -> ufl.core.expr.Expr: """Current density ``-kappa grad(phi) - F sum_j z_j D_j grad(c_j)``. + It is the charge the ions carry, ``F sum_j z_j N_j``, with the + Nernst-Planck fluxes ``N_j`` of every ion, also the one given by + electroneutrality: the conductivity and the migration use the same + concentrations, and advection carries no charge, as + ``sum_j z_j c_j = 0``. + Args: potential: The electrolyte potential. correction: See :meth:`conductivity`. diff --git a/src/sumfem/electrochemistry/half_cell.py b/src/sumfem/electrochemistry/half_cell.py index 235cf6d..3350732 100644 --- a/src/sumfem/electrochemistry/half_cell.py +++ b/src/sumfem/electrochemistry/half_cell.py @@ -814,7 +814,7 @@ def minimum_balanced_concentration(self, discretization: Discretization) -> floa assert e.balanced is not None space = discretization.function_space(next(iter(self._concentrations.values()))) balanced = dolfinx.fem.Function(space) - expression = discretization[e.concentration(e.balanced, positive=False)] + expression = discretization[e.concentration(e.balanced)] balanced.interpolate( dolfinx.fem.Expression(expression, space.element.interpolation_points) ) diff --git a/src/sumfem/electrochemistry/kinetics.py b/src/sumfem/electrochemistry/kinetics.py index bfda662..762556d 100644 --- a/src/sumfem/electrochemistry/kinetics.py +++ b/src/sumfem/electrochemistry/kinetics.py @@ -371,10 +371,7 @@ def __init__( ) 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() - } + return {electrolyte.concentration(ion): nu for ion, nu in ions.items()} super().__init__( name, @@ -433,8 +430,8 @@ def _bulk(self) -> tuple[ufl.core.expr.Expr, ufl.core.expr.Expr]: """Positive parts of the bulk concentrations of the reduced and oxidized ion.""" e = self._electrolyte return ( - e.concentration(self._reduced.name, positive=True), - e.concentration(self._oxidized.name, positive=True), + e.positive_part(e.concentration(self._reduced.name)), + e.positive_part(e.concentration(self._oxidized.name)), ) def rate_constants(self) -> tuple[ufl.core.expr.Expr, ufl.core.expr.Expr]: diff --git a/tests/electrochemistry/half_cell/conftest.py b/tests/electrochemistry/half_cell/conftest.py index d37e0e7..fde4214 100644 --- a/tests/electrochemistry/half_cell/conftest.py +++ b/tests/electrochemistry/half_cell/conftest.py @@ -366,8 +366,8 @@ def darcy_velocity(marked: MeshData, p): CONCENTRATION = syu.mol / syu.m**3 IONS = ( - Ion("V2", 2.0, 8.10e-12 * DIFFUSIVITY, positive=True), - Ion("V3", 3.0, 1.65e-11 * DIFFUSIVITY, positive=True), + Ion("V2", 2.0, 8.10e-12 * DIFFUSIVITY), + Ion("V3", 3.0, 1.65e-11 * DIFFUSIVITY), Ion("H", 1.0, 9.312e-9 * DIFFUSIVITY), Ion("HSO4", -1.0, 1.33e-9 * DIFFUSIVITY), Ion("SO4", -2.0, 1.065e-9 * DIFFUSIVITY), diff --git a/tests/electrochemistry/half_cell/test_models.py b/tests/electrochemistry/half_cell/test_models.py index 548f6c4..e238113 100644 --- a/tests/electrochemistry/half_cell/test_models.py +++ b/tests/electrochemistry/half_cell/test_models.py @@ -257,8 +257,8 @@ def test_the_reservoir_lacks_exactly_one_ion(reference, missing): 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("VO", 2, 3.9e-10 * DIFFUSIVITY_UNIT), + Ion("VO2", 1, 3.9e-10 * DIFFUSIVITY_UNIT), Ion("H", 1, 9.312e-9 * DIFFUSIVITY_UNIT), Ion("HSO4", -1, 1.33e-9 * DIFFUSIVITY_UNIT), Ion("H2O", 0, 2.3e-9 * DIFFUSIVITY_UNIT), diff --git a/tests/electrochemistry/test_electrochemistry_from_ufl.py b/tests/electrochemistry/test_electrochemistry_from_ufl.py index f32181a..86f7d6d 100644 --- a/tests/electrochemistry/test_electrochemistry_from_ufl.py +++ b/tests/electrochemistry/test_electrochemistry_from_ufl.py @@ -52,7 +52,7 @@ def test_conduction_and_transport_in_a_porous_region(cell, box, compare): electroneutrality, and Nernst-Planck transport of the three, with advection of V2+, in the pores of an electrode on the left half. - The conductivity takes the positive parts of the vanadium ions, and every + The conductivity takes the concentrations as they are, and every diffusivity Bruggeman's correction of the porosity. """ gdim = ufl.Cell(cell).topological_dimension @@ -61,19 +61,14 @@ def test_conduction_and_transport_in_a_porous_region(cell, box, compare): W = ufl.MixedFunctionSpace(V, V, V, V) phi, c2, c3, h = (ufl.Coefficient(V) for _ in range(4)) v_phi, w2, w3, wh = ufl.TestFunctions(W) - F, R, T, floor, eps, D2, D3, DH, DS = (ufl.Constant(mesh) for _ in range(9)) + F, R, T, eps, D2, D3, DH, DS = (ufl.Constant(mesh) for _ in range(8)) sulfate = -(2 * c2 + 3 * c3 + h) / -2 correction = eps**1.5 conductivity = ( F**2 / (R * T) * correction - * ( - 4 * D2 * positive_part(c2, floor) - + 9 * D3 * positive_part(c3, floor) - + DH * h - + 4 * DS * sulfate - ) + * (4 * D2 * c2 + 9 * D3 * c3 + DH * h + 4 * DS * sulfate) ) diffusion = ( F @@ -103,7 +98,7 @@ def flux(c, z, D, advected=False): names = { DOMAIN: mesh, **{n: c for n, c in zip(("phi", "c2", "c3", "h"), (phi, c2, c3, h))}, - **dict(F=F, R=R, T=T, floor=floor, eps=eps, D2=D2, D3=D3, DH=DH, DS=DS), + **dict(F=F, R=R, T=T, eps=eps, D2=D2, D3=D3, DH=DH, DS=DS), } plain = UFLResidual( "plain", @@ -117,12 +112,12 @@ def flux(c, z, D, advected=False): potential = Field("phi", domain, ("Lagrange", 1)) fields = [Field(name, domain, ("Lagrange", 1)) for name in ("c2", "c3", "h")] ions = [ - Ion("V2", 2, 2.4e-10, fields[0], positive=True), - Ion("V3", 3, 2.4e-10, fields[1], positive=True), + Ion("V2", 2, 2.4e-10, fields[0]), + Ion("V3", 3, 2.4e-10, fields[1]), Ion("H", 1, 9.3e-9, fields[2]), Ion("SO4", -2, 1.065e-9), ] - electrolyte = Electrolyte("e", domain, ions, concentration_floor=0.1) + electrolyte = Electrolyte("e", domain, ions) region = PorousElectrode("el", domain, ("LeftHalf",), porosity=0.8).region model = ElectrolyteConduction("charge", potential, electrolyte, region) model += NernstPlanck("V2", ions[0], electrolyte, potential, region, velocity=u) @@ -133,7 +128,6 @@ def flux(c, z, D, advected=False): "F": electrolyte.F.default_si, "R": electrolyte.R.default_si, "T": electrolyte.T.default_si, - "floor": 0.1, "eps": 0.8, "D2": 2.4e-10, "D3": 2.4e-10, @@ -198,8 +192,8 @@ def test_kinetics_of_a_redox_couple(cell, mass_transfer, box, compare): potentials = [Field(name, domain, ("Lagrange", 1)) for name in ("phi_s", "phi_l")] fields = [Field(name, domain, ("Lagrange", 1)) for name in ("c2", "c3")] ions = [ - Ion("V2", 2, 2.4e-10, fields[0], positive=True), - Ion("V3", 3, 2.4e-10, fields[1], positive=True), + Ion("V2", 2, 2.4e-10, fields[0]), + Ion("V3", 3, 2.4e-10, fields[1]), ] electrolyte = Electrolyte( "e", domain, ions, density=1.3e3, viscosity=5e-3, concentration_floor=0.1 @@ -286,8 +280,8 @@ def test_positive_couple_with_protons_water_and_piston_velocities(cell, box, com potentials = [Field(name, domain, ("Lagrange", 1)) for name in ("phi_s", "phi_l")] fields = [Field(name, domain, ("Lagrange", 1)) for name in ("c4", "c5", "ch", "cw")] ions = [ - Ion("VO", 2, 3.9e-10, fields[0], positive=True), - Ion("VO2", 1, 3.9e-10, fields[1], positive=True), + Ion("VO", 2, 3.9e-10, fields[0]), + Ion("VO2", 1, 3.9e-10, fields[1]), Ion("H", 1, 9.3e-9, fields[2]), Ion("H2O", 0, 2.3e-9, fields[3]), Ion("SO4", -2, 1.1e-9), diff --git a/tests/electrochemistry/test_electrolyte.py b/tests/electrochemistry/test_electrolyte.py index 40fb145..7e48904 100644 --- a/tests/electrochemistry/test_electrolyte.py +++ b/tests/electrochemistry/test_electrolyte.py @@ -3,6 +3,7 @@ import pytest import sympy.physics.units as syu +import ufl from sumfem.discrete import DomainMaps, discretize from sumfem.electrochemistry import ( @@ -31,8 +32,8 @@ def vanadium(domain): for name in ("c2", "c3", "h") ) return [ - Ion("V2", 2, 2.4e-10, c2, positive=True), - Ion("V3", 3, 2.4e-10, c3, positive=True), + Ion("V2", 2, 2.4e-10, c2), + Ion("V3", 3, 2.4e-10, c3), Ion("H", 1, 9.3e-9, h), Ion("SO4", -2, 1.065e-9), ] @@ -56,7 +57,7 @@ def uniform(domain, square, values): ions = vanadium(domain) electrolyte = Electrolyte("e", domain, ions) _, phi_l = potentials(domain) - # The conduction holds the parameters of every ion, and the floor + # The conduction holds the parameters of every ion model = ElectrolyteConduction("charge", phi_l, electrolyte) for ion in ions[:3]: model += NernstPlanck(ion.name, ion, electrolyte, phi_l) @@ -89,15 +90,64 @@ def test_the_balanced_ion_follows_from_electroneutrality(domain, square): assert sulfate == pytest.approx(-charge / electrolyte["SO4"].charge) -def test_positive_ions_take_the_positive_part(domain, square): +def test_concentrations_keep_their_sign(domain, square): + """The electrolyte gives a concentration that undershoots zero as it is, + and the balanced one from it.""" electrolyte, disc = uniform(domain, square, {**VALUES, "V2": -1.0}) - assert mean(disc, domain, electrolyte.concentration("V2")) == pytest.approx( - 0.0, abs=1e-12 + assert mean(disc, domain, electrolyte.concentration("V2")) == pytest.approx(-1.0) + sulfate = (2 * -1.0 + 3 * VALUES["V3"] + VALUES["H"]) / 2 + assert mean(disc, domain, electrolyte.concentration("SO4")) == pytest.approx( + sulfate ) - unclipped = electrolyte.concentration("V2", positive=False) - assert mean(disc, domain, unclipped) == pytest.approx(-1.0) - # An ion that is not positive keeps its sign unless asked - assert mean(disc, domain, electrolyte.concentration("H")) == pytest.approx(3.0) + + +def test_the_current_is_the_charge_the_ions_carry(domain, square): + """The electrolyte current equals ``F sum_j z_j N_j``, with the + Nernst-Planck fluxes of every ion, sulfate's from electroneutrality, in + a porous region with advection, also where V2+ undershoots zero. + + A regression test: the conductivity once took the positive part of the + vanadium concentrations while migration did not, so conduction did not + carry the charge the ions moved where a concentration was negative. + """ + ions = vanadium(domain) + electrolyte = Electrolyte("e", domain, ions) + _, phi_l = potentials(domain) + region = Region(porosity=0.6) + velocity = ufl.as_vector((1e-3, 2e-3)) + transport = [ + NernstPlanck(ion.name, ion, electrolyte, phi_l, region, velocity) + for ion in ions[:3] + ] + model = ElectrolyteConduction("charge", phi_l, electrolyte, region) + for summand in transport: + model += summand + disc = discretize([model], DomainMaps({domain: square()})) + profiles = { + ions[0].field: lambda x: -1.0 + 3.0 * x[0], + ions[1].field: lambda x: 2.0 + x[1], + ions[2].field: lambda x: 3.0 + x[0] * x[1], + phi_l: lambda x: x[0] ** 2 + 0.5 * x[1], + } + for field, profile in profiles.items(): + disc.function(field).interpolate(profile) + + e, sulfate = electrolyte, ions[3] + c = e.concentration(sulfate.name) + D = region.transport * e.diffusivity(sulfate.name) + migration = D * sulfate.charge * e.F * c / (e.R * e.T) * ufl.grad(phi_l) + fluxes = [m.flux() for m in transport] + fluxes.append(-D * ufl.grad(c) - migration + velocity * c) + charges = [e.F * ion.charge * N for ion, N in zip(ions, fluxes)] + current = e.current_density(phi_l, region.transport) + + # The charge fluxes of the ions largely cancel, so the round-off is + # relative to them, not to the current + difference = current - sum(charges) + error = mean(disc, domain, ufl.inner(difference, difference)) + size = sum(mean(disc, domain, ufl.inner(q, q)) for q in charges) + assert mean(disc, domain, ufl.inner(current, current)) > 0 + assert error <= 1e-24 * size @pytest.mark.parametrize( From 314e4145817987f73a3dafd59671483bc6d41a8f Mon Sep 17 00:00:00 2001 From: Jorgen Schartum Dokken Date: Thu, 12 Nov 2026 09:58:20 +0000 Subject: [PATCH 2/3] Temperature as a known or unknown expression/coefficient/field --- .../electrochemistry/vanadium_half_cell_1d.py | 10 +- src/sumfem/electrochemistry/charge.py | 27 ++--- src/sumfem/electrochemistry/electrolyte.py | 113 +++++++++++++----- src/sumfem/electrochemistry/membrane.py | 6 +- .../test_electrochemistry_from_ufl.py | 85 +++++++++++++ tests/electrochemistry/test_electrolyte.py | 24 +++- 6 files changed, 209 insertions(+), 56 deletions(-) diff --git a/examples/electrochemistry/vanadium_half_cell_1d.py b/examples/electrochemistry/vanadium_half_cell_1d.py index 18960b9..53ce2e5 100644 --- a/examples/electrochemistry/vanadium_half_cell_1d.py +++ b/examples/electrochemistry/vanadium_half_cell_1d.py @@ -140,17 +140,17 @@ # The ions of the electrolyte are given with their charges and diffusivities, # in their units, which sumfem converts to SI, see {ref}`units-checking`. The # kinetics take the positive part of the vanadium concentrations, as they may -# be depleted. The concentrations -# of the reservoir name the ions that are transported, each with a field and -# the equations of {eq}`half-cell`. Sulfate has none: electroneutrality, -# the last equation of {eq}`half-cell`, with the charges of the ions, +# be depleted. The concentrations of the reservoir name the ions that are +# transported, each with a field and the equations of {eq}`half-cell`. +# Sulfate has none: electroneutrality, the last equation of {eq}`half-cell`, +# {eq}`electroneutrality` with the charges of the ions, # # $$ # \begin{aligned} # 2 c_{\mathrm{V^{2+}}} + 3 c_{\mathrm{V^{3+}}} + c_{\mathrm{H^+}} # - c_{\mathrm{HSO_4^-}} - 2 c_{\mathrm{SO_4^{2-}}} &= 0, # \end{aligned} -# $$ (electroneutrality) +# $$ (vanadium-electroneutrality) # # is solved for its concentration, which is then no unknown, but an # expression of the others: diff --git a/src/sumfem/electrochemistry/charge.py b/src/sumfem/electrochemistry/charge.py index a19dad6..fccd567 100644 --- a/src/sumfem/electrochemistry/charge.py +++ b/src/sumfem/electrochemistry/charge.py @@ -15,24 +15,15 @@ class ElectrolyteConduction(Residual): r"""Charge conservation in the electrolyte, :math:`-(\mathbf{i}_l, \nabla v)`. - The current density of the electrolyte, - :meth:`~sumfem.electrochemistry.electrolyte.Electrolyte.current_density`, - 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`. + The weak form of :math:`\nabla \cdot \mathbf{i}_l = 0`, with the current + density :math:`\mathbf{i}_l` of :eq:`electrolyte-current`, the charge the + ions carry, computed by the electrolyte's + :meth:`~sumfem.electrochemistry.electrolyte.Electrolyte.current_density` + with the diffusivities corrected for the region, :math:`t` the region's + :attr:`~sumfem.electrochemistry.porous.Region.transport`, as in + :eq:`nernst-planck-flux`. Sources, e.g. of an + :class:`~sumfem.electrochemistry.ElectrodeReaction`, are added to it as + models of their own. Args: name: Name of the model. diff --git a/src/sumfem/electrochemistry/electrolyte.py b/src/sumfem/electrochemistry/electrolyte.py index 32bd24b..662f14e 100644 --- a/src/sumfem/electrochemistry/electrolyte.py +++ b/src/sumfem/electrochemistry/electrolyte.py @@ -61,29 +61,76 @@ def field(self) -> Field: class Electrolyte: - """A dilute electrolyte solution: its ions, temperature and constants. - - The electrolyte creates its parameters, such as the temperature - ``".T"``, named ``"."``: the constants, the - temperature and the properties of the solution, see its properties, and - the diffusivity of each ion, see :meth:`diffusivity`. Models that share - an electrolyte share these parameters. + r"""A dilute electrolyte solution: its ions, temperature and constants. + + An electrolyte is a solution of ions in a solvent, e.g. vanadium ions, + protons and sulfate in water. In a dilute solution the ions do not + interact with each other: each ion :math:`k`, of charge number + :math:`z_k`, concentration :math:`c_k` and diffusivity :math:`D_k`, + moves by diffusion, by migration in the electric field and with the + solvent, with the Nernst-Planck flux :math:`\mathbf{N}_k` of + :eq:`nernst-planck-flux`. The solution is electroneutral: + + .. math:: + :label: electroneutrality + + \begin{aligned} + \sum_k z_k c_k &= 0. + \end{aligned} + + The state of the solution is the concentrations of its ions and the + temperature :math:`T`. Its properties are the diffusivities of the ions + and the density and viscosity of the solution. + + The electrolyte gives the quantities of dilute-solution theory that the + models of a solution share. These are the concentration of each ion, + :meth:`concentration`, and the charge the ions carry, the current + density :meth:`current_density`: + + .. math:: + :label: electrolyte-current + + \begin{aligned} + \mathbf{i}_l &= F \sum_k z_k \mathbf{N}_k, \\ + \mathbf{i}_l &= -\kappa \nabla \phi + - F \sum_k z_k t D_k \nabla c_k, \\ + \kappa &= \frac{F^2}{R T} \sum_k z_k^2 t D_k c_k, + \end{aligned} + + with the electrolyte potential :math:`\phi`, the ionic + :meth:`conductivity` :math:`\kappa`, the current of diffusion, + :meth:`diffusion_current`, and a correction :math:`t` of every + diffusivity, e.g. Bruggeman's in a porous electrode. Advection carries + no charge, by :eq:`electroneutrality`. The transport, conduction and + reaction models, e.g. :class:`~sumfem.electrochemistry.NernstPlanck`, + :class:`~sumfem.electrochemistry.ElectrolyteConduction` and + :class:`~sumfem.electrochemistry.ButlerVolmer`, build their terms from + one electrolyte, so they agree on these quantities. At most one ion may have no field. Its concentration follows from - electroneutrality, ``sum_j z_j c_j = 0``, see :attr:`balanced`. Nothing - keeps that concentration positive, so choose the most abundant ion of - the supporting electrolyte, e.g. sulfate in a vanadium electrolyte, - never an ion of a reaction, and check its minimum, e.g. with + :eq:`electroneutrality`, see :attr:`balanced`. Nothing keeps that + concentration positive, so choose the most abundant ion of the + supporting electrolyte, e.g. sulfate in a vanadium electrolyte, never an + ion of a reaction, and check its minimum, e.g. with :meth:`~sumfem.electrochemistry.RedoxHalfCell.minimum_balanced_concentration`. Without ions, the electrolyte has a uniform composition, and only provides the constants, the temperature and the properties of the solution. + Its constants, temperature and properties of the solution, and the + diffusivity of each ion, see :meth:`diffusivity`, are parameters named + ``"."``, e.g. the temperature ``".T"``. Models + that share an electrolyte share these parameters. The temperature, + density and viscosity may instead be given as expressions, e.g. fields + of a heat equation or of a flow, which couples the models to them. + Args: name: Prefix of the parameter names. domain: The domain of the parameters. ions: The ions. - temperature: Default temperature in K. + temperature: The temperature: an expression in K, e.g. the field of + a heat equation, or the default value of the parameter the + electrolyte creates. density: Density of the solution: an expression in kg/m^3, e.g. a parameter shared with a flow model, or the default value of the parameter the electrolyte creates. ``None`` for no default. @@ -116,7 +163,7 @@ class Electrolyte: _R: Parameter """The gas constant.""" - _T: Parameter + _T: ufl.core.expr.Expr """The temperature.""" _floor: Parameter @@ -133,7 +180,7 @@ def __init__( name: str, domain: Domain, ions: typing.Sequence[Ion] = (), - temperature: float = 298.15, + temperature: ufl.core.expr.Expr | float = 298.15, density: ufl.core.expr.Expr | float | None = None, viscosity: ufl.core.expr.Expr | float | None = None, concentration_floor: float = 1e-14, @@ -172,8 +219,16 @@ def __init__( unit=syu.joule / (syu.mol * syu.kelvin), doc="Universal gas constant", ) - self._T = Parameter( - f"{name}.T", domain, default=temperature, unit=syu.kelvin, doc="Temperature" + self._T = ( + temperature + if isinstance(temperature, ufl.core.expr.Expr) + else Parameter( + f"{name}.T", + domain, + default=temperature, + unit=syu.kelvin, + doc="Temperature", + ) ) self._floor = Parameter( f"{name}.floor", @@ -251,8 +306,8 @@ def R(self) -> Parameter: return self._R @property - def T(self) -> Parameter: - """The temperature.""" + def T(self) -> ufl.core.expr.Expr: + """The temperature, the parameter ``".T"`` unless given.""" return self._T @property @@ -310,11 +365,11 @@ def diffusivity(self, name: str) -> Parameter: def conductivity( self, correction: ufl.core.expr.Expr | float = 1.0 ) -> ufl.core.expr.Expr: - """Ionic conductivity ``F^2 / (R T) sum_j z_j^2 D_j c_j``. + r"""The ionic conductivity :math:`\kappa` of :eq:`electrolyte-current`. Args: - correction: Factor applied to every diffusivity, such as a - Bruggeman correction in a porous electrode. + correction: The correction :math:`t` of every diffusivity, e.g. + Bruggeman's in a porous electrode. """ total = sum( s.charge**2 * self.diffusivity(s.name) * self.concentration(s.name) @@ -325,7 +380,9 @@ def conductivity( def diffusion_current( self, correction: ufl.core.expr.Expr | float = 1.0 ) -> ufl.core.expr.Expr: - """Current carried by diffusion, ``F sum_j z_j D_j grad(c_j)``. + r"""The current of diffusion, :math:`F \sum_k z_k t D_k \nabla c_k`. + + It is the current of :eq:`electrolyte-current` without migration. Args: correction: See :meth:`conductivity`. @@ -341,16 +398,14 @@ def current_density( potential: ufl.core.expr.Expr, correction: ufl.core.expr.Expr | float = 1.0, ) -> ufl.core.expr.Expr: - """Current density ``-kappa grad(phi) - F sum_j z_j D_j grad(c_j)``. + r"""The current density :math:`\mathbf{i}_l` of :eq:`electrolyte-current`. - It is the charge the ions carry, ``F sum_j z_j N_j``, with the - Nernst-Planck fluxes ``N_j`` of every ion, also the one given by - electroneutrality: the conductivity and the migration use the same - concentrations, and advection carries no charge, as - ``sum_j z_j c_j = 0``. + It is the charge the ions carry, with the Nernst-Planck fluxes of + every ion, also the one given by :eq:`electroneutrality`, as the + conductivity and the migration use the same concentrations. Args: - potential: The electrolyte potential. + potential: The electrolyte potential :math:`\phi`. correction: See :meth:`conductivity`. """ return -self.conductivity(correction) * ufl.grad( diff --git a/src/sumfem/electrochemistry/membrane.py b/src/sumfem/electrochemistry/membrane.py index 29d5257..1e81d27 100644 --- a/src/sumfem/electrochemistry/membrane.py +++ b/src/sumfem/electrochemistry/membrane.py @@ -216,7 +216,9 @@ class ProtonExchangeMembrane(ResidualSum): 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. + temperature: The temperature :math:`T`: an expression in K, e.g. + the field of a heat equation, or the default value of a + parameter. """ _potential: Field @@ -250,7 +252,7 @@ def __init__( electrokinetic_permeability: PROPERTY = None, hydraulic_permeability: PROPERTY = None, viscosity: PROPERTY = None, - temperature: float = 298.15, + temperature: ufl.core.expr.Expr | float = 298.15, ): domain = potential.domain self._potential = potential diff --git a/tests/electrochemistry/test_electrochemistry_from_ufl.py b/tests/electrochemistry/test_electrochemistry_from_ufl.py index 86f7d6d..dd801c3 100644 --- a/tests/electrochemistry/test_electrochemistry_from_ufl.py +++ b/tests/electrochemistry/test_electrochemistry_from_ufl.py @@ -22,6 +22,7 @@ PistonVelocityButlerVolmer, PorousElectrode, ) +from sumfem.models import Diffusion from sumfem.symbolic import Domain, Field, Residual, UFLResidual CELLS = ["interval", "triangle", "quadrilateral", "tetrahedron", "hexahedron"] @@ -144,6 +145,90 @@ def flux(c, z, D, advected=False): ) +@pytest.mark.parametrize("cell", CELLS) +def test_conduction_and_transport_at_the_temperature_of_a_heat_equation( + cell, box, compare +): + """Conduction and Nernst-Planck transport of V2+, V3+ and H+, with SO4 2- + by electroneutrality, at a temperature that is the unknown of a heat + equation, so the electrolyte's equations depend on 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, V, V) + phi, c2, c3, h, T = (ufl.Coefficient(V) for _ in range(5)) + v_phi, w2, w3, wh, wT = ufl.TestFunctions(W) + F, R, k, D2, D3, DH, DS = (ufl.Constant(mesh) for _ in range(7)) + sulfate = (2 * c2 + 3 * c3 + h) / 2 + conductivity = ( + F**2 / (R * T) * (4 * D2 * c2 + 9 * D3 * c3 + DH * h + 4 * DS * sulfate) + ) + diffusion = F * ( + 2 * D2 * ufl.grad(c2) + + 3 * D3 * ufl.grad(c3) + + DH * ufl.grad(h) + - 2 * DS * ufl.grad(sulfate) + ) + current = -conductivity * ufl.grad(phi) - diffusion + + def flux(c, z, D): + return -D * ufl.grad(c) - D * z * F * c / (R * T) * ufl.grad(phi) + + # Charge conservation, div(i) = 0, in the equation of the potential + charge = ufl.inner(-current, ufl.grad(v_phi)) * ufl.dx + # The balance of each transported ion, div(N) = 0 + vanadium2 = -ufl.inner(flux(c2, 2, D2), ufl.grad(w2)) * ufl.dx + vanadium3 = -ufl.inner(flux(c3, 3, D3), ufl.grad(w3)) * ufl.dx + protons = -ufl.inner(flux(h, 1, DH), ufl.grad(wh)) * ufl.dx + # The heat equation, -div(k grad(T)) = 0 + heat = k * ufl.inner(ufl.grad(T), ufl.grad(wT)) * ufl.dx + residual = charge + vanadium2 + vanadium3 + protons + heat + unknowns = dict(phi=phi, c2=c2, c3=c3, h=h, T=T) + constants_ = dict(F=F, R=R, k=k, D2=D2, D3=D3, DH=DH, DS=DS) + plain = UFLResidual( + "plain", + residual, + list(unknowns.values()), + names={DOMAIN: mesh, **unknowns, **constants_}, + ) + + domain = Domain(DOMAIN, cell) + potential, *fields, temperature = ( + Field(name, domain, ("Lagrange", 1)) for name in unknowns + ) + ions = [ + Ion("V2", 2, 2.4e-10, fields[0]), + Ion("V3", 3, 2.4e-10, fields[1]), + Ion("H", 1, 9.3e-9, fields[2]), + Ion("SO4", -2, 1.065e-9), + ] + electrolyte = Electrolyte("e", domain, ions, temperature=temperature) + assert electrolyte.T is temperature + model = ElectrolyteConduction("charge", potential, electrolyte) + for ion in ions[:3]: + model += NernstPlanck(ion.name, ion, electrolyte, potential) + model += Diffusion("heat", temperature, 0.6) + + values = { + "F": electrolyte.F.default_si, + "R": electrolyte.R.default_si, + "k": 0.6, + "D2": 2.4e-10, + "D3": 2.4e-10, + "DH": 9.3e-9, + "DS": 1.065e-9, + } + compare( + plain, + plain.unknowns, + model, + [potential, *fields, temperature], + box(cell), + [(constants(values), {})], + ) + + @pytest.mark.parametrize("mass_transfer", [False, True]) @pytest.mark.parametrize("cell", CELLS) def test_kinetics_of_a_redox_couple(cell, mass_transfer, box, compare): diff --git a/tests/electrochemistry/test_electrolyte.py b/tests/electrochemistry/test_electrolyte.py index 7e48904..37ce528 100644 --- a/tests/electrochemistry/test_electrolyte.py +++ b/tests/electrochemistry/test_electrolyte.py @@ -16,8 +16,8 @@ Region, TransientNernstPlanck, ) -from sumfem.models import Accumulation -from sumfem.symbolic import Domain, Field, build_system +from sumfem.models import Accumulation, Diffusion +from sumfem.symbolic import Domain, Field, Parameter, build_system CONCENTRATION = syu.mol / syu.m**3 @@ -228,3 +228,23 @@ def test_the_balanced_ion_needs_a_charge(): c = Field("c", domain, ("Lagrange", 1)) with pytest.raises(ValueError, match="has no charge"): Electrolyte("e", domain, [Ion("A", 1, 1e-9, c), Ion("N", 0, 1e-9)]) + + +def test_the_temperature_can_be_the_field_of_a_heat_equation(domain): + """A temperature given as a field replaces the parameter ``"e.T"``, and + the electrolyte's models keep consistent units with it.""" + T = Field("T", domain, ("Lagrange", 1), unit=syu.kelvin, scale=300.0) + ions = vanadium(domain) + electrolyte = Electrolyte("e", domain, ions, temperature=T) + assert electrolyte.T is T + _, phi_l = potentials(domain) + model = ElectrolyteConduction("charge", phi_l, electrolyte) + for ion in ions[:3]: + model += NernstPlanck(ion.name, ion, electrolyte, phi_l) + conductivity = Parameter( + "heat.k", domain, default=0.6, unit=syu.watt / (syu.m * syu.kelvin) + ) + model += Diffusion("heat", T, conductivity) + system = build_system(model, check_units=True) + assert "e.T" not in system.parameters + assert T in system.unknowns From fb8d4b3250e28d16fb351d0a4983b21ee18cc565 Mon Sep 17 00:00:00 2001 From: Jorgen Schartum Dokken Date: Thu, 12 Nov 2026 10:46:46 +0000 Subject: [PATCH 3/3] Add first full cell model. --- docs/api/electrochemistry.md | 16 + docs/api/symbolic.md | 2 + src/sumfem/electrochemistry/__init__.py | 13 +- src/sumfem/electrochemistry/electrolyte.py | 49 +- src/sumfem/electrochemistry/full_cell.py | 571 +++++++++++++++++++++ src/sumfem/electrochemistry/half_cell.py | 5 + src/sumfem/electrochemistry/membrane.py | 224 ++++++-- src/sumfem/models/reactions.py | 15 +- src/sumfem/symbolic/__init__.py | 2 + src/sumfem/symbolic/components.py | 36 ++ src/sumfem/symbolic/domains.py | 28 + tests/electrochemistry/test_full_cell.py | 297 +++++++++++ tests/symbolic/test_components.py | 18 + tests/symbolic/test_domains.py | 16 +- 14 files changed, 1213 insertions(+), 79 deletions(-) create mode 100644 src/sumfem/electrochemistry/full_cell.py create mode 100644 tests/electrochemistry/test_full_cell.py diff --git a/docs/api/electrochemistry.md b/docs/api/electrochemistry.md index 77ec833..0b2afce 100644 --- a/docs/api/electrochemistry.md +++ b/docs/api/electrochemistry.md @@ -64,6 +64,9 @@ .. autoclass:: SchloeglFlow :members: + +.. autoclass:: MembraneProtonFlux + :members: ``` ## Half-cells @@ -75,3 +78,16 @@ .. autoclass:: RedoxHalfCell :members: ``` + +## Full cells + +```{eval-rst} +.. autoclass:: FullCell + :members: + +.. autoclass:: HalfCellSpec + :members: + +.. autoclass:: MembraneSpec + :members: +``` diff --git a/docs/api/symbolic.md b/docs/api/symbolic.md index 3de36ed..53a02af 100644 --- a/docs/api/symbolic.md +++ b/docs/api/symbolic.md @@ -56,6 +56,8 @@ .. autofunction:: input_or_parameter +.. autofunction:: expression_or_parameter + .. autoclass:: Residual :members: diff --git a/src/sumfem/electrochemistry/__init__.py b/src/sumfem/electrochemistry/__init__.py index 79fd35a..1d985bf 100644 --- a/src/sumfem/electrochemistry/__init__.py +++ b/src/sumfem/electrochemistry/__init__.py @@ -11,9 +11,13 @@ on :class:`ElectrodeReaction`, with the Butler-Volmer equation :func:`~sumfem.models.butler_volmer`. - Membranes: :class:`ProtonExchangeMembrane`, with the flow of water of - :class:`SchloeglFlow`. + :class:`SchloeglFlow`, whose current the protons of an electrode carry + through its face, :class:`MembraneProtonFlux`. - Half-cells of a redox flow battery: :class:`RedoxHalfCell`, negative or positive, of any :class:`RedoxCouple` and electrolyte. +- Full cells: :class:`FullCell`, two half-cells, each a + :class:`HalfCellSpec`, and a membrane between them, a + :class:`MembraneSpec`. The package imports from the core of :mod:`sumfem`, and never the other way around. @@ -22,6 +26,7 @@ from .charge import ElectrolyteConduction from .couple import RedoxCouple from .electrolyte import Electrolyte, Ion +from .full_cell import FullCell, HalfCellSpec, MembraneSpec from .half_cell import RedoxHalfCell from .kinetics import ( ButlerVolmer, @@ -30,7 +35,7 @@ MassTransferButlerVolmer, PistonVelocityButlerVolmer, ) -from .membrane import ProtonExchangeMembrane, SchloeglFlow +from .membrane import MembraneProtonFlux, ProtonExchangeMembrane, SchloeglFlow from .porous import PorousElectrode, Region from .transport import NernstPlanck, TransientNernstPlanck @@ -44,8 +49,12 @@ "TransientNernstPlanck", "ProtonExchangeMembrane", "SchloeglFlow", + "MembraneProtonFlux", "RedoxCouple", "RedoxHalfCell", + "FullCell", + "HalfCellSpec", + "MembraneSpec", "ButlerVolmer", "ElectrodeReaction", "ExchangeCurrentButlerVolmer", diff --git a/src/sumfem/electrochemistry/electrolyte.py b/src/sumfem/electrochemistry/electrolyte.py index 662f14e..7c89534 100644 --- a/src/sumfem/electrochemistry/electrolyte.py +++ b/src/sumfem/electrochemistry/electrolyte.py @@ -10,7 +10,14 @@ import sympy.physics.units as syu import ufl -from sumfem.symbolic import FARADAY, GAS_CONSTANT, Domain, Field, Parameter +from sumfem.symbolic import ( + FARADAY, + GAS_CONSTANT, + Domain, + Field, + Parameter, + expression_or_parameter, +) from sumfem.symbolic.si import to_si __all__ = ["Electrolyte", "Ion"] @@ -219,16 +226,8 @@ def __init__( unit=syu.joule / (syu.mol * syu.kelvin), doc="Universal gas constant", ) - self._T = ( - temperature - if isinstance(temperature, ufl.core.expr.Expr) - else Parameter( - f"{name}.T", - domain, - default=temperature, - unit=syu.kelvin, - doc="Temperature", - ) + self._T = expression_or_parameter( + f"{name}.T", temperature, domain, unit=syu.kelvin, doc="Temperature" ) self._floor = Parameter( f"{name}.floor", @@ -237,27 +236,15 @@ def __init__( unit=CONCENTRATION, doc="Regularization of positive concentrations", ) - self._density = ( - density - if isinstance(density, ufl.core.expr.Expr) - else Parameter( - f"{name}.density", - domain, - default=density, - unit=syu.kg / syu.m**3, - doc="Density", - ) + self._density = expression_or_parameter( + f"{name}.density", density, domain, unit=syu.kg / syu.m**3, doc="Density" ) - self._viscosity = ( - viscosity - if isinstance(viscosity, ufl.core.expr.Expr) - else Parameter( - f"{name}.viscosity", - domain, - default=viscosity, - unit=syu.pascal * syu.second, - doc="Dynamic viscosity", - ) + self._viscosity = expression_or_parameter( + f"{name}.viscosity", + viscosity, + domain, + unit=syu.pascal * syu.second, + doc="Dynamic viscosity", ) self._ions = types.MappingProxyType({ion.name: ion for ion in ions}) self._diffusivities = { diff --git a/src/sumfem/electrochemistry/full_cell.py b/src/sumfem/electrochemistry/full_cell.py new file mode 100644 index 0000000..8851a6b --- /dev/null +++ b/src/sumfem/electrochemistry/full_cell.py @@ -0,0 +1,571 @@ +"""The full cell of a redox flow battery: two half-cells and a membrane. + +:class:`FullCell` builds a cell of two :class:`RedoxHalfCell` and a +:class:`ProtonExchangeMembrane` between them, with a current collector on +the outer side of each electrode, on one mesh, each unknown on the part of +the mesh where it lives. Each side is given by a :class:`HalfCellSpec`. +""" + +from __future__ import annotations + +import dataclasses +import types +import typing + +import numpy as np +import sympy +import sympy.physics.units as syu +import ufl + +from sumfem.discrete import Discretization, MeshData +from sumfem.models import BoundaryFlux, Conduction +from sumfem.symbolic import ( + DirichletBC, + Domain, + Field, + Parameter, + Residual, + ResidualSum, + Subdomain, + expression_or_parameter, +) + +from .couple import RedoxCouple +from .electrolyte import CONCENTRATION, Ion +from .half_cell import RedoxHalfCell +from .membrane import PROPERTY, ProtonExchangeMembrane +from .porous import PorousElectrode + +__all__ = ["FullCell", "HalfCellSpec", "MembraneSpec"] + +SIDES = ("negative", "positive") +"""The sides of the cell, in the order of the cell's layers.""" + +WATER = "H2O" +"""Name of the water of the electrolytes, shared by both sides.""" + +PROTONS = "H" +"""Name of the protons, which carry the current through the membrane.""" + + +@dataclasses.dataclass(frozen=True) +class HalfCellSpec: + """One side of a full cell: its electrolyte, its couple and its electrode. + + Attributes: + ions: The ions of the electrolyte, see :class:`RedoxHalfCell`. + Water, ``"H2O"``, if transported, must be an ion of both sides, + since the water moves through the membrane. + reservoir: The initial concentration of each transported ion, see + :class:`RedoxHalfCell`. The protons, ``"H"``, must be one of + them. + couple: The redox couple of the side. + porosity: The porosity of the electrode. + conductivity: The conductivity of the electrode's fibres, in S/m or + with units. + specific_area: The specific area of the electrode, in 1/m or with + units. + fiber_distance: The mean distance between the fibres, in m or with + units, which limits the kinetics by the piston velocities, see + :class:`~sumfem.electrochemistry.PistonVelocityButlerVolmer`. + ``None`` for kinetics without mass transfer. + """ + + ions: tuple[Ion, ...] + reservoir: typing.Mapping[str, float | sympy.Expr] + couple: RedoxCouple + porosity: float + conductivity: float | sympy.Expr + specific_area: float | sympy.Expr + fiber_distance: float | sympy.Expr | None = None + + +@dataclasses.dataclass(frozen=True) +class MembraneSpec: + """The properties of the membrane of a full cell. + + Each is an expression, used as given, or the default of the parameter + ``"."``, in SI units or with units, ``None`` for a + parameter without a default, see :class:`ProtonExchangeMembrane`. + + Attributes: + fixed_charge: The concentration of the fixed charges, in mol/m^3. + proton_diffusivity: The proton diffusivity in the membrane, in + m^2/s. + water_diffusivity: The water diffusivity in the membrane, in m^2/s. + """ + + fixed_charge: PROPERTY = None + proton_diffusivity: PROPERTY = None + water_diffusivity: PROPERTY = None + + +class FullCell: + r"""The full cell of a redox flow battery, without flow. + + The cell of `Shah et al., 2008 + `_, corrected, without + the flow of the electrolyte: from left to right, a current collector, the + negative porous electrode, a proton-exchange membrane, the positive + porous electrode and a current collector. Its mesh has these cell + markers: + + - ``"NegativeCollector"``, ``"NegativeElectrode"``, ``"Membrane"``, + ``"PositiveElectrode"`` and ``"PositiveCollector"``, the layers; + - ``"NegativeSolid"`` and ``"PositiveSolid"``, a collector and its + electrode, where the potential of that side's solid lives; + - ``"Electrolyte"``, the electrodes and the membrane, where the + electrolyte potential and the water live, + + and these facet markers: + + - ``"NegativeTerminal"`` and ``"PositiveTerminal"``, the outer faces of + the collectors; + - ``"NegativeFace"`` and ``"PositiveFace"``, the faces of the membrane + to the electrodes. + + :meth:`mesh_data` adds the unions to mesh data with the layers. + + Each side is a :class:`RedoxHalfCell`, whose ions live in its electrode. + The electrolyte potential :math:`\phi_l` and the water are one field + each, over the electrodes and the membrane, so the current and the water + flux are continuous across the faces. The membrane is a + :class:`ProtonExchangeMembrane`, and its protons carry the current + through the faces, :eq:`membrane-proton-flux`, so that no other ion + crosses them. In particular the ion that follows from electroneutrality, + e.g. sulfate, stays on its side. This corrects Shah et al., who hold + the protons at the concentration of the fixed charges at the faces, + which lets that ion carry part of the current through them. + + The collectors conduct with the conductivity + :math:`\sigma_\text{coll}`, the parameter + ``".collector_conductivity"``. The cell is galvanostatic: the + current density :math:`j` of the cell, the parameter + ``".current"``, positive on discharge, enters the solid of the + positive side at its terminal on charge and leaves it on discharge, and + the solid of the negative side is grounded at its terminal, + :meth:`ground`. + + The model, :attr:`model`, is the sum of the models of the two + half-cells, named ``"negative_"`` and ``"positive_"``, + of the membrane, ``"_membrane_"``, of the conduction in + the collectors, ``"_negative_collector"`` and + ``"_positive_collector"``, of the current at the positive + terminal, ``"_terminal_current"``, and of the proton fluxes at the + faces, ``"_negative_protons"`` and ``"_positive_protons"``. + + Args: + domain: The domain of the cell, with the markers above. + negative: The negative side. + positive: The positive side. + membrane: The properties of the membrane. Parameters without + defaults by default. + collector_conductivity: The conductivity of the collectors: an + expression, used as given, e.g. one of the temperature, or the + default of the parameter ``".collector_conductivity"``, in + S/m or with units. + degree: Degree of the Lagrange elements of the unknowns. + name: Prefix of the cell's own parameters and models. + + Raises: + ValueError: If a side does not transport protons, or only one side + transports water. + """ + + _domain: Domain + """The domain of the cell.""" + + _name: str + """Prefix of the cell's parameters and models.""" + + _subdomains: types.MappingProxyType[str, Subdomain] + """The subdomains of the cell's unknowns, by cell marker.""" + + _phi_l: Field + """Potential of the electrolyte, over the electrodes and the membrane.""" + + _water: Field | None + """Concentration of the water, over the electrodes and the membrane.""" + + _half_cells: types.MappingProxyType[str, RedoxHalfCell] + """The half-cells, by side.""" + + _membrane: ProtonExchangeMembrane + """The membrane.""" + + _collector_conductivity: ufl.core.expr.Expr + """Conductivity of the current collectors.""" + + _current: Parameter + """Current density of the cell, positive on discharge.""" + + _model: ResidualSum + """The sum of all models of the cell.""" + + _functionals: types.MappingProxyType[str, ufl.Form] + """The integrals of the voltage, created once.""" + + def __init__( + self, + domain: Domain, + negative: HalfCellSpec, + positive: HalfCellSpec, + membrane: MembraneSpec = MembraneSpec(), + collector_conductivity: PROPERTY = 1000.0, + degree: int = 1, + name: str = "cell", + ): + specs = dict(zip(SIDES, (negative, positive), strict=True)) + _check_sides(specs) + self._domain = domain + self._name = name + tdim = domain.tdim + markers = [ + *(f"{side.capitalize()}{part}" for side in SIDES for part in _PARTS), + "Membrane", + "Electrolyte", + ] + self._subdomains = types.MappingProxyType( + {m: Subdomain(f"{name}_{_snake(m)}", domain, tdim, m) for m in markers} + ) + element = ("Lagrange", degree) + electrolyte = self._subdomains["Electrolyte"] + self._phi_l = Field("phi_l", electrolyte, element, unit=syu.volt) + self._water = ( + Field(f"c_{WATER}", electrolyte, element, unit=CONCENTRATION, scale=1e4) + if WATER in negative.reservoir + else None + ) + self._half_cells = types.MappingProxyType( + { + side: self._create_half_cell(side, spec, element) + for side, spec in specs.items() + } + ) + self._membrane = ProtonExchangeMembrane( + f"{name}_membrane", + self._phi_l, + ("Membrane",), + water=self._water, + fixed_charge=membrane.fixed_charge, + proton_diffusivity=membrane.proton_diffusivity, + water_diffusivity=membrane.water_diffusivity, + ) + self._collector_conductivity = expression_or_parameter( + f"{name}.collector_conductivity", + collector_conductivity, + domain, + unit=syu.siemens / syu.m, + doc="Conductivity of the current collectors", + ) + self._current = Parameter( + f"{name}.current", + domain, + default=0.0, + unit=syu.ampere / syu.m**2, + doc="Current density of the cell, positive on discharge", + ) + self._model = ResidualSum(self._create_models(), name=name) + self._functionals = self._create_functionals() + + def _create_half_cell( + self, side: str, spec: HalfCellSpec, element: tuple[str, int] + ) -> RedoxHalfCell: + """Create the half-cell of a side, with its unknowns on their subdomains. + + Args: + side: ``"negative"`` or ``"positive"``. + spec: The side. + element: The element of the unknowns. + + Returns: + The half-cell, whose collector and membrane are resolved. + """ + prefix = side.capitalize() + electrode_domain = self._subdomains[f"{prefix}Electrode"] + phi_s = Field( + f"{side}_phi_s", + self._subdomains[f"{prefix}Solid"], + element, + unit=syu.volt, + ) + fields = { + ion: Field( + f"{side}_c_{ion}", + electrode_domain, + element, + unit=CONCENTRATION, + scale=1e3, + ) + for ion in spec.reservoir + if ion != WATER + } + if self._water is not None: + fields[WATER] = self._water + electrode = PorousElectrode( + f"{side}_electrode", + electrode_domain, + (f"{prefix}Electrode",), + porosity=spec.porosity, + conductivity=spec.conductivity, + specific_area=spec.specific_area, + ) + return RedoxHalfCell( + electrode_domain, + spec.ions, + spec.reservoir, + spec.couple, + electrode, + fiber_distance=spec.fiber_distance, + collector=None, + membrane=None, + name=side, + phi_s=phi_s, + phi_l=self._phi_l, + fields=fields, + electrolyte_name=f"{side}_electrolyte", + ) + + def _create_models(self) -> list[Residual]: + """Create the models of the cell, listed in the class documentation. + + Returns: + The models, to be summed. + """ + name, domain = self._name, self._domain + models: list[Residual] = [ + *self._half_cells["negative"].model.components, + *self._half_cells["positive"].model.components, + *self._membrane.components, + ] + for side, cell in self._half_cells.items(): + models.append( + Conduction( + f"{name}_{side}_collector", + cell.phi_s, + self._collector_conductivity, + markers=(f"{side.capitalize()}Collector",), + ) + ) + # On discharge the positive couple is reduced: the electrons it takes + # come in through the terminal, so the current of the cell leaves the + # solid there, and enters it on charge + positive = self._half_cells["positive"] + inflow = ( + self._current if positive.couple.oxidation_on_discharge else -self._current + ) + models.append( + BoundaryFlux( + f"{name}_terminal_current", + positive.phi_s, + inflow, + ("PositiveTerminal",), + ) + ) + membrane = self._subdomains["Membrane"] + for side, cell in self._half_cells.items(): + prefix = side.capitalize() + models.append( + self._membrane.proton_flux( + f"{name}_{side}_protons", + cell.concentrations[PROTONS], + domain.interface( + self._subdomains[f"{prefix}Electrode"], + membrane, + f"{prefix}Face", + ), + charge=cell.electrolyte[PROTONS].charge, + ) + ) + return models + + def _create_functionals(self) -> types.MappingProxyType[str, ufl.Form]: + """Create the integrals of the voltage. + + Returns: + The areas and the potential integrals of the terminals, by name. + """ + functionals: dict[str, ufl.Form] = {} + for side, cell in self._half_cells.items(): + ds = cell.phi_s.domain.ds(f"{side.capitalize()}Terminal") + functionals[f"{side}_area"] = 1.0 * ds + functionals[f"{side}_potential"] = cell.phi_s * ds + return types.MappingProxyType(functionals) + + @staticmethod + def mesh_data(marked: MeshData) -> MeshData: + """Mesh data with the unions of the layers the cell's unknowns live on. + + Args: + marked: Mesh data whose cell markers are the five layers, and + whose facet markers are the terminals and the faces, see + :class:`FullCell`. + + Returns: + The mesh data with the cell markers ``"NegativeSolid"``, + ``"PositiveSolid"`` and ``"Electrolyte"`` added. + """ + tdim = marked.mesh.topology.dim + cells = dict(marked.markers[tdim]) + for side in ("Negative", "Positive"): + cells[f"{side}Solid"] = ( + *cells[f"{side}Collector"], + *cells[f"{side}Electrode"], + ) + cells["Electrolyte"] = ( + *cells["NegativeElectrode"], + *cells["Membrane"], + *cells["PositiveElectrode"], + ) + return MeshData(marked.mesh, marked.tags, {**marked.markers, tdim: cells}) + + @property + def model(self) -> Residual: + """The sum of all models of the cell.""" + return self._model + + @property + def half_cells(self) -> typing.Mapping[str, RedoxHalfCell]: + """The half-cells, ``"negative"`` and ``"positive"``.""" + return self._half_cells + + @property + def membrane(self) -> ProtonExchangeMembrane: + """The membrane.""" + return self._membrane + + @property + def phi_l(self) -> Field: + """Potential of the electrolyte, over the electrodes and the membrane.""" + return self._phi_l + + @property + def water(self) -> Field | None: + """Concentration of the water, if transported.""" + return self._water + + @property + def current(self) -> Parameter: + """Current density of the cell, positive on discharge.""" + return self._current + + def subdomain(self, marker: str) -> Subdomain: + """The subdomain of a cell marker of the cell, e.g. ``"Membrane"``.""" + return self._subdomains[marker] + + def ground(self) -> list[DirichletBC]: + """The negative terminal at zero potential, which fixes the potentials.""" + return [ + DirichletBC(self._half_cells["negative"].phi_s, "NegativeTerminal", 0.0) + ] + + def initialize(self, discretization: Discretization) -> None: + """Set a state at rest, as Shah et al. start. + + The concentrations start at the reservoir values of the sides, the + defaults of their parameters ``".c0_"``, the electrolyte + potential at zero, and the solid of each side at the equilibrium + potential of its couple at those concentrations, by Nernst's + equation, so that no current flows. + + Raises: + ValueError: If the sides start with different water + concentrations, which one water field cannot hold. + """ + values: dict[str | ufl.Coefficient, float] = {self._phi_l: 0.0} + for cell in self._half_cells.values(): + for ion, c in cell.concentrations.items(): + values[c] = _initial(cell, ion) + values[cell.phi_s] = self._equilibrium_potential(cell, discretization) + if self._water is not None: + start = {_initial(c, WATER) for c in self._half_cells.values()} + if len(start) != 1: + raise ValueError( + "Both sides must start with the same water concentration, got " + f"{sorted(start)}." + ) + discretization.set(values) + + @staticmethod + def _equilibrium_potential( + cell: RedoxHalfCell, discretization: Discretization + ) -> float: + """The equilibrium potential of a half-cell's couple at its reservoir. + + Args: + cell: The half-cell. + discretization: The discretization, for the parameter values. + + Returns: + Nernst's potential ``E0 + R T / (n F) ln(c_ox / c_red)`` in V. + """ + value = discretization.value + couple, e = cell.couple, cell.electrolyte + # A half-cell's electrolyte has the temperature as a parameter + assert isinstance(e.T, Parameter) + E0 = value(cell.formal_potential) + thermal = value(e.R) * value(e.T) / (couple.electrons * value(e.F)) + ratio = _initial(cell, couple.oxidized) / _initial(cell, couple.reduced) + return float(E0 + thermal * np.log(ratio)) + + def voltage(self, discretization: Discretization) -> float: + """The cell voltage in V: the mean potential of the positive terminal + minus that of the negative terminal.""" + assemble = discretization.assemble_scalar + f = self._functionals + means = { + side: assemble(f[f"{side}_potential"]) / assemble(f[f"{side}_area"]) + for side in SIDES + } + return means["positive"] - means["negative"] + + +def _initial(cell: RedoxHalfCell, ion: str) -> float: + """The initial concentration of an ion of a half-cell, in mol/m^3. + + Args: + cell: The half-cell. + ion: The name of a transported ion. + + Returns: + The default of the reservoir concentration of the ion. + """ + default = cell.inlet[ion].default_si + assert default is not None + return float(default) + + +_PARTS = ("Collector", "Electrode", "Solid") +"""The cell markers of each side, after the side's name.""" + + +def _snake(marker: str) -> str: + """A cell marker in snake case, e.g. ``"negative_solid"``. + + Args: + marker: A cell marker in camel case, e.g. ``"NegativeSolid"``. + """ + return "".join(f"_{c.lower()}" if c.isupper() else c for c in marker).lstrip("_") + + +def _check_sides(specs: typing.Mapping[str, HalfCellSpec]) -> None: + """Check that the sides transport protons, and water on both or neither. + + Args: + specs: The sides, by name. + + Raises: + ValueError: If a side does not transport protons, or only one side + transports water. + """ + for side, spec in specs.items(): + if PROTONS not in spec.reservoir: + raise ValueError( + f"The {side} side must transport protons, {PROTONS!r}, which carry " + "the current through the membrane." + ) + if len({WATER in spec.reservoir for spec in specs.values()}) != 1: + raise ValueError( + f"Water, {WATER!r}, moves through the membrane, so both sides or " + "neither must transport it." + ) diff --git a/src/sumfem/electrochemistry/half_cell.py b/src/sumfem/electrochemistry/half_cell.py index 3350732..03df185 100644 --- a/src/sumfem/electrochemistry/half_cell.py +++ b/src/sumfem/electrochemistry/half_cell.py @@ -717,6 +717,11 @@ def current(self) -> Parameter: """Current density through the membrane, positive on discharge.""" return self._current + @property + def formal_potential(self) -> Parameter: + """Formal potential of the couple.""" + return self._formal_potential + @property def area_ratio(self) -> Parameter: """Ratio of the membrane area to the current collector area.""" diff --git a/src/sumfem/electrochemistry/membrane.py b/src/sumfem/electrochemistry/membrane.py index 1e81d27..105759e 100644 --- a/src/sumfem/electrochemistry/membrane.py +++ b/src/sumfem/electrochemistry/membrane.py @@ -10,19 +10,20 @@ import ufl from sumfem.models import Accumulation, Conduction, ConservativeAdvection, Diffusion +from sumfem.models.operators import restricted from sumfem.symbolic import ( DirichletBC, Domain, Field, - Parameter, Residual, ResidualSum, as_expr, + expression_or_parameter, ) from .electrolyte import CONCENTRATION, DIFFUSIVITY, Electrolyte -__all__ = ["ProtonExchangeMembrane", "SchloeglFlow"] +__all__ = ["MembraneProtonFlux", "ProtonExchangeMembrane", "SchloeglFlow"] INPUT: typing.TypeAlias = ufl.core.expr.Expr | npt.ArrayLike | None """A model input: an expression, a value, or ``None`` for a parameter.""" @@ -58,8 +59,14 @@ class SchloeglFlow(Residual): :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. + :class:`~sumfem.models.Darcy` flow. + + The equation is a bulk equation on the cells of the membrane, + ``markers``, and the velocity holds there only: :math:`\rho` is the + charge of the membrane's mobile ions, which the fixed charges set. With + one pressure over the membrane and the electrodes beside it, each with + a flow model of its own, the pressure is continuous across the faces, + and so, weakly, is the normal velocity. Args: name: Name of the model. @@ -183,6 +190,16 @@ class ProtonExchangeMembrane(ResidualSum): is continuous across its faces, and so is the flux of water with one water field. + Every equation is a bulk equation on the cells of the membrane, + ``markers``, whose thickness the mesh resolves. Its fields may extend + beyond them: in a :class:`~sumfem.electrochemistry.FullCell` the + potential and the water live on both electrodes and the membrane, and + each electrode's models add the equations there. The faces then need no + term of their own for these fields: their continuity, and that of their + fluxes, is natural. The protons are not a field of the membrane; the + faces pass them to the electrodes' own fields with + :meth:`proton_flux`, one condition per face. + The membrane is the sum of ``"_conduction"``, a :class:`~sumfem.models.Conduction`; with a water field, ``"_water_storage"``, an :class:`~sumfem.models.Accumulation`, and @@ -239,6 +256,9 @@ class ProtonExchangeMembrane(ResidualSum): _markers: tuple[str, ...] """Cell markers of the membrane.""" + _conduction: Conduction + """Ohmic conduction in the membrane.""" + def __init__( self, name: str, @@ -260,14 +280,14 @@ def __init__( self._markers = tuple(markers) self._constants = Electrolyte(name, domain, temperature=temperature) e = self._constants - self._fixed_charge = _property( + self._fixed_charge = expression_or_parameter( f"{name}.fixed_charge", fixed_charge, domain, unit=CONCENTRATION, doc="Concentration of the fixed charges", ) - proton = _property( + proton = expression_or_parameter( f"{name}.proton_diffusivity", proton_diffusivity, domain, @@ -275,9 +295,10 @@ def __init__( 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._conduction = Conduction( + f"{name}_conduction", potential, conductivity, self._markers + ) + models: list[Residual] = [self._conduction] self._flow = None if pressure is not None: self._flow = SchloeglFlow( @@ -285,21 +306,21 @@ def __init__( pressure, potential, charge_density=e.F * self._fixed_charge, - electrokinetic_permeability=_property( + electrokinetic_permeability=expression_or_parameter( f"{name}.electrokinetic_permeability", electrokinetic_permeability, domain, unit=PERMEABILITY, doc="Electrokinetic permeability", ), - hydraulic_permeability=_property( + hydraulic_permeability=expression_or_parameter( f"{name}.hydraulic_permeability", hydraulic_permeability, domain, unit=PERMEABILITY, doc="Hydraulic permeability", ), - viscosity=_property( + viscosity=expression_or_parameter( f"{name}.viscosity", viscosity, domain, @@ -310,7 +331,7 @@ def __init__( ) models.append(self._flow) if water is not None: - diffusivity = _property( + diffusivity = expression_or_parameter( f"{name}.water_diffusivity", water_diffusivity, domain, @@ -344,9 +365,50 @@ def fixed_charge(self) -> ufl.core.expr.Expr: """The concentration of the fixed charges.""" return self._fixed_charge + @property + def markers(self) -> tuple[str, ...]: + """Cell markers of the membrane.""" + return self._markers + + @property + def constants(self) -> Electrolyte: + """Faraday's constant, the gas constant and the temperature.""" + return self._constants + + def current_density(self) -> ufl.core.expr.Expr: + r"""The current density :math:`-\kappa_m \nabla \phi` of the membrane.""" + return self._conduction.flux() + + def proton_flux( + self, + name: str, + protons: Field, + measure: ufl.Measure, + charge: float = 1.0, + ) -> MembraneProtonFlux: + """The membrane current carried by the protons of an electrode. + + See :class:`MembraneProtonFlux`. + + Args: + name: Name of the model. + protons: The proton concentration of the electrode. + measure: The interface measure of the face, with the electrode on + its ``"+"`` side and the membrane on its ``"-"`` side, from + :meth:`~sumfem.symbolic.domains.Domain.interface`. + charge: The charge number of the protons. + + Returns: + The proton flux at the face. + """ + return MembraneProtonFlux(name, protons, self, measure, charge) + @property def velocity(self) -> ufl.core.expr.Expr: - """The velocity of the water, zero without a pressure.""" + """The velocity of the water in the membrane, zero without a pressure. + + It holds on the membrane's cells only, see :class:`SchloeglFlow`. + """ if self._flow is None: return ufl.zero(self._potential.domain.geometric_dimension) return self._flow.velocity @@ -370,28 +432,118 @@ def proton_conditions( 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. +class MembraneProtonFlux(Residual): + r"""The protons of an electrode carry the current of a membrane beside it. + + A cation-exchange membrane passes protons only, so at the face between + an electrode and the membrane the protons of the electrode carry the + whole current :math:`\mathbf{i}_m` of the membrane, + + .. math:: + :label: membrane-proton-flux + + \begin{aligned} + \mathbf{N}_\text{H} \cdot \mathbf{n} + &= \frac{\mathbf{i}_m \cdot \mathbf{n}}{z_\text{H} F}, + \end{aligned} + + with the normal :math:`\mathbf{n}` out of the electrode, and every other + ion has no flux. The residual is + :math:`(\mathbf{N}_\text{H} \cdot \mathbf{n}, w)` on the face, with the + test function :math:`w` of the protons. The current is that of the + membrane's side of the face, :meth:`ProtonExchangeMembrane.current_density`. + + With the ion of the electrolyte that follows from electroneutrality, e.g. + sulfate, the zero flux of that ion is no natural condition: it holds + when the protons carry the current, by this condition. The proton + concentration at the face is not held at the fixed charge, so the model + has no Donnan potential across the face. + + The condition is a surface term of the electrode's proton equation, on + one face of the membrane and from one side, the electrode's: each face + needs its own, with the protons of the electrode on that face. It is + not an equation of the membrane, which has no proton field. 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. + name: Name of the model. + protons: The proton concentration :math:`c_\text{H}` of the + electrode. + membrane: The membrane. + measure: The interface measure of the face, with the electrode on + its ``"+"`` side and the membrane on its ``"-"`` side, from + :meth:`~sumfem.symbolic.domains.Domain.interface`. + charge: The charge number :math:`z_\text{H}` of the protons. + + Raises: + ValueError: If the measure is not that of an interface with the + membrane on its ``"-"`` side and the protons' subdomain on its + ``"+"`` side. """ - if isinstance(value, ufl.core.expr.Expr): - return value - return Parameter(name, domain, default=value, unit=unit, doc=doc) + + _protons: Field + """The proton concentration of the electrode.""" + + _membrane: ProtonExchangeMembrane + """The membrane.""" + + _measure: ufl.Measure + """The interface measure of the face.""" + + _charge: float + """The charge number of the protons.""" + + def __init__( + self, + name: str, + protons: Field, + membrane: ProtonExchangeMembrane, + measure: ufl.Measure, + charge: float = 1.0, + ): + super().__init__(name) + self._protons = protons + self._membrane = membrane + self._measure = measure + self._charge = charge + self._check_sides() + + def _check_sides(self) -> None: + """Check that the face has the protons' electrode on its ``"+"`` side + and the membrane on its ``"-"`` side. + + Raises: + ValueError: If the measure is not that of such faces. + """ + domain = self._measure.ufl_domain() + interfaces = ( + domain.interfaces_of(self._measure) if isinstance(domain, Domain) else () + ) + if not interfaces: + raise ValueError( + f"{self.name!r} needs the measure of a face of the membrane, from " + "Domain.interface, to know its sides." + ) + electrode = self._protons.domain + for first, second, marker in interfaces: + if second.marker not in self._membrane.markers: + raise ValueError( + f"The face {marker!r} of {self.name!r} must have the membrane " + f"on its '-' side, got {second.marker!r}." + ) + if electrode is not first: + raise ValueError( + f"The face {marker!r} of {self.name!r} must have the electrode " + f"of {self._protons.name!r} on its '+' side, got " + f"{first.marker!r}." + ) + + def flux(self) -> ufl.core.expr.Expr: + r"""The proton flux :math:`\mathbf{N}_\text{H} \cdot \mathbf{n}` out of + the electrode, on its side of the face.""" + n = ufl.FacetNormal(self._measure.ufl_domain()) + current = restricted(self._membrane.current_density(), "-") + F = restricted(self._membrane.constants.F, "+") + return ufl.dot(current, n("+")) / (self._charge * F) + + def residual(self) -> ufl.Form: + return self.flux() * self._protons.v("+") * self._measure diff --git a/src/sumfem/models/reactions.py b/src/sumfem/models/reactions.py index d618db9..bbd5e56 100644 --- a/src/sumfem/models/reactions.py +++ b/src/sumfem/models/reactions.py @@ -269,16 +269,13 @@ def _sides(self) -> tuple[Subdomain, Subdomain] | None: if measure is None or measure.integral_type() != "interior_facet": return None domain = measure.ufl_domain() - ids = measure.subdomain_id() - if ( - not isinstance(domain, Domain) - or not isinstance(ids, tuple) - or len(ids) != 1 - ): + if not isinstance(domain, Domain): + return None + interfaces = domain.interfaces_of(measure) + if len(interfaces) != 1: return None - names = {i: name for name, i in domain.markers[domain.tdim - 1].items()} - sides = domain.interface_sides(names.get(ids[0], "")) - return None if sides is None else (sides[0], sides[1]) + ((first, second, _),) = interfaces + return first, second def _reaction_domain(self, fields: list[Field]) -> Domain: """The domain of lowest dimension among the participants'.""" diff --git a/src/sumfem/symbolic/__init__.py b/src/sumfem/symbolic/__init__.py index 2883bb5..fed517f 100644 --- a/src/sumfem/symbolic/__init__.py +++ b/src/sumfem/symbolic/__init__.py @@ -33,6 +33,7 @@ JacobianSum, Residual, ResidualSum, + expression_or_parameter, input_or_parameter, ) from .constants import BOLTZMANN, FARADAY, GAS_CONSTANT, boltzmann_constant @@ -77,6 +78,7 @@ "Clock", "Component", "Functional", + "expression_or_parameter", "input_or_parameter", "FunctionalSum", "Jacobian", diff --git a/src/sumfem/symbolic/components.py b/src/sumfem/symbolic/components.py index 92dba7a..c8f8e6b 100644 --- a/src/sumfem/symbolic/components.py +++ b/src/sumfem/symbolic/components.py @@ -27,6 +27,7 @@ import typing import numpy.typing as npt +import sympy import ufl from ufl.domain import extract_unique_domain @@ -53,6 +54,7 @@ "JacobianSum", "Residual", "ResidualSum", + "expression_or_parameter", "input_or_parameter", ] @@ -89,6 +91,40 @@ def input_or_parameter( return Parameter(name, domain, shape=shape, default=default, unit=unit, doc=doc) +def expression_or_parameter( + name: str, + value: ufl.core.expr.Expr | npt.ArrayLike | sympy.Expr | None, + domain: Domain, + *, + shape: tuple[int, ...] = (), + unit: Unit | None = None, + doc: str = "", +) -> ufl.core.expr.Expr: + """An expression, or else a new parameter with a value as its default. + + Unlike :func:`input_or_parameter`, which turns a number into a constant, + a number here becomes the default of a parameter named ``name``, so that + it has a unit and can be changed on a discretization. An expression, + e.g. a parameter shared with another model or a field, is used as given. + + Args: + name: The full name of the parameter, e.g. ``"."``. + value: An expression, used as given, or the default of the + parameter, in ``unit`` or a quantity with units; ``None`` for a + parameter without a default. + domain: The domain of the parameter. + shape: The shape 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, shape=shape, default=value, unit=unit, doc=doc) + + class Component: """What residuals and functionals share: a name, and named inputs. diff --git a/src/sumfem/symbolic/domains.py b/src/sumfem/symbolic/domains.py index 8dc1078..89bd178 100644 --- a/src/sumfem/symbolic/domains.py +++ b/src/sumfem/symbolic/domains.py @@ -417,6 +417,34 @@ def interface_sides(self, marker: str) -> tuple[Subdomain, Subdomain, str] | Non """ return self._interfaces.get(marker) + def interfaces_of( + self, measure: ufl.Measure + ) -> tuple[tuple[Subdomain, Subdomain, str], ...]: + """The interfaces an interior facet measure of this domain covers. + + See :meth:`interface` and :meth:`interface_sides`. + + Args: + measure: An interior facet measure. + + Returns: + The two subdomains and the facet marker of each interface the + measure covers, ``"+"`` side first. Empty if it covers any + facets that are not an interface, e.g. those of a plain + :meth:`measure`, or if it is not this domain's. + """ + ids = measure.subdomain_id() + if measure.ufl_domain() is not self or measure.integral_type() != ( + "interior_facet" + ): + return () + ids = ids if isinstance(ids, tuple) else (ids,) + names = {i: name for name, i in self._markers[self.tdim - 1].items()} + sides = [self.interface_sides(names.get(i, "")) for i in ids] + if any(side is None for side in sides): + return () + return tuple(typing.cast(list[tuple[Subdomain, Subdomain, str]], sides)) + def __str__(self) -> str: return f"<{type(self).__name__} {self._name!r} ({self.cell_name})>" diff --git a/tests/electrochemistry/test_full_cell.py b/tests/electrochemistry/test_full_cell.py new file mode 100644 index 0000000..5a7f05e --- /dev/null +++ b/tests/electrochemistry/test_full_cell.py @@ -0,0 +1,297 @@ +"""The full cell of Shah et al. (2008), corrected, in one dimension.""" + +from mpi4py import MPI + +import basix.ufl +import dolfinx.mesh +import numpy as np +import pytest +import sympy.physics.units as syu + +from sumfem.discrete import DomainMaps, MeshData, discretize +from sumfem.electrochemistry import ( + FullCell, + HalfCellSpec, + Ion, + MembraneSpec, + RedoxCouple, +) +from sumfem.solvers import BDF, Simulation +from sumfem.symbolic import Clock, Domain, build_system + +DIFFUSIVITY = syu.m**2 / syu.second +CONCENTRATION = syu.mol / syu.m**3 + +LAYERS = ( + ("NegativeCollector", 6e-3, 6), + ("NegativeElectrode", 4e-3, 40), + ("Membrane", 180e-6, 6), + ("PositiveElectrode", 4e-3, 40), + ("PositiveCollector", 6e-3, 6), +) +"""The layers of the cell of Shah et al., Table 3: name, thickness in m and +number of cells.""" + +ACID = { + "H": 1200.0 * CONCENTRATION, + "HSO4": 1200.0 * CONCENTRATION, + "H2O": 4.2e4 * CONCENTRATION, +} +"""The acid and the water of both electrolytes, Table 2, with the water +corrected from 4.2e3 to 4.2e4 mol/m^3.""" + + +def ions(reduced, oxidized, charges, diffusivity): + """The ions of a side: its couple, protons, bisulfate, water and sulfate.""" + return ( + Ion(reduced, charges[0], diffusivity * DIFFUSIVITY), + Ion(oxidized, charges[1], diffusivity * DIFFUSIVITY), + Ion("H", 1, 9.312e-9 * DIFFUSIVITY), + Ion("HSO4", -1, 1.33e-9 * DIFFUSIVITY), + Ion("H2O", 0, 2.3e-9 * DIFFUSIVITY), + Ion("SO4", -2, 1.065e-9 * DIFFUSIVITY), + ) + + +ELECTRODE = dict( + porosity=0.68, + conductivity=500.0, + specific_area=2e6, + fiber_distance=10 * syu.micrometer, +) +"""The carbon felt of both sides, Tables 3 and 5.""" + +NEGATIVE = HalfCellSpec( + ions=ions("V2", "V3", (2, 3), 2.4e-10), + reservoir={"V2": 27.0 * CONCENTRATION, "V3": 1053.0 * CONCENTRATION, **ACID}, + couple=RedoxCouple("V2", "V3", 0.68 * 1.75e-7, -0.255, True), + **ELECTRODE, +) +"""The negative side at 5% state of charge.""" + +POSITIVE = HalfCellSpec( + ions=ions("VO", "VO2", (2, 1), 3.9e-10), + reservoir={"VO": 1053.0 * CONCENTRATION, "VO2": 27.0 * CONCENTRATION, **ACID}, + couple=RedoxCouple( + "VO", + "VO2", + 0.68 * 3e-9, + 1.004, + False, + reactants={"VO": 1, "H2O": 1}, + products={"VO2": 1, "H": 2}, + ), + **ELECTRODE, +) +"""The positive side at 5% state of charge.""" + +MEMBRANE = MembraneSpec( + fixed_charge=1200.0, proton_diffusivity=1.4e-9, water_diffusivity=5.75e-10 +) +"""The membrane, Tables 4 and 5.""" + + +def cell_mesh(comm=MPI.COMM_WORLD) -> MeshData: + """The layers of the cell on an interval, with their markers.""" + edges = np.cumsum([0.0, *(width for _, width, _ in LAYERS)]) + if comm.rank == 0: + x = np.unique( + np.concatenate( + [ + np.linspace(a, b, n + 1) + for (a, b), (*_, n) in zip(zip(edges, edges[1:]), LAYERS) + ] + ) + ).reshape(-1, 1) + cells = np.stack([np.arange(len(x) - 1), np.arange(1, len(x))], axis=1) + else: + x, cells = np.empty((0, 1)), np.empty((0, 2), dtype=np.int64) + element = basix.ufl.element("Lagrange", "interval", 1, shape=(1,)) + mesh = dolfinx.mesh.create_mesh(comm, cells, element, x) + tolerance = 1e-9 + + def between(a, b): + return lambda x: (x[0] >= a - tolerance) & (x[0] <= b + tolerance) + + def at(a): + return lambda x: np.isclose(x[0], a, atol=tolerance) + + layers = {name: between(a, b) for (name, *_), a, b in zip(LAYERS, edges, edges[1:])} + facets = { + "NegativeTerminal": at(edges[0]), + "NegativeFace": at(edges[2]), + "PositiveFace": at(edges[3]), + "PositiveTerminal": at(edges[5]), + } + return FullCell.mesh_data(MeshData.from_locators(mesh, {1: layers, 0: facets})) + + +def full_cell(): + """The cell of Shah et al. on its domain.""" + domain = Domain("cell", "interval", length_scale=4e-3) + return FullCell(domain, NEGATIVE, POSITIVE, membrane=MEMBRANE), domain + + +def test_the_full_cell_has_consistent_units_and_shared_fields(): + cell, _ = full_cell() + system = build_system(cell.model, check_units=True) + names = sorted(f.name for f in system.unknowns) + assert names.count("phi_l") == 1 and names.count("c_H2O") == 1 + assert {"negative_phi_s", "positive_phi_s", "negative_c_H", "positive_c_H"} <= { + *names + } + assert "cell.collector_conductivity" in system.parameters + models = [m.name for m in cell.model.components] + assert {"cell_negative_protons", "cell_positive_protons"} <= {*models} + assert not any("membrane_current" in m for m in models) + + +@pytest.mark.parametrize( + "change, message", + [ + (lambda r: {k: v for k, v in r.items() if k != "H"}, "transport protons"), + (lambda r: {k: v for k, v in r.items() if k != "H2O"}, "both sides"), + ], +) +def test_the_sides_must_share_protons_and_water(change, message): + domain = Domain("cell", "interval") + ions = tuple( + i for i in POSITIVE.ions if i.name in (*change(POSITIVE.reservoir), "SO4") + ) + positive = HalfCellSpec( + ions=ions, + reservoir=change(POSITIVE.reservoir), + couple=POSITIVE.couple, + **ELECTRODE, + ) + with pytest.raises(ValueError, match=message): + FullCell(domain, NEGATIVE, positive, membrane=MEMBRANE) + + +def simulate(current, steps, dt=1.0): + """Run the cell at a current density from rest, with what the tests need.""" + cell, domain = full_cell() + marked = cell_mesh() + disc = discretize( + [cell.model, *cell.ground()], + DomainMaps({domain: marked}), + max_quadrature_degree=4, + ) + cell.initialize(disc) + clock = Clock(domain, dt=dt) + simulation = Simulation( + disc, + cell.model, + cell.ground(), + clock=clock, + scheme=BDF(1), + petsc_options={"snes_linesearch_type": "bt", "snes_max_it": 50}, + check_units=True, + factorize_equations=True, + time_scale=dt, + ) + disc.set({cell.current: current}) + return cell, disc, simulation + + +def inventories(cell, disc): + """The amounts per unit area of what the cell conserves, by name.""" + amounts = {} + for side, half in cell.half_cells.items(): + e = half.electrolyte + dx = half.concentrations["H"].domain.dx() + storage = half.electrode.porosity + for ion in e.ions: + if ion == "H2O": + continue + amounts[f"{side}_{ion}"] = disc.assemble_scalar( + storage * e.concentration(ion) * dx + ) + water = cell.water + electrolyte = cell.subdomain("Electrolyte") + amounts["water"] = sum( + disc.assemble_scalar(weight * water * electrolyte.dx(marker)) + for weight, marker in ( + (cell.half_cells["negative"].electrode.porosity, "NegativeElectrode"), + (1.0, "Membrane"), + (cell.half_cells["positive"].electrode.porosity, "PositiveElectrode"), + ) + ) + return amounts + + +def test_a_cell_at_rest_stays_at_rest(): + """Started at the equilibrium potentials of its couples, without current, + the cell does not change, and its voltage is the difference of those + potentials.""" + cell, disc, simulation = simulate(0.0, steps=2) + voltage = cell.voltage(disc) + before = inventories(cell, disc) + simulation.run(num_steps=2) + assert cell.voltage(disc) == pytest.approx(voltage, abs=1e-9) + after = inventories(cell, disc) + for name, amount in before.items(): + assert after[name] == pytest.approx(amount, rel=1e-10), name + RT_F = 8.314462618 * 298.15 / 96485.33212 + expected = 1.004 + 0.255 + 2 * RT_F * np.log(27.0 / 1053.0) + assert voltage == pytest.approx(expected, rel=1e-6) + + +def test_a_discharging_cell_conserves_its_species(): + """On discharge, V2+ is oxidized and VO2+ reduced, each by the charge + ``Q = j t`` per unit area; the protons carry the current through the + membrane, so sulfate and bisulfate stay on their sides, each side loses + ``Q / F`` protons, and the water gains what VO2+ loses. + + In one dimension the balances hold to the tolerance of Newton's method: + the electrolyte potential's test functions of an electrode and its face + sum to one in the electrode, so the current of the membrane cell at the + face is the electrode's total reaction current, which the protons carry. + """ + current, steps, dt = 100.0, 5, 1.0 + cell, disc, simulation = simulate(current, steps, dt) + before = inventories(cell, disc) + simulation.run(num_steps=steps) + after = inventories(cell, disc) + change = {name: after[name] - before[name] for name in before} + moles = current * steps * dt / 96485.33212 + + for side, reduced, oxidized, sign in ( + ("negative", "V2", "V3", 1.0), + ("positive", "VO", "VO2", -1.0), + ): + # Vanadium stays on its side, and the couple turns over Q / F + assert change[f"{side}_{reduced}"] == pytest.approx(-sign * moles, rel=1e-8) + assert change[f"{side}_{oxidized}"] == pytest.approx(sign * moles, rel=1e-8) + # Neither anion crosses the membrane, and each side loses Q / F protons + assert change[f"{side}_HSO4"] == pytest.approx(0.0, abs=1e-8 * moles) + assert change[f"{side}_SO4"] == pytest.approx(0.0, abs=1e-8 * moles) + assert change[f"{side}_H"] == pytest.approx(-moles, rel=1e-8) + # One water per VO2+ reduced + assert change["water"] == pytest.approx(moles, rel=1e-8) + assert cell.voltage(disc) < 1.004 + 0.255 + + +def test_the_protons_cross_a_face_from_their_own_side(): + """The proton flux is a condition of an electrode's protons on its face, + with the membrane on the other side; a face seen from the wrong side, or + a measure that knows no sides, is rejected.""" + cell, domain = full_cell() + negative = cell.half_cells["negative"].concentrations["H"] + electrode = cell.subdomain("NegativeElectrode") + membrane = cell.subdomain("Membrane") + proton_flux = cell.membrane.proton_flux + with pytest.raises(ValueError, match="membrane on its '-' side"): + proton_flux( + "wrong", negative, domain.interface(membrane, electrode, "NegativeFace") + ) + with pytest.raises(ValueError, match="'\\+' side"): + proton_flux( + "wrong", + negative, + domain.interface( + cell.subdomain("PositiveElectrode"), membrane, "PositiveFace" + ), + ) + with pytest.raises(ValueError, match="Domain.interface"): + proton_flux("wrong", negative, domain.dS("NegativeFace")) diff --git a/tests/symbolic/test_components.py b/tests/symbolic/test_components.py index 41e1aae..8a2251f 100644 --- a/tests/symbolic/test_components.py +++ b/tests/symbolic/test_components.py @@ -18,6 +18,7 @@ ResidualSum, SystemBuildError, build_system, + expression_or_parameter, ) @@ -467,3 +468,20 @@ def residual(self): with pytest.raises(ComponentError, match="ZeroBaseForm"): Vanishing("vanishing").fields + + +def test_an_expression_is_used_and_a_value_becomes_a_parameter(): + """``expression_or_parameter`` keeps an expression, and makes a number or + a quantity the default of a parameter with a unit.""" + domain = Domain("cell", "interval") + shared = Parameter("other.k", domain, default=2.0, unit=syu.siemens / syu.m) + assert expression_or_parameter("m.k", shared, domain) is shared + unit = syu.siemens / syu.m + number = expression_or_parameter("m.k", 1000.0, domain, unit=unit, doc="k") + assert isinstance(number, Parameter) and number.name == "m.k" + assert number.default_si == pytest.approx(1000.0) + quantity = expression_or_parameter( + "m.k", 10 * syu.siemens / syu.cm, domain, unit=unit + ) + assert quantity.default_si == pytest.approx(1000.0) + assert expression_or_parameter("m.k", None, domain, unit=unit).default_si is None diff --git a/tests/symbolic/test_domains.py b/tests/symbolic/test_domains.py index 1d27e6f..92ce3e7 100644 --- a/tests/symbolic/test_domains.py +++ b/tests/symbolic/test_domains.py @@ -1,7 +1,7 @@ import pytest import ufl -from sumfem.symbolic import ComponentError, Domain +from sumfem.symbolic import ComponentError, Domain, Subdomain from sumfem.symbolic.domains import MEASURE_ENTITY_TYPE @@ -139,3 +139,17 @@ def test_invalid_markers_and_integral_types(): domain.marker_ids("A", "edge") with pytest.raises(TypeError): domain.dx(3) + + +def test_an_interface_measure_knows_its_sides(): + """An interface measure gives its subdomains, ``"+"`` side first, and its + marker; other facet measures give none.""" + square = Domain("square", "triangle") + left = Subdomain("left", square, 2, "Left") + right = Subdomain("right", square, 2, "Right") + measure = square.interface(left, right, "Interface") + assert square.interfaces_of(measure) == ((left, right, "Interface"),) + assert square.interfaces_of(square.dS("Interface")) == () + assert square.interfaces_of(square.dS()) == () + assert square.interfaces_of(square.ds("Interface")) == () + assert left.interfaces_of(measure) == ()