diff --git a/_config.yml b/_config.yml index 5164b49..42ae3e3 100644 --- a/_config.yml +++ b/_config.yml @@ -61,6 +61,13 @@ sphinx: - fmt: py # A tutorial or demo that raises fails the build, instead of a warning nb_execution_raise_on_error: true + # Pages that cite references cited elsewhere too list them as footnotes, + # {footcite} and {footbibliography}. MyST checks footnotes before + # sphinxcontrib-bibtex creates them, and warns about each; the book has + # no MyST footnotes of its own + suppress_warnings: ["myst.footnote"] + # No rule before the footnotes, which end a page under References + myst_footnote_transition: false autodoc_member_order: bysource # The page tables of contents list objects by name, without their module toc_object_entries_show_parents: hide diff --git a/_toc.yml b/_toc.yml index 6eefd2f..4b1e5d9 100644 --- a/_toc.yml +++ b/_toc.yml @@ -24,6 +24,9 @@ parts: - caption: "Electrochemistry" chapters: - file: "examples/electrochemistry/porous_electrode.py" + - file: "examples/electrochemistry/vanadium_half_cell_1d.py" + - file: "examples/electrochemistry/vanadium_half_cell_2d.py" + - file: "examples/electrochemistry/vanadium_half_cell_3d.py" - caption: "Python API" chapters: - file: "docs/api/index" diff --git a/docs/api/electrochemistry.md b/docs/api/electrochemistry.md index 0078203..5021441 100644 --- a/docs/api/electrochemistry.md +++ b/docs/api/electrochemistry.md @@ -49,3 +49,10 @@ .. autoclass:: MassTransferButlerVolmer :members: ``` + +## Half-cells + +```{eval-rst} +.. autoclass:: RedoxHalfCell + :members: +``` diff --git a/docs/api/symbolic.md b/docs/api/symbolic.md index 3008296..3de36ed 100644 --- a/docs/api/symbolic.md +++ b/docs/api/symbolic.md @@ -44,6 +44,8 @@ :members: .. autofunction:: as_expr + +.. autofunction:: value_at_defaults ``` ## Components @@ -121,7 +123,11 @@ .. autoclass:: Interpolation :members: +.. currentmodule:: sumfem.symbolic.transfer + .. autodata:: Transferable + +.. currentmodule:: sumfem.symbolic ``` ## Units and dimensional analysis @@ -140,12 +146,20 @@ .. autoexception:: UnitError +.. currentmodule:: sumfem.symbolic.constants + .. autodata:: BOLTZMANN +.. currentmodule:: sumfem.symbolic.constants + .. autodata:: FARADAY +.. currentmodule:: sumfem.symbolic.constants + .. autodata:: GAS_CONSTANT +.. currentmodule:: sumfem.symbolic + .. autofunction:: boltzmann_constant ``` diff --git a/examples/electrochemistry/vanadium_half_cell_1d.py b/examples/electrochemistry/vanadium_half_cell_1d.py new file mode 100644 index 0000000..9eddc8d --- /dev/null +++ b/examples/electrochemistry/vanadium_half_cell_1d.py @@ -0,0 +1,399 @@ +# --- +# 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 half-cell in one dimension +# +# This demo charges, rests and discharges the negative half-cell of a +# vanadium redox flow battery across its porous electrode, as the 1D example +# of [RFBniCS](https://github.com/Ah-Aghabarari/RFBniCS) +# {footcite}`Aghabarari2026`. It shows how to: +# +# - Build a half-cell of any redox couple and electrolyte, negative or +# positive, with +# {py:class}`RedoxHalfCell ` +# - Run it through a galvanostatic schedule of charge, rest and discharge, +# with a {py:class}`Simulation ` +# - Follow the voltage and the state of charge of the half-cell +# +# ## The mathematical model +# +# The porous electrode holds the solid, of potential $\phi_s$, and in its pores +# of porosity $\varepsilon$ the electrolyte, of potential $\phi_l$, with the +# ions V²⁺, V³⁺, H⁺ and HSO₄⁻ of concentrations $c_k$, and sulfate SO₄²⁻ by +# electroneutrality. At the fibres of specific area $a$, the reaction +# $\mathrm{V^{2+} \rightleftharpoons V^{3+} + e^-}$ carries the current +# density $i$ from the electrolyte to the solid. Charge is conserved in both +# phases, and each ion is stored, transported and produced: +# +# $$ +# \begin{aligned} +# -\nabla \cdot (\sigma_{\text{eff}} \nabla \phi_s) + a i &= 0, \\ +# \nabla \cdot i_l - a i &= 0, \\ +# \varepsilon \frac{\partial c_k}{\partial t} + \nabla \cdot N_k &= s_k, +# \\ +# \sum_k z_k c_k &= 0, +# \end{aligned} +# $$ (half-cell) +# +# with the effective conductivity of the solid $\sigma_{\text{eff}}$, the +# current density of the electrolyte $i_l$ and the Nernst-Planck fluxes +# $N_k$, whose diffusivities are corrected for the porosity. The reaction +# consumes V²⁺ and produces V³⁺, $s_{\mathrm{V^{2+}}} = -a i / F$ and +# $s_{\mathrm{V^{3+}}} = a i / F$, at the Butler-Volmer current density +# +# $$ +# \begin{aligned} +# i &= F k_0 \left( +# e^{\alpha_a F \eta / (R T)} c_{\mathrm{V^{2+}}} +# - e^{-\alpha_c F \eta / (R T)} c_{\mathrm{V^{3+}}} \right), \\ +# \eta &= \phi_s - \phi_l - E_0, +# \end{aligned} +# $$ (half-cell-kinetics) +# +# with the formal potential $E_0$. The current density $j$ of the cell, +# positive on discharge, enters the solid at the current collector and leaves +# the electrolyte at the membrane, where the electrolyte has the +# concentrations of the reservoir. The voltage of the half-cell is the mean +# solid potential at the current collector less the mean electrolyte +# potential at the membrane. The model is the macro-homogeneous +# porous-electrode model of RFBniCS {footcite}`Aghabarari2026` (Section II), +# after {footcite}`NewmanTiedemann1975`, and its parameters those of Table I +# there, several from {footcite}`MunozPerales2022`. +# +# We import the relevant packages: + +# + tags=["hide-input"] +from mpi4py import MPI + +import dolfinx.mesh +import matplotlib.pyplot as plt +import numpy as np +import sympy.physics.units as syu + +from sumfem import discrete, solvers, symbolic +from sumfem.electrochemistry import Ion, PorousElectrode, RedoxHalfCell + +# - + +# ## The electrode +# +# Across the electrode, from the current collector at $x = 0$ to the +# membrane at $x = W$, the cell is an interval. Its ends are located with a +# tolerance relative to the width, see {ref}`units-length-scale`: + +# + +width, num_cells = 2.5e-3, 250 +dtype = dolfinx.default_real_type +interval = dolfinx.mesh.create_interval( + MPI.COMM_WORLD, num_cells, [0.0, width], dtype=dtype +) +tolerance = 100 * float(np.finfo(dtype).eps) * width +marked = discrete.MeshData.from_locators( + interval, + { + 0: { + "CurrentCollector": lambda x: np.isclose(x[0], 0.0, atol=tolerance), + "Membrane": lambda x: np.isclose(x[0], width, atol=tolerance), + } + }, + boundary=True, +) +domain = symbolic.Domain.from_mesh("cell", interval, length_scale=width) +# - + +# The porous electrode, which fills the interval, is a carbon felt, a mat of +# carbon fibres, with its porosity, the conductivity of its fibres, their +# specific area and diameter, see +# {py:class}`PorousElectrode `. +# The fibres of commercial carbon felts, papers and cloths share a diameter +# of several micrometres {footcite}`Barzegari2026`; RFBniCS +# {footcite}`Aghabarari2026` uses a diameter far below that range, kept in a +# comment. The diameter and the +# constants of the permeability of the pores and of the inertial drag matter +# only with flow, and this demo has none: + +electrode = PorousElectrode( + "electrode", + domain, + porosity=0.94, + conductivity=3402.1 * syu.siemens / syu.m, + specific_area=2.3e4 / syu.m, + fiber_diameter=8.0 * syu.micrometer, # RFBniCS: 0.5 micrometre + kozeny_carman_constant=5.55, + pore_length=100.0 * syu.micrometer, +) + + +# ## The half-cell +# +# The ions of the electrolyte are given with their charges and diffusivities, +# in their units, which sumfem converts to SI, see {ref}`units-checking`. The +# vanadium ions take the positive part of their concentrations in the +# kinetics and the conductivity, as they may be depleted. The concentrations +# of the reservoir name the ions that are transported, each with a field and +# the equations of {eq}`half-cell`. Sulfate has none: electroneutrality, +# the last equation of {eq}`half-cell`, with the charges of the ions, +# +# $$ +# \begin{aligned} +# 2 c_{\mathrm{V^{2+}}} + 3 c_{\mathrm{V^{3+}}} + c_{\mathrm{H^+}} +# - c_{\mathrm{HSO_4^-}} - 2 c_{\mathrm{SO_4^{2-}}} &= 0, +# \end{aligned} +# $$ (electroneutrality) +# +# is solved for its concentration, which is then no unknown, but an +# expression of the others: +# +# $$ +# \begin{aligned} +# c_{\mathrm{SO_4^{2-}}} &= \frac{1}{2} \left( 2 c_{\mathrm{V^{2+}}} +# + 3 c_{\mathrm{V^{3+}}} + c_{\mathrm{H^+}} +# - c_{\mathrm{HSO_4^-}} \right). +# \end{aligned} +# $$ (sulfate) +# +# It enters the conductivity and the diffusion current of the electrolyte +# like the concentrations of the other ions. The half-cell checks that the +# reservoir lacks exactly one ion, and that this ion is charged: + +# + +diffusivity = syu.m**2 / syu.second +concentration = syu.mol / syu.m**3 +ions = ( + Ion("V2", 2, 8.10e-12 * diffusivity, positive=True), + Ion("V3", 3, 1.65e-11 * diffusivity, positive=True), + Ion("H", 1, 9.312e-9 * diffusivity), + Ion("HSO4", -1, 1.33e-9 * diffusivity), + Ion("SO4", -2, 1.065e-9 * diffusivity), +) +reservoir = { + "V2": 300.0 * concentration, + "V3": 300.0 * concentration, + "H": 1000.0 * concentration, + "HSO4": 1000.0 * concentration, +} +# - + +# {py:class}`RedoxHalfCell ` sums the +# models of {eq}`half-cell` for a redox couple. The reaction +# $\mathrm{V^{2+} \rightleftharpoons V^{3+} + e^-}$ of {eq}`half-cell-kinetics` +# is given by three arguments: the ion `reduced`, V²⁺, the ion `oxidized`, +# V³⁺, and the number of `electrons` the reaction transfers, which the +# charges of the ions must match. The kinetics take the standard rate +# constant $k_0$ and the formal potential $E_0$ of the couple. +# The other ions take no part in the reaction; they carry current and keep +# the electrolyte neutral. Discharge oxidizes V²⁺ in this negative +# half-cell; the positive half-cell of the battery, with the couple +# VO²⁺/VO₂⁺, reduces its couple on discharge, which reverses the current. The +# concentrations of the reservoir are Dirichlet conditions at the membrane: + +cell = RedoxHalfCell( + domain, + ions, + reservoir, + reduced="V2", + oxidized="V3", + electrons=1, + electrode=electrode, + oxidation_on_discharge=True, + rate_constant=2.5e-7 * syu.m / syu.second, + formal_potential=-0.36 * syu.volt, +) +inlet = cell.inlet_conditions("Membrane") + +# The discretization bounds the quadrature degree of the nonlinear kinetics, +# as RFBniCS does. The initial state is at rest: the concentrations of the +# reservoir, the solid at zero and the electrolyte at minus the formal +# potential: + +disc = discrete.discretize( + [cell.model, *inlet], + discrete.DomainMaps({domain: marked}), + max_quadrature_degree=4, +) +cell.initialize(disc) + +# ## The charging and discharging schedule +# +# The cell is charged, rested and discharged, each with time steps that +# start short, while the cell responds fast, and then grow. Each operation +# gives its mode, the magnitude of the current density in A/m², and its time +# steps, as pairs of a number of steps and their size in seconds. The +# current of the cell is positive on discharge, so the mode gives its sign: + +# + +schedule = ( + ("charge", 0.5, ((20, 0.1), (10, 1.0), (10, 10.0))), + ("rest", 0.0, ((5, 0.1), (10, 1.0), (20, 2.0))), + ("discharge", 0.5, ((20, 0.1), (10, 1.0), (10, 10.0))), +) +sign = {"charge": -1.0, "rest": 0.0, "discharge": 1.0} +# - + +# ## Dimensionless groups +# +# The model is semi-discrete, so its units are checked on its system alone, +# see {ref}`units-checking`, with the length of the whole schedule as the +# reference time of the time derivatives. With the width of the electrode, +# the scales of the concentration of V²⁺ and of the solid potential, +# Faraday's constant, the temperature and the diffusivity of V²⁺ as the +# repeating quantities, {py:func}`pi_groups ` +# gives one group for every other quantity of the model: + +# + +duration = sum(count * dt for *_, steps in schedule for count, dt in steps) +semi_discrete = symbolic.build_system(cell.model, check_units=True, time_scale=duration) +electrolyte = cell.electrolyte +repeating = [ + domain, + cell.concentrations["V2"], + cell.phi_s, + electrolyte.F, + electrolyte.T, + electrolyte.diffusivity("V2"), +] +groups = symbolic.pi_groups(semi_discrete, time_scale=duration, repeating=repeating) +for group in groups: + print(f"{group.expression} = {group.value:.3g}") + + +def value_of(name: str) -> float: + """The value of the one group of the quantity of a name.""" + (value,) = [g.value for g in groups if name in str(g.expression)] + return value + + +# - + +# Most groups are ratios of like quantities: of the diffusivities and the +# concentrations of the ions, the porosity and the transfer coefficients. +# Four describe the behaviour of the cell: +# +# - The group of the reference time is the [mass Fourier +# number](https://dimless.org/docs/numbers/mass-fourier/) $D t / W^2$ of +# V²⁺ over the whole schedule. It is far below one: the ions diffuse only +# a small part of the way across the electrode, so the electrode exchanges +# little with the reservoir at the membrane. +# - The groups of the rate constant and the specific area multiply to the +# Damköhler number $a k_0 W^2 / D$ {footcite}`Fogler2020`, the rate of the +# reaction over that of diffusion across the electrode. Its square root is +# the Thiele modulus {footcite}`Thiele1939`, the width of the electrode over +# the depth a reactant supplied from the membrane would reach before it +# reacts. Both are far above one. +# - The group of the current of the cell is the current over the limiting +# current of diffusion across the electrode, $F D c / W$ +# {footcite}`NewmanBalsara2021`. Above one, the cell runs on the vanadium +# stored in its pores, not on what diffuses in from the reservoir. The +# current has no default, so the analysis takes a unit current density, +# and the group is scaled with the current of the schedule. +# - The groups of the conductivity of the solid and of the gas constant +# multiply to $\sigma R T / (F^2 D c)$, the conductivity of the solid over +# that V²⁺ alone would give the electrolyte. +# +# Neither the Damköhler number nor the Thiele modulus is yet on +# [dimless.org](https://dimless.org/docs/numbers/): + +# + +fourier = value_of("t_ref") +damkohler = value_of("rate_constant") * value_of("specific_area") +largest_current = max(current for _, current, _ in schedule) +current_ratio = value_of("cell.current") * largest_current +conductivity_ratio = value_of("electrode.conductivity") * value_of("electrolyte.R") +print(f"Mass Fourier number of the schedule: {fourier:.3g}") +print(f"Damkoehler number: {damkohler:.3g}, Thiele modulus: {np.sqrt(damkohler):.3g}") +print(f"Current over the limiting current of diffusion: {current_ratio:.3g}") +print(f"Conductivity of the solid over that of V2+: {conductivity_ratio:.3g}") +assert fourier < 1e-2 and damkohler > 1e2 and current_ratio > 1 +# - + +# A {py:class}`Simulation ` takes the steps with +# backward Euler, checks the units, and solves every time step in +# dimensionless form, see {ref}`units-scaling`. The steps of the schedule +# change in size, so the reference time of the time derivatives is the +# largest of them: a shorter step only makes Newton's absolute tolerance +# stricter, see {doc}`../cooling_plate`: + +# + +largest_step = max(dt for *_, steps in schedule for _, dt in steps) +clock = symbolic.Clock(domain, dt=schedule[0][2][0][1]) +simulation = solvers.Simulation( + disc, + cell.model, + inlet, + clock=clock, + scheme=solvers.BDF(1), + petsc_options={"snes_linesearch_type": "bt", "snes_max_it": 100}, + check_units=True, + factorize_equations=True, + time_scale=largest_step, +) +# - + +# ## Charge, rest and discharge +# +# Each operation sets the current of the cell, and each part of it the size +# of the time step, before the simulation takes its steps. After every step, +# a callback records the time, the mode, the voltage of the half-cell and its +# state of charge, the share of V²⁺ in the vanadium of the electrode: + +# + +records: list[tuple[float, str, float, float]] = [] +for mode, current, steps in schedule: + disc.set({cell.current: sign[mode] * current}) + for count, dt in steps: + disc.set({clock.dt: dt}) + simulation.run( + num_steps=count, + callback=lambda t, step, mode=mode: records.append( + (t, mode, cell.voltage(disc), cell.state_of_charge(disc)) + ), + ) +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]) +charge = np.array([r[3] for r in records]) +# - + +# Charging reduces V³⁺ to V²⁺, so the state of charge grows, and discharging +# turns it back. At rest only the exchange with the reservoir at the membrane +# changes the electrode, slowly. The voltage is lower on charge than on +# discharge, by the overpotentials that drive the reaction either way: + +# + +steps = {mode: np.diff(charge[modes == mode]) for mode in ("charge", "discharge")} +print(f"State of charge after charge: {charge[modes == 'charge'][-1]:.6f}") +print(f"State of charge at the end: {charge[-1]:.6f}") +assert np.all(steps["charge"] > 0) and np.all(steps["discharge"] < 0) +rest = charge[modes == "rest"] +assert np.abs(np.diff(rest)).max() < np.abs(steps["charge"]).min() +assert voltage[modes == "charge"].mean() < voltage[modes == "discharge"].mean() +assert np.isfinite(voltage).all() +# - + +# The voltage and the state of charge through the schedule: + +fig, (top, bottom) = plt.subplots(2, 1, sharex=True) +for mode in ("charge", "rest", "discharge"): + top.plot(times[modes == mode], 1e3 * voltage[modes == mode], ".", label=mode) + bottom.plot(times[modes == mode], charge[modes == mode], ".", label=mode) +top.set_ylabel("Voltage (mV)") +bottom.set_ylabel("State of charge") +bottom.set_xlabel("Time (s)") +top.legend() +plt.show() + +# ## References +# +# ```{footbibliography} +# ``` diff --git a/examples/electrochemistry/vanadium_half_cell_2d.py b/examples/electrochemistry/vanadium_half_cell_2d.py new file mode 100644 index 0000000..95ae2bd --- /dev/null +++ b/examples/electrochemistry/vanadium_half_cell_2d.py @@ -0,0 +1,421 @@ +# --- +# 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 half-cell in two dimensions, with flow +# +# This demo extends {doc}`vanadium_half_cell_1d` with the flow of the +# electrolyte through the porous electrode, as the 2D example of +# [RFBniCS](https://github.com/Ah-Aghabarari/RFBniCS) {footcite}`Aghabarari2026`. +# It shows how to: +# +# - Solve Darcy flow through the electrode with +# {py:class}`Darcy `, and check it against its exact +# solution +# - Couple the half-cell to the flow, which carries the ions in at the inlet +# and out at the outlet, and limits the kinetics by mass transfer to the +# fibres +# - Name the dimensionless numbers of the flow from the expressions of the +# model +# - Run the half-cell through a charging and discharging schedule, and plot +# its local state of charge +# +# ## The model +# +# The electrode spans the width $W$ from the current collector at $x = 0$ to +# the membrane at $x = W$, and the height $H$ from the inlet at $y = 0$ to +# the outlet at $y = H$. The electrolyte enters at the inlet with the +# superficial velocity $v_{\text{in}}$ and leaves at the outlet, held at the +# pressure $p_{\text{out}}$. In the pores, Darcy's law gives the +# superficial velocity $\mathbf{u}$ from the pressure $p$, with the +# permeability $K$ of the electrode and the viscosity $\mu$ of the +# electrolyte: +# +# $$ +# \begin{aligned} +# \mathbf{u} &= -\frac{K}{\mu} \nabla p, \\ +# \nabla \cdot \mathbf{u} &= 0, \\ +# -\mathbf{u} \cdot \mathbf{n} &= v_{\text{in}} \quad \text{at the inlet}, \\ +# p &= p_{\text{out}} \quad \text{at the outlet}, +# \end{aligned} +# $$ (darcy) +# +# and no flow through the other sides. Darcy's law is the balance of the +# pressure gradient and the Darcy drag in the volume-averaged Navier-Stokes +# equations, whose inertia, viscous diffusion and Forchheimer drag are +# negligible in the electrode {footcite}`Aghabarari2026` (Section II B 5). The +# permeability follows from the porosity and the fibre diameter by the +# Kozeny-Carman relation, see +# {py:func}`kozeny_carman `. +# +# The half-cell is that of {eq}`half-cell` in {doc}`vanadium_half_cell_1d`, +# with three changes. The Nernst-Planck fluxes carry the ions with +# $\mathbf{u}$. The concentrations of the reservoir hold at the inlet +# instead of the membrane, the ions leave with the flow at the outlet, and +# the membrane passes the current of the cell as protons. The reaction takes +# the concentrations at the surface of the fibres, which mass transfer from +# the pores, at the rate $k_m (c - c_s)$, ties to the concentrations $c$ of +# the pores, with a mass-transfer coefficient $k_m$ from the correlation of +# {footcite}`Krishnamurthy2011`, see +# {py:class}`~sumfem.electrochemistry.MassTransferButlerVolmer`. These are the +# boundary conditions and kinetics of {footcite}`Aghabarari2026` (Sections II B 4 +# and II C). +# The flow is steady and does not depend on the electrochemistry, so it is +# solved once, first. +# +# We import the relevant packages: + +# + tags=["hide-input"] +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 +import ufl + +from sumfem import discrete, solvers, symbolic +from sumfem.electrochemistry import Ion, PorousElectrode, RedoxHalfCell +from sumfem.models import BoundaryFlux, Darcy +from sumfem.plotting import gather_grid + +# - + +# ## The electrode +# +# The cell is a rectangle of quadrilaterals, with its sides marked. They are +# located with a tolerance relative to the height, see +# {ref}`units-length-scale`: + +# + +width, height, cells = 2.5e-3, 1e-2, (25, 100) +dtype = dolfinx.default_real_type +rectangle = dolfinx.mesh.create_rectangle( + MPI.COMM_WORLD, + np.array([[0.0, 0.0], [width, height]]), + list(cells), + cell_type=dolfinx.mesh.CellType.quadrilateral, + dtype=dtype, +) +tolerance = 100 * float(np.finfo(dtype).eps) * height +marked = discrete.MeshData.from_locators( + rectangle, + { + 1: { + "CurrentCollector": lambda x: np.isclose(x[0], 0.0, atol=tolerance), + "Membrane": lambda x: np.isclose(x[0], width, atol=tolerance), + "Inlet": lambda x: np.isclose(x[1], 0.0, atol=tolerance), + "Outlet": lambda x: np.isclose(x[1], height, atol=tolerance), + } + }, + boundary=True, +) +domain = symbolic.Domain.from_mesh("cell", rectangle, length_scale=width) +# - + +# The porous electrode is the carbon felt of {doc}`vanadium_half_cell_1d`, +# whose fibres of several micrometres {footcite}`Barzegari2026` now matter: their +# diameter sets the permeability and the mass transfer. RFBniCS +# {footcite}`Aghabarari2026` uses a diameter far below that range, kept in a +# comment: + +electrode = PorousElectrode( + "electrode", + domain, + porosity=0.94, + conductivity=3402.1 * syu.siemens / syu.m, + specific_area=2.3e4 / syu.m, + fiber_diameter=8.0 * syu.micrometer, # RFBniCS: 0.5 micrometre + kozeny_carman_constant=5.55, + pore_length=100.0 * syu.micrometer, +) + +# ## The flow +# +# The density and the viscosity of the electrolyte are parameters that the +# flow and the half-cell share, as are the inlet velocity and the outlet +# pressure. The flow depends on the gradient of the pressure alone, so the +# pressure is the gauge pressure, zero at the outlet: an offset as large as +# the atmospheric pressure would cost the small pressure drop its digits. + +# + +viscosity = symbolic.Parameter( + "electrolyte.viscosity", + domain, + default=6.963e-3 * syu.pascal * syu.second, + unit=syu.pascal * syu.second, + doc="Dynamic viscosity of the electrolyte", +) +density = symbolic.Parameter( + "electrolyte.density", + domain, + default=1399.8 * syu.kg / syu.m**3, + unit=syu.kg / syu.m**3, + doc="Density of the electrolyte", +) +inlet_velocity = symbolic.Parameter( + "flow.inlet_velocity", + domain, + default=1e-7 * syu.m / syu.second, + unit=syu.m / syu.second, + doc="Superficial velocity at the inlet", +) +outlet_pressure = symbolic.Parameter( + "flow.outlet_pressure", + domain, + default=0.0 * syu.pascal, + unit=syu.pascal, + doc="Gauge pressure at the outlet", +) +# - + +# The pressure drop $\mu v_{\text{in}} H / K$ of the uniform flow through the +# electrode is the scale of the pressure. +# {py:func}`value_at_defaults ` evaluates +# it, an expression of parameters, at their defaults: + +# + +pressure_drop = symbolic.value_at_defaults( + viscosity * inlet_velocity * height / electrode.permeability() +) +p = symbolic.Field("p", domain, ("Lagrange", 1), unit=syu.pascal, scale=pressure_drop) +# - + +# {py:class}`Darcy ` is the flow of {eq}`darcy`, and a +# {py:class}`BoundaryFlux ` its inflow. Its system +# is checked for units and solved in dimensionless form, see +# {ref}`units-scaling`: + +darcy = Darcy("darcy", p, electrode.permeability(), viscosity) +flow = darcy + BoundaryFlux("inflow", p, inlet_velocity, ("Inlet",)) +flow_system = symbolic.build_system(flow, check_units=True, factorize_equations=True) +outlet = symbolic.DirichletBC(p, "Outlet", outlet_pressure) + +# ## The half-cell +# +# The ions and the reservoir are those of {doc}`vanadium_half_cell_1d`, with +# sulfate from electroneutrality: + +# + +diffusivity = syu.m**2 / syu.second +concentration = syu.mol / syu.m**3 +ions = ( + Ion("V2", 2, 8.10e-12 * diffusivity, positive=True), + Ion("V3", 3, 1.65e-11 * diffusivity, positive=True), + Ion("H", 1, 9.312e-9 * diffusivity), + Ion("HSO4", -1, 1.33e-9 * diffusivity), + Ion("SO4", -2, 1.065e-9 * diffusivity), +) +reservoir = { + "V2": 300.0 * concentration, + "V3": 300.0 * concentration, + "H": 1000.0 * concentration, + "HSO4": 1000.0 * concentration, +} +# - + +# The half-cell is the negative one of {doc}`vanadium_half_cell_1d`. The +# superficial velocity of the flow carries the ions and drives the mass +# transfer to the fibres, which a velocity switches on. The flow carries the +# ions out at the `outlet`, the membrane passes the current as protons with +# `proton_flux`, and the reservoir concentrations hold at the inlet: + +cell = RedoxHalfCell( + domain, + ions, + reservoir, + reduced="V2", + oxidized="V3", + electrons=1, + electrode=electrode, + oxidation_on_discharge=True, + rate_constant=2.5e-7 * syu.m / syu.second, + formal_potential=-0.36 * syu.volt, + velocity=darcy.superficial_velocity, + outlet=("Outlet",), + proton_flux=True, + density=density, + viscosity=viscosity, +) +inlet = cell.inlet_conditions("Inlet") + +# The flow and the half-cell are discretized together, so they share their +# parameters, with the quadrature degree of the nonlinear kinetics bounded as +# in {doc}`vanadium_half_cell_1d`. A {py:class}`Solver ` +# solves the flow: + +disc = discrete.discretize( + [flow_system, outlet, cell.model, *inlet], + discrete.DomainMaps({domain: marked}), + max_quadrature_degree=4, +) +solvers.Solver(disc[flow_system], [outlet]).solve() + +# With a uniform permeability, the flow is uniform, and the pressure falls +# linearly from the inlet to the outlet by $\mu v_{\text{in}} H / K$. Linear +# elements hold it exactly, so the pressure drop and the outflow, the flux of +# the superficial velocity through the outlet, match it to round-off: + +# + +n = ufl.FacetNormal(domain) +ds_inlet, ds_outlet = domain.ds("Inlet"), domain.ds("Outlet") +area = disc.assemble_scalar(1.0 * domain.dx()) +v_in = float(disc.value(inlet_velocity)) +drop = disc.assemble_scalar((p - outlet_pressure) * ds_inlet) / width +outflow = disc.assemble_scalar(ufl.dot(darcy.superficial_velocity, n) * ds_outlet) +print(f"Pressure drop: {drop:.6g} Pa, exact {pressure_drop:.6g} Pa") +print(f"Outflow over inflow: {outflow / (v_in * width):.12f}") +assert np.isclose(drop, pressure_drop, rtol=1e-8) +assert np.isclose(outflow, v_in * width, rtol=1e-8) +# - + +# ## Dimensionless numbers of the flow +# +# The flow past the fibres has the [Reynolds +# number](https://dimless.org/docs/numbers/reynolds/) $\rho |\mathbf{u}| d / +# \mu$ of the fibre diameter $d$, which the kinetics compute, and sets the +# [Sherwood number](https://dimless.org/docs/numbers/sherwood/) $k_m d / D$ +# of the mass transfer by the correlation of +# {py:class}`~sumfem.electrochemistry.MassTransferButlerVolmer`. +# Along the electrode, the [mass Péclet +# number](https://dimless.org/docs/numbers/mass-peclet/) $v_{\text{in}} H / D$ +# of V²⁺ compares the transport by the flow with diffusion. Their means over +# the electrode: + +# + +reaction = cell.reaction +reynolds = disc.assemble_scalar(reaction.reynolds_number() * domain.dx()) / area +d_V2 = float(disc.value(cell.electrolyte.diffusivity("V2"))) +d_fibre = float(disc.value(electrode.fiber_diameter)) +(k_V2, _) = reaction.mass_transfer_coefficients() +sherwood = disc.assemble_scalar(k_V2 * domain.dx()) / area * d_fibre / d_V2 +peclet = v_in * height / d_V2 +print(f"Reynolds number of the fibres: {reynolds:.3g}") +print(f"Sherwood number of V2+: {sherwood:.3g}") +print(f"Mass Peclet number of V2+ along the electrode: {peclet:.3g}") +assert reynolds < 1 and peclet > 1 +# - + +# The flow is creeping, so Darcy's law holds, yet it carries the ions along +# the electrode faster than they diffuse. +# +# ## Charge, rest and discharge +# +# The schedule of {doc}`vanadium_half_cell_1d`, run by a +# {py:class}`Simulation ` in parts, one for each +# current and step size. The simulation's model is the half-cell alone: the +# flow does not depend on the electrochemistry and its inputs are fixed, so +# every step takes the pressure solved above as known, and the velocity from +# it, without solving the flow again. The initial state is at rest, and the +# reference time of the time derivatives is the largest step: + +# + +schedule = ( + ("charge", 0.5, ((20, 0.1), (10, 1.0), (10, 10.0))), + ("rest", 0.0, ((5, 0.1), (5, 1.0), (10, 10.0))), + ("discharge", 0.5, ((20, 0.1), (10, 1.0), (10, 10.0))), +) +sign = {"charge": -1.0, "rest": 0.0, "discharge": 1.0} +largest_step = max(dt for *_, steps in schedule for _, dt in steps) +cell.initialize(disc) +clock = symbolic.Clock(domain, dt=schedule[0][2][0][1]) +simulation = solvers.Simulation( + disc, + cell.model, + inlet, + clock=clock, + scheme=solvers.BDF(1), + petsc_options={"snes_linesearch_type": "bt", "snes_max_it": 100}, + check_units=True, + factorize_equations=True, + time_scale=largest_step, +) +# - + +# After every step, a callback records the time, the mode, the voltage and +# the state of charge. At the end of the charge, the local state of charge, +# the share of V²⁺ in the vanadium at each point, is kept for a plot: + +# + +local = cell.local_state_of_charge() +space = disc.function_space(cell.concentrations["V2"]) +charged = dolfinx.fem.Function(space, name="soc") +records: list[tuple[float, str, float, float]] = [] +for mode, current, steps in schedule: + disc.set({cell.current: sign[mode] * current}) + for count, dt in steps: + disc.set({clock.dt: dt}) + simulation.run( + num_steps=count, + callback=lambda t, step, mode=mode: records.append( + (t, mode, cell.voltage(disc), cell.state_of_charge(disc)) + ), + ) + if mode == "charge": + charged.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]) +charge = np.array([r[3] for r in records]) +# - + +# Charging raises the state of charge and discharging lowers it. The flow +# now renews the electrolyte, so at rest it carries the charged electrolyte +# out and brings the reservoir's in: the state of charge falls back towards +# that of the reservoir: + +# + +steps = {mode: np.diff(charge[modes == mode]) for mode in ("charge", "discharge")} +rest = charge[modes == "rest"] +print(f"State of charge after charge: {charge[modes == 'charge'][-1]:.6f}") +print(f"State of charge after rest: {rest[-1]:.6f}") +print(f"State of charge at the end: {charge[-1]:.6f}") +assert np.all(steps["charge"] > 0) and np.all(steps["discharge"] < 0) +assert np.all(np.diff(rest) < 0) +assert voltage[modes == "charge"].mean() < voltage[modes == "discharge"].mean() +assert np.isfinite(voltage).all() +# - + +# The voltage and the state of charge through the schedule: + +fig, (top, bottom) = plt.subplots(2, 1, sharex=True) +for mode in ("charge", "rest", "discharge"): + top.plot(times[modes == mode], 1e3 * voltage[modes == mode], ".", label=mode) + bottom.plot(times[modes == mode], charge[modes == mode], ".", label=mode) +top.set_ylabel("Voltage (mV)") +bottom.set_ylabel("State of charge") +bottom.set_xlabel("Time (s)") +top.legend() +plt.show() + +# The local state of charge at the end of the charge, gathered from all +# processes with {py:func}`gather_grid `. The +# electrode is drawn in millimetres, widened to show its width: + +grid = gather_grid(charged) +if MPI.COMM_WORLD.rank == 0: + grid.points[:, 0] *= 4e3 + 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_3d.py b/examples/electrochemistry/vanadium_half_cell_3d.py new file mode 100644 index 0000000..ef79e3c --- /dev/null +++ b/examples/electrochemistry/vanadium_half_cell_3d.py @@ -0,0 +1,488 @@ +# --- +# 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 half-cell in three dimensions, with flow channels +# +# This demo runs the half-cell of {doc}`vanadium_half_cell_2d` in a +# flow-through electrode with two channels, as the 3D example of +# [RFBniCS](https://github.com/Ah-Aghabarari/RFBniCS) {footcite}`Aghabarari2026`. +# It shows how to: +# +# - Build a mesh with physical groups in [Gmsh](https://gmsh.info), in +# memory, and mark it with them +# - Solve the flow of the electrolyte through the electrode and along the +# channels, Brinkman-Forchheimer in the electrode and Navier-Stokes in the +# channels, from a creeping-flow start +# - Couple the half-cell to the flow in the electrode and the channels, with +# the current of the cell balanced between the membrane and the smaller +# current collector +# +# ## The model +# +# The electrode is a thin plate. Two channels lie against its front: the +# electrolyte enters through the inlet channel, crosses the electrode along +# its height and leaves through the outlet channel. The front of the +# electrode that the channels leave free is the current collector, and its +# back the membrane. In the electrode, of porosity $\varepsilon$, the +# superficial velocity $\mathbf{u}$ and the pressure $p$ follow +# Brinkman-Forchheimer flow, with the Darcy drag of the permeability $K$ +# and the inertial drag of the Forchheimer coefficient $\beta$: +# +# $$ +# \begin{aligned} +# \frac{\rho}{\varepsilon^2} (\mathbf{u} \cdot \nabla) \mathbf{u} +# - \nabla \cdot \left( \frac{2 \mu}{\varepsilon} +# \boldsymbol{\epsilon}(\mathbf{u}) \right) + \nabla p +# + \frac{\mu}{K} \mathbf{u} + \rho \beta |\mathbf{u}| \mathbf{u} +# &= 0, \\ +# \nabla \cdot \mathbf{u} &= 0, +# \end{aligned} +# $$ (brinkman-forchheimer) +# +# with the symmetric gradient $\boldsymbol{\epsilon}$, see +# {py:class}`BrinkmanForchheimer `. In +# the channels, $\varepsilon = 1$ and the drag vanishes: the flow is +# Navier-Stokes flow. The velocity is given at the inlet, vanishes on the +# walls, the current collector and the membrane, and the outlet is free of +# stress. The Forchheimer coefficient is Ergun's, an estimate for fibres, +# see {py:func}`ergun_coefficient `. +# +# The half-cell is that of {doc}`vanadium_half_cell_2d`, with the +# electrolyte also in the channels, where it is stored and transported +# without porosity and without reaction. +# +# We import the relevant packages: + +# + tags=["hide-input"] +from mpi4py import MPI + +import dolfinx.fem +import dolfinx.io.gmsh +import gmsh +import matplotlib.pyplot as plt +import numpy as np +import pyvista +import sympy.physics.units as syu +import ufl + +from sumfem import discrete, solvers, symbolic +from sumfem.electrochemistry import Ion, PorousElectrode, RedoxHalfCell +from sumfem.models import Brinkman, BrinkmanForchheimer, GradDiv, NavierStokes, Stokes +from sumfem.plotting import gather_grid + +# - + +# ## The mesh +# +# The geometry is the flow-through flow field of RfbFoam +# {footcite}`Barzegari2026`, as RFBniCS {footcite}`Aghabarari2026` (Table III) builds +# it, from whose example (MIT licence) the code below is adapted, in +# millimetres: a plate electrode, and two channels with rounded ends that +# open at its front. A channel is a box with a cylinder at each end: + + +# + +def rounded_channel( + x_start: float, x_depth: float, y_center: float, z_start: float, z_end: float +) -> int: + """A channel of unit width along z with rounded ends, open along x.""" + radius = 0.5 + box = gmsh.model.occ.addBox( + x_start, y_center - radius, z_start, x_depth, 2 * radius, z_end - z_start + ) + caps = [ + gmsh.model.occ.addCylinder(x_start, y_center, z, x_depth, 0.0, 0.0, radius) + for z in (z_start, z_end) + ] + fused, _ = gmsh.model.occ.fuse([(3, box)], [(3, cap) for cap in caps]) + return fused[0][1] + + +# - + +# The electrode and the channels are fragmented into one conforming mesh, +# with the electrode and each channel a physical group of cells. The +# boundary surfaces are told apart by where they lie: the open ends of the +# channels are the inlet and the outlet, the free front of the electrode the +# current collector, its back the membrane, and the rest walls: + + +# + +def flow_through_electrode(sizes: tuple[float, float]) -> None: + """The Gmsh model of the cell, meshed with cells of the given sizes.""" + height, width, thickness, depth = 15.0, 17.0, 0.42, 1.0 + z_start, z_end, y_inlet, y_outlet = 1.5, 13.5, 1.0, 16.0 + electrode = gmsh.model.occ.addBox(depth, 0.0, 0.0, thickness, width, height) + channels = [ + rounded_channel(0.0, depth, y, z_start, z_end) for y in (y_inlet, y_outlet) + ] + gmsh.model.occ.synchronize() + _, fragments = gmsh.model.occ.fragment([(3, electrode)], [(3, c) for c in channels]) + gmsh.model.occ.synchronize() + for name, entities in zip( + ("Electrode", "InletChannel", "OutletChannel"), fragments, strict=True + ): + tags = [tag for dim, tag in entities if dim == 3] + gmsh.model.addPhysicalGroup(3, tags, name=name) + surfaces: dict[str, list[int]] = { + "Inlet": [], + "Outlet": [], + "CurrentCollector": [], + "Membrane": [], + "Walls": [], + } + tol = 1e-4 + for dim, tag in gmsh.model.getEntities(2): + cells, _ = gmsh.model.getAdjacencies(dim, tag) + if len(cells) == 2: + # An interface between the electrode and a channel + continue + xmin, ymin, _, xmax, ymax, _ = gmsh.model.getBoundingBox(dim, tag) + if np.isclose(xmin, 0.0, atol=tol) and np.isclose(xmax, 0.0, atol=tol): + y = 0.5 * (ymin + ymax) + near_inlet = abs(y - y_inlet) < abs(y - y_outlet) + surfaces["Inlet" if near_inlet else "Outlet"].append(tag) + elif np.isclose(xmin, depth, atol=tol) and np.isclose(xmax, depth, atol=tol): + surfaces["CurrentCollector"].append(tag) + elif np.isclose(xmin, depth + thickness, atol=tol) and np.isclose( + xmax, depth + thickness, atol=tol + ): + surfaces["Membrane"].append(tag) + else: + surfaces["Walls"].append(tag) + for name, tags in surfaces.items(): + gmsh.model.addPhysicalGroup(2, tags, name=name) + gmsh.option.setNumber("Mesh.MeshSizeMin", sizes[0]) + gmsh.option.setNumber("Mesh.MeshSizeMax", sizes[1]) + gmsh.model.mesh.generate(3) + + +# - + +# The first process builds the model, and +# {py:func}`model_to_mesh ` distributes it +# without a file. The mesh is coarse, a few cells across the electrode, to +# keep the demo short; RFBniCS meshes far finer. Gmsh scales coordinates only +# when it writes a file, so they are scaled from millimetres to metres here. +# {py:meth}`MeshData.from_gmsh ` names +# the markers after the physical groups: + +# + +gmsh.initialize() +gmsh.option.setNumber("General.Terminal", 0) +if MPI.COMM_WORLD.rank == 0: + flow_through_electrode(sizes=(0.2, 0.5)) +data = dolfinx.io.gmsh.model_to_mesh(gmsh.model, MPI.COMM_WORLD, 0, gdim=3) +gmsh.finalize() +data.mesh.geometry.x[:] *= 1e-3 +marked = discrete.MeshData.from_gmsh(data) +thickness = 0.42e-3 +domain = symbolic.Domain.from_mesh("cell", data.mesh, length_scale=thickness) +num_cells = data.mesh.topology.index_map(3).size_global +print(f"Cells: {num_cells}") +# - + +# ## The electrode and the electrolyte +# +# The electrode of {doc}`vanadium_half_cell_2d` fills the cells marked +# `Electrode`, the channels the others. The density and the viscosity of the +# electrolyte are shared by the flow and the half-cell: + +# + +electrode = PorousElectrode( + "electrode", + domain, + markers=("Electrode",), + porosity=0.94, + conductivity=3402.1 * syu.siemens / syu.m, + specific_area=2.3e4 / syu.m, + fiber_diameter=8.0 * syu.micrometer, # RFBniCS: 0.5 micrometre + kozeny_carman_constant=5.55, + pore_length=100.0 * syu.micrometer, +) +channels = ("InletChannel", "OutletChannel") +viscosity = symbolic.Parameter( + "electrolyte.viscosity", + domain, + default=6.963e-3 * syu.pascal * syu.second, + unit=syu.pascal * syu.second, + doc="Dynamic viscosity of the electrolyte", +) +density = symbolic.Parameter( + "electrolyte.density", + domain, + default=1399.8 * syu.kg / syu.m**3, + unit=syu.kg / syu.m**3, + doc="Density of the electrolyte", +) +inlet_velocity = symbolic.Parameter( + "flow.inlet_velocity", + domain, + default=1e-7 * syu.m / syu.second, + unit=syu.m / syu.second, + doc="Velocity at the inlet", +) +# - + +# ## The flow +# +# The velocity is of quadratic and the pressure of linear elements, the +# Taylor-Hood pair. The inlet velocity scales the velocity, and the pressure +# drop of Darcy flow at it across the height of the electrode, evaluated at +# the defaults of the parameters with +# {py:func}`value_at_defaults `, scales +# the pressure: + +# + +height = 15e-3 +pressure_scale = symbolic.value_at_defaults( + viscosity * inlet_velocity * height / electrode.permeability() +) +u = symbolic.Field( + "u", + domain, + ("Lagrange", 2, (3,)), + unit=syu.m / syu.second, + scale=symbolic.value_at_defaults(inlet_velocity), +) +p = symbolic.Field("p", domain, ("Lagrange", 1), unit=syu.pascal, scale=pressure_scale) +# - + +# The flow of {eq}`brinkman-forchheimer` is a sum of +# {py:class}`BrinkmanForchheimer ` in the +# electrode and {py:class}`NavierStokes ` in the +# channels, with grad-div stabilization, +# {py:class}`GradDiv `, whose weight RFBniCS takes as +# a small part of the viscosity. Its creeping-flow part, +# {py:class}`Brinkman ` and +# {py:class}`Stokes `, is linear, and its solution the +# start of Newton's method for the whole flow. Both are checked for units +# and solved in dimensionless form: + +# + +grad_div = GradDiv("grad_div", u, 1e-3 * viscosity) +creeping = ( + Brinkman( + "electrode_creeping_flow", + u, + p, + electrode.porosity, + electrode.permeability(), + viscosity, + markers=electrode.markers, + ) + + Stokes("channel_creeping_flow", u, p, viscosity, markers=channels) + + grad_div +) +flow = ( + BrinkmanForchheimer( + "electrode_flow", + u, + p, + electrode.porosity, + electrode.permeability(), + electrode.forchheimer(), + viscosity, + density, + markers=electrode.markers, + ) + + NavierStokes("channel_flow", u, p, viscosity, density, markers=channels) + + grad_div +) +flow_systems = [ + symbolic.build_system(model, check_units=True, factorize_equations=True) + for model in (creeping, flow) +] +flow_bcs = [ + symbolic.DirichletBC(u, "Inlet", ufl.as_vector([inlet_velocity, 0.0, 0.0])), + symbolic.DirichletBC(u, ("Walls", "CurrentCollector", "Membrane"), 0.0), +] +# - + +# ## The half-cell +# +# The ions, the reservoir and the half-cell are those of +# {doc}`vanadium_half_cell_2d`. The electrolyte fills the `channels` too. +# The membrane is larger than the current collector, which the channels +# partly cover, so the cell's current density, given at the membrane, is +# scaled at the collector by the ratio of their areas: + +# + +diffusivity = syu.m**2 / syu.second +concentration = syu.mol / syu.m**3 +ions = ( + Ion("V2", 2, 8.10e-12 * diffusivity, positive=True), + Ion("V3", 3, 1.65e-11 * diffusivity, positive=True), + Ion("H", 1, 9.312e-9 * diffusivity), + Ion("HSO4", -1, 1.33e-9 * diffusivity), + Ion("SO4", -2, 1.065e-9 * diffusivity), +) +reservoir = { + "V2": 300.0 * concentration, + "V3": 300.0 * concentration, + "H": 1000.0 * concentration, + "HSO4": 1000.0 * concentration, +} +cell = RedoxHalfCell( + domain, + ions, + reservoir, + reduced="V2", + oxidized="V3", + electrons=1, + electrode=electrode, + oxidation_on_discharge=True, + rate_constant=2.5e-7 * syu.m / syu.second, + formal_potential=-0.36 * syu.volt, + channels=channels, + velocity=u, + outlet=("Outlet",), + proton_flux=True, + density=density, + viscosity=viscosity, +) +inlet = cell.inlet_conditions("Inlet") +# - + +# ## Solving the flow +# +# Everything is discretized together, and the creeping flow and then the +# whole flow are solved: + +# + +disc = discrete.discretize( + [*flow_systems, *flow_bcs, cell.model, *inlet], + discrete.DomainMaps({domain: marked}), + max_quadrature_degree=4, +) +for name, system in zip(("Creeping flow", "Flow"), flow_systems, strict=True): + solver = solvers.Solver(disc[system], flow_bcs) + solver.solve() + print(f"{name}: {solver.iterations} Newton iterations") +# - + +# The flow is so slow that its inertia, and the Forchheimer drag, change it +# by less than Newton's tolerance: the creeping flow is already the flow. It +# is incompressible, so what enters at the inlet leaves at the outlet. The +# velocity of Taylor-Hood elements is not divergence-free cell by cell, but +# the constants are among the test functions of the pressure, so it +# conserves mass over the whole cell, and the outflow matches the inflow to +# round-off: + +# + +n = ufl.FacetNormal(domain) +inflow = -disc.assemble_scalar(ufl.dot(u, n) * domain.ds("Inlet")) +outflow = disc.assemble_scalar(ufl.dot(u, n) * domain.ds("Outlet")) +print(f"Outflow over inflow: {outflow / inflow:.8f}") +assert inflow > 0 and np.isclose(outflow, inflow, rtol=1e-8) +# - + +# ## Charge, rest and discharge +# +# The cell starts at rest, with the current balanced between the membrane +# and the collector, and follows the schedule of +# {doc}`vanadium_half_cell_2d`. As there, the simulation's model is the +# half-cell alone, and every step takes the velocity solved above as known: + +# + +cell.balance_currents(disc) +cell.initialize(disc) +schedule = ( + ("charge", 0.5, ((20, 0.1), (10, 1.0), (10, 10.0))), + ("rest", 0.0, ((5, 0.1), (5, 1.0), (10, 10.0))), + ("discharge", 0.5, ((20, 0.1), (10, 1.0), (10, 10.0))), +) +sign = {"charge": -1.0, "rest": 0.0, "discharge": 1.0} +largest_step = max(dt for *_, steps in schedule for _, dt in steps) +clock = symbolic.Clock(domain, dt=schedule[0][2][0][1]) +simulation = solvers.Simulation( + disc, + cell.model, + inlet, + clock=clock, + scheme=solvers.BDF(1), + petsc_options={"snes_linesearch_type": "bt", "snes_max_it": 100}, + check_units=True, + factorize_equations=True, + time_scale=largest_step, +) +# - + +# The simulation runs in parts, as in {doc}`vanadium_half_cell_2d`, and keeps +# the local state of charge at the end of the charge: + +# + +local = cell.local_state_of_charge() +space = disc.function_space(cell.concentrations["V2"]) +charged = dolfinx.fem.Function(space, name="soc") +records: list[tuple[float, str, float, float]] = [] +for mode, current, steps in schedule: + disc.set({cell.current: sign[mode] * current}) + for count, dt in steps: + disc.set({clock.dt: dt}) + simulation.run( + num_steps=count, + callback=lambda t, step, mode=mode: records.append( + (t, mode, cell.voltage(disc), cell.state_of_charge(disc)) + ), + ) + if mode == "charge": + charged.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]) +charge = np.array([r[3] for r in records]) +# - + +# Charging raises the state of charge and discharging lowers it, and the +# voltage is lower on charge than on discharge: + +# + +steps = {mode: np.diff(charge[modes == mode]) for mode in ("charge", "discharge")} +print(f"State of charge after charge: {charge[modes == 'charge'][-1]:.6f}") +print(f"State of charge at the end: {charge[-1]:.6f}") +assert np.all(steps["charge"] > 0) and np.all(steps["discharge"] < 0) +assert voltage[modes == "charge"].mean() < voltage[modes == "discharge"].mean() +assert np.isfinite(voltage).all() +# - + +# The voltage and the state of charge through the schedule: + +fig, (top, bottom) = plt.subplots(2, 1, sharex=True) +for mode in ("charge", "rest", "discharge"): + top.plot(times[modes == mode], 1e3 * voltage[modes == mode], ".", label=mode) + bottom.plot(times[modes == mode], charge[modes == mode], ".", label=mode) +top.set_ylabel("Voltage (mV)") +bottom.set_ylabel("State of charge") +bottom.set_xlabel("Time (s)") +top.legend() +plt.show() + +# The local state of charge at the end of the charge, on the membrane side +# of the electrode, in millimetres: + +grid = gather_grid(charged) +if MPI.COMM_WORLD.rank == 0: + grid.points[:] *= 1e3 + plotter = pyvista.Plotter(window_size=[640, 640]) + plotter.add_mesh(grid, scalars="soc", scalar_bar_args={"title": "SOC"}) + plotter.view_yz(negative=True) + plotter.show() + +# ## References +# +# ```{footbibliography} +# ``` diff --git a/references.bib b/references.bib index 1b06fae..4d7a89e 100644 --- a/references.bib +++ b/references.bib @@ -257,3 +257,72 @@ @book{CarslawJaeger1959 address = {Oxford}, year = {1959}, } + +@article{Barzegari2026, + author = {Barzegari, Mojtaba and de Waal, Martin and Bordignon, Davide and Rosinaite, Adele and Liu, Baichen and Forner-Cuenca, Antoni}, + title = {{RfbFoam}: An open-source framework for porous electrode modeling in redox flow batteries}, + journal = {ChemRxiv}, + year = {2026}, + doi = {10.26434/chemrxiv.10001658/v1}, +} + +@book{NewmanBalsara2021, + author = {Newman, John and Balsara, Nitash P.}, + title = {Electrochemical Systems}, + edition = {4}, + publisher = {Wiley}, + address = {Hoboken}, + year = {2021}, +} + +@book{Fogler2020, + author = {Fogler, H. Scott}, + title = {Elements of Chemical Reaction Engineering}, + edition = {6}, + publisher = {Pearson}, + address = {Boston}, + year = {2020}, +} + +@article{Thiele1939, + author = {Thiele, E. W.}, + title = {Relation between catalytic activity and size of particle}, + journal = {Industrial \& Engineering Chemistry}, + volume = {31}, + number = {7}, + pages = {916--920}, + year = {1939}, + doi = {10.1021/ie50355a027}, +} + +@misc{Aghabarari2026, + author = {Aghabarari, Amirhossein and Dokken, J{\o}rgen S. and Horsch, Martin T. and Valseth, Eirik and Janssen, Mathijs}, + title = {{RFBniCS}: An open-source simulation framework for redox flow batteries}, + year = {2026}, + eprint = {2610.01840}, + archivePrefix = {arXiv}, + primaryClass = {physics.comp-ph}, + url = {https://arxiv.org/abs/2610.01840}, +} + +@article{MunozPerales2022, + author = {Mu{\~n}oz-Perales, Vanesa and Berling, Sabrina and Garc{\'i}a-Quismondo, Enrique and Garc{\'i}a-Salaberri, Pablo A. and Palma, Jes{\'u}s and Vera, Marcos and Ib{\'a}{\~n}ez, Santiago E.}, + title = {Investigating the Effects of Operation Variables on {All-Vanadium} Redox Flow Batteries Through an Advanced Unit-Cell Model}, + journal = {Journal of The Electrochemical Society}, + volume = {169}, + number = {10}, + pages = {100522}, + year = {2022}, + doi = {10.1149/1945-7111/ac975f}, +} + +@article{Krishnamurthy2011, + author = {Krishnamurthy, Deepak and Johansson, Erik O. and Lee, Jin Wook and Kjeang, Erik}, + title = {Computational modeling of microfluidic fuel cells with flow-through porous electrodes}, + journal = {Journal of Power Sources}, + volume = {196}, + number = {23}, + pages = {10019--10031}, + year = {2011}, + doi = {10.1016/j.jpowsour.2011.08.024}, +} diff --git a/src/sumfem/electrochemistry/__init__.py b/src/sumfem/electrochemistry/__init__.py index 7b181b2..51853bb 100644 --- a/src/sumfem/electrochemistry/__init__.py +++ b/src/sumfem/electrochemistry/__init__.py @@ -1,4 +1,4 @@ -"""Models of electrochemistry, built on :mod:`sumfem.models`. +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. @@ -9,6 +9,8 @@ - Kinetics: :class:`ExchangeCurrentButlerVolmer`, :class:`ButlerVolmer` and :class:`MassTransferButlerVolmer`, on :class:`ElectrodeReaction`, with the Butler-Volmer equation :func:`~sumfem.models.butler_volmer`. +- Half-cells of a redox flow battery: :class:`RedoxHalfCell`, negative or + positive, of any couple and electrolyte. The package imports from the core of :mod:`sumfem`, and never the other way around. @@ -16,6 +18,7 @@ from .charge import ElectrolyteConduction from .electrolyte import Electrolyte, Ion +from .half_cell import RedoxHalfCell from .kinetics import ( ButlerVolmer, ElectrodeReaction, @@ -32,6 +35,7 @@ "Region", "ElectrolyteConduction", "NernstPlanck", + "RedoxHalfCell", "ButlerVolmer", "ElectrodeReaction", "ExchangeCurrentButlerVolmer", diff --git a/src/sumfem/electrochemistry/electrolyte.py b/src/sumfem/electrochemistry/electrolyte.py index dc4f106..f06b9d0 100644 --- a/src/sumfem/electrochemistry/electrolyte.py +++ b/src/sumfem/electrochemistry/electrolyte.py @@ -6,6 +6,7 @@ import types import typing +import sympy import sympy.physics.units as syu import ufl @@ -33,8 +34,9 @@ class Ion: Attributes: name: Name of the ion, e.g. ``"V2"``. charge: Charge number ``z``. - diffusivity: Diffusion coefficient in free electrolyte, in m^2/s. - It is the default value of the parameter the electrolyte creates. + diffusivity: Diffusion coefficient in free electrolyte, in m^2/s or a + quantity with units, e.g. ``1e-9 * syu.m**2 / syu.s``. It is the + default value of the parameter the electrolyte creates. concentration: The field of a transported ion. ``None`` for the one ion whose concentration follows from electroneutrality. positive: Use a smooth positive part of the concentration in the @@ -44,7 +46,7 @@ class Ion: name: str charge: float - diffusivity: float + diffusivity: float | sympy.Expr concentration: Field | None = None positive: bool = False @@ -91,7 +93,8 @@ class Electrolyte: Raises: ValueError: If there are ions but none has a field, more than one - ion has none, or two ions share a name. + ion has none, the ion without one has no charge, or two ions + share a name. """ _name: str @@ -143,6 +146,11 @@ def __init__( "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 " diff --git a/src/sumfem/electrochemistry/half_cell.py b/src/sumfem/electrochemistry/half_cell.py new file mode 100644 index 0000000..5e42a80 --- /dev/null +++ b/src/sumfem/electrochemistry/half_cell.py @@ -0,0 +1,566 @@ +"""Half-cells of a redox flow battery. + +:class:`RedoxHalfCell` is a half-cell of any redox couple of dissolved ions, +negative or positive: a porous electrode, optionally beside free channels, +filled with an electrolyte that flows from a reservoir, with charge +conservation in the solid and the electrolyte, transport of every ion, and +Butler-Volmer kinetics of the couple at the fibres. +""" + +from __future__ import annotations + +import dataclasses +import types +import typing + +import sympy +import sympy.physics.units as syu +import ufl + +from sumfem.discrete import Discretization +from sumfem.models import Accumulation, BoundaryFlux, Conduction, Outflow +from sumfem.symbolic import ( + DirichletBC, + Domain, + Field, + Parameter, + Residual, + ResidualSum, +) + +from .charge import ElectrolyteConduction +from .electrolyte import CONCENTRATION, Electrolyte, Ion +from .kinetics import ButlerVolmer, MassTransferButlerVolmer +from .porous import PorousElectrode, Region +from .transport import NernstPlanck + +__all__ = ["RedoxHalfCell"] + + +class RedoxHalfCell: + """A half-cell of a redox flow battery, of a couple of dissolved ions. + + The macro-homogeneous porous-electrode model of RFBniCS, `Aghabarari et + al., 2026 `_, Section II, for any + couple, electrolyte and sign of the half-cell. + + The unknowns are the potentials ``phi_s`` of the solid and ``phi_l`` of + the electrolyte, in V, and the concentration ``c_`` of every ion + with a reservoir concentration, in mol/m^3. The one ion without follows + from electroneutrality. In the porous electrode, the couple + ``reduced <-> oxidized + n e-`` follows Butler-Volmer kinetics, limited + by mass transfer to the fibres when there is flow. In free channels the + electrolyte is not corrected for porosity, and there is no reaction or + solid phase. + + The model, :attr:`model`, is the sum of: + + - ``"solid"``: Ohmic conduction in the solid, + :class:`~sumfem.models.Conduction` with Bruggeman's effective + conductivity of the electrode; + - ``"reaction"``: the kinetics of the couple, + :class:`~sumfem.electrochemistry.ButlerVolmer`, or + :class:`~sumfem.electrochemistry.MassTransferButlerVolmer` with mass + transfer, which moves the current between the phases and consumes and + produces the ions of the couple; + - ``"collector_current"`` and ``"membrane_current"``: the current of the + cell into the solid and out of the electrolyte, + :class:`~sumfem.models.BoundaryFlux`; + - for the electrode's pores and, with a suffix ``"_channels"``, the free + channels: charge conservation in the electrolyte, + :class:`~sumfem.electrochemistry.ElectrolyteConduction`, and for every + transported ion its storage, :class:`~sumfem.models.Accumulation` + weighted with the porosity, and its transport, + :class:`~sumfem.electrochemistry.NernstPlanck`, named + ``"storage_"`` and ``"transport_"``; + - with an outlet, ``"outflow_"``: advection out of the cell, + :class:`~sumfem.models.Outflow`; + - with a proton flux, ``"proton_flux"``: the membrane current carried by + protons, a :class:`~sumfem.models.BoundaryFlux` of the ion ``"H"``. + + The cell is galvanostatic: the parameter ``".current"`` is the + current density through the membrane, in A/m^2, positive on discharge. + On discharge the couple is oxidized in the negative half-cell and + reduced in the positive one, so the current enters the solid at the + current collector of the first, and leaves it at that of the second. + The current collector carries the current density ``current * + area_ratio``, so that the total currents match when the two areas + differ. + + Args: + domain: The domain of the cell. + ions: The ions of the electrolyte, without fields: the cell creates + the field of every ion that has a reservoir concentration. + reservoir: The concentration of each transported ion in the + reservoir, in mol/m^3 or a quantity with units, the default of + the parameter ``".c0_"``. + reduced: Name of the reduced ion of the couple. With ``oxidized`` + and ``electrons`` it gives the reaction + ``oxidized + n e- <-> reduced``, e.g. V3 + e- <-> V2. + oxidized: Name of the oxidized ion of the couple, of ``electrons`` + more charge than the reduced one. + electrode: The porous electrode. + oxidation_on_discharge: Whether discharge oxidizes the couple, as in + the negative half-cell. + rate_constant: Standard rate constant of the couple in m/s or with + units, the default of the parameter ``".rate_constant"``. + formal_potential: Formal potential of the couple in V or with + units, the default of the parameter ``".formal_potential"``. + anodic_transfer: Anodic transfer coefficient, the default of the + parameter ``".anodic_transfer"``. + cathodic_transfer: Cathodic transfer coefficient, the default of the + parameter ``".cathodic_transfer"``. + electrons: Number of electrons the couple transfers. + channels: Cell markers of free channels. Empty for none. + velocity: Superficial velocity of the electrolyte, or ``None`` + without flow. + mass_transfer: Limit the kinetics by mass transfer to the fibres, + which needs a velocity. Defaults to whether there is one. + collector: Facet markers of the current collector. + membrane: Facet markers of the membrane. + outlet: Facet markers of an outlet, through which ions leave by + advection. Empty for none. + proton_flux: Carry the membrane current by protons, the ion ``"H"``, + through the membrane. + density: Density of the electrolyte, see + :class:`~sumfem.electrochemistry.Electrolyte`. Pass the + parameter of a flow model to share it. + viscosity: Dynamic viscosity of the electrolyte, given like the + density. + degree: Degree of the Lagrange elements of the unknowns. + name: Prefix of the cell's own parameters. + + Raises: + ValueError: For mass transfer or an outlet without a velocity, a + reservoir concentration of an ion that the electrolyte does not + have, a reservoir that lacks not exactly one ion, an ion of the + couple without one, a couple that does not conserve charge, or a + balanced ion without charge. + """ + + _domain: Domain + """The domain of the cell.""" + + _electrode: PorousElectrode + """The porous electrode.""" + + _regions: tuple[tuple[str, Region], ...] + """The regions filled with electrolyte: the electrode's pores, and the + free channels if any. Each comes with the suffix of its models' names.""" + + _electrolyte: Electrolyte + """The electrolyte.""" + + _reduced: str + """Name of the reduced ion of the couple.""" + + _oxidized: str + """Name of the oxidized ion of the couple.""" + + _phi_s: Field + """Potential of the solid.""" + + _phi_l: Field + """Potential of the electrolyte.""" + + _concentrations: types.MappingProxyType[str, Field] + """Concentration of each transported ion, by ion name.""" + + _inlet: types.MappingProxyType[str, Parameter] + """Reservoir concentration of each transported ion, by ion name.""" + + _current: Parameter + """Current density through the membrane, positive on discharge.""" + + _area_ratio: Parameter + """Ratio of the membrane area to the current collector area.""" + + _formal_potential: Parameter + """Formal potential of the couple.""" + + _reaction: ButlerVolmer + """The reaction in the electrode.""" + + _model: ResidualSum + """The sum of all models of the cell.""" + + _collector: tuple[str, ...] + """Facet markers of the current collector.""" + + _membrane: tuple[str, ...] + """Facet markers of the membrane.""" + + _functionals: types.MappingProxyType[str, ufl.Form] + """The integrals behind the cell's diagnostics, created once so that a + discretization compiles each of them only once.""" + + def __init__( + self, + domain: Domain, + ions: typing.Sequence[Ion], + reservoir: typing.Mapping[str, float | sympy.Expr], + reduced: str, + oxidized: str, + electrode: PorousElectrode, + oxidation_on_discharge: bool, + rate_constant: float | sympy.Expr, + formal_potential: float | sympy.Expr, + anodic_transfer: float = 0.5, + cathodic_transfer: float = 0.5, + electrons: int = 1, + channels: tuple[str, ...] = (), + velocity: ufl.core.expr.Expr | None = None, + mass_transfer: bool | None = None, + collector: tuple[str, ...] = ("CurrentCollector",), + membrane: tuple[str, ...] = ("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, + degree: int = 1, + name: str = "cell", + ): + # Mass transfer to the fibres depends on the speed of the flow, so it + # is on by default whenever there is one + mass_transfer = velocity is not None if mass_transfer is None else mass_transfer + if velocity is None and (mass_transfer or outlet): + raise ValueError("Mass transfer and an outlet need a velocity.") + unknown = sorted(set(reservoir) - {ion.name for ion in ions}) + if unknown: + raise ValueError(f"The electrolyte has no ions {unknown}.") + # The one ion without a reservoir concentration has no field, and + # follows from electroneutrality, which closes the equations + balanced = [ion.name for ion in ions if ion.name not in reservoir] + if len(balanced) != 1: + raise ValueError( + "The reservoir must give the concentrations of all ions but one, " + "whose concentration follows from electroneutrality; it lacks " + f"{balanced}." + ) + # The couple is the reaction oxidized + n e- <-> reduced, which must + # conserve charge + charges = {ion.name: ion.charge for ion in ions} + for ion in (reduced, oxidized): + if ion not in reservoir: + raise ValueError( + f"The ion {ion!r} of the couple needs a reservoir concentration." + ) + if charges[oxidized] - charges[reduced] != electrons: + raise ValueError( + f"{oxidized} + {electrons} e- <-> {reduced} does not conserve " + f"charge: the charges are {charges[oxidized]} and " + f"{charges[reduced]}." + ) + self._domain = domain + self._electrode = electrode + self._reduced, self._oxidized = reduced, oxidized + # The regions filled with electrolyte, with the suffix of their models: + # the pores of the electrode, whose porosity weights storage and + # corrects transport, and free channels, of porosity one + regions = [("", self._electrode.region)] + if channels: + regions.append(("_channels", Region(tuple(channels)))) + self._regions = tuple(regions) + + # The unknowns. A concentration is of the order of a mole per litre, + # its scale for the dimensional analysis and the scaled equations + element = ("Lagrange", degree) + self._phi_s = Field("phi_s", domain, element, unit=syu.volt) + self._phi_l = Field("phi_l", domain, element, unit=syu.volt) + self._concentrations = types.MappingProxyType( + { + ion.name: Field( + f"c_{ion.name}", domain, element, unit=CONCENTRATION, scale=1e3 + ) + for ion in ions + if ion.name in reservoir + } + ) + # The electrolyte, with the field of every transported ion. The ion + # without a reservoir concentration, and so without a field, follows + # from electroneutrality. Its density and viscosity enter the mass + # transfer, and may be shared with a model of the flow + self._electrolyte = Electrolyte( + "electrolyte", + domain, + [ + dataclasses.replace( + ion, concentration=self._concentrations.get(ion.name) + ) + for ion in ions + ], + density=density, + viscosity=viscosity, + ) + # The concentrations of the reservoir: the initial state, and the + # values of the conditions at an inlet, see inlet_conditions + self._inlet = types.MappingProxyType( + { + ion: Parameter( + f"{name}.c0_{ion}", + domain, + default=value, + unit=CONCENTRATION, + doc=f"Reservoir concentration of {ion}", + ) + for ion, value in reservoir.items() + } + ) + # The current density of the cell through the membrane, and the + # ratio of the areas that makes the current through the collector + # equal, see balance_currents + self._current = Parameter( + f"{name}.current", + domain, + default=0.0, + unit=syu.ampere / syu.m**2, + doc="Current density through the membrane, positive on discharge", + ) + self._area_ratio = Parameter( + f"{name}.area_ratio", + domain, + default=1.0, + doc="Membrane area over current collector area", + ) + # The kinetics of the couple, parameters of the cell that the reaction + # shares, and the specific area of the electrode + self._formal_potential = Parameter( + f"{name}.formal_potential", + domain, + default=formal_potential, + unit=syu.volt, + doc=f"Formal potential of {reduced}/{oxidized}", + ) + kinetics: dict[str, typing.Any] = dict( + specific_area=self._electrode.specific_area, + formal_potential=self._formal_potential, + rate_constant=Parameter( + f"{name}.rate_constant", + domain, + default=rate_constant, + unit=syu.m / syu.second, + doc=f"Standard rate constant of {reduced}/{oxidized}", + ), + anodic_transfer=Parameter( + f"{name}.anodic_transfer", + domain, + default=anodic_transfer, + doc=f"Anodic transfer coefficient of {reduced}/{oxidized}", + ), + cathodic_transfer=Parameter( + f"{name}.cathodic_transfer", + domain, + default=cathodic_transfer, + doc=f"Cathodic transfer coefficient of {reduced}/{oxidized}", + ), + electrons=electrons, + markers=self._electrode.markers, + ) + # The reaction in the electrode: Butler-Volmer at the concentrations + # of the pores, or at the surface of the fibres with mass transfer + e = self._electrolyte + if mass_transfer: + assert velocity is not None + self._reaction = MassTransferButlerVolmer( + "reaction", + self._phi_s, + self._phi_l, + e, + reduced, + oxidized, + velocity, + fiber_diameter=self._electrode.fiber_diameter, + **kinetics, + ) + else: + self._reaction = ButlerVolmer( + "reaction", self._phi_s, self._phi_l, e, reduced, oxidized, **kinetics + ) + + # The current into the solid at the collector, and out of the + # electrolyte at the membrane, on discharge. Oxidation on discharge + # releases electrons into the solid, which leave through the + # collector, so the current of the cell enters there; reduction on + # discharge takes them, so the current leaves the solid there + inflow = self._current if oxidation_on_discharge else -self._current + models: list[Residual] = [ + Conduction( + "solid", + self._phi_s, + self._electrode.effective_conductivity(), + markers=self._electrode.markers, + ), + self._reaction, + BoundaryFlux( + "collector_current", self._phi_s, inflow * self._area_ratio, collector + ), + BoundaryFlux("membrane_current", self._phi_l, -inflow, membrane), + ] + # Every region filled with electrolyte conserves charge, and stores + # and transports every ion, with the porosity of the region + for suffix, region in self._regions: + models.append( + ElectrolyteConduction(f"electrolyte{suffix}", self._phi_l, e, region) + ) + for ion, c in self._concentrations.items(): + models.append( + Accumulation( + f"storage_{ion}{suffix}", + c, + weight=region.porosity, + markers=region.markers, + ) + ) + models.append( + NernstPlanck( + f"transport_{ion}{suffix}", + e[ion], + e, + self._phi_l, + region, + velocity, + ) + ) + # The ions that the flow carries out of the cell + if outlet: + assert velocity is not None + for ion, c in self._concentrations.items(): + models.append(Outflow(f"outflow_{ion}", c, velocity, outlet)) + # The membrane passes the current as protons: their flux is the + # current over their charge and Faraday's constant + if proton_flux: + models.append( + BoundaryFlux( + "proton_flux", + self._concentrations["H"], + -inflow / (e["H"].charge * e.F), + membrane, + ) + ) + self._model = ResidualSum(models, name=name) + + # The integrals of the diagnostics: the areas and mean potentials of + # the collector and the membrane, for the voltage, and the amounts of + # the couple in the electrolyte of every region, for the state of + # charge + self._collector = tuple(collector) + self._membrane = tuple(membrane) + ds_collector = domain.ds(*self._collector) + ds_membrane = domain.ds(*self._membrane) + functionals = { + "collector_area": 1.0 * ds_collector, + "membrane_area": 1.0 * ds_membrane, + "collector_potential": self._phi_s * ds_collector, + "membrane_potential": self._phi_l * ds_membrane, + } + for role, ion in (("reduced", reduced), ("oxidized", oxidized)): + c = self._concentrations[ion] + functionals[f"amount_{role}"] = sum( + region.storage * c * domain.dx(*region.markers) + for _, region in self._regions + ) + self._functionals = types.MappingProxyType(functionals) + + @property + def model(self) -> Residual: + """The sum of all models of the cell.""" + return self._model + + @property + def electrode(self) -> PorousElectrode: + """The porous electrode.""" + return self._electrode + + @property + def electrolyte(self) -> Electrolyte: + """The electrolyte.""" + return self._electrolyte + + @property + def reaction(self) -> ButlerVolmer: + """The reaction in the electrode.""" + return self._reaction + + @property + def phi_s(self) -> Field: + """Potential of the solid.""" + return self._phi_s + + @property + def phi_l(self) -> Field: + """Potential of the electrolyte.""" + return self._phi_l + + @property + def concentrations(self) -> typing.Mapping[str, Field]: + """Concentration of each transported ion, by ion name.""" + return self._concentrations + + @property + def inlet(self) -> typing.Mapping[str, Parameter]: + """Reservoir concentration of each transported ion, by ion name.""" + return self._inlet + + @property + def current(self) -> Parameter: + """Current density through the membrane, positive on discharge.""" + return self._current + + @property + def area_ratio(self) -> Parameter: + """Ratio of the membrane area to the current collector area.""" + return self._area_ratio + + def inlet_conditions(self, markers: str | tuple[str, ...]) -> list[DirichletBC]: + """The reservoir concentration of every transported ion on some facets.""" + return [ + DirichletBC(c, markers, self._inlet[ion]) + for ion, c in self._concentrations.items() + ] + + 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. + """ + values: dict[str | ufl.Coefficient, typing.Any] = { + c: discretization.value(self._inlet[ion]) + for ion, c in self._concentrations.items() + } + values[self._phi_s] = 0.0 + values[self._phi_l] = -discretization.value(self._formal_potential) + discretization.set(values) + + def balance_currents(self, discretization: Discretization) -> None: + """Set the area ratio from the areas of membrane and current collector.""" + assemble = discretization.assemble_scalar + ratio = assemble(self._functionals["membrane_area"]) / assemble( + self._functionals["collector_area"] + ) + discretization.set({self._area_ratio: ratio}) + + def voltage(self, discretization: Discretization) -> float: + """The half-cell voltage in V. + + That is the mean solid potential on the current collector minus the + mean electrolyte potential on the membrane. + """ + assemble = discretization.assemble_scalar + f = self._functionals + return assemble(f["collector_potential"]) / assemble( + f["collector_area"] + ) - assemble(f["membrane_potential"]) / assemble(f["membrane_area"]) + + def state_of_charge(self, discretization: Discretization) -> float: + """The share of the reduced ion in the couple, over the cell.""" + assemble = discretization.assemble_scalar + reduced = assemble(self._functionals["amount_reduced"]) + oxidized = assemble(self._functionals["amount_oxidized"]) + return reduced / max(reduced + oxidized, 1e-30) + + def local_state_of_charge(self) -> ufl.core.expr.Expr: + """The share of the reduced ion in the couple, ``c_red / (c_red + c_ox)``.""" + c_red = self._concentrations[self._reduced] + c_ox = self._concentrations[self._oxidized] + return c_red / (c_red + c_ox + self._electrolyte.floor) diff --git a/src/sumfem/electrochemistry/kinetics.py b/src/sumfem/electrochemistry/kinetics.py index cf1400b..ed09660 100644 --- a/src/sumfem/electrochemistry/kinetics.py +++ b/src/sumfem/electrochemistry/kinetics.py @@ -442,7 +442,11 @@ class MassTransferButlerVolmer(ButlerVolmer): ``k_m (c - c_s)`` from the bulk. The mass-transfer coefficient of each ion is ``k_m = Sh D / d`` with the Sherwood number ``Sh = a Re^b`` and the Reynolds number ``Re = rho |u| d / mu`` of the flow past fibers of - diameter ``d``. ``D`` is the diffusivity in free electrolyte. + diameter ``d``. ``D`` is the diffusivity in free electrolyte. The + defaults ``a = 7`` and ``b = 0.4`` are the correlation of `Krishnamurthy + et al., 2011 `_, and the + surface concentrations those of equation 20 of `Aghabarari et al., 2026 + `_. Args: name: Name of the model. diff --git a/src/sumfem/electrochemistry/porous.py b/src/sumfem/electrochemistry/porous.py index 2fe4b0d..3386866 100644 --- a/src/sumfem/electrochemistry/porous.py +++ b/src/sumfem/electrochemistry/porous.py @@ -4,6 +4,7 @@ import dataclasses +import sympy import sympy.physics.units as syu import ufl @@ -53,7 +54,8 @@ class PorousElectrode: 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. + 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. @@ -90,11 +92,11 @@ def __init__( domain: Domain, markers: tuple[str, ...] = (), porosity: float | None = None, - conductivity: float | None = None, - specific_area: float | None = None, - fiber_diameter: 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 | None = None, + pore_length: float | sympy.Expr | None = None, ): self._markers = tuple(markers) self._porosity = Parameter( @@ -165,11 +167,13 @@ def effective_conductivity(self) -> ufl.core.expr.Expr: return bruggeman(self._conductivity, 1.0 - self._porosity) def permeability(self) -> ufl.core.expr.Expr: - """Kozeny-Carman permeability of the fibers.""" + """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.""" + """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/models/porous.py b/src/sumfem/models/porous.py index 7435b8f..45f0e75 100644 --- a/src/sumfem/models/porous.py +++ b/src/sumfem/models/porous.py @@ -19,6 +19,15 @@ def bruggeman( ) -> ufl.core.expr.Expr: """Bruggeman's effective property ``value * fraction**exponent``. + It corrects a transport property of a phase for the volume the other + phase takes and for the tortuosity of the paths through it. The exponent + ``3/2`` is that of spherical obstructions, see `Tjaden et al., 2016 + `_, where the tortuosity + scales with the fraction to the power ``-1/2``. It suits macroporous + electrodes; where the tortuosity is measured, the property + ``value * fraction / tortuosity`` is the alternative of `Barzegari et + al., 2026 `_. + Args: value: The property of the phase, e.g. a diffusivity or the conductivity of the solid. @@ -36,6 +45,14 @@ def kozeny_carman( ) -> ufl.core.expr.Expr: """Kozeny-Carman permeability ``d^2 eps^3 / (K (1 - eps)^2)`` of a fibrous medium. + An estimate of the permeability of a fibrous porous electrode, such as a + carbon felt, from the diameter of its fibres, with an empirical constant + ``K``. The carbon fibres of commercial cloths, papers and felts have + diameters of 7 to 9 micrometres, see `Barzegari et al., 2026 + `_. The estimate is of + limited use for complex media, whose permeability is best fitted to + measured pressure drops. + Args: fiber_diameter: Fiber diameter ``d``. porosity: Porosity ``eps``. @@ -70,6 +87,17 @@ def ergun_coefficient( ) -> ufl.core.expr.Expr: """Ergun's inertial drag coefficient ``1.75 (1 - eps) / (L eps^3)``. + The coefficient ``beta``, in 1/m, of Forchheimer's drag ``rho beta |u| u`` + on the superficial velocity ``u``, from Ergun's relation for randomly + packed beds of spheres of diameter ``L`` (`Ergun and Orning, 1949 + `_). For a fibrous electrode, such + as a carbon felt, it is an estimate only: the model of `Tomadakis and + Robertson, 2005 `_ is made for + random fibre structures, and fits to experiments give the coefficient of + a given electrode, e.g. of the order of Ergun's for the 3D printed + electrodes of `van der Heijden et al., 2023 + `_. + Args: porosity: Porosity ``eps``. length: Characteristic pore length ``L``. diff --git a/src/sumfem/symbolic/__init__.py b/src/sumfem/symbolic/__init__.py index 41af21f..2883bb5 100644 --- a/src/sumfem/symbolic/__init__.py +++ b/src/sumfem/symbolic/__init__.py @@ -48,7 +48,7 @@ from .domains import Domain, Subdomain from .expressions import as_expr from .system import System, build_system -from .terminals import Field, KnownCoefficient, Parameter +from .terminals import Field, KnownCoefficient, Parameter, value_at_defaults from .time import Clock, Dt from .transfer import Assignment, Interpolation, Transfer, Transferable from .ufl_forms import UFLResidual @@ -105,4 +105,5 @@ "equation_dimensions", "equation_scales", "pi_groups", + "value_at_defaults", ] diff --git a/src/sumfem/symbolic/terminals.py b/src/sumfem/symbolic/terminals.py index 6d04ac7..534bddf 100644 --- a/src/sumfem/symbolic/terminals.py +++ b/src/sumfem/symbolic/terminals.py @@ -39,7 +39,8 @@ from ..error_handling import ComponentError, check_name, qualified_name from .domains import Domain, domain_of from .elements import ElementSpec, as_element -from .si import Unit, to_si +from .expressions import evaluate_constant +from .si import Unit, has_units, si_factor, to_si _PIOLA = ( ufl.pullback.contravariant_piola, @@ -58,6 +59,7 @@ "Parameter", "TimeDerivative", "field_of", + "value_at_defaults", ] # Map from the part number of a field's test and trial functions to the @@ -442,14 +444,20 @@ class Parameter(ufl.Coefficient, QuantityMixin): named ``"."``. domain: The domain the parameter is defined on. shape: Value shape, ``()`` for a scalar. - default: Default value in ``unit``, used when no other value is set. - A scalar sets every component of a vector or tensor parameter. + default: Default value in ``unit``, used when no other value is set, + or a sympy quantity carrying its own units, e.g. ``0.5 * + syu.millimeter``, converted to ``unit``. A scalar sets every + component of a vector or tensor parameter. unit: Unit of the value, a sympy unit expression. ``None`` means dimensionless. The discretization stores values in SI base units. doc: Description, shown in documentation and archived metadata. The name must be passed positionally: ``ufl.Coefficient.__new__`` reads the first positional argument. + + Raises: + ValueError: If a default with units does not have the dimension of + ``unit``. """ _name: str @@ -467,7 +475,7 @@ def __init__( /, domain: Domain, shape: tuple[int, ...] = (), - default: npt.ArrayLike | None = None, + default: npt.ArrayLike | sympy.Expr | None = None, unit: Unit | None = None, doc: str = "", ): @@ -484,6 +492,12 @@ def __init__( space = ufl.FunctionSpace(domain, element, label=f"{name}#parameter") ufl.Coefficient.__init__(self, space) self._name = name + if has_units(default): + # A quantity with its own units, as a value in the declared unit + default = to_si(default, unit) / si_factor(unit) + elif isinstance(default, sympy.Expr) and default.is_number: + # A zero quantity, e.g. 0 * syu.pascal, which sympy reduces to 0 + default = float(default) self._default = ( None if default is None else _validate_value(name, default, shape) ) @@ -522,6 +536,51 @@ def doc(self) -> str: return self._doc +def value_at_defaults(expr: ufl.core.expr.Expr) -> float | npt.NDArray[np.float64]: + """The value of an expression of parameters at their defaults, in SI units. + + Use it for a value that follows from the parameters of a model before + anything is discretized, e.g. the permeability of an electrode or the + pressure drop of a flow, as the scale of a field. + + Args: + expr: An expression of parameters and numbers. + + Raises: + ComponentError: If the expression depends on anything but parameters + with defaults, e.g. a field or the spatial coordinate. + + Example: + The scale of a pressure from a viscosity, a velocity, a length and a + permeability, all parameters:: + + p = Field("p", domain, ("Lagrange", 1), unit=syu.pascal, + scale=value_at_defaults(mu * v * L / K)) + """ + defaults: dict[ufl.Coefficient, ufl.core.expr.Expr] = {} + for c in ufl.algorithms.extract_coefficients(expr): + default = c.default_si if isinstance(c, Parameter) else None + if default is None: + name = getattr(c, "name", str(c)) + raise ComponentError( + f"The expression depends on {name!r}, which is not a parameter " + "with a default." + ) + defaults[c] = ( + ufl.as_ufl(float(default)) + if default.ndim == 0 + else (ufl.as_tensor(default.tolist())) + ) + replaced = ufl.replace(expr, defaults) if defaults else expr + try: + return evaluate_constant(replaced) + except TypeError as error: + raise ComponentError( + "The expression depends on more than parameters, e.g. on the " + "spatial coordinate." + ) from error + + class KnownCoefficient(ufl.Coefficient, QuantityMixin): """Named, spatially varying data that a model reads but does not solve for. diff --git a/tests/electrochemistry/half_cell/conftest.py b/tests/electrochemistry/half_cell/conftest.py new file mode 100644 index 0000000..76f9e95 --- /dev/null +++ b/tests/electrochemistry/half_cell/conftest.py @@ -0,0 +1,435 @@ +"""Meshes and the original RFBniCS forms, for comparison with the models. + +The reference residuals are transcribed from the examples of RFBniCS +(https://github.com/Ah-Aghabarari/RFBniCS, MIT licence; Aghabarari et al., +2026, https://arxiv.org/abs/2610.01840), with one change: +each unknown is its own function, where RFBniCS splits one blocked +function. The tests pass the functions of a discretization, so both sides +assemble on the same spaces. + +The tests use pytest's importlib import mode, so these are handed to tests +through the ``reference`` fixture instead of being imported. +""" + +import types + +from mpi4py import MPI + +import dolfinx.fem +import dolfinx.mesh +import numpy as np +import pytest +import sympy.physics.units as syu +import ufl + +from sumfem.discrete import MeshData +from sumfem.electrochemistry import Ion, PorousElectrode, RedoxHalfCell +from sumfem.symbolic import Domain + +W = 2.5e-3 +"""Width of the electrode in the examples, in m.""" + +H = 1e-2 +"""Height of the electrode in the 2D example, in m.""" + +PARAMETERS = { + "sigma": 3402.1, + "a": 2.3e4, + "alpha_a": 0.5, + "alpha_c": 0.5, + "k0": 2.5e-7, + "E_formal": -0.36, + "df": 0.5e-6, + "F": 96485.33212, + "R": 8.314462618, + "T": 298.15, + "D_V2": 8.10e-12, + "D_V3": 1.65e-11, + "D_H": 9.312e-9, + "D_HSO4": 1.33e-9, + "D_SO4": 1.065e-9, + "c0_V2": 300.0, + "c0_V3": 300.0, + "c0_H": 1000.0, + "c0_HSO4": 1000.0, + "z_V2": 2.0, + "z_V3": 3.0, + "z_H": 1.0, + "z_HSO4": -1.0, + "z_SO4": -2.0, + "mu": 6.963e-3, + "rho": 1399.8, + "Kc": 5.55, + "epsilon": 0.94, + "u_in": 1.0e-7, + "p_out": 3.0e5, + "Lp": 100.0e-6, +} +"""Default parameters of the RFBniCS examples.""" + +QUADRATURE_DEGREE = 4 +"""The quadrature degree of the examples for degree-one elements.""" + + +def interval(num_cells: int = 20) -> MeshData: + """The 1D example's electrode, with its current collector and membrane.""" + mesh = dolfinx.mesh.create_interval(MPI.COMM_WORLD, num_cells, [0.0, W]) + return MeshData.from_locators( + mesh, + { + 0: { + "CurrentCollector": lambda x: np.isclose(x[0], 0.0), + "Membrane": lambda x: np.isclose(x[0], W), + } + }, + boundary=True, + ) + + +def rectangle(nx: int = 5, ny: int = 8) -> MeshData: + """The 2D example's electrode, with flow from the inlet at the bottom.""" + mesh = dolfinx.mesh.create_rectangle( + MPI.COMM_WORLD, + np.array([[0.0, 0.0], [W, H]]), + [nx, ny], + cell_type=dolfinx.mesh.CellType.quadrilateral, + ) + return MeshData.from_locators( + mesh, + { + 1: { + "CurrentCollector": lambda x: np.isclose(x[0], 0.0), + "Membrane": lambda x: np.isclose(x[0], W), + "Inlet": lambda x: np.isclose(x[1], 0.0), + "Outlet": lambda x: np.isclose(x[1], H), + } + }, + boundary=True, + ) + + +def box(n: int = 2) -> MeshData: + """A unit cube with an electrode and two channels, as in the 3D example. + + The electrode fills ``x > 1/2``. The inlet and outlet channels fill the + rest, below and above ``y = 1/2``, and open at ``x = 0``. + """ + mesh = dolfinx.mesh.create_unit_cube(MPI.COMM_WORLD, 2 * n, 2 * n, 2 * n) + half = 0.5 + 1e-10 + cells = MeshData.from_locators( + mesh, + { + 3: { + "InletChannel": lambda x: (x[0] <= half) & (x[1] <= half), + "OutletChannel": lambda x: (x[0] <= half) & (x[1] >= 0.5 - 1e-10), + "Electrode": lambda x: x[0] >= 0.5 - 1e-10, + } + }, + ) + sides = { + "Inlet": lambda x: np.isclose(x[0], 0.0) & (x[1] <= half), + "Outlet": lambda x: np.isclose(x[0], 0.0) & (x[1] >= 0.5 - 1e-10), + "CurrentCollector": lambda x: np.isclose(x[2], 0.0) & (x[0] >= 0.5 - 1e-10), + "Membrane": lambda x: np.isclose(x[0], 1.0), + } + + def walls(x): + """The rest of the boundary: markers of facets do not overlap.""" + return ~np.logical_or.reduce([side(x) for side in sides.values()]) + + facets = MeshData.from_locators(mesh, {2: {"Walls": walls, **sides}}, boundary=True) + return MeshData( + mesh, {**cells.tags, **facets.tags}, {**cells.markers, **facets.markers} + ) + + +def measures(marked: MeshData, degree: int | None = QUADRATURE_DEGREE): + """Cell and facet measures of a marked mesh, and a lookup of marker ids.""" + mesh = marked.mesh + tdim = mesh.topology.dim + metadata = None if degree is None else {"quadrature_degree": degree} + dx = ufl.Measure( + "dx", domain=mesh, subdomain_data=marked.tags.get(tdim), metadata=metadata + ) + ds = ufl.Measure( + "ds", domain=mesh, subdomain_data=marked.tags.get(tdim - 1), metadata=metadata + ) + + def ids(dim, *names): + return tuple(marked.markers[dim][name][0] for name in names) + + return dx, ds, ids + + +def _constants(mesh, values): + return { + name: dolfinx.fem.Constant(mesh, np.float64(value)) + for name, value in values.items() + } + + +def electrochemistry( + marked: MeshData, + u, + u_old, + v, + dt, + current, + *, + velocity=None, + regions: bool = False, + area_ratio: float = 1.0, + membrane_flux: bool = False, + outlet: bool = False, +): + """The residual of the RFBniCS examples, in one form. + + Without a velocity it is the 1D example's, with a velocity the 2D + example's, and with regions the 3D example's, where the electrode and + channel cells differ. + + Args: + marked: The mesh. + u: Functions of phi_s, phi_l, c_V2, c_V3, c_H and c_HSO4. + u_old: Functions of the previous concentrations of V2, V3, H, HSO4. + v: Test functions of the six unknowns. + dt: Time step, a number or a constant. + current: Membrane current density, signed, a number or a constant. + velocity: The flow velocity, if any. + regions: Whether the mesh has electrode and channel cells. + area_ratio: Membrane area over current collector area. + membrane_flux: Whether protons carry the membrane current. + outlet: Whether ions leave through an outlet. + """ + mesh = marked.mesh + c = _constants(mesh, PARAMETERS) + tdim = mesh.topology.dim + dx, ds, ids = measures(marked) + if regions: + dxE = dx(ids(tdim, "Electrode")) + dxAll = dxE + dx(ids(tdim, "InletChannel", "OutletChannel")) + Q0 = dolfinx.fem.functionspace(mesh, ("DG", 0)) + mask_E = dolfinx.fem.Function(Q0) + mask_E.x.array[marked.entities(tdim, "Electrode")] = 1.0 + mask_F = dolfinx.fem.Function(Q0) + mask_F.x.array[:] = 1.0 - mask_E.x.array + else: + dxE = dxAll = dx + mask_E = mask_F = None + ds_cc = ds(ids(tdim - 1, "CurrentCollector")) + ds_mem = ds(ids(tdim - 1, "Membrane")) + + phi_s, phi_l, c_V2, c_V3, c_H, c_HSO4 = u + c_V2_old, c_V3_old, c_H_old, c_HSO4_old = u_old + v_s, v_l, v_c_V2, v_c_V3, v_c_H, v_c_HSO4 = v + _F, _R, _T, _eps = c["F"], c["R"], c["T"], c["epsilon"] + dt_c = ( + dt + if isinstance(dt, dolfinx.fem.Constant) + else dolfinx.fem.Constant(mesh, np.float64(dt)) + ) + _i_mem = ( + current + if isinstance(current, dolfinx.fem.Constant) + else dolfinx.fem.Constant(mesh, np.float64(current)) + ) + _i_cc = _i_mem * area_ratio + + # Electroneutrality + c_SO4 = ( + -(c["z_V2"] * c_V2 + c["z_V3"] * c_V3 + c["z_H"] * c_H + c["z_HSO4"] * c_HSO4) + / c["z_SO4"] + ) + eps_c = dolfinx.fem.Constant(mesh, 1e-14) + c_V2_eff = 0.5 * (c_V2 + ufl.sqrt(c_V2 * c_V2 + eps_c * eps_c)) + c_V3_eff = 0.5 * (c_V3 + ufl.sqrt(c_V3 * c_V3 + eps_c * eps_c)) + + delta_phi = phi_s - phi_l - c["E_formal"] + A = c["k0"] * ufl.exp(c["alpha_a"] * _F * delta_phi / (_R * _T)) + B = c["k0"] * ufl.exp(-c["alpha_c"] * _F * delta_phi / (_R * _T)) + if velocity is None: + c_V2_s, c_V3_s = c_V2_eff, c_V3_eff + else: + eps_speed = dolfinx.fem.Constant(mesh, 1e-16) + u_mag = ufl.sqrt(ufl.dot(velocity, velocity) + eps_speed**2) + Re = c["rho"] * u_mag * c["df"] / c["mu"] + k_mR = 7.0 * c["D_V2"] / c["df"] * (Re**0.4) + k_mO = 7.0 * c["D_V3"] / c["df"] * (Re**0.4) + den = 1.0 + A / (k_mR + eps_c) + B / (k_mO + eps_c) + c_V3_s = ( + (1.0 + A / (k_mR + eps_c)) * c_V3_eff + (A / (k_mO + eps_c)) * c_V2_eff + ) / den + c_V2_s = ( + (1.0 + B / (k_mO + eps_c)) * c_V2_eff + (B / (k_mR + eps_c)) * c_V3_eff + ) / den + i_F = _F * (A * c_V2_s - B * c_V3_s) + if mask_E is not None: + i_F = mask_E * i_F + f = c["a"] * i_F + + def effective(D): + if mask_E is None: + return D * _eps**1.5 + return mask_E * (D * _eps**1.5) + mask_F * D + + D_V2, D_V3, D_H, D_HSO4, D_SO4 = ( + effective(c[f"D_{ion}"]) for ion in ("V2", "V3", "H", "HSO4", "SO4") + ) + storage = _eps if mask_E is None else mask_E * _eps + mask_F * 1.0 + + kappa = (_F**2 / (_R * _T)) * ( + c["z_V2"] ** 2 * D_V2 * c_V2_eff + + c["z_V3"] ** 2 * D_V3 * c_V3_eff + + c["z_H"] ** 2 * D_H * c_H + + c["z_HSO4"] ** 2 * D_HSO4 * c_HSO4 + + c["z_SO4"] ** 2 * D_SO4 * c_SO4 + ) + i_l = -kappa * ufl.grad(phi_l) - _F * ( + c["z_V2"] * D_V2 * ufl.grad(c_V2) + + c["z_V3"] * D_V3 * ufl.grad(c_V3) + + c["z_H"] * D_H * ufl.grad(c_H) + + c["z_HSO4"] * D_HSO4 * ufl.grad(c_HSO4) + + c["z_SO4"] * D_SO4 * ufl.grad(c_SO4) + ) + sigma = (1.0 - _eps) ** 1.5 * c["sigma"] + residual = ufl.inner(sigma * ufl.grad(phi_s), ufl.grad(v_s)) * dxE + f * v_s * dxE + residual += ufl.inner(-i_l, ufl.grad(v_l)) * dxAll - ufl.inner(f, v_l) * dxE + + def flux(D, z, conc): + N = -D * ufl.grad(conc) - D * z * _F * conc / (_R * _T) * ufl.grad(phi_l) + return N if velocity is None else N + velocity * conc + + N_V2 = flux(D_V2, c["z_V2"], c_V2) + N_V3 = flux(D_V3, c["z_V3"], c_V3) + N_H = flux(D_H, c["z_H"], c_H) + N_HSO4 = flux(D_HSO4, c["z_HSO4"], c_HSO4) + residual += ( + storage * (c_V2 - c_V2_old) / dt_c * v_c_V2 * dxAll + - ufl.inner(N_V2, ufl.grad(v_c_V2)) * dxAll + + (1.0 / _F) * f * v_c_V2 * dxE + ) + residual += ( + storage * (c_V3 - c_V3_old) / dt_c * v_c_V3 * dxAll + - ufl.inner(N_V3, ufl.grad(v_c_V3)) * dxAll + - (1.0 / _F) * f * v_c_V3 * dxE + ) + residual += ( + storage * (c_H - c_H_old) / dt_c * v_c_H * dxAll + - ufl.inner(N_H, ufl.grad(v_c_H)) * dxAll + ) + residual += ( + storage * (c_HSO4 - c_HSO4_old) / dt_c * v_c_HSO4 * dxAll + - ufl.inner(N_HSO4, ufl.grad(v_c_HSO4)) * dxAll + ) + residual -= _i_cc * v_s * ds_cc + residual += _i_mem * v_l * ds_mem + if membrane_flux: + residual += _i_mem / (c["z_H"] * _F) * v_c_H * ds_mem + if outlet: + n = ufl.FacetNormal(mesh) + ds_out = ds(ids(tdim - 1, "Outlet")) + for conc, test in ( + (c_V2, v_c_V2), + (c_V3, v_c_V3), + (c_H, v_c_H), + (c_HSO4, v_c_HSO4), + ): + residual += ufl.dot(velocity, n) * conc * test * ds_out + return residual + + +def _darcy_coefficient(mesh): + """The Darcy coefficient, permeability over viscosity, of the 2D example.""" + c = _constants(mesh, PARAMETERS) + return ( + (c["df"] ** 2 / (c["Kc"] * c["mu"])) + * (c["epsilon"] ** 3) + / ((1.0 - c["epsilon"]) ** 2) + ) + + +def darcy(marked: MeshData, p, q): + """The pressure equation of the 2D example, with the inflow.""" + mesh = marked.mesh + dx, ds, ids = measures(marked) + u_in = dolfinx.fem.Constant(mesh, np.float64(PARAMETERS["u_in"])) + a = ufl.inner(_darcy_coefficient(mesh) * ufl.grad(p), ufl.grad(q)) * dx + return a - u_in * q * ds(ids(mesh.topology.dim - 1, "Inlet")) + + +def darcy_velocity(marked: MeshData, p): + """The Darcy velocity of the 2D example.""" + return -_darcy_coefficient(marked.mesh) * ufl.grad(p) + + +DIFFUSIVITY = syu.m**2 / syu.second +CONCENTRATION = syu.mol / syu.m**3 + +IONS = ( + Ion("V2", 2.0, 8.10e-12 * DIFFUSIVITY, positive=True), + Ion("V3", 3.0, 1.65e-11 * DIFFUSIVITY, positive=True), + Ion("H", 1.0, 9.312e-9 * DIFFUSIVITY), + Ion("HSO4", -1.0, 1.33e-9 * DIFFUSIVITY), + Ion("SO4", -2.0, 1.065e-9 * DIFFUSIVITY), +) +"""The ions of the examples, sulfate by electroneutrality.""" + +RESERVOIR = { + "V2": 300.0 * CONCENTRATION, + "V3": 300.0 * CONCENTRATION, + "H": 1000.0 * CONCENTRATION, + "HSO4": 1000.0 * CONCENTRATION, +} +"""The concentrations of the reservoir of the examples.""" + + +def electrode(domain: Domain, markers: tuple[str, ...] = ()) -> PorousElectrode: + """The porous electrode of the examples.""" + return PorousElectrode( + "electrode", + domain, + markers, + porosity=0.94, + conductivity=3402.1 * syu.siemens / syu.m, + specific_area=2.3e4 / syu.m, + fiber_diameter=0.5 * syu.micrometer, + kozeny_carman_constant=5.55, + pore_length=100.0 * syu.micrometer, + ) + + +def half_cell(domain: Domain, porous=None, **kwargs) -> RedoxHalfCell: + """The negative half-cell of the examples, V2+ oxidized on discharge.""" + return RedoxHalfCell( + domain, + IONS, + RESERVOIR, + "V2", + "V3", + electrode(domain) if porous is None else porous, + oxidation_on_discharge=True, + rate_constant=2.5e-7 * syu.m / syu.second, + formal_potential=-0.36 * syu.volt, + **kwargs, + ) + + +@pytest.fixture +def reference(): + return types.SimpleNamespace( + PARAMETERS=PARAMETERS, + QUADRATURE_DEGREE=QUADRATURE_DEGREE, + W=W, + H=H, + interval=interval, + rectangle=rectangle, + box=box, + measures=measures, + IONS=IONS, + RESERVOIR=RESERVOIR, + electrode=electrode, + half_cell=half_cell, + electrochemistry=electrochemistry, + darcy=darcy, + darcy_velocity=darcy_velocity, + ) diff --git a/tests/electrochemistry/half_cell/test_models.py b/tests/electrochemistry/half_cell/test_models.py new file mode 100644 index 0000000..045f0f9 --- /dev/null +++ b/tests/electrochemistry/half_cell/test_models.py @@ -0,0 +1,259 @@ +"""The RFBniCS models: units, structure, schedules and post-processing.""" + +from mpi4py import MPI + +import dolfinx.fem +import dolfinx.la +import numpy as np +import pytest +import sympy.physics.units as syu + +from sumfem.discrete import DomainMaps, discretize +from sumfem.electrochemistry import RedoxHalfCell +from sumfem.models import ( + BoundaryFlux, + BrinkmanForchheimer, + Darcy, + GradDiv, + NavierStokes, +) +from sumfem.post_processing import evaluate_over_line, function_range +from sumfem.solvers import BDF +from sumfem.symbolic import Domain, Field, Parameter, build_system + + +def fluid(domain): + return ( + Parameter( + "electrolyte.viscosity", + domain, + default=6.963e-3, + unit=syu.pascal * syu.second, + ), + Parameter( + "electrolyte.density", domain, default=1399.8, unit=syu.kg / syu.m**3 + ), + ) + + +def assert_consistent(model): + """Units are consistent with and without time stepping, and when factorized.""" + build_system(model, check_units=True) + build_system(model, substitutions=BDF(1).substitutions(model), check_units=True) + scaled = build_system(model, factorize_equations=True) + assert all(scale > 0 for scale in scaled.equation_scales.values()) + + +def test_half_cells_have_consistent_units(reference): + line = Domain("cell", "interval") + cell = reference.half_cell(line) + assert [f.name for f in cell.model.equations] == [ + "phi_s", + "phi_l", + "c_V2", + "c_V3", + "c_H", + "c_HSO4", + ] + assert_consistent(cell.model) + + rect = Domain("cell", "quadrilateral") + viscosity, density = fluid(rect) + electrode = reference.electrode(rect) + p = Field("p", rect, ("Lagrange", 1), unit=syu.pascal) + inflow = Parameter("u_in", rect, default=1e-7, unit=syu.m / syu.second) + darcy = Darcy("darcy", p, electrode.permeability(), viscosity) + assert_consistent(darcy + BoundaryFlux("inflow", p, inflow, ("Inlet",))) + cell = reference.half_cell( + rect, + electrode, + velocity=darcy.superficial_velocity, + outlet=("Outlet",), + proton_flux=True, + density=density, + viscosity=viscosity, + ) + assert_consistent(cell.model) + # The flow is data for the electrochemistry + assert list(build_system(cell.model).known_coefficients) == ["p"] + + box = Domain("cell", "tetrahedron") + viscosity, density = fluid(box) + electrode = reference.electrode(box, ("Electrode",)) + channels = ("InletChannel", "OutletChannel") + u = Field("u", box, ("Lagrange", 2, (3,)), unit=syu.m / syu.second, scale=1e-7) + p = Field("p", box, ("Lagrange", 1), unit=syu.pascal) + flow = BrinkmanForchheimer( + "porous_flow", + u, + p, + electrode.porosity, + electrode.permeability(), + electrode.forchheimer(), + viscosity, + density, + markers=("Electrode",), + ) + NavierStokes("channel_flow", u, p, viscosity, density, markers=channels) + flow += GradDiv("grad_div", u, 1e-3 * viscosity) + assert_consistent(flow) + cell = reference.half_cell( + box, + electrode, + channels=channels, + velocity=u, + outlet=("Outlet",), + proton_flux=True, + density=density, + viscosity=viscosity, + ) + assert_consistent(cell.model) + names = [m.name for m in cell.model.components] + assert "transport_V2_channels" in names and "storage_H_channels" in names + assert "electrolyte_channels" in names + + +def test_shared_parameters_are_used_by_flow_and_kinetics(reference): + rect = Domain("cell", "quadrilateral") + viscosity, density = fluid(rect) + p = Field("p", rect, ("Lagrange", 1), unit=syu.pascal) + electrode = reference.electrode(rect) + darcy = Darcy("darcy", p, electrode.permeability(), viscosity) + cell = reference.half_cell( + rect, + electrode, + velocity=darcy.superficial_velocity, + density=density, + viscosity=viscosity, + ) + parameters = build_system(cell.model + darcy).parameters + assert parameters["electrolyte.viscosity"] is viscosity + assert parameters["electrode.porosity"] is electrode.porosity + assert "darcy.viscosity" not in parameters + + +def test_half_cell_errors(reference): + domain = Domain("cell", "interval") + with pytest.raises(ValueError, match="need a velocity"): + reference.half_cell(domain, mass_transfer=True) + with pytest.raises(ValueError, match="need a velocity"): + reference.half_cell(domain, outlet=("Outlet",)) + + +def test_the_couple_is_oxidized_or_reduced_on_discharge(reference): + """On discharge the current enters the solid at the current collector of a + half-cell that oxidizes its couple, and leaves it at that of one that + reduces it: the equations of the potentials, with only the current in + them, are opposite.""" + marked = reference.interval(num_cells=4) + residuals = [] + for oxidation in (True, False): + domain = Domain.from_mesh("cell", marked.mesh) + cell = RedoxHalfCell( + domain, + reference.IONS, + reference.RESERVOIR, + "V2", + "V3", + reference.electrode(domain), + oxidation_on_discharge=oxidation, + rate_constant=0.0, + formal_potential=0.0, + ) + scheme = BDF(1) + system = build_system( + cell.model, substitutions=scheme.substitutions(cell.model) + ) + disc = discretize([system], DomainMaps({domain: marked})) + scheme.set_step(1, 1.0, disc) + # Without gradients, rate or change, only the currents are left in + # the equations of the potentials + disc.set({cell.current: 0.5}) + problem = disc[system] + residuals.append([owned_residual(form) for form in problem.residual[:2]]) + # Each process may hold only one of the boundaries, so the checks are + # over all processes, which all reach the same verdict + comm = marked.mesh.comm + for oxidizing, reducing in zip(*residuals, strict=True): + largest = comm.allreduce(np.abs(oxidizing).max(initial=0.0), op=MPI.MAX) + mismatch = comm.allreduce( + np.abs(oxidizing + reducing).max(initial=0.0), op=MPI.MAX + ) + assert largest > 0 + assert mismatch <= 1e-12 * largest + + +def owned_residual(form: dolfinx.fem.Form) -> np.ndarray: + """A residual vector with its ghost contributions added to their owners, + restricted to the entries this process owns.""" + vector = dolfinx.fem.assemble_vector(form) + vector.scatter_reverse(dolfinx.la.InsertMode.add) + size = vector.index_map.size_local * vector.block_size + return vector.array[:size].copy() + + +def test_the_reservoir_names_ions_of_the_electrolyte(reference): + domain = Domain("cell", "interval") + with pytest.raises(ValueError, match="no ions"): + RedoxHalfCell( + domain, + reference.IONS, + {**reference.RESERVOIR, "Cl": 1.0}, + "V2", + "V3", + reference.electrode(domain), + oxidation_on_discharge=True, + rate_constant=2.5e-7, + formal_potential=-0.36, + ) + + +def test_post_processing(reference): + marked = reference.interval(num_cells=10) + domain = Domain.from_mesh("cell", marked.mesh) + c = Field("c", domain, ("Lagrange", 1)) + disc = discretize([c], DomainMaps({domain: marked})) + disc.set({c: lambda x: 1.0 + x[0] / reference.W}) + function = disc.function(c) + assert function_range(function) == pytest.approx((1.0, 2.0)) + points, values = evaluate_over_line(function, [0.0], [reference.W], 5) + np.testing.assert_allclose(values[:, 0], 1.0 + points[:, 0] / reference.W) + assert values.shape[1] == 1 + + +def test_the_couple_conserves_charge(reference): + """V3 + 2 e- <-> V2 does not conserve charge, nor V3 + e- <-> SO4.""" + domain = Domain("cell", "interval") + with pytest.raises(ValueError, match="does not conserve charge"): + reference.half_cell(domain, electrons=2) + with pytest.raises(ValueError, match="needs a reservoir"): + RedoxHalfCell( + domain, + reference.IONS, + reference.RESERVOIR, + "SO4", + "V3", + reference.electrode(domain), + oxidation_on_discharge=True, + rate_constant=2.5e-7, + formal_potential=-0.36, + ) + + +@pytest.mark.parametrize("missing", [(), ("H", "SO4")], ids=["none", "two"]) +def test_the_reservoir_lacks_exactly_one_ion(reference, missing): + """One ion, and only one, follows from electroneutrality.""" + domain = Domain("cell", "interval") + reservoir = {**reference.RESERVOIR, "SO4": 100.0} + reservoir = {k: v for k, v in reservoir.items() if k not in missing} + with pytest.raises(ValueError, match="all ions but one"): + RedoxHalfCell( + domain, + reference.IONS, + reservoir, + "V2", + "V3", + reference.electrode(domain), + oxidation_on_discharge=True, + rate_constant=2.5e-7, + formal_potential=-0.36, + ) diff --git a/tests/electrochemistry/half_cell/test_reference.py b/tests/electrochemistry/half_cell/test_reference.py new file mode 100644 index 0000000..8e1ab0c --- /dev/null +++ b/tests/electrochemistry/half_cell/test_reference.py @@ -0,0 +1,353 @@ +"""The models reproduce the forms and solutions of the RFBniCS examples.""" + +import basix.ufl +import dolfinx.fem +import dolfinx.fem.petsc +import numpy as np +import pytest +import sympy.physics.units as syu +import ufl +import ufl.algorithms + +from sumfem.discrete import DomainMaps, discretize +from sumfem.models import BoundaryFlux, Darcy +from sumfem.solvers import BDF, Simulation +from sumfem.symbolic import Clock, Domain, Field, Parameter, build_system + +IONS = ("V2", "V3", "H", "HSO4") + + +def fluid(domain): + """Viscosity and density of the examples, shared by flow and kinetics.""" + return ( + Parameter( + "electrolyte.viscosity", + domain, + default=6.963e-3, + unit=syu.pascal * syu.second, + ), + Parameter( + "electrolyte.density", domain, default=1399.8, unit=syu.kg / syu.m**3 + ), + ) + + +def unknowns(cell): + return [cell.phi_s, cell.phi_l, *(cell.concentrations[ion] for ion in IONS)] + + +def set_state(disc, cell, reference): + """A smooth state away from equilibrium, with a different previous step.""" + W = reference.W + disc.set( + { + cell.phi_s: lambda x: 0.01 * np.sin(np.pi * x[0] / W), + cell.phi_l: lambda x: 0.36 + 0.005 * x[0] / W, + cell.concentrations["V2"]: lambda x: ( + 300 * (1 + 0.1 * np.sin(2 * np.pi * x[0] / W)) + ), + cell.concentrations["V3"]: lambda x: ( + 300 * (1 - 0.1 * np.cos(np.pi * x[0] / W)) + ), + cell.concentrations["H"]: lambda x: 1000 * (1 + 0.05 * x[0] / W), + cell.concentrations["HSO4"]: lambda x: 1000 * (1 - 0.02 * x[0] / W), + } + ) + for ion in IONS: + field = cell.concentrations[ion] + previous = disc.function(field.previous(1)) + previous.x.array[:] = 0.99 * disc.function(field).x.array + + +def compare(disc, cell, problem, reference, dt, current, **options): + """Assemble both residuals, block by block, and compare them.""" + fields = unknowns(cell) + spaces = [disc.function_space(f) for f in fields] + tests = ufl.TestFunctions(ufl.MixedFunctionSpace(*spaces)) + marked = disc.maps.mesh_data(cell.phi_s.domain) + residual = reference.electrochemistry( + marked, + [disc.function(f) for f in fields], + [disc.function(cell.concentrations[ion].previous(1)) for ion in IONS], + tests, + dt, + current, + **options, + ) + expected = ufl.extract_blocks(residual) + for field, ours, theirs in zip(fields, problem.residual, expected, strict=True): + a = dolfinx.fem.assemble_vector(ours).array + b = dolfinx.fem.assemble_vector(dolfinx.fem.form(theirs)).array + scale = np.abs(b).max() + assert scale > 0, field.name + np.testing.assert_allclose( + a, b, rtol=1e-9, atol=1e-11 * scale, err_msg=field.name + ) + + +def backward_euler(cell): + """The scheme, and the system of one backward Euler step of the cell's model.""" + scheme = BDF(1) + return scheme, build_system( + cell.model, substitutions=scheme.substitutions(cell.model) + ) + + +def transient(disc, scheme, system, cell, dt, current): + """Bind the step's system, and set the step and current.""" + problem = disc[system] + scheme.set_step(1, dt, disc) + disc.set({cell.current: current}) + return problem + + +def test_1d_residual_matches_rfbnics(reference): + marked = reference.interval() + domain = Domain.from_mesh("cell", marked.mesh) + cell = reference.half_cell(domain) + scheme, system = backward_euler(cell) + disc = discretize( + [system], + DomainMaps({domain: marked}), + max_quadrature_degree=reference.QUADRATURE_DEGREE, + ) + set_state(disc, cell, reference) + problem = transient(disc, scheme, system, cell, dt=0.1, current=-0.5) + compare(disc, cell, problem, reference, dt=0.1, current=-0.5) + + +def test_2d_residual_matches_rfbnics(reference): + marked = reference.rectangle() + domain = Domain.from_mesh("cell", marked.mesh) + viscosity, density = fluid(domain) + electrode = reference.electrode(domain) + p = Field("p", domain, ("Lagrange", 1), unit=syu.pascal) + flow = Darcy("darcy", p, electrode.permeability(), viscosity) + cell = reference.half_cell( + domain, + electrode, + velocity=flow.superficial_velocity, + outlet=("Outlet",), + proton_flux=True, + density=density, + viscosity=viscosity, + ) + scheme, system = backward_euler(cell) + disc = discretize( + [system, p], + DomainMaps({domain: marked}), + max_quadrature_degree=reference.QUADRATURE_DEGREE, + ) + set_state(disc, cell, reference) + # A pressure drop of 1 kPa from inlet to outlet + disc.set({p: lambda x: 3e5 + 1e3 * (1 - x[1] / reference.H)}) + problem = transient(disc, scheme, system, cell, dt=0.5, current=0.5) + velocity = reference.darcy_velocity(marked, disc.function(p)) + compare( + disc, + cell, + problem, + reference, + dt=0.5, + current=0.5, + velocity=velocity, + membrane_flux=True, + outlet=True, + ) + + +def test_2d_pressure_matches_rfbnics(reference): + marked = reference.rectangle() + domain = Domain.from_mesh("cell", marked.mesh) + viscosity, _ = fluid(domain) + electrode = reference.electrode(domain) + p = Field("p", domain, ("Lagrange", 1), unit=syu.pascal) + inflow = Parameter("u_in", domain, default=1e-7, unit=syu.m / syu.second) + model = Darcy("darcy", p, electrode.permeability(), viscosity) + BoundaryFlux( + "inflow", p, inflow, ("Inlet",) + ) + disc = discretize( + [model], + DomainMaps({domain: marked}), + max_quadrature_degree=reference.QUADRATURE_DEGREE, + ) + disc.set({p: lambda x: 3e5 + 1e3 * x[0] * x[1] / (reference.W * reference.H)}) + problem = disc[model] + q = ufl.TestFunction(disc.function_space(p)) + expected = reference.darcy(marked, disc.function(p), q) + a = dolfinx.fem.assemble_vector(problem.residual).array + b = dolfinx.fem.assemble_vector(dolfinx.fem.form(expected)).array + np.testing.assert_allclose(a, b, rtol=1e-12, atol=1e-12 * np.abs(b).max()) + + +def test_3d_residual_matches_rfbnics(reference): + marked = reference.box() + domain = Domain.from_mesh("cell", marked.mesh) + viscosity, density = fluid(domain) + electrode = reference.electrode(domain, ("Electrode",)) + u = Field("u", domain, ("Lagrange", 2, (3,)), unit=syu.m / syu.second) + cell = reference.half_cell( + domain, + electrode, + channels=("InletChannel", "OutletChannel"), + velocity=u, + outlet=("Outlet",), + proton_flux=True, + density=density, + viscosity=viscosity, + ) + scheme, system = backward_euler(cell) + disc = discretize( + [system], + DomainMaps({domain: marked}), + max_quadrature_degree=reference.QUADRATURE_DEGREE, + ) + # The cube is 1 m wide, so scale the state of the 2.5 mm electrode + set_state(disc, cell, type(reference)(**{**vars(reference), "W": 1.0})) + disc.set({u: lambda x: (1e-4 * (1 + x[1]), 2e-5 * x[2], -1e-5 * x[0])}) + cell.balance_currents(disc) + ratio = float(disc.value(cell.area_ratio)) + assert ratio == pytest.approx(2.0) + problem = transient(disc, scheme, system, cell, dt=0.2, current=-0.5) + compare( + disc, + cell, + problem, + reference, + dt=0.2, + current=-0.5, + velocity=disc.function(u), + regions=True, + area_ratio=ratio, + membrane_flux=True, + outlet=True, + ) + + +def rfbnics_1d(marked, reference, schedule): + """The time loop of the 1D example of RFBniCS, on a blocked function.""" + mesh = marked.mesh + el = basix.ufl.element("Lagrange", mesh.basix_cell(), 1) + V = dolfinx.fem.functionspace(mesh, basix.ufl.blocked_element(el, shape=(6,))) + X = dolfinx.fem.Function(V) + X_old = dolfinx.fem.Function(V) + parameters = reference.PARAMETERS + for k, ion in enumerate(IONS, start=2): + X.sub(k).interpolate( + lambda x, c=parameters[f"c0_{ion}"]: np.full(x.shape[1], c) + ) + X.sub(1).interpolate(lambda x: np.full(x.shape[1], -parameters["E_formal"])) + X_old.x.array[:] = X.x.array + u = ufl.split(X) + dt = dolfinx.fem.Constant(mesh, 0.1) + current = dolfinx.fem.Constant(mesh, 0.0) + residual = reference.electrochemistry( + marked, u, ufl.split(X_old)[2:], ufl.TestFunctions(V), dt, current + ) + tdim = mesh.topology.dim + membrane = marked.entities(tdim - 1, "Membrane") + bcs = [ + dolfinx.fem.dirichletbc( + np.float64(parameters[f"c0_{ion}"]), + dolfinx.fem.locate_dofs_topological(V.sub(k), tdim - 1, membrane), + V.sub(k), + ) + for k, ion in enumerate(IONS, start=2) + ] + problem = dolfinx.fem.petsc.NonlinearProblem( + residual, + X, + bcs=bcs, + petsc_options_prefix="rfbnics_", + petsc_options={ + "snes_type": "newtonls", + "snes_linesearch_type": "bt", + "snes_atol": 1e-8, + "snes_rtol": 1e-12, + "snes_stol": 1e-14, + "snes_error_if_not_converged": True, + "ksp_type": "preonly", + "pc_type": "lu", + "pc_factor_mat_solver_type": "mumps", + # No potential is fixed, so the Jacobian is singular + "mat_mumps_icntl_24": 1, + "mat_mumps_icntl_25": 0, + }, + ) + for signed_current, steps in schedule: + current.value = signed_current + for count, step in steps: + dt.value = step + for _ in range(count): + problem.solve() + X_old.x.array[:] = X.x.array + return X, V + + +def test_1d_charge_and_discharge_match_rfbnics(reference): + marked = reference.interval(num_cells=30) + # The signed current density, and the number and size of the time steps + schedule = ( + (-0.5, ((3, 0.1), (2, 1.0))), + (0.0, ((2, 1.0),)), + (0.5, ((3, 1.0),)), + ) + X, V = rfbnics_1d(marked, reference, schedule) + + domain = Domain.from_mesh("cell", marked.mesh) + cell = reference.half_cell(domain) + inlet = cell.inlet_conditions("Membrane") + disc = discretize( + [cell.model, *inlet], + DomainMaps({domain: marked}), + max_quadrature_degree=reference.QUADRATURE_DEGREE, + ) + cell.initialize(disc) + clock = Clock(domain) + simulation = Simulation( + disc, + cell.model, + inlet, + clock=clock, + scheme=BDF(1), + petsc_options={"snes_atol": 1e-8, "snes_rtol": 1e-12, "snes_stol": 1e-14}, + ) + voltages, charge = [], [] + + def monitor(t, step): + voltages.append(cell.voltage(disc)) + charge.append(cell.state_of_charge(disc)) + + for signed_current, steps in schedule: + disc.set({cell.current: signed_current}) + for count, dt in steps: + disc.set({clock.dt: dt}) + simulation.run(num_steps=count, callback=monitor) + assert len(voltages) == 10 + assert float(disc.value(clock.t)) == pytest.approx(0.3 + 2.0 + 2.0 + 3.0) + values = {} + for k, field in enumerate(unknowns(cell)): + space, dofs = V.sub(k).collapse() + expected = dolfinx.fem.Function(space) + expected.x.array[:] = X.x.array[dofs] + ours = disc.function(field) + # Both are degree-one Lagrange on one interval, so order by vertex + order_ours = np.argsort(ours.function_space.tabulate_dof_coordinates()[:, 0]) + order_theirs = np.argsort( + expected.function_space.tabulate_dof_coordinates()[:, 0] + ) + values[field.name] = ours.x.array[order_ours], expected.x.array[order_theirs] + # Only differences of potentials are defined, so measure both potentials + # from the electrolyte potential at the membrane + for name in ("phi_s", "phi_l"): + values[name] = tuple(values[name][i] - values["phi_l"][i][-1] for i in range(2)) + for name, (ours, theirs) in values.items(): + scale = max(np.abs(theirs).max(), 1e-3) + np.testing.assert_allclose( + ours, theirs, rtol=1e-7, atol=1e-9 * scale, err_msg=name + ) + # Charging turns V3 into V2, and discharging turns it back + assert charge[4] > 0.5 + assert charge[-1] < charge[6] + assert charge[6] == pytest.approx(charge[4], rel=1e-6) + assert np.isfinite(voltages).all() diff --git a/tests/electrochemistry/test_electrolyte.py b/tests/electrochemistry/test_electrolyte.py index 7c9a790..0990f15 100644 --- a/tests/electrochemistry/test_electrolyte.py +++ b/tests/electrochemistry/test_electrolyte.py @@ -13,7 +13,7 @@ NernstPlanck, PorousElectrode, ) -from sumfem.symbolic import Field, build_system +from sumfem.symbolic import Domain, Field, build_system CONCENTRATION = syu.mol / syu.m**3 @@ -152,3 +152,10 @@ def test_a_vanadium_electrode_has_consistent_units(domain): ) system = build_system(model, check_units=True) assert {f.name for f in system.unknowns} >= {"phi_s", "phi_l", "c2", "c3", "h"} + + +def test_the_balanced_ion_needs_a_charge(): + domain = Domain("cell", "interval") + c = Field("c", domain, ("Lagrange", 1)) + with pytest.raises(ValueError, match="has no charge"): + Electrolyte("e", domain, [Ion("A", 1, 1e-9, c), Ion("N", 0, 1e-9)]) diff --git a/tests/symbolic/test_terminals.py b/tests/symbolic/test_terminals.py index 9efe556..ee4cb3e 100644 --- a/tests/symbolic/test_terminals.py +++ b/tests/symbolic/test_terminals.py @@ -11,6 +11,7 @@ KnownCoefficient, Parameter, as_expr, + value_at_defaults, ) from sumfem.symbolic.elements import as_element from sumfem.symbolic.terminals import FieldRole, FieldTerminal, field_of @@ -164,6 +165,19 @@ def test_the_default_of_a_parameter_in_si(domain): np.testing.assert_allclose(ratio.default_si, [1.0, 3.0]) +def test_a_default_with_units_is_converted_to_the_unit(domain): + """A quantity becomes a value in the declared unit, so that constants can + carry their units; another dimension is rejected.""" + length = Parameter("length", domain, default=5 * syu.centimeter, unit=syu.mm) + assert length.default == pytest.approx(50.0) + assert length.default_si == pytest.approx(0.05) + diffusivity = syu.m**2 / syu.s + rate = Parameter("D", domain, default=0.36 * syu.cm**2 / syu.hour, unit=diffusivity) + assert rate.default_si == pytest.approx(1e-8) + with pytest.raises(ValueError, match="does not have the dimension"): + Parameter("bad", domain, default=1.0 * syu.m / syu.s, unit=diffusivity) + + def test_scalar_default_sets_every_component(domain): np.testing.assert_allclose( Parameter("g", domain, shape=(2,), default=3.0).default, [3, 3] @@ -222,3 +236,30 @@ def test_as_expr(domain): def test_parameter_names_may_be_qualified(domain): assert Parameter("solid.sigma", domain).name == "solid.sigma" + + +def test_a_zero_quantity_is_a_default(domain): + """Sympy reduces a zero quantity to the number zero, without units.""" + parameter = Parameter("p", domain, default=0.0 * syu.pascal, unit=syu.pascal) + assert parameter.default_si == 0.0 + + +def test_an_expression_of_parameters_has_a_value_at_their_defaults(domain): + mu = Parameter("mu", domain, default=2.0) + v = Parameter("v", domain, default=[1.0, 3.0], shape=(2,)) + assert value_at_defaults(mu * ufl.sqrt(ufl.dot(v, v))) == pytest.approx( + 2.0 * np.sqrt(10.0) + ) + np.testing.assert_allclose(value_at_defaults(mu * v), [2.0, 6.0]) + + +def test_a_value_at_defaults_needs_parameters_with_defaults(domain): + free = Parameter("free", domain) + with pytest.raises(ComponentError, match="'free'"): + value_at_defaults(2.0 * free) + u = Field("u", domain, ("Lagrange", 1)) + with pytest.raises(ComponentError, match="'u'"): + value_at_defaults(u + 1.0) + x = ufl.SpatialCoordinate(domain) + with pytest.raises(ComponentError, match="spatial coordinate"): + value_at_defaults(x[0] * Parameter("a", domain, default=1.0)) diff --git a/tests/test_api_docs.py b/tests/test_api_docs.py index 77fb05d..0ed9a82 100644 --- a/tests/test_api_docs.py +++ b/tests/test_api_docs.py @@ -23,6 +23,15 @@ def documented(page: pathlib.Path) -> set[str]: return names +def resolve(name: str) -> object: + """The object a documented full name refers to, or ``None``.""" + module, _, attribute = name.rpartition(".") + try: + return getattr(importlib.import_module(module), attribute, None) + except ImportError: + return None + + @pytest.mark.parametrize( "module, page", [ @@ -39,12 +48,18 @@ def documented(page: pathlib.Path) -> set[str]: ], ) def test_every_public_name_is_on_its_api_page(module, page): - public = importlib.import_module(module).__all__ + imported = importlib.import_module(module) names = documented(API / f"{page}.md") - # A name may be documented under its module or a package that re-exports it + # A name may be documented under its module, a package that re-exports + # it, or the module that defines it, as a constant must be, whose + # docstring autodoc finds only in the source of the module that defines it parts = module.split(".") owners = [".".join(parts[: k + 1]) for k in range(1, len(parts))] + objects = [resolve(name) for name in names] missing = sorted( - n for n in public if not any(f"{owner}.{n}" in names for owner in owners) + n + for n in imported.__all__ + if not any(f"{owner}.{n}" in names for owner in owners) + and not any(o is getattr(imported, n) for o in objects if o is not None) ) assert not missing, f"{missing} of {module} are missing from docs/api/{page}.md"