From eecc650c322b7ac96185aad8a3ed0bae974b23f5 Mon Sep 17 00:00:00 2001 From: Jorgen Schartum Dokken Date: Thu, 12 Nov 2026 12:28:29 +0000 Subject: [PATCH 1/7] Start on v2 of a full cell model. Move units declarations to si units. --- _toc.yml | 1 + docs/api/electrochemistry.md | 3 + references.bib | 22 + src/sumfem/electrochemistry/__init__.py | 6 +- src/sumfem/electrochemistry/full_cell.py | 506 +++++++++++++++++++---- src/sumfem/symbolic/si.py | 28 +- tests/electrochemistry/test_full_cell.py | 264 +++++++++++- tests/electrochemistry/test_membrane.py | 22 +- 8 files changed, 738 insertions(+), 114 deletions(-) diff --git a/_toc.yml b/_toc.yml index 4b1e5d9..dea3e38 100644 --- a/_toc.yml +++ b/_toc.yml @@ -27,6 +27,7 @@ parts: - file: "examples/electrochemistry/vanadium_half_cell_1d.py" - file: "examples/electrochemistry/vanadium_half_cell_2d.py" - file: "examples/electrochemistry/vanadium_half_cell_3d.py" + - file: "examples/electrochemistry/vanadium_full_cell_2d.py" - caption: "Python API" chapters: - file: "docs/api/index" diff --git a/docs/api/electrochemistry.md b/docs/api/electrochemistry.md index 0b2afce..88e0fb2 100644 --- a/docs/api/electrochemistry.md +++ b/docs/api/electrochemistry.md @@ -90,4 +90,7 @@ .. autoclass:: MembraneSpec :members: + +.. autoclass:: FlowSpec + :members: ``` diff --git a/references.bib b/references.bib index 4d7a89e..0ea0f10 100644 --- a/references.bib +++ b/references.bib @@ -326,3 +326,25 @@ @article{Krishnamurthy2011 year = {2011}, doi = {10.1016/j.jpowsour.2011.08.024}, } + +@article{Shah2008, + author = {Shah, A. A. and Watt-Smith, M. J. and Walsh, F. C.}, + title = {A dynamic performance model for redox-flow batteries involving soluble species}, + journal = {Electrochimica Acta}, + volume = {53}, + number = {27}, + pages = {8087--8100}, + year = {2008}, + doi = {10.1016/j.electacta.2008.05.067}, +} + +@article{Knehr2012, + author = {Knehr, K. W. and Agar, Ertan and Dennison, C. R. and Kalidindi, A. R. and Kumbur, E. C.}, + title = {A Transient Vanadium Flow Battery Model Incorporating Vanadium Crossover and Water Transport through the Membrane}, + journal = {Journal of The Electrochemical Society}, + volume = {159}, + number = {9}, + pages = {A1446--A1459}, + year = {2012}, + doi = {10.1149/2.017209jes}, +} diff --git a/src/sumfem/electrochemistry/__init__.py b/src/sumfem/electrochemistry/__init__.py index 1d985bf..fc9712d 100644 --- a/src/sumfem/electrochemistry/__init__.py +++ b/src/sumfem/electrochemistry/__init__.py @@ -17,7 +17,8 @@ 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`. + :class:`MembraneSpec`, with the flow of the electrolytes from their + tanks, a :class:`FlowSpec`. The package imports from the core of :mod:`sumfem`, and never the other way around. @@ -26,7 +27,7 @@ from .charge import ElectrolyteConduction from .couple import RedoxCouple from .electrolyte import Electrolyte, Ion -from .full_cell import FullCell, HalfCellSpec, MembraneSpec +from .full_cell import FlowSpec, FullCell, HalfCellSpec, MembraneSpec from .half_cell import RedoxHalfCell from .kinetics import ( ButlerVolmer, @@ -55,6 +56,7 @@ "FullCell", "HalfCellSpec", "MembraneSpec", + "FlowSpec", "ButlerVolmer", "ElectrodeReaction", "ExchangeCurrentButlerVolmer", diff --git a/src/sumfem/electrochemistry/full_cell.py b/src/sumfem/electrochemistry/full_cell.py index 8851a6b..73071ea 100644 --- a/src/sumfem/electrochemistry/full_cell.py +++ b/src/sumfem/electrochemistry/full_cell.py @@ -3,7 +3,9 @@ :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`. +the mesh where it lives. Each side is given by a :class:`HalfCellSpec`, the +membrane by a :class:`MembraneSpec`, and the flow of the electrolytes from +their tanks, if any, by a :class:`FlowSpec`. """ from __future__ import annotations @@ -18,7 +20,8 @@ import ufl from sumfem.discrete import Discretization, MeshData -from sumfem.models import BoundaryFlux, Conduction +from sumfem.error_handling import ComponentError +from sumfem.models import BoundaryFlux, Conduction, Darcy from sumfem.symbolic import ( DirichletBC, Domain, @@ -28,6 +31,7 @@ ResidualSum, Subdomain, expression_or_parameter, + value_at_defaults, ) from .couple import RedoxCouple @@ -36,7 +40,7 @@ from .membrane import PROPERTY, ProtonExchangeMembrane from .porous import PorousElectrode -__all__ = ["FullCell", "HalfCellSpec", "MembraneSpec"] +__all__ = ["FlowSpec", "FullCell", "HalfCellSpec", "MembraneSpec"] SIDES = ("negative", "positive") """The sides of the cell, in the order of the cell's layers.""" @@ -47,6 +51,9 @@ PROTONS = "H" """Name of the protons, which carry the current through the membrane.""" +VISCOSITY = syu.pascal * syu.second +"""Unit of dynamic viscosities.""" + @dataclasses.dataclass(frozen=True) class HalfCellSpec: @@ -57,8 +64,8 @@ class HalfCellSpec: 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. + :class:`RedoxHalfCell`, and with flow, that of its tank. 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 @@ -68,7 +75,17 @@ class HalfCellSpec: 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. + ``None`` for kinetics without them, limited by the flow's mass + transfer with flow. + fiber_diameter: The diameter of the electrode's fibres, in m or with + units, for the permeability with flow, and the mass transfer + without a fibre distance. + kozeny_carman_constant: The Kozeny-Carman constant of the + permeability with flow, see + :func:`~sumfem.models.kozeny_carman`. + density: The density of the electrolyte, in kg/m^3 or with units. + viscosity: The dynamic viscosity of the electrolyte, in Pa s or with + units, which the flow and the mass transfer share. """ ions: tuple[Ion, ...] @@ -78,6 +95,10 @@ class HalfCellSpec: conductivity: float | sympy.Expr specific_area: float | sympy.Expr fiber_distance: float | sympy.Expr | None = None + fiber_diameter: float | sympy.Expr | None = None + kozeny_carman_constant: float | None = None + density: float | sympy.Expr = 1399.8 + viscosity: float | sympy.Expr = 6.963e-3 @dataclasses.dataclass(frozen=True) @@ -86,76 +107,108 @@ class MembraneSpec: 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`. + parameter without a default, see :class:`ProtonExchangeMembrane`. The + last three are those of the flow of water, with a :class:`FlowSpec`. 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. + electrokinetic_permeability: The electrokinetic permeability, in + m^2. + hydraulic_permeability: The hydraulic permeability, in m^2. + viscosity: The viscosity of the water in the membrane, in Pa s. """ fixed_charge: PROPERTY = None proton_diffusivity: PROPERTY = None water_diffusivity: PROPERTY = None + electrokinetic_permeability: PROPERTY = None + hydraulic_permeability: PROPERTY = None + viscosity: PROPERTY = None + + +@dataclasses.dataclass(frozen=True) +class FlowSpec: + r"""The flow of the electrolytes through the electrodes, from their tanks. + + One pressure :math:`p` lives on both electrodes and the membrane: + :class:`~sumfem.models.Darcy` flow in each electrode, with the + Kozeny-Carman permeability of its fibres and the viscosity of its + electrolyte, and :class:`~sumfem.electrochemistry.SchloeglFlow` in the + membrane, so the pressure is continuous across the faces. The + electrolyte of each side enters its electrode at the inlet with the + superficial velocity :math:`v_\text{in}` and the concentrations + :math:`c_k^\text{in}` of its tank, and leaves at the outlet, held at the + pressure :math:`p_\text{out}`, with the flow, + + .. math:: + :label: full-cell-flow + + \begin{aligned} + -\mathbf{u} \cdot \mathbf{n} &= v_\text{in} + && \text{at the inlet}, \\ + -\mathbf{N}_k \cdot \mathbf{n} &= v_\text{in} c_k^\text{in} + && \text{at the inlet}, \\ + p &= p_\text{out} + && \text{at the outlet}, \\ + (\mathbf{N}_k - \mathbf{u} c_k) \cdot \mathbf{n} &= 0 + && \text{at the outlet}, + \end{aligned} + + with the superficial velocity :math:`\mathbf{u}` of the electrode and + the Nernst-Planck flux :math:`\mathbf{N}_k` of each ion. The inlet holds + the flux of the ions, the flow's, rather than their concentration as + Shah et al. do, so that the cell takes in exactly what its tank gives. + + Attributes: + inlet_velocity: The superficial velocity at the inlets of both + electrodes: an expression, or the default of the parameter + ``".inlet_velocity"``, in m/s or with units. + outlet_pressure: The gauge pressure at the outlets, given like the + inlet velocity, in Pa. + tank_ratio: The volume of each tank over the volume of electrolyte in + the pores of its electrode, ``None`` for tanks so large that the + electrolyte at the inlets does not change. + """ + + inlet_velocity: PROPERTY = None + outlet_pressure: PROPERTY = 0.0 + tank_ratio: float | None = None class FullCell: - r"""The full cell of a redox flow battery, without flow. + r"""The full cell of a redox flow battery. 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"``. + `_, corrected: from + left to right, a current collector, the negative porous electrode, a + proton-exchange membrane, the positive porous electrode and a current + collector, each a layer of the mesh, see :meth:`mesh_data`. Each unknown + lives where it is defined: the potential of each side's solid on its + collector and electrode, the electrolyte potential, the water and, with + flow, the pressure on both electrodes and the membrane, and the other + ions of each side in its electrode. + + The cell is made of: + + - a :class:`RedoxHalfCell` for each side, whose collector and membrane + are layers of the mesh; + - the :class:`ProtonExchangeMembrane`, whose protons carry the current + through each face, :class:`MembraneProtonFlux`, so the ion that + follows from electroneutrality, e.g. sulfate, stays on its side. Shah + et al. hold the protons at the fixed charge at the faces instead, + which lets that ion carry part of the current through them; + - Ohmic conduction in the collectors, :class:`~sumfem.models.Conduction`; + - the current density of the cell at the positive terminal, with the + negative terminal grounded, see :meth:`conditions`; + - with flow, :class:`~sumfem.models.Darcy` flow in the electrodes from + their inlets to their outlets, :eq:`full-cell-flow`, and tanks, + :meth:`update_tanks`. Args: - domain: The domain of the cell, with the markers above. + domain: The domain of the cell, with the markers of :meth:`mesh_data`. negative: The negative side. positive: The positive side. membrane: The properties of the membrane. Parameters without @@ -164,6 +217,7 @@ class FullCell: expression, used as given, e.g. one of the temperature, or the default of the parameter ``".collector_conductivity"``, in S/m or with units. + flow: The flow of the electrolytes, ``None`` for none. degree: Degree of the Lagrange elements of the unknowns. name: Prefix of the cell's own parameters and models. @@ -187,6 +241,21 @@ class FullCell: _water: Field | None """Concentration of the water, over the electrodes and the membrane.""" + _flow: FlowSpec | None + """The flow of the electrolytes, if any.""" + + _pressure: Field | None + """Pressure over the electrodes and the membrane, with flow.""" + + _inlet_velocity: ufl.core.expr.Expr | None + """Superficial velocity at the inlets, with flow.""" + + _outlet_pressure: ufl.core.expr.Expr | None + """Gauge pressure at the outlets, with flow.""" + + _darcy: types.MappingProxyType[str, Darcy] + """The flow in each electrode, by side, empty without flow.""" + _half_cells: types.MappingProxyType[str, RedoxHalfCell] """The half-cells, by side.""" @@ -203,7 +272,7 @@ class FullCell: """The sum of all models of the cell.""" _functionals: types.MappingProxyType[str, ufl.Form] - """The integrals of the voltage, created once.""" + """The integrals of the voltage and the tanks, created once.""" def __init__( self, @@ -212,6 +281,7 @@ def __init__( positive: HalfCellSpec, membrane: MembraneSpec = MembraneSpec(), collector_conductivity: PROPERTY = 1000.0, + flow: FlowSpec | None = None, degree: int = 1, name: str = "cell", ): @@ -219,6 +289,7 @@ def __init__( _check_sides(specs) self._domain = domain self._name = name + self._flow = flow tdim = domain.tdim markers = [ *(f"{side.capitalize()}{part}" for side in SIDES for part in _PARTS), @@ -236,9 +307,25 @@ def __init__( if WATER in negative.reservoir else None ) + electrodes = { + side: self._create_electrode(side, spec) for side, spec in specs.items() + } + viscosities = { + side: expression_or_parameter( + f"{side}_electrolyte.viscosity", + spec.viscosity, + self._subdomains[f"{side.capitalize()}Electrode"], + unit=VISCOSITY, + doc="Dynamic viscosity", + ) + for side, spec in specs.items() + } + self._create_flow(electrodes, viscosities, element) self._half_cells = types.MappingProxyType( { - side: self._create_half_cell(side, spec, element) + side: self._create_half_cell( + side, spec, element, electrodes[side], viscosities[side] + ) for side, spec in specs.items() } ) @@ -247,9 +334,13 @@ def __init__( self._phi_l, ("Membrane",), water=self._water, + pressure=self._pressure, fixed_charge=membrane.fixed_charge, proton_diffusivity=membrane.proton_diffusivity, water_diffusivity=membrane.water_diffusivity, + electrokinetic_permeability=membrane.electrokinetic_permeability, + hydraulic_permeability=membrane.hydraulic_permeability, + viscosity=membrane.viscosity, ) self._collector_conductivity = expression_or_parameter( f"{name}.collector_conductivity", @@ -268,8 +359,104 @@ def __init__( self._model = ResidualSum(self._create_models(), name=name) self._functionals = self._create_functionals() + def _create_electrode(self, side: str, spec: HalfCellSpec) -> PorousElectrode: + """Create the porous electrode of a side, with its parameters. + + Args: + side: ``"negative"`` or ``"positive"``. + spec: The side. + + Returns: + The electrode, on the cells of the side's electrode. + """ + prefix = side.capitalize() + return PorousElectrode( + f"{side}_electrode", + self._subdomains[f"{prefix}Electrode"], + (f"{prefix}Electrode",), + porosity=spec.porosity, + conductivity=spec.conductivity, + specific_area=spec.specific_area, + fiber_diameter=spec.fiber_diameter, + kozeny_carman_constant=spec.kozeny_carman_constant, + ) + + def _create_flow( + self, + electrodes: typing.Mapping[str, PorousElectrode], + viscosities: typing.Mapping[str, ufl.core.expr.Expr], + element: tuple[str, int], + ) -> None: + """Create the pressure, its parameters and the flow in the electrodes. + + The pressure's scale is the pressure drop of a uniform flow at the + inlet velocity through the negative electrode over the domain's + length scale, if the defaults give it, else one pascal. + + Args: + electrodes: The electrodes, by side. + viscosities: The viscosities of the electrolytes, by side. + element: The element of the pressure. + """ + flow, name, domain = self._flow, self._name, self._domain + if flow is None: + self._pressure = None + self._inlet_velocity = self._outlet_pressure = None + self._darcy = types.MappingProxyType({}) + return + self._inlet_velocity = expression_or_parameter( + f"{name}.inlet_velocity", + flow.inlet_velocity, + domain, + unit=syu.m / syu.second, + doc="Superficial velocity at the inlets", + ) + self._outlet_pressure = expression_or_parameter( + f"{name}.outlet_pressure", + flow.outlet_pressure, + domain, + unit=syu.pascal, + doc="Gauge pressure at the outlets", + ) + negative = electrodes["negative"] + try: + scale = float( + value_at_defaults( + viscosities["negative"] + * self._inlet_velocity + * domain.length_scale + / negative.permeability() + ) + ) + except ComponentError: + scale = 1.0 + self._pressure = Field( + "p", + self._subdomains["Electrolyte"], + element, + unit=syu.pascal, + scale=scale, + ) + self._darcy = types.MappingProxyType( + { + side: Darcy( + f"{name}_{side}_darcy", + self._pressure, + electrode.permeability(), + viscosities[side], + markers=electrode.markers, + ) + for side, electrode in electrodes.items() + } + ) + def _create_half_cell( - self, side: str, spec: HalfCellSpec, element: tuple[str, int] + self, + side: str, + spec: HalfCellSpec, + element: tuple[str, int], + electrode: PorousElectrode, + viscosity: ufl.core.expr.Expr, ) -> RedoxHalfCell: """Create the half-cell of a side, with its unknowns on their subdomains. @@ -277,9 +464,12 @@ def _create_half_cell( side: ``"negative"`` or ``"positive"``. spec: The side. element: The element of the unknowns. + electrode: The electrode of the side. + viscosity: The viscosity of the side's electrolyte. Returns: - The half-cell, whose collector and membrane are resolved. + The half-cell, whose collector and membrane are resolved, with the + flow of its electrode and its outlet, if the cell has flow. """ prefix = side.capitalize() electrode_domain = self._subdomains[f"{prefix}Electrode"] @@ -302,23 +492,20 @@ def _create_half_cell( } 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, - ) + darcy = self._darcy.get(side) return RedoxHalfCell( electrode_domain, spec.ions, spec.reservoir, spec.couple, electrode, + velocity=None if darcy is None else darcy.superficial_velocity, fiber_distance=spec.fiber_distance, collector=None, membrane=None, + outlet=() if darcy is None else (f"{prefix}Outlet",), + density=spec.density, + viscosity=viscosity, name=side, phi_s=phi_s, phi_l=self._phi_l, @@ -377,29 +564,86 @@ def _create_models(self) -> list[Residual]: charge=cell.electrolyte[PROTONS].charge, ) ) + models.extend(self._create_inflow()) + return models + + def _create_inflow(self) -> list[Residual]: + """Create the flow in the electrodes and the inflow of the inlets. + + Returns: + The models of :eq:`full-cell-flow` but the outlets': the Darcy + flow and the inflow of each electrode, and the inflow of each of + its ions, empty without flow. + """ + if self._pressure is None: + return [] + assert self._inlet_velocity is not None + name, velocity = self._name, self._inlet_velocity + models: list[Residual] = [] + for side, cell in self._half_cells.items(): + inlet = (f"{side.capitalize()}Inlet",) + models.append(self._darcy[side]) + models.append( + BoundaryFlux(f"{name}_{side}_inflow", self._pressure, velocity, inlet) + ) + for ion, c in cell.concentrations.items(): + models.append( + BoundaryFlux( + f"{name}_{side}_inflow_{ion}", + c, + velocity * cell.inlet[ion], + inlet, + ) + ) return models def _create_functionals(self) -> types.MappingProxyType[str, ufl.Form]: - """Create the integrals of the voltage. + """Create the integrals of the voltage and of the tanks. Returns: - The areas and the potential integrals of the terminals, by name. + The areas and the potential integrals of the terminals, and with + flow, for each side, the volume of electrolyte in the pores of + its electrode, the area of its inlet and the outflow of each ion, + 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") + prefix = side.capitalize() + ds = cell.phi_s.domain.ds(f"{prefix}Terminal") functionals[f"{side}_area"] = 1.0 * ds functionals[f"{side}_potential"] = cell.phi_s * ds + if side not in self._darcy: + continue + electrode = self._subdomains[f"{prefix}Electrode"] + functionals[f"{side}_pore_volume"] = ( + cell.electrode.porosity * electrode.dx() + ) + functionals[f"{side}_inlet_area"] = 1.0 * electrode.ds(f"{prefix}Inlet") + u = self._darcy[side].superficial_velocity + for ion, c in cell.concentrations.items(): + n = ufl.FacetNormal(c.domain) + functionals[f"{side}_outflow_{ion}"] = ( + ufl.dot(u, n) * c * c.domain.ds(f"{prefix}Outlet") + ) 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. + The cell's mesh has the cell markers of its layers, + ``"NegativeCollector"``, ``"NegativeElectrode"``, ``"Membrane"``, + ``"PositiveElectrode"`` and ``"PositiveCollector"``, and the facet + markers ``"NegativeTerminal"`` and ``"PositiveTerminal"``, the outer + faces of the collectors, ``"NegativeFace"`` and ``"PositiveFace"``, + the faces of the membrane, and with flow, ``"NegativeInlet"``, + ``"NegativeOutlet"``, ``"PositiveInlet"`` and ``"PositiveOutlet"``, + the ends of the electrodes. This adds ``"NegativeSolid"`` and + ``"PositiveSolid"``, a collector and its electrode, and + ``"Electrolyte"``, the electrodes and the membrane. + Args: - marked: Mesh data whose cell markers are the five layers, and - whose facet markers are the terminals and the faces, see - :class:`FullCell`. + marked: Mesh data with the layers and the facet markers. Returns: The mesh data with the cell markers ``"NegativeSolid"``, @@ -421,7 +665,17 @@ def mesh_data(marked: MeshData) -> MeshData: @property def model(self) -> Residual: - """The sum of all models of the cell.""" + """The sum of all models of the cell. + + These are the models of the half-cells, ``"negative_"`` and + ``"positive_"``, of the membrane, ``"_membrane_"``, + the conduction in the collectors, ``"__collector"``, the + current at the positive terminal, ``"_terminal_current"``, the + proton fluxes at the faces, ``"__protons"``, and with + flow, the flow in each electrode, ``"__darcy"``, its + inflow, ``"__inflow"``, and that of each of its ions, + ``"__inflow_"``. + """ return self._model @property @@ -444,6 +698,16 @@ def water(self) -> Field | None: """Concentration of the water, if transported.""" return self._water + @property + def pressure(self) -> Field | None: + """Pressure over the electrodes and the membrane, with flow.""" + return self._pressure + + @property + def darcy(self) -> typing.Mapping[str, Darcy]: + """The flow in each electrode, by side, empty without flow.""" + return self._darcy + @property def current(self) -> Parameter: """Current density of the cell, positive on discharge.""" @@ -453,20 +717,35 @@ 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 [ + def conditions(self) -> list[DirichletBC]: + """The Dirichlet conditions of the cell. + + These are the negative terminal at zero potential, which fixes the + potentials, and with flow, the outlet pressure. + """ + conditions = [ DirichletBC(self._half_cells["negative"].phi_s, "NegativeTerminal", 0.0) ] + if self._pressure is not None: + assert self._outlet_pressure is not None + conditions.append( + DirichletBC( + self._pressure, + ("NegativeOutlet", "PositiveOutlet"), + self._outlet_pressure, + ) + ) + return conditions 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 at zero, the solid of each side at the equilibrium potential of its couple at those concentrations, by Nernst's - equation, so that no current flows. + equation, so that no current flows, and the pressure at that of the + outlets. Raises: ValueError: If the sides start with different water @@ -484,6 +763,9 @@ def initialize(self, discretization: Discretization) -> None: "Both sides must start with the same water concentration, got " f"{sorted(start)}." ) + if self._pressure is not None: + assert isinstance(self._outlet_pressure, Parameter) + values[self._pressure] = float(discretization.value(self._outlet_pressure)) discretization.set(values) @staticmethod @@ -519,6 +801,66 @@ def voltage(self, discretization: Discretization) -> float: } return means["positive"] - means["negative"] + def tank_volume(self, discretization: Discretization, side: str) -> float: + """The volume of a side's tank, per unit of the dimensions the mesh + lacks, e.g. per metre of depth for a 2D cell. + + Args: + discretization: The discretization of the cell. + side: ``"negative"`` or ``"positive"``. + + Raises: + ValueError: For a cell without flow, or without tanks. + """ + if self._flow is None or self._flow.tank_ratio is None: + raise ValueError(f"The cell {self._name!r} has no tanks.") + pores = discretization.assemble_scalar(self._functionals[f"{side}_pore_volume"]) + return self._flow.tank_ratio * pores + + def update_tanks(self, discretization: Discretization, dt: float) -> None: + r"""Take a step of the tanks after a time step of the cell. + + Each tank, of volume :math:`V`, mixes what the outlet of its + electrode returns, Shah et al.'s eq. 33 with the fluxes of the + outlet and the inlet, + + .. math:: + :label: full-cell-tank + + \begin{aligned} + V \frac{\mathrm{d} c_k^\text{in}}{\mathrm{d} t} + &= \int_\text{outlet} \mathbf{u} \cdot \mathbf{n}\, c_k + \,\mathrm{d}s + - v_\text{in} \int_\text{inlet} c_k^\text{in} \,\mathrm{d}s, + \end{aligned} + + see :eq:`full-cell-flow`. The step takes the outflow of the cell's + new state and the inflow of the inlet concentrations the cell took + its step with, so the cell and its tanks together conserve every + ion. The concentrations of the tanks are the reservoir parameters of + the half-cells, ``".c0_"``. A cell without tanks keeps + them. + + Args: + discretization: The discretization of the cell, after the step. + dt: The size of the step, in s. + """ + if self._flow is None or self._flow.tank_ratio is None: + return + assemble, value = discretization.assemble_scalar, discretization.value + assert isinstance(self._inlet_velocity, Parameter) + velocity = float(value(self._inlet_velocity)) + updates: dict[str | ufl.Coefficient, float] = {} + for side, cell in self._half_cells.items(): + volume = self.tank_volume(discretization, side) + area = assemble(self._functionals[f"{side}_inlet_area"]) + for ion in cell.concentrations: + tank = float(value(cell.inlet[ion])) + outflow = assemble(self._functionals[f"{side}_outflow_{ion}"]) + inflow = velocity * area * tank + updates[cell.inlet[ion]] = tank + dt * (outflow - inflow) / volume + discretization.set(updates) + def _initial(cell: RedoxHalfCell, ion: str) -> float: """The initial concentration of an ion of a half-cell, in mol/m^3. diff --git a/src/sumfem/symbolic/si.py b/src/sumfem/symbolic/si.py index d7da669..16d860f 100644 --- a/src/sumfem/symbolic/si.py +++ b/src/sumfem/symbolic/si.py @@ -18,7 +18,18 @@ from sympy.physics.units.systems.si import SI from ufl_units import get_factor -__all__ = ["Unit", "has_units", "molecule", "si_factor", "to_si"] +__all__ = [ + "CONCENTRATION", + "DENSITY", + "DIFFUSIVITY", + "PERMEABILITY", + "VISCOSITY", + "Unit", + "has_units", + "molecule", + "si_factor", + "to_si", +] Unit: typing.TypeAlias = sympy.Expr """A sympy unit expression, e.g. ``syu.meter / syu.second``.""" @@ -29,6 +40,21 @@ SI.set_quantity_dimension(molecule, amount_of_substance) SI.set_quantity_scale_factor(molecule, syu.mol / syu.avogadro_number) +CONCENTRATION = syu.mol / syu.m**3 +"""Unit of concentrations.""" + +DENSITY = syu.kg / syu.m**3 +"""Unit of densities.""" + +DIFFUSIVITY = syu.m**2 / syu.second +"""Unit of diffusivities.""" + +PERMEABILITY = syu.m**2 +"""Unit of permeabilities.""" + +VISCOSITY = syu.pascal * syu.second +"""Unit of dynamic viscosities.""" + def _as_unit(unit: Unit | None) -> sympy.Expr: return sympy.sympify(1) if unit is None else unit diff --git a/tests/electrochemistry/test_full_cell.py b/tests/electrochemistry/test_full_cell.py index 5a7f05e..a5207a4 100644 --- a/tests/electrochemistry/test_full_cell.py +++ b/tests/electrochemistry/test_full_cell.py @@ -1,4 +1,6 @@ -"""The full cell of Shah et al. (2008), corrected, in one dimension.""" +"""The full cell of Shah et al. (2008), corrected, in one and two dimensions.""" + +import dataclasses from mpi4py import MPI @@ -10,14 +12,23 @@ from sumfem.discrete import DomainMaps, MeshData, discretize from sumfem.electrochemistry import ( + DarcyFlow, + FlowSpec, FullCell, HalfCellSpec, Ion, MembraneSpec, RedoxCouple, ) +from sumfem.models import Accumulation from sumfem.solvers import BDF, Simulation -from sumfem.symbolic import Clock, Domain, build_system +from sumfem.symbolic import ( + Clock, + Domain, + Parameter, + build_system, + value_at_defaults, +) DIFFUSIVITY = syu.m**2 / syu.second CONCENTRATION = syu.mol / syu.m**3 @@ -41,7 +52,7 @@ corrected from 4.2e3 to 4.2e4 mol/m^3.""" -def ions(reduced, oxidized, charges, diffusivity): +def create_ions(reduced, oxidized, charges, diffusivity): """The ions of a side: its couple, protons, bisulfate, water and sulfate.""" return ( Ion(reduced, charges[0], diffusivity * DIFFUSIVITY), @@ -62,7 +73,7 @@ def ions(reduced, oxidized, charges, diffusivity): """The carbon felt of both sides, Tables 3 and 5.""" NEGATIVE = HalfCellSpec( - ions=ions("V2", "V3", (2, 3), 2.4e-10), + ions=create_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, @@ -70,7 +81,7 @@ def ions(reduced, oxidized, charges, diffusivity): """The negative side at 5% state of charge.""" POSITIVE = HalfCellSpec( - ions=ions("VO", "VO2", (2, 1), 3.9e-10), + ions=create_ions("VO", "VO2", (2, 1), 3.9e-10), reservoir={"VO": 1053.0 * CONCENTRATION, "VO2": 27.0 * CONCENTRATION, **ACID}, couple=RedoxCouple( "VO", @@ -91,7 +102,7 @@ def ions(reduced, oxidized, charges, diffusivity): """The membrane, Tables 4 and 5.""" -def cell_mesh(comm=MPI.COMM_WORLD) -> MeshData: +def create_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: @@ -126,14 +137,14 @@ def at(a): return FullCell.mesh_data(MeshData.from_locators(mesh, {1: layers, 0: facets})) -def full_cell(): +def create_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() + cell, _ = create_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 @@ -158,22 +169,22 @@ def test_the_sides_must_share_protons_and_water(change, message): 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): + positive = HalfCellSpec( + ions=ions, + reservoir=change(POSITIVE.reservoir), + couple=POSITIVE.couple, + **ELECTRODE, + ) 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() + cell, domain = create_full_cell() + marked = create_cell_mesh() disc = discretize( - [cell.model, *cell.ground()], + [cell.model, *cell.conditions()], DomainMaps({domain: marked}), max_quadrature_degree=4, ) @@ -182,7 +193,7 @@ def simulate(current, steps, dt=1.0): simulation = Simulation( disc, cell.model, - cell.ground(), + cell.conditions(), clock=clock, scheme=BDF(1), petsc_options={"snes_linesearch_type": "bt", "snes_max_it": 50}, @@ -276,7 +287,7 @@ 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() + cell, domain = create_full_cell() negative = cell.half_cells["negative"].concentrations["H"] electrode = cell.subdomain("NegativeElectrode") membrane = cell.subdomain("Membrane") @@ -295,3 +306,218 @@ def test_the_protons_cross_a_face_from_their_own_side(): ) with pytest.raises(ValueError, match="Domain.interface"): proton_flux("wrong", negative, domain.dS("NegativeFace")) + + +HEIGHT = 1e-2 +"""The height of the 2D cell, from the inlets to the outlets, in m.""" + +FLOW = FlowSpec(inlet_velocity=2.5e-3, tank_ratio=8.0) +"""Shah et al.'s 1 mL/s through an electrode of 4 mm by 10 cm, Table 2, +and tanks of about the paper's 250 mL each, Table 3.""" + +FLOWING = dict(fiber_diameter=10 * syu.micrometer, kozeny_carman_constant=5.55) +"""What the flow needs of the felt of both sides, Tables 3 and 5.""" + +MEMBRANE_FLOW = MembraneSpec( + fixed_charge=1200.0, + proton_diffusivity=1.4e-9, + water_diffusivity=5.75e-10, + electrokinetic_permeability=1.13e-19, + hydraulic_permeability=1.58e-18, + viscosity=1e-3, +) +"""The membrane with the flow of its water, Table 5.""" + + +def create_flowing_mesh( + comm=MPI.COMM_WORLD, columns=(2, 6, 2, 6, 2), rows=4 +) -> MeshData: + """The layers of the cell on a rectangle of the cell's height, with the + inlets of the electrodes at the bottom and their outlets at the top.""" + edges = np.cumsum([0.0, *(width for _, width, _ in LAYERS)]) + xs = np.unique( + np.concatenate( + [np.linspace(a, b, n + 1) for a, b, n in zip(edges, edges[1:], columns)] + ) + ) + mesh = dolfinx.mesh.create_rectangle( + comm, [[0.0, 0.0], [1.0, HEIGHT]], [len(xs) - 1, rows] + ) + # Move the vertices of the uniform grid onto the edges of the layers + index = np.rint(mesh.geometry.x[:, 0] * (len(xs) - 1)).astype(int) + mesh.geometry.x[:, 0] = xs[index] + tolerance = 1e-9 + + def between(a, b): + return lambda x: (x[0] >= a - tolerance) & (x[0] <= b + tolerance) + + def end(a, b, y): + return lambda x: between(a, b)(x) & np.isclose(x[1], y, atol=tolerance) + + layers = {name: between(a, b) for (name, *_), a, b in zip(LAYERS, edges, edges[1:])} + facets = { + "NegativeTerminal": lambda x: np.isclose(x[0], edges[0], atol=tolerance), + "NegativeFace": lambda x: np.isclose(x[0], edges[2], atol=tolerance), + "PositiveFace": lambda x: np.isclose(x[0], edges[3], atol=tolerance), + "PositiveTerminal": lambda x: np.isclose(x[0], edges[5], atol=tolerance), + "NegativeInlet": end(edges[1], edges[2], 0.0), + "NegativeOutlet": end(edges[1], edges[2], HEIGHT), + "PositiveInlet": end(edges[3], edges[4], 0.0), + "PositiveOutlet": end(edges[3], edges[4], HEIGHT), + } + marked = MeshData.from_locators(mesh, {2: layers, 1: facets}) + return FullCell.mesh_data(marked) + + +def create_flowing_cell(current, dt=1.0, flow=FLOW): + """The cell of Shah et al. in 2D, with flow and tanks, from rest.""" + domain = Domain("cell", "triangle", length_scale=HEIGHT) + sides = [dataclasses.replace(s, **FLOWING) for s in (NEGATIVE, POSITIVE)] + cell = FullCell(domain, *sides, membrane=MEMBRANE_FLOW, flow=flow) + disc = discretize( + [cell.model, *cell.conditions()], + DomainMaps({domain: create_flowing_mesh()}), + max_quadrature_degree=4, + ) + cell.initialize(disc) + simulation = Simulation( + disc, + cell.model, + cell.conditions(), + clock=Clock(domain, dt=dt), + 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, domain + + +def totals(cell, disc): + """The amounts per unit depth in the cell and its tanks, by name.""" + amounts = inventories(cell, disc) + water = 0.0 + for side, half in cell.half_cells.items(): + volume = cell.tank_volume(disc, side) + for ion in half.concentrations: + tank = volume * float(disc.value(half.inlet[ion])) + if ion == "H2O": + water += tank + else: + amounts[f"{side}_{ion}"] += tank + amounts["water"] += water + return amounts + + +def test_a_cell_at_rest_with_flow_stays_at_rest(): + """The tanks hold the electrolyte the cell starts with, so without current + the flow changes nothing, and the pressure falls through each electrode + by the drop ``mu v H / K`` of a uniform flow.""" + cell, disc, simulation, domain = create_flowing_cell(0.0) + before, voltage = totals(cell, disc), cell.voltage(disc) + simulation.run(num_steps=2, callback=lambda t, step: cell.update_tanks(disc, 1.0)) + after = totals(cell, disc) + for name, amount in before.items(): + assert after[name] == pytest.approx(amount, rel=1e-8), name + assert cell.voltage(disc) == pytest.approx(voltage, abs=1e-8) + + electrode = cell.half_cells["negative"].electrode + viscosity = cell.half_cells["negative"].electrolyte.viscosity + drop = value_at_defaults(viscosity * 2.5e-3 * HEIGHT / electrode.permeability()) + inlet = cell.subdomain("NegativeElectrode").ds("NegativeInlet") + mean = disc.assemble_scalar(cell.pressure * inlet) / disc.assemble_scalar( + 1.0 * inlet + ) + assert mean == pytest.approx(drop, rel=1e-6) + + +def test_a_discharging_cell_and_its_tanks_conserve_every_ion(): + """With flow and tanks, the cell and its tanks together keep the vanadium + and the bisulfate of each side, turn over ``Q / F`` of each couple, with + ``Q`` the charge through the terminals of the cell's height, and gain the + water that the positive couple makes. + + Each side loses the protons that carry the current through the membrane, + ``Q / F``, only to first order in the mesh size in two dimensions: the + current of the membrane's cells at a face differs from the current the + electrode's reactions send to it by the discretization error. It halves + as the mesh is refined. + """ + current, steps, dt = 100.0, 3, 1.0 + cell, disc, simulation, _ = create_flowing_cell(current, dt) + before = totals(cell, disc) + simulation.run( + num_steps=steps, callback=lambda t, step: cell.update_tanks(disc, dt) + ) + after = totals(cell, disc) + change = {name: after[name] - before[name] for name in before} + moles = current * HEIGHT * steps * dt / 96485.33212 + for side, reduced, oxidized, sign in ( + ("negative", "V2", "V3", 1.0), + ("positive", "VO", "VO2", -1.0), + ): + assert change[f"{side}_{reduced}"] == pytest.approx(-sign * moles, rel=1e-8) + assert change[f"{side}_{oxidized}"] == pytest.approx(sign * moles, rel=1e-8) + assert change[f"{side}_HSO4"] == pytest.approx(0.0, abs=1e-8 * moles) + assert change[f"{side}_H"] == pytest.approx(-moles, rel=2e-2) + assert change["water"] == pytest.approx(moles, rel=1e-8) + + +class StoredDarcyFlow(DarcyFlow): + """Darcy flow with a storage of the pressure in the electrodes: a flow + that changes in time, plugged into the cell in a few lines.""" + + @property + def models(self): + storage = Accumulation( + f"{self.setup.name}_storage", + self.pressure, + weight=Parameter( + "storage.storativity", + self.setup.domain, + default=1e-3, + unit=1 / syu.pascal, + ), + markers=("NegativeElectrode", "PositiveElectrode"), + ) + return [*super().models, storage] + + +def test_any_flow_plugs_into_the_cell(): + """A flow of the cell's own making, here one that changes in time, takes + the place of Darcy flow. Its storage keeps the pressure below the steady + drop of Darcy flow over the steps, so its velocity is not divergence + free; the cell and its tanks still conserve the vanadium and turn over + ``Q / F`` of each couple, as the species balances hold the flux of + whatever velocity the flow gives.""" + current, steps, dt = 100.0, 3, 1.0 + flow = dataclasses.replace(FLOW, model=StoredDarcyFlow) + cell, disc, simulation, _ = create_flowing_cell(current, dt, flow=flow) + assert isinstance(cell.flow, StoredDarcyFlow) + assert "cell_storage" in [m.name for m in cell.model.components] + before = totals(cell, disc) + simulation.run( + num_steps=steps, callback=lambda t, step: cell.update_tanks(disc, dt) + ) + after = totals(cell, disc) + moles = current * HEIGHT * steps * dt / 96485.33212 + for side, reduced, oxidized, sign in ( + ("negative", "V2", "V3", 1.0), + ("positive", "VO", "VO2", -1.0), + ): + change = after[f"{side}_{oxidized}"] - before[f"{side}_{oxidized}"] + assert change == pytest.approx(sign * moles, rel=1e-8) + vanadium = (after[f"{side}_{reduced}"] + after[f"{side}_{oxidized}"]) - ( + before[f"{side}_{reduced}"] + before[f"{side}_{oxidized}"] + ) + assert vanadium == pytest.approx(0.0, abs=1e-8 * moles) + electrode = cell.half_cells["negative"].electrode + viscosity = cell.half_cells["negative"].electrolyte.viscosity + drop = value_at_defaults(viscosity * 2.5e-3 * HEIGHT / electrode.permeability()) + inlet = cell.subdomain("NegativeElectrode").ds("NegativeInlet") + mean = disc.assemble_scalar(cell.pressure * inlet) / disc.assemble_scalar( + 1.0 * inlet + ) + assert 0 < mean < 0.99 * drop diff --git a/tests/electrochemistry/test_membrane.py b/tests/electrochemistry/test_membrane.py index 5a37956..45b53ee 100644 --- a/tests/electrochemistry/test_membrane.py +++ b/tests/electrochemistry/test_membrane.py @@ -44,7 +44,7 @@ T = 298.15 -def steady(model): +def create_steady_membrane(model): """The membrane without the storage of water.""" return ResidualSum( [m for m in model.components if not m.name.endswith("storage")], @@ -52,7 +52,7 @@ def steady(model): ) -def membrane(domain, water=True, pressure=True): +def create_membrane(domain, water=True, pressure=True): phi = Field("phi", domain, ("Lagrange", 1), unit=syu.volt) c_w = Field("c_w", domain, ("Lagrange", 1), unit=syu.mol / syu.m**3, scale=4e4) p = Field("p", domain, ("Lagrange", 1), unit=syu.pascal, scale=1e3) @@ -93,7 +93,7 @@ def test_the_membrane_against_plain_ufl(cell, box, compare): plain = UFLResidual("plain", residual, [phi, c_w, p], names=names) domain = Domain(DOMAIN, cell) - model, *fields = membrane(domain) + model, *fields = create_membrane(domain) values = { "plain.F": F, "plain.R": R, @@ -108,7 +108,7 @@ def test_the_membrane_against_plain_ufl(cell, box, compare): compare( plain, plain.unknowns, - steady(model), + create_steady_membrane(model), fields, box(cell), [(values, {})], @@ -117,7 +117,7 @@ def test_the_membrane_against_plain_ufl(cell, box, compare): def test_the_membrane_has_consistent_units(): domain = Domain("cell", "interval") - model, *_ = membrane(domain) + model, *_ = create_membrane(domain) build_system(model, check_units=True) build_system(model, substitutions=BDF(1).substitutions(model), check_units=True) names = [m.name for m in model.components] @@ -128,13 +128,13 @@ def test_the_membrane_has_consistent_units(): "membrane_water_diffusion", "membrane_water_advection", ] - still, *_ = membrane(domain, pressure=False) + still, *_ = create_membrane(domain, pressure=False) assert [m.name for m in still.components] == [ "membrane_conduction", "membrane_water_storage", "membrane_water_diffusion", ] - bare, *_ = membrane(domain, water=False, pressure=False) + bare, *_ = create_membrane(domain, water=False, pressure=False) assert [m.name for m in bare.components] == ["membrane_conduction"] @@ -158,7 +158,7 @@ def test_a_membrane_between_two_potentials_and_pressures(): boundary=True, ) domain = Domain.from_mesh("membrane", mesh, length_scale=L) - model, phi, c_w, p = membrane(domain) + model, phi, c_w, p = create_membrane(domain) dphi, dp, c0, cL = 0.01, -5e4, 4.2e4, 4.1e4 bcs = [ DirichletBC(phi, "Negative", 0.0), @@ -170,7 +170,9 @@ def test_a_membrane_between_two_potentials_and_pressures(): ] # The pressure's equation is some fifteen orders of magnitude smaller # than the potential's, so the equations are solved in dimensionless form - system = build_system(steady(model), check_units=True, factorize_equations=True) + system = build_system( + create_steady_membrane(model), check_units=True, factorize_equations=True + ) disc = discretize([system, *bcs], DomainMaps({domain: marked})) Solver(disc[system], bcs).solve() @@ -198,7 +200,7 @@ def test_a_membrane_between_two_potentials_and_pressures(): def test_the_membrane_holds_the_protons_at_its_faces(): domain = Domain("cell", "interval") - model, *_ = membrane(domain) + model, *_ = create_membrane(domain) protons = [Field(f"c_H_{side}", domain, ("Lagrange", 1)) for side in "ab"] bcs = model.proton_conditions(protons, ("Negative", "Positive")) assert [bc.field for bc in bcs] == protons From 20a5de1a2c6c1a84a1f2da6c7a3238fe98a462f3 Mon Sep 17 00:00:00 2001 From: Jorgen Schartum Dokken Date: Thu, 12 Nov 2026 12:28:52 +0000 Subject: [PATCH 2/7] Add type-hint helper for expression or value --- src/sumfem/symbolic/components.py | 19 ++++++++++++++++--- 1 file changed, 16 insertions(+), 3 deletions(-) diff --git a/src/sumfem/symbolic/components.py b/src/sumfem/symbolic/components.py index c8f8e6b..cd9bfce 100644 --- a/src/sumfem/symbolic/components.py +++ b/src/sumfem/symbolic/components.py @@ -54,14 +54,27 @@ "JacobianSum", "Residual", "ResidualSum", + "ExpressionOrDefault", + "ExpressionOrValue", "expression_or_parameter", "input_or_parameter", ] +ExpressionOrValue: typing.TypeAlias = ufl.core.expr.Expr | npt.ArrayLike | None +"""A model input: an expression, a value converted with :func:`as_expr`, +or ``None`` for a parameter, see :func:`input_or_parameter`.""" + +ExpressionOrDefault: typing.TypeAlias = ( + ufl.core.expr.Expr | npt.ArrayLike | sympy.Expr | None +) +"""A property: an expression, used as given, or the default of a parameter, +``None`` for none, see :func:`expression_or_parameter`.""" + + def input_or_parameter( name: str, - value: ufl.core.expr.Expr | npt.ArrayLike | None, + value: ExpressionOrValue, domain: Domain, *, shape: tuple[int, ...] = (), @@ -93,7 +106,7 @@ def input_or_parameter( def expression_or_parameter( name: str, - value: ufl.core.expr.Expr | npt.ArrayLike | sympy.Expr | None, + value: ExpressionOrDefault, domain: Domain, *, shape: tuple[int, ...] = (), @@ -216,7 +229,7 @@ def known_coefficient( def input( self, name: str, - value: ufl.core.expr.Expr | npt.ArrayLike | None = None, + value: ExpressionOrValue = None, *, domain: Domain, element: ElementSpec | None = None, From f7e2694ffe00bab362c38a02dbe25c9e482631f6 Mon Sep 17 00:00:00 2001 From: Jorgen Schartum Dokken Date: Thu, 12 Nov 2026 12:28:59 +0000 Subject: [PATCH 3/7] Expose helpers --- src/sumfem/symbolic/__init__.py | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/src/sumfem/symbolic/__init__.py b/src/sumfem/symbolic/__init__.py index fed517f..3f9f317 100644 --- a/src/sumfem/symbolic/__init__.py +++ b/src/sumfem/symbolic/__init__.py @@ -27,6 +27,8 @@ from .bcs import DirichletBC from .components import ( Component, + ExpressionOrDefault, + ExpressionOrValue, Functional, FunctionalSum, Jacobian, @@ -78,6 +80,8 @@ "Clock", "Component", "Functional", + "ExpressionOrDefault", + "ExpressionOrValue", "expression_or_parameter", "input_or_parameter", "FunctionalSum", From c49112bf443163509829bec32154c1441f7d2ae3 Mon Sep 17 00:00:00 2001 From: Jorgen Schartum Dokken Date: Thu, 12 Nov 2026 12:29:15 +0000 Subject: [PATCH 4/7] Use helpers --- src/sumfem/models/flow.py | 53 ++++++++++++++++++++------------------- 1 file changed, 27 insertions(+), 26 deletions(-) diff --git a/src/sumfem/models/flow.py b/src/sumfem/models/flow.py index 9419b0b..928a3b2 100644 --- a/src/sumfem/models/flow.py +++ b/src/sumfem/models/flow.py @@ -11,11 +11,17 @@ from __future__ import annotations -import numpy.typing as npt import sympy.physics.units as syu import ufl -from sumfem.symbolic import Field, Residual, ResidualSum, input_or_parameter +from sumfem.symbolic import ( + ExpressionOrValue, + Field, + Residual, + ResidualSum, + input_or_parameter, +) +from sumfem.symbolic.si import DENSITY, VISCOSITY __all__ = [ "Brinkman", @@ -29,15 +35,10 @@ "Stokes", ] -VISCOSITY = syu.pascal * syu.second """Unit of dynamic viscosities.""" -DENSITY = syu.kg / syu.m**3 """Unit of densities.""" -INPUT = ufl.core.expr.Expr | npt.ArrayLike | None -"""A model input: an expression, a value, or ``None`` for a parameter.""" - class Darcy(Residual): """Darcy flow through a porous medium, ``(K / mu grad(p), grad(q))``. @@ -73,8 +74,8 @@ def __init__( self, name: str, pressure: Field, - permeability: INPUT = None, - viscosity: INPUT = None, + permeability: ExpressionOrValue = None, + viscosity: ExpressionOrValue = None, markers: tuple[str, ...] = (), ): super().__init__(name) @@ -143,7 +144,7 @@ def __init__( name: str, velocity: Field, pressure: Field, - viscosity: INPUT = None, + viscosity: ExpressionOrValue = None, markers: tuple[str, ...] = (), ): super().__init__(name) @@ -243,7 +244,7 @@ def __init__( self, name: str, velocity: Field, - density: INPUT = None, + density: ExpressionOrValue = None, markers: tuple[str, ...] = (), ): super().__init__(name, velocity, markers) @@ -280,8 +281,8 @@ def __init__( self, name: str, velocity: Field, - viscosity: INPUT = None, - permeability: INPUT = None, + viscosity: ExpressionOrValue = None, + permeability: ExpressionOrValue = None, markers: tuple[str, ...] = (), ): super().__init__(name, velocity, markers) @@ -334,8 +335,8 @@ def __init__( self, name: str, velocity: Field, - density: INPUT = None, - coefficient: INPUT = None, + density: ExpressionOrValue = None, + coefficient: ExpressionOrValue = None, markers: tuple[str, ...] = (), speed_floor: ufl.core.expr.Expr | float = 1e-12, ): @@ -399,7 +400,7 @@ def __init__( self, name: str, velocity: Field, - weight: INPUT = None, + weight: ExpressionOrValue = None, markers: tuple[str, ...] = (), ): super().__init__(name) @@ -474,8 +475,8 @@ def __init__( name: str, velocity: Field, pressure: Field, - viscosity: INPUT = None, - density: INPUT = None, + viscosity: ExpressionOrValue = None, + density: ExpressionOrValue = None, markers: tuple[str, ...] = (), ): domain = velocity.domain @@ -521,9 +522,9 @@ def __init__( name: str, velocity: Field, pressure: Field, - porosity: INPUT = None, - permeability: INPUT = None, - viscosity: INPUT = None, + porosity: ExpressionOrValue = None, + permeability: ExpressionOrValue = None, + viscosity: ExpressionOrValue = None, markers: tuple[str, ...] = (), ): domain = velocity.domain @@ -600,11 +601,11 @@ def __init__( name: str, velocity: Field, pressure: Field, - porosity: INPUT = None, - permeability: INPUT = None, - forchheimer: INPUT = None, - viscosity: INPUT = None, - density: INPUT = None, + porosity: ExpressionOrValue = None, + permeability: ExpressionOrValue = None, + forchheimer: ExpressionOrValue = None, + viscosity: ExpressionOrValue = None, + density: ExpressionOrValue = None, markers: tuple[str, ...] = (), speed_floor: float = 1e-12, ): From 8c9f57739609069bc106ca5a74243cad39c64412 Mon Sep 17 00:00:00 2001 From: Jorgen Schartum Dokken Date: Thu, 12 Nov 2026 12:31:24 +0000 Subject: [PATCH 5/7] Second attempt at structuring. Not great.. --- docs/api/electrochemistry.md | 12 +- docs/api/symbolic.md | 18 + .../electrochemistry/vanadium_full_cell_2d.py | 641 ++++++++++++++++++ .../electrochemistry/vanadium_half_cell_1d.py | 1 + src/sumfem/electrochemistry/__init__.py | 12 +- src/sumfem/electrochemistry/electrode.py | 424 ++++++++++++ src/sumfem/electrochemistry/electrolyte.py | 8 +- src/sumfem/electrochemistry/full_cell.py | 553 ++++++--------- src/sumfem/electrochemistry/half_cell.py | 43 +- src/sumfem/electrochemistry/kinetics.py | 46 +- src/sumfem/electrochemistry/membrane.py | 40 +- src/sumfem/electrochemistry/porous.py | 152 +---- src/sumfem/electrochemistry/specs.py | 163 +++++ 13 files changed, 1550 insertions(+), 563 deletions(-) create mode 100644 examples/electrochemistry/vanadium_full_cell_2d.py create mode 100644 src/sumfem/electrochemistry/electrode.py create mode 100644 src/sumfem/electrochemistry/specs.py diff --git a/docs/api/electrochemistry.md b/docs/api/electrochemistry.md index 88e0fb2..e8e7cdb 100644 --- a/docs/api/electrochemistry.md +++ b/docs/api/electrochemistry.md @@ -14,12 +14,21 @@ :members: ``` -## Porous electrodes +## Porous electrodes and their flow ```{eval-rst} .. autoclass:: PorousElectrode :members: +.. autoclass:: ElectrodeFlow + :members: + +.. autoclass:: DarcyFlow + :members: + +.. autoclass:: FlowSetup + :members: + .. autoclass:: Region :members: ``` @@ -93,4 +102,5 @@ .. autoclass:: FlowSpec :members: + :exclude-members: model ``` diff --git a/docs/api/symbolic.md b/docs/api/symbolic.md index 53a02af..d31bbd1 100644 --- a/docs/api/symbolic.md +++ b/docs/api/symbolic.md @@ -58,6 +58,10 @@ .. autofunction:: expression_or_parameter +.. autodata:: ExpressionOrValue + +.. autodata:: ExpressionOrDefault + .. autoclass:: Residual :members: @@ -163,6 +167,20 @@ .. currentmodule:: sumfem.symbolic .. autofunction:: boltzmann_constant + +.. currentmodule:: sumfem.symbolic.si + +.. autodata:: CONCENTRATION + +.. autodata:: DENSITY + +.. autodata:: DIFFUSIVITY + +.. autodata:: PERMEABILITY + +.. autodata:: VISCOSITY + +.. currentmodule:: sumfem.symbolic ``` ## Errors diff --git a/examples/electrochemistry/vanadium_full_cell_2d.py b/examples/electrochemistry/vanadium_full_cell_2d.py new file mode 100644 index 0000000..ec94e88 --- /dev/null +++ b/examples/electrochemistry/vanadium_full_cell_2d.py @@ -0,0 +1,641 @@ +# --- +# jupyter: +# jupytext: +# text_representation: +# extension: .py +# format_name: light +# format_version: '1.5' +# jupytext_version: 1.13.6 +# kernelspec: +# display_name: Python 3 (ipykernel) +# language: python +# name: python3 +# --- + +# # A vanadium redox flow cell in two dimensions +# +# This demo builds a full cell model of Shah, Watt-Smith and Walsh +# {footcite}`Shah2008`, with a slight correction, see +# {ref}`full-cell-corrections`: two half-cells of +# {doc}`vanadium_half_cell_2d`, a proton-exchange membrane between them and +# a current collector on either side, with the electrolytes flowing through +# the electrodes from their tanks. It shows how to: +# +# - Describe each side of the cell with a +# {py:class}`HalfCellSpec `, the +# membrane with a +# {py:class}`MembraneSpec ` and the +# flow with a {py:class}`FlowSpec ` +# - Build the cell with {py:class}`FullCell ` +# on a mesh of its layers, each unknown on the layers where it lives +# - Charge, rest and discharge the cell, with the tanks mixing what the +# electrodes return +# - Check that the cell and its tanks conserve the vanadium, and that each +# couple turns over the charge that passes +# +# ## The model +# +# The cell has five layers, from left to right: the negative current +# collector, the negative porous electrode, the membrane, the positive +# porous electrode and the positive current collector. Each half-cell is +# that of {doc}`vanadium_half_cell_1d`, {eq}`redox-half-cell`, for its couple: +# V²⁺/V³⁺ on the negative side and VO²⁺/VO₂⁺, which takes up water and gives +# off protons, {eq}`butler-volmer-couple`, on the positive side, with the +# piston velocities of Shah et al., {eq}`piston-velocity`. The membrane +# conducts protons only, {eq}`proton-exchange-membrane`, and lets water +# through by Schlögl's flow, {eq}`schloegl-velocity`. The electrolytes +# enter at the bottom of the electrodes and leave at the top, +# {eq}`full-cell-flow`, and return to their tanks, {eq}`full-cell-tank`. +# +# (full-cell-corrections)= +# ### Corrections to Shah et al. +# +# The model is that of Shah et al., with one correction to the model and +# some to the paper. The correction to the model is at the membrane: +# +# - At the faces of the membrane the protons of the electrodes carry its +# current, {eq}`membrane-proton-flux`, and no other ion crosses. Shah et +# al. hold the protons at the fixed charge of the membrane there, which +# lets sulfate, whose concentration follows from electroneutrality, carry +# part of the current through the faces. +# +# The corrections to the paper are: +# +# - Their Table 1 puts the protons that the positive couple makes in the +# negative electrode. +# - Their eq. 22 has the decimal logarithm in place of the natural one, and +# the ratio of V(IV) to itself in the positive potential, which should be +# that of V(V) to V(IV). The text after it gives the formal potentials in +# their Table 3; they are in their Table 4. +# - Their eq. 23 pairs the rate constant of each couple with the formal +# potential of the other; their eq. 24 has them right. +# - Their Table 3 gives the volume of electrolyte of a half-cell in metres. +# Millilitres match the charge time and the state of charge at the end of +# the charge of their Figure 2. +# - Their Table 2 gives an initial water concentration ten times smaller +# than the one their Figure 7 starts from, and their Figure 6 shows. +# - Their Table 5 gives the diffusivities of bisulfate and sulfate in the +# membrane, which holds no anions in their model, and none of the +# protons, bisulfate and sulfate in the electrolyte, which their eqs. 5 +# and 10 need. Those are taken from Knehr et al. {footcite}`Knehr2012`. +# +# The paper leaves two choices open. It does not say how the electrolyte +# potential, the water and the pressure meet across the faces of the +# membrane: they are continuous here. It relates its inlet velocity to the +# flow rate through the porosity, as a velocity in the pores, but uses it in +# Darcy's law, whose velocity is superficial: here it is the flow rate over +# the area of the inlet. The inlets also take in the flux of the +# electrolyte, rather than holding its concentration, so that the cell and +# its tanks exchange exactly what they hold. +# +# The model leaves out the Donnan potentials at the faces of the membrane +# and the protons in the potential of the positive couple, which Shah et al. +# make up for with a fitted contact resistance. Their measured curves are +# not compared with. + +# + +from mpi4py import MPI + +import dolfinx.fem +import dolfinx.mesh +import matplotlib.pyplot as plt +import numpy as np +import pyvista +import sympy.physics.units as syu + +from sumfem import discrete, solvers, symbolic +from sumfem.electrochemistry import ( + DarcyFlow, + FlowSpec, + FullCell, + HalfCellSpec, + Ion, + MembraneSpec, + RedoxCouple, +) +from sumfem.plotting import gather_grid +from sumfem.symbolic.si import to_si + +# - + +# ## The cell +# +# The cell lives on a symbolic domain of quadrilaterals. Its models name the +# parts of the cell by markers, the layers and their faces, so any mesh with +# those markers discretizes it, see +# {py:meth}`mesh_data `. Its +# length scale is the thickness of an electrode, Shah et al.'s Table 3, +# which with the height and the width of the cell also sets the flow below: + +# + +thickness = 4.0 * syu.millimeter +height, width = 10.0 * syu.centimeter, 10.0 * syu.centimeter +domain = symbolic.Domain("cell", "quadrilateral", length_scale=to_si(thickness, syu.m)) +# - + +# ## The sides +# +# Each side has a +# {py:class}`RedoxCouple `, as the +# half-cell of {ref}`vanadium_half_cell_1d `: +# those of Shah et al.'s eqs. 1 and 2, with the formal potentials and the +# rate constants of their Table 4. Their eq. 20 has the porosity times the +# rate constant where the half-cell has its rate constant +# {eq}`butler-volmer-couple`, with the same specific area, so the rate +# constants are multiplied by the porosity. +# +# The positive couple, $\mathrm{VO^{2+} + H_2O \rightleftharpoons VO_2^+ + +# 2 H^+ + e^-}$, also takes up water and gives off protons when it is +# oxidized, on charge. Its `reactants` and `products` give every species of +# the reaction with its coefficient, the ions of the couple included, and +# the couple checks that the reaction conserves charge. The coefficients +# are those of the sources of the ions in each half-cell, +# {eq}`redox-half-cell`: + +# + +porosity = 0.68 +negative_couple = RedoxCouple( + reduced="V2", + oxidized="V3", + rate_constant=porosity * 1.75e-7 * syu.m / syu.second, + formal_potential=-0.255 * syu.volt, + oxidation_on_discharge=True, +) +positive_couple = RedoxCouple( + reduced="VO", + oxidized="VO2", + rate_constant=porosity * 3e-9 * syu.m / syu.second, + formal_potential=1.004 * syu.volt, + oxidation_on_discharge=False, + reactants={"VO": 1, "H2O": 1}, + products={"VO2": 1, "H": 2}, +) +# - + +# Each electrolyte holds its couple, protons, bisulfate and water, with +# sulfate by electroneutrality. The diffusivities of the vanadium ions and +# the water are those of Shah et al.'s Table 5, and those of the protons, +# bisulfate and sulfate those of Knehr et al.'s Table II +# {footcite}`Knehr2012`, see {ref}`full-cell-corrections`: + +# + +diffusivity = syu.m**2 / syu.second + + +def create_ions(reduced, oxidized, charges, couple_diffusivity): + """The ions of a side: its couple, protons, bisulfate, water and sulfate. + + Args: + reduced: The name of the couple's reduced ion. + oxidized: The name of the couple's oxidized ion. + charges: The charge numbers of the reduced and oxidized ion. + couple_diffusivity: The diffusivity of both ions of the couple. + """ + return ( + Ion(reduced, charges[0], couple_diffusivity), + Ion(oxidized, charges[1], couple_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), + ) + + +# - + +# Both sides start at the state of charge of Shah et al.'s Table 2: the +# concentrations of the discharged and the charged ion of each couple, and +# those of the acid and the water, the latter corrected: + +# + +concentration = syu.mol / syu.m**3 +discharged, charged = 1053.0 * concentration, 27.0 * concentration +acid = { + "H": 1200.0 * concentration, + "HSO4": 1200.0 * concentration, + "H2O": 4.2e4 * concentration, # Shah et al., Table 2: 4.2e3 +} +# - + +# A {py:class}`HalfCellSpec ` gathers +# what a side of the cell is made of: the ions of its electrolyte and their +# initial concentrations, its couple, and the properties of its electrode. +# The full cell makes a +# {py:class}`RedoxHalfCell ` of each +# side, see {doc}`vanadium_half_cell_1d`. Both electrodes are the carbon +# felt of Shah et al.'s Tables 3 and 5: its porosity, the conductivity of +# its fibres and their specific area, the distance between the fibres, +# which sets the piston velocities, {eq}`piston-velocity`, and their +# diameter and the Kozeny-Carman constant, which set the permeability of the +# flow, {py:func}`kozeny_carman `. The electrolytes have +# the viscosity of water, as theirs: + +# + +carbon_felt = dict( + porosity=porosity, + conductivity=500.0 * syu.siemens / syu.m, + specific_area=2e6 / syu.m, + fiber_distance=10.0 * syu.micrometer, + fiber_diameter=10.0 * syu.micrometer, + kozeny_carman_constant=5.55, + viscosity=1e-3 * syu.pascal * syu.second, +) +negative = HalfCellSpec( + ions=create_ions("V2", "V3", (2, 3), 2.4e-10 * diffusivity), + reservoir={"V2": charged, "V3": discharged, **acid}, + couple=negative_couple, + **carbon_felt, +) +positive = HalfCellSpec( + ions=create_ions("VO", "VO2", (2, 1), 3.9e-10 * diffusivity), + reservoir={"VO": discharged, "VO2": charged, **acid}, + couple=positive_couple, + **carbon_felt, +) +# - + +# ## The membrane and the flow +# +# A {py:class}`MembraneSpec ` gives the +# properties of the +# {py:class}`ProtonExchangeMembrane `: +# the concentration of its fixed charges and the diffusivity of its +# protons, which set its conductivity, {eq}`proton-exchange-membrane`, the +# diffusivity of its water, and the permeabilities and the viscosity of +# Schlögl's flow through it, {eq}`schloegl-velocity`. Those of Shah et al.'s +# Tables 4 and 5: + +membrane = MembraneSpec( + fixed_charge=1200.0 * concentration, + proton_diffusivity=1.4e-9 * diffusivity, + water_diffusivity=5.75e-10 * diffusivity, + electrokinetic_permeability=1.13e-19 * syu.m**2, + hydraulic_permeability=1.58e-18 * syu.m**2, + viscosity=1e-3 * syu.pascal * syu.second, +) + +# A {py:class}`FlowSpec ` adds the flow of +# the electrolytes to the cell. Its `model` is the flow, an +# {py:class}`ElectrodeFlow `, which +# gives the superficial velocity of each electrode: any flow plugs in, e.g. +# Brinkman-Forchheimer or transient Stokes flow. Shah et al.'s is +# {py:class}`DarcyFlow `: one pressure +# lives on the electrodes and the membrane, with +# {py:class}`Darcy ` flow in the electrodes, +# {eq}`darcy`, and +# {py:class}`Schlögl's ` in the +# membrane, {eq}`schloegl-velocity`. +# The electrolyte of each side enters its electrode at the inlet, at the +# inlet velocity and with the concentrations of its tank, and leaves at the +# outlet, held at the outlet pressure, {eq}`full-cell-flow`. Each tank, +# whose volume is given relative to the pores of its electrode, mixes what +# the outlet returns, {eq}`full-cell-tank`. +# +# Shah et al. pump each electrolyte at the flow rate of their Table 2. The +# inlet velocity is the superficial velocity of that flow rate, over the +# inlet, the thickness of the electrode times the width of the cell. Each +# tank holds the volume of electrolyte of a half-cell, their Table 3, but +# for what fills the pores of the electrode: + +# + +flow_rate, half_cell_volume = 1.0 * syu.milliliter / syu.second, 250 * syu.milliliter +inlet = thickness * width +pores = to_si(porosity * height * inlet, syu.m**3) +tank = to_si(half_cell_volume, syu.m**3) - pores +flow = FlowSpec( + inlet_velocity=flow_rate / inlet, tank_ratio=tank / pores, model=DarcyFlow +) +# - + +# {py:class}`FullCell ` builds the cell +# from its sides, its membrane and its flow, with the conductivity of the +# collectors of Shah et al.'s Table 5. Its Dirichlet conditions, +# {py:meth}`conditions `, ground +# the negative terminal and hold the pressure at the outlets: + +# + +cell = FullCell( + domain, + negative, + positive, + membrane=membrane, + collector_conductivity=1000.0 * syu.siemens / syu.m, + flow=flow, +) +conditions = cell.conditions() +# - + +# ## The mesh +# +# The cell is discretized on the geometry of Shah et al., their Table 3: the +# layers with their thicknesses and the number of columns of cells across +# each, and the rows of cells over the height of the cell, from the inlets +# to the outlets: + +# + +layers = ( + ("NegativeCollector", 6.0 * syu.millimeter, 3), + ("NegativeElectrode", thickness, 16), + ("Membrane", 180.0 * syu.micrometer, 3), + ("PositiveElectrode", thickness, 16), + ("PositiveCollector", 6.0 * syu.millimeter, 3), +) +rows = 40 +# - + +# A rectangle of quadrilaterals, whose columns are moved onto the edges of +# the layers, with the markers of the cell: the cells of each layer, the +# outer faces of the collectors, the faces of the membrane and the ends of +# the electrodes. They are located with a tolerance relative to the thinnest layer, +# see {ref}`units-length-scale`. The cell's unknowns live on unions of the +# layers, which +# {py:meth}`mesh_data ` adds: + +# + +edges = np.cumsum([0.0, *(to_si(thickness, syu.m) for _, thickness, _ in layers)]) +columns = np.unique( + np.concatenate( + [ + np.linspace(a, b, n + 1) + for a, b, (*_, n) in zip(edges[:-1], edges[1:], layers, strict=True) + ] + ) +) +dtype = dolfinx.default_real_type +rectangle = dolfinx.mesh.create_rectangle( + MPI.COMM_WORLD, + np.array([[0.0, 0.0], [1.0, to_si(height, syu.m)]]), + [len(columns) - 1, rows], + cell_type=dolfinx.mesh.CellType.quadrilateral, + dtype=dtype, +) +index = np.rint(rectangle.geometry.x[:, 0] * (len(columns) - 1)).astype(int) +rectangle.geometry.x[:, 0] = columns[index] +tolerance = 100 * float(np.finfo(dtype).eps) * to_si(layers[2][1], syu.m) + + +def between(a, b): + """The points from ``x = a`` to ``x = b``. + + Args: + a: The left edge, in m. + b: The right edge, in m. + """ + return lambda x: (x[0] >= a - tolerance) & (x[0] <= b + tolerance) + + +def at(a): + """The points at ``x = a``. + + Args: + a: The position, in m. + """ + return lambda x: np.isclose(x[0], a, atol=tolerance) + + +def end(a, b, y): + """The points from ``x = a`` to ``x = b`` at the height ``y``. + + Args: + a: The left edge, in m. + b: The right edge, in m. + y: The height, in m. + """ + return lambda x: between(a, b)(x) & np.isclose(x[1], y, atol=tolerance) + + +top = to_si(height, syu.m) +marked = FullCell.mesh_data( + discrete.MeshData.from_locators( + rectangle, + { + 2: { + name: between(a, b) + for (name, *_), a, b in zip(layers, edges[:-1], edges[1:], strict=True) + }, + 1: { + "NegativeTerminal": at(edges[0]), + "NegativeFace": at(edges[2]), + "PositiveFace": at(edges[3]), + "PositiveTerminal": at(edges[5]), + "NegativeInlet": end(edges[1], edges[2], 0.0), + "NegativeOutlet": end(edges[1], edges[2], top), + "PositiveInlet": end(edges[3], edges[4], 0.0), + "PositiveOutlet": end(edges[3], edges[4], top), + }, + }, + ) +) +# - + +# ## The discrete problem +# +# The discretization pairs the symbolic domain with the marked mesh. It +# bounds the quadrature degree of the nonlinear kinetics, as in +# {doc}`vanadium_half_cell_1d`, whose integrals no degree makes exact; a +# higher bound changes the voltage and the states of charge only at +# round-off. The initial state is at rest, Shah et al.'s eq. 32: the +# concentrations of the tanks, the electrolyte at zero and each solid at the +# equilibrium potential of its couple: + +disc = discrete.discretize( + [cell.model, *conditions], + discrete.DomainMaps({domain: marked}), + max_quadrature_degree=4, +) +cell.initialize(disc) + +# ## Charge, rest and discharge +# +# The current of Shah et al.'s Table 2 enters and leaves through the outer +# faces of the collectors, their eq. 31, whose area is the height times the +# width of the cell. The cell is charged for the time of +# their Figure 2, rested, and discharged until its voltage falls to a +# cut-off or the time of the charge has passed. Each operation gives its +# mode and its time steps, as pairs of a number of steps and their size in +# seconds: + +# + +current = 10.0 * syu.ampere / (height * width) +density = to_si(current, syu.ampere / syu.m**2) +charging = ((7, 1.0), (10, 5.0), (98, 20.0)) +schedule = ( + ("charge", -density, charging), + ("rest", 0.0, ((12, 10.0),)), + ("discharge", density, charging), +) +cutoff = 1.0 +largest_step = max(dt for *_, steps in schedule for _, dt in steps) +# - + +# A {py:class}`Simulation ` takes the steps with +# backward Euler and solves each in dimensionless form, see +# {ref}`units-scaling`. After every step, a callback takes a step of the +# tanks, {py:meth}`update_tanks +# `, and records the time, +# the mode, the voltage and the state of charge of each side, tank included, +# {py:meth}`state_of_charge +# `: + +# + +clock = symbolic.Clock(domain, dt=schedule[0][2][0][1]) +simulation = solvers.Simulation( + disc, + cell.model, + conditions, + clock=clock, + scheme=solvers.BDF(1), + petsc_options={"snes_linesearch_type": "bt", "snes_max_it": 50}, + check_units=True, + factorize_equations=True, + time_scale=largest_step, +) +records: list[tuple[float, str, float, float, float]] = [] + + +def record(t: float, mode: str, dt: float) -> None: + """Take a step of the tanks, and record the state of the cell. + + Args: + t: The time, in s. + mode: The operation. + dt: The size of the step, in s. + """ + cell.update_tanks(disc, dt) + records.append( + ( + t, + mode, + cell.voltage(disc), + cell.state_of_charge(disc, "negative"), + cell.state_of_charge(disc, "positive"), + ) + ) + + +# - + +# The amounts of vanadium on each side, in the electrode and its tank, before +# the schedule, and the charge that passes through the cell, per unit of its +# width, are kept for the checks below. At the end of the charge, the local +# state of charge of the negative electrode is kept for a plot: + + +# + +def vanadium(side: str) -> float: + """The amount of vanadium of a side, in its electrode and its tank. + + Args: + side: ``"negative"`` or ``"positive"``. + """ + couple = cell.half_cells[side].couple + return sum( + cell.amount(disc, side, ion) for ion in (couple.reduced, couple.oxidized) + ) + + +start = {side: vanadium(side) for side in ("negative", "positive")} +reduced = cell.amount(disc, "negative", "V2") +negative_cell = cell.half_cells["negative"] +local = negative_cell.local_state_of_charge() +space = disc.function_space(negative_cell.concentrations["V2"]) +local_charge = dolfinx.fem.Function(space, name="soc") +passed = {"charge": 0.0, "discharge": 0.0} +for mode, density_of_mode, steps in schedule: + disc.set({cell.current: density_of_mode}) + stop = False + for count, dt in steps: + disc.set({clock.dt: dt}) + for _ in range(count): + simulation.run( + num_steps=1, + callback=lambda t, step, mode=mode, dt=dt: record(t, mode, dt), + ) + if mode in passed: + passed[mode] += abs(density_of_mode) * dt * top + if mode == "discharge" and records[-1][2] < cutoff: + stop = True + break + if stop: + break + if mode == "charge": + local_charge.interpolate( + dolfinx.fem.Expression(disc[local], space.element.interpolation_points) + ) +times = np.array([r[0] for r in records]) +modes = np.array([r[1] for r in records]) +voltage = np.array([r[2] for r in records]) +soc = { + "negative": np.array([r[3] for r in records]), + "positive": np.array([r[4] for r in records]), +} +# - + +# ## Checks +# +# The cell and its tanks keep the vanadium of each side, and the negative +# couple turns over the charge that passed, divided by Faraday's constant. +# At rest the states of charge, tanks included, stay as they are, while the +# flow mixes the electrolytes of the electrodes and the tanks. Sulfate, +# which follows from electroneutrality, stays positive on both sides, see +# {py:meth}`minimum_balanced_concentration +# `: + +# + +faraday = float(disc.value(negative_cell.electrolyte.F)) +turned_over = (passed["charge"] - passed["discharge"]) / faraday +for side in ("negative", "positive"): + print( + f"Vanadium of the {side} side: {vanadium(side) / start[side]:.12f} of the start" + ) + assert np.isclose(vanadium(side), start[side], rtol=1e-9) +change = cell.amount(disc, "negative", "V2") - reduced +print(f"V2+ made over the charge passed: {change / turned_over:.9f}") +assert np.isclose(change, turned_over, rtol=1e-6) +rest = soc["negative"][modes == "rest"] +assert np.ptp(rest) < 1e-9 * rest.mean() +assert np.all(np.diff(soc["negative"][modes == "charge"]) > 0) +assert np.all(np.diff(soc["negative"][modes == "discharge"]) < 0) +assert voltage[modes == "charge"].mean() > voltage[modes == "discharge"].mean() +for half in cell.half_cells.values(): + assert half.minimum_balanced_concentration(disc) > 0 +# - + +# The voltage and the states of charge through the schedule. Charging raises +# the voltage above the difference of the formal potentials of the couples, +# and discharging lowers it below, by the overpotentials and the Ohmic +# losses of the cell. The states of charge of the two sides follow each +# other, since each electron that one couple takes the other gives: + +fig, (upper, lower) = plt.subplots(2, 1, sharex=True) +for mode in ("charge", "rest", "discharge"): + upper.plot(times[modes == mode], voltage[modes == mode], ".", label=mode) + lower.plot(times[modes == mode], soc["negative"][modes == mode], ".", label=mode) +lower.plot(times, soc["positive"], "k--", label="positive side") +upper.set_ylabel("Voltage (V)") +lower.set_ylabel("State of charge") +lower.set_xlabel("Time (s)") +upper.legend() +lower.legend() +plt.show() + +# The local state of charge of the negative electrode at the end of the +# charge, gathered from all processes with +# {py:func}`gather_grid `. The electrolyte +# flows up from the inlet at the bottom, so it is charged more towards the +# outlet. The electrode is drawn in millimetres, widened to show its +# thickness: + +grid = gather_grid(local_charge) +if MPI.COMM_WORLD.rank == 0: + grid.points[:, 0] *= 1e4 + grid.points[:, 1] *= 1e3 + plotter = pyvista.Plotter(window_size=[480, 640]) + plotter.add_mesh(grid, scalars="soc", scalar_bar_args={"title": "SOC"}) + plotter.view_xy() + plotter.show() + +# ## References +# +# ```{footbibliography} +# ``` diff --git a/examples/electrochemistry/vanadium_half_cell_1d.py b/examples/electrochemistry/vanadium_half_cell_1d.py index 53ce2e5..812cdbe 100644 --- a/examples/electrochemistry/vanadium_half_cell_1d.py +++ b/examples/electrochemistry/vanadium_half_cell_1d.py @@ -185,6 +185,7 @@ } # - +# (vanadium-half-cell-couple)= # A {py:class}`RedoxCouple ` gives the # reaction $\mathrm{V^{2+} \rightleftharpoons V^{3+} + e^-}$ of # {eq}`half-cell-kinetics`: the ion `reduced`, V²⁺, the ion `oxidized`, V³⁺, diff --git a/src/sumfem/electrochemistry/__init__.py b/src/sumfem/electrochemistry/__init__.py index fc9712d..179e742 100644 --- a/src/sumfem/electrochemistry/__init__.py +++ b/src/sumfem/electrochemistry/__init__.py @@ -18,7 +18,8 @@ - Full cells: :class:`FullCell`, two half-cells, each a :class:`HalfCellSpec`, and a membrane between them, a :class:`MembraneSpec`, with the flow of the electrolytes from their - tanks, a :class:`FlowSpec`. + tanks, a :class:`FlowSpec`, by any :class:`ElectrodeFlow`, e.g. + :class:`DarcyFlow`. The package imports from the core of :mod:`sumfem`, and never the other way around. @@ -26,8 +27,9 @@ from .charge import ElectrolyteConduction from .couple import RedoxCouple +from .electrode import DarcyFlow, ElectrodeFlow, FlowSetup, PorousElectrode from .electrolyte import Electrolyte, Ion -from .full_cell import FlowSpec, FullCell, HalfCellSpec, MembraneSpec +from .full_cell import FullCell from .half_cell import RedoxHalfCell from .kinetics import ( ButlerVolmer, @@ -37,7 +39,8 @@ PistonVelocityButlerVolmer, ) from .membrane import MembraneProtonFlux, ProtonExchangeMembrane, SchloeglFlow -from .porous import PorousElectrode, Region +from .porous import Region +from .specs import FlowSpec, HalfCellSpec, MembraneSpec from .transport import NernstPlanck, TransientNernstPlanck __all__ = [ @@ -57,6 +60,9 @@ "HalfCellSpec", "MembraneSpec", "FlowSpec", + "ElectrodeFlow", + "DarcyFlow", + "FlowSetup", "ButlerVolmer", "ElectrodeReaction", "ExchangeCurrentButlerVolmer", diff --git a/src/sumfem/electrochemistry/electrode.py b/src/sumfem/electrochemistry/electrode.py new file mode 100644 index 0000000..b0f440c --- /dev/null +++ b/src/sumfem/electrochemistry/electrode.py @@ -0,0 +1,424 @@ +"""Porous electrodes, and the flow of the electrolyte through them. + +:class:`PorousElectrode` is an electrode of carbon fibres filled with +electrolyte: its porosity, the conductivity of its fibres, its specific +area and the permeability of its fibres to the flow. The flow of the +electrolytes through the electrodes of a +:class:`~sumfem.electrochemistry.FullCell` is a component, an +:class:`ElectrodeFlow`, which the cell creates from a :class:`FlowSetup` of +what it knows. :class:`DarcyFlow` is Darcy flow through the electrodes; +another flow, e.g. Brinkman-Forchheimer flow or transient Stokes flow, is +another subclass, given to the cell's +:class:`~sumfem.electrochemistry.FlowSpec`. +""" + +from __future__ import annotations + +import abc +import dataclasses +import typing + +import sympy +import sympy.physics.units as syu +import ufl + +from sumfem.discrete import Discretization +from sumfem.error_handling import ComponentError +from sumfem.models import ( + BoundaryFlux, + Darcy, + bruggeman, + ergun_coefficient, + kozeny_carman, +) +from sumfem.symbolic import ( + DirichletBC, + Domain, + Field, + Parameter, + Residual, + Subdomain, + value_at_defaults, +) + +from .porous import Region + +__all__ = ["DarcyFlow", "ElectrodeFlow", "FlowSetup", "PorousElectrode"] + + +class PorousElectrode: + """A porous electrode of carbon fibers, filled with electrolyte. + + The electrode creates a parameter for each argument below from the + porosity on, named ``"."``, such as + ``".porosity"``, with the argument as its default value. Models + that share an electrode share these parameters. A parameter without a + default value must be set on the discretization before it is used. + + Args: + name: Prefix of the parameter names. + domain: The domain of the parameters. + markers: Cell markers of the electrode. Empty for the whole domain. + porosity: Default porosity, if any. Every default with a unit is in + SI units, or a quantity with units, e.g. ``0.5 * syu.micrometer``. + conductivity: Default conductivity of the fiber material in S/m, if + any. + specific_area: Default specific area in 1/m, if any. + fiber_diameter: Default fiber diameter in m, if any. + kozeny_carman_constant: Default Kozeny-Carman constant, if any. + pore_length: Default characteristic pore length of Ergun's drag in + m, if any. + """ + + _markers: tuple[str, ...] + """Cell markers of the electrode.""" + + _porosity: Parameter + """Porosity.""" + + _conductivity: Parameter + """Conductivity of the fiber material.""" + + _specific_area: Parameter + """Specific area.""" + + _fiber_diameter: Parameter + """Fiber diameter.""" + + _kozeny_carman_constant: Parameter + """Kozeny-Carman constant.""" + + _pore_length: Parameter + """Characteristic pore length of Ergun's drag.""" + + def __init__( + self, + name: str, + domain: Domain, + markers: tuple[str, ...] = (), + porosity: float | None = None, + conductivity: float | sympy.Expr | None = None, + specific_area: float | sympy.Expr | None = None, + fiber_diameter: float | sympy.Expr | None = None, + kozeny_carman_constant: float | None = None, + pore_length: float | sympy.Expr | None = None, + ): + self._markers = tuple(markers) + self._porosity = Parameter( + f"{name}.porosity", domain, default=porosity, doc="Porosity" + ) + self._conductivity = Parameter( + f"{name}.conductivity", + domain, + default=conductivity, + unit=syu.siemens / syu.m, + doc="Conductivity of the fiber material", + ) + self._specific_area = Parameter( + f"{name}.specific_area", + domain, + default=specific_area, + unit=1 / syu.m, + doc="Specific area", + ) + self._fiber_diameter = Parameter( + f"{name}.fiber_diameter", + domain, + default=fiber_diameter, + unit=syu.m, + doc="Fiber diameter", + ) + self._kozeny_carman_constant = Parameter( + f"{name}.kozeny_carman_constant", + domain, + default=kozeny_carman_constant, + doc="Kozeny-Carman constant", + ) + self._pore_length = Parameter( + f"{name}.pore_length", + domain, + default=pore_length, + unit=syu.m, + doc="Characteristic pore length of Ergun's drag", + ) + + @property + def markers(self) -> tuple[str, ...]: + """Cell markers of the electrode.""" + return self._markers + + @property + def region(self) -> Region: + """The electrolyte-filled pores, as a region.""" + return Region(self._markers, self._porosity) + + @property + def porosity(self) -> Parameter: + """Porosity.""" + return self._porosity + + @property + def specific_area(self) -> Parameter: + """Specific area.""" + return self._specific_area + + @property + def fiber_diameter(self) -> Parameter: + """Fiber diameter.""" + return self._fiber_diameter + + def effective_conductivity(self) -> ufl.core.expr.Expr: + r"""Bruggeman's conductivity of the solid. + + That is :math:`\sigma_\text{eff} = (1 - \varepsilon)^{3/2} \sigma`, with + the conductivity :math:`\sigma` of the fibres and the porosity + :math:`\varepsilon`. + """ + return bruggeman(self._conductivity, 1.0 - self._porosity) + + def permeability(self) -> ufl.core.expr.Expr: + """Kozeny-Carman permeability of the fibers, see + :func:`~sumfem.models.kozeny_carman`.""" + return kozeny_carman( + self._fiber_diameter, self._porosity, self._kozeny_carman_constant + ) + + def forchheimer(self) -> ufl.core.expr.Expr: + """Ergun's inertial drag coefficient, an estimate for fibres, see + :func:`~sumfem.models.ergun_coefficient`.""" + return ergun_coefficient(self._porosity, self._pore_length) + + +@dataclasses.dataclass(frozen=True) +class FlowSetup: + """What a full cell gives the flow of its electrolytes. + + Attributes: + name: Prefix of the flow's models and parameters, the cell's name. + domain: The domain of the cell. + electrolyte: The subdomain of the electrolyte, the electrodes and + what lies between them, where the pressure lives. + electrodes: The porous electrodes, by side, ``"negative"`` and + ``"positive"``, with their cell markers and properties, e.g. the + permeability of their fibres. + subdomains: The subdomain of each electrode, by side, e.g. for a + velocity field of the electrode. + viscosities: The dynamic viscosity of each side's electrolyte. + densities: The density of each side's electrolyte. + inlets: The facet marker of each electrode's inlet, by side. + outlets: The facet marker of each electrode's outlet, by side. + inlet_velocity: The superficial velocity at the inlets. + outlet_pressure: The gauge pressure at the outlets. + degree: The degree of the Lagrange elements of the cell's unknowns. + """ + + name: str + domain: Domain + electrolyte: Subdomain + electrodes: typing.Mapping[str, PorousElectrode] + subdomains: typing.Mapping[str, Subdomain] + viscosities: typing.Mapping[str, ufl.core.expr.Expr] + densities: typing.Mapping[str, ufl.core.expr.Expr] + inlets: typing.Mapping[str, str] + outlets: typing.Mapping[str, str] + inlet_velocity: ufl.core.expr.Expr + outlet_pressure: ufl.core.expr.Expr + degree: int + + +class ElectrodeFlow(abc.ABC): + """The flow of the electrolytes through the electrodes of a full cell. + + A flow gives the cell the superficial velocity :math:`\\mathbf{u}` of + each electrode, :meth:`velocity`, which carries the ions in the + electrode, out of its outlet and to its tank, and its own models and + conditions, which the cell adds to its own. The cell holds the ions' + inflow at the inlets itself, :eq:`full-cell-flow`, so the flow needs to + take in the electrolyte there at the inlet velocity of the + :class:`FlowSetup`. + + A subclass creates its unknowns and models from the setup in its + constructor. Its models may hold time derivatives, e.g. of the velocity + of transient Stokes flow, which the cell's simulation steps in time + with the rest. Its :attr:`pressure` lives on + :attr:`FlowSetup.electrolyte`, which the cell shares with what lies + between the electrodes, see :class:`~sumfem.electrochemistry.FullCell`. + + Example: + Darcy flow with a storage of the pressure, a flow that changes in + time, is :class:`DarcyFlow` with one more model, given to the cell + by its :class:`~sumfem.electrochemistry.FlowSpec`:: + + class StoredDarcyFlow(DarcyFlow): + @property + def models(self): + storage = Accumulation( + f"{self.setup.name}_storage", + self.pressure, + weight=storativity, + markers=("NegativeElectrode", "PositiveElectrode"), + ) + return [*super().models, storage] + + flow = FlowSpec(inlet_velocity=v_in, model=StoredDarcyFlow) + + A Brinkman-Forchheimer or a transient Stokes flow creates a velocity + field on each electrode's subdomain, :attr:`FlowSetup.subdomains`, + and its pressure on the electrolyte's, + :attr:`FlowSetup.electrolyte`, with + :class:`~sumfem.models.BrinkmanForchheimer` or + :class:`~sumfem.models.Stokes` and an + :class:`~sumfem.models.Accumulation` of the velocity, returns them + from :meth:`velocity` and :attr:`pressure`, and holds the inlet + velocity on the velocity in :meth:`conditions`. + + Args: + setup: What the cell gives the flow. + """ + + _setup: FlowSetup + """What the cell gives the flow.""" + + def __init__(self, setup: FlowSetup): + self._setup = setup + + @property + def setup(self) -> FlowSetup: + """What the cell gives the flow.""" + return self._setup + + @abc.abstractmethod + def velocity(self, side: str) -> ufl.core.expr.Expr: + """The superficial velocity in a side's electrode. + + Args: + side: ``"negative"`` or ``"positive"``. + """ + + @property + @abc.abstractmethod + def models(self) -> list[Residual]: + """The models of the flow, which the cell adds to its own.""" + + @property + @abc.abstractmethod + def pressure(self) -> Field: + """The pressure of the flow, on :attr:`FlowSetup.electrolyte`.""" + + def conditions(self) -> list[DirichletBC]: + """The Dirichlet conditions of the flow, none by default.""" + return [] + + def initialize(self, discretization: Discretization) -> None: + """Set the initial state of the flow's unknowns, none by default. + + Args: + discretization: The discretization of the cell. + """ + + +class DarcyFlow(ElectrodeFlow): + """Darcy flow through the electrodes, with one pressure. + + The pressure :math:`p` lives on :attr:`FlowSetup.electrolyte`, see + :attr:`ElectrodeFlow.pressure`. In each electrode it follows + :class:`~sumfem.models.Darcy` flow, ``"__darcy"``, with the + Kozeny-Carman permeability of its fibres, + :meth:`PorousElectrode.permeability`, and the viscosity of its + electrolyte. + + The electrolyte enters each electrode at its inlet at the inlet + velocity, ``"__inflow"``, a + :class:`~sumfem.models.BoundaryFlux`, and leaves at its outlet, held at + the outlet pressure, :meth:`conditions`. + + The pressure's scale is the pressure drop of a uniform flow at the + inlet velocity through the negative electrode over the domain's length + scale, if the defaults give it, else one pascal. + + Args: + setup: What the cell gives the flow. + """ + + _pressure: Field + """The pressure, on the electrolyte's subdomain.""" + + _darcy: dict[str, Darcy] + """The flow in each electrode, by side.""" + + def __init__(self, setup: FlowSetup): + super().__init__(setup) + negative = setup.electrodes["negative"] + try: + scale = float( + value_at_defaults( + setup.viscosities["negative"] + * setup.inlet_velocity + * setup.domain.length_scale + / negative.permeability() + ) + ) + except ComponentError: + scale = 1.0 + self._pressure = Field( + "p", + setup.electrolyte, + ("Lagrange", setup.degree), + unit=syu.pascal, + scale=scale, + ) + self._darcy = { + side: Darcy( + f"{setup.name}_{side}_darcy", + self._pressure, + electrode.permeability(), + setup.viscosities[side], + markers=electrode.markers, + ) + for side, electrode in setup.electrodes.items() + } + + def velocity(self, side: str) -> ufl.core.expr.Expr: + """The superficial velocity of the side's Darcy flow. + + Args: + side: ``"negative"`` or ``"positive"``. + """ + return self._darcy[side].superficial_velocity + + @property + def models(self) -> list[Residual]: + """The Darcy flow and the inflow of each electrode.""" + setup = self._setup + models: list[Residual] = [] + for side, darcy in self._darcy.items(): + models.append(darcy) + models.append( + BoundaryFlux( + f"{setup.name}_{side}_inflow", + self._pressure, + setup.inlet_velocity, + (setup.inlets[side],), + ) + ) + return models + + @property + def pressure(self) -> Field: + """The pressure, on :attr:`FlowSetup.electrolyte`.""" + return self._pressure + + def conditions(self) -> list[DirichletBC]: + """The outlet pressure at both outlets.""" + outlets = tuple(self._setup.outlets.values()) + return [DirichletBC(self._pressure, outlets, self._setup.outlet_pressure)] + + def initialize(self, discretization: Discretization) -> None: + """Set the pressure to the outlet pressure. + + Args: + discretization: The discretization of the cell. + """ + outlet = self._setup.outlet_pressure + assert isinstance(outlet, Parameter) + discretization.set({self._pressure: float(discretization.value(outlet))}) diff --git a/src/sumfem/electrochemistry/electrolyte.py b/src/sumfem/electrochemistry/electrolyte.py index 7c89534..0fc57d4 100644 --- a/src/sumfem/electrochemistry/electrolyte.py +++ b/src/sumfem/electrochemistry/electrolyte.py @@ -18,16 +18,10 @@ Parameter, expression_or_parameter, ) -from sumfem.symbolic.si import to_si +from sumfem.symbolic.si import CONCENTRATION, DIFFUSIVITY, to_si __all__ = ["Electrolyte", "Ion"] -CONCENTRATION = syu.mol / syu.m**3 -"""Unit of concentrations.""" - -DIFFUSIVITY = syu.m**2 / syu.second -"""Unit of diffusion coefficients.""" - @dataclasses.dataclass(frozen=True) class Ion: diff --git a/src/sumfem/electrochemistry/full_cell.py b/src/sumfem/electrochemistry/full_cell.py index 73071ea..3e34a7e 100644 --- a/src/sumfem/electrochemistry/full_cell.py +++ b/src/sumfem/electrochemistry/full_cell.py @@ -10,173 +10,37 @@ 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.error_handling import ComponentError -from sumfem.models import BoundaryFlux, Conduction, Darcy +from sumfem.models import BoundaryFlux, Conduction from sumfem.symbolic import ( DirichletBC, Domain, + ExpressionOrDefault, Field, Parameter, Residual, ResidualSum, Subdomain, expression_or_parameter, - value_at_defaults, ) +from sumfem.symbolic.si import CONCENTRATION, DENSITY, VISCOSITY -from .couple import RedoxCouple -from .electrolyte import CONCENTRATION, Ion +from .electrode import ElectrodeFlow, FlowSetup, PorousElectrode from .half_cell import RedoxHalfCell -from .membrane import PROPERTY, ProtonExchangeMembrane -from .porous import PorousElectrode +from .membrane import ProtonExchangeMembrane +from .specs import PROTONS, WATER, FlowSpec, HalfCellSpec, MembraneSpec -__all__ = ["FlowSpec", "FullCell", "HalfCellSpec", "MembraneSpec"] +__all__ = ["FullCell"] 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.""" - -VISCOSITY = syu.pascal * syu.second -"""Unit of dynamic viscosities.""" - - -@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`, and with flow, that of its tank. 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 them, limited by the flow's mass - transfer with flow. - fiber_diameter: The diameter of the electrode's fibres, in m or with - units, for the permeability with flow, and the mass transfer - without a fibre distance. - kozeny_carman_constant: The Kozeny-Carman constant of the - permeability with flow, see - :func:`~sumfem.models.kozeny_carman`. - density: The density of the electrolyte, in kg/m^3 or with units. - viscosity: The dynamic viscosity of the electrolyte, in Pa s or with - units, which the flow and the mass transfer share. - """ - - 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 - fiber_diameter: float | sympy.Expr | None = None - kozeny_carman_constant: float | None = None - density: float | sympy.Expr = 1399.8 - viscosity: float | sympy.Expr = 6.963e-3 - - -@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`. The - last three are those of the flow of water, with a :class:`FlowSpec`. - - 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. - electrokinetic_permeability: The electrokinetic permeability, in - m^2. - hydraulic_permeability: The hydraulic permeability, in m^2. - viscosity: The viscosity of the water in the membrane, in Pa s. - """ - - fixed_charge: PROPERTY = None - proton_diffusivity: PROPERTY = None - water_diffusivity: PROPERTY = None - electrokinetic_permeability: PROPERTY = None - hydraulic_permeability: PROPERTY = None - viscosity: PROPERTY = None - - -@dataclasses.dataclass(frozen=True) -class FlowSpec: - r"""The flow of the electrolytes through the electrodes, from their tanks. - - One pressure :math:`p` lives on both electrodes and the membrane: - :class:`~sumfem.models.Darcy` flow in each electrode, with the - Kozeny-Carman permeability of its fibres and the viscosity of its - electrolyte, and :class:`~sumfem.electrochemistry.SchloeglFlow` in the - membrane, so the pressure is continuous across the faces. The - electrolyte of each side enters its electrode at the inlet with the - superficial velocity :math:`v_\text{in}` and the concentrations - :math:`c_k^\text{in}` of its tank, and leaves at the outlet, held at the - pressure :math:`p_\text{out}`, with the flow, - - .. math:: - :label: full-cell-flow - - \begin{aligned} - -\mathbf{u} \cdot \mathbf{n} &= v_\text{in} - && \text{at the inlet}, \\ - -\mathbf{N}_k \cdot \mathbf{n} &= v_\text{in} c_k^\text{in} - && \text{at the inlet}, \\ - p &= p_\text{out} - && \text{at the outlet}, \\ - (\mathbf{N}_k - \mathbf{u} c_k) \cdot \mathbf{n} &= 0 - && \text{at the outlet}, - \end{aligned} - - with the superficial velocity :math:`\mathbf{u}` of the electrode and - the Nernst-Planck flux :math:`\mathbf{N}_k` of each ion. The inlet holds - the flux of the ions, the flow's, rather than their concentration as - Shah et al. do, so that the cell takes in exactly what its tank gives. - - Attributes: - inlet_velocity: The superficial velocity at the inlets of both - electrodes: an expression, or the default of the parameter - ``".inlet_velocity"``, in m/s or with units. - outlet_pressure: The gauge pressure at the outlets, given like the - inlet velocity, in Pa. - tank_ratio: The volume of each tank over the volume of electrolyte in - the pores of its electrode, ``None`` for tanks so large that the - electrolyte at the inlets does not change. - """ - - inlet_velocity: PROPERTY = None - outlet_pressure: PROPERTY = 0.0 - tank_ratio: float | None = None - class FullCell: r"""The full cell of a redox flow battery. @@ -203,9 +67,11 @@ class FullCell: - Ohmic conduction in the collectors, :class:`~sumfem.models.Conduction`; - the current density of the cell at the positive terminal, with the negative terminal grounded, see :meth:`conditions`; - - with flow, :class:`~sumfem.models.Darcy` flow in the electrodes from - their inlets to their outlets, :eq:`full-cell-flow`, and tanks, - :meth:`update_tanks`. + - with flow, the flow of the electrolytes through the electrodes from + their inlets to their outlets, :eq:`full-cell-flow`, an + :class:`~sumfem.electrochemistry.ElectrodeFlow` such as + :class:`~sumfem.electrochemistry.DarcyFlow`, see :attr:`flow`, and + tanks, :meth:`update_tanks`. Args: domain: The domain of the cell, with the markers of :meth:`mesh_data`. @@ -222,8 +88,7 @@ class FullCell: 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. + ValueError: If only one side transports water. """ _domain: Domain @@ -241,21 +106,15 @@ class FullCell: _water: Field | None """Concentration of the water, over the electrodes and the membrane.""" - _flow: FlowSpec | None - """The flow of the electrolytes, if any.""" + _spec: FlowSpec | None + """The specification of the flow of the electrolytes, if any.""" - _pressure: Field | None - """Pressure over the electrodes and the membrane, with flow.""" + _flow: ElectrodeFlow | None + """The flow of the electrolytes, if any.""" _inlet_velocity: ufl.core.expr.Expr | None """Superficial velocity at the inlets, with flow.""" - _outlet_pressure: ufl.core.expr.Expr | None - """Gauge pressure at the outlets, with flow.""" - - _darcy: types.MappingProxyType[str, Darcy] - """The flow in each electrode, by side, empty without flow.""" - _half_cells: types.MappingProxyType[str, RedoxHalfCell] """The half-cells, by side.""" @@ -280,24 +139,32 @@ def __init__( negative: HalfCellSpec, positive: HalfCellSpec, membrane: MembraneSpec = MembraneSpec(), - collector_conductivity: PROPERTY = 1000.0, + collector_conductivity: ExpressionOrDefault = 1000.0, flow: FlowSpec | None = None, degree: int = 1, name: str = "cell", ): specs = dict(zip(SIDES, (negative, positive), strict=True)) - _check_sides(specs) + if (WATER in negative.reservoir) != (WATER in positive.reservoir): + raise ValueError( + f"Water, {WATER!r}, moves through the membrane, so both sides or " + "neither must transport it." + ) self._domain = domain self._name = name - self._flow = flow + self._spec = flow tdim = domain.tdim markers = [ - *(f"{side.capitalize()}{part}" for side in SIDES for part in _PARTS), + *( + f"{side.capitalize()}{part}" + for side in SIDES + for part in ("Collector", "Electrode", "Solid") + ), "Membrane", "Electrolyte", ] self._subdomains = types.MappingProxyType( - {m: Subdomain(f"{name}_{_snake(m)}", domain, tdim, m) for m in markers} + {m: Subdomain(f"{name}_{m}", domain, tdim, m) for m in markers} ) element = ("Lagrange", degree) electrolyte = self._subdomains["Electrolyte"] @@ -310,21 +177,14 @@ def __init__( electrodes = { side: self._create_electrode(side, spec) for side, spec in specs.items() } - viscosities = { - side: expression_or_parameter( - f"{side}_electrolyte.viscosity", - spec.viscosity, - self._subdomains[f"{side.capitalize()}Electrode"], - unit=VISCOSITY, - doc="Dynamic viscosity", - ) - for side, spec in specs.items() + properties = { + side: self._create_properties(side, spec) for side, spec in specs.items() } - self._create_flow(electrodes, viscosities, element) + self._create_flow(electrodes, properties, degree) self._half_cells = types.MappingProxyType( { side: self._create_half_cell( - side, spec, element, electrodes[side], viscosities[side] + side, spec, element, electrodes[side], *properties[side] ) for side, spec in specs.items() } @@ -334,7 +194,7 @@ def __init__( self._phi_l, ("Membrane",), water=self._water, - pressure=self._pressure, + pressure=None if self._flow is None else self._flow.pressure, fixed_charge=membrane.fixed_charge, proton_diffusivity=membrane.proton_diffusivity, water_diffusivity=membrane.water_diffusivity, @@ -381,74 +241,98 @@ def _create_electrode(self, side: str, spec: HalfCellSpec) -> PorousElectrode: kozeny_carman_constant=spec.kozeny_carman_constant, ) + def _create_properties( + self, side: str, spec: HalfCellSpec + ) -> tuple[ufl.core.expr.Expr, ufl.core.expr.Expr]: + """Create the density and the viscosity of a side's electrolyte. + + They are parameters of the side's electrolyte, which its half-cell + and the flow share. + + Args: + side: ``"negative"`` or ``"positive"``. + spec: The side. + + Returns: + The density and the viscosity. + """ + domain = self._subdomains[f"{side.capitalize()}Electrode"] + density = expression_or_parameter( + f"{side}_electrolyte.density", + spec.density, + domain, + unit=DENSITY, + doc="Density", + ) + viscosity = expression_or_parameter( + f"{side}_electrolyte.viscosity", + spec.viscosity, + domain, + unit=VISCOSITY, + doc="Dynamic viscosity", + ) + return density, viscosity + def _create_flow( self, electrodes: typing.Mapping[str, PorousElectrode], - viscosities: typing.Mapping[str, ufl.core.expr.Expr], - element: tuple[str, int], + properties: typing.Mapping[str, tuple[ufl.core.expr.Expr, ufl.core.expr.Expr]], + degree: int, ) -> None: - """Create the pressure, its parameters and the flow in the electrodes. - - The pressure's scale is the pressure drop of a uniform flow at the - inlet velocity through the negative electrode over the domain's - length scale, if the defaults give it, else one pascal. + """Create the inlet velocity and the flow of the flow's specification. Args: electrodes: The electrodes, by side. - viscosities: The viscosities of the electrolytes, by side. - element: The element of the pressure. + properties: The density and the viscosity of each side's + electrolyte. + degree: The degree of the Lagrange elements of the unknowns. """ - flow, name, domain = self._flow, self._name, self._domain - if flow is None: - self._pressure = None - self._inlet_velocity = self._outlet_pressure = None - self._darcy = types.MappingProxyType({}) + spec, name, domain = self._spec, self._name, self._domain + if spec is None: + self._flow = self._inlet_velocity = None return self._inlet_velocity = expression_or_parameter( f"{name}.inlet_velocity", - flow.inlet_velocity, + spec.inlet_velocity, domain, unit=syu.m / syu.second, doc="Superficial velocity at the inlets", ) - self._outlet_pressure = expression_or_parameter( + outlet_pressure = expression_or_parameter( f"{name}.outlet_pressure", - flow.outlet_pressure, + spec.outlet_pressure, domain, unit=syu.pascal, doc="Gauge pressure at the outlets", ) - negative = electrodes["negative"] - try: - scale = float( - value_at_defaults( - viscosities["negative"] - * self._inlet_velocity - * domain.length_scale - / negative.permeability() - ) - ) - except ComponentError: - scale = 1.0 - self._pressure = Field( - "p", - self._subdomains["Electrolyte"], - element, - unit=syu.pascal, - scale=scale, - ) - self._darcy = types.MappingProxyType( - { - side: Darcy( - f"{name}_{side}_darcy", - self._pressure, - electrode.permeability(), - viscosities[side], - markers=electrode.markers, - ) - for side, electrode in electrodes.items() - } + setup = FlowSetup( + name=name, + domain=domain, + electrolyte=self._subdomains["Electrolyte"], + electrodes=types.MappingProxyType(dict(electrodes)), + subdomains=types.MappingProxyType( + { + side: self._subdomains[f"{side.capitalize()}Electrode"] + for side in SIDES + } + ), + viscosities=types.MappingProxyType( + {side: viscosity for side, (_, viscosity) in properties.items()} + ), + densities=types.MappingProxyType( + {side: density for side, (density, _) in properties.items()} + ), + inlets=types.MappingProxyType( + {side: f"{side.capitalize()}Inlet" for side in SIDES} + ), + outlets=types.MappingProxyType( + {side: f"{side.capitalize()}Outlet" for side in SIDES} + ), + inlet_velocity=self._inlet_velocity, + outlet_pressure=outlet_pressure, + degree=degree, ) + self._flow = spec.model(setup) def _create_half_cell( self, @@ -456,6 +340,7 @@ def _create_half_cell( spec: HalfCellSpec, element: tuple[str, int], electrode: PorousElectrode, + density: ufl.core.expr.Expr, viscosity: ufl.core.expr.Expr, ) -> RedoxHalfCell: """Create the half-cell of a side, with its unknowns on their subdomains. @@ -465,6 +350,7 @@ def _create_half_cell( spec: The side. element: The element of the unknowns. electrode: The electrode of the side. + density: The density of the side's electrolyte. viscosity: The viscosity of the side's electrolyte. Returns: @@ -492,19 +378,19 @@ def _create_half_cell( } if self._water is not None: fields[WATER] = self._water - darcy = self._darcy.get(side) + flow = self._flow return RedoxHalfCell( electrode_domain, spec.ions, spec.reservoir, spec.couple, electrode, - velocity=None if darcy is None else darcy.superficial_velocity, + velocity=None if flow is None else flow.velocity(side), fiber_distance=spec.fiber_distance, collector=None, membrane=None, - outlet=() if darcy is None else (f"{prefix}Outlet",), - density=spec.density, + outlet=() if flow is None else (f"{prefix}Outlet",), + density=density, viscosity=viscosity, name=side, phi_s=phi_s, @@ -568,24 +454,19 @@ def _create_models(self) -> list[Residual]: return models def _create_inflow(self) -> list[Residual]: - """Create the flow in the electrodes and the inflow of the inlets. + """Create the models of the flow and the inflow of the ions. Returns: - The models of :eq:`full-cell-flow` but the outlets': the Darcy - flow and the inflow of each electrode, and the inflow of each of - its ions, empty without flow. + The models of the flow, and the inflow of each ion of each + electrode, :eq:`full-cell-flow`, empty without flow. """ - if self._pressure is None: + if self._flow is None: return [] assert self._inlet_velocity is not None name, velocity = self._name, self._inlet_velocity - models: list[Residual] = [] + models: list[Residual] = list(self._flow.models) for side, cell in self._half_cells.items(): inlet = (f"{side.capitalize()}Inlet",) - models.append(self._darcy[side]) - models.append( - BoundaryFlux(f"{name}_{side}_inflow", self._pressure, velocity, inlet) - ) for ion, c in cell.concentrations.items(): models.append( BoundaryFlux( @@ -598,13 +479,13 @@ def _create_inflow(self) -> list[Residual]: return models def _create_functionals(self) -> types.MappingProxyType[str, ufl.Form]: - """Create the integrals of the voltage and of the tanks. + """Create the integrals of the voltage, the amounts and the tanks. Returns: - The areas and the potential integrals of the terminals, and with - flow, for each side, the volume of electrolyte in the pores of - its electrode, the area of its inlet and the outflow of each ion, - by name. + The areas and the potential integrals of the terminals, the + amount of each ion in each electrode, the volume of electrolyte + in its pores, and with flow, the area of its inlet and the + outflow of each ion, by name. """ functionals: dict[str, ufl.Form] = {} for side, cell in self._half_cells.items(): @@ -612,14 +493,17 @@ def _create_functionals(self) -> types.MappingProxyType[str, ufl.Form]: ds = cell.phi_s.domain.ds(f"{prefix}Terminal") functionals[f"{side}_area"] = 1.0 * ds functionals[f"{side}_potential"] = cell.phi_s * ds - if side not in self._darcy: - continue electrode = self._subdomains[f"{prefix}Electrode"] - functionals[f"{side}_pore_volume"] = ( - cell.electrode.porosity * electrode.dx() - ) + porosity = cell.electrode.porosity + functionals[f"{side}_pore_volume"] = porosity * electrode.dx() + for ion, c in cell.concentrations.items(): + functionals[f"{side}_amount_{ion}"] = ( + porosity * c * c.domain.dx(f"{prefix}Electrode") + ) + if self._flow is None: + continue functionals[f"{side}_inlet_area"] = 1.0 * electrode.ds(f"{prefix}Inlet") - u = self._darcy[side].superficial_velocity + u = self._flow.velocity(side) for ion, c in cell.concentrations.items(): n = ufl.FacetNormal(c.domain) functionals[f"{side}_outflow_{ion}"] = ( @@ -672,9 +556,9 @@ def model(self) -> Residual: the conduction in the collectors, ``"__collector"``, the current at the positive terminal, ``"_terminal_current"``, the proton fluxes at the faces, ``"__protons"``, and with - flow, the flow in each electrode, ``"__darcy"``, its - inflow, ``"__inflow"``, and that of each of its ions, - ``"__inflow_"``. + flow, the models of the flow, e.g. those of :class:`DarcyFlow + `, and the inflow of each ion of + each electrode, ``"__inflow_"``. """ return self._model @@ -699,14 +583,15 @@ def water(self) -> Field | None: return self._water @property - def pressure(self) -> Field | None: - """Pressure over the electrodes and the membrane, with flow.""" - return self._pressure + def flow(self) -> ElectrodeFlow | None: + """The flow of the electrolytes, if any.""" + return self._flow @property - def darcy(self) -> typing.Mapping[str, Darcy]: - """The flow in each electrode, by side, empty without flow.""" - return self._darcy + def pressure(self) -> Field | None: + """The pressure of the flow over the electrodes and the membrane, + ``None`` without flow.""" + return None if self._flow is None else self._flow.pressure @property def current(self) -> Parameter: @@ -721,20 +606,14 @@ def conditions(self) -> list[DirichletBC]: """The Dirichlet conditions of the cell. These are the negative terminal at zero potential, which fixes the - potentials, and with flow, the outlet pressure. + potentials, and with flow, those of the flow, e.g. the outlet + pressure of :class:`DarcyFlow `. """ conditions = [ DirichletBC(self._half_cells["negative"].phi_s, "NegativeTerminal", 0.0) ] - if self._pressure is not None: - assert self._outlet_pressure is not None - conditions.append( - DirichletBC( - self._pressure, - ("NegativeOutlet", "PositiveOutlet"), - self._outlet_pressure, - ) - ) + if self._flow is not None: + conditions.extend(self._flow.conditions()) return conditions def initialize(self, discretization: Discretization) -> None: @@ -744,8 +623,9 @@ def initialize(self, discretization: Discretization) -> None: defaults of their parameters ``".c0_"``, the electrolyte potential at zero, the solid of each side at the equilibrium potential of its couple at those concentrations, by Nernst's - equation, so that no current flows, and the pressure at that of the - outlets. + equation, so that no current flows, and the flow as it initializes + itself, see :meth:`ElectrodeFlow.initialize + `. Raises: ValueError: If the sides start with different water @@ -754,41 +634,20 @@ def initialize(self, discretization: Discretization) -> None: 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) + values[c] = cell.reservoir_concentration(ion) + values[cell.phi_s] = cell.equilibrium_potential(discretization) if self._water is not None: - start = {_initial(c, WATER) for c in self._half_cells.values()} + start = { + c.reservoir_concentration(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)}." ) - if self._pressure is not None: - assert isinstance(self._outlet_pressure, Parameter) - values[self._pressure] = float(discretization.value(self._outlet_pressure)) 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)) + if self._flow is not None: + self._flow.initialize(discretization) def voltage(self, discretization: Discretization) -> float: """The cell voltage in V: the mean potential of the positive terminal @@ -812,10 +671,55 @@ def tank_volume(self, discretization: Discretization, side: str) -> float: Raises: ValueError: For a cell without flow, or without tanks. """ - if self._flow is None or self._flow.tank_ratio is None: + if self._spec is None or self._spec.tank_ratio is None: raise ValueError(f"The cell {self._name!r} has no tanks.") pores = discretization.assemble_scalar(self._functionals[f"{side}_pore_volume"]) - return self._flow.tank_ratio * pores + return self._spec.tank_ratio * pores + + def amount(self, discretization: Discretization, side: str, ion: str) -> float: + """The amount of an ion in a side's electrode and its tank. + + The amounts are per unit of the dimensions the mesh lacks, e.g. per + metre of depth for a 2D cell. + + Args: + discretization: The discretization of the cell. + side: ``"negative"`` or ``"positive"``. + ion: The name of an ion of the side with a field. The water of + the membrane is in neither side's amount. + + Returns: + The amount in mol, that of the tank included with tanks. + """ + cell = self._half_cells[side] + amount = discretization.assemble_scalar( + self._functionals[f"{side}_amount_{ion}"] + ) + if self._spec is not None and self._spec.tank_ratio is not None: + tank = float(discretization.value(cell.inlet[ion])) + amount += self.tank_volume(discretization, side) * tank + return amount + + def state_of_charge(self, discretization: Discretization, side: str) -> float: + """The state of charge of a side's couple, its tank included. + + Args: + discretization: The discretization of the cell. + side: ``"negative"`` or ``"positive"``. + + Returns: + The share of the ion that charging makes: the reduced ion of a + couple oxidized on discharge, the oxidized ion of one reduced on + discharge. + """ + couple = self._half_cells[side].couple + charged, other = ( + (couple.reduced, couple.oxidized) + if couple.oxidation_on_discharge + else (couple.oxidized, couple.reduced) + ) + made = self.amount(discretization, side, charged) + return made / (made + self.amount(discretization, side, other)) def update_tanks(self, discretization: Discretization, dt: float) -> None: r"""Take a step of the tanks after a time step of the cell. @@ -845,7 +749,7 @@ def update_tanks(self, discretization: Discretization, dt: float) -> None: discretization: The discretization of the cell, after the step. dt: The size of the step, in s. """ - if self._flow is None or self._flow.tank_ratio is None: + if self._spec is None or self._spec.tank_ratio is None: return assemble, value = discretization.assemble_scalar, discretization.value assert isinstance(self._inlet_velocity, Parameter) @@ -860,54 +764,3 @@ def update_tanks(self, discretization: Discretization, dt: float) -> None: inflow = velocity * area * tank updates[cell.inlet[ion]] = tank + dt * (outflow - inflow) / volume discretization.set(updates) - - -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 03df185..7d3f94b 100644 --- a/src/sumfem/electrochemistry/half_cell.py +++ b/src/sumfem/electrochemistry/half_cell.py @@ -31,16 +31,18 @@ Residual, ResidualSum, ) +from sumfem.symbolic.si import CONCENTRATION from .charge import ElectrolyteConduction from .couple import RedoxCouple -from .electrolyte import CONCENTRATION, Electrolyte, Ion +from .electrode import PorousElectrode +from .electrolyte import Electrolyte, Ion from .kinetics import ( ButlerVolmer, MassTransferButlerVolmer, PistonVelocityButlerVolmer, ) -from .porous import PorousElectrode, Region +from .porous import Region from .transport import TransientNernstPlanck __all__ = ["RedoxHalfCell"] @@ -734,6 +736,43 @@ def inlet_conditions(self, markers: str | tuple[str, ...]) -> list[DirichletBC]: for ion, c in self._concentrations.items() ] + def reservoir_concentration(self, ion: str) -> float: + """The reservoir concentration of a transported ion, as given. + + Args: + ion: The name of the ion. + + Returns: + The default of its parameter ``".c0_"``, in mol/m^3. + """ + default = self._inlet[ion].default_si + assert default is not None + return float(default) + + def equilibrium_potential(self, discretization: Discretization) -> float: + """Nernst's potential of the couple at its reservoir concentrations. + + That is ``E0 + R T / (n F) ln(c_ox / c_red)``, see + :meth:`~sumfem.electrochemistry.ButlerVolmer.equilibrium_potential`, + with the reservoir concentrations as given, + :meth:`reservoir_concentration`. + + Args: + discretization: The discretization, for the values of the + formal potential and the constants. + + Returns: + The potential in V. + """ + value, couple, e = discretization.value, self._couple, self._electrolyte + # The half-cell's electrolyte has the temperature as a parameter + assert isinstance(e.T, Parameter) + thermal = value(e.R) * value(e.T) / (couple.electrons * value(e.F)) + ratio = self.reservoir_concentration( + couple.oxidized + ) / self.reservoir_concentration(couple.reduced) + return float(value(self._formal_potential) + thermal * np.log(ratio)) + def initialize(self, discretization: Discretization) -> None: """Set a state at rest, as RFBniCS starts. diff --git a/src/sumfem/electrochemistry/kinetics.py b/src/sumfem/electrochemistry/kinetics.py index 762556d..cdbffab 100644 --- a/src/sumfem/electrochemistry/kinetics.py +++ b/src/sumfem/electrochemistry/kinetics.py @@ -9,12 +9,11 @@ import typing -import numpy.typing as npt import sympy.physics.units as syu import ufl from sumfem.models import RateReaction, butler_volmer -from sumfem.symbolic import Field, as_expr +from sumfem.symbolic import ExpressionOrValue, Field, as_expr from .electrolyte import Electrolyte, Ion @@ -26,9 +25,6 @@ "PistonVelocityButlerVolmer", ] -INPUT = ufl.core.expr.Expr | npt.ArrayLike | None -"""A model input: an expression, a value, or ``None`` for a parameter.""" - class ElectrodeReaction(RateReaction): """An electrode reaction ``reactants -> products + n e-`` in a porous electrode. @@ -84,7 +80,7 @@ def __init__( electrolyte: Electrolyte, reactants: typing.Mapping[ufl.core.expr.Expr, int] | None = None, products: typing.Mapping[ufl.core.expr.Expr, int] | None = None, - specific_area: INPUT = None, + specific_area: ExpressionOrValue = None, electrons: int = 1, markers: tuple[str, ...] = (), ): @@ -189,13 +185,13 @@ def __init__( solid_potential: Field, electrolyte_potential: Field, electrolyte: Electrolyte, - exchange_current: INPUT = None, - equilibrium_potential: INPUT = None, - anodic_transfer: INPUT = None, - cathodic_transfer: INPUT = None, + exchange_current: ExpressionOrValue = None, + equilibrium_potential: ExpressionOrValue = None, + anodic_transfer: ExpressionOrValue = None, + cathodic_transfer: ExpressionOrValue = None, reduced: ufl.core.expr.Expr | None = None, oxidized: ufl.core.expr.Expr | None = None, - specific_area: INPUT = None, + specific_area: ExpressionOrValue = None, electrons: int = 1, markers: tuple[str, ...] = (), ): @@ -345,11 +341,11 @@ def __init__( electrolyte: Electrolyte, reduced: str, oxidized: str, - specific_area: INPUT = None, - rate_constant: INPUT = None, - formal_potential: INPUT = None, - anodic_transfer: INPUT = None, - cathodic_transfer: INPUT = None, + specific_area: ExpressionOrValue = None, + rate_constant: ExpressionOrValue = None, + formal_potential: ExpressionOrValue = None, + anodic_transfer: ExpressionOrValue = None, + cathodic_transfer: ExpressionOrValue = None, electrons: int = 1, markers: tuple[str, ...] = (), reactants: typing.Mapping[str, int] | None = None, @@ -581,18 +577,18 @@ def __init__( reduced: str, oxidized: str, velocity: ufl.core.expr.Expr, - fiber_diameter: INPUT = None, - density: INPUT = None, - viscosity: INPUT = None, + fiber_diameter: ExpressionOrValue = None, + density: ExpressionOrValue = None, + viscosity: ExpressionOrValue = None, sherwood_factor: float = 7.0, reynolds_exponent: float = 0.4, speed_floor: float = 1e-16, transfer_floor: float = 1e-14, - specific_area: INPUT = None, - rate_constant: INPUT = None, - formal_potential: INPUT = None, - anodic_transfer: INPUT = None, - cathodic_transfer: INPUT = None, + specific_area: ExpressionOrValue = None, + rate_constant: ExpressionOrValue = None, + formal_potential: ExpressionOrValue = None, + anodic_transfer: ExpressionOrValue = None, + cathodic_transfer: ExpressionOrValue = None, electrons: int = 1, markers: tuple[str, ...] = (), reactants: typing.Mapping[str, int] | None = None, @@ -771,7 +767,7 @@ def __init__( electrolyte: Electrolyte, reduced: str, oxidized: str, - fiber_distance: INPUT = None, + fiber_distance: ExpressionOrValue = None, transfer_floor: float = 1e-14, **kwargs: typing.Any, ): diff --git a/src/sumfem/electrochemistry/membrane.py b/src/sumfem/electrochemistry/membrane.py index 105759e..4325379 100644 --- a/src/sumfem/electrochemistry/membrane.py +++ b/src/sumfem/electrochemistry/membrane.py @@ -4,8 +4,6 @@ import typing -import numpy.typing as npt -import sympy import sympy.physics.units as syu import ufl @@ -14,30 +12,20 @@ from sumfem.symbolic import ( DirichletBC, Domain, + ExpressionOrDefault, + ExpressionOrValue, Field, Residual, ResidualSum, as_expr, expression_or_parameter, ) +from sumfem.symbolic.si import CONCENTRATION, DIFFUSIVITY, PERMEABILITY, VISCOSITY -from .electrolyte import CONCENTRATION, DIFFUSIVITY, Electrolyte +from .electrolyte import Electrolyte __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.""" - -PROPERTY: typing.TypeAlias = ufl.core.expr.Expr | float | sympy.Expr | None -"""A property of the membrane: an expression, used as given, or the default -of the parameter the membrane creates, ``None`` for none.""" - -PERMEABILITY = syu.m**2 -"""Unit of permeabilities.""" - -VISCOSITY = syu.pascal * syu.second -"""Unit of dynamic viscosities.""" - class SchloeglFlow(Residual): r"""Flow of water through a charged membrane, by Schlögl's equation. @@ -109,10 +97,10 @@ def __init__( name: str, pressure: Field, potential: ufl.core.expr.Expr, - charge_density: INPUT = None, - electrokinetic_permeability: INPUT = None, - hydraulic_permeability: INPUT = None, - viscosity: INPUT = None, + charge_density: ExpressionOrValue = None, + electrokinetic_permeability: ExpressionOrValue = None, + hydraulic_permeability: ExpressionOrValue = None, + viscosity: ExpressionOrValue = None, markers: tuple[str, ...] = (), ): super().__init__(name) @@ -266,12 +254,12 @@ def __init__( markers: tuple[str, ...], water: Field | None = None, pressure: Field | None = None, - fixed_charge: PROPERTY = None, - proton_diffusivity: PROPERTY = None, - water_diffusivity: PROPERTY = None, - electrokinetic_permeability: PROPERTY = None, - hydraulic_permeability: PROPERTY = None, - viscosity: PROPERTY = None, + fixed_charge: ExpressionOrDefault = None, + proton_diffusivity: ExpressionOrDefault = None, + water_diffusivity: ExpressionOrDefault = None, + electrokinetic_permeability: ExpressionOrDefault = None, + hydraulic_permeability: ExpressionOrDefault = None, + viscosity: ExpressionOrDefault = None, temperature: ufl.core.expr.Expr | float = 298.15, ): domain = potential.domain diff --git a/src/sumfem/electrochemistry/porous.py b/src/sumfem/electrochemistry/porous.py index ef29544..9ee46c4 100644 --- a/src/sumfem/electrochemistry/porous.py +++ b/src/sumfem/electrochemistry/porous.py @@ -1,17 +1,14 @@ -"""Porous electrodes and the regions of a cell filled with electrolyte.""" +"""The regions of a cell filled with electrolyte, free or in pores.""" from __future__ import annotations import dataclasses -import sympy -import sympy.physics.units as syu import ufl -from sumfem.models import bruggeman, ergun_coefficient, kozeny_carman -from sumfem.symbolic import Domain, Parameter +from sumfem.models import bruggeman -__all__ = ["PorousElectrode", "Region"] +__all__ = ["Region"] @dataclasses.dataclass(frozen=True) @@ -42,146 +39,3 @@ def transport(self) -> ufl.core.expr.Expr | float: It is one in free electrolyte, with porosity one. """ return bruggeman(1.0, self.porosity, self.bruggeman_exponent) - - -class PorousElectrode: - """A porous electrode of carbon fibers, filled with electrolyte. - - The electrode creates a parameter for each argument below from the - porosity on, named ``"."``, such as - ``".porosity"``, with the argument as its default value. Models - that share an electrode share these parameters. A parameter without a - default value must be set on the discretization before it is used. - - Args: - name: Prefix of the parameter names. - domain: The domain of the parameters. - markers: Cell markers of the electrode. Empty for the whole domain. - porosity: Default porosity, if any. Every default with a unit is in - SI units, or a quantity with units, e.g. ``0.5 * syu.micrometer``. - conductivity: Default conductivity of the fiber material in S/m, if - any. - specific_area: Default specific area in 1/m, if any. - fiber_diameter: Default fiber diameter in m, if any. - kozeny_carman_constant: Default Kozeny-Carman constant, if any. - pore_length: Default characteristic pore length of Ergun's drag in - m, if any. - """ - - _markers: tuple[str, ...] - """Cell markers of the electrode.""" - - _porosity: Parameter - """Porosity.""" - - _conductivity: Parameter - """Conductivity of the fiber material.""" - - _specific_area: Parameter - """Specific area.""" - - _fiber_diameter: Parameter - """Fiber diameter.""" - - _kozeny_carman_constant: Parameter - """Kozeny-Carman constant.""" - - _pore_length: Parameter - """Characteristic pore length of Ergun's drag.""" - - def __init__( - self, - name: str, - domain: Domain, - markers: tuple[str, ...] = (), - porosity: float | None = None, - conductivity: float | sympy.Expr | None = None, - specific_area: float | sympy.Expr | None = None, - fiber_diameter: float | sympy.Expr | None = None, - kozeny_carman_constant: float | None = None, - pore_length: float | sympy.Expr | None = None, - ): - self._markers = tuple(markers) - self._porosity = Parameter( - f"{name}.porosity", domain, default=porosity, doc="Porosity" - ) - self._conductivity = Parameter( - f"{name}.conductivity", - domain, - default=conductivity, - unit=syu.siemens / syu.m, - doc="Conductivity of the fiber material", - ) - self._specific_area = Parameter( - f"{name}.specific_area", - domain, - default=specific_area, - unit=1 / syu.m, - doc="Specific area", - ) - self._fiber_diameter = Parameter( - f"{name}.fiber_diameter", - domain, - default=fiber_diameter, - unit=syu.m, - doc="Fiber diameter", - ) - self._kozeny_carman_constant = Parameter( - f"{name}.kozeny_carman_constant", - domain, - default=kozeny_carman_constant, - doc="Kozeny-Carman constant", - ) - self._pore_length = Parameter( - f"{name}.pore_length", - domain, - default=pore_length, - unit=syu.m, - doc="Characteristic pore length of Ergun's drag", - ) - - @property - def markers(self) -> tuple[str, ...]: - """Cell markers of the electrode.""" - return self._markers - - @property - def region(self) -> Region: - """The electrolyte-filled pores, as a region.""" - return Region(self._markers, self._porosity) - - @property - def porosity(self) -> Parameter: - """Porosity.""" - return self._porosity - - @property - def specific_area(self) -> Parameter: - """Specific area.""" - return self._specific_area - - @property - def fiber_diameter(self) -> Parameter: - """Fiber diameter.""" - return self._fiber_diameter - - def effective_conductivity(self) -> ufl.core.expr.Expr: - r"""Bruggeman's conductivity of the solid. - - That is :math:`\sigma_\text{eff} = (1 - \varepsilon)^{3/2} \sigma`, with - the conductivity :math:`\sigma` of the fibres and the porosity - :math:`\varepsilon`. - """ - return bruggeman(self._conductivity, 1.0 - self._porosity) - - def permeability(self) -> ufl.core.expr.Expr: - """Kozeny-Carman permeability of the fibers, see - :func:`~sumfem.models.kozeny_carman`.""" - return kozeny_carman( - self._fiber_diameter, self._porosity, self._kozeny_carman_constant - ) - - def forchheimer(self) -> ufl.core.expr.Expr: - """Ergun's inertial drag coefficient, an estimate for fibres, see - :func:`~sumfem.models.ergun_coefficient`.""" - return ergun_coefficient(self._porosity, self._pore_length) diff --git a/src/sumfem/electrochemistry/specs.py b/src/sumfem/electrochemistry/specs.py new file mode 100644 index 0000000..cd475a3 --- /dev/null +++ b/src/sumfem/electrochemistry/specs.py @@ -0,0 +1,163 @@ +"""What a full cell is made of: its sides, its membrane and its flow. + +The specifications of a :class:`~sumfem.electrochemistry.FullCell`: a +:class:`HalfCellSpec` for each side, a :class:`MembraneSpec` and a +:class:`FlowSpec`. +""" + +from __future__ import annotations + +import dataclasses +import typing + +import sympy + +from sumfem.symbolic import ExpressionOrDefault + +from .couple import RedoxCouple +from .electrode import DarcyFlow, ElectrodeFlow, FlowSetup +from .electrolyte import Ion + +__all__ = ["FlowSpec", "HalfCellSpec", "MembraneSpec"] + +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`, and with flow, that of its tank. The + protons, ``"H"``, must be one of them. + + Raises: + ValueError: If the side does not transport protons. + 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 them, limited by the flow's mass + transfer with flow. + fiber_diameter: The diameter of the electrode's fibres, in m or with + units, for the permeability with flow, and the mass transfer + without a fibre distance. + kozeny_carman_constant: The Kozeny-Carman constant of the + permeability with flow, see + :func:`~sumfem.models.kozeny_carman`. + density: The density of the electrolyte, in kg/m^3 or with units. + viscosity: The dynamic viscosity of the electrolyte, in Pa s or with + units, which the flow and the mass transfer share. + """ + + 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 + fiber_diameter: float | sympy.Expr | None = None + kozeny_carman_constant: float | None = None + density: float | sympy.Expr = 1399.8 + viscosity: float | sympy.Expr = 6.963e-3 + + def __post_init__(self) -> None: + if PROTONS not in self.reservoir: + raise ValueError( + f"A side of a full cell must transport protons, {PROTONS!r}, which " + "carry the current through the membrane." + ) + + +@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`. The + last three are those of the flow of water, with a :class:`FlowSpec`. + + 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. + electrokinetic_permeability: The electrokinetic permeability, in + m^2. + hydraulic_permeability: The hydraulic permeability, in m^2. + viscosity: The viscosity of the water in the membrane, in Pa s. + """ + + fixed_charge: ExpressionOrDefault = None + proton_diffusivity: ExpressionOrDefault = None + water_diffusivity: ExpressionOrDefault = None + electrokinetic_permeability: ExpressionOrDefault = None + hydraulic_permeability: ExpressionOrDefault = None + viscosity: ExpressionOrDefault = None + + +@dataclasses.dataclass(frozen=True) +class FlowSpec: + r"""The flow of the electrolytes through the electrodes, from their tanks. + + The flow of the electrolytes is a :class:`ElectrodeFlow`, by default + :class:`DarcyFlow`, which gives the superficial velocity + :math:`\mathbf{u}` of each electrode. The electrolyte of each side + enters its electrode at the inlet with the superficial velocity + :math:`v_\text{in}` and the concentrations :math:`c_k^\text{in}` of its + tank, and leaves at the outlet, held at the pressure + :math:`p_\text{out}`, with the flow, + + .. math:: + :label: full-cell-flow + + \begin{aligned} + -\mathbf{u} \cdot \mathbf{n} &= v_\text{in} + && \text{at the inlet}, \\ + -\mathbf{N}_k \cdot \mathbf{n} &= v_\text{in} c_k^\text{in} + && \text{at the inlet}, \\ + p &= p_\text{out} + && \text{at the outlet}, \\ + (\mathbf{N}_k - \mathbf{u} c_k) \cdot \mathbf{n} &= 0 + && \text{at the outlet}, + \end{aligned} + + with the superficial velocity :math:`\mathbf{u}` of the electrode and + the Nernst-Planck flux :math:`\mathbf{N}_k` of each ion. The inlet holds + the flux of the ions, the flow's, rather than their concentration as + Shah et al. do, so that the cell takes in exactly what its tank gives. + + Attributes: + inlet_velocity: The superficial velocity at the inlets of both + electrodes: an expression, or the default of the parameter + ``".inlet_velocity"``, in m/s or with units. + outlet_pressure: The gauge pressure at the outlets, given like the + inlet velocity, in Pa. + tank_ratio: The volume of each tank over the volume of electrolyte in + the pores of its electrode, ``None`` for tanks so large that the + electrolyte at the inlets does not change. + model: The flow, created from the :class:`FlowSetup` the cell gives + it: an :class:`ElectrodeFlow` subclass, or any callable of the + setup, e.g. a :func:`functools.partial` of one with options of + its own. + """ + + inlet_velocity: ExpressionOrDefault = None + outlet_pressure: ExpressionOrDefault = 0.0 + tank_ratio: float | None = None + model: typing.Callable[[FlowSetup], ElectrodeFlow] = DarcyFlow From 96c06df5e2ca327e1271a32538f5f8bf8192da75 Mon Sep 17 00:00:00 2001 From: Jorgen Schartum Dokken Date: Thu, 12 Nov 2026 12:38:28 +0000 Subject: [PATCH 6/7] Refactor electrochemistry into cells (what each bigger component of our problem describes), materials (electrolyte, electrode and what couples them) and models (the various UFL forms) --- docs/api/electrochemistry.md | 62 ++++-- src/sumfem/electrochemistry/__init__.py | 60 +++--- src/sumfem/electrochemistry/cells/__init__.py | 15 ++ .../electrochemistry/{ => cells}/full_cell.py | 4 +- .../electrochemistry/{ => cells}/half_cell.py | 105 ++++++---- .../electrochemistry/{ => cells}/specs.py | 34 +-- .../electrochemistry/materials/__init__.py | 19 ++ .../{ => materials}/couple.py | 16 +- .../electrochemistry/materials/electrode.py | 194 ++++++++++++++++++ .../{ => materials}/electrolyte.py | 7 +- .../electrochemistry/models/__init__.py | 46 +++++ .../electrochemistry/{ => models}/charge.py | 7 +- .../electrode_flow.py} | 175 +--------------- .../electrochemistry/{ => models}/kinetics.py | 4 +- .../electrochemistry/{ => models}/membrane.py | 2 +- .../{ => models}/transport.py | 7 +- src/sumfem/electrochemistry/porous.py | 41 ---- tests/electrochemistry/test_full_cell.py | 76 ++++++- 18 files changed, 542 insertions(+), 332 deletions(-) create mode 100644 src/sumfem/electrochemistry/cells/__init__.py rename src/sumfem/electrochemistry/{ => cells}/full_cell.py (99%) rename src/sumfem/electrochemistry/{ => cells}/half_cell.py (92%) rename src/sumfem/electrochemistry/{ => cells}/specs.py (87%) create mode 100644 src/sumfem/electrochemistry/materials/__init__.py rename src/sumfem/electrochemistry/{ => materials}/couple.py (88%) create mode 100644 src/sumfem/electrochemistry/materials/electrode.py rename src/sumfem/electrochemistry/{ => materials}/electrolyte.py (98%) create mode 100644 src/sumfem/electrochemistry/models/__init__.py rename src/sumfem/electrochemistry/{ => models}/charge.py (89%) rename src/sumfem/electrochemistry/{electrode.py => models/electrode_flow.py} (61%) rename src/sumfem/electrochemistry/{ => models}/kinetics.py (99%) rename src/sumfem/electrochemistry/{ => models}/membrane.py (99%) rename src/sumfem/electrochemistry/{ => models}/transport.py (96%) delete mode 100644 src/sumfem/electrochemistry/porous.py diff --git a/docs/api/electrochemistry.md b/docs/api/electrochemistry.md index e8e7cdb..53d8f53 100644 --- a/docs/api/electrochemistry.md +++ b/docs/api/electrochemistry.md @@ -4,7 +4,15 @@ .. automodule:: sumfem.electrochemistry ``` -## Electrolytes +## Materials + +```{eval-rst} +.. automodule:: sumfem.electrochemistry.materials + +.. currentmodule:: sumfem.electrochemistry +``` + +### Electrolytes ```{eval-rst} .. autoclass:: Electrolyte @@ -14,26 +22,32 @@ :members: ``` -## Porous electrodes and their flow +### Redox couples ```{eval-rst} -.. autoclass:: PorousElectrode - :members: - -.. autoclass:: ElectrodeFlow +.. autoclass:: RedoxCouple :members: +``` -.. autoclass:: DarcyFlow - :members: +### Porous electrodes -.. autoclass:: FlowSetup +```{eval-rst} +.. autoclass:: PorousElectrode :members: .. autoclass:: Region :members: ``` -## Charge and transport +## Models + +```{eval-rst} +.. automodule:: sumfem.electrochemistry.models + +.. currentmodule:: sumfem.electrochemistry +``` + +### Charge and transport ```{eval-rst} .. autoclass:: ElectrolyteConduction @@ -46,7 +60,7 @@ :members: ``` -## Kinetics +### Kinetics ```{eval-rst} .. autoclass:: ElectrodeReaction @@ -65,7 +79,7 @@ :members: ``` -## Membranes +### Membranes ```{eval-rst} .. autoclass:: ProtonExchangeMembrane @@ -78,17 +92,35 @@ :members: ``` -## Half-cells +### Flow through the electrodes ```{eval-rst} -.. autoclass:: RedoxCouple +.. autoclass:: ElectrodeFlow + :members: + +.. autoclass:: DarcyFlow + :members: + +.. autoclass:: FlowSetup :members: +``` + +## Cells + +```{eval-rst} +.. automodule:: sumfem.electrochemistry.cells +.. currentmodule:: sumfem.electrochemistry +``` + +### Half-cells + +```{eval-rst} .. autoclass:: RedoxHalfCell :members: ``` -## Full cells +### Full cells ```{eval-rst} .. autoclass:: FullCell diff --git a/src/sumfem/electrochemistry/__init__.py b/src/sumfem/electrochemistry/__init__.py index 179e742..7f9db2b 100644 --- a/src/sumfem/electrochemistry/__init__.py +++ b/src/sumfem/electrochemistry/__init__.py @@ -1,47 +1,45 @@ r"""Models of electrochemistry, built on :mod:`sumfem.models`. -- Electrolytes: :class:`Electrolyte`, a dilute solution of :class:`Ion` - and the constants and properties its models share. -- Porous electrodes: :class:`PorousElectrode`, and the :class:`Region` of - a cell filled with electrolyte. -- Charge and transport in the electrolyte: :class:`ElectrolyteConduction` - and :class:`NernstPlanck`, with storage :class:`TransientNernstPlanck`. -- Kinetics: :class:`ExchangeCurrentButlerVolmer`, :class:`ButlerVolmer`, - :class:`MassTransferButlerVolmer` and :class:`PistonVelocityButlerVolmer`, - on :class:`ElectrodeReaction`, with the Butler-Volmer equation - :func:`~sumfem.models.butler_volmer`. -- Membranes: :class:`ProtonExchangeMembrane`, with the flow of water of - :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`, with the flow of the electrolytes from their - tanks, a :class:`FlowSpec`, by any :class:`ElectrodeFlow`, e.g. - :class:`DarcyFlow`. +The package has three layers, each importing only from those before it: + +- :mod:`~sumfem.electrochemistry.materials`, what the cells are made of: + :class:`Electrolyte` of :class:`Ion`, :class:`RedoxCouple`, + :class:`PorousElectrode` and :class:`Region`. +- :mod:`~sumfem.electrochemistry.models`, the equations: charge and + transport, :class:`ElectrolyteConduction`, :class:`NernstPlanck` and + :class:`TransientNernstPlanck`; kinetics, :class:`ButlerVolmer` and its + variants on :class:`ElectrodeReaction`; membranes, + :class:`ProtonExchangeMembrane`, :class:`SchloeglFlow` and + :class:`MembraneProtonFlux`; and the flow through the electrodes, + :class:`ElectrodeFlow` and :class:`DarcyFlow`. +- :mod:`~sumfem.electrochemistry.cells`, assembled from them: + :class:`RedoxHalfCell` and :class:`FullCell`, with the specifications + :class:`HalfCellSpec`, :class:`MembraneSpec` and :class:`FlowSpec`. + +Every name is also imported from here. The package imports from the core of :mod:`sumfem`, and never the other way around. """ -from .charge import ElectrolyteConduction -from .couple import RedoxCouple -from .electrode import DarcyFlow, ElectrodeFlow, FlowSetup, PorousElectrode -from .electrolyte import Electrolyte, Ion -from .full_cell import FullCell -from .half_cell import RedoxHalfCell -from .kinetics import ( +from .cells import FlowSpec, FullCell, HalfCellSpec, MembraneSpec, RedoxHalfCell +from .materials import Electrolyte, Ion, PorousElectrode, RedoxCouple, Region +from .models import ( ButlerVolmer, + DarcyFlow, + ElectrodeFlow, ElectrodeReaction, + ElectrolyteConduction, ExchangeCurrentButlerVolmer, + FlowSetup, MassTransferButlerVolmer, + MembraneProtonFlux, + NernstPlanck, PistonVelocityButlerVolmer, + ProtonExchangeMembrane, + SchloeglFlow, + TransientNernstPlanck, ) -from .membrane import MembraneProtonFlux, ProtonExchangeMembrane, SchloeglFlow -from .porous import Region -from .specs import FlowSpec, HalfCellSpec, MembraneSpec -from .transport import NernstPlanck, TransientNernstPlanck __all__ = [ "Electrolyte", diff --git a/src/sumfem/electrochemistry/cells/__init__.py b/src/sumfem/electrochemistry/cells/__init__.py new file mode 100644 index 0000000..78107a9 --- /dev/null +++ b/src/sumfem/electrochemistry/cells/__init__.py @@ -0,0 +1,15 @@ +"""Cells of a redox flow battery, assembled from the materials and models. + +- Half-cells: :class:`RedoxHalfCell`, negative or positive, of any redox + couple and electrolyte. +- Full cells: :class:`FullCell`, two half-cells, each a + :class:`HalfCellSpec`, and a membrane between them, a + :class:`MembraneSpec`, with the flow of the electrolytes from their + tanks, a :class:`FlowSpec`. +""" + +from .full_cell import FullCell +from .half_cell import RedoxHalfCell +from .specs import FlowSpec, HalfCellSpec, MembraneSpec + +__all__ = ["FlowSpec", "FullCell", "HalfCellSpec", "MembraneSpec", "RedoxHalfCell"] diff --git a/src/sumfem/electrochemistry/full_cell.py b/src/sumfem/electrochemistry/cells/full_cell.py similarity index 99% rename from src/sumfem/electrochemistry/full_cell.py rename to src/sumfem/electrochemistry/cells/full_cell.py index 3e34a7e..6188c55 100644 --- a/src/sumfem/electrochemistry/full_cell.py +++ b/src/sumfem/electrochemistry/cells/full_cell.py @@ -31,9 +31,9 @@ ) from sumfem.symbolic.si import CONCENTRATION, DENSITY, VISCOSITY -from .electrode import ElectrodeFlow, FlowSetup, PorousElectrode +from ..materials import PorousElectrode +from ..models import ElectrodeFlow, FlowSetup, ProtonExchangeMembrane from .half_cell import RedoxHalfCell -from .membrane import ProtonExchangeMembrane from .specs import PROTONS, WATER, FlowSpec, HalfCellSpec, MembraneSpec __all__ = ["FullCell"] diff --git a/src/sumfem/electrochemistry/half_cell.py b/src/sumfem/electrochemistry/cells/half_cell.py similarity index 92% rename from src/sumfem/electrochemistry/half_cell.py rename to src/sumfem/electrochemistry/cells/half_cell.py index 7d3f94b..c46df01 100644 --- a/src/sumfem/electrochemistry/half_cell.py +++ b/src/sumfem/electrochemistry/cells/half_cell.py @@ -10,6 +10,7 @@ from __future__ import annotations import dataclasses +import math import types import typing @@ -26,24 +27,23 @@ from sumfem.symbolic import ( DirichletBC, Domain, + ExpressionOrDefault, Field, Parameter, Residual, ResidualSum, + expression_or_parameter, ) from sumfem.symbolic.si import CONCENTRATION -from .charge import ElectrolyteConduction -from .couple import RedoxCouple -from .electrode import PorousElectrode -from .electrolyte import Electrolyte, Ion -from .kinetics import ( +from ..materials import Electrolyte, Ion, PorousElectrode, RedoxCouple, Region +from ..models import ( ButlerVolmer, + ElectrolyteConduction, MassTransferButlerVolmer, PistonVelocityButlerVolmer, + TransientNernstPlanck, ) -from .porous import Region -from .transport import TransientNernstPlanck __all__ = ["RedoxHalfCell"] @@ -251,19 +251,19 @@ class RedoxHalfCell: _area_ratio: Parameter """Ratio of the membrane area to the current collector area.""" - _formal_potential: Parameter + _formal_potential: ufl.core.expr.Expr """Formal potential of the couple.""" - _rate_constant: Parameter + _rate_constant: ufl.core.expr.Expr """Standard rate constant of the couple.""" - _anodic_transfer: Parameter + _anodic_transfer: ufl.core.expr.Expr """Anodic transfer coefficient of the couple.""" - _cathodic_transfer: Parameter + _cathodic_transfer: ufl.core.expr.Expr """Cathodic transfer coefficient of the couple.""" - _fiber_distance: Parameter | None + _fiber_distance: ufl.core.expr.Expr | None """Mean distance between the fibres, if the piston velocities limit the kinetics.""" @@ -293,13 +293,13 @@ def __init__( channels: tuple[str, ...] = (), velocity: ufl.core.expr.Expr | None = None, mass_transfer: bool | None = None, - fiber_distance: float | sympy.Expr | None = None, + fiber_distance: ExpressionOrDefault = None, collector: tuple[str, ...] | None = ("CurrentCollector",), membrane: tuple[str, ...] | None = ("Membrane",), outlet: tuple[str, ...] = (), proton_flux: bool = False, - density: ufl.core.expr.Expr | float = 1399.8, - viscosity: ufl.core.expr.Expr | float = 6.963e-3, + density: ExpressionOrDefault = 1399.8, + viscosity: ExpressionOrDefault = 6.963e-3, degree: int = 1, name: str = "cell", phi_s: Field | None = None, @@ -390,8 +390,8 @@ def _create_electrolyte( self, ions: typing.Sequence[Ion], name: str, - density: ufl.core.expr.Expr | float, - viscosity: ufl.core.expr.Expr | float, + density: ExpressionOrDefault, + viscosity: ExpressionOrDefault, ) -> Electrolyte: """Create the electrolyte, with the fields of the transported ions. @@ -423,7 +423,7 @@ def _create_electrolyte( def _create_parameters( self, reservoir: typing.Mapping[str, float | sympy.Expr], - fiber_distance: float | sympy.Expr | None, + fiber_distance: ExpressionOrDefault, ) -> None: """Create the parameters of the cell. @@ -465,39 +465,39 @@ def _create_parameters( default=1.0, doc="Membrane area over current collector area", ) - self._formal_potential = Parameter( + self._formal_potential = expression_or_parameter( f"{name}.formal_potential", + couple.formal_potential, domain, - default=couple.formal_potential, unit=syu.volt, doc=f"Formal potential of {pair}", ) - self._rate_constant = Parameter( + self._rate_constant = expression_or_parameter( f"{name}.rate_constant", + couple.rate_constant, domain, - default=couple.rate_constant, unit=syu.m / syu.second, doc=f"Standard rate constant of {pair}", ) - self._anodic_transfer = Parameter( + self._anodic_transfer = expression_or_parameter( f"{name}.anodic_transfer", + couple.anodic_transfer, domain, - default=couple.anodic_transfer, doc=f"Anodic transfer coefficient of {pair}", ) - self._cathodic_transfer = Parameter( + self._cathodic_transfer = expression_or_parameter( f"{name}.cathodic_transfer", + couple.cathodic_transfer, domain, - default=couple.cathodic_transfer, doc=f"Cathodic transfer coefficient of {pair}", ) self._fiber_distance = ( None if fiber_distance is None - else Parameter( + else expression_or_parameter( f"{name}.fiber_distance", + fiber_distance, domain, - default=fiber_distance, unit=syu.m, doc="Mean distance between the fibres", ) @@ -641,12 +641,14 @@ def _create_models( return models def _create_functionals(self) -> types.MappingProxyType[str, ufl.Form]: - """Create the integrals of the voltage and the state of charge. + """Create the integrals of the voltage, the potentials and the state + of charge. Returns: The areas and potential integrals of the collector and the - membrane, when they have markers, and the amounts of the reduced - and oxidized ion, by name. + membrane, when they have markers, the volume of the electrode and + the integrals over it of the formal and the equilibrium potential, + and the amounts of the reduced and oxidized ion, by name. """ domain = self._domain functionals: dict[str, ufl.Form] = {} @@ -658,6 +660,17 @@ def _create_functionals(self) -> types.MappingProxyType[str, ufl.Form]: ds_membrane = domain.ds(*self._membrane) functionals["membrane_area"] = 1.0 * ds_membrane functionals["membrane_potential"] = self._phi_l * ds_membrane + couple, e = self._couple, self._electrolyte + dx = domain.dx(*self._electrode.markers) + ratio = self.reservoir_concentration( + couple.oxidized + ) / self.reservoir_concentration(couple.reduced) + thermal = e.R * e.T / (couple.electrons * e.F) + functionals["electrode_volume"] = 1.0 * dx + functionals["formal_potential"] = self._formal_potential * dx + functionals["equilibrium_potential"] = ( + self._formal_potential + thermal * math.log(ratio) + ) * dx for role, ion in ( ("reduced", self._couple.reduced), ("oxidized", self._couple.oxidized), @@ -720,7 +733,7 @@ def current(self) -> Parameter: return self._current @property - def formal_potential(self) -> Parameter: + def formal_potential(self) -> ufl.core.expr.Expr: """Formal potential of the couple.""" return self._formal_potential @@ -755,36 +768,38 @@ def equilibrium_potential(self, discretization: Discretization) -> float: That is ``E0 + R T / (n F) ln(c_ox / c_red)``, see :meth:`~sumfem.electrochemistry.ButlerVolmer.equilibrium_potential`, with the reservoir concentrations as given, - :meth:`reservoir_concentration`. + :meth:`reservoir_concentration`, averaged over the electrode, so that + a formal potential or a temperature given as an expression, e.g. of + a field, counts with its values on the discretization. Args: discretization: The discretization, for the values of the - formal potential and the constants. + formal potential, the constants and any fields they hold. Returns: - The potential in V. + The mean potential in V. """ - value, couple, e = discretization.value, self._couple, self._electrolyte - # The half-cell's electrolyte has the temperature as a parameter - assert isinstance(e.T, Parameter) - thermal = value(e.R) * value(e.T) / (couple.electrons * value(e.F)) - ratio = self.reservoir_concentration( - couple.oxidized - ) / self.reservoir_concentration(couple.reduced) - return float(value(self._formal_potential) + thermal * np.log(ratio)) + assemble = discretization.assemble_scalar + f = self._functionals + return assemble(f["equilibrium_potential"]) / assemble(f["electrode_volume"]) def initialize(self, discretization: Discretization) -> None: """Set a state at rest, as RFBniCS starts. The concentrations start at their reservoir values, the solid at zero - and the electrolyte at minus the formal potential. + and the electrolyte at minus the formal potential, its mean over the + electrode if it is an expression. """ values: dict[str | ufl.Coefficient, typing.Any] = { c: discretization.value(self._inlet[ion]) for ion, c in self._concentrations.items() } + assemble = discretization.assemble_scalar + f = self._functionals values[self._phi_s] = 0.0 - values[self._phi_l] = -discretization.value(self._formal_potential) + values[self._phi_l] = -assemble(f["formal_potential"]) / assemble( + f["electrode_volume"] + ) discretization.set(values) def _check_boundaries(self, what: str) -> None: diff --git a/src/sumfem/electrochemistry/specs.py b/src/sumfem/electrochemistry/cells/specs.py similarity index 87% rename from src/sumfem/electrochemistry/specs.py rename to src/sumfem/electrochemistry/cells/specs.py index cd475a3..78053e7 100644 --- a/src/sumfem/electrochemistry/specs.py +++ b/src/sumfem/electrochemistry/cells/specs.py @@ -14,9 +14,8 @@ from sumfem.symbolic import ExpressionOrDefault -from .couple import RedoxCouple -from .electrode import DarcyFlow, ElectrodeFlow, FlowSetup -from .electrolyte import Ion +from ..materials import Ion, RedoxCouple +from ..models import DarcyFlow, ElectrodeFlow, FlowSetup __all__ = ["FlowSpec", "HalfCellSpec", "MembraneSpec"] @@ -31,6 +30,13 @@ class HalfCellSpec: """One side of a full cell: its electrolyte, its couple and its electrode. + The properties of the electrode and the electrolyte, from the porosity + on, are each an expression, used as given, e.g. a field of another model + such as the porosity of a compressed felt or a viscosity of the + temperature, or the default of the parameter the cell creates, see + :func:`~sumfem.symbolic.expression_or_parameter`. The reservoir + concentrations are values: with flow, the tanks change them. + Attributes: ions: The ions of the electrolyte, see :class:`RedoxHalfCell`. Water, ``"H2O"``, if transported, must be an ion of both sides, @@ -38,9 +44,6 @@ class HalfCellSpec: reservoir: The initial concentration of each transported ion, see :class:`RedoxHalfCell`, and with flow, that of its tank. The protons, ``"H"``, must be one of them. - - Raises: - ValueError: If the side does not transport protons. 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 @@ -61,19 +64,22 @@ class HalfCellSpec: density: The density of the electrolyte, in kg/m^3 or with units. viscosity: The dynamic viscosity of the electrolyte, in Pa s or with units, which the flow and the mass transfer share. + + Raises: + ValueError: If the side does not transport protons. """ 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 - fiber_diameter: float | sympy.Expr | None = None - kozeny_carman_constant: float | None = None - density: float | sympy.Expr = 1399.8 - viscosity: float | sympy.Expr = 6.963e-3 + porosity: ExpressionOrDefault + conductivity: ExpressionOrDefault + specific_area: ExpressionOrDefault + fiber_distance: ExpressionOrDefault = None + fiber_diameter: ExpressionOrDefault = None + kozeny_carman_constant: ExpressionOrDefault = None + density: ExpressionOrDefault = 1399.8 + viscosity: ExpressionOrDefault = 6.963e-3 def __post_init__(self) -> None: if PROTONS not in self.reservoir: diff --git a/src/sumfem/electrochemistry/materials/__init__.py b/src/sumfem/electrochemistry/materials/__init__.py new file mode 100644 index 0000000..bc771de --- /dev/null +++ b/src/sumfem/electrochemistry/materials/__init__.py @@ -0,0 +1,19 @@ +"""What the cells are made of: ions, electrolytes, couples and electrodes. + +These describe materials and their properties, as parameters and the +expressions of them that the models share, and hold no equations of their +own. + +- Electrolytes: :class:`Electrolyte`, a dilute solution of :class:`Ion`. +- Redox couples: :class:`RedoxCouple`, the reaction of a half-cell. +- Porous electrodes: :class:`PorousElectrode`, and the :class:`Region` of + a cell filled with electrolyte. + +The materials import nothing from the models or the cells. +""" + +from .couple import RedoxCouple +from .electrode import PorousElectrode, Region +from .electrolyte import Electrolyte, Ion + +__all__ = ["Electrolyte", "Ion", "PorousElectrode", "RedoxCouple", "Region"] diff --git a/src/sumfem/electrochemistry/couple.py b/src/sumfem/electrochemistry/materials/couple.py similarity index 88% rename from src/sumfem/electrochemistry/couple.py rename to src/sumfem/electrochemistry/materials/couple.py index eb93950..8736b7e 100644 --- a/src/sumfem/electrochemistry/couple.py +++ b/src/sumfem/electrochemistry/materials/couple.py @@ -6,7 +6,7 @@ import types import typing -import sympy +from sumfem.symbolic import ExpressionOrDefault from .electrolyte import Ion @@ -26,6 +26,12 @@ class RedoxCouple: :class:`~sumfem.electrochemistry.ButlerVolmer`, in the concentrations of the couple. + The rate constant, the formal potential and the transfer coefficients + are each an expression, used as given, e.g. a rate constant of a + temperature field by Arrhenius' law, or the default of the parameter of + the half-cell of the couple, see + :func:`~sumfem.symbolic.expression_or_parameter`. + Attributes: reduced: Name of the reduced ion, e.g. ``"V2"``. oxidized: Name of the oxidized ion, e.g. ``"V3"``. @@ -50,12 +56,12 @@ class RedoxCouple: reduced: str oxidized: str - rate_constant: float | sympy.Expr - formal_potential: float | sympy.Expr + rate_constant: ExpressionOrDefault + formal_potential: ExpressionOrDefault oxidation_on_discharge: bool electrons: int = 1 - anodic_transfer: float = 0.5 - cathodic_transfer: float = 0.5 + anodic_transfer: ExpressionOrDefault = 0.5 + cathodic_transfer: ExpressionOrDefault = 0.5 reactants: typing.Mapping[str, int] | None = None products: typing.Mapping[str, int] | None = None diff --git a/src/sumfem/electrochemistry/materials/electrode.py b/src/sumfem/electrochemistry/materials/electrode.py new file mode 100644 index 0000000..6dea808 --- /dev/null +++ b/src/sumfem/electrochemistry/materials/electrode.py @@ -0,0 +1,194 @@ +"""Porous electrodes, and the regions of a cell filled with electrolyte. + +:class:`PorousElectrode` is an electrode of carbon fibres filled with +electrolyte: its porosity, the conductivity of its fibres, its specific +area and the permeability of its fibres to the flow. A :class:`Region` is +a part of a cell filled with electrolyte, the pores of an electrode or a +free channel, with the porosity that weights storage and corrects +transport. +""" + +from __future__ import annotations + +import dataclasses + +import sympy.physics.units as syu +import ufl + +from sumfem.models import bruggeman, ergun_coefficient, kozeny_carman +from sumfem.symbolic import Domain, ExpressionOrDefault, expression_or_parameter + +__all__ = ["PorousElectrode", "Region"] + + +@dataclasses.dataclass(frozen=True) +class Region: + """Cells filled with electrolyte: free electrolyte, or an electrode's pores. + + Attributes: + markers: Cell markers of the region. Empty for the whole domain. + porosity: Porosity of a porous electrode. Defaults to one, for free + electrolyte. + bruggeman_exponent: Exponent of the Bruggeman correction of + transport in the pores. + """ + + markers: tuple[str, ...] = () + porosity: ufl.core.expr.Expr | float = 1.0 + bruggeman_exponent: float = 1.5 + + @property + def storage(self) -> ufl.core.expr.Expr | float: + """Volume fraction of electrolyte, the porosity.""" + return self.porosity + + @property + def transport(self) -> ufl.core.expr.Expr | float: + r"""Correction of diffusivities, Bruggeman's :math:`\varepsilon^b`. + + It is one in free electrolyte, with porosity one. + """ + return bruggeman(1.0, self.porosity, self.bruggeman_exponent) + + +class PorousElectrode: + """A porous electrode of carbon fibers, filled with electrolyte. + + Each property is an expression, used as given, e.g. a field of another + model such as the porosity of a compressed felt or a conductivity of the + temperature, or the default of the parameter ``"."``, + such as ``".porosity"``, in SI units or a quantity with units, + e.g. ``0.5 * syu.micrometer``, see + :func:`~sumfem.symbolic.expression_or_parameter`. Models that share an + electrode share these. A parameter without a default value must be set + on the discretization before it is used. + + Args: + name: Prefix of the parameter names. + domain: The domain of the parameters. + markers: Cell markers of the electrode. Empty for the whole domain. + porosity: The porosity. + conductivity: The conductivity of the fiber material, in S/m. + specific_area: The specific area, in 1/m. + fiber_diameter: The fiber diameter, in m. + kozeny_carman_constant: The Kozeny-Carman constant. + pore_length: The characteristic pore length of Ergun's drag, in m. + """ + + _markers: tuple[str, ...] + """Cell markers of the electrode.""" + + _porosity: ufl.core.expr.Expr + """Porosity.""" + + _conductivity: ufl.core.expr.Expr + """Conductivity of the fiber material.""" + + _specific_area: ufl.core.expr.Expr + """Specific area.""" + + _fiber_diameter: ufl.core.expr.Expr + """Fiber diameter.""" + + _kozeny_carman_constant: ufl.core.expr.Expr + """Kozeny-Carman constant.""" + + _pore_length: ufl.core.expr.Expr + """Characteristic pore length of Ergun's drag.""" + + def __init__( + self, + name: str, + domain: Domain, + markers: tuple[str, ...] = (), + porosity: ExpressionOrDefault = None, + conductivity: ExpressionOrDefault = None, + specific_area: ExpressionOrDefault = None, + fiber_diameter: ExpressionOrDefault = None, + kozeny_carman_constant: ExpressionOrDefault = None, + pore_length: ExpressionOrDefault = None, + ): + self._markers = tuple(markers) + self._porosity = expression_or_parameter( + f"{name}.porosity", porosity, domain, doc="Porosity" + ) + self._conductivity = expression_or_parameter( + f"{name}.conductivity", + conductivity, + domain, + unit=syu.siemens / syu.m, + doc="Conductivity of the fiber material", + ) + self._specific_area = expression_or_parameter( + f"{name}.specific_area", + specific_area, + domain, + unit=1 / syu.m, + doc="Specific area", + ) + self._fiber_diameter = expression_or_parameter( + f"{name}.fiber_diameter", + fiber_diameter, + domain, + unit=syu.m, + doc="Fiber diameter", + ) + self._kozeny_carman_constant = expression_or_parameter( + f"{name}.kozeny_carman_constant", + kozeny_carman_constant, + domain, + doc="Kozeny-Carman constant", + ) + self._pore_length = expression_or_parameter( + f"{name}.pore_length", + pore_length, + domain, + unit=syu.m, + doc="Characteristic pore length of Ergun's drag", + ) + + @property + def markers(self) -> tuple[str, ...]: + """Cell markers of the electrode.""" + return self._markers + + @property + def region(self) -> Region: + """The electrolyte-filled pores, as a region.""" + return Region(self._markers, self._porosity) + + @property + def porosity(self) -> ufl.core.expr.Expr: + """Porosity.""" + return self._porosity + + @property + def specific_area(self) -> ufl.core.expr.Expr: + """Specific area.""" + return self._specific_area + + @property + def fiber_diameter(self) -> ufl.core.expr.Expr: + """Fiber diameter.""" + return self._fiber_diameter + + def effective_conductivity(self) -> ufl.core.expr.Expr: + r"""Bruggeman's conductivity of the solid. + + That is :math:`\sigma_\text{eff} = (1 - \varepsilon)^{3/2} \sigma`, with + the conductivity :math:`\sigma` of the fibres and the porosity + :math:`\varepsilon`. + """ + return bruggeman(self._conductivity, 1.0 - self._porosity) + + def permeability(self) -> ufl.core.expr.Expr: + """Kozeny-Carman permeability of the fibers, see + :func:`~sumfem.models.kozeny_carman`.""" + return kozeny_carman( + self._fiber_diameter, self._porosity, self._kozeny_carman_constant + ) + + def forchheimer(self) -> ufl.core.expr.Expr: + """Ergun's inertial drag coefficient, an estimate for fibres, see + :func:`~sumfem.models.ergun_coefficient`.""" + return ergun_coefficient(self._porosity, self._pore_length) diff --git a/src/sumfem/electrochemistry/electrolyte.py b/src/sumfem/electrochemistry/materials/electrolyte.py similarity index 98% rename from src/sumfem/electrochemistry/electrolyte.py rename to src/sumfem/electrochemistry/materials/electrolyte.py index 0fc57d4..0a89c51 100644 --- a/src/sumfem/electrochemistry/electrolyte.py +++ b/src/sumfem/electrochemistry/materials/electrolyte.py @@ -14,6 +14,7 @@ FARADAY, GAS_CONSTANT, Domain, + ExpressionOrDefault, Field, Parameter, expression_or_parameter, @@ -181,9 +182,9 @@ def __init__( name: str, domain: Domain, ions: typing.Sequence[Ion] = (), - temperature: ufl.core.expr.Expr | float = 298.15, - density: ufl.core.expr.Expr | float | None = None, - viscosity: ufl.core.expr.Expr | float | None = None, + temperature: ExpressionOrDefault = 298.15, + density: ExpressionOrDefault = None, + viscosity: ExpressionOrDefault = None, concentration_floor: float = 1e-14, ): names = [ion.name for ion in ions] diff --git a/src/sumfem/electrochemistry/models/__init__.py b/src/sumfem/electrochemistry/models/__init__.py new file mode 100644 index 0000000..d1ee51d --- /dev/null +++ b/src/sumfem/electrochemistry/models/__init__.py @@ -0,0 +1,46 @@ +r"""The equations of electrochemistry: residuals of the materials' fields. + +- Charge and transport in the electrolyte: :class:`ElectrolyteConduction` + and :class:`NernstPlanck`, with storage :class:`TransientNernstPlanck`. +- Kinetics: :class:`ExchangeCurrentButlerVolmer`, :class:`ButlerVolmer`, + :class:`MassTransferButlerVolmer` and :class:`PistonVelocityButlerVolmer`, + on :class:`ElectrodeReaction`, with the Butler-Volmer equation + :func:`~sumfem.models.butler_volmer`. +- Membranes: :class:`ProtonExchangeMembrane`, with the flow of water of + :class:`SchloeglFlow`, whose current the protons of an electrode carry + through its face, :class:`MembraneProtonFlux`. +- The flow of the electrolyte through the electrodes: any + :class:`ElectrodeFlow`, e.g. :class:`DarcyFlow`, created from a + :class:`FlowSetup`. + +The models import from the materials, and nothing from the cells. +""" + +from .charge import ElectrolyteConduction +from .electrode_flow import DarcyFlow, ElectrodeFlow, FlowSetup +from .kinetics import ( + ButlerVolmer, + ElectrodeReaction, + ExchangeCurrentButlerVolmer, + MassTransferButlerVolmer, + PistonVelocityButlerVolmer, +) +from .membrane import MembraneProtonFlux, ProtonExchangeMembrane, SchloeglFlow +from .transport import NernstPlanck, TransientNernstPlanck + +__all__ = [ + "ButlerVolmer", + "DarcyFlow", + "ElectrodeFlow", + "ElectrodeReaction", + "ElectrolyteConduction", + "ExchangeCurrentButlerVolmer", + "FlowSetup", + "MassTransferButlerVolmer", + "MembraneProtonFlux", + "NernstPlanck", + "PistonVelocityButlerVolmer", + "ProtonExchangeMembrane", + "SchloeglFlow", + "TransientNernstPlanck", +] diff --git a/src/sumfem/electrochemistry/charge.py b/src/sumfem/electrochemistry/models/charge.py similarity index 89% rename from src/sumfem/electrochemistry/charge.py rename to src/sumfem/electrochemistry/models/charge.py index fccd567..7572904 100644 --- a/src/sumfem/electrochemistry/charge.py +++ b/src/sumfem/electrochemistry/models/charge.py @@ -6,8 +6,7 @@ from sumfem.symbolic import Field, Residual -from .electrolyte import Electrolyte -from .porous import Region +from ..materials import Electrolyte, Region __all__ = ["ElectrolyteConduction"] @@ -18,9 +17,9 @@ class ElectrolyteConduction(Residual): 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` + :meth:`~sumfem.electrochemistry.Electrolyte.current_density` with the diffusivities corrected for the region, :math:`t` the region's - :attr:`~sumfem.electrochemistry.porous.Region.transport`, as in + :attr:`~sumfem.electrochemistry.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. diff --git a/src/sumfem/electrochemistry/electrode.py b/src/sumfem/electrochemistry/models/electrode_flow.py similarity index 61% rename from src/sumfem/electrochemistry/electrode.py rename to src/sumfem/electrochemistry/models/electrode_flow.py index b0f440c..d72f24e 100644 --- a/src/sumfem/electrochemistry/electrode.py +++ b/src/sumfem/electrochemistry/models/electrode_flow.py @@ -1,14 +1,11 @@ -"""Porous electrodes, and the flow of the electrolyte through them. - -:class:`PorousElectrode` is an electrode of carbon fibres filled with -electrolyte: its porosity, the conductivity of its fibres, its specific -area and the permeability of its fibres to the flow. The flow of the -electrolytes through the electrodes of a -:class:`~sumfem.electrochemistry.FullCell` is a component, an -:class:`ElectrodeFlow`, which the cell creates from a :class:`FlowSetup` of -what it knows. :class:`DarcyFlow` is Darcy flow through the electrodes; -another flow, e.g. Brinkman-Forchheimer flow or transient Stokes flow, is -another subclass, given to the cell's +"""The flow of the electrolyte through the electrodes of a full cell. + +A :class:`~sumfem.electrochemistry.FullCell` takes its flow as a component: +an :class:`ElectrodeFlow`, which the cell creates from a :class:`FlowSetup` +of what it knows, e.g. the electrodes and their subdomains, the +electrolytes' viscosities and the inlet velocity. :class:`DarcyFlow` is +Darcy flow through the electrodes; another flow, e.g. Brinkman-Forchheimer +flow or transient Stokes flow, is another subclass, given to the cell's :class:`~sumfem.electrochemistry.FlowSpec`. """ @@ -18,19 +15,12 @@ import dataclasses import typing -import sympy import sympy.physics.units as syu import ufl from sumfem.discrete import Discretization from sumfem.error_handling import ComponentError -from sumfem.models import ( - BoundaryFlux, - Darcy, - bruggeman, - ergun_coefficient, - kozeny_carman, -) +from sumfem.models import BoundaryFlux, Darcy from sumfem.symbolic import ( DirichletBC, Domain, @@ -41,152 +31,9 @@ value_at_defaults, ) -from .porous import Region - -__all__ = ["DarcyFlow", "ElectrodeFlow", "FlowSetup", "PorousElectrode"] - - -class PorousElectrode: - """A porous electrode of carbon fibers, filled with electrolyte. - - The electrode creates a parameter for each argument below from the - porosity on, named ``"."``, such as - ``".porosity"``, with the argument as its default value. Models - that share an electrode share these parameters. A parameter without a - default value must be set on the discretization before it is used. - - Args: - name: Prefix of the parameter names. - domain: The domain of the parameters. - markers: Cell markers of the electrode. Empty for the whole domain. - porosity: Default porosity, if any. Every default with a unit is in - SI units, or a quantity with units, e.g. ``0.5 * syu.micrometer``. - conductivity: Default conductivity of the fiber material in S/m, if - any. - specific_area: Default specific area in 1/m, if any. - fiber_diameter: Default fiber diameter in m, if any. - kozeny_carman_constant: Default Kozeny-Carman constant, if any. - pore_length: Default characteristic pore length of Ergun's drag in - m, if any. - """ - - _markers: tuple[str, ...] - """Cell markers of the electrode.""" - - _porosity: Parameter - """Porosity.""" - - _conductivity: Parameter - """Conductivity of the fiber material.""" - - _specific_area: Parameter - """Specific area.""" - - _fiber_diameter: Parameter - """Fiber diameter.""" - - _kozeny_carman_constant: Parameter - """Kozeny-Carman constant.""" - - _pore_length: Parameter - """Characteristic pore length of Ergun's drag.""" - - def __init__( - self, - name: str, - domain: Domain, - markers: tuple[str, ...] = (), - porosity: float | None = None, - conductivity: float | sympy.Expr | None = None, - specific_area: float | sympy.Expr | None = None, - fiber_diameter: float | sympy.Expr | None = None, - kozeny_carman_constant: float | None = None, - pore_length: float | sympy.Expr | None = None, - ): - self._markers = tuple(markers) - self._porosity = Parameter( - f"{name}.porosity", domain, default=porosity, doc="Porosity" - ) - self._conductivity = Parameter( - f"{name}.conductivity", - domain, - default=conductivity, - unit=syu.siemens / syu.m, - doc="Conductivity of the fiber material", - ) - self._specific_area = Parameter( - f"{name}.specific_area", - domain, - default=specific_area, - unit=1 / syu.m, - doc="Specific area", - ) - self._fiber_diameter = Parameter( - f"{name}.fiber_diameter", - domain, - default=fiber_diameter, - unit=syu.m, - doc="Fiber diameter", - ) - self._kozeny_carman_constant = Parameter( - f"{name}.kozeny_carman_constant", - domain, - default=kozeny_carman_constant, - doc="Kozeny-Carman constant", - ) - self._pore_length = Parameter( - f"{name}.pore_length", - domain, - default=pore_length, - unit=syu.m, - doc="Characteristic pore length of Ergun's drag", - ) - - @property - def markers(self) -> tuple[str, ...]: - """Cell markers of the electrode.""" - return self._markers - - @property - def region(self) -> Region: - """The electrolyte-filled pores, as a region.""" - return Region(self._markers, self._porosity) - - @property - def porosity(self) -> Parameter: - """Porosity.""" - return self._porosity - - @property - def specific_area(self) -> Parameter: - """Specific area.""" - return self._specific_area - - @property - def fiber_diameter(self) -> Parameter: - """Fiber diameter.""" - return self._fiber_diameter - - def effective_conductivity(self) -> ufl.core.expr.Expr: - r"""Bruggeman's conductivity of the solid. - - That is :math:`\sigma_\text{eff} = (1 - \varepsilon)^{3/2} \sigma`, with - the conductivity :math:`\sigma` of the fibres and the porosity - :math:`\varepsilon`. - """ - return bruggeman(self._conductivity, 1.0 - self._porosity) - - def permeability(self) -> ufl.core.expr.Expr: - """Kozeny-Carman permeability of the fibers, see - :func:`~sumfem.models.kozeny_carman`.""" - return kozeny_carman( - self._fiber_diameter, self._porosity, self._kozeny_carman_constant - ) +from ..materials import PorousElectrode - def forchheimer(self) -> ufl.core.expr.Expr: - """Ergun's inertial drag coefficient, an estimate for fibres, see - :func:`~sumfem.models.ergun_coefficient`.""" - return ergun_coefficient(self._porosity, self._pore_length) +__all__ = ["DarcyFlow", "ElectrodeFlow", "FlowSetup"] @dataclasses.dataclass(frozen=True) diff --git a/src/sumfem/electrochemistry/kinetics.py b/src/sumfem/electrochemistry/models/kinetics.py similarity index 99% rename from src/sumfem/electrochemistry/kinetics.py rename to src/sumfem/electrochemistry/models/kinetics.py index cdbffab..cc18cf5 100644 --- a/src/sumfem/electrochemistry/kinetics.py +++ b/src/sumfem/electrochemistry/models/kinetics.py @@ -15,7 +15,7 @@ from sumfem.models import RateReaction, butler_volmer from sumfem.symbolic import ExpressionOrValue, Field, as_expr -from .electrolyte import Electrolyte, Ion +from ..materials import Electrolyte, Ion __all__ = [ "ButlerVolmer", @@ -277,7 +277,7 @@ class ButlerVolmer(ElectrodeReaction): of the reduced and oxidized ion at the electrode surface. These are, by default, the smooth positive parts of the bulk concentrations, :math:`c^+ = (c + \sqrt{c^2 + c_\text{floor}^2}) / 2`, see - :meth:`~sumfem.electrochemistry.electrolyte.Electrolyte.positive_part`. + :meth:`~sumfem.electrochemistry.Electrolyte.positive_part`. The reduced and oxidized ions are ions in the electrolyte, given by name. Subclasses change the kinetics by overriding diff --git a/src/sumfem/electrochemistry/membrane.py b/src/sumfem/electrochemistry/models/membrane.py similarity index 99% rename from src/sumfem/electrochemistry/membrane.py rename to src/sumfem/electrochemistry/models/membrane.py index 4325379..baf82c6 100644 --- a/src/sumfem/electrochemistry/membrane.py +++ b/src/sumfem/electrochemistry/models/membrane.py @@ -22,7 +22,7 @@ ) from sumfem.symbolic.si import CONCENTRATION, DIFFUSIVITY, PERMEABILITY, VISCOSITY -from .electrolyte import Electrolyte +from ..materials import Electrolyte __all__ = ["MembraneProtonFlux", "ProtonExchangeMembrane", "SchloeglFlow"] diff --git a/src/sumfem/electrochemistry/transport.py b/src/sumfem/electrochemistry/models/transport.py similarity index 96% rename from src/sumfem/electrochemistry/transport.py rename to src/sumfem/electrochemistry/models/transport.py index 694bd87..8bdb322 100644 --- a/src/sumfem/electrochemistry/transport.py +++ b/src/sumfem/electrochemistry/models/transport.py @@ -7,8 +7,7 @@ from sumfem.models import Accumulation from sumfem.symbolic import Field, Residual, ResidualSum, as_expr -from .electrolyte import Electrolyte, Ion -from .porous import Region +from ..materials import Electrolyte, Ion, Region __all__ = ["NernstPlanck", "TransientNernstPlanck"] @@ -29,7 +28,7 @@ class NernstPlanck(Residual): with its diffusivity :math:`D` corrected for the region, :math:`D_\text{eff} = t D` with :math:`t` the region's - :attr:`~sumfem.electrochemistry.porous.Region.transport`, the electrolyte + :attr:`~sumfem.electrochemistry.Region.transport`, the electrolyte potential :math:`\phi` and the superficial velocity :math:`\mathbf{u}`. For transient transport, see :class:`TransientNernstPlanck`. @@ -126,7 +125,7 @@ class TransientNernstPlanck(ResidualSum): e.g. of an :class:`~sumfem.electrochemistry.ElectrodeReaction`, are added to it as models of their own. It is the sum of ``"_storage"``, an :class:`~sumfem.models.Accumulation` weighted - with the region's :attr:`~sumfem.electrochemistry.porous.Region.storage`, + with the region's :attr:`~sumfem.electrochemistry.Region.storage`, and ``"_transport"``, a :class:`NernstPlanck`. Args: diff --git a/src/sumfem/electrochemistry/porous.py b/src/sumfem/electrochemistry/porous.py deleted file mode 100644 index 9ee46c4..0000000 --- a/src/sumfem/electrochemistry/porous.py +++ /dev/null @@ -1,41 +0,0 @@ -"""The regions of a cell filled with electrolyte, free or in pores.""" - -from __future__ import annotations - -import dataclasses - -import ufl - -from sumfem.models import bruggeman - -__all__ = ["Region"] - - -@dataclasses.dataclass(frozen=True) -class Region: - """Cells filled with electrolyte: free electrolyte, or an electrode's pores. - - Attributes: - markers: Cell markers of the region. Empty for the whole domain. - porosity: Porosity of a porous electrode. Defaults to one, for free - electrolyte. - bruggeman_exponent: Exponent of the Bruggeman correction of - transport in the pores. - """ - - markers: tuple[str, ...] = () - porosity: ufl.core.expr.Expr | float = 1.0 - bruggeman_exponent: float = 1.5 - - @property - def storage(self) -> ufl.core.expr.Expr | float: - """Volume fraction of electrolyte, the porosity.""" - return self.porosity - - @property - def transport(self) -> ufl.core.expr.Expr | float: - r"""Correction of diffusivities, Bruggeman's :math:`\varepsilon^b`. - - It is one in free electrolyte, with porosity one. - """ - return bruggeman(1.0, self.porosity, self.bruggeman_exponent) diff --git a/tests/electrochemistry/test_full_cell.py b/tests/electrochemistry/test_full_cell.py index a5207a4..e7393b8 100644 --- a/tests/electrochemistry/test_full_cell.py +++ b/tests/electrochemistry/test_full_cell.py @@ -9,6 +9,7 @@ import numpy as np import pytest import sympy.physics.units as syu +import ufl from sumfem.discrete import DomainMaps, MeshData, discretize from sumfem.electrochemistry import ( @@ -20,11 +21,12 @@ MembraneSpec, RedoxCouple, ) -from sumfem.models import Accumulation +from sumfem.models import Accumulation, Diffusion from sumfem.solvers import BDF, Simulation from sumfem.symbolic import ( Clock, Domain, + Field, Parameter, build_system, value_at_defaults, @@ -521,3 +523,75 @@ def test_any_flow_plugs_into_the_cell(): 1.0 * inlet ) assert 0 < mean < 0.99 * drop + + +def test_inputs_can_be_fields_of_other_models(): + """The properties of the electrodes and the couples may be expressions of + the fields of other models: here a temperature, which sets the + conductivity of the fibres, an Arrhenius rate constant and a formal + potential that changes with it, and a compression, which sets the + porosity. Their parameters are then not created, the units are + consistent, and a cell at rest starts and stays at the Nernst voltage of + that temperature.""" + domain = Domain("cell", "interval", length_scale=4e-3) + T = Field("T", domain, ("Lagrange", 1), unit=syu.kelvin, scale=300.0) + strain = Field("strain", domain, ("Lagrange", 1)) + reference = 298.15 * syu.kelvin + slope = 1e-3 * syu.volt / syu.kelvin + activation = 4e3 * syu.kelvin + sides = [] + for spec in (NEGATIVE, POSITIVE): + couple = dataclasses.replace( + spec.couple, + rate_constant=spec.couple.rate_constant + * ufl.exp(-activation * (1 / T - 1 / reference)), + formal_potential=spec.couple.formal_potential + slope * (T - reference), + reactants=None if spec is NEGATIVE else spec.couple.reactants, + products=None if spec is NEGATIVE else spec.couple.products, + ) + sides.append( + dataclasses.replace( + spec, + couple=couple, + porosity=0.68 * (1 - 0.1 * strain), + conductivity=500.0 * syu.siemens / syu.m * (T / reference), + ) + ) + cell = FullCell(domain, *sides, membrane=MEMBRANE) + heat = Diffusion( + "heat", + T, + Parameter("heat.k", domain, default=0.6, unit=syu.watt / (syu.m * syu.kelvin)), + ) + compression = Diffusion("compression", strain, 1.0) + model = cell.model + heat + compression + system = build_system(model, check_units=True) + assert "negative_electrode.porosity" not in system.parameters + assert "positive.rate_constant" not in system.parameters + assert {T, strain} <= set(system.unknowns) + + disc = discretize( + [model, *cell.conditions()], + DomainMaps({domain: create_cell_mesh()}), + max_quadrature_degree=4, + ) + temperature = 308.15 + disc.set({T: temperature, strain: 0.1}) + cell.initialize(disc) + simulation = Simulation( + disc, + model, + cell.conditions(), + clock=Clock(domain, dt=1.0), + scheme=BDF(1), + petsc_options={"snes_linesearch_type": "bt", "snes_max_it": 50}, + check_units=True, + factorize_equations=True, + time_scale=1.0, + ) + RT_F = 8.314462618 * temperature / 96485.33212 + shift = 2 * 1e-3 * (temperature - 298.15) + expected = 1.004 + 0.255 + shift + 2 * RT_F * np.log(27.0 / 1053.0) + assert cell.voltage(disc) == pytest.approx(expected, rel=1e-9) + simulation.run(num_steps=1) + assert cell.voltage(disc) == pytest.approx(expected, rel=1e-8) From 023b031d6e02e85bd5a908074782e76824c03baa Mon Sep 17 00:00:00 2001 From: Jorgen Schartum Dokken Date: Thu, 12 Nov 2026 14:29:33 +0000 Subject: [PATCH 7/7] Add temperature as input to cell --- .../electrochemistry/cells/full_cell.py | 2 + .../electrochemistry/cells/half_cell.py | 10 +++- src/sumfem/electrochemistry/cells/specs.py | 6 ++ .../electrochemistry/materials/electrolyte.py | 57 +++++++++++++------ .../electrochemistry/models/membrane.py | 2 +- tests/electrochemistry/test_full_cell.py | 55 +++++++++++------- 6 files changed, 92 insertions(+), 40 deletions(-) diff --git a/src/sumfem/electrochemistry/cells/full_cell.py b/src/sumfem/electrochemistry/cells/full_cell.py index 6188c55..e910b48 100644 --- a/src/sumfem/electrochemistry/cells/full_cell.py +++ b/src/sumfem/electrochemistry/cells/full_cell.py @@ -201,6 +201,7 @@ def __init__( electrokinetic_permeability=membrane.electrokinetic_permeability, hydraulic_permeability=membrane.hydraulic_permeability, viscosity=membrane.viscosity, + temperature=membrane.temperature, ) self._collector_conductivity = expression_or_parameter( f"{name}.collector_conductivity", @@ -392,6 +393,7 @@ def _create_half_cell( outlet=() if flow is None else (f"{prefix}Outlet",), density=density, viscosity=viscosity, + temperature=spec.temperature, name=side, phi_s=phi_s, phi_l=self._phi_l, diff --git a/src/sumfem/electrochemistry/cells/half_cell.py b/src/sumfem/electrochemistry/cells/half_cell.py index c46df01..365cfd9 100644 --- a/src/sumfem/electrochemistry/cells/half_cell.py +++ b/src/sumfem/electrochemistry/cells/half_cell.py @@ -190,6 +190,9 @@ class RedoxHalfCell: parameter of a flow model to share it. viscosity: Dynamic viscosity of the electrolyte, given like the density. + temperature: The temperature of the electrolyte, an expression, e.g. + the field of a heat equation, or the default of its parameter, in + K, see :class:`~sumfem.electrochemistry.Electrolyte`. degree: Degree of the Lagrange elements of the unknowns. name: Prefix of the cell's own parameters, ``"."``, and of the names of its models, ``"_"``. The @@ -300,6 +303,7 @@ def __init__( proton_flux: bool = False, density: ExpressionOrDefault = 1399.8, viscosity: ExpressionOrDefault = 6.963e-3, + temperature: ExpressionOrDefault = 298.15, degree: int = 1, name: str = "cell", phi_s: Field | None = None, @@ -336,7 +340,7 @@ def __init__( self._create_fields(ions, reservoir, degree, phi_s, phi_l, fields) self._electrolyte = self._create_electrolyte( - ions, electrolyte_name, density, viscosity + ions, electrolyte_name, density, viscosity, temperature ) self._create_parameters(reservoir, fiber_distance) self._reaction = self._create_reaction(velocity, mass_transfer) @@ -392,6 +396,7 @@ def _create_electrolyte( name: str, density: ExpressionOrDefault, viscosity: ExpressionOrDefault, + temperature: ExpressionOrDefault, ) -> Electrolyte: """Create the electrolyte, with the fields of the transported ions. @@ -403,6 +408,8 @@ def _create_electrolyte( expression. viscosity: Dynamic viscosity of the electrolyte, a value or a shared expression. + temperature: Temperature of the electrolyte, a value or an + expression. Returns: The electrolyte. @@ -418,6 +425,7 @@ def _create_electrolyte( ], density=density, viscosity=viscosity, + temperature=temperature, ) def _create_parameters( diff --git a/src/sumfem/electrochemistry/cells/specs.py b/src/sumfem/electrochemistry/cells/specs.py index 78053e7..2d8023f 100644 --- a/src/sumfem/electrochemistry/cells/specs.py +++ b/src/sumfem/electrochemistry/cells/specs.py @@ -64,6 +64,8 @@ class HalfCellSpec: density: The density of the electrolyte, in kg/m^3 or with units. viscosity: The dynamic viscosity of the electrolyte, in Pa s or with units, which the flow and the mass transfer share. + temperature: The temperature of the electrolyte, in K, e.g. the field + of a heat equation. Raises: ValueError: If the side does not transport protons. @@ -80,6 +82,7 @@ class HalfCellSpec: kozeny_carman_constant: ExpressionOrDefault = None density: ExpressionOrDefault = 1399.8 viscosity: ExpressionOrDefault = 6.963e-3 + temperature: ExpressionOrDefault = 298.15 def __post_init__(self) -> None: if PROTONS not in self.reservoir: @@ -107,6 +110,8 @@ class MembraneSpec: m^2. hydraulic_permeability: The hydraulic permeability, in m^2. viscosity: The viscosity of the water in the membrane, in Pa s. + temperature: The temperature of the membrane, in K, e.g. the field of + a heat equation. """ fixed_charge: ExpressionOrDefault = None @@ -115,6 +120,7 @@ class MembraneSpec: electrokinetic_permeability: ExpressionOrDefault = None hydraulic_permeability: ExpressionOrDefault = None viscosity: ExpressionOrDefault = None + temperature: ExpressionOrDefault = 298.15 @dataclasses.dataclass(frozen=True) diff --git a/src/sumfem/electrochemistry/materials/electrolyte.py b/src/sumfem/electrochemistry/materials/electrolyte.py index 0a89c51..3f170f1 100644 --- a/src/sumfem/electrochemistry/materials/electrolyte.py +++ b/src/sumfem/electrochemistry/materials/electrolyte.py @@ -190,22 +190,7 @@ def __init__( names = [ion.name for ion in ions] if len(set(names)) != len(names): raise ValueError(f"Ion names must be unique, got {names}.") - balanced = [ion for ion in ions if ion.concentration is None] - if len(balanced) > 1: - raise ValueError( - "At most one ion can follow from electroneutrality, got " - f"{[ion.name for ion in balanced]}." - ) - if balanced and balanced[0].charge == 0: - raise ValueError( - f"Ion {balanced[0].name!r} has no charge, so electroneutrality " - "cannot give its concentration. Give it a field." - ) - if balanced and len(balanced) == len(ions): - raise ValueError( - "Electroneutrality needs at least one transported ion, got only " - f"{balanced[0].name!r}." - ) + balanced = _check_balance(ions) self._name = name self._F = Parameter( f"{name}.F", @@ -256,8 +241,8 @@ def __init__( ion.name: ion.concentration for ion in ions if ion.concentration is not None } charge = sum(self._ions[n].charge * c for n, c in self._concentrations.items()) - for ion in balanced: - self._concentrations[ion.name] = -charge / ion.charge + if balanced is not None: + self._concentrations[balanced.name] = -charge / balanced.charge @property def name(self) -> str: @@ -393,3 +378,39 @@ def current_density( return -self.conductivity(correction) * ufl.grad( potential ) - self.diffusion_current(correction) + + +def _check_balance(ions: typing.Sequence[Ion]) -> Ion | None: + """Check that electroneutrality can give the ion without a field, if any. + + Args: + ions: The ions of an electrolyte. + + Returns: + The one ion without a field, whose concentration follows from + electroneutrality, or ``None`` if every ion has one. + + Raises: + ValueError: If more than one ion has no field, the ion without one + has no charge, or it is the only ion. + """ + balanced = [ion for ion in ions if ion.concentration is None] + if len(balanced) > 1: + raise ValueError( + "At most one ion can follow from electroneutrality, got " + f"{[ion.name for ion in balanced]}." + ) + if not balanced: + return None + (ion,) = balanced + if ion.charge == 0: + raise ValueError( + f"Ion {ion.name!r} has no charge, so electroneutrality cannot give " + "its concentration. Give it a field." + ) + if len(ions) == 1: + raise ValueError( + "Electroneutrality needs at least one transported ion, got only " + f"{ion.name!r}." + ) + return ion diff --git a/src/sumfem/electrochemistry/models/membrane.py b/src/sumfem/electrochemistry/models/membrane.py index baf82c6..70155e5 100644 --- a/src/sumfem/electrochemistry/models/membrane.py +++ b/src/sumfem/electrochemistry/models/membrane.py @@ -260,7 +260,7 @@ def __init__( electrokinetic_permeability: ExpressionOrDefault = None, hydraulic_permeability: ExpressionOrDefault = None, viscosity: ExpressionOrDefault = None, - temperature: ufl.core.expr.Expr | float = 298.15, + temperature: ExpressionOrDefault = 298.15, ): domain = potential.domain self._potential = potential diff --git a/tests/electrochemistry/test_full_cell.py b/tests/electrochemistry/test_full_cell.py index e7393b8..cdbabbd 100644 --- a/tests/electrochemistry/test_full_cell.py +++ b/tests/electrochemistry/test_full_cell.py @@ -526,38 +526,53 @@ def test_any_flow_plugs_into_the_cell(): def test_inputs_can_be_fields_of_other_models(): - """The properties of the electrodes and the couples may be expressions of - the fields of other models: here a temperature, which sets the - conductivity of the fibres, an Arrhenius rate constant and a formal - potential that changes with it, and a compression, which sets the - porosity. Their parameters are then not created, the units are - consistent, and a cell at rest starts and stays at the Nernst voltage of - that temperature.""" + """The inputs of the cell may be expressions of the fields of other + models: here a temperature, the temperature of the electrolytes and the + membrane, which also sets the conductivity of the fibres, an Arrhenius + rate constant and a formal potential that changes with it, and a + compression, which sets the porosity. The expressions hold parameters + with units, so that the units stay consistent. Their own parameters are + then not created, and a cell at rest starts and stays at the Nernst + voltage of that temperature; the change of the formal potentials, the + same on both sides, cancels in it.""" domain = Domain("cell", "interval", length_scale=4e-3) T = Field("T", domain, ("Lagrange", 1), unit=syu.kelvin, scale=300.0) strain = Field("strain", domain, ("Lagrange", 1)) - reference = 298.15 * syu.kelvin - slope = 1e-3 * syu.volt / syu.kelvin - activation = 4e3 * syu.kelvin + + def create_parameter(name, default, unit=None): + """A parameter of the thermal model, with its unit.""" + return Parameter(f"thermal.{name}", domain, default=default, unit=unit) + + reference = create_parameter("reference", 298.15, syu.kelvin) + slope = create_parameter("slope", 1e-3, syu.volt / syu.kelvin) + activation = create_parameter("activation", 4e3, syu.kelvin) + conductivity = create_parameter("conductivity", 500.0, syu.siemens / syu.m) sides = [] - for spec in (NEGATIVE, POSITIVE): + for side, spec in (("negative", NEGATIVE), ("positive", POSITIVE)): + k0 = create_parameter( + f"{side}_rate", spec.couple.rate_constant, syu.m / syu.second + ) + E0 = create_parameter( + f"{side}_potential", spec.couple.formal_potential, syu.volt + ) couple = dataclasses.replace( spec.couple, - rate_constant=spec.couple.rate_constant - * ufl.exp(-activation * (1 / T - 1 / reference)), - formal_potential=spec.couple.formal_potential + slope * (T - reference), - reactants=None if spec is NEGATIVE else spec.couple.reactants, - products=None if spec is NEGATIVE else spec.couple.products, + rate_constant=k0 * ufl.exp(-activation * (1 / T - 1 / reference)), + formal_potential=E0 + slope * (T - reference), + reactants=None if side == "negative" else spec.couple.reactants, + products=None if side == "negative" else spec.couple.products, ) sides.append( dataclasses.replace( spec, couple=couple, porosity=0.68 * (1 - 0.1 * strain), - conductivity=500.0 * syu.siemens / syu.m * (T / reference), + conductivity=conductivity * (T / reference), + temperature=T, ) ) - cell = FullCell(domain, *sides, membrane=MEMBRANE) + membrane = dataclasses.replace(MEMBRANE, temperature=T) + cell = FullCell(domain, *sides, membrane=membrane) heat = Diffusion( "heat", T, @@ -568,6 +583,7 @@ def test_inputs_can_be_fields_of_other_models(): system = build_system(model, check_units=True) assert "negative_electrode.porosity" not in system.parameters assert "positive.rate_constant" not in system.parameters + assert "negative_electrolyte.T" not in system.parameters assert {T, strain} <= set(system.unknowns) disc = discretize( @@ -590,8 +606,7 @@ def test_inputs_can_be_fields_of_other_models(): time_scale=1.0, ) RT_F = 8.314462618 * temperature / 96485.33212 - shift = 2 * 1e-3 * (temperature - 298.15) - expected = 1.004 + 0.255 + shift + 2 * RT_F * np.log(27.0 / 1053.0) + expected = 1.004 + 0.255 + 2 * RT_F * np.log(27.0 / 1053.0) assert cell.voltage(disc) == pytest.approx(expected, rel=1e-9) simulation.run(num_steps=1) assert cell.voltage(disc) == pytest.approx(expected, rel=1e-8)