diff --git a/examples/electrochemistry/vanadium_half_cell_1d.py b/examples/electrochemistry/vanadium_half_cell_1d.py index c20f353..34f5513 100644 --- a/examples/electrochemistry/vanadium_half_cell_1d.py +++ b/examples/electrochemistry/vanadium_half_cell_1d.py @@ -367,7 +367,10 @@ def value_of(name: str) -> float: # Charging reduces V³⁺ to V²⁺, so the state of charge grows, and discharging # turns it back. At rest only the exchange with the reservoir at the membrane # changes the electrode, slowly. The voltage is lower on charge than on -# discharge, by the overpotentials that drive the reaction either way: +# discharge, by the overpotentials that drive the reaction either way. +# Sulfate, which follows from electroneutrality, stays positive, see +# {py:meth}`minimum_balanced_concentration +# `: # + steps = {mode: np.diff(charge[modes == mode]) for mode in ("charge", "discharge")} @@ -378,6 +381,7 @@ def value_of(name: str) -> float: assert np.abs(np.diff(rest)).max() < np.abs(steps["charge"]).min() assert voltage[modes == "charge"].mean() < voltage[modes == "discharge"].mean() assert np.isfinite(voltage).all() +assert cell.minimum_balanced_concentration(disc) > 0 # - # The voltage and the state of charge through the schedule: diff --git a/examples/electrochemistry/vanadium_half_cell_2d.py b/examples/electrochemistry/vanadium_half_cell_2d.py index 93d7be8..1c7609e 100644 --- a/examples/electrochemistry/vanadium_half_cell_2d.py +++ b/examples/electrochemistry/vanadium_half_cell_2d.py @@ -380,7 +380,10 @@ # Charging raises the state of charge and discharging lowers it. The flow # now renews the electrolyte, so at rest it carries the charged electrolyte # out and brings the reservoir's in: the state of charge falls back towards -# that of the reservoir: +# that of the reservoir. Sulfate, which follows from electroneutrality, stays +# positive, see +# {py:meth}`minimum_balanced_concentration +# `: # + steps = {mode: np.diff(charge[modes == mode]) for mode in ("charge", "discharge")} @@ -392,6 +395,7 @@ assert np.all(np.diff(rest) < 0) assert voltage[modes == "charge"].mean() < voltage[modes == "discharge"].mean() assert np.isfinite(voltage).all() +assert cell.minimum_balanced_concentration(disc) > 0 # - # The voltage and the state of charge through the schedule: diff --git a/examples/electrochemistry/vanadium_half_cell_3d.py b/examples/electrochemistry/vanadium_half_cell_3d.py index 0a3b99a..e088b40 100644 --- a/examples/electrochemistry/vanadium_half_cell_3d.py +++ b/examples/electrochemistry/vanadium_half_cell_3d.py @@ -450,7 +450,10 @@ def flow_through_electrode(sizes: tuple[float, float]) -> None: # - # Charging raises the state of charge and discharging lowers it, and the -# voltage is lower on charge than on discharge: +# voltage is lower on charge than on discharge. Sulfate, which follows from +# electroneutrality, stays positive, see +# {py:meth}`minimum_balanced_concentration +# `: # + steps = {mode: np.diff(charge[modes == mode]) for mode in ("charge", "discharge")} @@ -459,6 +462,7 @@ def flow_through_electrode(sizes: tuple[float, float]) -> None: assert np.all(steps["charge"] > 0) and np.all(steps["discharge"] < 0) assert voltage[modes == "charge"].mean() < voltage[modes == "discharge"].mean() assert np.isfinite(voltage).all() +assert cell.minimum_balanced_concentration(disc) > 0 # - # The voltage and the state of charge through the schedule: diff --git a/src/sumfem/electrochemistry/electrolyte.py b/src/sumfem/electrochemistry/electrolyte.py index f06b9d0..ff3929a 100644 --- a/src/sumfem/electrochemistry/electrolyte.py +++ b/src/sumfem/electrochemistry/electrolyte.py @@ -74,9 +74,14 @@ class Electrolyte: an electrolyte share these parameters. At most one ion may have no field. Its concentration follows from - electroneutrality, ``sum_j z_j c_j = 0``. Without ions, the electrolyte - has a uniform composition, and only provides the constants, the - temperature and the properties of the solution. + 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 + :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. Args: name: Prefix of the parameter names. @@ -231,6 +236,14 @@ def ions(self) -> typing.Mapping[str, Ion]: """The ions by name, in the given order, as given.""" return self._ions + @property + def balanced(self) -> str | None: + """Name of the ion that follows from electroneutrality, if any.""" + for ion in self._ions.values(): + if ion.concentration is None: + return ion.name + return None + @property def F(self) -> Parameter: """Faraday's constant.""" diff --git a/src/sumfem/electrochemistry/half_cell.py b/src/sumfem/electrochemistry/half_cell.py index 504e632..235cf6d 100644 --- a/src/sumfem/electrochemistry/half_cell.py +++ b/src/sumfem/electrochemistry/half_cell.py @@ -13,6 +13,10 @@ import types import typing +from mpi4py import MPI + +import dolfinx.fem +import numpy as np import sympy import sympy.physics.units as syu import ufl @@ -790,6 +794,34 @@ def state_of_charge(self, discretization: Discretization) -> float: oxidized = assemble(self._functionals["amount_oxidized"]) return reduced / max(reduced + oxidized, 1e-30) + def minimum_balanced_concentration(self, discretization: Discretization) -> float: + """The smallest concentration of the ion given by electroneutrality. + + Nothing keeps that concentration positive. A negative minimum means + the electrolyte cannot stay neutral with the computed + concentrations, e.g. after a time step that is too long, or because + the balanced ion is not abundant enough, see + :class:`~sumfem.electrochemistry.Electrolyte`. + + Args: + discretization: The discretization of the cell. + + Returns: + The minimum over the nodes of the concentrations' space, over all + processes. + """ + e = self._electrolyte + 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)] + balanced.interpolate( + dolfinx.fem.Expression(expression, space.element.interpolation_points) + ) + owned = space.dofmap.index_map.size_local * space.dofmap.index_map_bs + local = balanced.x.array[:owned].min(initial=np.inf) + return float(space.mesh.comm.allreduce(local, op=MPI.MIN)) + def local_state_of_charge(self) -> ufl.core.expr.Expr: """The share of the reduced ion in the couple, ``c_red / (c_red + c_ox)``.""" c_red = self._concentrations[self._couple.reduced] diff --git a/tests/electrochemistry/half_cell/test_models.py b/tests/electrochemistry/half_cell/test_models.py index 6bd4020..548f6c4 100644 --- a/tests/electrochemistry/half_cell/test_models.py +++ b/tests/electrochemistry/half_cell/test_models.py @@ -401,3 +401,17 @@ def test_two_half_cells_share_the_electrolyte_potential_and_water(reference): assert {"negative_c_H", "positive_c_H", "c_VO", "c_V2"} <= set(names) assert "negative.rate_constant" in system.parameters assert "positive_electrolyte.D_H" in system.parameters + + +def test_the_balanced_ion_is_checked_for_its_minimum(reference): + """At rest the vanadium cell holds sulfate at the concentration that + neutralizes the others, (2 c_V2 + 3 c_V3 + c_H - c_HSO4) / 2.""" + marked = reference.interval(num_cells=8) + domain = Domain.from_mesh("cell", marked.mesh) + cell = reference.half_cell(domain) + inlet = cell.inlet_conditions("Membrane") + disc = discretize([cell.model, *inlet], DomainMaps({domain: marked})) + cell.initialize(disc) + assert cell.electrolyte.balanced == "SO4" + expected = (2 * 300.0 + 3 * 300.0 + 1000.0 - 1000.0) / 2 + assert cell.minimum_balanced_concentration(disc) == pytest.approx(expected)