Temperature as a known or unknown expression/coefficient/field

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:

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.

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
-    ``"<name>.T"``, named ``"<name>.<parameter>"``: 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
+    ``"<name>.<parameter>"``, e.g. the temperature ``"<name>.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 ``"<name>.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(

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

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):

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