From bc911ecc22b5047baea4e9f2c5512edbbbd284cb Mon Sep 17 00:00:00 2001 From: Jorgen Schartum Dokken Date: Thu, 12 Nov 2026 09:58:20 +0000 Subject: [PATCH] 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