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(