[1]Thomas Roy

1]Lawrence Livermore National Laboratory - 7000 East Ave, Livermore, CA 94550,
USA

# Topology Optimization for the Full-Cell Design of Porous Electrodes in
Electrochemical Energy Storage Devices

Hanyu Li Giovanna Bucci Nicholas W. Brady Nicholas R. Cross Victoria M.
Ehlinger Tiras Y. Lin Miguel Salazar de Troya Daniel Tortorelli Marcus A.
Worsley<EMAIL_ADDRESS>[

###### Abstract

In this paper, we introduce a density-based topology optimization framework to
design porous electrodes for maximum energy storage. We simulate the full cell
with a model that incorporates electronic potential, ionic potential, and
electrolyte concentration. The system consists of three materials, namely pure
liquid electrolyte and the porous solids of the anode and cathode, for which
we determine the optimal placement. We use separate electronic potentials to
model each electrode, which allow interdigitated designs. As the result, a
penalization is required to ensure that the anode and cathode do not touch,
i.e. causing a short circuit. We compare multiple 2D designs generated for
different fixed conditions, e.g. material properties. A 3D design with complex
channel and interlocking structure is also created. All optimized designs are
far superior to the traditional monolithic electrode design with respect to
energy storage metrics. We observe up to 750% increase in energy storage for
cases with slow effective ionic diffusion within the porous electrode.

###### keywords:

Topology optimization, Electrochemistry, batteries, supercapacitors

##  1 Introduction

Electrochemical energy storage devices provide a shift away from fossil fuels
by enabling electric vehicles and supporting the adoption of intermittent
renewable energy sources (Chu and Majumdar, 2012; Chu et al, 2016; Gür, 2018).
Batteries and capacitors are examples of such devices that are ubiquitous in
modern technologies. And thus improving their performance is crucial for the
green energy transition.

The electrochemical charge storage mechanisms can be broadly grouped into two
types: charge separation and Faradaic charge-transfer reactions (Simon et al,
2014). In the former, electrostatic double-layer capacitors (EDLCs) store
energy by forming an electric double layer at the interface between the
electrodes and electrolyte. This is different from the latter wherein
electrons are transferred via reduction and oxidation (i.e. _redox_)
reactions. Such reactions are the principal energy storage mechanism of
batteries, which utilize slow kinetics and therefore lower power density but
enable much higher energy density than EDLCs. As a middle ground,
pseudocapacitors use fast redox kinetics, resulting in smaller energy density
than batteries but larger than EDLCs while also demonstrating high power
density (Simon et al, 2014).

Batteries and capacitors share a common design: two porous electrodes, namely
the anode and the cathode, are immersed in a liquid electrolyte and separated
by a membrane or simply an electrolyte-filled gap. The traditional monolithic
porous electrode cell illustrated in Fig. 1 consists of two slabs described by
a uniform porosity and surface area to volume ratio. The electrodes are
typically assemblies of micron-scale conductive particles. Electrochemical
reactions occur over the interface between the electrode-particles and the
electrolyte. As such, cell-design criteria primarily focus on maximizing the
surface area to volume ratio while minimizing resistances related to electron
or ion transport (Wang et al, 2008), as these criteria enable higher currents,
power outputs, and energy storage.

Current collector (positive)Current collector
(negative)CathodeElectrolyte/separatorAnode Power
---
source /
load
$e^{-}$$e^{-}$ Figure 1: Diagram of a traditional monolithic electrochemical
energy storage device.

Simply maximizing electrochemical surface area, in the attempt to increase
energy density, often deteriorates ionic transport within the electrode by
replacing the pore network for ion transport with an extremely tortuous path
(Wang et al, 2008). Such paths restrict ionic penetration into the monolithic
electrodes and thus reduce their effectiveness. Fortunately, promising
alternatives to these monolithic designs are emerging; e.g. electrode designs
with geometric features, multiple lengthscales and/or materials, and spatially
varying porosity (Beck et al, 2021a; Trembacki et al, 2019; Zekoll et al,
2018; Gao et al, 2023; Yang et al, 2022). These new designs are made possible
thanks to new advanced manufacturing processes (Chu et al, 2021; Soares et al,
2021; Chandrasekaran et al, 2023). It is this nearly complete control over the
electrode design that raises the question of what these designs should be.

One approach to increase the energy density consists in increasing the
thickness of conventional slab electrodes and reducing the number of cells
needed in a stack. This approach allows to reduce the weight associated with
non-active materials such as current collectors, foils and separators (Heubner
et al, 2021). However, the electrode areal capacity is limited by the ion
transport within the tortuous electrode structure, and typical commercial-
electrodes have a thickness in the range of 100-150 . (Heubner et al, 2021; Xu
et al, 2019; Kuang et al, 2019). Further increasing the electrode thickness
(and the areal capacity) leads to a point of diminishing returns in volumetric
energy density for charge/discharge rates relevant to electric-vehicle
applications (Singh et al, 2015; Zheng et al, 2012; Gallagher et al, 2015).

Thicker electrodes with three-dimensional features, in contrast to
conventional slabs, have been previously fabricated and tested (Long et al,
2004, 2020), demonstrating the importance of providing structure and low-
tortuosity transport pathways that facilitate ion transport in the
electrolyte. A categorization of the various three-dimensional designs has
been proposed by Hung et al (2022): the electrodes may be individually shaped
and separated by a planar separator, which is either solid or porous, or
completely intertwined with one another. In particular, an interdigitated
design, where alternating planar fins or rods protruding from the cathode and
anode separated by the electrolyte, has demonstrated increased power and
energy performance (Hung et al, 2021; Zadin et al, 2011). Similar three-
dimensional interpenetrating structures have been investigated for
microbatteries and wearable devices applications (Lyu et al, 2021; Zeng et al,
2014; Wang et al, 2017). The non-planar structures improve performance via
increased surface area to volume ratios to promote reaction and increased
ionic diffusion through complex pathways to promote mass transport (Trembacki
et al, 2019; Chen et al, 2020), therefore enabling high power and fast
charging lithium-ion batteries (Chen et al, 2020; Zekoll et al, 2018).
However, because the current distribution and concentration profiles in these
interdigitated and interpenetrating electrodes are highly non-uniform (Hart et
al, 2003; Trembacki et al, 2019), their performance becomes strongly dependent
on the electrode and electrolyte properties, e.g. specific surface area and
conductivity, thus motivating the need for a systematic design optimization
paradigm.

Computational techniques such as topology optimization offer a promising
opportunity to design optimal porous electrodes. Density-based topology
optimization was initially formulated as a mass distribution problem in which
the volume fraction field is optimized to maximize the stiffness of a linear
elastic structure subject to a mass constraint Bendsøe (1989). Topology
optimization has since been adapted to design electrochemical devices. Yaji et
al (2018); Chen et al (2019); Lin et al (2022) design the channels that
transport the electrolyte fluid to the porous electrodes in redox flow
batteries. Behrou et al (2019); Wang et al (2023); Qu et al (2023) similarly
design channels in fuel cells to minimize electrical and pumping power losses
while maximizing power density. Heterogeneous porosity fields have been
designed for lithium-ion batteries (Ramadesigan et al, 2010; Golmon et al,
2012; Xue et al, 2015), and redox flow batteries (Beck et al, 2021b).

Topology optimization has been used in Roy et al (2022) to design redox porous
electrodes and EDLC electrodes. This study, which only considers a half cell,
i.e. a single electrode, generates designs for a wide range of fixed
dimensionless groups encapsulating material parameters, electrode length
scale, and operating conditions. The best performance improvement over
monolithic designs was achieved when ionic transport within the porous
material is slow. Batista et al (2023) used a similar optimization framework
to design a supercapacitor electrode, which was then fabricated using additive
manufacturing. The optimized electrode exhibited nearly twice the volumetric
capacitance over a lattice structured electrode. The models used to predict
the response in these previous studies assumed no concentration gradients; the
model in a more recent optimization study by Alizadeh et al (2023) included
mass accumulation and entropy effects. In this work, we optimize both
electrodes simultaneously using a full-cell model that includes mass
transport.

A fundamental flaw with density-based topology optimization is the use of the
material volume fraction field as the control. As such, regions in the
optimized designs contain mixtures of materials, which is not possible, i.e.
we must have so-called $0-1$ designs in which a distinct material resides at
each location. The size of these mixed “intermediate” material regions is
minimized by using clever penalization schemes in the topology optimization
problem formulation. For structural optimization, the penalization is
relatively straightforward because of the explicit trade-off between the
structural stiffness and mass cf. e.g. the Simple Isotropic Material with
Penalization (SIMP) (Bendsøe, 1989). For electrochemical energy storage
systems, the trade-off between ion/electron transport and chemical activity is
less obvious, thus special care is required when devising penalization
schemes. Schemes that concurrently increase ohmic loss and reduce energy
storage in intermediate material regions have proven effective (Roy et al,
2022; Batista et al, 2023). Developing effective penalization schemes for
topology optimization can be very tedious. However those that conquer this
battle are rewarded with a very effective design tool.

In this manuscript, we use topology optimization to design full-cell
electrochemical energy storage devices. In Sect. 2, we review topology
optimization concepts, and describe the boundary propagation method required
for the two-electrode optimization. Next, in Sect. 3, we introduce the
simulation model, i.e. the Nernst-Planck equation for a binary electrolyte
with redox reactions and electrostatic double layer capacitance.
Implementation and solver details are provided in Sect. 4. Following that, in
Sect. 5, the optimization problem is formulated and solved to generate designs
for various key dimensionless groups encapsulating device length scale and
material parameters such as conductivity and reaction rate. We also compare
the performances of two optimized designs against traditional monolithic
electrodes to demonstrate effectiveness of the optimization. Next, we showcase
a 3D design, possessing complex channels and interlocking features. We
conclude in Sect. 6 with some final discussion.

##  2 Topology optimization

Current collector $\Gamma_{\mathrm{a}}$Current collector
$\Gamma_{\mathrm{c}}$Cathode ($\chi=1$)Electrolyte/ Separator($\chi=0$)Anode
($\chi=1$) Figure 2: Two-dimensional design space of a full cell with porous
electrodes.

Our topology optimization problem concerns the material indicator field $\chi$
in cell domain $\Omega$ defined such that $\chi(\mathbf{x})=1$ indicates that
the solid porous electrode material occupies $\mathbf{x}$, whereas
$\chi(\mathbf{x})=0$ indicates that the separator or pure liquid electrolyte
occupies $\mathbf{x}$, cf. Fig. 2. Over this domain we solve the model
equations $F(\chi,\Phi_{1},\Phi_{2},c)=0$ to be introduced in Sect. 3, for the
electronic potential $\Phi_{1}$, ionic potential $\Phi_{2}$, and salt
concentration $c$. These pieces are used to define the porous electrodes
topology optimization problem, i.e.

$\displaystyle\min_{\chi\colon\Omega\to\\{0,1\\}}\theta_{0}(\chi)=\int_{\Omega}\pi(\chi,\Phi_{1},\Phi_{2},c)\mathop{}\\!\mathrm{d}V,$
(1) $\displaystyle\text{s.t.
}\theta_{i}(\chi)=\int_{\Omega}g_{i}(\chi,\Phi_{1},\Phi_{2},c)\mathop{}\\!\mathrm{d}V\leq
0,\;i=1,\dots,m,$ $\displaystyle\text{and}\ F(\chi,\Phi_{1},\Phi_{2},c)=0.$

Here, $\theta_{0}$ is the cost function to be minimized and $\theta_{i}$ are
the $m$ design constraint functions, with local functions $\pi$ and $g_{i}$ to
be defined in Sect. 5.1.

Since the electrode phase is comprised of separate anode and cathode
materials, a three-phase (three-material) optimization formulation would be
natural (Sigmund and Torquato, 1997; Watts and Tortorelli, 2016). This three-
phase approach requires a second material indicator field that would delineate
the anode and cathode phases. However, since cathode and anode are in contact
with different boundaries $\Gamma_{\mathrm{c}}$ and $\Gamma_{\mathrm{a}}$ (see
Fig. 2), we instead use a boundary propagation method (Sá et al, 2022) for
this delineation, cf. Sect. 2.2. Thus, we optimize with respect to a single
field $\chi$.

In the context of density-based topology optimization (Bendsøe, 1989), the
discrete nature of the design field $\chi$ prevents the application of
gradient-based optimization algorithms. Because of this, and at the risk of
not obtaining $0-1$ designs, we replace $\chi\colon\Omega\to\\{0,1\\}$ with a
continuous volume fraction variable $\rho\colon\Omega\to[0,1]$. Additionally,
special care is required to ensure that optimization problem is well-posed.
These concerns are remedied using filtering detailed in Sect. 2.1 and material
penalization described in Sect. 5.1.

###  2.1 Filtering

Topology optimization problems are often ill-posed. Their “solutions” are
characterized by non-converging sequences of structures with highly
oscillatory material distributions. One option to obtain a well-posed problem
imposes a feature length scale on the design by filtering $\rho$ to obtain the
smoothed field $\tilde{\rho}$, which is used to define the structures of
material distribution (Bruns and Tortorelli, 2001; Bourdin, 2001). In our
problem, we use the diffusion-reaction PDE filter described in Lazarov and
Sigmund (2011) wherein $\tilde{\rho}$ solves

$\displaystyle-r^{2}\nabla^{2}\tilde{\rho}+\tilde{\rho}$
$\displaystyle=\rho\quad$ $\displaystyle\text{in }\Omega,$ (2) $\displaystyle
r^{2}\nabla\tilde{\rho}\cdot\mathbf{n}$ $\displaystyle=0\quad$
$\displaystyle\text{on }\partial\Omega,$

for the given design $\rho$. Fine scale oscillations in $\rho$ are mollified
in $\tilde{\rho}$; more mollification occurs as the filter radius $r$
increases.

Unfortunately, the filtered design $\tilde{\rho}$ contains large transitional
regions where $\tilde{\rho}\neq\\{0,1\\}$ and our original design problem only
allows for $\chi(\mathbf{x})\in\\{0,1\\}$ regions, i.e. $0-1$ designs. To
remedy this, the filtered distribution $\tilde{\rho}$ is thresholded to define
the field $\bar{\rho}$ such that $\bar{\rho}(\mathbf{x})=\\{0,1\\}$ throughout
$\Omega$. We can do this by equating $\bar{\rho}=H(\tilde{\rho}-0.5)$ where
$H$ is the unit Heaviside function. However, the lack of differentiability of
$H$ is not conducive to our gradient-based optimization. As such, we replace
$H$ with a smoothed approximation so that

$\bar{\rho}=H_{b,k}(\tilde{\rho})=\frac{\tanh(bk)+\tanh\left(b(\tilde{\rho}-k)\right)}{\tanh(bk)+\tanh\left(b(1-k)\right)}.$
(3)

Here, $b$ controls the sharpness of the $0-1$ transition and $k$ controls the
cutoff such that

$\lim_{b\rightarrow\infty}H_{b,k}(\tilde{\rho})=H(\tilde{\rho}-k)$

(Guest et al, 2004; Sigmund, 2007; Sigmund and Maute, 2013; Xu et al, 2010).
In this way, we obtain “nearly” $0-1$ designs such that
$\bar{\rho}(\mathbf{x})\approx\\{0,1\\}$ throughout $\Omega$. Here, we use
$b=4$ and $k=0.5$ for a mild transition. Although mild, it combines well with
the continuous boundary propagation method to identify the extent of the
cathode and anode regions. Aggressive thresholding, i.e. larger $b$ values,
leads to premature convergence to poor local optima. Because of this,
continuation techniques on $b$ are often employed to improve optimization
convergence and obtain nearly $0-1$ designs. In this work, we do not use
continuation approaches, since the combination of mild thresholding and the
interpolation described in 5.1 lead to nearly $0-1$ designs.

###  2.2 Continuous boundary propagation

Figure 3: Physical density $\bar{\rho}$ (top) with its corresponding $\beta$
(middle) and anode indicator $\mathcal{I}_{a}$ (bottom) for the initial design
(left) and optimized design (right).

We can use $\bar{\rho}(\mathbf{x})=\\{0,1\\}$ to determine if electrolyte or
electrode resides at $\mathbf{x}$, but we do not know if the electrode is the
anode or cathode. This is important because the anode and cathode may be
manufactured with different materials and because we need to separate the
anode and cathode regions to avoid a short circuit. Both of these concerns are
resolved by using the continuous boundary propagation approach (Sá et al,
2022).

The boundary propagation approach consists of solving the following diffusion-
reaction problem cf. Fig. 3 for $\beta$:

$\displaystyle\bar{\rho}\nabla\cdot\nabla\beta-(1-\bar{\rho})\beta$
$\displaystyle=0\quad$ $\displaystyle\text{in }\Omega,$ (4)
$\displaystyle\beta$ $\displaystyle=1\quad$ $\displaystyle\text{on
}\Gamma_{\mathrm{a}},$ $\displaystyle\beta$ $\displaystyle=-1\quad$
$\displaystyle\text{on }\Gamma_{\mathrm{c}},$
$\displaystyle\nabla\beta\cdot\mathbf{n}$ $\displaystyle=0\quad$
$\displaystyle\text{on
}\partial\Omega/\Gamma_{\mathrm{c}}\cup\Gamma_{\mathrm{a}}.$

In electrode regions ($\bar{\rho}=1$), the equation reduces to
$\nabla\cdot\nabla\beta=0$, and in electrolyte regions ($\bar{\rho}=0$), it
reduces to $\beta=0$. By prescribing $\beta=1$ and $\beta=-1$ on anode and
cathode current collector boundaries, namely $\Gamma_{a}$ and $\Gamma_{c}$, we
are assured that $0\leq\beta\leq 1$ in the anode and $-1\leq\beta\leq 0$ in
the cathode. The resulting $\beta$ is transformed into differentiable
indicator functions for the anode and cathode regions such that

$\displaystyle\mathcal{I}_{\mathrm{a}}(\beta)=$ $\displaystyle
H_{b,k}\left(\frac{\beta+1}{2}\right),$ (5)
$\displaystyle\mathcal{I}_{\mathrm{c}}(\beta)=$ $\displaystyle
H_{b,k}\left(\frac{-\beta+1}{2}\right),$

with $k=0.5$.

Examples of $\bar{\rho}$ distributions, their corresponding $\beta$
propagation, and anode indicator computed with $b=100$ are illustrated in Fig.
3. To summarize, $\bar{\rho}$ distinguishes between solid electrode and liquid
electrolyte regions, and $\beta$ produces the indicator fields
$\mathcal{I}_{\mathrm{a}}$ and $\mathcal{I}_{\mathrm{c}}$ to distinguish
between anode and cathode regions. The indicator fields are also used to
ensure model consistency as discussed in Sect. 3.3.

##  3 Governing equations

In this section, we introduce the governing equations to model electrochemical
energy storage devices with both redox and capacitive charging effects
(Johnson and Newman, 1971; Doyle et al, 1993; Newman and Tobias, 1962; Newman
and Tiedemann, 1975). Extending the two-potential system described in our
previous study (Roy et al, 2022), we solve the Nernst–Planck equation for the
ion concentration fields, which affect the charge-transfer and adsorption
kinetics, and the electrolyte conductivity. Therefore, accounting for its mass
transport improves the predictability of our designs. Battery models may also
include ionic diffusion in the solid particles (Doyle et al, 1993; Brosa
Planella et al, 2022), however, for simplicity we ignore solid diffusion,
which is a reasonable assumption when the electrode particles are of the order
of 1 $\mu$m in size. Additionally, we assume that the anode and cathode have
identical material properties.

###  3.1 Physical model

We consider a binary electrolyte in which the cation participates in the redox
reaction that produces current density $i_{n}$. In the battery, this reaction
has the form

$s_{+}A^{z_{+}}+ne^{-}\rightleftharpoons s_{+}A,$ (6)

where $n$ is the number of electrons transferred, $z_{+}$ is the cation charge
number, and $s_{+}=n/z_{+}$ is the cation stoichiometric number. Both cations
and anions are adsorbed by double-layer charging, which produces current
density $i_{c}$. The Nernst–Planck equation for the positive ($c_{+}$) and
negative ($c_{-}$) charged species concentrations are as follows:

$\displaystyle\frac{\partial\left(\epsilon c_{+}\right)}{\partial t}=$
$\displaystyle-\nabla\cdot\mathbf{N}_{+}+\frac{dq_{+}}{dq}\frac{1}{z_{+}F}ai_{c}+\frac{s_{+}}{nF}ai_{n}$
$\displaystyle\text{in }\Omega,$ (7)
$\displaystyle\frac{\partial\left(\epsilon c_{-}\right)}{\partial t}=$
$\displaystyle-\nabla\cdot\mathbf{N}_{-}+\frac{dq_{-}}{dq}\frac{1}{z_{-}F}ai_{c}$
$\displaystyle\text{in }\Omega,$

where $\mathbf{N}_{i}$ is the ionic flux for species $\textit{i}=\pm$,
$\epsilon$ is porosity, $F$ is the Faraday constant (96485), and $a$ is the
specific surface area, which equates zero in the liquid electrolyte region.
The variation in the amount of adsorbed cations/anions is captured by
$\dfrac{dq_{i}}{dq}$, where $q_{i}$ is the adsorbed cation/anion surface
charge density in the double layer, and $q$ is the total surface charge
density in the double layer. Assuming the absence of cosorption of anions and
cations, $\frac{dq_{+}}{dq}$ = 1 and $\frac{dq_{-}}{dq}$ = 0 for a cation
responsive electrode, and $\frac{dq_{+}}{dq}$ = 0 and $\frac{dq_{-}}{dq}$= 1
for an anion responsive electrode (Johnson and Newman, 1971).

We now define constitutive equations for $i_{n}$, $i_{c}$, and
$\mathbf{N}_{i}$. The redox and capacitive current densities are defined as

$\displaystyle i_{n}$
$\displaystyle=i_{0}\left[\exp\left(\frac{\alpha_{a}F}{RT}\eta\right)-\exp\left(\frac{-\alpha_{c}F}{RT}\eta\right)\right],$
(8) $\displaystyle i_{c}$
$\displaystyle=C_{d}\frac{\partial\left(\Phi_{1}-\Phi_{2}\right)}{\partial
t},$

where $\alpha_{a}$ and $\alpha_{c}$ are the anodic and cathodic charge
transfer coefficients, respectively. For simplicity, we assume
$\alpha=\alpha_{a}=\alpha_{c}$. Here, the exchange current density is
$i_{0}=Fk_{\mathrm{rxn}}c^{\alpha}$, where $k_{\mathrm{rxn}}$ is the reaction
rate and $c$ is the salt concentration defined below. The overpotential
$\eta=\Phi_{1}-\Phi_{2}-U_{0}$ is obtained from the electronic potential
$\Phi_{1}$, the ionic potential $\Phi_{2}$, and the equilibrium potential
$U_{0}$. $R$ (8.3145J.K^-1.mol^-1) is the ideal gas constant, $T$ is the
temperature, and $C_{d}$ is the specific capacitance. Note that battery models
typically treat $U_{0}$ as a function of solid surface concentration (Doyle et
al, 1993); for simplicity, we assume $U_{0}$ to be constant. Therefore, the
Nernstian loss, which is defined by the change in the equilibrium potential as
a function of the concentration, is not modeled. To account for the effect of
concentration on ion adsorption, we define
$C_{d}=C_{0}\left(\frac{c}{c_{0}}\right)^{\alpha}$. This nonstandard addition
to the model prevents nonphysical extreme capacitive current density in the
designs where high reaction rates could lead to ion depletion.

The ionic flux for an individual species is calculated by

$\mathbf{N}_{i}=-z_{i}u_{i}Fc_{i}\nabla\Phi_{2}-D_{i}\nabla c_{i}.$ (9)

In the above, $u_{i}$ is the effective ion mobility and $D_{i}$ is the
effective molecular diffusion coefficient of the species, which satisfy the
Nernst-Einstein relationship as

$D_{i}=RTu_{i}.$ (10)

The stoichiometric coefficient and ion charge number are related by

$z_{i}=\frac{n}{s_{i}}.$ (11)

Both relations are useful for nondimensionalization of the governing system.
Assuming electroneutrality, we have the constraint

$\sum_{i}z_{i}c_{i}=0,$ (12)

which is also equivalent to

$\sum_{i}z_{i}\nu_{i}=0,$ (13)

where $\nu_{i}$ is the number of moles of cations or anions per mole of
electrolyte. For a binary electrolyte, we can then define the salt
concentration as

$c=\frac{c_{+}}{\nu_{+}}=\frac{c_{-}}{\nu_{-}}.$ (14)

Considering electroneutrality, equations (7) are combined into one
conservation equation for the salt concentration (see Newman and Thomas-Alyea
(2012) for the derivation of the salt concentration equation). Combining that
with the charge conservation in the electronic and ionic potentials, we obtain
the governing field equations that solve for $\Phi_{1}$, $\Phi_{2}$, and $c$
(Newman and Tiedemann, 1975):

$\displaystyle-\nabla\cdot(\sigma\nabla\Phi_{1})$
$\displaystyle=-a\left(i_{n}+i_{c}\right)$ $\displaystyle\text{in }\Omega,$
(15a) $\displaystyle-\nabla\cdot(\kappa\nabla\Phi_{2})$ $\displaystyle-
z_{+}\nu_{+}F\nabla\cdot\left((D_{+}-D_{-})\nabla c\right)$
$\displaystyle=a\left(i_{n}+i_{c}\right)$ $\displaystyle\text{in }\Omega,$
(15b) $\displaystyle\frac{\partial(\epsilon c)}{\partial t}-\nabla\cdot($
$\displaystyle D\nabla c)=\frac{s_{+}}{Fn\nu_{+}}t_{-}ai_{n}$
$\displaystyle+\frac{1}{Fz_{+}\nu_{+}}$
$\displaystyle\left(t_{-}\frac{dq_{+}}{dq}-t_{+}\frac{dq_{-}}{dq}\right)ai_{c}$
$\displaystyle\text{in }\Omega,$ (15c) subject to the boundary and initial
conditions: $\displaystyle\Phi_{1}$ $\displaystyle=0\quad$
$\displaystyle\text{on }\Gamma_{\mathrm{a}},$ (15d) $\displaystyle\Phi_{1}$
$\displaystyle=\zeta t\quad$ $\displaystyle\text{on }\Gamma_{\mathrm{c}},$
(15e) $\displaystyle\nabla\Phi_{1}\cdot\mathbf{n}$ $\displaystyle=0\quad$
$\displaystyle\text{on
}\partial\Omega\setminus(\Gamma_{\mathrm{a}}\cup\Gamma_{\mathrm{c}}),$ (15f)
$\displaystyle\nabla\Phi_{2}\cdot\mathbf{n}$ $\displaystyle=0\quad$
$\displaystyle\text{on }\partial\Omega,$ (15g) $\displaystyle\nabla
c\cdot\mathbf{n}$ $\displaystyle=0\quad$ $\displaystyle\text{on
}\partial\Omega,$ (15h) $\displaystyle\Phi_{1}(\mathbf{x},t=0)$
$\displaystyle=0\quad$ $\displaystyle\text{in }\Omega,$ (15i)
$\displaystyle\Phi_{2}(\mathbf{x},t=0)$ $\displaystyle=0\quad$
$\displaystyle\text{in }\Omega,$ (15j) $\displaystyle c(\mathbf{x},t=0)$
$\displaystyle=c_{0}\quad$ $\displaystyle\text{in }\Omega,$ (15k)

for $t\in(0,T]$, where $\sigma$ is the effective electronic conductivity,
$\kappa$ is the effective ionic conductivity, and $c_{0}$ is the initial
concentration. The dimensionless quantity $t_{+}$ denotes the transference
number of the cation; the transference number of the anion is related by
$t_{-}=1-t_{+}$. The boundary conditions describe a linearly increasing
potential in time on the cathode current collector, with scan rate $\zeta$.

In the electrode region, the effective transport properties, namely effective
electronic conductivity $\sigma$, effective ionic conductivity $\kappa$, and
the effective salt diffusion coefficient or ambipolar diffusivity $D$ are
influenced by the tortuosity of the pores. The Bruggeman correlation
(Bruggeman, 1935) quantifies tortuosity as a power law of porosity with
exponent $\frac{3}{2}$ for packed spheres, and thus the effective transport
properties are calculated as

$\sigma=(1-\epsilon)^{\frac{3}{2}}\sigma_{0},\quad\kappa=\epsilon^{\frac{3}{2}}\kappa_{0},\quad
D=\epsilon^{\frac{3}{2}}D_{0},$ (16)

where the $0$ subscript represents the bulk value. Specifically, $\kappa_{0}$
and $D_{0}$ refer to the ionic conductivity and diffusion coefficient in pure
electrolyte ($\epsilon=1$), and $\sigma_{0}$ refers to the electronic
conductivity in non-porous solid ($\epsilon=0$). Note that $D_{0}$ and
$\kappa_{0}$ are computed from the individual ion species as

$D_{0}=D_{0_{+}}D_{0_{-}}\frac{z_{+}-z_{-}}{z_{+}D_{0_{+}}-z_{-}D_{0_{-}}},$
(17) $\kappa_{0}=F^{2}c(z^{2}_{+}\nu_{+}u_{0_{+}}+z^{2}_{-}\nu_{-}u_{0_{-}}).$
(18)

The values for bulk conductivities, i.e. $\sigma_{0}$ and $\kappa_{0}$, and
porosity are readily available from experimental data. The bulk diffusion
coefficient, transference number, and specific surface area are more difficult
to measure.

Comparison with experimental data indicates that the Bruggeman correlation
overestimates ionic conductivity in porous electrodes with certain particle
arrangements (Tjaden et al, 2016, 2018). Therefore, we also consider a
modified correlation that reduces ionic conductivity and diffusion coefficient
inside the porous material by equating
$\kappa=f_{m}\epsilon^{\frac{3}{2}}\kappa_{0}$ and
$D_{i}=f_{m}\epsilon^{\frac{3}{2}}D_{0,i}$, where $f_{m}$ is a correction
coefficient that takes values in the range $[0.02,0.79]$, depending on the
electrode microstructure (Madabattula and Kumar, 2020).

###  3.2 Nondimensionalization

To facilitate the analysis, we nondimensionalize the system described by (15)
by defining the dimensionless quantities:

$\displaystyle\hat{\Phi}_{1}=\frac{F}{RT}\Phi_{1},$ $\displaystyle\hat{c}$
$\displaystyle=\frac{c}{c_{0}},$ $\displaystyle\hat{t}$
$\displaystyle=\frac{D_{0}}{L^{2}}t,$ (19)
$\displaystyle\hat{\Phi}_{2}=\frac{F}{RT}(\Phi_{2}+U_{0}),$
$\displaystyle\hat{\mathbf{x}}$ $\displaystyle=\frac{\mathbf{x}}{L}.$

The dimensionless location $\hat{\mathbf{x}}$ leads to the nondimensional
domain $\hat{\Omega}$ and boundaries $\hat{\Gamma}_{\mathrm{a}}$ and
$\hat{\Gamma}_{\mathrm{c}}$, as well as gradient $\hat{\nabla}$. We also
introduce the dimensionless specific surface area $\hat{a}$, effective
diffusion coefficient $\hat{D}$, and effective electronic conductivity
$\hat{\sigma}$, which are design-dependent quantities we will describe in
Sect. 3.4. Using (10) and (11) the dimensionless field equations read

$\displaystyle-\hat{\nabla}\cdot\left(\hat{\sigma}\hat{\nabla}\hat{\Phi}_{1}\right)$
$\displaystyle=-\lambda\hat{a}\left(\delta_{r}\hat{i}_{n}+\delta_{c}\hat{i}_{c}\right)$
$\displaystyle\text{in }\hat{\Omega},$ (20a)
$\displaystyle-\hat{\nabla}\cdot\left(\hat{D}\hat{c}\hat{\nabla}\hat{\Phi}_{2}\right)-$
$\displaystyle\left(\frac{t_{+}}{z_{+}}+\frac{t_{-}}{z_{-}}\right)\hat{\nabla}\cdot\left(\hat{D}\hat{\nabla}\hat{c}\right)$
$\displaystyle=(1-\lambda)\hat{a}\left(\delta_{r}\hat{i}_{n}+\delta_{c}\hat{i}_{c}\right)$
$\displaystyle\text{in }\hat{\Omega},$ (20b)
$\displaystyle\frac{\partial(\epsilon\hat{c})}{\partial\hat{t}}-\hat{\nabla}\cdot(\hat{D}\hat{\nabla}\hat{c})$
$\displaystyle=\frac{z_{+}z_{-}}{t_{+}t_{-}(z_{-}-z_{+})}(1-\lambda)$
$\displaystyle\hat{a}t_{-}\bigg{(}\delta_{r}\hat{i}_{n}+$
$\displaystyle\left(\frac{dq_{+}}{dq}-\frac{t_{+}}{t_{-}}\frac{dq_{-}}{dq}\right)\delta_{c}\hat{i}_{c}\bigg{)}$
$\displaystyle\text{in }\hat{\Omega},$ (20c) subject to the boundary and
initial conditions $\displaystyle\hat{\Phi}_{1}$ $\displaystyle=0\quad$
$\displaystyle\text{on }\hat{\Gamma}_{\mathrm{a}},$ (20d)
$\displaystyle\hat{\Phi}_{1}$ $\displaystyle=\xi\hat{t}\quad$
$\displaystyle\text{on }\hat{\Gamma}_{\mathrm{c}},$ (20e)
$\displaystyle\hat{\nabla}\hat{\Phi}_{1}\cdot\mathbf{n}$
$\displaystyle=0\quad$ $\displaystyle\text{on
}\partial\hat{\Omega}\setminus(\hat{\Gamma}_{\mathrm{a}}\cup\hat{\Gamma}_{\mathrm{c}}),$
(20f) $\displaystyle\hat{\nabla}\hat{\Phi}_{2}\cdot\mathbf{n}$
$\displaystyle=0\quad$ $\displaystyle\text{on }\partial\hat{\Omega},$ (20g)
$\displaystyle\hat{\nabla}\hat{c}\cdot\mathbf{n}$ $\displaystyle=0\quad$
$\displaystyle\text{on }\partial\hat{\Omega},$ (20h)
$\displaystyle\hat{\Phi}_{1}(\hat{\mathbf{x}},\hat{t}=0)$
$\displaystyle=0\quad$ $\displaystyle\text{in }\hat{\Omega},$ (20i)
$\displaystyle\hat{\Phi}_{2}(\hat{\mathbf{x}},\hat{t}=0)$
$\displaystyle=0\quad$ $\displaystyle\text{in }\hat{\Omega},$ (20j)
$\displaystyle\hat{c}(\hat{\mathbf{x}},\hat{t}=0)$
$\displaystyle=\hat{c}_{0}\quad$ $\displaystyle\text{in }\hat{\Omega},$ (20k)

for $\hat{t}\in(0,\hat{T}]$. In the above we introduce the dimensionless
current densities

$\displaystyle\hat{i}_{n}=$
$\displaystyle\hat{c}^{\alpha}\left[e^{\alpha(\hat{\Phi}_{1}-\hat{\Phi}_{2})}-e^{-\alpha(\hat{\Phi}_{1}-\hat{\Phi}_{2})}\right],$
(21) $\displaystyle\hat{i}_{c}=$
$\displaystyle\hat{c}^{\alpha}\frac{\partial\left(\hat{\Phi}_{1}-\hat{\Phi}_{2}\right)}{\partial\hat{t}}.$

and dimensionless ionic conductivity $\hat{\kappa}=\hat{c}\hat{D}$.

For future reference we also introduce the independent dimensionless groups as

$\displaystyle\delta_{r}=\frac{a_{0}i_{0}L^{2}F}{RT}\left(\frac{1}{\sigma_{0}}+\frac{1}{\kappa_{0}}\right),$
$\displaystyle t_{+}$
$\displaystyle=\frac{z_{+}u_{+}}{z_{+}u_{+}-z_{-}u_{-}},$ (22)
$\displaystyle\delta_{c}=a_{0}C_{d}D_{0}\left(\frac{1}{\sigma_{0}}+\frac{1}{\kappa_{0}}\right),$
$\displaystyle\lambda$
$\displaystyle=\frac{\kappa_{0}}{\sigma_{0}+\kappa_{0}},$
$\displaystyle\xi=\frac{\zeta L^{2}F}{D_{0}RT}.$

with the following meanings

  * •

$\delta_{r}$ and $\delta_{c}$ represent the ratios between the kinetic
resistance and ohmic resistance, for redox kinetics and capacitive charging,
respectively;

  * •

$\lambda$ is the ratio between the ionic conductivity and the total
conductivity;

  * •

$t_{+}$ denotes the transference number of cation; the anion transference
number is related by $t_{-}=1-t_{+}$;

  * •

$\xi$ is the scan rate to diffusion timescale ratio.

Here, the ionic conductivity $\kappa_{0}$, the exchange current density
$i_{0}$, and specific capacitance $C_{d}$ are computed with the initial
concentration $c_{0}$.

We also define $\delta=\delta_{r}+\delta_{c}$ as the ratio of total kinetic
resistance and ohmic resistance, and $\gamma=\delta_{r}/\delta$ as the
fractional contribution from the redox reaction, so that

$\displaystyle\delta$
$\displaystyle=\frac{a_{0}D_{0}F}{RT}\left(\frac{1}{\sigma_{0}}+\frac{1}{\kappa_{0}}\right)\left(\frac{i_{0}L^{2}}{D_{0}}+\frac{C_{d}RT}{F}\right),$
(23) $\displaystyle\gamma$
$\displaystyle=\frac{i_{0}L^{2}/D_{0}}{i_{0}L^{2}/D_{0}+C_{d}RT/F}.$

Hence, $\gamma$ can be varied while keeping $\delta$ constant to have
different reaction compositions for the same total kinetic resistance. In the
literature, the Wagner number, i.e. $1/\delta$, is often used to analyze the
limiting kinetics and current distribution.

###  3.3 Electrode-specific electronic potential

In a full-cell model, the electronic potential only exists in the electrode
regions (where $\bar{\rho}\to 1$). System (20) models $\hat{\Phi}_{1}$ on the
entire domain $\Omega$, by assigning the dimensionless electronic conductivity
to a very small positive value in the separator region (where $\bar{\rho}\to
0$) so that $\hat{\sigma}\hat{\nabla}\hat{\Phi}_{1}\approx 0$ in the
electrolyte. However, because of filtering, $\bar{\rho}$ is not exactly zero
the electrolyte region and hence the interpolated conductivity may not be
negligible, which allows some electronic current to bypass the ionic phase,
creating a short circuit between the electrodes. To minimize the short circuit
effect, the optimized designs will position the electrodes far apart, so that
$\bar{\rho}$ approaches zero in the center of the electrolyte region.
Consequently, the optimization produces suboptimal designs, with no
possibility for creating interdigitated structures that are known to be
effective (Hung et al, 2022; Long et al, 2004).

To prevent a short circuit and allow for interdigitation, we represent
$\hat{\Phi}_{1}$ via two fields: $\hat{\Phi}_{1,\mathrm{a}}$ in the anode and
$\hat{\Phi}_{1,\mathrm{c}}$ in the cathode and use the indicator functions
from (5) to define their effective subdomains. The electronic potential
equation (15a) is thus replaced by the two equations

$\displaystyle-\hat{\nabla}\cdot\Big{(}\mathcal{I}_{\mathrm{a}}$
$\displaystyle\hat{\sigma}\hat{\nabla}\hat{\Phi}_{1,\mathrm{a}}\Big{)}$
$\displaystyle=-\lambda\hat{a}\left(\delta_{r}\hat{i}_{n,\mathrm{a}}+\delta_{c}\hat{i}_{c,\mathrm{a}}\right)\quad$
$\displaystyle\text{in }\hat{\Omega},$ (24a)
$\displaystyle-\hat{\nabla}\cdot\Big{(}\mathcal{I}_{\mathrm{c}}$
$\displaystyle\hat{\sigma}\hat{\nabla}\hat{\Phi}_{1,\mathrm{c}}\Big{)}$
$\displaystyle=-\lambda\hat{a}\left(\delta_{r}\hat{i}_{n,\mathrm{c}}+\delta_{c}\hat{i}_{c,\mathrm{c}}\right)\quad$
$\displaystyle\text{in }\hat{\Omega},$ (24b) with associated boundary
conditions $\displaystyle\hat{\Phi}_{1,\mathrm{a}}$ $\displaystyle=0\quad$
$\displaystyle\text{on }\hat{\Gamma}_{\mathrm{a}},$ (24c)
$\displaystyle\hat{\nabla}\hat{\Phi}_{1,\mathrm{a}}\cdot\mathbf{n}$
$\displaystyle=0\quad$ $\displaystyle\text{on
}\partial\hat{\Omega}\setminus\hat{\Gamma}_{\mathrm{a}},$ (24d)
$\displaystyle\hat{\Phi}_{1,\mathrm{c}}$ $\displaystyle=\xi\hat{t}\quad$
$\displaystyle\text{on }\hat{\Gamma}_{\mathrm{c}},$ (24e)
$\displaystyle\hat{\nabla}\hat{\Phi}_{1,\mathrm{c}}\cdot\mathbf{n}$
$\displaystyle=0\quad$ $\displaystyle\text{on
}\partial\hat{\Omega}\setminus\hat{\Gamma}_{\mathrm{c}},$ (24f)

where the anodic currents are computed as

$\displaystyle\hat{i}_{n,\mathrm{a}}=$
$\displaystyle\mathcal{I}_{\mathrm{a}}\hat{c}^{\alpha}\left[e^{\alpha(\hat{\Phi}_{1,\mathrm{a}}-\hat{\Phi}_{2})}-e^{-\alpha(\hat{\Phi}_{1,\mathrm{a}}-\hat{\Phi}_{2})}\right],$
(25) $\displaystyle\hat{i}_{c,\mathrm{a}}=$
$\displaystyle\mathcal{I}_{\mathrm{a}}\hat{c}^{\alpha}\frac{\partial\left(\hat{\Phi}_{1,\mathrm{a}}-\hat{\Phi}_{2}\right)}{\partial\hat{t}},$

and cathodic currents $\hat{i}_{n,\mathrm{c}}$ and $\hat{i}_{c,\mathrm{c}}$
are computed similarly by replacing $\mathcal{I}_{\mathrm{a}}$ and
$\hat{\Phi}_{1,\mathrm{a}}$ with $\mathcal{I}_{\mathrm{c}}$ and
$\hat{\Phi}_{1,\mathrm{c}}$. For all other terms defined above, the total
current densities are evaluated as
$\hat{i}_{n}=\hat{i}_{n,\mathrm{a}}+\hat{i}_{n,\mathrm{c}}$, and
$\hat{i}_{c}=\hat{i}_{c,\mathrm{a}}+\hat{i}_{c,\mathrm{c}}$, and the
electronic potential as
$\hat{\Phi}_{1}=\mathcal{I}_{\mathrm{a}}\hat{\Phi}_{1,\mathrm{a}}+\mathcal{I}_{\mathrm{c}}\hat{\Phi}_{1,\mathrm{c}}$.

###  3.4 Design-dependent material parameters

In our topology optimization, the model parameters such as the porosity, the
dimensionless surface area, electronic conductivity, and diffusivity, vary
spatially as a function of the material indicator field $\bar{\rho}$. The
porosity and dimensionless surface area are given by interpolations as

$\displaystyle\epsilon=$
$\displaystyle\epsilon_{M}+\bar{\rho}(\epsilon_{N}-\epsilon_{M}),$ (26)
$\displaystyle\hat{a}=$ $\displaystyle\bar{\rho}^{q},$

where $\epsilon_{M}$ is the porosity in the electrolyte/separator,
$\epsilon_{N}$ is the porosity in the electrode, and $q\geq 1$ is a
penalization parameter. The porosity $\epsilon$ is volume averaged while the
specific area $\hat{a}$ is penalized. As such, intermediate regions where
$\bar{\rho}(\hat{\mathbf{x}})\neq\\{0,1\\}$ have artificially smaller
$\hat{a}$, thus reducing current densities. Because of this, the optimized
designs predominantly have $\bar{\rho}(\hat{\mathbf{x}})=\\{0,1\\}$, which
approximates our original problem in (1), i.e. $\chi\approx\bar{\rho}$.
Meanwhile, the dimensionless electronic conductivity $\hat{\sigma}$ and
diffusion coefficients $\hat{D}$ are computed from the Bruggeman correlation
(16) and interpolations

$\displaystyle\hat{\sigma}=$
$\displaystyle\bar{\rho}^{p}(1-\epsilon_{N})^{\frac{3}{2}},$ (27)
$\displaystyle\hat{D}=$
$\displaystyle\epsilon_{M}^{\frac{3}{2}}+\bar{\rho}^{p}\left(\epsilon_{N}^{\frac{3}{2}}-\epsilon_{M}^{\frac{3}{2}}\right),$

where the parameter $p\geq 1$ is similarly used to penalize ionic and
electronic transport in intermediate material regions. Note that
$\epsilon_{M}$ is not involved in the interpolation for $\hat{\sigma}$ since
the electrolyte has negligible electronic conductivity. For the modified
Bruggeman correlation, the $\epsilon_{N}^{\frac{3}{2}}$ is replaced by
$f_{m}\epsilon_{N}^{\frac{3}{2}}$ in the formula for $\hat{D}$. In this paper,
we choose $\epsilon_{M}=1$, $\epsilon_{N}=0.5$, and $f_{m}=0.02$. The choice
of the penalization parameters $p$ and $q$ is discussed in Sect. 5.1.

###  3.5 Energy metrics

We list quantities used to construct the optimization cost and constraint
functions, namely energy input, kinetic energy, and ohmic loss,

$\displaystyle E_{\mathrm{in}}=$
$\displaystyle\frac{1}{\lambda}\int_{0}^{\hat{T}}\int_{\partial\hat{\Omega}}\hat{\sigma}\hat{\nabla}\hat{\Phi}_{1}\cdot\mathbf{n}\hat{\Phi}_{1}\mathop{}\\!\mathrm{d}\hat{s}\mathop{}\\!\mathrm{d}\hat{t},$
(28a) $\displaystyle E_{\mathrm{kin}}=$
$\displaystyle\delta_{r}\int_{0}^{\hat{T}}\int_{\hat{\Omega}}\hat{a}\hat{i}_{n}(\hat{\Phi}_{1}-\hat{\Phi}_{2})\mathop{}\\!\mathrm{d}\hat{x}\mathop{}\\!\mathrm{d}\hat{t}$
$\displaystyle+$
$\displaystyle\delta_{c}\int_{0}^{\hat{T}}\int_{\hat{\Omega}}\hat{a}\hat{i}_{c}(\hat{\Phi}_{1}-\hat{\Phi}_{2})\mathop{}\\!\mathrm{d}\hat{x}\mathop{}\\!\mathrm{d}\hat{t},$
(28b) $\displaystyle E_{\mathrm{ohm}}=$
$\displaystyle\frac{1}{\lambda}\int_{0}^{\hat{T}}\int_{\hat{\Omega}}\hat{\sigma}\hat{\nabla}\hat{\Phi}_{1}\cdot\hat{\nabla}\hat{\Phi}_{1}\mathop{}\\!\mathrm{d}\hat{x}\mathop{}\\!\mathrm{d}\hat{t}$
$\displaystyle+$
$\displaystyle\frac{1}{1-\lambda}\int_{0}^{\hat{T}}\int_{\hat{\Omega}}\hat{D}\hat{c}\hat{\nabla}\hat{\Phi}_{2}\cdot\hat{\nabla}\hat{\Phi}_{2}\mathop{}\\!\mathrm{d}\hat{x}\mathop{}\\!\mathrm{d}\hat{t}$
$\displaystyle+$
$\displaystyle\frac{1}{1-\lambda}\left(\frac{t_{+}}{z_{+}}+\frac{t_{-}}{z_{-}}\right)\int_{0}^{\hat{T}}\int_{\hat{\Omega}}\hat{D}\hat{\nabla}\hat{c}\cdot\hat{\nabla}\hat{\Phi}_{2}\mathop{}\\!\mathrm{d}\hat{x}\mathop{}\\!\mathrm{d}\hat{t}.$
(28c)

The detailed derivation of these expressions is presented in Appendix A. The
kinetic energy, henceforth referred to as the stored energy, includes both the
energy stored in the chemical reaction and the dissipation due to finite
charge-transfer kinetics. The energy storage depends on the thermodynamic
potential, which is characteristic of a specific material system. To keep the
analysis general, we do not specify the thermodynamic potential and instead
infer its value via the dimensionless ionic potential, which naturally groups
kinetic loss with kinetic storage.

The optimization problem focuses on reducing the ohmic losses, because this
can be lowered by shaping the electrodes to facilitate ion and electron
transport across the cell. By reducing the ohmic losses, a larger current
output can be extracted from the cell, for a given applied voltage. While a
larger current also causes larger kinetic losses (due to the activation
barrier for redox reaction) we expect the cell efficiency to be dominated by
transport losses.

##  4 Implementation

We resolve the governing equations (20), (24), filtering equation (2), and
boundary propagation equation (4) with the finite element library Firedrake
(Ham et al, 2023), which uses PETSc (Balay et al, 2023a, b) as the backend for
the iterative solvers. The sensitivities of the optimization cost and
constraint functions are computed by pyadjoint (Mitusch et al, 2019). A Python
implementation of the MMA algorithm (Svanberg, 1987), pyMMAopt (Salazar de
Troya, 2021) is used to solve the optimization problem. The designs are
considered converged after 350 iterations, where the change in the normalized
objective function is typically in the magnitude of $10^{-3}$.

Piecewise linear finite elements are used for spatial discretization of the
governing physics equations (20) and (24). We use a mixed finite element
method to solve the filtering equation (2) and the continuous boundary
propagation problem (4), as detailed in Appendix B. The domain is either
partitioned into triangular (2D) or tetrahedral (3D) elements. The two-
dimensional domain is a rectangle with distance between the current collectors
of nondimensional length 1 and the width of 2. The additional dimension of the
three-dimensional domain is a depth of 2. We simulate half of the total domain
in 2D, and a quarter in 3D, due to the domain symmetry. The total number of
degrees of freedom (DoF) for 2D and 3D models are $\sim$60k and $\sim$4M,
respectively. Backward Euler is used for temporal discretization with 20 total
timesteps for a dimensionless simulation time of $1$, which is the smallest
number of timesteps that is stable for all studies. We parallelize the
computation via an approximately 20k DoFs/CPU load.

With a transient problem, computing sensitivities with the adjoint method
requires storing the primal response for all timesteps. In this work, we have
not encountered any memory bottleneck, but it is an inevitability for larger
simulations. In response, checkpointing schemes such as Griewank and Walther
(2000) can be utilized.

The finite element system corresponding to (20) with the two-electronic
potential formulation (24) is solved using Newton’s method with linesearch.
The 2D linearized systems are solved using a direct method. For the 3D cases,
however, an iterative method with a scalable preconditioner is essential
(Wathen, 2015). The linearized systems have the following block form:

$\begin{bmatrix}A_{11}&A_{12}&A_{13}&A_{14}\\\ A_{21}&A_{22}&A_{23}&A_{24}\\\
A_{31}&A_{32}&A_{33}&A_{34}\\\ A_{41}&A_{42}&A_{43}&A_{44}\\\
\end{bmatrix}\begin{bmatrix}\Delta\mathbf{\Phi}_{1,\mathrm{c}}\\\
\Delta\mathbf{\Phi}_{1,\mathrm{a}}\\\ \Delta\mathbf{\Phi}_{2}\\\
\Delta\mathbf{c}\end{bmatrix}=\begin{bmatrix}\mathbf{r}_{1}\\\
\mathbf{r}_{2}\\\ \mathbf{r}_{3}\\\ \mathbf{r}_{4}\\\ \end{bmatrix},$ (29)

where $\mathbf{r}_{i},\ i\in\\{1,2,3,4\\}$, are the residual vectors
corresponding to Eqn. (24b), (24a), (20b), and (20c) that are expressed in
terms of the $j\in\\{1,2,3,4\\}$ response vectors $\Delta\mathbf{\Phi}_{1,c}$,
$\Delta\mathbf{\Phi}_{1,a}$, $\Delta\mathbf{\Phi}_{2}$, and $\mathbf{c}$. Each
block $A_{ij}$ represents the linearization of the finite element residual
$\mathbf{r}_{i}$ with respect to vector $j$.

The diagonal blocks $A_{ii}$ result from the discretization of elliptic or
parabolic operators, for which multigrid methods (Brandt, 1977) are often
ideal preconditioners. We thus use the preconditioned GMRES (Saad and Schultz,
1986), with a block triangular preconditioner of the form

$P=\begin{bmatrix}A_{11}&0&0&0\\\ A_{21}&A_{22}&0&0\\\
A_{31}&A_{32}&A_{33}&0\\\ A_{41}&A_{42}&A_{43}&A_{44}\\\ \end{bmatrix}.$ (30)

Evaluating the inverse of this matrix requires the inverses of the diagonal
blocks $A_{ii}$, which is done approximately using a single AMG V-cycle (Ruge
and Stüben, 1987) BoomerAMG (Henson and Yang, 2002) from the hypre library
(Falgout and Yang, 2002). The PETSc options for solving the field equations,
the filter problem, and boundary propagation problem can be found in Appendix
C.

Figure 4: Strong scaling for the 3D solver. The timings show the total
solution time and the total linear solver time where the number of GMRES
iterations appears in the parentheses.

In Fig. 4, we investigate the scalability of this preconditioning approach.
For this investigation, we simulate the first timestep for the initial uniform
$\rho=0.5$ design. We observe that the linear solver of (29) exhibits almost
linear scalability, which indicates that our preconditioning is effective.

##  5 Results and discussion

Herein, we introduce the optimization strategy and optimize designs for
various choices of dimensionless parameters. The fixed dimensionless
parameters among all designs are listed in Table 1.

Table 1: Dimensionless parameters that are fixed amongst all designs. Parameter | Quantity | Value
---|---|---
$\alpha$ | charge transfer coefficient | 0.5
$t_{+}$ | cation transference number | 0.5
$\xi$ | potential scan rate | 1.0
$\hat{T}$ | final time | 1.0
$\epsilon_{M}$ | electrolyte porosity | 1.0
$\epsilon_{N}$ | electrode porosity | 0.5

Table 2 lists the three dimensionless parameters that are varied amongst the
designs. The results presented in Sect. 5.2 show how different combinations of
$(\lambda,\delta,\gamma)$ lead to different optimal electrode designs.

Table 2: Dimensionless parameters that are varied amongst the designs. Parameter | Quantity
---|---
$\lambda$ | ratio between ionic and total conductivity
$\delta$ | ratio between ohmic and kinetic resistance
$\gamma$ | redox contribution to total reaction

###  5.1 Optimization problem setup

In this section, we introduce several strategies devised to improve the
convergence and effectiveness of the optimization algorithm. Given the
complexity of the physical problem, we need to carefully define the cost and
constraint functions and the initial design in order to effectively traverse
the design space and find optimized designs.

####  5.1.1 Parameter penalizations

To promote binary $\bar{\rho}(\hat{\mathbf{x}})\approx\\{0,1\\}$ designs, we
use different interpolations for the design-dependent quantities $\hat{a}$,
$\hat{\sigma}$, and $\hat{D}$ in the cost and constraint functions, as
compared to the primal model, cf. (20) and (24) in Sect. 3.4.

Figure 5: Different penalization for specific surface area and ionic
conductivity for the primal analysis and cost/constraint computation to
produce binary designs.

In the cost and constraint computations, $\hat{a}$ is interpolated with $q=3$,
whereas in the primal analysis $q=1$, cf. (26). Reversely, $\hat{\sigma}$ and
$\hat{D}$ use $p=1$ in the cost and constraint computations and $p=1.5$ in the
primal analysis. The different interpolations adopted for the primal analysis
and the cost/constraint computations are visualized in Fig. 5. The $q=3$
penalization for the specific surface area and $p=1.5$ for the conductivity
disproportionately reduces the energy storage and disproportionately increases
the ohmic loss in the intermediate material regions. These penalizations lead
to inefficient use of intermediate material regions and hence the size of such
regions in the optimized designs is small.

####  5.1.2 Cost and Constraint

Maximizing the kinetic energy $E_{\mathrm{kin}}$ defined in (28) is an
intuitive choice for the (negative) cost function, as a proxy for maximizing
energy stored. Unfortunately, this choice does not consider the ohmic loss.
Most importantly, this cost function would not benefit from the (27)
diffusivity/conductivity interpolation scheme. In response, we formulate the
optimization problem as follows

$\max_{\rho\colon\Omega\to[0,1]}\theta_{0}=\frac{1}{2}\left(E_{\mathrm{kin}}+E_{\mathrm{in}}-E_{\mathrm{ohm}}\right)\\\
\text{subject to
}\theta_{1}=\frac{E_{\mathrm{ohm}}}{E_{\mathrm{in}}}\leq\eta_{\max},$ (31)

i.e. we maximize energy storage and ensure a minimum energy efficiency of
$(1-\eta_{\max})$. According to (31), the design can be optimized by either
increasing $E_{\mathrm{kin}}$ or decreasing $E_{\mathrm{ohm}}$. And hence the
interpolations of $\hat{a}$, $\hat{\sigma}$, and $\hat{D}$ produce
disproportionately lower $E_{\mathrm{kin}}$ and disproportionately higher
$E_{\mathrm{ohm}}$ from the intermediate material regions for the cost and
constraint functions, relative to the primal analysis. As such the optimizer
generates predominantly $0-1$ designs.

Figure 6: Interdigitated electrodes with near uniform $\beta\approx\pm 1$
distribution bounded by a thin separator layer with large $(1-|\beta|)^{3}$
values.

According to (24), the only relation between the two electronic potentials is
through the reaction terms. Consequently, a physical connection between the
anode and cathode is not simulated as a short circuit (electronic current
bypassing the reaction). To prevent such connections, we define the short
circuit intensity

$\mathcal{I}_{\mathrm{SC}}=\int_{\Omega}\left(1-|\beta|\right)^{3}\bar{\rho}\mathop{}\\!\mathrm{d}\hat{x}$
(32)

The term $(1-|\beta|)^{3}$ takes its maximum value when $\beta=0$. In the
absence of a short circuit, in the electrodes where $\bar{\rho}\approx 1$, we
have a near uniform $\beta=\pm 1$ distribution. The exception occurs in the
separator layer where $\beta$ experiences a sharp transition as illustrated in
Fig. 6. A “small” short circuit intensity requires $\bar{\rho}\approx 0$ when
$\beta\neq\pm 1$, thus ensuring the electrolyte fills the gap between the
electrodes. We weakly enforce this short circuit intensity constraint via the
penalty method, i.e. we now solve the optimization problem

$\min_{\rho\colon\Omega\to[0,1]}\tilde{\theta}_{0}=\frac{1}{\theta_{0}}+w_{\mathrm{SC}}\mathcal{I}_{\mathrm{SC}}\\\
\text{subject to
}\theta_{1}=\frac{E_{\mathrm{ohm}}}{E_{\mathrm{in}}}\leq\eta_{\max},$ (33)

where $w_{\mathrm{SC}}=1$ is the penalty weight used in our study. Note that
the minimization of $-\theta_{0}$, is less robust than $\frac{1}{\theta_{0}}$
since it requires case-specific scaling values for $\theta_{0}$. However, in
order for (33) to work properly,
$E_{\mathrm{ohm}}<E_{\mathrm{in}}+E_{\mathrm{kin}}$ must be satisfied at
initial design, so that the cost $\tilde{\theta}_{0}$ is always positive, and
hence minimized by maximizing $\theta_{0}$. If $\theta_{0}<0$, then
$\tilde{\theta}_{0}$ is minimized as $\theta_{0}\rightarrow 0^{-}$, which
completely changes the optimization trajectory and hence designs.

It is well known that the optimized designs may depend on the initial design.
To this end, we solve the optimization problem using three uniform initial
designs $\rho(\mathbf{x})\in\\{0.45,0.5,0.55\\}$ and choose the best overall
design. Let $\theta_{1}^{(0)}$ denote the ohmic loss ratio of the initial
design. We set the constraint bound $\eta_{\max}=\theta_{1}^{(0)}*\Sigma$,
where $\Sigma\in(0.0,1.0]$. We generate designs for different combinations of
dimensionless parameters, with both original and modified Bruggeman
correlation. $\Sigma=0.4$ is used for all $\delta=0.5$, $\lambda=0.1$, and the
three-dimensional designs; $\Sigma=0.5$ is used for all other designs.

Figure 7: Cost and constraint function history for $\delta=2,\gamma=0.5$,
$\lambda=0.01$ with modified Bruggeman correlation. The inserts are the
designs at iteration (a) 8, (b) 20, (c) 50, and (d) 150. Only one symmetric
half of the design is illustrated to save space.

A typical optimization convergence history is illustrated in Fig. 7. During
early design iterations, the constraint drives the optimization to satisfy the
$\theta_{1}\leq\eta_{\mathrm{max}}$ constraint by separating the electrodes.
Also, more electrode material leads to additional reaction current and
increased ohmic loss. Because of this, the constraint limits the maximum
amount of electrode material and thus prevents short circuits from appearing
in the design as the optimization continues. While eliminating short circuits
is necessary, a tight constraint hinders the development of structures as the
optimization continues. Therefore, we deactivate this constraint after the
bulk design features have formed, which for our study is after iteration 150.
The ensuing iterations refine the shape rather than the overall layout. To
summarize, $\Sigma$ strongly affects the optimization during early iterations
to find the interdigitation topology but it has no effect during later
iterations.

####  5.1.3 Case specific treatments

Certain cases within our parameter sweep require extra care. 1) For cases with
the original Bruggeman correlation, i.e. $f_{m}=1$, large effective ionic
conductivity/diffusivity within the porous electrode allows ions to travel
deeper into the pore network, resulting in better material utilization as
discussed in Sect. 5.3. Accordingly, during early design iterations, the
optimizer positions a considerable amount of electrode material in the middle
region to benefit from the largest potential drop. Unfortunately this causes a
strong connection between the electrodes that cannot be removed in the later
design iterations. To remedy this, the optimization adopts a continuation
strategy wherein the initial design uses the modified Bruggemann correlation
with $f_{m}=0.02$ which increases to the target $f_{m}=1$ after iteration 50.
2) For cases with large $\delta$, we observe $\theta_{0}<0$ from the initial
design, which is undesirable as stated above. We remedy this by invoking yet
another continuation strategy where the initial design uses a smaller $\delta$
such that $E_{\mathrm{ohm}}<E_{\mathrm{in}}+E_{\mathrm{kin}}$ is satisfied. In
our study, the $\delta=5$ cases start with $\delta=2$ and increase to its
target value after iteration 50 (see Fig. 7).

In some extreme situations, namely the $\delta=5$ with modified Bruggeman
correlation cases, large regions of intermediate material appear in the
optimized designs. The small values for $E_{\mathrm{kin}}$ and
$E_{\mathrm{ohm}}$ in these regions are overshadowed by the other regions
where the current density is substantial. As such, the penalization scheme is
not effective. To counter this effect, we scale $E_{\mathrm{kin}}$ and
$E_{\mathrm{ohm}}$ after the bulk structures are established (i.e. after
iteration 150 when the efficiency constraint is deactivated). To motivate our
scaling scheme, refer to the typical $\beta$ distribution shown in Fig. 3
(right-hand side) where $|\beta|\approx 0$ around the electrode tips and
$|\beta|\approx 1$ near the current collectors. Accordingly, we consider three
different types of scaling, namely $|\beta|$, $1-|\beta|$, and
$|\beta|(1-|\beta|)$, to scale down the energy integrands around the electrode
tip, the current collector, and both, respectively. The scaling is readily
incorporated into our strategy with a continuation scheme. With $|\beta|$
scaling as an example, the energy storage is computed by

$\tilde{E}_{\mathrm{kin}}=\delta_{r}\int_{0}^{\hat{T}}\int_{\hat{\Omega}}|\beta|^{s}\hat{a}\hat{i}_{n}(\hat{\Phi}_{1}-\hat{\Phi}_{2})\mathop{}\\!\mathrm{d}\hat{x}\mathop{}\\!\mathrm{d}\hat{t}\\\
+\frac{\delta_{c}}{2}\int_{\hat{\Omega}}|\beta|^{s}\hat{a}(\hat{\Phi}_{1}-\hat{\Phi}_{2})^{2}\mathop{}\\!\mathrm{d}\hat{x},$
(34)

where $s=0$ for the first 150 iterations where after it is increased to $s=1$.
We use the $|\beta|$ scaling for $E_{\mathrm{kin}}$ and $E_{\mathrm{ohm}}$ in
the $\gamma=\\{0.5,1\\}$ cases. For the $\gamma=0$ case, we use $|\beta|$
scaling for $E_{\mathrm{kin}}$ and $|\beta|(1-|\beta|)$ scaling for
$E_{\mathrm{ohm}}$.

For the 3D optimization case, a slight connection of the electrodes appears
near a wall. Some of the strategies described above could be used to prevent
this short circuit design: different initial designs, scaling, etc. Instead,
motivated by Lazarov et al (2016), we implement a simple dilation/erosion
strategy, which allows us to remove the short circuit. This strategy applies a
one-time adjustment to the design field $\rho$, which is localized to the
connection region. Namely, we use shifted versions of the electrode indicator
functions (5) such that

$\displaystyle\tilde{\mathcal{I}}_{\mathrm{a}}$
$\displaystyle=H_{100,0.52}\left(\frac{\beta+1}{2}\right),$ (35)
$\displaystyle\tilde{\mathcal{I}}_{\mathrm{c}}$
$\displaystyle=H_{100,0.52}\left(\frac{-\beta+1}{2}\right).$

In this way, regions where $\tilde{\mathcal{I}}_{\mathrm{a}}=1$ and
$\tilde{\mathcal{I}}_{\mathrm{c}}=1$ are reduced in size as compared to
$\mathcal{I}_{\mathrm{a}}$ and $\mathcal{I}_{\mathrm{c}}$ computed by
$H_{100,0.5}$. Then there exist a subdomain $\Omega_{\mathrm{s}}$ such that
$\tilde{\mathcal{I}}_{\mathrm{a}}(\mathbf{x})+\tilde{\mathcal{I}}_{\mathrm{c}}(\mathbf{x})=0$
for all $\mathbf{x}\in\Omega_{\mathrm{s}}$. Subsequently we equate
$\rho\coloneqq\rho(\tilde{\mathcal{I}}_{\mathrm{a}}+\tilde{\mathcal{I}}_{\mathrm{c}})$
to ensure $\rho=0$ in $\Omega_{\mathrm{s}}$. Once the connection is removed,
it is unlikely to reconnect due to the short circuit intensity penalization.
We perform the shifted indicator adjustment after iteration 150. Note that the
remaining iterations are still completed with the original indicator functions
$\mathcal{I}_{\mathrm{a}}$ and $\mathcal{I}_{\mathrm{c}}$. This disconnection
strategy is only appropriate when the electrode connections are of a small
“grazing” type, and if so, has very little effect on the trajectory of the
optimization. Meanwhile, “head-on collision” connections should be avoided
earlier in the optimization using strategies mentioned above.

###  5.2 Discussion of generated designs

We generate optimized designs with various combinations of the dimensionless
parameters $\delta$, $\lambda$, and $\gamma$. We choose energy storage and
ohmic loss as performance metrics and compare the optimized designs to a
monolithic design with electrode separation distance of 0.05, as this is the
average gap among all optimized designs. The optimized designs are shown in
Fig. 8 \- 12 with the percentage improvement over the monolithic design. In
the following subsections, we discuss the designs corresponding to various
combinations of the three nondimensional parameters listed in Table 2, and
Bruggeman correlation functions (which approximates the tortuosity of the ion-
diffusion pathway within the porous electrodes.)

Here, $\lambda$, i.e. the ratio between ionic and total conductivity (defined
in (22)), affects the macroscale localization of the reaction front. In cases
where ionic conductivity is limiting ($\lambda<0.5$), the reaction front
localizes at the electrode/separator interface. Conversely, when the
electronic conductivity is the limiting factor, ($\lambda>0.5$), the reaction
front shifts towards the current collector.

$\delta$ 0.5 5 $\gamma$ 0.0 $E_{\mathrm{kin}}=.0427\,(110\%),\
E_{\mathrm{ohm}}=.0038\,(39\%)$ $E_{\mathrm{kin}}=.1098\,(231\%),\
E_{\mathrm{ohm}}=.0403\,(65\%)$ (a) (b) 0.5 $E_{\mathrm{kin}}=.0383\,(107\%),\
E_{\mathrm{ohm}}=.0024\,(41\%)$ $E_{\mathrm{kin}}=.1395\,(207\%),\
E_{\mathrm{ohm}}=.0404\,(56\%)$ (c) (d) 1.0 $E_{\mathrm{kin}}=.0334\,(103\%),\
E_{\mathrm{ohm}}=.0015\,(42\%)$ $E_{\mathrm{kin}}=.2254\,(192\%),\
E_{\mathrm{ohm}}=.0497\,(59\%)$ (e) (f)

Figure 8: Optimized porous electrode designs modeled with the original
Bruggeman correlation, $\lambda=0.01$, and varying $\delta$/$\gamma$. Black is
$\bar{\rho}=1$; white is $\bar{\rho}=0$.

The parameter $\delta$, i.e. the ratio between ohmic and kinetic resistance
(defined in (22)), controls the electrode penetration depth. In the absence of
concentration gradients, the relative reaction rate between the back and the
front of electrodes can be quantified by (Fuller and Harb, 2018)

$\text{Relative rate}\approx\frac{1}{\cosh\left(\sqrt{\delta}\right)}.$ (36)

For small $\delta$ values, the system is dominated by the kinetic resistance,
and results in a uniform or near-uniform distribution of current density,
which leads to a uniform reaction throughout the electrode. In this case, the
macroscale geometry of the electrode has less effect on the electrode material
utilization, i.e. monolithic designs perform sufficiently well. On the other
hand, for large $\delta$ values, the system is dominated by ohmic resistance
and results in a non-uniform (localized) distribution of current density. As
mentioned in the introduction, thick monolothic electrodes perform poorly due
to small penetration depth. Therefore, this case favors designs with complex
electrode geometries.

####  5.2.1 Original Bruggeman correlation

Monolithic electrodes

(a) $\hat{c}\in[0.7,1.3],\ \hat{c}_{\mathrm{var}}=0.2161$

(b) $\hat{\Phi}_{2}\in[0.15,0.85],\
|\hat{\nabla}\hat{\Phi}_{2}|_{\mathrm{avg}}=0.8280$

Optimized electrodes

(c) $\hat{c}\in[0.8,1.2],\ \hat{c}_{\mathrm{var}}=0.1015$

(d) $\hat{\Phi}_{2}\in[0.3,0.7],\
|\hat{\nabla}\hat{\Phi}_{2}|_{\mathrm{avg}}=0.4932$

Figure 9: Concentration (top) and ionic potential (bottom) distributions of
monolithic and optimized electrodes with $\delta=5$, $\gamma=1$,
$\lambda=0.01$, and the original Bruggeman correlation.

$\lambda$ 0.01 0.1 $\gamma$ 0.0 $E_{\mathrm{kin}}=.0910\,(179\%),\
E_{\mathrm{ohm}}=.0189\,(40\%)$ $E_{\mathrm{kin}}=.0731\,(136\%),\
E_{\mathrm{ohm}}=.0337\,(66\%)$ (a) (b) 0.5 $E_{\mathrm{kin}}=.1040\,(163\%),\
E_{\mathrm{ohm}}=.0153\,(39\%)$ $E_{\mathrm{kin}}=.0775\,(118\%),\
E_{\mathrm{ohm}}=.0283\,(68\%)$ (c) (d) 1.0 $E_{\mathrm{kin}}=.1092\,(135\%),\
E_{\mathrm{ohm}}=.0145\,(47\%)$ $E_{\mathrm{kin}}=.0846\,(106\%),\
E_{\mathrm{ohm}}=.0250\,(78\%)$ (e) (f)

Figure 10: Optimized porous electrode designs modeled with the original
Bruggeman correlation, $\delta=2$, and varying $\lambda$/$\gamma$. Black is
$\bar{\rho}=1$; white is $\bar{\rho}=0$.

Shown in Fig. 8 is the first series of designs generated using the original
Bruggeman correlation, a fixed $\lambda=0.01$, and varying $\delta$ and
$\gamma$. All optimized designs possess an interdigitated structure, which was
found by the optimization algorithm thanks to the continuation schemes
described in Sect. 5.1. The interdigitation allows reactions to happen closer
to the current collectors, boosting current density. The structure also
reduces ohmic losses by promoting more efficient ionic transport. It is
important to note that, if we simply force electrode separation by increasing
the weight $w_{SC}$ on the short circuit intensity penalty in (33), we obtain
near monolithic designs with a large distance between electrodes.

The first benefit of interdigitation is a more uniform concentration
distribution. We quantify the concentration variation by

$\hat{c}_{\mathrm{var}}=|\hat{\Omega}|^{-\frac{1}{2}}\|\hat{c}-\bar{c}\|_{L^{2}(\hat{\Omega})},$
(37)

where $\bar{c}=1$ is the average dimensionless concentration in the system. As
illustrated in Fig. 9, using the $\delta=5$ and $\gamma=1$ case as an example,
the interdigitated design decreases the variation in concentration by more
than 50% over the monolithic design. A more uniform concentration distribution
ensures better ion adsorption and redox reaction, since small concentrations
hinder the adsorption/reaction rate, while large localized concentrations
indicate poor ion utilization.

In the case of redox-active materials ($\gamma=1$), the interdigitation also
increases the electrode area in contact with the pure electrolyte region,
facilitating ion transport between electrolyte and electrode regions. For an
applied voltage boundary condition, interdigitation reduces the cell
resistance and hence increases the current output. Summarizing, by
facilitating ion transport in the electrolyte, the interdigitated design
reduces the energy dissipation and improves cell efficiency.

For small $\lambda$, energy dissipation mainly results from ion transport and
is proportional to the magnitude of the ionic-potential gradient. We compute
the average gradient magnitude of ionic potential by

$|\hat{\nabla}\hat{\Phi}_{2}|_{\mathrm{avg}}=|\hat{\Omega}|^{-1}\|\hat{\nabla}\hat{\Phi}_{2}\|_{L^{2}(\hat{\Omega})}.$
(38)

Fig. 9 shows a more uniform ionic potential distribution in the optimized
design compared to that in the monolithic design, and a $40\%$ decrease in the
average gradient.

In kinetic resistance limited cases ($\delta=0.5$), there are no notable
differences between the designs obtained with different $\gamma$ values (see
left column in Fig. 8.). In other words, the bulk geometry is not sensitive to
the partitioning between capacitive and redox reaction due to reactions being
more uniform. Indeed, the designs are driven by the heterogeneous ohmic losses
$E_{\mathrm{ohm}}$. In these cases, the interdigitated optimized designs
achieve a maximum energy storage increase of 10%, and a maximum ohmic loss
decrease of 60% over the monolithic design.

For the ohmic resistance limited cases ($\delta=5$), we notice obvious
variations among the designs (see the right column in Fig. 8). The periodicity
of the designs remains the same for varying $\gamma$ values. Since the current
model does not incorporate surface concentration that limits the maximum ion
adsorption, the overall concentration decreases significantly for capacitor
electrodes. Therefore, for capacitive-dominated reactions (small $\gamma$),
more electrolyte is needed to support the reaction and prevent reactant
starvation and dramatic energy dissipation. This explains the wider white
separation region between the electrodes in the $\gamma=0$ design. As the
charging becomes more redox dominated ($\gamma\rightarrow 1$), the designs’
interdigitation depth increases, reducing the separation distance between the
electrodes. Such an effect is only obvious in these ohmic resistance limited
cases. A smaller separation distance results in a faster ion exchange between
the electrodes and hence sustains high reaction rates. Additionally, more
electrode material is placed close to the current collectors to support the
reaction on the electrode tips. For these $\delta=5$ designs, we observe
approximately 100% increase in energy storage and 40% reduction in ohmic loss
over the monolithic design.

Next, we study the effect of varying $\lambda$ and $\gamma$ for a fixed
$\delta=2$. A small $\lambda$ indicates larger relative electronic
conductivity, which benefits designs with more intensive interdigitation. A
small $\lambda$ is characteristic of Graphite-based anodes with large
electronic conductivity, in the order of 1000mScm, and in porous cathodes that
incorporate carbon additives to enhance the average electric conductivity. On
the other hand, the average room-temperature conductivity of liquid
electrolyte materials is on the order of 10mScm. Overall, it is reasonable to
assume that the electronic conductivity is at least one or two orders of
magnitude larger than the ionic conductivity. This corresponds to $\lambda$
values ranging from 0.01 to 0.1.

As shown in Fig. 10, we observe long and bulky columns in the $\lambda=0.01$
designs. These designs achieve substantial current densities over the
electrode tip due to the larger over-potential. The ohmic loss associated with
transporting the electrons to and from the current collector is negligible,
and thus results in more efficient designs.

$\delta$ 0.5 5 $\gamma$ 0.0 $E_{\mathrm{kin}}=.0328\,(656\%),\
E_{\mathrm{ohm}}=.0101\,(204\%)$ $E_{\mathrm{kin}}=.0962\,(786\%),\
E_{\mathrm{ohm}}=.0465\,(312\%)$ (a) (b) 0.5
$E_{\mathrm{kin}}=.0291\,(557\%),\ E_{\mathrm{ohm}}=.0069\,(120\%)$
$E_{\mathrm{kin}}=.1152\,(855\%),\ E_{\mathrm{ohm}}=.0400\,(267\%)$ (c) (d)
1.0 $E_{\mathrm{kin}}=.0248\,(452\%),\ E_{\mathrm{ohm}}=.0048\,(72\%)$
$E_{\mathrm{kin}}=.1330\,(851\%),\ E_{\mathrm{ohm}}=.0414\,(265\%)$ (e) (f)

Figure 11: Optimized porous electrode designs modeled with the modified
Bruggeman correlation, $\lambda=0.01$, and varying $\delta$/$\gamma$. Black is
$\bar{\rho}=1$; white is $\bar{\rho}=0$.

Comparing the $\lambda=0.1$ cases (right column in Fig. 10) with the
$\lambda=0.01$ cases (left column in Fig. 10), we notice that an increase in
$\lambda$ favors designs with more complexity. This can be attributed to the
fact that as $\lambda$ approaches $0.5$, electron and ion transport become
equally limiting. Consequently, smaller electrode columns are needed to
shorten electron transport distance. This explain why the designs have more
thinner columns. And unlike the sharp teeth-like structures observed in Roy et
al (2022), the features here have round corners since smaller curvature
enhances ion diffusion, and prevents electric field congestion. Deep
interdigitation provides smaller benefits as compared to the previous cases.
However, maximizing the interfacial area with the pure electrolyte region is
still desirable as it introduces additional pathways for ionic current to
enter and leave the pure electrolyte region. Therefore, the complex bulk
geometry enhances the ionic transport within the porous electrode region, and
thus reduces local ion depletion.

The optimization mainly tries to reduce energy dissipated due to ohmic
resistance in these $\delta=2$ designs. Therefore, the $\lambda=0.1$ designs
with increased electronic resistance generally store less energy. Although all
optimized designs outperform the monolithic design, the designs generated with
$\lambda=0.01$ provide an additional $40\%$ increase in energy storage and
$30\%$ reduction in ohmic loss, as compared to the corresponding $\lambda=0.1$
cases.

####  5.2.2 Modified Bruggeman correlation

$\lambda$ 0.01 0.1 $\gamma$ 0.0 $E_{\mathrm{kin}}=.0803\,(858\%),\
E_{\mathrm{ohm}}=.0189\,(165\%)$ $E_{\mathrm{kin}}=.0601\,(617\%),\
E_{\mathrm{ohm}}=.0167\,(136\%)$ (a) (b) 0.5
$E_{\mathrm{kin}}=.0761\,(798\%),\ E_{\mathrm{ohm}}=.0169\,(159\%)$
$E_{\mathrm{kin}}=.0592\,(602\%),\ E_{\mathrm{ohm}}=.0195\,(173\%)$ (c) (d)
1.0 $E_{\mathrm{kin}}=.0734\,(741\%),\ E_{\mathrm{ohm}}=.0153\,(155\%)$
$E_{\mathrm{kin}}=.0572\,(594\%),\ E_{\mathrm{ohm}}=.0131\,(160\%)$ (e) (f)

Figure 12: Optimized porous electrode designs modeled with the modified
Bruggeman correlation, $\delta=2$, and varying $\lambda$/$\gamma$. Black is
$\bar{\rho}=1$; white is $\bar{\rho}=0$.

We now compare the results generated using the modified Bruggeman correlation
with $f_{m}=0.02$, which is equivalent to reducing ionic diffusion inside the
porous electrode by 50 times. As a result, the penetration depth of the
reaction front in the monolithic design is further diminished. We first
compare the designs from different reaction regimes. The results are shown in
Fig. 11. To compensate for extremely slow ionic diffusion inside the electrode
phase, all designs form spikey teeth-like structures to increase the
interfacial area with the pure electrolyte.

Similar to the observations made from the original Bruggeman cases, among the
$\delta=0.5$ designs, changing $\gamma$ does not noticeably change the design
because it does not improve the uniformity of current distribution. For the
$\delta=5$ designs, capacitor electrodes ($\gamma=0$) are thinner with shorter
spikes as compared to redox electrodes ($\gamma=1$) due to the depletion of
ions in the electrolyte that causes substantial increases in ohmic loss.

By increasing $\delta$ (see right column of Fig. 11), we observe a marked
increase in the number of spikes and the teeth features protrude deeper into
the electrode columns. The limiting kinetics of localized reactions is
compensated by the increased interfacial area with the pure electrolyte
region, allowing for better ion transport.

Conversely, in monolithic designs, electrode utilization is poor. Charging
occurs only in the regions immediately adjacent to the electrolyte interface,
leading to minimal energy storage, as further discussed in Sect. 5.3.
Therefore, we observe exceptional improvements on energy storage across all
cases, ranging from 350% to 750%. And despite significantly larger charging
currents, the ohmic loss is only increased by a maximum of 200%. Considering
that the energy input is much larger due to a lower overall cell resistance,
the energy efficiency for the optimized designs improves significantly.

The second group of results obtained with the modified Bruggeman correlation
is illustrated in Fig. 12 for two values of $\lambda$. Similar to the
observations made in the corresponding original Bruggeman designs (Fig. 10),
going from $\lambda=0.01$ to $\lambda=0.1$, the design complexity increases.
Interestingly, the frequency of the microscale teeth remains roughly the same
between the $\lambda=0.01$ and $\lambda=0.1$ designs with the same $\gamma$.
The teeth in the $\lambda=0.01$ designs are enlarged proportionally with the
macroscopic structures. We observe $500\%$ to $750\%$ improvement in energy
storage with only $\sim 50\%$ increase in ohmic loss.

####  5.2.3 Summary

To summarize the effect of different dimensionless groups on the performance
and optimized design:

  * •

The energy gain is more significant in systems limited by ohmic resistance
(mainly mass transport induced resistance), corresponding to large $\delta$
values, with $\delta$ being the inverse Wagner number.

  * •

Macroscale features are most effective in enhancing overall ionic transport
for porous electrodes that have high electronic conductivity corresponding to
small $\lambda$ values.

  * •

In extreme tortuosity systems characterized by the modified Bruggeman case,
features with spikey teeth are advantageous because they provide secondary
paths for ion transport.

  * •

Interdigitated electrodes reduce the transport distance and increase active-
material utilization.

###  5.3 Performance improvements

Monolithic electrodes

(a) $E_{\mathrm{kin}}=1.1712\times 10^{-1},\ E_{\mathrm{var}}=0.8077$

(b) $E_{\mathrm{ohm}}=8.4885\times 10^{-2},\ R_{\mathrm{avg}}=161.29$

Optimized electrodes

(c) $E_{\mathrm{kin}}=2.2538\times 10^{-1},\ E_{\mathrm{var}}=0.3623$

(d) $E_{\mathrm{ohm}}=4.9681\times 10^{-2},\ R_{\mathrm{avg}}=121.98$

Figure 13: Energy storage density (top) and ohmic loss density (bottom)
distributions of monolithic (left) and optimized (right) electrodes with
$\delta=5$, $\gamma=1$, $\lambda=0.01$, and the original Bruggeman
correlation.

Monolithic electrodes

(a) $E_{\mathrm{kin}}=1.5622\times 10^{-2},\ E_{\mathrm{var}}=3.6411$

(b) $E_{\mathrm{ohm}}=1.5616\times 10^{-2},\ R_{\mathrm{avg}}=920.62$

Optimized electrodes

(c) $E_{\mathrm{kin}}=1.3296\times 10^{-1},\ E_{\mathrm{var}}=0.3164$

(d) $E_{\mathrm{ohm}}=4.1372\times 10^{-2},\ R_{\mathrm{avg}}=192.77$

Figure 14: Energy storage density (top) and ohmic loss density (bottom)
between monolithic (left) and optimized (right) electrodes with $\delta=5$,
$\gamma=1$, $\lambda=0.01$, and the modified Bruggeman correlation.

In this section, we compare the energy storage and ohmic loss between the
optimized and monolithic designs to illustrate the improved material
utilization and ionic transport of the optimized cells. The material
utilization is quantified by the normalized variation in energy storage
density

$E_{\mathrm{var}}=\left\|\bar{\rho}\left(\frac{E(\hat{x})}{\bar{E}}-1\right)\right\|_{L^{2}(\hat{\Omega})},$
(39)

where

$\bar{E}=\left(\int_{\Omega}\bar{\rho}\mathop{}\\!\mathrm{d}\hat{x}\right)^{-1}E_{\mathrm{kin}}$
(40)

is the average energy density. The effectiveness of ionic transport is
quantified by the average cell resistance

$R_{\mathrm{avg}}=\hat{T}^{-1}\int_{0}^{\hat{T}}\frac{\xi\hat{t}}{\int_{\hat{\Gamma}_{c}}\hat{\sigma}\hat{\nabla}\cdot\hat{\Phi}_{1}\cdot\mathbf{n}\mathop{}\\!\mathrm{d}\hat{x}}\mathop{}\\!\mathrm{d}\hat{t}.$
(41)

With a fixed boundary potential across all designs, a smaller cell resistance
results in a stronger charging current and thus more energy storage.

$E_{\mathrm{kin}}$ and $E_{\mathrm{ohm}}$ distributions for the $\delta=5$,
$\gamma=1$, and $\lambda=0.01$ design with either the original or modified
Bruggeman correlations are illustrated in Figs. 13 and 14. First, looking at
the original Bruggeman cases in Fig. 13, notably better electrode utilization
is seen in the optimized designs. Regions near the current collector in the
monolithic electrodes have little energy storage, as opposed to the optimized
electrodes. The optimized designs also achieve a larger energy density at the
electrode tips due to their increased over-potential, which improves their
material utilization by approximately 50%. The monolithic design also exhibits
significant ohmic loss, in the electrode region near the separator, since the
ionic currents must pass these regions to cross the gap. The ohmic loss
density in the optimized electrodes is similar, however, its net effect is
less since the interdigitated structure provides more than twice as much
electrode-electrolyte interface area than the monolithic design. The increased
area provides additional pathways for the currents to cross the gap.
Consequently, the ionic current is more uniformly distributed, resulting in
less ohmic loss. This is also reflected in the 25% decrease in average cell
resistance.

Analogous $E_{\mathrm{kin}}$ and $E_{\mathrm{ohm}}$ plots for designs obtained
with modified Bruggeman correlation are provided in Fig. 14. We observe poor
material utilization by the monolithic electrodes wherein most regions are
considered “dead zones” with negligible energy storage. On the other hand, the
spikey teethed optimized structures promote ion exchange, thereby the energy
density is more uniformly distributed throughout the electrodes. Therefore,
the overall storage is enhanced by an order of magnitude and the material
utilization is enhanced by $\sim 90\%$. The ohmic loss distribution is similar
to the original Bruggeman case. Although the maximum loss per unit area is
significantly less for the optimized design, the total ohmic loss increases by
$\sim 150\%$ due to the substantial increase in energy input. Nevertheless, we
observe 5x reduction in average cell resistance. Also, the energy efficiency
reaches 76%, a 26% improvement over the monolithic design.

###  5.4 3D example

Figure 15: Optimized anode and its energy storage distribution.

Figure 16: Cross sections of both electrodes for 3D optimized designs. The
blue (top) electrode is the cathode, and the red/pink (bottom) electrode is
the anode.

Figure 17: Zoomed in view of the interlocking structure formed by the
electrode arches. Figure 18: Schematic of electrode structures. Both 2D and 3D
cases form interdigitated columns. In addition, 3D designs also possess
interlocked structure.

(a) Monolithic design

(b) Optimized design

Figure 19: Cross sections of the energy storage distribution in the 3D
monolithic (a) and optimized (b) designs.

We optimize a 3D design with $\delta=2$, $\gamma=0.5$, $\lambda=0.01$, and
modified Bruggeman correlation. The two optimized electrodes are similar and
so we only illustrate the anode and its energy storage distribution, cf. Fig.
15. Similar protruding features of the cathode are interdigitated with these
of the anode. Throughout the design of both electrodes, we observe small
penetrating holes, which increase interfacial area similar to the spikey teeth
from the 2D designs. From the cross sections of both electrodes plotted in
Fig. 16, we observe the interdigitation similar to the 2D designs. In
addition, connections between the main pillar and the wall are observed in
certain locations, resulting in an arch shape. The space below the arch is
filled with the arch from the other electrode, creating an interlocking
structure as illustrated by the zoomed in view of Fig. 17. A simplified
schematic summarizing the different design features is presented in Fig. 18

We also compare the performance between the optimized and monolithic
electrodes. The energy storage density distribution over a cross section is
presented in Fig. 19 to demonstrate superior material utilization of the
optimized design. Similar to the 2D observations, only the regions neighboring
the separator have noticeable energy storage in the monolithic design. Better
material utilization is observed in the optimized design, wherein the stored
energy is more uniformly distributed. The total energy storage and ohmic loss
is $814\%$ and $145\%$ to the ones from the monolith, which is slightly better
than the corresponding 2D case, possibly due to three-dimensional effects such
as interlocking.

##  6 Conclusion

In this work, we present a density-based topology optimization strategy for
the design of porous electrodes in electrochemical energy storage devices with
Faradaic reactions and capacitive storage. A full-cell model is utilized to
simultaneously optimize the cathode and anode.

We present various 2D optimized designs to analyze the changes in their
features as we vary the value of dimensionless groups. Interdigitation is
observed in all optimized designs, which unlocks additional energy storage by
allowing reactions to occur closer to the opposite current collectors and by
enhancing ionic transport between the electrodes. In general, more energy
storage is acquired through electrode structure for $\delta=5$ cases where
reactions are localized near the separator interface. Among these cases, the
designs also vary noticeably when changing between Faradaic ($\gamma=1$) and
capacitive storage-mechanisms ($\gamma=0$). $\lambda$, i.e. the ratio between
ionic and total conductivity, controls the interdigitation depth. We observe
more shorter protrusions in designs with lower electronic conduction, i.e. as
$\lambda\rightarrow 0.1$. Designs generated using the standard Bruggeman
correlation for effective transport possess smooth surfaces. The designs
generated with modified Bruggeman correlation, where effective ionic transport
is much slower, possess spiky teeth cutting into the main pillars to further
boost the ion exchange between the electrode and the pure electrolyte.

The key performance metrics of the optimized designs are compared against the
ones from monolithic electrodes to quantify the improvements. We observe
increased energy storage for all cases with the best improvement exceeding
700% from the modified Bruggeman correlation cases. Comparisons of energy
storage and ohmic loss distributions suggest superior material utilization in
the optimized designs. In addition, the ohmic loss is more evenly distributed
and thus does not grow significantly even when the total charging current
increases by an order of magnitude.

Finally, a 3D design is generated to further demonstrate the utility of the
proposed optimization strategy. More complex features, such as channels,
holes, and interlocking arches between the electrodes are observed in this
design.

As mentioned above, our current work does not track surface or solid
concentrations. To address this, our future work may incorporate diffusion in
the solid particles as in the Doyle-Fuller-Newman model (Doyle et al, 1993).
Typically used for a 1D geometry, this model includes diffusion inside
spherical particles, making it a pseudo-2D model. For 3D geometries, this
becomes a pseudo-4D model, further exacerbating the computational effort.

## Funding information

This work was performed under the auspices of the U.S. Department of Energy by
Lawrence Livermore National Laboratory under Contract DE-AC52-07NA27344 and
was supported by the LLNL-LDRD program under project numbers 20-ERD-019 and
23-SI-002. LLNL Release Number LLNL-JRNL-862205-DRAFT.

##  Appendix A Energy balance derivation

We derive the energy balance for the dimensionless system with energy input,
energy storage, and ohmic loss. We first scale the electric and ionic
potential equation such that the source terms are the same:

$\displaystyle-\frac{1}{\lambda}\hat{\nabla}\cdot\big{(}\hat{\sigma}\hat{\nabla}\hat{\Phi}_{1}\big{)}=$
$\displaystyle-\delta_{r}\hat{a}\hat{i}_{n}-\delta_{c}\hat{a}\hat{i}_{c},$
(42)
$\displaystyle-\frac{1}{1-\lambda}\hat{\nabla}\cdot(\hat{D}\hat{c}\hat{\nabla}\hat{\Phi}_{2})-\frac{1}{1-\lambda}$
$\displaystyle\left(\frac{t_{0,+}}{z_{+}}+\frac{t_{0,-}}{z_{-}}\right)$
$\displaystyle\hat{\nabla}\cdot(\hat{D}\hat{\nabla}\hat{c})=$
$\displaystyle\delta_{r}\hat{a}\hat{i}_{n}+\delta_{c}\hat{a}\hat{i}_{c}.$

Since $\hat{\Phi}_{1}=0$ on $\hat{\Gamma}_{a}$ and
$\nabla\hat{\Phi}_{1}\cdot\mathbf{n}=0$ on
$\partial\hat{\Omega}\setminus(\hat{\Gamma}_{a}\bigcup\hat{\Gamma}_{c})$ the
dimensionless energy input can be expressed as

$\displaystyle E_{\mathrm{in}}$
$\displaystyle=\frac{1}{\lambda}\int_{0}^{\hat{T}}\int_{\hat{\Gamma}_{c}}\hat{\sigma}\hat{\nabla}\hat{\Phi}_{1}\cdot\mathbf{n}\hat{\Phi}_{1}\mathop{}\\!\mathrm{d}\hat{s}\mathop{}\\!\mathrm{d}\hat{t}$
(43)
$\displaystyle=\frac{1}{\lambda}\int_{0}^{\hat{T}}\int_{\partial\hat{\Omega}}\hat{\sigma}\hat{\nabla}\hat{\Phi}_{1}\cdot\mathbf{n}\hat{\Phi}_{1}\mathop{}\\!\mathrm{d}\hat{s}\mathop{}\\!\mathrm{d}\hat{t},$

which extends integral from the cathode current collector $\hat{\Gamma}_{c}$
to the entire boundary $\partial\hat{\Omega}$. With
$\hat{\nabla}\hat{\Phi}_{2}\cdot\mathbf{n}=0$ and
$\hat{\nabla}\hat{c}\cdot\mathbf{n}=0$ on $\partial\hat{\Omega}$, the
dimensionless energy input can be further extended as

$\displaystyle
E_{\mathrm{in}}=\frac{1}{\lambda}\int_{0}^{\hat{T}}\int_{\partial\hat{\Omega}}\hat{\sigma}\hat{\nabla}\hat{\Phi}_{1}\cdot\mathbf{n}\hat{\Phi}_{1}\mathop{}\\!\mathrm{d}\hat{s}\mathop{}\\!\mathrm{d}\hat{t}$
(44)
$\displaystyle+\frac{1}{1-\lambda}\int_{0}^{\hat{T}}\int_{\partial\hat{\Omega}}\hat{D}\hat{c}\hat{\nabla}\hat{\Phi}_{2}\cdot\mathbf{n}\hat{\Phi}_{2}\mathop{}\\!\mathrm{d}\hat{s}\mathop{}\\!\mathrm{d}\hat{t}$
$\displaystyle+\frac{1}{1-\lambda}\left(\frac{t_{0,+}}{z_{+}}+\frac{t_{0,-}}{z_{-}}\right)\int_{0}^{\hat{T}}\int_{\partial\hat{\Omega}}\hat{D}\hat{\nabla}\hat{c}\cdot\mathbf{n}\hat{\Phi}_{2}\mathop{}\\!\mathrm{d}\hat{s}\mathop{}\\!\mathrm{d}\hat{t},$

Applying the divergence theorem we obtain

$\displaystyle
E_{\mathrm{in}}=\frac{1}{\lambda}\int_{0}^{\hat{T}}\int_{\hat{\Omega}}\hat{\nabla}\cdot\left(\hat{\sigma}\hat{\nabla}\hat{\Phi}_{1}\hat{\Phi}_{1}\right)\mathop{}\\!\mathrm{d}\hat{x}\mathop{}\\!\mathrm{d}\hat{t}$
(45)
$\displaystyle+\frac{1}{1-\lambda}\int_{0}^{\hat{T}}\int_{\hat{\Omega}}\hat{\nabla}\cdot\left(\hat{D}\hat{c}\hat{\nabla}\hat{\Phi}_{2}\hat{\Phi}_{2}\right)\mathop{}\\!\mathrm{d}\hat{x}\mathop{}\\!\mathrm{d}\hat{t}$
$\displaystyle+\frac{1}{1-\lambda}\left(\frac{t_{0,+}}{z_{+}}+\frac{t_{0,-}}{z_{-}}\right)\int_{0}^{\hat{T}}\int_{\hat{\Omega}}\hat{\nabla}\cdot\left(\hat{D}\hat{\nabla}\hat{c}\hat{\Phi}_{2}\right)\mathop{}\\!\mathrm{d}\hat{x}\mathop{}\\!\mathrm{d}\hat{t}.$

An application of the product rule subsequently gives

$\displaystyle
E_{\mathrm{in}}=\frac{1}{\lambda}\int_{0}^{\hat{T}}\int_{\hat{\Omega}}\hat{\sigma}\hat{\nabla}\hat{\Phi}_{1}\cdot\hat{\nabla}\hat{\Phi}_{1}\mathop{}\\!\mathrm{d}\hat{x}\mathop{}\\!\mathrm{d}\hat{t}$
(46)
$\displaystyle+\frac{1}{\lambda}\int_{0}^{\hat{T}}\int_{\hat{\Omega}}\hat{\nabla}\cdot(\hat{\sigma}\hat{\nabla}\hat{\Phi}_{1})\hat{\Phi}_{1}\mathop{}\\!\mathrm{d}\hat{x}\mathop{}\\!\mathrm{d}\hat{t}$
$\displaystyle+\frac{1}{1-\lambda}\int_{0}^{\hat{T}}\int_{\hat{\Omega}}\hat{D}\hat{c}\hat{\nabla}\hat{\Phi}_{2}\cdot\hat{\nabla}\hat{\Phi}_{2}\mathop{}\\!\mathrm{d}\hat{x}\mathop{}\\!\mathrm{d}\hat{t}$
$\displaystyle+\frac{1}{1-\lambda}\int_{0}^{\hat{T}}\int_{\hat{\Omega}}\hat{\nabla}\cdot(\hat{D}\hat{c}\hat{\nabla}\hat{\Phi}_{2})\hat{\Phi}_{2}\mathop{}\\!\mathrm{d}\hat{x}\mathop{}\\!\mathrm{d}\hat{t}$
$\displaystyle+\frac{1}{1-\lambda}\left(\frac{t_{0,+}}{z_{+}}+\frac{t_{0,-}}{z_{-}}\right)\int_{0}^{\hat{T}}\int_{\hat{\Omega}}\hat{D}\hat{\nabla}\hat{c}\cdot\hat{\nabla}\hat{\Phi}_{2}\mathop{}\\!\mathrm{d}\hat{x}\mathop{}\\!\mathrm{d}\hat{t}$
$\displaystyle+\frac{1}{1-\lambda}\left(\frac{t_{0,+}}{z_{+}}+\frac{t_{0,-}}{z_{-}}\right)\int_{0}^{\hat{T}}\int_{\hat{\Omega}}\hat{\nabla}\cdot(\hat{D}\hat{\nabla}\hat{c})\hat{\Phi}_{2}\mathop{}\\!\mathrm{d}\hat{x}\mathop{}\\!\mathrm{d}\hat{t}.$

Combining (42) and (46) yields

$\displaystyle
E_{\mathrm{in}}=\frac{1}{\lambda}\int_{0}^{\hat{T}}\int_{\hat{\Omega}}\hat{\sigma}\hat{\nabla}\hat{\Phi}_{1}\cdot\hat{\nabla}\hat{\Phi}_{1}\mathop{}\\!\mathrm{d}\hat{x}\mathop{}\\!\mathrm{d}\hat{t}$
(47)
$\displaystyle+\frac{1}{1-\lambda}\int_{0}^{\hat{T}}\int_{\hat{\Omega}}\hat{D}\hat{c}\hat{\nabla}\hat{\Phi}_{2}\cdot\hat{\nabla}\hat{\Phi}_{2}\mathop{}\\!\mathrm{d}\hat{x}\mathop{}\\!\mathrm{d}\hat{t}$
$\displaystyle+\frac{1}{1-\lambda}\left(\frac{t_{+}}{z_{+}}+\frac{t_{-}}{z_{-}}\right)\int_{0}^{\hat{T}}\int_{\hat{\Omega}}\hat{D}\hat{\nabla}\hat{c}\cdot\hat{\nabla}\hat{\Phi}_{2}\mathop{}\\!\mathrm{d}\hat{x}\mathop{}\\!\mathrm{d}\hat{t}$
$\displaystyle+\delta_{r}\int_{0}^{\hat{T}}\int_{\hat{\Omega}}\hat{a}\hat{i}_{n}(\hat{\Phi}_{1}-\hat{\Phi}_{2})\mathop{}\\!\mathrm{d}\hat{x}\mathop{}\\!\mathrm{d}\hat{t}$
$\displaystyle+\delta_{c}\int_{0}^{\hat{T}}\int_{\hat{\Omega}}\hat{a}\hat{i}_{c}(\hat{\Phi}_{1}-\hat{\Phi}_{2})\mathop{}\\!\mathrm{d}\hat{x}\mathop{}\\!\mathrm{d}\hat{t},$

which is equivalent to $E_{\mathrm{in}}=E_{\mathrm{kin}}+E_{\mathrm{ohm}}$.

##  Appendix B Mixed discretizations

A mixed formulation is utilized to solve the filter problem (2) to obtain
piecewise constant $\tilde{\rho}$ over triangular and tetrahedral meshes. In
this way, we are sure that $\tilde{\rho}(\mathbf{x})\in[0,1]$. Indeed, the
mixed method enforces flux continuity across the element interfaces, which
prevents oscillations in the filtered field as the control field goes to
$\\{0,1\\}$ (Salazar de Troya and Tortorelli, 2020). As opposed to the two-
point flux approximation scheme from Roy et al (2022), no orthogonality
condition is required on the mesh regularity.

With $r^{2}>0$, we define the auxiliary flux term as

$r^{-2}\bm{u}_{f}=-\nabla\tilde{\rho}.$ (48)

Then the filter problem is rewritten as

$\displaystyle r^{-2}\bm{u}_{f}+\nabla\tilde{\rho}$ $\displaystyle=0\quad$
$\displaystyle\text{in }\Omega,$ (49)
$\displaystyle-\nabla\cdot\bm{u}_{f}-\tilde{\rho}$ $\displaystyle=-\rho\quad$
$\displaystyle\text{in }\Omega,$ $\displaystyle\bm{u}_{f}\cdot\mathbf{n}$
$\displaystyle=0\quad$ $\displaystyle\text{on }\partial\Omega.$

The function spaces for the mixed weak variational form are

$\bm{V}=\left\\{\bm{v}\in
H(\mathrm{div},\Omega)\;\middle|\;\bm{v}\cdot\mathbf{n}=0\text{ on
}\partial\Omega\right\\},\;\;W=L^{2}(\Omega),$ (50)

where the Neumann flux boundary condition are imposed strongly. We seek to
find $\left(\bm{u}_{f},\tilde{\rho}\right)\in\bm{V}\times W$ such that

$\displaystyle(r^{-2}\bm{u}_{f},\bm{v})_{\Omega}$
$\displaystyle-(\tilde{\rho},\nabla\cdot\bm{v})_{\Omega}=0,$ (51)
$\displaystyle-(\nabla\cdot\bm{u}_{f},w)_{\Omega}$
$\displaystyle-(\tilde{\rho},w)_{\Omega}=-(\rho,w)_{\Omega},$

$\forall\left(\bm{v},w\right)\in\bm{V}\times W$.

A similar mixed formulation is used for solving the boundary propagation PDE
(4). From (4), since $\bar{\rho}$ is strictly positive, we define the
auxiliary flux term as

$\bar{\rho}^{-1}\bm{u}_{b}=-\nabla\beta.$ (52)

Then the continuous boundary propagation problem can be reformulated as

$\displaystyle\bar{\rho}^{-1}\bm{u}_{b}+\nabla\beta$ $\displaystyle=0\quad$
$\displaystyle\text{in }\Omega,$ (53)
$\displaystyle-\nabla\cdot\bm{u}_{b}-(1-\bar{\rho})\beta$
$\displaystyle=0\quad$ $\displaystyle\text{in }\Omega,$ $\displaystyle\beta$
$\displaystyle=1\quad$ $\displaystyle\text{on }\Gamma_{\mathrm{a}},$
$\displaystyle\beta$ $\displaystyle=-1\quad$ $\displaystyle\text{on
}\Gamma_{\mathrm{c}},$ $\displaystyle\bm{u}_{b}\cdot\mathbf{n}$
$\displaystyle=0\quad$ $\displaystyle\text{on
}\partial\Omega/\Gamma_{\mathrm{c}}\cup\Gamma_{\mathrm{a}}.$

Let $\Gamma_{D}=\Gamma_{a}\cup\Gamma_{c}$. We introduce the function space

$\bm{V}_{D}=\left\\{\bm{v}\in
H(\mathrm{div},\Omega)\middle|\bm{v}\cdot\mathbf{n}=0\text{ on
}\Gamma_{D}\right\\}.$ (54)

We seek to find $\left(\bm{u}_{b},\beta\right)\in\bm{V}_{D}\times W$ such that

$\displaystyle(\bar{\rho}^{-1}\bm{u}_{b},\bm{v})_{\Omega}-(\beta,\nabla\cdot\bm{v})_{\Omega}$
$\displaystyle=$ (55) $\displaystyle-\langle
1,\bm{v}\cdot\mathbf{n}\rangle_{\Gamma_{\mathrm{a}}}-$
$\displaystyle\langle-1,\bm{v}\cdot\mathbf{n}\rangle_{\Gamma_{\mathrm{c}}},$
$\displaystyle-(\nabla\cdot\bm{u}_{b},w)_{\Omega}-((1-\bar{\rho})\beta,w)_{\Omega}$
$\displaystyle=0,$

$\forall\left(\bm{v},w\right)\in\bm{V}_{D}\times W$. Both problems use
Raviart-Thomas basis functions with the lowest order to obtain the discretized
solution, which results in piecewise constant $\tilde{\rho}$/$\beta$ and
piecewise linear $\bm{u}_{f}$/$\bm{u}_{b}$ with continuous normal flux across
the element interfaces.

For the three-dimensional case, we use an iterative method to solve the linear
systems resulting from the weak forms described by (51) and (55). These
systems have symmetric block matrices of the form

$\begin{bmatrix}A&B\\\ B^{\top}&D\\\ \end{bmatrix}$ (56)

where $A$ represents a weighted mass matrix for the flux, $D$ is a weighted
mass matrix for $\tilde{\rho}$ or $\beta$, and $B$ is the coupling term. The
approach introduced in Silvester and Wathen (1994) is used to solve the block
system efficiently. The block matrix (56) admits a factorization of the type

$\begin{bmatrix}A&B\\\ B^{\top}&D\\\ \end{bmatrix}=\begin{bmatrix}I&0\\\
B^{\top}A^{-1}&I\end{bmatrix}\begin{bmatrix}A&0\\\
0&S\end{bmatrix}\begin{bmatrix}I&A^{-1}B\\\ 0&I\end{bmatrix},$ (57)

where $S=D-CA^{-1}B$ is called the _Schur complement_. The inverse of the
block matrix is

$\begin{bmatrix}A&B\\\ B^{\top}&D\\\ \end{bmatrix}^{-1}\\\
=\begin{bmatrix}I&-A^{-1}B\\\ 0&I\end{bmatrix}\begin{bmatrix}A^{-1}&0\\\
0&S^{-1}\end{bmatrix}\begin{bmatrix}I&0\\\ -B^{\top}A^{-1}&I\end{bmatrix},$
(58)

which requires inverting the blocks $A$ and $S$. Our preconditioner
approximates the blocks of this factorization. We approximate the Schur
complements using the diagonal of $A$, i.e. we use the sparse
$\tilde{S}=D-C\mathrm{diag}(A)^{-1}B$. The inverse of $\tilde{S}$ is
approximated by Algebraic Multigrid (Ruge and Stüben, 1987), and the inverse
of $A$ is approximated by a block incomplete LU factorization with 0 fill
(ILU0) (Meijerink and Van Der Vorst, 1977). The PETSc solver options used for
the three-dimensional simulations are given in Appendix C.

##  Appendix C PETSc solver options

Here we provide he PETSc solver options used for the three-dimensional case.
The options for solving the governing physical system are:

⬇

1 "snes_rtol": 0,

2 "snes_atol": 1e-4,

3 "mat_type": "aij",

4 "ksp_type": "gmres",

5 "ksp_rtol": 1e-4,

6 "ksp_gmres_restart": 30,

7 "pc_type": "fieldsplit",

8 "pc_fieldsplit_type": "multiplicative",

9 "fieldsplit_pc_type": "hypre",

10 "fieldsplit_pc_hypre_boomeramg": {

11 "strong_threshold": 0.7,

12 "coarsen_type": "HMIS",

13 "agg_nl": 3,

14 "interp_type": "ext+i",

15 "agg_num_paths": 5,

16 },

The options for the mixed formulation used for the filtering and boundary
propagation problems are:

⬇

1 "mat_type": "nest",

2 "ksp_type": "gmres",

3 "ksp_rtol": 1e-4,

4 "pc_type": "fieldsplit",

5 "pc_fieldsplit_type": "schur",

6 "pc_fieldsplit_schur_fact_type": "full",

7 "fieldsplit_0_ksp_type": "preonly",

8 "fieldsplit_0_pc_type": "bjacobi",

9 "fieldsplit_0_sub_pc_type": "ilu",

10 "fieldsplit_1_ksp_type": "preonly",

11 "pc_fieldsplit_schur_precondition": "selfp",

12 "fieldsplit_1_pc_type": "hypre",

13 "fieldsplit_1_pc_hypre_boomeramg": {

14 "strong_threshold": 0.7,

15 "coarsen_type": "HMIS",

16 "agg_nl": 3,

17 "interp_type": "ext+i",

18 "agg_num_paths": 5,

19 }.

## References

  *   * Alizadeh et al (2023) Alizadeh M, Charoen-amornkitt P, Suzuki T, et al (2023) A numerical simulation of evolution processes and entropy generation for optimal architecture of an electrochemical reaction-diffusion system: Comparison of two optimization strategies. Journal of The Electrochemical Society 10.1149/1945-7111/ad0a7c, URL http://iopscience.iop.org/article/10.1149/1945-7111/ad0a7c
  * Balay et al (2023a) Balay S, Abhyankar S, Adams MF, et al (2023a) PETSc/TAO users manual. Tech. Rep. ANL-21/39 - Revision 3.20, Argonne National Laboratory, 10.2172/1968587
  * Balay et al (2023b) Balay S, Abhyankar S, Adams MF, et al (2023b) PETSc Web page. https://petsc.org/, URL https://petsc.org/
  * Batista et al (2023) Batista MDR, Chandrasekaran S, Moran BD, et al (2023) Design and additive manufacturing of optimized electrodes for energy storage applications. Carbon 205:262–269. 10.1016/j.carbon.2023.01.044
  * Beck et al (2021a) Beck VA, Ivanovskaya AN, Chandrasekaran S, et al (2021a) Inertially enhanced mass transport using 3D-printed porous flow-through electrodes with periodic lattice structures. Proceedings of the National Academy of Sciences 118(32):e2025562,118. 10.1073/pnas.2025562118
  * Beck et al (2021b) Beck VA, Wong JJ, Jekel CF, et al (2021b) Computational design of microarchitected porous electrodes for redox flow batteries. Journal of Power Sources 512:230,453. 10.1016/j.jpowsour.2021.230453
  * Behrou et al (2019) Behrou R, Pizzolato A, Forner-Cuenca A (2019) Topology optimization as a powerful tool to design advanced PEMFCs flow fields. International Journal of Heat and Mass Transfer 135:72–92. 10.1016/j.ijheatmasstransfer.2019.01.050
  * Bendsøe (1989) Bendsøe MP (1989) Optimal shape design as a material distribution problem. Structural optimization 1(4):193–202. 10.1007/bf01650949
  * Bourdin (2001) Bourdin B (2001) Filters in topology optimization. International Journal for Numerical Methods in Engineering 50(9):2143–2158
  * Brandt (1977) Brandt A (1977) Multi-level adaptive solutions to boundary-value problems. Mathematics of Computation 31(138):333–390. 10.1090/s0025-5718-1977-0431719-x
  * Brosa Planella et al (2022) Brosa Planella F, Ai W, Boyce AM, et al (2022) A continuum of physics-based lithium-ion battery models reviewed. Progress in Energy 4(4):042,003. 10.1088/2516-1083/ac7d31
  * Bruggeman (1935) Bruggeman VD (1935) Berechnung verschiedener physikalischer konstanten von heterogenen substanzen. i. dielektrizitätskonstanten und leitfähigkeiten der mischkörper aus isotropen substanzen. Annalen der physik 416(7):636–664. 10.1002/andp.19374210205
  * Bruns and Tortorelli (2001) Bruns TE, Tortorelli DA (2001) Topology optimization of non-linear elastic structures and compliant mechanisms. Computer Methods in Applied Mechanics and Engineering 190(26-27):3443–3459
  * Chandrasekaran et al (2023) Chandrasekaran S, Lin D, Li Y, et al (2023) Aerogels, additive manufacturing, and energy storage. Joule 7(5):866–883
  * Chen et al (2019) Chen CH, Yaji K, Yamasaki S, et al (2019) Computational design of flow fields for vanadium redox flow batteries via topology optimization. Journal of Energy Storage 26:100,990. 10.1016/j.est.2019.100990
  * Chen et al (2020) Chen KH, Namkoong MJ, Goel V, et al (2020) Efficient fast-charging of lithium-ion batteries enabled by laser-patterned three-dimensional graphite anode architectures. Journal of Power Sources 471:228,475
  * Chu and Majumdar (2012) Chu S, Majumdar A (2012) Opportunities and challenges for a sustainable energy future. Nature 488(7411):294–303. 10.1038/nature11475
  * Chu et al (2016) Chu S, Cui Y, Liu N (2016) The path towards sustainable energy. Nature Materials 16(1):16–22. 10.1038/nmat4834
  * Chu et al (2021) Chu T, Park S, Fu KK (2021) 3D printing-enabled advanced electrode architecture design. Carbon Energy 3(3):424–439
  * Doyle et al (1993) Doyle M, Fuller TF, Newman J (1993) Modeling of galvanostatic charge and discharge of the lithium/polymer/insertion cell. Journal of The Electrochemical Society 140(6):1526. 10.1149/1.2221597
  * Falgout and Yang (2002) Falgout RD, Yang UM (2002) hypre: A library of high performance preconditioners. In: International Conference on Computational Science, Springer, pp 632–641, 10.1007/3-540-47789-6_66
  * Fuller and Harb (2018) Fuller TF, Harb JN (2018) Electrochemical Engineering. John Wiley & Sons, Hoboken, NJ USA
  * Gallagher et al (2015) Gallagher KG, Trask SE, Bauer C, et al (2015) Optimizing areal capacities through understanding the limitations of lithium-ion electrodes. Journal of The Electrochemical Society 163(2):A138
  * Gao et al (2023) Gao Y, Cao Q, Pu J, et al (2023) Stable Zn anodes with triple gradients. Advanced Materials 35(6):2207,573
  * Golmon et al (2012) Golmon S, Maute K, Dunn ML (2012) Multiscale design optimization of lithium ion batteries using adjoint sensitivity analysis. International Journal for Numerical Methods in Engineering 92(5):475–494. 10.1002/nme.4347
  * Griewank and Walther (2000) Griewank A, Walther A (2000) Algorithm 799: Revolve: An implementation of checkpointing for the reverse or adjoint mode of computational differentiation. ACM Trans Math Softw 26(1):19–45. 10.1145/347837.347846
  * Guest et al (2004) Guest JK, Prévost JH, Belytschko T (2004) Achieving minimum length scale in topology optimization using nodal design variables and projection functions. International journal for numerical methods in engineering 61(2):238–254. 10.1002/nme.1064
  * Gür (2018) Gür TM (2018) Review of electrical energy storage technologies, materials and systems: challenges and prospects for large-scale grid storage. Energy & Environmental Science 11(10):2696–2767. 10.1039/c8ee01419a
  * Ham et al (2023) Ham DA, Kelly PHJ, Mitchell L, et al (2023) Firedrake User Manual. Imperial College London and University of Oxford and Baylor University and University of Washington, first edition edn., 10.25561/104839
  * Hart et al (2003) Hart RW, White HS, Dunn B, et al (2003) 3-D microbatteries. Electrochemistry Communications 5(2):120–123
  * Henson and Yang (2002) Henson VE, Yang UM (2002) BoomerAMG: A parallel algebraic multigrid solver and preconditioner. Applied Numerical Mathematics 41(1):155–177. 10.1016/s0168-9274(01)00115-5
  * Heubner et al (2021) Heubner C, Voigt K, Marcinkowski P, et al (2021) From active materials to battery cells: A straightforward tool to determine performance metrics and support developments at an application-relevant level. Advanced Energy Materials 11(46):2102,647
  * Hung et al (2021) Hung CH, Allu S, Cobb CL (2021) Modeling current density non-uniformities to understand high-rate limitations in 3D interdigitated lithium-ion batteries. Journal of the Electrochemical Society 168(10):100,512
  * Hung et al (2022) Hung CH, Huynh P, Teo K, et al (2022) Are three-dimensional batteries beneficial? analyzing historical data to elucidate performance advantages. ACS Energy Letters 8(1):296–305
  * Johnson and Newman (1971) Johnson A, Newman J (1971) Desalting by means of porous carbon electrodes. Journal of The Electrochemical Society 118(3):150. 10.1149/1.2408094
  * Kuang et al (2019) Kuang Y, Chen C, Kirsch D, et al (2019) Thick electrode batteries: Principles, opportunities, and challenges. Advanced Energy Materials 9(33):1901,457. https://doi.org/10.1002/aenm.201901457, URL https://onlinelibrary.wiley.com/doi/abs/10.1002/aenm.201901457
  * Lazarov and Sigmund (2011) Lazarov BS, Sigmund O (2011) Filters in topology optimization based on Helmholtz-type differential equations. International Journal for Numerical Methods in Engineering 86(6):765–781. 10.1002/nme.3072
  * Lazarov et al (2016) Lazarov BS, Wang F, Sigmund O (2016) Length scale and manufacturability in density-based topology optimization. Archive of Applied Mechanics 86:189–218. 10.1007/s00419-015-1106-4
  * Lin et al (2022) Lin TY, Baker SE, Duoss EB, et al (2022) Topology optimization of 3D flow fields for flow batteries. Journal of The Electrochemical Society 169(5):050,540. 10.1149/1945-7111/ac716d
  * Long et al (2004) Long JW, Dunn B, Rolison DR, et al (2004) Three-dimensional battery architectures. Chemical Reviews 104(10):4463–4492
  * Long et al (2020) Long JW, Dunn B, Rolison DR, et al (2020) 3D architectures for batteries and electrodes. Advanced Energy Materials 10(46):2002,457
  * Lyu et al (2021) Lyu Z, Lim G, Koh JJ, et al (2021) Design and manufacture of 3D-printed batteries. Joule 5:89–114
  * Madabattula and Kumar (2020) Madabattula G, Kumar S (2020) Model and measurement based insights into double layer capacitors: Voltage-dependent capacitance and low ionic conductivity in pores. Journal of The Electrochemical Society 167(8):080,535. 10.1149/1945-7111/ab90aa
  * Meijerink and Van Der Vorst (1977) Meijerink JA, Van Der Vorst HA (1977) An iterative solution method for linear systems of which the coefficient matrix is a symmetric _M_ -matrix. Mathematics of computation 31(137):148–162. 10.2307/2005786
  * Mitusch et al (2019) Mitusch SK, Funke SW, Dokken JS (2019) dolfin-adjoint 2018.1: automated adjoints for FEniCS and Firedrake. Journal of Open Source Software 4(38):1292. 10.21105/joss.01292
  * Newman and Thomas-Alyea (2012) Newman J, Thomas-Alyea KE (2012) Electrochemical systems. John Wiley & Sons
  * Newman and Tiedemann (1975) Newman J, Tiedemann W (1975) Porous-electrode theory with battery applications. AlChE Journal 21(1):25–41. 10.1002/aic.690210103
  * Newman and Tobias (1962) Newman JS, Tobias CW (1962) Theoretical analysis of current distribution in porous electrodes. Journal of The Electrochemical Society 109(12):1183. 10.1149/1.2425269
  * Qu et al (2023) Qu C, Zhang Y, Zhang Z, et al (2023) Design of radial flow channel proton exchange membrane fuel cell based on topology optimization. Processes 11(8):2482. 10.3390/pr11082482
  * Ramadesigan et al (2010) Ramadesigan V, Methekar RN, Latinwo F, et al (2010) Optimal porosity distribution for minimized ohmic drop across a porous electrode. Journal of The Electrochemical Society 157(12):A1328. 10.1149/1.3495992
  * Roy et al (2022) Roy T, Salazar de Troya MA, Worsley MA, et al (2022) Topology optimization for the design of porous electrodes. Structural and Multidisciplinary Optimization 65(171). 10.1007/s00158-022-03249-2
  * Ruge and Stüben (1987) Ruge JW, Stüben K (1987) Algebraic multigrid. In: Multigrid methods, vol 3 of Frontiers in Applied Mathematics. SIAM, Philadelphia, chap 4, p 73–130, 10.1137/1.9781611971057.ch4
  * Sá et al (2022) Sá LF, Okubo Jr. CM, Sá AN, et al (2022) Continuous boundary condition propagation model for topology optimization. Structural and Multidisciplinary Optimization 65(15). 10.1007/s00158-021-03148-y
  * Saad and Schultz (1986) Saad Y, Schultz MH (1986) Gmres: A generalized minimal residual algorithm for solving nonsymmetric linear systems. SIAM Journal on scientific and statistical computing 7(3):856–869. 10.1137/0907058
  * Sigmund (2007) Sigmund O (2007) Morphology-based black and white filters for topology optimization. Structural and Multidisciplinary Optimization 33:401–424. 10.1007/s00158-006-0087-x
  * Sigmund and Maute (2013) Sigmund O, Maute K (2013) Topology optimization approaches: A comparative review. Structural and multidisciplinary optimization 48(6):1031–1055. 10.1007/s00158-013-0978-6
  * Sigmund and Torquato (1997) Sigmund O, Torquato S (1997) Design of materials with extreme thermal expansion using a three-phase topology optimization method. Journal of the Mechanics and Physics of Solids 45(6):1037–1067. 10.1016/S0022-5096(96)00114-7
  * Silvester and Wathen (1994) Silvester D, Wathen A (1994) Fast iterative solution of stabilised stokes systems part II: Using general block preconditioners. SIAM Journal on Numerical Analysis 31(5):1352–1367. 10.1137/0731070
  * Simon et al (2014) Simon P, Gogotsi Y, Dunn B (2014) Where do batteries end and supercapacitors begin? Science 343(6176):1210–1211
  * Singh et al (2015) Singh M, Kaiser J, Hahn H (2015) Thick electrodes for high energy lithium ion batteries. Journal of The Electrochemical Society 162(7):A1196
  * Soares et al (2021) Soares DM, Ren Z, Mujib SB, et al (2021) Additive manufacturing of electrochemical energy storage systems electrodes. Advanced Energy and Sustainability Research 2(5):2000,111
  * Svanberg (1987) Svanberg K (1987) The method of moving asymptotes – a new method for structural optimization. International journal for numerical methods in engineering 24(2):359–373. 10.1002/nme.1620240207
  * Tjaden et al (2016) Tjaden B, Cooper SJ, Brett DJ, et al (2016) On the origin and application of the Bruggeman correlation for analysing transport phenomena in electrochemical systems. Current Opinion in Chemical Engineering 12:44–51. 10.1016/j.coche.2016.02.006
  * Tjaden et al (2018) Tjaden B, Brett DJ, Shearing PR (2018) Tortuosity in electrochemical devices: a review of calculation approaches. International Materials Reviews 63(2):47–67. 10.1080/09506608.2016.1249995
  * Trembacki et al (2019) Trembacki B, Duoss E, Oxberry G, et al (2019) Mesoscale electrochemical performance simulation of 3D interpenetrating lithium-ion battery electrodes. Journal of The Electrochemical Society 166(6):A923—A934
  * Salazar de Troya (2021) Salazar de Troya MA (2021) LLNL/pyMMAopt: Method of Moving Asymptotes for Firedrake. 10.5281/zenodo.5524961
  * Salazar de Troya and Tortorelli (2020) Salazar de Troya MA, Tortorelli DA (2020) Three-dimensional adaptive mesh refinement in stress-constrained topology optimization. Structural and Multidisciplinary Optimization 62(5):2467–2479. 10.1007/s00158-020-02618-z
  * Wang et al (2008) Wang DW, Li F, Liu M, et al (2008) 3D aperiodic hierarchical porous graphitic carbon material for high-rate electrochemical capacitive energy storage. Angewandte Chemie International Edition 47(2):373–376. 10.1002/ange.200702721
  * Wang et al (2017) Wang Y, Chen C, Xie H, et al (2017) 3D-printed all-fiber Li-ion battery toward wearable energy storage. Advanced Functional Materials 27(43):1703,140
  * Wang et al (2023) Wang Z, Wang H, Xiao H, et al (2023) Enhancing heat dissipation and mass transfer of oxygen gas flow channel in a proton exchange membrane fuel cell using multiobjective topology optimization. International Journal of Hydrogen Energy 48(83):32,495–32,511. https://doi.org/10.1016/j.ijhydene.2023.05.023
  * Wathen (2015) Wathen AJ (2015) Preconditioning. Acta Numerica 24:329–376. 10.1017/s0962492915000021
  * Watts and Tortorelli (2016) Watts S, Tortorelli DA (2016) An $n$-material thresholding method for improving integerness of solutions in topology optimization. International Journal for Numerical Methods in Engineering 108(12):1498–1524. 10.1002/nme.5265
  * Xu et al (2019) Xu M, Reichman B, Wang X (2019) Modeling the effect of electrode thickness on the performance of lithium-ion batteries with experimental validation. Energy 186:115,864. https://doi.org/10.1016/j.energy.2019.115864, URL https://www.sciencedirect.com/science/article/pii/S0360544219315361
  * Xu et al (2010) Xu S, Cai Y, Cheng G (2010) Volume preserving nonlinear density filter based on heaviside functions. Structural and Multidisciplinary Optimization 41:495–505. 10.1007/s00158-009-0452-7
  * Xue et al (2015) Xue N, Du W, Martins JR, et al (2015) Lithium-ion batteries: Thermomechanics, performance, and design optimization. In: Handbook of Clean Energy Systems, John Wiley & Sons Ltd., vol. 6. Wiley Online Library, p 2849–2864, 10.1002/9781118991978.hces225
  * Yaji et al (2018) Yaji K, Yamasaki S, Tsushima S, et al (2018) Topology optimization for the design of flow fields in a redox flow battery. Structural and multidisciplinary optimization 57(2):535–546. 10.1007/s00158-017-1763-8
  * Yang et al (2022) Yang J, Li Y, Mijailovic A, et al (2022) Gradient porosity electrodes for fast charging lithium-ion batteries. Journal of Materials Chemistry A 10(22):12,114–12,124
  * Zadin et al (2011) Zadin V, Brandell D, Kasemägi H, et al (2011) Finite element modelling of ion transport in the electrolyte of a 3D-microbattery. Solid State Ionics 192(1):279–283
  * Zekoll et al (2018) Zekoll S, Marriner-Edwards C, Hekselman AKO, et al (2018) Hybrid electrolytes with 3D bicontinuous ordered ceramic and polymer microchannels for all-solid-state batteries. Energy Environ Sci 11:185–201
  * Zeng et al (2014) Zeng W, Shu L, Li Q, et al (2014) Fiber-based wearable electronics: A review of materials, fabrication, devices, and applications. Advanced Materials 26(31):5310–5336
  * Zheng et al (2012) Zheng H, Li J, Song X, et al (2012) A comprehensive understanding of electrode thickness effects on the electrochemical performances of Li-ion battery cathodes. Electrochimica Acta 71:258–265